This special issue is a nice opportunity to honor Professor Radim Blaheta, a well-known Czech numerical mathematician. It was supported by his former collaborators, colleagues, friends, and students. Some of them have also contributed to this issue.

This paper is devoted to the extension of a quasi-Newton/variable preconditioning (QNVP) method to non-smooth problems, motivated by elasto-plastic models. Two approaches are discussed: the first one is carried out via regularized approximations of the nonsmooth problem, and the second one gives an extension to nonsmooth operators in order to be applied directly. Convergence analysis is presented for both variants. Then these abstract methods are applied to elasto-plasticity where two different variants of QNVP are investigated and combined with the deflated conjugate gradient and aggregation-based algebraic multigrid methods. The convergence results are illustrated on numerical examples in 3D inspired by real-life problems, and they demonstrate that the suggested QNVP methods are competitive with the standard Newton method. Well-documented Matlab codes on elasto-plasticity are used and enriched by the suggested methods.

This paper is motivated by the limit load, limit analysis and shear strength reduction methods, which are commonly employed in geotechnical stability analysis or similar applications. The aim is to make these methods more approachable by introducing a unified framework based on abstract convex optimization and its parametric studies. We establish suitable assumptions on the abstract problems that capture the selected features of these methods and facilitate rigorous theoretical investigation. Further, we propose continuation techniques tailored to the resulting parametric problem formulations and show that the developed abstract framework could also be useful outside the domain of geotechnical stability analysis. The main results are illustrated with analytical and numerical examples. The numerical example deals with a 3D slope stability problem.

This paper addresses the solution of slope stability problems in 3D using the finite element method and incremental procedures like the shear strength reduction or limit load methods. We build on Mohr-Coulomb plasticity, Davis' modifications of the non-associated plastic flow rule and recent mathematical results, which relate the factor of safety (FoS) with convex optimization. A complex solution concept is presented in detail and completed with in-house developed, publicly available open-source MATLAB codes. The concept consists of a combination of indirect continuation techniques, inexact Newton-like solvers and deflated Krylov methods with preconditioners. Further, mesh adaptivity is used to reduce overestimation of FoS and determine failure zones more accurately. The solution concept is tested on slope stability benchmarks in 3D and its efficiency is demonstrated. Numerical results are validated against either literature or software COMSOL Multiphysics.

This work continues an earlier study, aimed at determining the threshold of elastic behaviour for specimens at the microscale. Given the size-dependent nature of the response the theory makes use of a model of rate-independent strain-gradient plasticity to establish conditions for lower and upper bounds to the threshold. The model is based on a family of yield or dissipation functions. In the earlier work bounds were established for one member of the family, in which the dissipation depends on the root mean square of the magnitude of plastic strain and its gradient. In this work bounds are derived based on a second alternative, in which the dependence is linear. These two alternatives are most relevant in determining length scales through fits to experimental results. The second aim of this work is to carry out comparisons between the bounds obtained theoretically, with experimental data on initial yielding for torsion tests. The correlations are in some cases reasonable, while in others there are discrepancies. Generally they illustrate the complexities - for example in accurate experimental determination of the elastic threshold, and of the length scale - in predicting aspects of behaviour at the microscale.

This contribution is focused on severely nonlinear systems of equations with nonsmooth operators. A continuation Newton method with a smoothing operator is suggested and its convergence analyzed. Then the method is applied to elasto-plasticity with hardening. Finally, another continuation method convenient for an elastic-perfectly plastic problem is introduced and used for finding the so-called limit load. Dedication: This paper is dedicated to the memory of Professor Owe Axelsson. The presented results arise from joint work with Owe.

The Strategic Research Agenda (SRA; https://www.ejp-eurad.eu/publications/eurad-sra) of the European Joint Programme on Radioactive Waste Management (EURAD; https://www.ejp-eurad.eu/) describes the scientific and technical domains and sub-domains and knowledge management needs of common interest between EURAD participant organizations. Theme number 7 is entitled “Performance assessment, safety case development and safety analyses.” A list of research and development priorities and activities of common interest to be addressed within EURAD for theme 7 have been established. Amongst others, the Understanding and modelling of multi-physical Thermo-Hydro-Mechanical-Chemical coupled processes (THMC) occurring in radioactive waste disposal is a major and permanent issue to support optimization of design and safety case abstraction. To tackle this challenge a research work package entitled “DONUT: Development and improvement of numerical methods and tools for modelling coupled processes” has been conducted within the EURAD join programming initiative. The purpose of this work package is to improve/develop methods or numerical tools in order to go a step further in development of (i) relevant, performant and cutting-edge numerical methods that can easily be implemented in existing or new tools, in order to carry out high-performance computing to facilitate the study of highly coupled processes in large systems, (ii) numerical scale transition schemes for coupled processes, (iii) innovative numerical methods to carry out uncertainty and sensitivity analyses. In this paper the work carried out within the DONUT work package is put in perspective regarding the existing concept and literature on the field. It does not pretend to be exhaustive but rather to put emphasis on particular issues tackled during the project.

Quasi-Newton iterations are constructed for the finite element solution of small-strain nonlinear elasticity systems in 3D. The linearizations are based on spectral equivalence and hence considered as variable preconditioners arising from proper simplifications in the differential operator. Convergence is proved, providing bounds uniformly w.r.t. the FEM discretization. Convenient iterative solvers for linearized systems are also proposed. Numerical experiments in 3D confirm that the suggested quasi-Newton methods are competitive with Newton's method.

This paper focuses on the analysis and the solution of the saddle-point problem arising from a three-field formulation of Biot's model of poroelasticity, discretized in time by the implicit Euler method. A block diagonal-preconditioner, based on the Schur complement, is analyzed on a functional level and compared with two other block-diagonal preconditioners having a similar structure. The problem is discretized in space using mixed finite elements and solved with appropriate iterative solvers, incorporating the investigated preconditioners. The solvers are tested on numerical examples inspired by geotechnical practice, with particular attention devoted to the solvers' robustness concerning strong heterogeneity in permeability.

This paper is focused on continuation techniques and Newton-like methods suggested for numerical determination of safety factors within stability assessment. Especially, we are interested in the stability of slopes and related limit load and shear strength reduction methods. We build on computational plasticity and the finite element method, but we mainly work on an algebraic level to be the topic understandable for broader class of scientists and our algorithms more transparent. The presented algorithms are based on the associated plasticity to be more robust. For non-associated models, we use Davis-type approximations enabling us to apply the associated approach. A particular attention is devoted to the Mohr-Coulomb elastic-perfectly plastic constitutive problem. On this example, we explain some important features of the presented methods which are beyond the algebraic settings of the problems. We also summarize the Mohr-Coulomb constitutive solution and some implementation details.

In this paper, an optimization variant of the shear strength reduction method is introduced and used for the solution of embankment stability problems with unconfined seepage. The optimization framework is based on approximations of non-associated Mohr-Coulomb plastic models with associated ones, especially by using various Davis' approaches. Next, the finite element method is considered and mesh adaptive solution concepts are developed for both the unconfined seepage and stability problems. In-house codes in Matlab are used for their implementation. Finally, two numerical examples inspired by geotechnical practice are investigated in order to demonstrate the accuracy of the optimization framework and to evaluate three different Davis' approaches. The results are compared with commercial codes in Plaxis and Comsol Multiphysics. (c) 2023 Elsevier Ltd. All rights reserved.

The shear strength reduction method (SSRM) is a standard method in slope stability enabling to determine the factor of safety and related failure zones. In case of the non-associated Mohr-Coulomb model, the method can oscillate with respect to the refinement of a finite element mesh. To suppress this drawback, the non-associated model is approximated by the associated one such that the strength parameters are reduced by using a function depending on a scalar factor and on the effective friction and dilatancy angles. This modification (MSSRM) can be easily implemented in commercial codes like Plaxis or Comsol Multiphysics. Next, an optimization approach to the modified SSRM (OPT-MSSRM) is introduced. It is shown that the optimization problem is well-defined and can be analyzed by variational principles. For its solution, a regularization method is combined with mesh adaptivity and implemented in Matlab. The SSRM, MSSRM and OPT-MSSRM methods are compared on numerical examples representing a case study of a real heterogeneous slope.

This paper is focused on a new optimization variant of the shear strength reduction (OPT-SSR) in non-associated Mohr-Coulomb plasticity. The OPT-SSR method mimics the limit analysis problem and enables to compute the factor of safety without performing an elasto-plastic analysis. It is shown that this optimization problem is well defined and closely related to recently developed Davis approaches used in combination with the standard SSR method. Next, the duality between the static and kinematic principles of OPT-SSR is introduced. For the numerical solution, a regularization method is suggested.

This paper is focused on the definition, analysis and numerical solution of a new optimization variant (OPT) of the shear strength reduction (SSR) problem with applications to slope stability problems. This new variant is derived on the basis of recent results by Tschuchnigg et al. 2015, where limit analysis and a modified Davis approach were used for approximation of the standard SSR method. The OPT-SSR method computes the factor of safety without performing an elasto-plastic analysis, similarly as in limit analysis. It is shown that this optimization problem is well-defined. Next, the duality between the static and kinematic principles of OPT-SSR is derived. For the numerical solution, a regularization method is introduced and analyzed. This method is combined with the finite element method, mesh adaptivity and a damped Newton method. In-house codes (Matlab) are used for the implementation of this solution concept. Finally, two slope stability problems are considered, one of which follows from analysis of a real slope. The softwares packages Plaxis and Comsol Multiphysics are used for comparison of the results.

Longest-edge bisection algorithms are often used for local mesh refinements within the finite element method in 2D. In this paper, we discuss and describe their conforming variant. A particular attention is devoted to the so-called Right-Triangulated Irregular Network (RTIN) based on isosceles right triangles and its tranformation to more general domains. We suggest to combine RTIN with a balanced quadrant tree (QuadTree) decomposition. This combination does not produce hanging nodes within the mesh refinements and could be extended to tetrahedral meshes in 3D.

This work is concerned with an abstract inf-sup problem generated by a bilinear Lagrangian and convex constraints. We study the conditions that guarantee no gap between the inf-sup and related sup-inf problems. The key assumption introduced in the paper generalizes the well-known Babuska-Brezzi condition. It is based on an inf-sup condition defined for convex cones in function spaces. We also apply a regularization method convenient for solving the inf-sup problem and derive a computable majorant of the critical (inf-sup) value, which can be used in a posteriori error analysis of numerical results. Results obtained for the abstract problem are applied to continuum mechanics. In particular, examples of limit load problems and similar ones arising in classical plasticity, gradient plasticity and delamination are introduced.

This work is concerned with the purely dissipative version of a well-established model of rate-independent strain-gradient plasticity. In the conventional theory of plasticity the approach to determining plastic flow is local, and based on the stress distribution in the body. For the dissipative problem of strain-gradient plasticity such an approach is not valid as the yield function depends on microstresses that are not known in the elastic region. Instead, yield and plastic flow must be considered at the global level. This work addresses the problem of determining the elastic threshold by formulating primal and dual versions of the global problem and, motivated by techniques used in limit analysis for perfect plasticity, establishing conditions for lower and upper bounds to the threshold. The general approach is applied to two examples: of a plate under plane stress, and subjected to a prescribed displacement; and of a bar subjected to torsion.