The inspection of the stiffness matrix eigenvalues for a single element is important to assess the differences between formulations. It is now known that different formulations can produce exactly the same results, and new derivations are often clones of long-standing ones (cf. [23]). Based on experience, we use four distinct configurations (square, trapezoidal, parallelogram and irregular) and three distinct Poisson coefficient values (ν = 0, ν = 0.3 and ν = 0.49995). Besides our formulation, we also inspect Simo and Armero Q1E4 element (cf. [50]), which entails the cost of additional degrees-of-freedom. Figures 1↓ and 2↓ summarize the results.

Figure 1 Eigenvalue results for the tested elements (α − δ) and Simo and Armero Q1E4. α and δ are used in SIMPLAS. (Square and trapezoid elements).

Figure 2 Eigenvalue results for the tested elements (α − δ) and Simo and Armero Q1E4. α and δ are used in SIMPLAS. (Parallelogram and distorted elements).

A set of five linear verification tests is depicted in Figure 3↓. Besides the Patch-test from T.J.R. Hughes textbook [27] (which is a requirement of SIMPLAS), we use the linear version of Cook’s membrane (as described in Zienkiewicz and Taylor [58]) to assess the behavior under bending combined with the incompressibility condition (Poisson coefficients ν1 = 1 ⁄ 3 and ν2 = 0.499995 are tested). Distortion locking (see, e.g. [23] for this nomenclature) is assessed by the beam problems with two and five distorted elements as well as the MacNeal slender beam with both parallel and trapezoid distortions (cf. [35]). We implemented the element by Simo and Armero formulation (cf. [50]) for comparison. We compare with elements from the literature that satisfy the Patch-test: Abaqus [1] CPS4 , CPS4I, CPE4H and CPE4I elements, Pian-Sumihara hybrid [37], Simo-Rifai enhanced assumed strain element (EAS) [52], Piltner-Taylor consistent EAS [38]. In addition, the recent “area-coordinate” quadrilaterals (e.g. [16]), presenting often exact solutions for benchmark problems, fail (cf. [39]) the Patch-test (the so-called weak Patch test was devised) and show grossly wrong solutions (cf. [39]) in problems other than the traditional benchmarks. Results for the Patch test and Cook’s membrane are shown in Figure 4↓ from which the following conclusions can be drawn:

(c) Normalized results for the MacNeal beam for two load cases and three mesh types.

Figure 5 Results for of the linear verification tests. Best results are highlighted.

Concerning the examples with distorted meshes, the following conclusions can be drawn:

For the two-element beam, element δ shows superior results to elements in the literature. Element α presents poor accuracy, similar to Abaqus CPS4. Element β shows results very similar to the Pian-Sumihara and Piltner-Taylor elements.

For the five-element beam, element δ shows excellent results, and with exception of the exactly integrated Korelc-Wriggers element [29] (not shown here since that version is limited to linear elastic problems), we could not find better results.

For the MacNeal beam, element δ is better in 3 out of 6 cases, tied with the Pian-Sumihara element.

Concluding, we found the element δ the only appropriate candidate for the use in linear and nonlinear problems when bending is present. Element α can also be used in problems where bending behavior is not essential. These two elements are used in SIMPLAS. The α is the default for plane strain and δ is the default for plane stress and optional for plane strain.

The well-known upsetting problem shows hourglassing patterns when the Q1E4 (the element of Simo and Armero [50]) is employed, see the paper by Crisfield et al.[19]. In that paper, the presence of hourglass instabilities was detected, despite the use of the 5 − point integration rule proposed by Simo, Armero and Taylor [49]. Therefore, it is an excellent candidate for inspection of the δ element in terms of hourglass patterns. Schonauer et al.[47] managed to reduce the reported instabilities and also reported the excessive coarseness of Crisfield’s mesh to model the shear band zone. The problem is defined in Figure 6↓ where consistent units are used. Pressure and effective plastic strain contour plots are shown, reproducing a diffuse shear band region. A comparison with Abaqus 6.8 element CPE4H is made for the upsetting/load curve (element CPE4I failed to converge at an early stage). Both proposed elements show slightly lower loads and present no signs of the hourglassing behavior reported by Crisfield.

Element assessment in localization in a plane strain bar has been performed by Glaser and Armero [24]. A 200 element mesh is used, locally refined near the predictable necking zone (see Figure 7↓). This necking zone has a width 1.8% smaller than the pulling region. Taking advantage of symmetry, only one quarter of the bar is actually discretized. In Glaser and Armero, the context was the use of a special five-point quadrature aimed at attenuating the hourglass instabilities present in the original Q1E4 element. The test and results for elements α and δ are depicted in Figure 7↓. No hourglass was observed and results were compared with Abaqus 6.8 CPE4H and CPE4I elements. It is noticeable that the latter failed to converge. We can observe that our elements α and δ present somehow delayed localization when compared with Abaqus, but peak load is lower.

Two classical versions of the asymmetric block indentation problem inspired by Crisfield work (cf. [19]) are assessed here: one with quasi-incompressibility elasticity (Quasi-incompressible Neo-Hookean material) and another with J2 plasticity/relative Kirchhoff Saint Venant elasticity. The problem data is shown in Figure 8↓. Results for both versions are shown in Figure 9↓ where the displacement vector maps and contour plots are shown. Conclusions are drawn from this test:

No hourglass patterns are present in the results;

Both elements α and δ are immune to locking in the incompressible limit;

Results are less stiff than both elements CPE4H and CPE4I in Abaqus 6.8 (see also the convergence curve 10↓);

Element CPE4I fails to finish the run due to lack of convergence in the Newton iteration.

We use a nonlinear version of the Cook’s membrane shown in Figure 3↑ with J2 plasticity as discussed by Simo and Armero [50]. Other comparisons were made by César de Sá, Areias and Jorge (cf. [17]) where improved results were shown. We here are concerned with the convergence of the α and δ elements with the number of elements, as well as oscillations in the pressure distribution for the δ element. Figure 11↓ describes the problem and shows the contour plots. Neither hourglassing nor pressure oscillations are present. In terms of convergence of the tip displacement, the results by Simo and Armero are shown for comparison in Figure 12↓ where the excellent behavior of the δ element can be observed.

Plate bending tests for verifying the implementation in the linear case are detailed in Figure 13↓ (where the relevant data for the tests is shown) and in Figures 14↓ and 15↓ (comparisons with results from the literature). In Figure 14↓ comparisons are shown with the results presented in references [37, 12, 4, 33, 38, 54, 32, 44, 7] and the Abaqus S4 element. In all examples we use standard Gaussian quadrature with 2 × 2 Gauss points in the plane and 2 points along the thickness. Noting that although no enrichment of the displacement and rotation fields are performed, the present element is competitive with the best available formulations.

Linear shell tests are important for assessing competitiveness in both linear and nonlinear codes. Results in nonlinear problems are conditioned by the inherent linear formulation. However, due to the presence of interfering constitutive laws, stress integration and other interfering ingredients, results from nonlinear problems are seldom conclusive concerning the accuracy of the underlying interpolation approach. Therefore, we first focus on the five classical shell tests depicted in Figure 16↓. This Figure contains the relevant data. Results are compared with known element formulations in Figures 17↓ and 18↓. Test #1 is the well-known pinched cylinder problem (a reference value of 1.82488 × 10 − 5 m is used in the normalization, in agreement with Simo, Fox and Rifai [51]). For this test, results shown in Figure 17↓ attest the good performance of the present element. Test #2 is the Scordelis-Lo roof (see [35]) in which our element was found to have similar performance to the S4 element from Abaqus. Our implementation of Dvorkin and Bathe [20] MITC4 element results in slightly stiffer results in both tests #1 and #2. Test #3 is the closed pinched hemisphere [A] [A] The open 18○ variant does not provide additional information. with a target value of the radial displacement of 0.0924[51]. In this test, we found that our element converges faster than the S4 and MITC4 elements for coarse meshes. The test #4 is the so-called twisted beam and it has been adopted in many references with several variants. Our variant is the most demanding (with a thickness of 0.0032 m) and was proposed by Belytschko and Wong [13]. Two load cases are inspected, as indicated in Figure 16↓. Excellent results were obtained by the present formulation in case A. Test #5 is the so-called Raasch’s hook (cf. [48]) and involves bending and torsion. The target value adopted here is 5.02 consistent units (also employed in the first Author’s Ph.D. thesis and the Abaqus/Standard benchmark suite). Considered meshes are (elements × elements): 1 × 9, 3 × 18 , 5 × 36, 10 × 72, 20 × 144 and 40 × 288 elements. Results appear somehow worse than Abaqus S4 but the target value is well approximated by the present element. An inspection shows very good results for coarse meshes. The present results are better than the previous high-performance element HIS (cf. [7]) and QBM (cf. [10]). Overall, SIMPLAS shows one of the best performing elements four-node shell elements.

Two classical geometrically non-linear problems are solved, with the corresponding data summarized in Figure 19↑. Results for the clamped pinched cylinder (cf. [22]) are shown in Figures 20↓ and 21↓. These attest the robustness of the proposed formulation. Specifically, Newton-Raphson convergence is asymptotically quadratic (cf. Figure 21↓) since an exact linearization is performed [B] [B] We restrict the Newton step size to avoid iterating outside the convergence radius (5% of the geometrical dimensions for the displacement step and 0.05 Radians for the rotation degrees-of-freedom)..

Two finite strain plasticity problems are solved using our recent algorithm (cf. [9]): the elasto-plastic plate under deformation-dependent pressure and the pinched cylinder (both with plane stress von-Mises plasticity). Relevant data and deformed meshes (and contour plots) for these problems are shown in Figure 22↓. The plate under pressure was described by Hauptmann et al.[26] and further inspected in [7]. It consists of a plate with dimensions 508 × 508 × 2.54 consistent units. Elasticity modulus is E = 6.9 × 104 consistent units, Poisson coefficient ν = 0.3 and the yield stress is σy = 248 consistent units. One-fourth of the plate is analyzed with 32 × 32 elements with a sequence of deformed meshes and the effective plastic strain contour plot shown in Sub-Figure a↓. Results are presented in Figure 23↓. As for the other problem (pinched cylinder test with soft support), it is depicted in Sub-Figure b↓ and was also tested in [7]. It consists of a cylinder with 300 mm of radius, with E = 3000 MPa, ν = 0.3 and a linear hardening law: σy = 24.3 + 300εp MPa. Load-deflection results for the 16 × 16 are shown in Figure 23↓, are close to the ones by Wagner, Klinkel and Gruttmann [56] for high values of pinching displacement.

The Muscat-Fenech and Atkins plate (studied by these Authors and reported in [36]) was already tested numerically using the extended finite element method (XFEM) in Areias and Belytschko (see [5]). It was found that reasonable results are obtained combining J2 plasticity with a cohesive law. Only a slight deviation from the experiments was observed near the end of the analysis. Relevant data for this problem is shown in Figure 24↑. A pre-cracked steel plate is subject to two vertical forces F as indicated in that Figure. The crack advance criterion is the principal maximum strain (ε = 1.45714 × 10 − 3) corresponding to the maximum normal stress of 306 MPa. A fracture energy of JR = 255 N/m is used, the upper limit of what is reported in [36]. The Sutton criterion [34] is used for the crack path. Results are shown in Figure 25↑ where reasonable agreement can be observed.

A slender beam with a C profile made from ductile aluminum was clamped to our laboratory frame as shown in Figure 26↓. The loading strategy was very simple: 2 N, 2.5 N and 5 N disks were added to a hook hanged at a drilling in the beam and depth comparators were set above the beam. We employed a non-linear least-square code to determine the elasticity modulus. The hardening law was, in contrast, estimated by trial-and-error. Elasto-plastic buckling occurs at the flanges as shown in Figure 27↓. Comparative results for the proposed element are also shown in Figure 26↓.

Figure 27 Cantilever beam: comparative pictures of the elasto-plastic buckling of the flange and numerical prediction.

The predictive capability was found to be very good. The channel section cantilever beam by Eberlein and Wriggers [21] is reproduced here. Recently, Klinkel, Gruttmann and Wagner [28] tested this problem with a hybrid quadrilateral, although to a lower deformation than in the original paper. Relevant data is shown in Figure 28↓ where deformed and undeformed meshes are shown for a sequence of steps. We subdivide each quadrilateral element in [28] with four triangles. The load-deflection results are shown in Figure 29↓ and compared with reference [28]. Reasonable agreement with the reported results in [28] can be observed.

The first smooth plate is the ring test first proposed by Başar and Ding [11] to test formulations for finite rotations in shells. The relevant geometry, boundary conditions and elastic properties are depicted in Figure 30↓. We extend the ring up to 28 consistent units. Figure 30↓ depicts this extension. A comparison with the results by Sansour et al.[45] is also shown in Figure 31↓. We can observe that with our method we attain a higher value of ring extension and slightly more flexible results. Between the present element and a combination of DKT and Allman’s element, no differences were noted, in contrast with what was observed in the linear tests.

This benchmark has been used to advocate the use of mortar discretizations in contact problems, cf. Yang, Laursen and Meng [57]. We solve the ironing problem introduced by to support the use of mortar contact elements. Figure 32↓ shows the relevant geometric and constitutive data, agreeing with the references [57] and [25]. We compare the present approach with the results of these References in Figure 33↓. Differences exist in the magnitude of forces, and we concluded that this is due to the continuum finite element technology. We use finite strain B-Bar elements (cf. [53]), known to be more compliant than the standard Q4 isoparametric elements.

Figure 32 Ironing problem: relevant data and deformed mesh snapshots.

Figure 33 Ironing problem: results for the load in terms of pseudo-time, compared with the values reported by Yang et al.[57] and Hartmann et al.[25].

This classical elasto-plastic axisymmetric problem, relevant for crashworthiness studies was introduced by Laursen (cf. [31]). In his work, Laursen used the operator split method with the augmented Lagrangian method to solve the crushing of a cylinder. He used μ = 0 and μ = 0.2 to assess the implementation. We here reproduce this problem, and solve it with μ = 0, μ = 0.2 and μ = 1. We use the consistent elasto-plastic integration algorithm discussed in [9]. Results are shown in Figures 34↓ and 35↓ allow the following conclusions:

For μ = 0, our algorithm predicts the post-buckling behavior earlier than Laursen. In addition, higher imposed displacements are possible with our algorithm. We note that a slightly lower penetration in the rigid fixture when compared with Laursen.

For μ = 0.2, our algorithm produces very similar results to the ones reported by Laursen.

We also tested μ = 1 with visible differences with respect to the case μ = 0.2 in the fifth wrinkle.

We also show, for μ = 0.2, the error evolution as a function of iteration count in Figure 45↓.

Figure 34 Post buckling of an elasto-plastic cylinder. A picture corresponding to the 5 wrinkles (cf. 4 were reported by Laursen [31]) is also shown.

Figure 35 Post buckling of an elasto-plastic cylinder: compression/reaction results. A comparison with the results reported by Laursen [31] is performed.

These are other examples typically adopted within the mortar formulation contexts (cf. Yang et al.[57]). We here use these examples to assess the (quadratic) cone projection method for problems typically solved by mortar discretization techniques. Figure 37↓ shows the relevant data for the two tests. In test I, a Neo-Hookean constitutive law is used and in test II J2 plasticity is used with our recently proposed algorithm [9]. A sequence of results produced by the proposed algorithm is shown in Figure 38↓ for both tests. The effective plastic strain contour plot for test II is very similar to the one obtained by Yang et al.[57]. The evolution of reactions is also compared with published results in Figure 39↓. Good agreement with the mortar technique can be observed. We also show results for μ = 1 and μ = 2 and for these high values, slight oscillations appear.

We reproduce the problem of Puso and Laursen [41] who detected locking with node-to-segment approach. Relevant problem data is shown in Figure 41↓. Two stages of displacement of the deformable die (whose center is located at 2.5 units from the slab end) are used: in the first stage, t ∈ [0, 0.2] the die lowers into the slab a total of 1.4 units. In the second stage, t ∈ ]0.2, 1.5], the die moves 4 units in the longitudinal direction of the slab. Reaction results are shown in Figure 42↓ where good agreement with Reference [41] can be observed for the normal reactions. However, the tangential reactions are slightly lower in magnitude for our case.

Figure 41 3D ironing problem: relevant data and deformed meshes.

Figure 42 Reactions and comparison with the results of Puso and Laursen (cf. [41]).

The paper by Puso and Laursen [40] presents the concentric sphere problem indicating premature failure of the node-to-face algorithm (even for the frictionless case). Meshes are purposively incompatible to force gap discontinuities to occur during sliding of sphere faces. Figure 43↓ presents the relevant data for this problem (we explicitly represent the two obstacles). We solve this problem to a slightly lower imposed displacement and inspect the effect of friction in the interface between the two hollow spheres (contact with the obstacles has zero friction coefficient). Note that, in [40], no friction was present. Figure 44↓ shows the reaction results, compared with the values reported by Puso and Laursen [40]. The evolution of the energy error during Newton iteration is shown in Figure 45↓.

We now reproduce the benchmark proposed by Temizer [55] which introduces sharp corners and edges in 3D and finite strains (with the neo-Hookean constitutive model). We represent two compressed blocks using one-quarter of the geometry, as depicted in Figure 46↓. Temizer uses two distinct meshes (a coarse mesh with 5 × 5 × 3 elements for the lower block and a fine mesh with 8 × 8 × 4 elements). The present method requires a finer mesh so that severe spurious inter-penetration (edge/edge crossing) does not occur. Often, contact problems are apparently successfully solved with the node-to-face algorithm but in close inspection there are large edge crossing or wrong target face selections. Also, this fact is known to cause oscillations in the load response. Although there are ad-hoc remedies to this, we usually employ two sufficiently refined meshes. Our coarser mesh consists of:

Top block with 6 × 6 × 2 elements.

Bottom block with 9 × 9 × 5 elements.

A finer mesh with a bottom block containing a 12 × 12 × 8 arrangement is also adopted for comparison purposes. Although Temizer considered frictionless contact, we here adopt μ = 0.2 for application with our projection algorithm. The load/deflection results obtained are compared with the results of Temizer (which are scaled by 1 ⁄ 100) and shown in Figure 47↓, where a good agreement can be observed.

Figure 48 Verification with the three-point-bending test: relevant data and forces in the cohesive region.

In Figure 48↑, we can also observe two steps with the injected softening elements, for v = 0.6 mm and v = 1 mm. Geometry is 100 × magnified. For comparison, we use both numerical and experimental results provided, respectively, by Claudia Comi [18] and co-workers (numerical results) and also the experimental data gathered by Alfaiate and co-workers [2]. Figure 49↓ shows the effect of band thickness (in the non-dimensional form θ) in the load/displacement results. Subsequent Figures 50↓ and 51↓ show the effect of the step and mesh size, respectively.

This version of the single-edge-notched beam problem was introduced by Schlangen (cf. [46]) and is now modeled. It consists of a pre-notched beam loaded in two points and supported in other two points. Figure 52↓ contains the relevant data for this problem. It is appropriate for the assessment of the injection method since the experimentally observed result is a curved crack emanating from the right corner of the notch. Curved cracks are difficult to reproduce with smeared models. Since \bmP cannot be monotonous when crack propagation occurs, a control equation is used, which increases the shear mouth relative displacement of the initial notch (cf. Figure 52↓). This is called crack mouth sliding displacement (CMSD) and is imposed at the global level.

Figure 52 Schlangen’s SEN specimen: relevant data.

We test three triangular meshes containing 2587, 3729 and 5573 nodes. The resulting crack path is not far from the experimental envelope, as can be observed in Figure 53↓; even near the support the experimental observations are reasonably reproduced. A comparison with the experimental results and the DSDA method [2, 6], along with a study of mesh size influence is performed. As can be observed in Figure 55↓, after the peak load is reached, the numerical results are more brittle than the experimental results. According to Alfaiate, Wells and Sluys [3], this is due to the fact that an isotropic mode-I traction-jump law is used. Since crack paths are nearly insensitive to the value of θ, we fix the value θ = 0.2 for the mesh size effect study.

The bi-notched concrete beam proposed by Bocca et al.[15] is tested. The beam is simply supported in two points and subjected to two point loads. The effect of beam size in brittleness is assessed by using two different specimen dimensions. The corresponding experimental setting is described in detail in the original work [15]. From the set of specimens studied by Bocca et al. we retain the specimens, both with c ⁄ b = 0.8: one with b = 50 and another with b = 200 mm. These have reported experimental measurements in [15]. We are also concerned with the crack paths that were reported in [15]. Using the well-known cracking particle method, Rabczuk and Belytschko [42] obtained very good results for the crack path prediction, although with a slightly higher load than the experimental one. However, with the particle methods, there is an ambiguity in assigning the support dimension for the crack region. We use a single uniform initial mesh, with 7488 nodes and 14536 triangular elements. All relevant data is shown in Figure 56↓. For anti-symmetry reasons, we force the same mouth horizontal displacement at the edge of notches A and B: ΔuB = ΔuA. It has been debated if quasi-static simulations allow propagation of more than one crack (concerning this topic, see the excellent thesis by Chaves [110]), and the imposition of same relative displacement forces both cracks to evolve simultaneously. We obtain an excellent agreement with the experimental crack paths, as shown in Figure 57↓. The relatively wide spread of experimental crack paths is typical and is a consequence of the use of 6 specimens of reference [15]. Experimentally, some residual crack evolution in the opposite direction of the final path was observed and we also obtained that effect. Load-displacement results are shown in Figure 58↓ where a comparison with the measurements of Bocca et al.[15] and the cracking particle method of Rabczuk and Belytschko [42] is made. For the smaller specimen there is a somehow longer and lower curve than the observed one.

A tension test producing cup and cone fracture was described by Besson [14] and is here reproduced. We include explicit crack propagation by element injection subsequently to the satisfaction of f ≥ ff. Figure 59↓ shows relevant data for this problem. The stress/displacement results are compared with results obtained by Besson et al.[14]. We are also concerned with the cup and cone formation, which requires a sufficiently refined mesh and elongated elements in the radial direction. Figure 61↓ shows the crack formed in the necking region and propagating toward the outer surface, forming the cup and cone. This type of result, combining the localization in ductile materials and a physically meaningful ductile crack formation is very rare in the literature. In addition, robustness is an advantage when compared with GFEM/XFEM. The stress/displacement results are compared with results obtained by Besson et al.[14]. Crack path formation is compared between the two meshes in Figure 62↓. Some differences can be observed in the crack path, but these results are in-line with what is expected from crack propagation analysis with ductile materials.

The compact tension specimen described by Samal et al.[43] is reproduced. In that paper, the Rousselier yield function was adopted to model ductile fracture of a pre-cracked (a0 = 0.0161 m) specimen. Experimental results were also reported. The Authors used a gradient model to attenuate the mesh dependence (equivalent to the one by Areias et al.[8]). Relevant data for this test is shown in Figure 63↓. Here, the initial crack is explicitly represented (which was not in the original paper). Hardening law is inserted as a set of ordered pairs, as shown in Table 1↓. We monitor the force F as a function of the imposed displacement v and compare with the experimental results reported in Samal’s paper. This comparison is presented in Figure 64↓ where good agreement can be observed, despite the slightly higher values of reaction obtained here. Note that higher numerical values were also reported by Samal et al.[43]. A sequence of contour plots for the void fraction and effective plastic strain is shown in Figure 65↓. Very high deformations are possible without convergence problems.

Figure 63 Compact tension specimen. Lengths given in meters. The hardening function is provided in [43].

Table 1 Hardening law used in the compact tension specimen (Figure 63↑).

εp

y [MPa]

0

405

0.0568

569

0.0710

584

0.0994

621

0.1207

634

0.1917

675

0.2864

721

0.5089

796

0.8260

881

1.1030

935

Figure 64 Compact tension specimen: comparison between experimental (reported in [43]) and numerical results for three different initial meshes.