A Selection of Benchmark Problems in Solid Mechanics and Applied Mathematics

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’.


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, computeraided 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.

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.

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 2 P 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 1 P 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].

Hyperelastic Material Model
In the following we use two polyconvex strain-energy functions. The first function reads where and constitute the Lamé parameters, J = det F the determinant of the deformation gradient F and I C = trC the first principal invariant of the right Cauchy-Green tensor C = F T F , see e.g. Ciarlet [5]. The 2nd Piola Kirchhoff stress tensor S = 2 C can be computed from the strain-energy function 1 as follows  here denotes the second-order identity tensor. The material parameters ( , ) as well as the load p 0 are given in Table 1.
The second strain-energy function under consideration is based on a split into an isochoric U(J) and unimodular W(IC, IIC) part see Hartmann and Neff [6] and Netz et al. [7]. The principal invariants of the unimodular right Cauchy-Green tensor are given with Following Hartmann and Neff [6] we choose the isochoric and unimodular parts of the strain energy as and The 2nd Piola Kirchhoff stress tensor for 2 can be computed as follows The material parameters used for the second strain-energy function ( , c 10 , c 01 , ) as well as the load p 0 are given in Table 2. Furthermore, a Fortran77 code for the computation of the 2nd Piola Kirchhoff stress tensor and the material tangent ℂ = 2 C S will be provided online under (https ://www.doi. org/10.5281/zenod o.40148 97 (see also Korelc and Stupkiewicz [14])) with the subroutines: Here, the input parameters are the deformation gradient (F as 3 × 3 field) as well as the material parameters. The output parameters are the 2nd Piola Kirchhoff stress tensor (S as 3 × 3 field) and the material tangent (Cmat as 3 × 3 × 3 × 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].

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 1 a mesh consisting of 13 hexahedral elements with two refinements towards the singularity is used and for 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 mm ) and reverse side ( z = −0.5 mm ) of the hexahedral meshes. In the following, we will consider both strain-energy functions 1 and 2 . An uniform and isotropic increase of the polynomial degree p = 1, 2, 3, … 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 2 P 0 for triangular elements. Furthermore, the results of the vertical displacement of point A are given in Tables 3 and 4. Figure 4 shows the contour plot of the xx stresses for the strain-energy function 1 . Here, the results of the T 2 and the T 2 P 0 element are depicted in Fig. 4a and b, respectively. Figure 4c shows the contour plot of the 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 2 is depicted. It shows the xx stresses for the T 2 and the T 2 P 0 element in Fig. 5a and b, respectively. To this end, the stress values are linearly     6 Cook's membrane: a convergence of displacement u y at point A for 1 b convergence of displacement u y at point A for 2 extrapolated from the integration points to the nodal points. Stress oscillations can be recognized in the case of the displacement based element. The 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.

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 1 and the hexahedral mesh with 15 elements for 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 2 P 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.
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 for the first free energy function 1 and the bulk modulus for the second free energy function 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 1 is computed as and for 2 as whereas according to Hartmann and Neff [6] the shear modulus is defined as 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 1 and 15 elements for 2 . The corresponding values are given in Tables 9 and 10.

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.

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 F is decomposed into an elastic F e and a plastic F p part by applying a multiplicative split 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 F e = FF −1 p , thus, the elastic right Cauchy-Green tensor C e can be formulated as and the elastic left Cauchy-Green tensor b e as where C p = F T p 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 In Eq. (2.16), and define the Lamé constants and I C e and III C e denote the first and third invariant of the elastic right Cauchy-Green tensor C e or the elastic left Cauchy-Green tensor b e as follows (2.14) The 1st and 2nd Piola-Kirchhoff as well as the Cauchy and Kirchhoff stress tensors read The yield criterion of the plasticity model reads In Eq. (2.19) ̄ denotes the hardening variable, s defines the deviatoric part of the Kirchhoff stress tensor , and K(̄) describes the hardening which consists of a linear and an exponential part as follows see e.g. Simo [12]. Here, Y 0 denotes the initial yield stress, H the linear hardening parameter, Y ∞ the saturation stress, and the hardening exponent. In accordance to Simo and Miehe [13], the evolution of the plastic flow as well as the internal variable reads Here, L v b e denotes the Lie derivative of the elastic left Cauchy-Green tensor b e and ≥ 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 ( Φ = 0 ), depending on the Dissipation inequality and the yield criterion in relation with . Furthermore, following Korelc and Stupkiewicz [14] Eq. (2.15) can be inserted into the plastic flow formula in Eq. (2.21) to obtain, (2.20) .       p an exponential map algorithm within an implicit time integration scheme is applied as follows which was introduced by Weber and Anand [15] and Eterovic and Bathe [16]. Here, Δ t 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.
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 isotropic increase of the polynomial degree p = 1, … , 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 2 P 0 , H 1 and H 1 P 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 = F P will be provided online under (https :// www.doi.org/10.5281/zenod o.40148 97 (see also Korelc and Stupkiewicz [14])) with the subroutine:   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 1 P 0 formulations is depicted in Fig. 8b, whereas in case of the T 2 and T 2 P 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 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 vM and Fig. 10b the equivalent plastic strain ̄ 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 SUBROUTINE mate(v,niter,MatIO,FIO,h1,SubIterationTolerance,P,Amat,h2,error) at point A and the von Mises stress 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.

Remark
The software used in this study are AceGen and AceFEM (Korelc [17,18], Korelc and Wriggers [19]) as well as inhouse 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.

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 (highorder) displacement-based and mixed finite elements.

Geometry and Boundary Conditions
The system is depicted in Fig. 12.
Its dimensions are given in Table 15.
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].

Hyperelastic Material Model
We consider a deformation (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: (X, t) = X − u(X, t) . Using this deformation map, the deformation gradient can be computed as: where H(X, t) is the displacement gradient. The volume change J is given by the determinant of the deformation gradient: J = det F and the right Cauchy-Green-tensor is defined by: C = F T F . Based on these quantities, the following isotropic strain energy function W( ) is used in this benchmark The material parameters used in this benchmark are given in Table 16. They lead to a nearly incompressible material with a Poisson ratio of = 0.4983.

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.
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 × 4 × 4, 8 × 8 × 8, 16 × 16 × 16 and 32 × 32 × 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:    with quadratic ansatz and 3 × 3 × 3 nodes. • H1/E9: Hexahedral mixed type ansatz with assumed enhanced strain fields and linear displacements. The element uses 2 × 2 × 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 hourglassing 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.
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 × 4 × 4 elements. The results are depicted in Fig. 14 and Table 18.
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 × 4 × 4 elements.
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 × 8 × 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.

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. The corresponding boundary conditions read as follows: where u, v and w are the displacements in x, y and z-directions, respectively.

Material Model
In the present work, a Neo-Hooke model for isotropic hyperelastic material behavior is applied. The material parameters are Λ = 144.2307692 MPa and = 96.1538 MPa (corresponding to E = 250 MPa and = 0.3 ). The strain energy function W is given by where = T is the right Cauchy-Green-tensor with 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.

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.

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.

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. 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 | × 100 ) lies around 1% for having only five load steps. Nonetheless, both element formulations converge quadratically in each load step as expected.

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   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 ≤ 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.

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 where U ex and U denote the exact and approximated energy, respectively. For a linear elastic problem U = 1 2 ∫ Ω ∶ dΩ and for a heat conduction problem U = 1 2 ∫ Ω ∇u ⋅ ∇udΩ. Two types of discretization methods, Isogeometric Analysis (IGA) [37] and p-FEM [38], are used to provide reference solutions for the energy U 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.

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.

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: where Ω (1) , Ω (2) ⊂ ℝ 2 are defined in Fig. 24, (i) is the stress tensor, ℂ (i) the plane strain elastic material tensor, the strain tensor, and u (i) is the displacement vector of the sub-domain Ω (i) . The circular inclusion Ω (1) has material parameters E 1 , 1 and a radius a that is embedded in an isotropic and linear elastic ( E 2 , 2 ) disc Ω (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] with and the Lamé parameters and . The strains are given by r = u r,r and = u r r and the stresses by r = 2 r + ( r + ) and = 2 + ( r + ). 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 (5.2c) (u (i) ) = 1 2 ∇u (i) + ∇u (i) ⊤ ∀x ∈ Ω (i) , i = 1, … , n,   The dimensions are given in Table 20 and the material parameters in Table 21. The analytically computed strain energy is

Discretization and Results
An overall impression of the stress state 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 U ex = 1.239852433865801 × 10 6 .  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.
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.

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.

System
The model problem is again defined by Eq. (5.2), where ℂ (i) is the plane stress material tensor, Ω (2) is an elliptical inclusion with material parameters E 1 , 1 and aspect ratio r y ∕r x that is embedded into a quadratic plate Ω (1) with material parameters E 2 , 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 t is applied to the right side of the disc.

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.
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   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 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 As shown in Table 25 and illustrated in Fig. 29, both discretization methods converge to the extrapolated reference energy.

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.

System
This benchmark is governed by the Poisson equation: (5.6) (i) ∇ 2 (i) = −1 ∀x ∈ Ω (i) Fig. 30 Inclusion with a corner: model problem with domain and boundary condition as well as heat flux magnitude is illustrated   Table 27 Inclusion with a corner: material parameters and traction load

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 where (i) denotes the temperature, (i) the thermal diffusivity, Ω (i) and Γ 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 Γ 12 has a sharp corner, inducing a vertex singularity. Moreover, the intersection of the material interface with the Dirichlet boundary Γ D introduces two additional weak singularities where the solution exhibits reduced continuity.
The exact solution to this problem is given in polar coordinates (r, ) by [43]: where A 1 and A 2 are scalar constants, h 1 ( ) and h 2 ( ) are smooth sinusoidal functions and

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.
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 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.

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 which depends on the strain field 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 G c and describe the cracking resistance (related to the fracture toughness) and the width of the regularization zone. The parameter 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( ) 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 and by U ex = 1.016844315974886 × 10 −1 .  where is the linearized strain, and u represents the displacement field. The " : " in (6.2) indicates the scalar product between two second order tensors, and " tr " is the trace of a second order tensor. The free energy density (6.1) defines the stress via where 1 is the second order identity tensor. The stress has to satisfy the equilibrium condition 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 with prescribed displacements u * on the displacement boundary B u , or Neumann conditions (6.3) = Ψ = (s 2 + )( tr 1 + 2 ) (6.4) div = 0, . 34 Geometry of the 3D benchmark problem  with prescribed traction t * on the traction boundary B . The vector 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: where ̇s is the rate of s with respect to time t, and " Δ " 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 with prescribed fracture field s * on the Dirichlet boundary B s , or Neumann conditions with prescribed flux q * on the flux boundary 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(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: 8) s = s * on B s (6.9) s n = ∇s ⋅ n = q * on B q (6.10) s(x, 0) = s 0 (x).
If the material is unbroken and no preexisting cracks are treated, s 0 (x) = 1 everywhere.
To conclude the presentation of the mathematical model some remarks are given: • the fracture field s is bounded s ∈ [0, 1] • different evolution strategies are possible, see Kuhn [47] -damage-like behavior: ̇s ≤ 0 -indicator-like behavior: ̇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.

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.
The following dimensionless parameters are used: Geometry parameters: Material and model parameters: 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 (6.14) 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 Δt = 0.01 . The time step Δ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 Δ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 (E) is reported with respect to the prescribed displacement w in Fig. 36.

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    Table 30. The monolithic solution strategy was used. As a quantitative result the overall reaction force F in y-direction on the segment (D) is reported with respect to the displacement v in Fig. 39.

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 yy dA 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 'pressuredriven 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.

Equations
In this section, we introduce the used notation and state the governing equations. In the following, the L 2 product is denoted as (⋅, ⋅): and for vector-valued quantities by is an open, connected and bounded set. The crack C is a onedimensional set contained in Ω.

The Energy Formulation
The Francfort-Marigo functional [44] describes the energy of a crack in an elastic medium as where u ∶ Ω → ℝ d is the vector-valued displacement function and = (u) the classical stress tensor of linearized elasticity defined as with the Lamé parameters , > 0 . The symmetric strain tensor e(u) is defined as 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 . The critical energy release rate is denoted by G C > 0 . The term H 1 (C) stands for the Hausdorff measure denoting the length of the crack. Combining this notation with the poro-elastic stress (u) ∶= 2 e(u) + tr(e(u))I, including the phase-field variable , we get the following energy functional Mikelić et al. [52]: where p ∶ Ω → ℝ is a given pressure.

Remark 1
The considered benchmark setting in this study is posed in an elastic medium. For this reason, we set = 0 . Moreover, a constant pressure is applied, and therefore the pressure gradient is zero, i.e., ∇p = 0 . This yields the simplified energy functional which is used in the following.
Following Bourdin et al. [45], the crack C is approximated by a continuous phase-field variable ∶ Ω → [0, 1] , which is zero in the crack and one in the unbroken solid with a transition zone with size > 0 . The regularization parameter > 0 allows control of the diffusive transition zone. So the crack surface energy G C 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) where is a positive regularization parameter for the elastic energy with ≪ .

Formulation 2 (Fracture energy minimization under constraint) Find functions u and for almost all times t with
The constraint realizes the crack irreversibility. Physicallyspeaking: the crack cannot heal. To derive an incremental version, the constraint is discretized in time via:

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.

The Weak Formulation
Our weak formulation is based on [51]. With we obtain the following weak formulation of Eq. (7.2).

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.  Table 31.
If it is of interest to compute the exact values of the Lamé coefficients and , the converting formulae are: = E (1+ )(1−2 ) and = E 2(1+ ) .

Quantities of Interest
As benchmark comparisons, we propose the following quantities of interest:    [48]) is given by:  where p is the applied pressure, l 0 is the half crack length and E � ∶= E 1− 2 .  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] × {0} , one obtains where 0 ± denotes the respective limit from above or below, and u y denotes the second (y) component of the displacement. y).
Using the exact representation of u y (cf. Sneddon and Lowengrub [48], page 29) we obtain: Applied to our parameter settings, we consequently obtain the reference value for an infinite domain as: In the following, the formula to compute bulk and crack energy as two further numerical quantities of interest are given. 7   (e) ∶= tr(e(u) 2 ) + 1 2 tr(e(u)) 2 .

Solution Approaches
We use H 1 conforming finite elements on quadrilaterals (2D). Specifically, we use bilinear elements Q c 1 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 primaldual 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 1 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 ://githu b.com/tjhei /crack s/tree/ snedd on_spp_bench mark including all necessary parameter files.

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 = √ 8 ≈ 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 = 0 and for an infinite domain Ω ∞ . While we can expect convergence of the results presented here against the reference values with → 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.  Table 32 the values of COD 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 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.

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.