A Selection of Benchmark Problems in Solid Mechanics and Applied Mathematics

Abstract

In this contribution we provide benchmark problems in the field of computational solid mechanics. In detail, we address classical fields as elasticity, incompressibility, material interfaces, thin structures and plasticity at finite deformations. For this we describe explicit setups of the benchmarks and introduce the numerical schemes. For the computations the various participating groups use different (mixed) Galerkin finite element and isogeometric analysis formulations. Some programming codes are available open-source. The output is measured in terms of carefully designed quantities of interest that allow for a comparison of other models, discretizations, and implementations. Furthermore, computational robustness is shown in terms of mesh refinement studies. This paper presents benchmarks, which were developed within the Priority Programme of the German Research Foundation ‘SPP 1748 Reliable Simulation Techniques in Solid Mechanics—Development of Non-Standard Discretisation Methods, Mechanical and Mathematical Analysis’.

1 Introduction

Solving partial differential equations on complex geometries is perhaps one of the most important scientific achievement of the last decades. Analytical or manufactured solutions of such differential equations, e.g. from engineering or economics, is in most cases not available. Therefore, computer-aided numerical algorithms play an important role. At this point we mention that the rapid development of the available computing power contributes significantly to the success of these numerical methods. An end to this progress is currently not conceivable and becomes even more important due to the capabilities to address challenging multiphysics applications. Popular numerical methods are, for example, mesh-free methods, the finite volume method, finite differences, isogeometric analysis and—with possibly the greatest impact—the finite element method.

In this field, the German Priority Program 1748 (DFG SPP 1748) Reliable Simulation Techniques in Solid Mechanics. Development of Non-standard Discretization Methods, Mechanical and Mathematical Analysis is located. The main goal of this priority program was the development of modern non-conventional discretization methods based e.g. on mixed finite elements or discontinuous Galerkin formulations, including mathematical analysis for geometrically as well as physically non-linear problems in the areas of incompressibility, anisotropies and discontinuities to name a few. Typical problems especially in the field of geometric and material nonlinearities are for instance insufficient stress approximations due to unsuitable approximation spaces as well as weak convergence behavior due to stiffening effects or mesh distortions. Similar problems arise in crack or contact problems, where the local determination of discontinuities as well as their development play an important role in many fields of application.

This paper presents benchmark collections, which were compiled within the SPP 1748. The benchmarks were designed and shall serve as future reference for comparisons with other discretizations, nonlinear and linear solution algorithms. In the first two benchmarks we show results on hyperelastic and elasto-plastic problems at finite strains. Different (mixed) finite element technologies and p-FEM are used to obtain convergent results for certain displacement and stress values for Cook’s membrane and an incompressible block problem. In the third benchmark problem (a thin as well as a thick plate) the convergence behavior was investigated. Here, both shell and continuum finite element technologies are used. The fourth benchmark concentrates on material interfaces that lead to high gradients or singularities. Therein, two discretization methods, finite elements and isogeometric analysis are compared. In the fifth benchmark, a phase-field method is applied for two- and three-dimensional propagating fractures. Lastly, again a phase-field fracture benchmark is provided. Therein, the fracture is subject to a constant pressure. For this setting, even manufactured solutions are known and taken as reference values. In all configurations, all necessary data are provided to reproduce our findings. These data include reference values, computations on mesh hierarchies and various quantities of interest.

2 Cook’s Membrane: Hyperelasticity and Finite Elastoplasticity

In the first part of this benchmark, the Cook’s membrane is investigated by applying a hyperelastic material model (Sect. 2.3) while in the second part a finite $J_2$ elastoplastic material model is used (Sect. 2.4). To this end, different element formulations are investigated by showing displacement and stress convergence studies. Furthermore, the locking effect due to the incompressibilty contraint is examined.

2.1 Geometry and Boundary Conditions

Cook’s membrane, a tapered cantilever see Cook [1], combines bending and shearing. The left-hand side of the domain is clamped and a constant shear load in vertical direction is applied on the right-hand side. The thickness of the Cook’s membrane is chosen to be 1 mm. The geometry and boundary conditions are illustrated in Fig. 1.

Fig. 1

Cook’s membrane (all dimensions in mm)

2.2 Element Formulations

Five different element formulations are used in this study. Two formulations are based on triangular (2D) or tetrahedral (3D) shaped elements. The first one applies a pure displacement based formulation with quadratic shape functions (in the following denoted as $T_2$) and the second one is a mixed formulation with an additional constant field for the pressure and volume dilatation (denoted with $T_2P_0$), see e.g. Boffi et al. [2]. The third and fourth formulation are based on hexahedral shaped elements with a first order polynomial interpolation of the displacements (represented by $H_1$) and a mixed formulation with a piecewise constant pressure field (denoted with $H_1P_0$), see e.g. Boffi et al. [2]. Furthermore, the three-dimensional problem is discretized also with hexahedral elements, which are used within a p-extension (p-FEM) based on the trunk space utilizing hierarchic shape functions, see Szabó and Babuška [3] and Düster et al. [4].

2.3 Cook’s Membrane for Hyperelasticity

2.3.1 Hyperelastic Material Model

In the following we use two polyconvex strain-energy functions. The first function reads

$$\begin{aligned} \psi _1=\frac{\mu }{2} (I_{{{{\varvec{C}}}}} - 3)+ \frac{\lambda }{4} (J^2 - 1) - \left(\frac{\lambda }{2}+\mu \right)\ln {J} , \end{aligned}$$

(2.1)

where $\lambda$ and $\mu$ constitute the Lamé parameters, $J=\det {{{\varvec{F}}}}$ the determinant of the deformation gradient ${{{\varvec{F}}}}$ and $I_{{{{\varvec{C}}}}}={{\text {tr}}}{{{\varvec{C}}}}$ the first principal invariant of the right Cauchy-Green tensor ${{{\varvec{C}}}}={{{\varvec{F}}}}^T {{{\varvec{F}}}}$, see e.g. Ciarlet [5]. The 2nd Piola Kirchhoff stress tensor ${{{\varvec{S}}}}=2\partial _{{{\varvec{C}}}}\psi$ can be computed from the strain-energy function $\psi _1$ as follows

$$\begin{aligned} {{{\varvec{S}}}}= \frac{\lambda }{2} (J^2 - 1 ){{{\varvec{C}}}}^{-1} + \mu (\mathbf{1}- {{{\varvec{C}}}}^{-1}) , \end{aligned}$$

(2.2)

here $\mathbf{1}$ denotes the second-order identity tensor. The material parameters ($\lambda$, $\mu$) as well as the load $p_0$ are given in Table 1.

Table 1 Material parameters and applied load for $\psi _1$

The second strain-energy function under consideration is based on a split into an isochoric U(J) and unimodular $W(I_{{\widetilde{{{{\varvec{C}}}}}}},II_{{\widetilde{{{{\varvec{C}}}}}}})$ part

$$\begin{aligned} \psi _2=U(J)+W(I_{{\widetilde{{{{\varvec{C}}}}}}},II_{{\widetilde{{{{\varvec{C}}}}}}}) , \end{aligned}$$

(2.3)

see Hartmann and Neff [6] and Netz et al. [7]. The principal invariants of the unimodular right Cauchy-Green tensor

$$\begin{aligned} {\widetilde{{{{\varvec{C}}}}}}=(\det {{{\varvec{C}}}})^{-1/3}{{{\varvec{C}}}}\end{aligned}$$

(2.4)

are given with

$$\begin{aligned} I_{{\widetilde{{{{\varvec{C}}}}}}}={{\text {tr}}}{\widetilde{{{{\varvec{C}}}}}} \quad \hbox {and} \quad II_{{\widetilde{{{{\varvec{C}}}}}}}={{\text {tr}}}({{\text {Cof}}}{\widetilde{{{{\varvec{C}}}}}}) . \end{aligned}$$

(2.5)

Following Hartmann and Neff [6] we choose the isochoric and unimodular parts of the strain energy as

$$\begin{aligned} U(J)=\frac{\kappa }{50}(J^5+J^{-5}-2) \end{aligned}$$

(2.6)

and

$$\begin{aligned} W(I_{{\widetilde{{{{\varvec{C}}}}}}},II_{{\widetilde{{{{\varvec{C}}}}}}})= & {} \alpha (I_{{\widetilde{{{{\varvec{C}}}}}}}^3-27)+c_{10}(I_{{\widetilde{{{{\varvec{C}}}}}}}-3) \nonumber \\&+\,c_{01}(II_{{\widetilde{{{{\varvec{C}}}}}}}^{3/2}-3 \sqrt{3}) . \end{aligned}$$

(2.7)

The 2nd Piola Kirchhoff stress tensor for $\psi _2$ can be computed as follows

$$\begin{aligned} {{{\varvec{S}}}}= & {} J U'(J) {{{\varvec{C}}}}^{-1} + 2J^{-2/3} \nonumber \\&\left( (W_{,1}+W_{,2} I_{{\widetilde{{{{\varvec{C}}}}}}}){{{\varvec{I}}}}- W_{,2} {\widetilde{{{{\varvec{C}}}}}} \right. \nonumber \\&\qquad \left. -\, \frac{1}{3}(W_{,1} I_{{\widetilde{{{{\varvec{C}}}}}}}+ 2W_{,2} II_{{\widetilde{{{{\varvec{C}}}}}}}) {\widetilde{{{{\varvec{C}}}}}}^{-1} \right) \end{aligned}$$

(2.8)

with

$$\begin{aligned}&W_{,1}=\partial W/\partial I_{{\widetilde{{{{\varvec{C}}}}}}}=3\ \alpha \ I_{{\widetilde{{{{\varvec{C}}}}}}}^2 \quad \hbox {and} \nonumber \\&W_{,2}=\partial W/\partial II_{{\widetilde{{{{\varvec{C}}}}}}}=c_{01}\ (3/2)\ II_{{\widetilde{{{{\varvec{C}}}}}}}^{(1/2)} . \end{aligned}$$

(2.9)

The material parameters used for the second strain-energy function ($\alpha ,c_{10},c_{01},\kappa$) as well as the load $p_0$ are given in Table 2.

Table 2 Material parameters and applied load for $\psi _2$

Furthermore, a Fortran77 code for the computation of the 2nd Piola Kirchhoff stress tensor and the material tangent ${\mathbb {C}}=2\partial _{{{\varvec{C}}}}{{{\varvec{S}}}}$ will be provided online under (https://www.doi.org/10.5281/zenodo.4014897 (see also Korelc and Stupkiewicz [14])) with the subroutines:

Here, the input parameters are the deformation gradient (F as $3\times 3$ field) as well as the material parameters. The output parameters are the 2nd Piola Kirchhoff stress tensor (S as $3\times 3$ field) and the material tangent (Cmat as $3\times 3\times 3\times 3$ field). For further details on the material model and the solution of hyperelastic problems, the reader is referred to Hartmann and Neff [6], Düster et al. [8] and Wriggers [9].

2.3.2 Analysis for the 2D Plane Strain Case

In the following the convergence of the spatial discretization in terms of the vertical tip displacement $u_y$ of the upper right node (point A, see Fig. 1) is compared. The meshes are uniformly refined by increasing the number of elements per edge. An exemplary mesh with triangular elements is depicted in Fig. 2a. Furthermore, for $\psi _1$ a mesh consisting of 13 hexahedral elements with two refinements towards the singularity is used and for $\psi _2$ a mesh with 15 hexahedral elements with three refinements, as depicted in Fig. 2b and c, respectively. The discretization is based on the trunk space utilizing 3D hierarchic shape functions. The plane strain conditions are enforced by suppressing the displacements in z-direction at the front ($z=0.5\,\hbox {mm}$) and reverse side ($z=-0.5\,\hbox {mm}$) of the hexahedral meshes. In the following, we will consider both strain-energy functions $\psi _1$ and $\psi _2$. An uniform and isotropic increase of the polynomial degree $p=1,2,3,\ldots$ of the shape functions yields the results depicted in Fig. 3a and b, which are plotted together with the displacements computed with the formulations $T_2$ and $T_2P_0$ for triangular elements. Furthermore, the results of the vertical displacement of point A are given in Tables 3 and 4.

Fig. 2

Spatial discretization: a triangular mesh $T_2$ and $T_2P_0$ with 169 nodes b hexahedral mesh consisting of 13 elements with two refinements towards the singularity c hexahedral mesh consisting of 15 elements with three refinements towards the singularity

Fig. 3

Cook’s membrane: a convergence of displacement $u_y$ at point A for $\psi _1$ b convergence of displacement $u_y$ at point A for $\psi _2$

Table 3 Convergence of the displacement component $u_y$ of point A for $\psi _1$—for the p-FEM the hexahedral mesh with 13 elements was used with $p=1,\ldots ,14$

Table 4 Convergence of the displacement component $u_y$ of point A for $\psi _2$—for the p-FEM the hexahedral mesh with 15 elements was used with $p=1,\ldots ,9$

Figure 4 shows the contour plot of the $\sigma _{xx}$ stresses for the strain-energy function $\psi _1$. Here, the results of the $T_2$ and the $T_2P_0$ element are depicted in Fig. 4a and b, respectively. Figure 4c shows the contour plot of the $\sigma _{xx}$ stresses for the p-FEM analysis using 13 hexahedral elements and the order of the Ansatz $p=11$. In Fig. 5 the stress distribution for the strain-energy function $\psi _2$ is depicted. It shows the $\sigma _{xx}$ stresses for the $T_2$ and the $T_2P_0$ element in Fig. 5a and b, respectively. To this end, the stress values are linearly extrapolated from the integration points to the nodal points. Stress oscillations can be recognized in the case of the displacement based element. The $\sigma _{xx}$ stress distribution of the p-FEM solution is plotted in Fig. 5c. Here, a mesh with 1026 hexahedral elements is used and $p=9$ to achieve a smooth solution of the stress distribution.

Fig. 4

Contour plot of the $\sigma _{xx}$ stresses for $\psi _1$: a $T_2$ and b $T_2P_0$ with 4704 DOFs c p-FEM solution with 104 hexahedral elements and $p=5$

Fig. 5

Contour plot of the $\sigma _{xx}$ stresses for $\psi _2$: a $T_2$ and b $T_2P_0$ with 4704 DOFs c p-FEM solution with 1026 hexahedral elements and $p=9$

2.3.3 Analysis for the 3D Case

For the three-dimensional problem we choose the same constitutive model, geometry and boundary conditions analogous to the 2D setting. However, instead of plane strain conditions we consider stress-free out-of-plane surfaces, i.e. the front and reverse side of Cook’s membrane are not constrained. Again, an uniform and isotropic p-refinement is performed utilizing the hexahedral mesh with 13 elements for $\psi _1$ and the hexahedral mesh with 15 elements for $\psi _2$ as depicted in Fig. 2b and c, respectively. The symmetry in z-direction is exploited by choosing one element layer of thickness 0.5 mm. In Fig. 6a and b the convergence of the displacement component $u_y$ at point A is depicted together with the results obtained with the $T_2$ and $T_2P_0$ formulations for tetrahedral elements applying also a mesh layer thickness of 0.5 mm. The corresponding numbers are also depicted in Tables 5 and 6.

Fig. 6

Cook’s membrane: a convergence of displacement $u_y$ at point A for $\psi _1$ b convergence of displacement $u_y$ at point A for $\psi _2$

Table 5 Convergence of the displacement component $u_y$ at point A for $\psi _1$—for the p-FEM the hexahedral mesh with 13 elements was used with $p=1,\ldots ,10$

Table 6 Convergence of the displacement component $u_y$ at point A for $\psi _2$—for the p-FEM the hexahedral mesh with 15 elements was used with $p=1,\ldots ,12$

Influence of incompressibility: In this section, we investigate the influence of the incompressibility level on the convergence of the different discretizations. To this end, we keep the same material parameters and load as before and only change the Lamé’s parameter $\lambda$ for the first free energy function $\psi _1$ and the bulk modulus $\kappa$ for the second free energy function $\psi _2$. The values of Poisson’s ratio are approximated considering the corresponding relations for homogeneous isotropic linear materials, see Tables 7 and 8. Consequently, Poisson’s ratio for $\psi _1$ is computed as

$$\begin{aligned} \nu = \frac{\lambda }{2 (\lambda + \mu )} \end{aligned}$$

(2.10)

and for $\psi _2$ as

$$\begin{aligned} \nu = \frac{3 \kappa - 2G}{6 \kappa + 2G} , \end{aligned}$$

(2.11)

whereas according to Hartmann and Neff [6] the shear modulus is defined as

$$\begin{aligned} G = 54 \alpha + 2 c_{10} + 3 \sqrt{3} c_{01} . \end{aligned}$$

(2.12)

Table 7 Lamé’s parameter $\lambda$ and Poisson’s ratio $\nu$ for $\psi _1$

Table 8 Bulk modulus $\kappa$ and Poisson’s ratio $\nu$ for $\psi _2$

In order to investigate the influence of the increasing incompressibility level, we study the convergence of the tip displacements. The reference solutions of the tip displacement $u_y$ at point A are computed by the p-FEM discretization using 13 elements for $\psi _1$ and 15 elements for $\psi _2$. The corresponding values are given in Tables 9 and 10. Figure 7 depicts the convergence of the normalized displacement $u_y / u_{y,\mathrm {ref}}$ for both strain-energy functions $\psi _1$ and $\psi _2$ computed with the different finite element formulations.

Fig. 7

Locking behavior considering the normalized displacement $u_y / u_{y,\mathrm {ref}}$ of point A computed with the different finite element formulations

Table 9 Reference tip displacement considering strain energy $\psi _1$

Table 10 Reference tip displacement considering strain energy $\psi _2$

2.4 Cook’s Membrane for Elastoplasticity

In this section we use the same setup as depicted in Fig. 1. Furthermore, we consider only the 3D case.

2.4.1 Finite $J_2$ Elastoplastic Material Model

First, a short summary of the governing equations for the underlying material model will be given which is based on the $J_2$ flow theory of elastoplasticity for finite strain including a nonlinear isotropic hardening. For a detailed overview about the theory the interested reader is referred to Simo [10,11,12] and Simo and Miehe [13].

For the material model the deformation gradient $\varvec{F}$ is decomposed into an elastic $\varvec{F}_e$ and a plastic $\varvec{F}_p$ part by applying a multiplicative split

$$\begin{aligned} \varvec{F} = \varvec{F}_e \varvec{F}_p. \end{aligned}$$

(2.13)

Assuming $J_2$ von Mises plasticity the plastic deformation only depends on the deviatoric part of the stress tensor. Based on the multiplicative split in Eq. (2.13) the elastic deformation gradient reads $\varvec{F}_e = \varvec{F} \varvec{F}_p^{-1}$, thus, the elastic right Cauchy-Green tensor $\varvec{C}_e$ can be formulated as

$$\begin{aligned} \varvec{C}_e = \varvec{F}_p^{-T} \varvec{C} \varvec{F}_p^{-1}, \end{aligned}$$

(2.14)

and the elastic left Cauchy-Green tensor $\varvec{b}_e$ as

$$\begin{aligned} \varvec{b}_e = \varvec{F} \varvec{C}_p^{-1} \varvec{F}^T, \end{aligned}$$

(2.15)

where $\varvec{C}_p = \varvec{F}_p^T \varvec{F}_p$ defines the plastic right Cauchy-Green tensor.

An isotropic compressible Neo-Hookean material behavior is assumed for the elastic part of the deformation utilizing the following strain-energy function

$$\begin{aligned} \psi _e= & {} \frac{\lambda }{4} \left( III_{{{{\varvec{C}}}}_e} - 1 - \ln III_{{{{\varvec{C}}}}_e} \right) \nonumber \\&+ \,\frac{\mu }{2} \left( I_{{{{\varvec{C}}}}_e} - 3 - \ln III_{{{{\varvec{C}}}}_e} \right) . \end{aligned}$$

(2.16)

In Eq. (2.16), $\lambda$ and $\mu$ define the Lamé constants and $I_{{{{\varvec{C}}}}_e}$ and $III_{{{{\varvec{C}}}}_e}$ denote the first and third invariant of the elastic right Cauchy-Green tensor $\varvec{C}_e$ or the elastic left Cauchy-Green tensor $\varvec{b}_e$ as follows

$$\begin{aligned}&I_{{{{\varvec{C}}}}_e} = \mathrm {tr} \varvec{C}_e \quad \text {and} \quad III_{{{{\varvec{C}}}}_e} = \mathrm {det} \varvec{C}_e \quad \text {or} \nonumber \\&I_{{{{\varvec{b}}}}_e} = \mathrm {tr} \varvec{b}_e \quad \text {and} \quad III_{{{{\varvec{b}}}}_e} = \mathrm {det} \varvec{b}_e. \end{aligned}$$

(2.17)

The 1st and 2nd Piola-Kirchhoff as well as the Cauchy and Kirchhoff stress tensors read

$$\begin{aligned}&\varvec{P} = \frac{\partial \psi _e}{\partial \varvec{F}_e} , \quad \varvec{S} = 2 \frac{\partial \psi _e}{\partial \varvec{C}_e} , \nonumber \\&\varvec{\sigma } = \frac{2}{J_e} \varvec{F}_e \frac{\partial \psi _e}{\partial \varvec{C}_e} \varvec{F}_e^T , \quad \varvec{\tau } = 2 \frac{\partial \psi _e}{\partial \varvec{b}_e} \varvec{b}_e. \end{aligned}$$

(2.18)

The yield criterion of the plasticity model reads

$$\begin{aligned} \Phi \left( \varvec{\tau }, {\bar{\alpha }} \right) = \sqrt{\frac{3}{2} \varvec{s} : \varvec{s}} - K ({\bar{\alpha }} ) \quad \text {with} \quad \varvec{s} = \varvec{\tau } - \frac{1}{3} \mathrm {tr} (\varvec{\tau }) \varvec{1}. \end{aligned}$$

(2.19)

In Eq. (2.19) ${\bar{\alpha }}$ denotes the hardening variable, $\varvec{s}$ defines the deviatoric part of the Kirchhoff stress tensor $\varvec{\tau }$, and $K({\bar{\alpha }})$ describes the hardening which consists of a linear and an exponential part as follows

$$\begin{aligned} K({\bar{\alpha }}) = Y_0 + H {\bar{\alpha }} + \left( Y_\infty - Y_0 \right) \left( 1 - e^{-\delta {\bar{\alpha }}} \right) , \end{aligned}$$

(2.20)

see e.g. Simo [12]. Here, $Y_0$ denotes the initial yield stress, H the linear hardening parameter, $Y_\infty$ the saturation stress, and $\delta$ the hardening exponent. In accordance to Simo and Miehe [13], the evolution of the plastic flow as well as the internal variable reads

$$\begin{aligned} -\frac{1}{2} \mathcal {L}_v \varvec{b}_e = {\gamma } \frac{\partial \Phi }{\partial \varvec{\tau }} \varvec{b}_e \quad \hbox {and} \quad \dot{{\bar{\alpha }}} = - {\gamma } \frac{\partial \Phi }{\partial K({\bar{\alpha }})}. \end{aligned}$$

(2.21)

Here, $\mathcal {L}_v \varvec{b}_e$ denotes the Lie derivative of the elastic left Cauchy-Green tensor $\varvec{b}_e$ and $\gamma \ge 0$ is the non-negative plastic multiplier. Based on the principle of maximum dissipation, this follows by enforcing stationarity conditions for the Lagrange functional, which is constructed as an optimization problem with the constraint condition ($\Phi = 0$), depending on the Dissipation inequality and the yield criterion in relation with $\gamma$.

Furthermore, following Korelc and Stupkiewicz [14] Eq. (2.15) can be inserted into the plastic flow formula in Eq. (2.21) to obtain,

$$\begin{aligned} \dot{\varvec{C}}_p^{-1} = -2 {\gamma } \varvec{F}^{-1} \frac{\partial \Phi }{\partial \varvec{\tau }} \varvec{F} \varvec{C}_p^{-1} , \end{aligned}$$

(2.22)

which is then used for the algorithmic treatment of the model. Finally, the material model is complemented by the Kuhn-Tucker conditions

$$\begin{aligned} \Phi \le 0 , \quad {\gamma } \ge 0 , \quad \Phi {\gamma } = 0 , \end{aligned}$$

(2.23)

describing loading and unloading. For the integration of the internal variable $\dot{{\bar{\alpha }}}$, a backward Euler scheme is used. For the plastic variable $\dot{\varvec{C}}_p$ an exponential map algorithm within an implicit time integration scheme is applied as follows

$$\begin{aligned} {\bar{\alpha }}_{n+1}= & {} {\bar{\alpha }}_{n} + \Delta _t\gamma , \end{aligned}$$

(2.24)

$$\begin{aligned} (\varvec{C}_{p}^{-1})_{n+1}= & {} \varvec{F}_{n+1}^{-1} \; \text {exp} \left[ -2 {\Delta _t\gamma } \frac{\partial \Phi _{n+1}}{\partial \varvec{\tau }_{n+1}} \right] \varvec{F}_{n+1} (\varvec{C}_{p}^{-1})_{n}, \end{aligned}$$

(2.25)

which was introduced by Weber and Anand [15] and Eterovic and Bathe [16]. Here, $\Delta _t\gamma$ denotes the increment of the plastic multiplier. The material parameters used in this benchmark as well as the load $p_0$ are given in Table 11.

Table 11 Material parameters and applied load

As in the hyperelastic case the convergence of the spatial discretization in terms of the vertical tip displacement $u_y$ of the upper right node (point A, see Fig. 1) is studied. As initial discretization a mesh consisting of 104 hexahedral elements with two refinements towards the singularity at the top left corner is used, as depicted in Fig. 8a. Again, we take advantage of symmetry and apply just one layer of elements of thickness 0.5mm. The starting mesh for the $H_1$ and $H_1P_0$ formulations is depicted in Fig. 8b, whereas in case of the $T_2$ and $T_2P_0$ each hexahedral element is subdivided into five tetrahedral elements, taking also advantage of symmetry in thickness direction. In the case of the h-FEM approach, the elements are uniformly refined in in-plane direction. In contrast, for the p-FEM the discretization is based on the trunk space utilizing hierarchic shape functions. An uniform and isotropic increase of the polynomial degree $p=1,\ldots ,8$ of the shape functions yields the results depicted in Fig. 9a, which are plotted together with the displacements computed with the formulations $T_2$, $T_2P_0$, $H_1$ and $H_1P_0$. Furthermore, the results of the vertical displacement of point A are given in Table 12. Furthermore, a Fortran77 code for the computation of the 1st Piola Kirchhoff stress tensor and the material tangent ${\mathbb {A}}=\partial _{{{\varvec{F}}}}{{{\varvec{P}}}}$ will be provided online under (https://www.doi.org/10.5281/zenodo.4014897 (see also Korelc and Stupkiewicz [14])) with the subroutine:

Fig. 8

Spatial discretization: a a mesh consisting of 104 hexahedral elements with two refinements towards the singularity used for the p-FEM computations b the starting mesh for the $H_1$ and $H_1P_0$ formulations

Table 12 Convergence of the displacement component $u_y$ at point A

Von Mises stress convergence: In the following, we study the convergence behavior of the von Mises stress. In doing so, we use the same meshes as in the previous section. In Fig. 9b the von Mises stress at point B is plotted over the number of degrees of freedom. The corresponding numbers are also depicted in Table 13. Figure 10a shows the contour plot for the von Mises stress $\sigma _{vM}$ and Fig. 10b the equivalent plastic strain ${\bar{\alpha }}$ for $p=6$ at the last load step.

Table 13 Convergence of the von Mises stress $\sigma _{vM}$ at point B

Fig. 9

Cook’s membrane: a convergence of displacement $u_y$ at point A b convergence of von Mises stress $\sigma _{vM}$ at point B

Fig. 10

Cook’s membrane contour plots for p-FEM a von Mises stress $\sigma _{vM}$ b equivalent plastic strain ${\bar{\alpha }}$, for $p = 6$ at the last load step

Convergence study using different load steps: In the following we study the influence of the total number of load steps used to apply the load. To this end, the displacement $u_y$ at point A and the von Mises stress $\sigma _{vM}$ at point B are investigated. In doing so, we use the same mesh like in the previous section, as depicted in Fig. 8a applying a polynomial degree of $p = 5$ based on the trunk space utilizing hierarchic shape functions. In Fig. 11a and b, the displacements as well as the von Mises stress are plotted over the number of load steps chosen to apply the prescribed load. The corresponding numbers are also depicted in Table 14.

Fig. 11

Convergence study for different load steps using p-FEM, $p = 5$ a displacement $u_y$ at point A b von Mises stress $\sigma _{vM}$ at point B

Table 14 Convergence of the displacement component $u_y$ at point A and the von Mises stress $\sigma _{vM}$ at point B for different numbers of load steps using $p = 5$

2.4.2 Remark

The software used in this study are AceGen and AceFEM (Korelc [17, 18], Korelc and Wriggers [19]) as well as in-house codes AdhoC$^4$ (Düster and Kollmannsberger [20]) and AdhoC++ which are jointly developed at the Technical University of Munich and the Hamburg University of Technology.

3 Incompressible Block Under Constant Partial Load

In this benchmark, a cube is subject to a vertical load on part of its upper surface. This benchmark poses two challenges for a FEM simulation in a geometrically non-linear setting:

• two line-singularities induced by a jump in the Neumann boundary condition which intersect in a point,

• almost incompressible material.

The benchmark allows to asses the performance of (high-order) displacement-based and mixed finite elements.

3.1 Geometry and Boundary Conditions

The system is depicted in Fig. 12.

Fig. 12

Quasi-incompressible cube with partial load

Its dimensions are given in Table 15.

Table 15 Dimensions of the problem

At the upper surface of the block, i.e. at $z=h$, displacements are fixed in both x- and y-direction. The bottom of the block, i.e. at $z=0$, is fixed in z-direction.

Due to axial symmetry of loads and boundary conditions, computations were only carried out on the quarter of the block marked in light gray. To this end, homogeneous Dirichlet boundary conditions were additionally applied on the symmetry planes such that at $x = 0.5w$, displacements were fixed in x-direction and at $y = 0.5l$ displacements were fixed in y-direction.

Furthermore, the block is partially loaded with a constant surface load q on the area defined by a and b. The area in dark gray, thus, bears a mixed Dirichlet-Neumann boundary condition. These boundary conditions are chosen according to a similar test presented in Reese et al. [21].

3.2 Hyperelastic Material Model

We consider a deformation ${{\varvec{\varphi }}}({{{\varvec{X}}}},t)$ which maps points in the initial configuration to the current or deformed configuration. This deformation can be computed using the coordinates of the initial configuration and the displacement field: ${{\varvec{\varphi }}}({{{\varvec{X}}}},t)={{{\varvec{X}}}}-{{{\varvec{u}}}}({{{\varvec{X}}}},t)$. Using this deformation map, the deformation gradient can be computed as:

$$\begin{aligned} {{{\varvec{F}}}}= & {} {\text {Grad}}{{\varvec{\varphi }}}({{{\varvec{X}}}},t)={\text {Grad}}( {{{\varvec{X}}}}-{{{\varvec{u}}}}({{{\varvec{X}}}},t ) ) \nonumber \\= & {} 1+{{{\varvec{H}}}}({{{\varvec{X}}}},t) \end{aligned}$$

(3.1)

where ${{{\varvec{H}}}}({{{\varvec{X}}}},t)$ is the displacement gradient. The volume change J is given by the determinant of the deformation gradient: $J=\det {{{\varvec{F}}}}$ and the right Cauchy-Green-tensor is defined by: ${{{\varvec{C}}}}={{{\varvec{F}}}}^T {{{\varvec{F}}}}$. Based on these quantities, the following isotropic strain energy function $W({{\varvec{\varphi }}})$ is used in this benchmark

$$\begin{aligned} W({{\varvec{\varphi }}}) = \frac{\mu }{2} ( {{\text {tr}}}{{{\varvec{C}}}}-3-2\,\text {ln}\,J) + \frac{\lambda }{2} (J-1)^2 \end{aligned}$$

(3.2)

The material parameters used in this benchmark are given in Table 16. They lead to a nearly incompressible material with a Poisson ratio of $\nu = 0.4983$.

Table 16 Material parameters

3.3 Discretization and Results

An overall impression of the deformation is provided in Fig. 13 which depicts the displacement magnitude in the deformed configuration along with the wireframe of the corresponding hexahedral discretization in its initial configuration.

Fig. 13

Overall impression of results left: displacement magnitude on a fixed mesh of $4 \times 4 \times 4$ elements; right: von Mises stress on a $8\times 8\times 8$ mesh of $p=8$, geometrically graded towards the traction boundary

The convergence of the vertical displacement in z-direction at point P, $u_z(P)$, versus the degrees of freedom (DOF) of the discretization was investigated for the finite elements H1/EI9, H1, H2, H1/P0, O2/P1, H1/E9 (see e.g. Wriggers [9]) on regular meshes. Here $4 \times 4 \times 4, 8 \times 8 \times 8, 16 \times 16 \times 16$ and $32 \times 32 \times 32$ elements were used for the elements H1, H1/EI9, H1/P0 and H1/E9 with linear ansatz functions. For the hexahedral element H2 (27 nodes) and the tetrahedral element O2/P1with quadratic ansatz functions the number of elements was reduced such that the number of nodes was the same as for the elements with linear ansatz. The used elements are based on following ideas:

• H1: Classical hexahedral displacement based element with linear ansatz and $2\times 2\times 2$ nodes.

• H2: Classical hexahedral displacement based element with quadratic ansatz and $3\times 3\times 3$ nodes.

• H1/E9: Hexahedral mixed type ansatz with assumed enhanced strain fields and linear displacements. The element uses $2\times 2\times 2$ nodes and has 9 internal degrees of freedom related to the enhanced modes. These mode can be eliminated at element level using Schur complement techniques, details can be found in Simo and Armero [22].

• H1/EI9: Hexahedral mixed type ansatz with assumed enhanced strain fields and linear displacements, like H1/E9. However the enhanced strain modes are included here as a stabilization of a constant strain element, see Mueller-Hoeppe et al. [23].

• H1/TSCG: Hexahedral mixed type ansatz with 12 assumed enhanced strain fields and linear displacements. This element has a special stabilization to avoid hour-glassing of H1/E9, see Korelc et al. [24].

• H1/P0: Hexahedral mixed element with linear ansatz for the displacement field and discontinuous constant pressure field at element level.

• O2/P1: Tetrahedral mixed element with quadratic ansatz for the displacement field and continuous linear pressure field.

Table 17 $\left| u_z(P)\right|$ for given element types, Dirichlet BC with DOF elimination

The results are presented in Table 17 and Fig. 14. Since all these elements are developed using AceGen, see Korelc and Wriggers [19], they exhibit quadratic convergence within the Newton-Raphson solution algorithm. Both enhanced strain elements and the mixed O2/P1 and the H1/P0 elements are softer than the H2 element. All elements except the H1/E9 element converge practically to the same solution. For the H1/E9 element, solutions can only be obtained for the two coarsest meshes. For finer mesh resolutions the H1/E9 element depicts nonphysical hourglass instabilities.

Furthermore, a purely displacement based high-order finite element using the integrated Legendre basis functions described in Szabó et al. [25] (trunk space) was used and a p-extension was carried out on a fixed mesh with $4 \times 4 \times 4$ elements. The results are depicted in Fig. 14 and Table 18.

Fig. 14

Displacement at point P in z-direction ($u_z(P)$) for different element formulations and discretizations

It can be observed that the pure displacement formulation locks for $p=1$ on the h-refined meshes. This is expected for this nearly incompressible problem. However, locking may be controlled by increasing the order the finite elements as demonstrated by the curve given by a p-extension on a fixed mesh with $4\times 4\times 4$ elements.

Table 18 $\left| u_z(P)\right|$ for p-extension on a $4\times 4\times 4$ mesh, isotropic trunk space, Dirichlet conditions implemented with a penalty method

For the high-order finite element using the integrated Legendre basis functions, additionally the stress component $P_{zz}$ and the von Mises stress are measured along the diagonal line $A-P$ at $z=h$. The results on a discretization with $8\times 8\times 8$ pure displacement based elements of $p=10$ are given in Figs. 15 and 16. The expected stress is zero at the homogeneous part of the Neumann boundary until it jumps to the value of the load q under the load surface. The discretization is only able to capture this jump approximately even though the discretization was geometrically propagated with a factor of 0.15 towards the singularity for this particular case.

Fig. 15

Von Mises stress along diagonal from point P to point A on top face for the isotropic trunk space of order $p=10$ on a $8\times 8\times 8$ mesh, geometrically graded towards the traction boundary and for a $64\times 64\times 64$ mesh of TSCG elements

Fig. 16

$1 {\rm{st}}$ Piola-Kirchhoff stress component $P_{zz}$ along diagonal from point P to point A on top face for the isotropic trunk space of order $p=10$ on a $8\times 8\times 8$ mesh, geometrically graded towards the traction boundary and for a $64\times 64\times 64$ mesh of TSCG elements

4 Analysis of a Thin Plate

For the analysis of the complex thin geometries, such as shell and plate structures, various finite element methods have been developed to deliver accurate and efficient simulations. In this section, a study on the performance of three different element formulations is established for a thin plate. To this end, the thickness of the plate is varied from a thick to a very thin plate. Accordingly, two solid finite element methodologies as well as a solid-shell finite element technology are exploited.

4.1 Problem Definition

A square thin plate is to be investigated in this benchmark (see Korelc et al. [24]). It is clamped at all four sides (refer to Eq. 4.1) and is loaded by a distributed load of $q=0.0002$ MPa. Geometry, boundary conditions and loading of the plate in 3D are illustrated in Fig. 17. The thickness of the plate h is a variable parameter for later investigations, whereas the side lengths a of the square are constant and equal to 1 mm.

Fig. 17

Geometry, boundary conditions and loading of the thin plate

The corresponding boundary conditions read as follows:

$$\begin{aligned} \begin{aligned} x = 0 : u = v = w = 0 , \\ x = a : u = v = w = 0 , \\ y = 0 : u = v = w = 0 , \\ y = a : u = v = w = 0 , \end{aligned} \end{aligned}$$

(4.1)

where u, v and w are the displacements in x, y and z-directions, respectively.

4.2 Material Model

In the present work, a Neo-Hooke model for isotropic hyper-elastic material behavior is applied. The material parameters are $\Lambda = 144.2307692$ MPa and $\mu = 96.1538$ MPa (corresponding to $E=250$ MPa and $\nu =0.3$). The strain energy function W is given by

$$\begin{aligned} W= & {}\, \frac{\mu }{2} (\hbox {tr}\mathbf{C }-3) -\mu \, \hbox {ln}(\sqrt{\hbox {det}\mathbf{C }}) \nonumber \\&+\, \frac{\Lambda }{4} (\hbox {det} \mathbf{C }-1-2 \, \hbox {ln}(\sqrt{\hbox {det}\mathbf{C }}), \end{aligned}$$

(4.2)

where $\mathbf{C }=\mathbf{F }^T \mathbf{F }$ is the right Cauchy-Green-tensor with $\mathbf{F }$ denoting the deformation gradient (see [26]).

The same plate is examined for the case of linear elasticity with the thickness $h= 0.01$ mm ($a/h = 100$) in the work of Bayat et al. [27] and Bayat et al. [28]. In the aforementioned studies, different conforming and non-conforming element formulations are compared in terms of rate of convergence with respect to the mesh refinement level.

4.3 Finite Element Technologies

Three different finite element methods are applied in this work. The first solid element formulation is the standard eight-node tri-linear Q1 brick element with eight Gauss points. Furthermore, an eight-node tri-linear solid element (Q1SP) with reduced integration (one Gauss point in the center of the element) and hourglass stabilization technique is applied. This element formulation is well-known for its performance in overcoming the problems concerned with shear as well as volumetric locking (please refer to Reese and Wriggers [26] and Reese [29]). Finally, an eight-node solid-shell (Q1STs) finite element with reduced integration (one integration point) within the shell plane and at least two integration points over the thickness is used (Schwarze and Reese [30, 31]). Similar to the Q1SP element formulation, this element benefits from enhanced assumed stains (EAS). However, in addition to the Q1SP element, the transverse shear and curvature thickness locking are additionally addressed by the employment of the assumed natural strain (ANS) concept.

4.4 Convergence Study

In order to investigate the convergence rate of different element formulations, the vertical displacement $w_P$ at the point P on the upper side of the middle of the surface is investigated, see Fig. 17. To this end, different thicknesses h, different number of elements $n_e$ and $n_z$ in the planar directions (x and y) as well as in the out-of-plane direction (z) are considered. The ratio between the length and thickness of the plate (a/h) is varied by 20, 50, 100, 200, 500 and 1000 to cover a wide range of thick to thin plates.

Figure 18 shows the deformed plate with 10 times magnified displacements.

Fig. 18

a Displacement contour and b stress distribution of the deformed plate in z-direction (10 times magnified displacements). The geometrical aspect ratio $a/h=100$ and the mesh size $n_e \times n_z =16 \times 8$ are set

4.4.1 Influence of the Number of Elements in Thickness Direction

First, the influence of the number of elements in thickness direction on the displacement convergence rate is studied for the two finite element technologies Q1SP and Q1STs for the given aspect ratio $a/h=100$. Figure 19 shows the influence of the number of elements ($n_z$) in thickness direction on the deflection $w_P$ of the plate at the point P for different number of elements ($n_e$) in planar directions. In addition, the vertical displacements of the point P are given for the specific planar discretization $n_e=32$ in Table 19.

Fig. 19

Deflection $w_P$ at the central point of the plate for different number of elements in planar and thickness directions ($a/h=100$). a Q1SP element formulation, b Q1STs element formulation with 3 integration points in thickness direction

Table 19 Deflection $w_P$ at the central point of the plate ($a/h=100$) for different number of elements in the thickness direction with 32 elements in the planar directions ($n_e=32$)

According to Fig. 19, both solid and solid-shell element formulations converge with very few number of elements in thickness direction as long as enough number of elements in the planar directions are used. The difference between the performance of these two element formulations lies mainly in the application of only one element in the thickness direction. For a single element in z-direction, the Q1SP element formulation shows still a deviation from the converged solution. This is expected because one uses only one Gauss point in the thickness direction in this shell structure. The error of the Q1STs element technique for $n_z=1$ is small. This is due to the fact that the latter approach benefits from possessing two or more integration points in the thickness direction.

The hourglass stabilization can be chosen constant. In this way, the deformation of the element is not considered. Another option is to take the deformation of the element into account and change the hourglass stabilization based on the deformed element configuration. To maintain quadratic convergence in the Newton iterations, the stabilization must be kept constant within the load step. This action could bring about a dependency of the solution on the number of load steps (Reese et al. [32]). To examine this dependency, the convergence behavior of two element formulations with respect to the number of the load steps within a pseudo time range (1 sec) of a calculation is plotted in Fig. 20. The sensitivity of the results to the size of the load step is negligible as the estimated relative error ($|(w_p-w_{P,conv})/w_{P,conv}| \times 100$) lies around $1\%$ for having only five load steps. Nonetheless, both element formulations converge quadratically in each load step as expected.

Fig. 20

Estimated relative error of the point P with $a/h = 100$ against different number of load steps in one pseudo second

4.4.2 Influence of the Geometrical Aspect Ratio

The relative displacement $w_P/w_{P,conv}$ at the point P with respect to the converged solution $w_{P,conv}$ is investigated for thick as well as thin plates in Figs. 21 and 22, respectively. The x-axis of theses figures describes the mesh refinement in the planar directions as follows: $4 \times 4$, $8 \times 8$, $16 \times 16$ and $32 \times 32$, indicating the number of elements in x and y-direction, respectively.

Fig. 21

Deflection at the central point of the thick plate for different number of elements ($n_e$) in planar directions. a $a/h=100$, b $a/h=50$ and c $a/h=20$

Fig. 22

Deflection at the central point of the thin plate for different number of elements ($n_e$) in planar directions. a $a/h=200$, b $a/h=500$ and c $a/h=1000$

Figure 21 pictures the behavior of different element formulations for thick plates. In our examples, this implies the geometrical aspect ratio $a/h \ge 50$ with $a/h=100$ lying between the thick and thin plate regimes. The Q1STs element formulation overcomes the problem of transverse shear locking with very few number of elements (almost one) in z-direction due to the application of assumed natural strain (ANS) concept. On the other hand, The Q1SP element needs more elements in the thickness direction to correctly capture the shear deformations. Only few number of elements in the planar directions are used for these two elements since these elements eliminate shear locking effects. On the contrary, the Q1 element observes severe locking due to the presence of shear locking in spite of the refinement of the mesh.

For much thinner structures (here, $a/h \le 200$ illustrated in Fig. 22), the transverse shear deformations become negligible and thus transverse shear locking effects tend to disappear. In other words, for too high aspect ratios ($a/h = 1000$), the shell becomes a membrane in our example. This means that the need for the ANS concept in the solid-shell element formulation becomes less pronounced. As a result, both solid (Q1SP) and solid-shell (Q1STs) element formulations tend to perform similarly. Unlike the outstanding convergence rate of these two element formulations, the Q1 element shows still severe shear locking effects. However, by making the structures excessively thinner, the convergence rate of the Q1 element shows a slight increase. This can be due to the fact that by changing the geometries from shell to membrane, the bending stresses start to decrease.

Finally, the vertical displacement $w_P$ at the point P is plotted for different geometrical aspect ratios in Fig. 23. The size of the mesh is set to the coarsest level which is necessary for the Q1SP and Q1STs element formulations to converge. For these mesh refinement levels, the Q1 elements show severe locking effects.

Fig. 23

Vertical displacement $w_P$ at the point P with different thicknesses h

5 Material Interfaces

This collection of benchmarks pose the challenge of resolving material interfaces as well as high gradients or singularities induced by them. Especially discretization methods that avoid geometry conforming meshes, as presented e.g. in [33,34,35] will find interesting challenges in the presented two-dimensional benchmarks. Some extensions towards three dimensions can be found in [36].

In the following benchmarks, the convergence is measured by the relative error in the total energy with respect to the reference solution

$$\begin{aligned} \Vert e\Vert _{{\mathcal {U}}} = \sqrt{\frac{|{\mathcal {U}}_\text {ex}-{\mathcal {U}} |}{{\mathcal {U}}_\text {ex}} } \times 100\%, \end{aligned}$$

(5.1)

where ${\mathcal {U}}_\text {ex}$ and ${\mathcal {U}}$ denote the exact and approximated energy, respectively. For a linear elastic problem ${\mathcal {U}}=\frac{1}{2}\int _\Omega \varvec{\sigma }: \varvec{\varepsilon }\text {d}\Omega$ and for a heat conduction problem ${\mathcal {U}}=\frac{1}{2}\int _\Omega \kappa \nabla u \cdot \nabla u \text {d}\Omega$.

Two types of discretization methods, Isogeometric Analysis (IGA) [37] and p-FEM [38], are used to provide reference solutions for the energy ${\mathcal {U}}_\text {ex}$ on locally refined conforming meshes. Curved geometries are represented in IGA using NURBS or B-splines while the blending function method [38] is used in the p-version of the finite element method. The unknown field variable is approximated in IGA using Bézier extraction of truncated hierarchical B-splines [39] in combination with a safe refinement strategy [40]. In p-FEM geometrically graded meshes are used. In those cases where no analytic solution is available, improved energy values are obtained by the extrapolation technique described in [38].

In the following, all material and geometric parameters are given in consistent units.

5.1 Circular Inclusion

This benchmark is to test the representation of stress and strain fields across a material interface. The solution on the rest of the domain is smooth and optimal convergence rates should be achieved easily once the discontinuity introduced by the material interface is resolved properly. An analytical solution is available.

5.1.1 System

The two-dimensional plane strain problem is governed by the partial differential equations of linear elastostatics without body load. Given in the strong form of the boundary value problem it reads:

$$\begin{aligned} \nabla \cdot {\varvec{{\sigma }}}^{(i)}&= {\varvec{{0}}}&\forall {\varvec{{x}}}&\in \Omega ^{(i)},\quad i=1,\ldots ,n, \end{aligned}$$

(5.2a)

$$\begin{aligned} {\varvec{{\sigma }}}^{(i)}&= {\varvec{{{\mathbb {C}}}}}^{(i)}: {\varvec{{\varepsilon }}}({\varvec{{u}}}^{(i)})&\forall {\varvec{{x}}}&\in \Omega ^{(i)},\quad i=1,\ldots ,n, \end{aligned}$$

(5.2b)

$$\begin{aligned} {\varvec{{\varepsilon }}}({\varvec{{u}}}^{(i)})&= \frac{1}{2} \left( \nabla {\varvec{{u}}}^{(i)} +\nabla {{\varvec{{u}}}^{(i)}}^{\top }\right)&\forall {\varvec{{x}}}&\in \Omega ^{(i)},\quad i=1,\ldots ,n, \end{aligned}$$

(5.2c)

where $\Omega ^{(1)}$, $\Omega ^{(2)} \; \subset {\mathbb {R}}^{2}$ are defined in Fig. 24, ${\varvec{{\sigma }}}^{(i)}$ is the stress tensor, ${\varvec{{{\mathbb {C}}}}}^{(i)}$ the plane strain elastic material tensor, ${\varvec{{\varepsilon }}}$ the strain tensor, and ${\varvec{{u}}}^{(i)}$ is the displacement vector of the sub-domain $\Omega ^{(i)}$. The circular inclusion $\Omega ^{(1)}$ has material parameters $E_1$, $\nu _1$ and a radius a that is embedded in an isotropic and linear elastic ($E_2$, $\nu _2$) disc $\Omega ^{(2)}$ with radius b. A radial displacement is applied on the outer boundary as depicted in Fig. 24. The problem can be solved analytically under plane strain conditions leading to the displacement field in polar coordinates [41]

$$\begin{aligned} u_r&= {\left\{ \begin{array}{ll} \left[ \left( 1-\frac{b^2}{a^2} \right) \alpha + \frac{b^2}{a^2} \right] r &{} \text {for }\; 0 \le r \le a \\ \left( r- \frac{b^2}{r} \right) \alpha + \frac{b^2}{r} &{} \text {for }\; a \le r \le b \end{array}\right. } \end{aligned}$$

(5.3)

$$\begin{aligned} u_\theta&= 0 \end{aligned}$$

(5.4)

with

$$\begin{aligned} \alpha = \frac{(\lambda _1+\mu _1+\mu _2)b^2}{(\lambda _2+\mu _2)a^2+(\lambda _1+\mu _1)(b^2-a^2)+\mu _2 b^2} \end{aligned}$$

(5.5)

and the Lamé parameters $\lambda$ and $\mu$. The strains are given by $\epsilon _r = u_{r,r}$ and $\epsilon _\theta =\frac{u_r}{r}$ and the stresses by $\sigma _r=2\mu \epsilon _r + \lambda (\epsilon _r+\epsilon _\theta )$ and $\sigma _\theta =2\mu \epsilon _\theta + \lambda (\epsilon _r+\epsilon _\theta )$.

Fig. 24

Circular inclusion problem: model problem with domain and boundary condition as well as stress solution for $\sigma _{yy}$ is illustrated in marked part of the domain

To reduce the computational effort, the problem is only solved on a square with dimension c cut out of the cylinder. This setting eases the application of boundary conditions and thereby avoids additional error sources, especially if embedded domain methods with a regular background mesh are used. Homogeneous Dirichlet boundary conditions are applied at the lower and left edge of the square in y and x direction, respectively. Neumann boundary conditions are applied at the upper and right edge of the square using the analytical stresses.

The dimensions are given in Table 20 and the material parameters in Table 21. The analytically computed strain energy is

$$\begin{aligned} {\mathcal {U}}_{\text {ex}} = 1.239852433865801 \times 10^{6}. \end{aligned}$$

Table 20 Circular inclusion problem: dimensions of the domain in Fig. 24

Table 21 Circular inclusion problem: material parameters

5.1.2 Discretization and Results

An overall impression of the stress state $\sigma _{yy}$ is provided in Fig. 24. The jump of the stresses at the material interface is clearly visible.

The p-FEM discretization uses 5 elements. The circular shape at their boundary is represented exactly using the blending function method, c.f. Fig. 25. The solution is obtained by keeping the element mesh fixed and uniformly increasing the polynomial degree of the basis functions of each element (uniform p-refinement).

The IGA discretization is based on 7 NURBS-patches using at least a bi-quadratic basis that allows for an exact geometric representation of the circular inclusion, c.f. Fig. 25. The basis is $C^0$-continuous along the material interface. Due to the absence of singularities and the mild stress concentration, an uniform mesh refinement under constant polynomial degrees was carried out in the convergence studies.

Fig. 25

Circular inclusion problem: initial mesh for p-FEM (left) and IGA (right) discretization are illustrated. Please note that the NURBS mesh represents the circle exactly. Linear approximation results from post processing only

Fig. 26

Circular inclusion problem: convergence of total energy

Table 22 Circular inclusion problem: the strain energy ${\mathcal {U}}$ and corresponding degrees of freedom DOF are given for the p-FEM and IGA ($p=4$) computation

The results of the convergence study are shown in Table 22 and illustrated in Fig. 26. It is clearly visible that due to the $C^0$-continuous basis along the material interface and the absence of singularities, all methods are able to reproduce optimal convergence rates. For both methods, p-FEM and IGA with $p=4$, only 1000 degrees of freedom are needed to compute a strain energy that equals the analytically computed strain energy except for the last digit.

5.2 Elliptical Inclusion

The second benchmark is a linear elastic plane stress problem of an embedded elliptical inclusion. The high aspect ratio of the ellipse leads stress concentrations that have to be resolved.

Fig. 27

Elliptical inclusion problem: model problem with domain and boundary condition as well as stress solution for $\sigma _{yy}$ is illustrated

5.2.1 System

The model problem is again defined by Eq. (5.2), where ${\varvec{{{\mathbb {C}}}}}^{(i)}$ is the plane stress material tensor, $\Omega ^{(2)}$ is an elliptical inclusion with material parameters $E_1$, $\nu _1$ and aspect ratio $r_y/r_x$ that is embedded into a quadratic plate $\Omega ^{(1)}$ with material parameters $E_2$, $\nu _2$ and dimension L. The domain and boundary conditions are given in Fig. 27. The dimensions and material properties are depicted in Tables 23 and 24, respectively. The disc has a loose bearing on the lower and left side and a traction load $\hat{{\varvec{t}}}$ is applied to the right side of the disc.

Table 23 Elliptical inclusion problem: dimensions of the domain

Table 24 Elliptical inclusion problem: material parameters and traction load

5.2.2 Discretization and Results

To validate the extrapolated reference energy, the convergence of the solution is measured in the energy norm for p-refinement in a FE analysis and for h-refinement in an isogeometric analysis.

Fig. 28

Elliptical inclusion problem: initial mesh for P-FEM (left) and locally refined mesh for IGA (right) discretization are illustrated. In IGA the mesh is refined towards the stress concentration at the material interface

In the case of IGA, 12 patches are used to discretise the computational domain. The elliptical shape of the material interface is represented exactly using NURBS. To properly resolve the stress concentration, the mesh is refined using Bézier extraction of truncated hierarchical B-splines [39] as illustrated in Fig. 28 on the right. To allow for well graded meshes, a safe refinement strategy is used [40]. Starting from an initial mesh of 768 elements, seven refinement steps are performed for three different polynomial degrees. To also refine the vicinity of the stress concentration, the maximum number of hierarchical levels in a computation is restricted to four.

In the case of p-FEM, the domain is discretized by a mesh composed of 26 elements with blended edges conforming to the material interface. The mesh is spatially graded towards the stress concentration due to the high curvature of the interface at the major axis of the ellipse. The solution $\sigma _{yy}$ and the graded meshes are depicted in Figs. 27, and 28. To obtain a reference solution, the energy computed for a discretization with polynomial degrees $p=30$ and $p=47$ is extrapolated. The extrapolated value is

$$\begin{aligned} {\mathcal {U}}_\text {ex} =9.10131116644128\times 10^{-2}. \end{aligned}$$

Fig. 29

Elliptical inclusion problem: convergence in energy norm

Table 25 Elliptical inclusion problem: the strain energy ${\mathcal {U}}$ and corresponding degrees of freedom DOF are given for the p-FEM and IGA ($p=6$) computation

As shown in Table 25 and illustrated in Fig. 29, both discretization methods converge to the extrapolated reference energy.

5.3 Inclusion with a Corner

The third benchmark is a linear Poisson problem with a sharp inclusion adapted from [42]. The re-entrant corner at the material interface leads to a singularity in the solution.

Fig. 30

Inclusion with a corner: model problem with domain and boundary condition as well as heat flux magnitude is illustrated

5.3.1 System

This benchmark is governed by the Poisson equation:

$$\begin{aligned} \kappa ^{(i)}~\nabla ^2 \phi ^{(i)}&= -1&\forall {\varvec{x}}&\in \Omega ^{(i)} \end{aligned}$$

(5.6)

$$\begin{aligned} \phi&= 0&\forall {\varvec{x}}&\in \Gamma _D, \end{aligned}$$

(5.7)

where $\phi ^{(i)}$ denotes the temperature, $\kappa ^{(i)}$ the thermal diffusivity, $\Omega ^{(i)}$ and $\Gamma _D$ are defined as shown in Fig. 30. The dimensions and material properties are given in Tables 26 and 27, respectively. Here, the material interface $\Gamma _{12}$ has a sharp corner, inducing a vertex singularity. Moreover, the intersection of the material interface with the Dirichlet boundary $\Gamma _D$ introduces two additional weak singularities where the solution exhibits reduced continuity.

Table 26 Inclusion with a corner: dimensions of the domain in Fig. 30

Table 27 Inclusion with a corner: material parameters and traction load

The exact solution to this problem is given in polar coordinates $(r,\theta )$ by [43]:

$$\begin{aligned} \phi (r,\theta ) = A_1 r^{\lambda _1} h_1(\theta ) + A_2 r^{\lambda _2} h_2(\theta ) + {\mathcal {O}}(r^2), \end{aligned}$$

(5.8)

where $A_1$ and $A_2$ are scalar constants, $h_1(\theta )$ and $h_2(\theta )$ are smooth sinusoidal functions and

$$\begin{aligned} \lambda _1 = 0.731691779, \quad \lambda _2 = 1.268308221. \end{aligned}$$

(5.9)

5.3.2 Discretization and Results

The domain is discretized by 28 NURBS patches with at least quadratic polynomial degree to exactly represent the circular outer boundary of the domain, see Fig. 31. The initial mesh is locally refined towards the singularity at the re-entrant corner. Due to the local pre-refinement, uniform refinement was applied in the convergence study of the IGA discretization.

Fig. 31

Inclusion with a corner: the initial mesh for IGA (left) consists of 28 patches and is locally refined towards the singularity (zoom on the right). The NURBS mesh represents the outer circle exactly. The visible linear approximation results from post processing only

Fig. 32

Inclusion with a corner: the initial mesh for p-FEM (left) is locally refined towards the singularity (zoom on the right)

Fig. 33

Inclusion with a corner: convergence in energy norm

Table 28 Inclusion with a corner: the energy ${\mathcal {U}}$ and corresponding degrees of freedom DOF are given for the p-FEM and IGA ($p=5$) computation

The p-FEM discretization is depicted in Fig. 32. As in the previous examples, the mesh is conforming to material interface and the domain boundary is exactly represented by means of the blending function method. The mesh is graded geometrically towards the singularity located at the re-entrant corner, see Fig. 32. The problem is solved by keeping the mesh fixed and performing a uniform p-extension.

For this example, the extrapolated reference energy is given by

$$\begin{aligned} {\mathcal {U}}_\text {ex} = 1.016844315974886\times 10^{-1}. \end{aligned}$$

The convergence of the p-FEM and the IGA discretizations shown in Table 28 and illustrated in Fig. 33. In this case the grading of the p-FEM mesh towards the singularities considerably improves the accuracy per degree of freedom.

6 Phase Field Model for Brittle Fracture

6.1 Mathematical Model

Phase field models for fracture were introduced by the seminal works of Francfort and Marigo [44], Bourdin et al. [45], Bourdin [46]. Initially developed as a regularized version of the fracture mechanical treatment of cracks in the spirit of Griffith, recent extensions treat dynamic and inelastic fracture processes. In order to provide a robust setup for different numerical schemes and discretization methods, the proposed benchmarks are straightforward and disregard special problems concerning the distinction between tension and compression, as well as subtleties in the evolution equation. Details on these issues can for example be found in Kuhn [47].

With these preliminary remarks, we state the governing differential equations, which follow the original work of Bourdin [46]. Starting point is the free energy density of the system

$$\begin{aligned} \Psi (\varvec{\varepsilon },s)= & {} (s^2+\eta ) W(\varvec{\varepsilon }) \nonumber \\&+\,{\mathcal {G}}_{\mathrm{c}} \left( \frac{1}{4\epsilon } (1-s)^2+\epsilon |\nabla s|^2\right) \end{aligned}$$

(6.1)

which depends on the strain field $\varvec{\varepsilon }$ and the fracture field s. The meaning of s is the following: If $s=1$ the material is intact, while $s=0$ represents the fully broken material. The parameters ${\mathcal {G}}_{\mathrm{c}}$ and $\epsilon$ describe the cracking resistance (related to the fracture toughness) and the width of the regularization zone. The parameter $\eta$ models the residual stiffness of the fully broken material ($s=0$). The residual stiffness is required to ensure a non vanishing stiffness in a static analysis. In (6.1) $W(\varvec{\varepsilon })$ represents the elastic strain energy density. If an isotropic material behavior is assumed the strain energy density can be expressed with the help of the two Lamé constants $\lambda$ and $\mu$ by

$$\begin{aligned}&W(\varvec{\varepsilon } )=\frac{1}{2} \lambda \left( \text {tr} \, \varvec{\varepsilon } \right) ^2 + \mu \, \varvec{\varepsilon } : \varvec{\varepsilon } \quad \text {with} \nonumber \\&\varvec{\varepsilon } = \frac{1}{2} \left( \nabla {\varvec{u}} + \nabla ^\text {T} {\varvec{u}} \right) \end{aligned}$$

(6.2)

where $\varvec{\varepsilon }$ is the linearized strain, and ${\varvec{u}}$ represents the displacement field. The “ : ” in (6.2) indicates the scalar product between two second order tensors, and “$\text {tr}$” is the trace of a second order tensor. The free energy density (6.1) defines the stress $\varvec{\sigma }$ via

$$\begin{aligned} \varvec{\sigma } = \dfrac{\partial \Psi }{\partial \varvec{\varepsilon }} = (s^2+\eta ) \left( \lambda \, \text {tr} \, \varvec{\varepsilon } \, {\varvec{1}} + 2 \mu \, \varvec{\varepsilon } \right) \end{aligned}$$

(6.3)

where ${\varvec{1}}$ is the second order identity tensor. The stress has to satisfy the equilibrium condition

$$\begin{aligned} \text {div} \varvec{\sigma } = {\varvec{0}} , \end{aligned}$$

(6.4)

with “div” as the divergence. Thus volume forces are neglected. The static mechanical problem has to be supplied with proper boundary conditions, either Dirichlet conditions

$$\begin{aligned} {\varvec{u}} = {\varvec{u}}^* \quad \hbox {on} \, \partial {\mathcal {B}}_u \end{aligned}$$

(6.5)

with prescribed displacements ${\varvec{u}}^*$ on the displacement boundary ${\mathcal {B}}_u$, or Neumann conditions

$$\begin{aligned} \varvec{\sigma } {\varvec{n}} = {\varvec{t}}^* \quad \hbox {on} \, \partial {\mathcal {B}}_\sigma \end{aligned}$$

(6.6)

with prescribed traction ${\varvec{t}}^*$ on the traction boundary ${\mathcal {B}}_\sigma$. The vector ${\varvec{n}}$ denotes the outward unit normal. The fracture field s is governed by an nonlocal evolution equation, which is driven by the variational derivative of the free energy with respect to s. In the terminology of phase field models this evolution equation is also referred to as the time dependent Ginzburg-Landau equation:

$$\begin{aligned} \dot{s}= & {} -M\,\frac{\delta \Psi }{\delta s} \nonumber \\= & {} M\left[ 2{\mathcal {G}}_{\mathrm{c}}\epsilon \,\Delta s- 2s\, W(\varvec{\varepsilon }) + \frac{{\mathcal {G}}_{\mathrm{c}}}{2\epsilon }(1-s) \right] , \end{aligned}$$

(6.7)

where ${\dot{s}}$ is the rate of s with respect to time t, and “$\Delta$” denotes the Laplace operator. The mobility constant M has to chosen sufficiently large to approximate quasi-static crack growth conditions. The fracture field is also associated with either Dirichlet conditions

$$\begin{aligned} s= s^* \quad \hbox {on} \, \partial {\mathcal {B}}_s \end{aligned}$$

(6.8)

with prescribed fracture field $s^*$ on the Dirichlet boundary ${\mathcal {B}}_s$, or Neumann conditions

$$\begin{aligned} \partial s _{n} = \nabla s \cdot {\varvec{n}}= q^* \quad \hbox {on} \, \partial {\mathcal {B}}_q \end{aligned}$$

(6.9)

with prescribed flux $q^*$ on the flux boundary $\partial {\mathcal {B}}_q$. Frequently, homogeneous Neumann boundary conditions are set on the fracture field s and s is set to zero at certain locations to model initial cracks. Additionally, the fracture field $s({\varvec{x}},t)$ has to be supplied by initial conditions, as (6.7) is a first order differential equation in time t. Thus, we must provide:

$$\begin{aligned} s({\varvec{x}},0) = s_0 ({\varvec{x}}) . \end{aligned}$$

(6.10)

If the material is unbroken and no preexisting cracks are treated, $s_0 ({\varvec{x}}) = 1$ everywhere.

To conclude the presentation of the mathematical model some remarks are given:

• the fracture field s is bounded $s \in [0,1]$

• different evolution strategies are possible, see Kuhn [47]

  • damage-like behavior: ${\dot{s}} \le 0$

  • indicator-like behavior: ${\dot{s}}=0$ if $s=0$

• different solution strategies are possible

  • monolithic solution of (6.4) and (6.7)

  • staggered solution of (6.4) and (6.7)

  using the stress from (6.3) and strain energy density from (6.2) in both strategies.

6.2 3D Setup

In Fig. 34a 3D initial-boundary value problem is defined on a block-like geometry. The block is loaded on a part of its top surface by a linearly increasing displacement load, while being fixed on the lower surface.

Fig. 34

Geometry of the 3D benchmark problem

The following dimensionless parameters are used:

Geometry parameters:

$$\begin{aligned} a=2 , b=2 , c=1 \end{aligned}$$

(6.11)

Material and model parameters:

$$\begin{aligned}&\lambda = 100000 , \mu = 100000 , \epsilon = 0.1 ,\nonumber \\&\eta = 0.0001 , {\mathcal {G}}_{\mathrm{c}} = 1 , M = 10 \end{aligned}$$

(6.12)

Initial conditions at $t=0$:

$$\begin{aligned} s_0=1 \end{aligned}$$

(6.13)

Boundary conditions on the surfaces (A)-(E):

$$\begin{aligned} \begin{array}{ll} \text {(A)}: u = 0 , v=0 , w=0 , &{}\partial _n s = 0 \\ \text {(B)}: \sigma _{xx} = 0 , \sigma _{xy} = 0 , \sigma _{xz} = 0 , &{} \partial _n s = 0 \\ \text {(C)}: \sigma _{yy} = 0 , \sigma _{xy} = 0 , \sigma _{yz} = 0 , &{} \partial _n s = 0 \\ \text {(D)}: \sigma _{zz} = 0 , \sigma _{xz} = 0 , \sigma _{yz} = 0 , &{} \partial _n s = 0 \\ \text {(E)}: u = 0 , v = 0 , w = 0.01 t , &{} s = 1 \\ \end{array} \end{aligned}$$

(6.14)

This setup leads to a crack initiation in the bulk. The crack then propagates to the top surface resulting in a shell like fracture surface. Intermediate steps of the computation are reported in Fig. 35. The finest mesh of Table 29 is used for the plots in Fig. 35.

The results are obtained by standard 8-node finite elements with tri-linear shape functions. The same spatial interpolation is used for u, v, w and s. The time integration of (6.7) is done by an implicit/backward Euler scheme using adaptive time step control. The initial time step is set to $\Delta t = 0.01$. The time step $\Delta t$ is halved if the Newton method for the increments does not converge in 8 iterations. This time step halving is repeated at most 10 times. The time step $\Delta t$ is increased to twice its value if 4 or less Newton iterations are required. Using n elements along the edges a, b and c, the number of degrees of freedom (d.o.f.) is reported in Table 29. The monolithic solution strategy was used. As a quantitative result, the reaction force F in z-direction on the segment $\text {(E)}$

$$\begin{aligned} F = \int _{A(\text {E})} \sigma _{zz} \, \text {d}A \end{aligned}$$

(6.15)

is reported with respect to the prescribed displacement w in Fig. 36.

Fig. 35

Evolution of crack configuration and normal stress $\sigma _{zz}$ at four time steps ($t=0.2/0.367755/0.460529/0.92381$). Elements with $s<0.05$ are suppressed in the plots to indicate fracture

Table 29 Number of degrees of freedom for different spatial discretisations

Fig. 36

Displacement-reaction force curve for different spatial discretisations

6.3 2D Setup

In Fig. 37 a 2D initial-boundary value problem is defined on a block geometry. The block is loaded on a half of its top surface by a linearly increasing displacement load, while being fixed on the lower surface.

Fig. 37

Geometry of the 2D benchmark problem

The following dimensionless parameters are used:

Geometry parameters:

$$\begin{aligned} a=2 , b=1 \end{aligned}$$

(6.16)

Material and model parameters:

$$\begin{aligned}&\lambda = 100000 , \mu = 100000 , \epsilon = 0.1, \nonumber \\&\eta = 0.0001 , {\mathcal {G}}_{\mathrm{c}} = 1 , M = 10 \end{aligned}$$

(6.17)

Initial conditions at $t=0$:

$$\begin{aligned} s_0=1 \end{aligned}$$

(6.18)

Boundary conditions on the surfaces (A)-(D):

$$\begin{aligned} \begin{array}{lll} \text {(A)}: &{} u = 0 , v=0 , &{}\partial _n s = 0 \\ \text {(B)}: &{} \sigma _{xx} = 0 , \sigma _{xy} = 0 , &{} \partial _n s = 0 \\ \text {(C)}: &{} \sigma _{yy} = 0 , \sigma _{xy} = 0 , &{} \partial _n s = 0 \\ \text {(D)}: &{} u = 0 , v = 0.01 t , &{} s = 1 \\ \end{array} \end{aligned}$$

(6.19)

This setup leads to a crack initiation in the bulk. The crack then propagates to the top surface resulting in an arch like fracture line. Intermediate steps of the computation are reported in Fig. 38. The plots in Fig. 38 are computed by a discretization of $100\times 100$ elements. The results are obtained by standard 4-node finite elements with bi-linear shape functions. The same spatial interpolation is used for u, v and s. The time integration of (6.7) is done by an implicit/backward Euler scheme using adaptive time step control. The initial time step is set to $\Delta t = 0.01$. The time step $\Delta t$ is halved if the Newton method for the increments does not converge in 8 iterations. This time step halving is repeated at most 10 times. The time step $\Delta t$ is increased to twice its value if 4 or less Newton iterations are required.Using n elements along the edges a and b, the number of degrees of freedom (d.o.f.) is reported in Table 30. The monolithic solution strategy was used. As a quantitative result the overall reaction force F in y-direction on the segment $\text {(D)}$

$$\begin{aligned} F = \int _{A(\text {D})} \sigma _{yy} \, \text {d}A \end{aligned}$$

(6.20)

is reported with respect to the displacement v in Fig. 39.

Fig. 38

Evolution of crack configuration and normal stress $\sigma _{zz}$ at four time steps ($t=0.3 / 0.386697 / 0.426888 / 0.94779$). Elements with $s<0.05$ are suppressed in the plots to indicate fracture

Table 30 Number of degrees of freedom for different plane discretizartions

Fig. 39

Displacement-reaction force curve for different plane discretizations

7 Quasi-Static Pressure-Driven Cavity

7.1 Introduction

In this benchmark, we consider a lower-dimensional fracture in a d-dimensional domain. For the time being we restrict ourselves to $d=2$. This fracture has a constant length and varying width (also known in the literature as aperture, crack opening displacement or COD). The driving force is a given constant pressure prescribed in the fracture. The setting is motivated by the book of Sneddon and Lowengrub [48] and therefore known as ‘Sneddon’ benchmark or ‘pressure-driven cavity’. Analytical solutions are derived in Sneddon and Lowengrub [48] and are also discussed in [49]. Subsequently, [50, 51] coin the proposed benchmark, and provide numerical results.

7.2 Equations

In this section, we introduce the used notation and state the governing equations. In the following, the $L^2$ product is denoted as $(\cdot ,\cdot )$:

$$\begin{aligned} (x,y) := \int _{\Omega } x\cdot y \, d\Omega , \end{aligned}$$

and for vector-valued quantities by

$$\begin{aligned} (X,Y) := \int _{\Omega } X: Y \, d\Omega . \end{aligned}$$

The domain $\Omega \subset {\mathbb {R}}^d$ ($d=2$ for the 2D benchmark problem) is an open, connected and bounded set. The crack C is a one-dimensional set contained in $\Omega$.

7.2.1 The Energy Formulation

The Francfort-Marigo functional [44] describes the energy of a crack in an elastic medium as

$$\begin{aligned} E(u,C)=\frac{1}{2} (\sigma ,e(u))-\int _C \tau \cdot u\ \mathrm {d}s + G_C {\mathcal {H}}^1(C), \end{aligned}$$

(7.1)

where $u: \Omega \rightarrow {\mathbb {R}}^d$ is the vector-valued displacement function and $\sigma =\sigma (u)$ the classical stress tensor of linearized elasticity defined as

$$\begin{aligned} \sigma (u) :=2\mu e(u) +\lambda {{\text {tr}}}(e(u))I, \end{aligned}$$

with the Lamé parameters $\mu ,\lambda > 0$. The symmetric strain tensor e(u) is defined as

$$\begin{aligned} e(u):=\frac{1}{2} (\nabla u+\nabla u^T). \end{aligned}$$

The energy functional E consists of three terms: a bulk energy term, a traction energy, and a crack energy contribution. Specifically, traction forces are denoted by $\tau$. The critical energy release rate is denoted by $G_C>0$. The term ${\mathcal {H}}^1(C)$ stands for the Hausdorff measure denoting the length of the crack.

Combining this notation with the poro-elastic stress $\sigma _{\text {poro}}=\sigma - \alpha p I$ with the Biot coefficient $\alpha \in [0,1]$ and including the phase-field variable $\varphi$, we get the following energy functional Mikelić et al. [52]:

$$\begin{aligned} E(u,C)&= \frac{1}{2} (\sigma , e(u))-(\alpha -1)(\varphi ^2 p,{{\text {Div}}}u) \\&+ (\varphi ^2 \nabla p, u) + G_C {\mathcal {H}}^1(C), \end{aligned}$$

where $p:\Omega \rightarrow {\mathbb {R}}$ is a given pressure.

Remark 1

The considered benchmark setting in this study is posed in an elastic medium. For this reason, we set $\alpha =0$. Moreover, a constant pressure is applied, and therefore the pressure gradient is zero, i.e., $\nabla p= 0$. This yields the simplified energy functional

$$\begin{aligned} E(u,C)=\frac{1}{2} (\sigma , e(u))+(\varphi ^2 p,{{\text {Div}}}u) + G_C {\mathcal {H}}^1(C), \end{aligned}$$

which is used in the following.

Following Bourdin et al. [45], the crack C is approximated by a continuous phase-field variable $\varphi :\Omega \rightarrow [0,1]$, which is zero in the crack and one in the unbroken solid with a transition zone with size $\epsilon > 0$. The regularization parameter $\epsilon >0$ allows control of the diffusive transition zone. So the crack surface energy $G_C {\mathcal {H}}^1(C)$ in Eq. (7.1) can be reformulated to an elliptic Ambrosio-Tortorelli functional, which yields:

Formulation 1

(Regularized energy functional for pressurized fractures)

$$\begin{aligned} \begin{aligned} E_{\epsilon }(u,\varphi )=&\ \frac{1}{2}\left( \left( (1-\kappa )\varphi ^2+\kappa \right) \sigma (u),e(u)\right) +(\varphi ^2 p,{{\text {Div}}}u)\\&+ \ G_C\left( \frac{1}{2\epsilon }\Vert 1-\varphi \Vert ^2 + \frac{\epsilon }{2}\Vert \nabla \varphi \Vert ^2\right) , \end{aligned} \end{aligned}$$

where $\kappa$ is a positive regularization parameter for the elastic energy with $\kappa \ll \epsilon$.

Formulation 2

(Fracture energy minimization under constraint) Find functions u and $\varphi$ for almost all times t with

$$\begin{aligned} \text {min}\ E_{\epsilon }(u,\varphi )\quad \text {s. t.}\ \partial _t \varphi \le 0. \end{aligned}$$

(7.2)

The constraint realizes the crack irreversibility. Physically-speaking: the crack cannot heal. To derive an incremental version, the constraint is discretized in time via:

$$\begin{aligned} \frac{\varphi (t^{n+1})-\varphi (t^n)}{t^{n+1}-t^n} \le 0. \end{aligned}$$

Remark 2

We notice that this benchmark setting is a stationary test case. Due to the crack irreversibility constraint, we compute a few time step solutions in order to converge to this stationary limit.

7.2.2 The Weak Formulation

Our weak formulation is based on [51]. With

$$\begin{aligned}&V:=H_0^1(\Omega ), \\&W_{\text {in}}:=\{w \in H^1(\Omega )| w\le \varphi ^{n-1} \le 1\text { a.e. on }\Omega \}, \quad \text {and} \\&W:=H^1(\Omega ), \end{aligned}$$

we obtain the following weak formulation of Eq. (7.2).

Formulation 3

(Euler–Lagrange System of Formulation 2) Find $(u,\varphi ) \in V \times W$ with

$$\begin{aligned}&\left( \left( (1-\kappa )\varphi ^2 + \kappa \right) \sigma (u),e(w)\right) +(\varphi ^2 p,{{\text {Div}}}w)=0 \\&\quad \forall w \in V, \end{aligned}$$

and

$$\begin{aligned}&(1- \kappa )(\varphi \sigma (u):e(u),\psi -\varphi ) + 2(\varphi p {{\text {Div}}}u,\psi -\varphi )\\&\quad +\ G_C\left( -\frac{1}{\epsilon }(1-\varphi ,\psi -\varphi )+\epsilon \left( \nabla \varphi ,\nabla (\psi -\varphi )\right) \right) \ge 0 \\&\qquad \forall \psi \in W_{\text {in}} \cap L^{\infty }(\Omega ). \end{aligned}$$

7.3 A 2D Pressurized Fracture

The theoretical calculations of Sneddon [53] and Sneddon and Lowengrub [48] form the basis of this example. All settings to simulate the benchmark of the given pressurized fracture problem are listed in the following sections.

7.3.1 Setup

7.3.1.1 Domain Definition The two-dimensional domain $\Omega = (-10,10)^2$ is sketched in Fig. 40.

Fig. 40

Domain $\Omega$ (in 2D) with Dirichlet boundaries $\partial \Omega$, an initial crack C of length $2l_0$ and a zone of width $\epsilon$, where the phase-field function $\varphi$ is defined

7.3.1.2 Initial Conditions An initial crack with length $2l_0 = 2.0$ and thickness d of two cells on $\Omega _c=[-1,1] \times [-d, d] \subset \Omega$ is prescribed by help of the phase-field function $\varphi$, i.e., $\varphi = 0$ in $\Omega _c$ and $\varphi = 1$ in $\Omega \setminus \Omega _c$. Note that the thickness of 2d corresponds to $2h/\sqrt{2}$.

Fig. 41

Zoom-in to the center of the domain $\Omega$. The lower-dimensional crack with length $2l_0 = 2.0$ (middle line in black) is approximated as a volume by extending it by one cell in normal up and down direction

7.3.1.3 Boundary Conditions As boundary conditions, the displacements u are set to zero on $\partial \Omega$. For the phase-field variable, we use homogeneous Neumann conditions (traction free), i.e., $\epsilon \partial _n \varphi = 0$ on $\partial \Omega$.

7.3.1.4 Parameter Settings The determined parameter values are listed in Table 31.

Table 31 Setting of the benchmark and numerical parameters in 2D

If it is of interest to compute the exact values of the Lamé coefficients $\lambda$ and $\mu$, the converting formulae are: $\lambda = \frac{\nu E}{(1+\nu )(1-2\nu )}$ and $\mu = \frac{E}{2(1+\nu )}$.

7.3.2 Quantities of Interest

As benchmark comparisons, we propose the following quantities of interest:

• 1. Section 7.3.2.1: Crack opening displacement (COD) as a function of the x coordinate;

• 2. Section 7.3.2.2: The maximal $\text {COD}_{\text {max}}$ in the middle of the fracture;

• 3. Section 7.3.2.3: Total volume of the crack (TCV);

• 4. Section 7.3.2.4: Bulk energy;

• 5. Section 7.3.2.5: Crack energy.

For the COD, $\hbox {COD}_{\text {max}}$ and TCV, manufactured reference values can be computed for a infinite domain from the formulae presented in [48, Section 2.4]. We restate the corresponding formulae below and provide the reference values for our proposed material, and model parameters. In our comparisons, the reference values on an infinite domain are denoted by reference Sneddon. Additionally, we provide values based on fine computations on adaptively refined meshes for the finite domain. See Heister and Wick [54] for more details about the used method. These reference values are denoted by reference adaptive. Here, we also calculate the errors between the reference adaptive solution and our numerical results. These errors are ‘relative’ errors computed via

$$\begin{aligned} \frac{|{\texttt {numerical~result}} - {\texttt {reference~adaptive}}|}{|{\texttt {reference~adaptive}}|}. \end{aligned}$$

7.3.2.1 Crack Opening Displacement The width or crack opening displacement (COD) under spatial mesh refinement can be measured in the numerical tests. The width is defined as

$$\begin{aligned} \text {COD}(x):=[u \cdot n](x) \approx \int _{{-\infty }}^{{\infty }} u(x,y) \cdot \nabla \varphi (x,y)\ \mathrm {d}y, \end{aligned}$$

(7.3)

where $x\in [-1,1]$ the x-coordinate along the integral in y direction. The reference value for a infinite domain (cf. Sneddon and Lowengrub [48]) is given by:

$$\begin{aligned} \text {COD}_{\text {ref}} (x)&= 2 \frac{p l_0}{E^{\prime }} \left( 1-\frac{x^2}{l_0^2}\right) ^{1/2} \nonumber \\&= 1.92 \times 10^{-3} \left( 1-\frac{x^2}{l_0^2}\right) ^{1/2}, \end{aligned}$$

(7.4)

where p is the applied pressure, $l_0$ is the half crack length and $E^{\prime }:=\frac{E}{1-\nu ^2}$.

7.3.2.2 Maximum Crack Opening Displacement We evaluate the COD in the middle of the fracture (due to symmetry the maximum is attained at $x=0$):

$$\begin{aligned} \text {COD}_{\text {max}}:=[u \cdot n]({0}) \approx \int _{{-\infty }}^{{\infty }} u({0},y) \cdot \nabla \varphi ({0},y)\ \mathrm {d}y. \end{aligned}$$

7.3.2.3 Total Crack Volume The total crack volume (TCV) can be computed numerically using

$$\begin{aligned} \text {TCV}_{h,\epsilon }=\int _{\Omega } u(x,y) \cdot \nabla \varphi (x,y)\ \mathrm {d}{(x,y)}. \end{aligned}$$

(7.5)

A formula for the limit can be obtained using Sneddon and Lowengrub [48]. Using symmetry of the configuration, i.e., $u_y(x,0^+) = - u_y(x,0^-)$, and the known crack location $[-1,1]\times \{0\}$, one obtains

$$\begin{aligned} \text {TCV}_{2D}=2 \int _{-\infty }^{\infty } u_y (x,0^+)\ \mathrm {d}x, \end{aligned}$$

where $0^\pm$ denotes the respective limit from above or below, and $u_y$ denotes the second (y) component of the displacement.

Using the exact representation of $u_y$ (cf. Sneddon and Lowengrub [48], page 29)

$$\begin{aligned} u_y(x,0^+) = \frac{p l_0}{E^{\prime }} \left( 1-\frac{x^2}{l_0^2}\right) ^{1/2} \end{aligned}$$

we obtain:

$$\begin{aligned} \text {TCV}_{2D} = \int _{-\infty }^{\infty } 2 u_y (x,0^+)\ \mathrm {d}x = \frac{2\pi p l_0^2}{E^{\prime }}. \end{aligned}$$

(7.6)

Applied to our parameter settings, we consequently obtain the reference value for an infinite domain as:

$$\begin{aligned} \text {TCV}_{2D} \approx 6.03186\times 10^{-3}. \end{aligned}$$

In the following, the formula to compute bulk and crack energy as two further numerical quantities of interest are given.

7.3.2.4 Bulk Energy Next, we compute the bulk energy $E_B$ given by

$$\begin{aligned} E_B= \int _{\Omega } ((1-\kappa )\varphi ^2 + \kappa ) \psi (e)\, \mathrm {d}{(x,y)}. \end{aligned}$$

(7.7)

The strain energy functional $\psi (e)$ in Eq. (7.7) is defined as

$$\begin{aligned} \psi (e):= \mu {{\text {tr}}}(e(u)^2) + \frac{1}{2} \lambda {{\text {tr}}}(e(u))^2. \end{aligned}$$

Here, no manufactured reference values are provided and we only present values computed numerically.

7.3.2.5 Crack Energy Finally, we compute the crack energy

$$\begin{aligned} E_C = \frac{G_C}{2} \int _{\Omega } \left( \frac{(\varphi -1)^2}{\epsilon } + \epsilon |\nabla \varphi |^2\right) \, \mathrm {d}{(x,y)}. \end{aligned}$$

(7.8)

Again, no manufactured reference values are provided and we only present values computed numerically.

7.4 Numerical Results for the 2D Benchmark

7.4.1 Solution Approaches

We use $H^1$ conforming finite elements on quadrilaterals (2D). Specifically, we use bilinear elements $Q_1^c$ Ciarlet [55] for both, the displacements and the phase-field variable. In accordance with Heister et al. [51], a monolithic approach with an extrapolation of the phase-field variable in the displacement equation is used.

To treat the variational inequality, we employ a primal-dual active set method (again following [51]) with the recent open-source code published in Heister and Wick [54] based on deal.II Arndt et al. [56]. Furthermore, we emphasize that in this benchmark setting we exclusively use uniform mesh refinement.

The benchmark is sensitive to how the phase field variable is initialized from the given initial condition in 7.3.1.2. We use nodal interpolation¹ of the piecewise bi-linear finite element variable as shown in Fig. 41. Using an $L^2$ projection (or any other form of projection) instead, yields slightly different benchmark results.

The overall nonlinear discrete problem is then solved with a Newton method in which (when necessary), globalization is achieved with a simple backtracking line search method.

The code for the benchmarks performed here, again based on the Finite Element library deal.II Arndt et al. [56], is available at https://github.com/tjhei/cracks/tree/sneddon_spp_benchmark including all necessary parameter files.

7.4.2 Numerical Results

In the following sections all numerical results for the determined quantities of interest are given and compared to reference values.

We start the computation with a structured mesh of 10 by 10 cells, where each cell has an extent of 2 by 2 units and a diameter $h=\sqrt{8}\approx 2.828$. From there, we refine the mesh globally, cutting the size of each cell by a factor of 2 in each step.

Remark 3

Keep in mind that the reference values of COD and TCV in the previous section are valid for $\epsilon = 0$ and for an infinite domain $\Omega _\infty$. While we can expect convergence of the results presented here against the reference values with $\epsilon \rightarrow 0$, the problem is still given on a finite domain. Therefore, convergence against these reference values can not be expected.. Instead, we compare with with numerical results on the same finite domain but on very fine, adaptively refined meshes from Heister and Wick [54]. The reference values from Sneddon/Lowengrub are denoted as Reference Sneddon, while the numerical reference values are denoted as Reference adaptive.

7.4.2.1 COD (Quantity of Interest No. 1) Figure 42 shows the crack opening displacements (COD), computed with Eq. (7.3), over the width of the crack for several tests with different sizes of the mesh size parameter h. The computations correspond to 4 to 8 global refinements of the initial mesh. The dotted black line gives the exact COD values from Eq. (7.4).

Fig. 42

Benchmark 6: crack opening displacements (COD) for various refinement levels

7.4.2.2 Maximum COD (Quantity of Interest No. 2) To compare the COD values at one fixed point, in Table 32 the values of $\text {COD}_{\text {max}}$ using Eq. (7.3) are listed and the (relative) error of the numerical COD value in comparison to the computed value on an adaptively refined mesh (Reference adaptive). A zoom on the pressurized fracture is given in Fig. 43.

Fig. 43

Benchmark 6: a zoom of the left half of the pressurized fracture based on the setup in 2D with 7 refinements: 4,922,883 degrees of freedom, $h = 0.0220971$

Table 32 Benchmark 6: numerical $\text {COD}_{\text {max}}$ values for five refinement levels in 2D and respectively the error to reference values based on an adaptively refined mesh Heister and Wick [54]

7.4.2.3 Total Crack Volume (Quantity of Interest No. 3) In Table 33, the computed values of the total crack volume Eq. (7.5) are provided for five mesh refinement levels and compared to the exact TCV value Eq. (7.6).

Table 33 Benchmark 6: total crack volume for five refinement levels in 2D and respectively the error to reference values based on an adaptively refined mesh (Reference adaptive) Heister and Wick [54]

7.4.2.4 Bulk and Crack Energies (Quantity of Interest No. 4 and No. 5) For the bulk and crack energy no reference values are provided. Here, in Table 34, the bulk and crack energy values are shown for five mesh refinement levels using the definitions of Eqs. (7.7) and (7.8), respectively.

Table 34 Benchmark 6: bulk and crack energy for five refinement levels in 2D

8 Conclusion

In the first two benchmarks we provided results on hyperelastic and elasto-plastic problems at finite strains. Different (mixed) finite element technologies as well as p-FEM was used in order to obtain convergent results for certain displacement and stress values for the Cook’s membrane and an incompressible block problem. In benchmark problem in Sect. 4, the convergence behavior of a thin as well as a thick plate with a distributed load was investigated. Two different finite element technologies, namely solid and solid-shell were applied in addition to the standard Q1 finite element formulation. It was shown that due to the shell structure of the problem, the standard finite element method observes severe locking whereas the other methods show an outstanding convergence rate with respect to the mesh refinement. In the fourth benchmark, material interfaces were considered. Therein, two discretization methods, finite elements and isogeometric analysis are compared, which show very good agreements for various inclusion configurations. In the fifth benchmark, a phase-field fracture method was applied to two settings. First a three dimensional configuration and then a two-dimensional setting. As quantity of interest, displacement-reaction curves were adopted. In the last benchmark, a pressurized fracture in elasticity was considered for which numerical and manufactured reference values were generated. Moreover, some codes are open-source to reproduce the presented findings. All benchmarks shall serve as future reference for other discretization methods and other numerical solution algorithms.