## 3 Model Problems and Numerical Results
### Model Problem 2
#### 3.3.2 Finite Element Method  
Because of the nonlinear stress-strain relation in the plastic domain, a Newton-Raphson iterative solver is used. The plastic region is governed by the Drucker-Prager yield criterion with a small amount of isotropic linear hardening for numerical stability reasons. An incremental load approach is used starting with a force of 0 N until the downward force resulting from the slope's weight is reached. The methods used to solve the slope stability problem are based on ([46], Chapter 2 Section 4, and Chapter 7 Sections 3 and 4). For this case, the system equation takes the following form:  
\[\mathbf{K}\Delta\underline{\mathbf{u}}=\underline{\mathbf{r}} \tag{34}\]  
where \(\Delta\underline{\mathbf{u}}\) stands for the displacement increment. The vector \(\underline{\mathbf{r}}\) is the residual,  
\[\underline{\mathbf{r}}=\mathbf{f}+\Delta\mathbf{f}-\underline{\mathbf{q}}, \tag{35}\]  
where \(\mathbf{f}\) stands for the sum of the external force increments applied in the previous steps, \(\Delta\underline{\mathbf{f}}\) for the applied load increment of the current step and \(\underline{\mathbf{q}}\) for the internal force resulting from the stresses  
\[\underline{\mathbf{q}}=\int_{\Omega}\mathbf{B}^{T}\sigma d\Omega. \tag{36}\]  
First, the displacement increment of all the nodes is computed according to Equation (34), with an initial system stiffness matrix \(\mathbf{K}\) resulting from the assembly of the element stiffness matrix \(\mathbf{K}^{\mathbf{e}}\), computed by means of a Gauss quadrature, i.e.,  
\[\mathbf{K}^{\mathbf{e}}=\int_{\Omega}\mathbf{B}^{T}\mathbf{D}^{\mathrm{ep}} \mathbf{B}d\Omega. \tag{37}\]  
Here \(\mathbf{D}^{\mathrm{ep}}\) denotes the elastoplastic constitutive matrix. The initial state of \(\mathbf{D}^{\mathrm{ep}}\) is defined as  
\[\mathbf{D}^{\mathrm{ep}}=\frac{E(1-\
u)}{(1+\
u)(1-2\
u)}\begin{bmatrix}1& \frac{\
u}{(1-\
u)}&0\\ \frac{\
u}{(1-\
u)}&1&0\\ 0&0&\frac{(1-2\
u)}{2(1-\
u)}\end{bmatrix}. \tag{38}\]  
Secondly, the strain increment \(\Delta\varepsilon\) is computed as  
\[\Delta\varepsilon=\mathbf{B}\Delta\underline{\mathbf{u}}. \tag{39}\]  
Finally, the nonlinear stress-strain relationship,  
\[d\sigma=\mathbf{D}^{\mathrm{ep}}d\varepsilon, \tag{40}\]  
is integrated by means of a backward Euler method. The backward Euler method essentially acts as an elastic predictor-plastic corrector; an initial stress state that is purely elastic is computed and then projected in the direction of the yield surface to obtain the plastic stress state. Due to the implicit nature of the integrated stress-strain relation, this equation must be supplemented with the integrated form of the hardening rule and the yield condition. This system of nonlinear equations is then solved with an iterative Newton-Raphson method. Afterwards, the consistent tangent stiffness matrix is computed, see [47]. This matrix is then used to compute the updated element stiffness matrix, Equation (37),resulting in an updated system stiffness matrix \(\mathbf{K}\). The inner iteration step of solving the stress-strain relation and the updated system stiffness matrix is repeated for each outer iteration step which solves Equation (34). The outer step consists in balancing the internal forces with the external ones as to satisfy the residual, which in our case equals \(10^{-4}\) multiplied with the load increment. The procedure used is incremental-iterative, relying on the iterative Newton-Raphson method. This process is repeated for each load increment.