An enhanced corotational Virtual Element Method for large displacements in plane elasticity

Computational Mechanics (2024) 74:379–392 https://doi.org/10.1007/s00466-023-02437-1 ORIGINAL PAPER An enhanced corotational Virtual Element Method for large displacements in plane elasticity Marco Nale1 ·Cristina Gatta2 ·Daniela Addessi2 ·Elena Benvenuti1 ·Elio Sacco3 Received: 20 July 2023 / Accepted: 12 December 2023 / Published online: 13 January 2024 © The Author(s) 2024 Abstract An enhanced virtual element formulation for large displacement analyses is presented. Relying on the corotational approach, the nonlinear geometric effects are introduced by assuming nodal large displacements but small strains in the element. The element deformable behavior is analyzed with reference to the local system, corotating with the element during its motion. Then, the large displacement-induced nonlinearity is accounted for through the transformation matrices relating the local and global quantities. At the local level, the Virtual Element Method is adopted, proposing an enhanced procedure for strain interpolation within the element. The reliability of the proposed approach is explored through several benchmark tests by comparing the results with those evaluated by standard virtual elements, finite element formulations, and analytical solutions. The results prove that: (i) the corotational formulation can be efficiently used within the virtual element framework to account for geometric nonlinearity in the presence of large displacements and small strains; (ii) the adoption of enhanced polynomial approximation for the strain field in the virtual element avoids, in many cases, the need for ad-hoc stabilization procedures also in the nonlinear geometric framework. Keywords Virtual element method ·Corotational approach ·Stabilization-free 1 Introduction The structural analysis in the presence of large displacements can be carried out by different approaches, including more established Finite Element (FE) methods integrating the Total BMarco Nale [email protected] Cristina Gatta [email protected] Daniela Addessi [email protected] Elena Benvenuti elena.benv[email protected] Elio Sacco [email protected] 1Department of Engineering, University of Ferrara, via Saragat 1, Ferrara 44122, Italy 2Department of Structural and Geotechnical Engineering, Sapienza University of Rome, via Eudossiana 18, 00184 Rome, Italy 3Department of Structures for Engineering and Architecture, University of Naples Federico II, via Claudio 21, 80125 Naples, Italy and Updated Lagrangian formulation, and the Corotational (CR) approach. Pioneered by the seminal papers by Wempner [1], Belytschko and Hsieh [2], and Argyris and his collaborators [3], the key idea of the CR approach is to account for the effects of large displacements and rotations by decomposing the kinematics into a contribution provided by the large rigid-body displacements and the small-strain part. The CR approach has gained an ever-increasing interest in the computational mechanics community for it may reach results comparable to those obtainable through finite strain formulations. It offers indeed higher accuracy with respect to finite strain formulations and other approaches based on the small strain assumption and considering large displacements [4]. Another advantage is that the CR formulation obeys the principle of material frame indifference [5,6] and can deal very simply with material nonlinear laws [7]. CR formulations were developed for one, twoand threedimensional finite elements, such as beams [8,9], plates [10, 11], shells [12], and bricks [13,14], while unified CR formulations for general elements are also available [13,15]. The CR approach has been exploited in various fields of structural analysis, including applications to masonry walls both 123 380 Computational Mechanics (2024) 74:379–392 using beams [16,17] and bricks elements [14], shells made of shape-memory alloys [18], thin and moderately thick laminated composite structures [19], and soft biological tissues [20]. Variants are represented by the Implicit Corotational Method [21] to generate accurate geometrically nonlinear models, while Yaw and Sukumar [22] have extended the CR approach to the Maximum Entropy Meshfree Method [23]. However, the CR approach has been not yet extended to new generation generalized Galerkin finite element methods, such as the Mimetic Element Method [24], the Polygonal Finite Element Method [25–27], and the Virtual Element Method (VEM). Recently derived by Brezzi, Beirao, and coworkers [28] from the Mimetic Element Method, since its appearance, VEM has captured the general interest for it allows using general polygonal elements while completely releasing the meshing process from the requirement of generating regular undistorted meshes and avoiding hanging nodes. VEM was applied to a broad range of structural problems in linear [29], nonlinear [30], and finite elasticity and elastoplasticity [31], in elastodynamics, and fracture [32– 34]. Classical VEM requires stabilization, which can be carried out through one of the various stabilization techniques proposed in the literature, examples being the stabilization technique proposed by Beirão da Veiga et al. [28], Chi et al. [35], and the locking-free technique devised by Wriggers et al. [36]. As the unique identification of the stabilization parameters may be a hard task, while being occasionally problem-dependent, stabilization-free VEM formulations have been recently developed [37] which combine low-order displacement interpolations and high-order strain descriptions. Indeed, it is known that a meaningful advance in terms of computational efficiency and high accuracy in bendingdominated situations and even in the nearly-incompressible limit can be reached by utilizing enhanced formulations for low-order finite element interpolations. The relevant seminal enhanced formulations can be retraced back to the assumed strain-based approach [38,39], the B-bar methods [40–42], and the mode-decomposition and Hu–Washizu methods of Belytschko and co-workers [43]. Owing to their versatility, soon after their development, enhanced strain FE formulations were extended to geometrically non-linear problems [44], including a seminal CR-like approach [10]. Recently, it has been shown that the enhanced strain approach offers analogous advantages in the context of VEM [45], including the ability to deal with irregular meshes and improve accuracy, even for nearly incompressible materials. Of particular interest are the recent extensions of the enhanced strain approach to self-stabilized VE formulations. The key idea is to keep low the interpolation order even in the presence of many displacement degrees of freedom (DOFs) on the element boundary, as typical of polygonal VEs [45]. Another stabilization-free VE formulation, that does not require any additional internal DOFs and admits higherorder strain approximation, has been recently derived using Serendipity elements [46]. Similarly, a VEM based on the Hu–Washizu mixed variational principle has been recently proposed for 2D linear elastostatic framework [47]. The most notable advantage of this formulation is the possibility to prescribe the strain model a priori so that self-stabilized and locking-free low-order VEs with an accurate strain description can be obtained. In the present contribution, for the first time, the CR formulation is implemented within an enhanced, stabilizationfree, VEM framework for plane elasticity. In particular, the stabilization-free enhanced VEM developed by D’Altri et al. [45] is extended by exploiting a divergence-free polynomial approximation. The paper is organized as follows. Section2details the strong and weak formulations of the problem. Section 3is dedicated to presenting the CR approach. Section4describes the adopted enhanced virtual element formulation (EVE). Section5shows some numerical examples with a comprehensive discussion of the results. Finally, the salient conclusions and possible areas for further development are drawn in Sect.6. 2 Problem formulation The fundamental equations governing the two-dimensional (2D) continuum problem are described in the following. A 2D domain ⊂R2is considered with its boundary ∂ constituted by u∪t=∂ and u∩t=∅, where uand tare the portions where displacements uand tractions tare imposed, respectively. The problem is governed by the following strong form of compatibility and constitutive equations: ε=Du σ=Cε (1) where εand σare the strain and stress organized in vectors according to the Voigt notation, uis the displacement vector, Dis the compatibility operator, and Cthe material constitutive matrix. The elastic isotropic constitutive behavior is herein assumed. The above equations are complemented by the following Dirichlet and Neumann boundary conditions, that is: u=uon u(2) Nσ=ton t(3) where Ncontains the outward to unit normal vector. To derive the approximated VEM procedure, the weak form of 123 Computational Mechanics (2024) 74:379–392 381 the equilibrium equations is required, which is written by resorting to the virtual work principle as: Find u∈Usuch that  (Dδu)TCDu d= δuTbd +t δuTtd∀δu∈U0(4) with: U=u∈H1()2:u=uon u(5) U0=δu∈H1()2:δu=0on u(6) and band tdenoting the body force and traction vectors. For the considered 2D problem, u={uv}T,b= bxbyT,u={uv}T,t=txtyT,ε=εxεyγxyT, σ=σxσyτxyT. 3 Corotational formulation for the enhanced virtual element The 2D solid occupying the body is discretized into nonoverlapping polygons. Each polygon Eis characterized by the area Eand boundary ∂E, which is made up of nV vertices, nEedges and nnodes required for the displacement interpolation on the element boundary. Large displacement and small strain assumptions are considered to develop the enhanced VEM presented. The corotational formulation is introduced to this end by making reference to polygonal elements. In what follows, the procedure proposed by Battini et al. [11] for quadrilateral elements is extended to more general polygons. The overall motion of the polygonal virtual element Eis expressed as the composition of the rigid (Fig.1a) and deformation part (Fig.1b). With reference to the global coordinate system (O,X,Y)(Fig.1), the current location of the node iis denoted by Xi={XiYi}T, and that referred to the initial undeformed configuration by X0 i=X0 iY0 iT. Thus, the global displacement at node iis defined as: Ui=Xi−X0 i(7) and that of the element centroid Cresults: UC=XC−X0 C(8) The rigid motion can be expressed as the composition of the rigid translation, expressed by vector UC, and the rigid rotation θ(Fig.1a). The varied position of the element, after rigid body motion, is considering as the ‘corotational’ reference configuration. Referring to this latter, that is the local corotated system (C,x,y), the deformation displacement at node iis expressed as: ˜ Vi=xi−x0 i(9) where xiindicates the local position of node ias referred to the corotated reference system after the deformation process and x0 iits initial position in the same system. In Eq.9the current local position xican be evaluated as: xi=R(Xi−XC)(10) where Ris the matrix ruling the rigid rotation tranformation from the initial global reference system to the local one, expressed as: R=cos θsin θ −sin θcos θ(11) By introducing Eq.10 in 9, and using Eqs.7and 8,itfollows: ˜ Vi=RUi+X0 i−UC−X0 C−x0 i(12) The rigid rotation θis evaluated by minimizing the square of the Euclidean norm of the nodal deformation displacements and results: tan θ=n i=1x0 i(Yi−YC)−y0 i(Xi−XC) n i=1x0T i(Xi−XC)(13) The global and local deformation displacement vectors at all the nodes of each element are collected in vectors Uand ˜ V, respectively. The following relation links their variations: δ˜ V=BδU(14) where Bis the compatibility operator defined in the following. Considering the local deformation displacement at node i, it results: δ˜ Vi=R(δXi−δXC)+δR(Xi−XC)−δx0 i(15) and, after some manipulations, it is obtained: δ˜ Vi=R(δUi−δUC)−Yiδθ (16) with: Yi=−yi xi(17) 123 382 Computational Mechanics (2024) 74:379–392 Fig. 1 Kinematics of the virtual element: (a) rigid and (b) deformation parts The variation of the rigid rotation δθ is evaluated as: δθ = n  i=1 1 n i=1x0T ixi Y0T iR(δUi−δUC)(18) where: Y0 i=−y0 i x0 i(19) Introducing Eq.18 in Eq.16, it results: δ˜ V=(I−AG)ETδU(20) where matrix Eis the 2n×2noperator collecting the rotation matrices, i.e.: E= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ RT0... 0 0R T. . . . . ..... . . 0... ... RT ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ (21) Matrix A=YT 1YT 2...YT nTcollects vectors Yifor all the nodes located on the element boundary, and: G=1 n i=1x0T ixi A0(22) with A0=Y0T 1Y0T 2...Y0T n. The transformation matrix Bin Eq.14 is then defined as follows: B=(I−AG)ET(23) with Idenoting the identity matrix. The nodal forces workconjugate to displacement global and local vectors, Uand ˜ V, are collected in vectors Pand ˜ Pand are linked by the following relation: δP=BTδ˜ P(24) derived by invoking the equivalence of the element virtual work evaluated at the global and local level, respectively. The incremental relation between the global force and displacement vectors is: δP=∂P ∂UδU(25) 123 Computational Mechanics (2024) 74:379–392 383 Considering the incremental force-displacement relation at the element local level in the form: δ˜ P=˜ Kδ˜ V(26) and making use of Eqs.14 and 24,itfollows: ∂P ∂U=∂(BT˜ P) ∂U=BT∂˜ P ∂˜ V ∂˜ V ∂U+∂BT ∂U ˜ P=BT˜ KB+Kg (27) where the first term at the RHS indicates the material tangent matrix, while the second term is the geometric tangent matrix defined as: Kg=δE(I−AG)T˜ P+Eδ(I−AG)T˜ P(28) 4 Enhanced virtual element method The enhanced VEM formulation presented in [45] is here adopted and the main steps are recalled in the following. The virtual displacement field vat the generic point in the element domain E, according to the corotational approach, is the approximated deformation displacement as referred to the local corotated system (C,x,y), defined as (see Eq. 12): v=RU+X0−UC−X0 C−x0(29) According to the classical VEM assumptions, the displacement field vis not explicitly interpolated. Conversely, the approximated displacement field ˜ von the element boundary ∂Eis expressed as function of the displacement degrees of freedom defined at the nnodes lying on it by means of prescribed shape functions collected in matrix NV. The assumed interpolation of ˜ vis written as: ˜ v=NV˜ V(30) where vector ˜ Vhas 2nVkcomponents, with kdenoting the order of the polynomial interpolation adopted. To obtain a suitable representation of the approximated strain associated to the unknown virtual displacement, v, in the element interior, the projected strain εPis defined as the unique function that minimizes the energy norm, that is the solution of the following equation: E (δεP)TCεP−DvdE=0∀δεP∈Pq(E) (31) where Pq(E)indicates the polynomial space up to order q. Denoting with σPthe stress vector associated to the strain εP,Eq.31 can be expressed as follows: E (C−1δσP)TσP−CDvdE=0∀δσP∈Pq(E) (32) Then, the stress vector σPis interpolated on the basis of the parameters contained in vector ˆ σand the polynomial functions collected in the matrix ˇ N, that is: σP=ˇ Nˆ σ(33) By introducing the interpolation expressed by Eq.33 in Eq.32 and integrating by parts this latter term, the following expression results for the stress parameters: ˆ σ=E ˇ NTC−1ˇ NdE−1 ∂E (NT Eˇ N)TNVdE˜ V−E (DTˇ N)TvdE(34) where NEis the matrix containing the direction cosines of the outward unit normal vectors on E. Assuming the following definitions: G=E ˇ NTC−1ˇ NdE ˜ B˜ V=∂E (NT Eˇ N)TNVdE˜ V ˆ Bˆ V=−E (DTˇ N)TvdE (35) Equation34 can be written in compact form as follows: ˆ σ=G−1˜ B˜ V+ˆ Bˆ V(36) with vector ˆ Vcollecting the moments of the virtual displacement vrepresenting additional internal degrees of freedom. After the evaluation of the stress parameters by Eq.36,the projected strain εPis computed as: εP=C−1ˇ Nˆ σ(37) The LHS term in the virtual work equation in Eq. 4, expressing the internal work, can now be written in each element domain as a function of the projected strain εP, resulting as: E [εP(δv)]TCεP(v)dE =EG−1BδVTˇ NTC−1ˇ NG−1BVdE(38) 123 384 Computational Mechanics (2024) 74:379–392 having collected: B=˜ Bˆ B,V=˜ Vˆ VT (39) Equation38 can be manipulated as: δVTBTG−TE ˇ NTC−1ˇ NdEG−1BV=δVT˜ KcV (40) where the consistent definition of the element stiffness matrix ˜ Kcis introduced. Among the different possible procedures to select the polynomial representation for ˇ N,the divergence-free formulation is here adopted, corresponding to a self-equilibrated representation of σP. This is a special choice that leads to no internal degrees of freedom, as the last term at the RHS of Eq.34 vanishes being DTˇ N=0. Results reported in [45] proved that this assumption allows to obtain satisfactory solutions, especially if no distributed volume forces are considered. The slightly worse results obtained in the presence of volume forces could be improved by selecting other types of polynomial representation for ˇ N, which lead to more complex formulations involving also the internal degrees of freedom. In the presented enhanced formulation, differently from the standard VEM, the degree qof the polynomial representation in ˇ Nis not linked to the order kused in NV,but it is selected so that the number of modes considered in ˇ Nis greater than or equal to the element degrees of freedom purged by the number of rigid motions (i.e., modes in ˇ N≥nDOF −3 in the presented 2D formulation). Appendix A shows the matrix ˇ Nadopted in cases of q=1, q=2 and q=3. In such a way, self-stabilized elements are derived for which the stabilization term is not required and, consequently, the stiffness matrix ˜ Kis equal to ˜ Kcin Eq.40. To clarify, Table 1reports the number of modes in ˇ N and the polynomial degree qto be taken for the single enhanced virtual element with nVvertices and nDOF degrees of freedom, assuming k=1. For instance, if the element is characterized by nV=4 and nDOF =8, at least 5 parameters are needed to describe the strain field. Thus, the three strain components can be assumed linear (q=1) into the element leading to 9 strain (or stress) parameters, that reduce to 7 by enforcing the free-divergence. 5 Numerical examples This section aims to validate both the efficiency of the numerical procedure and the accuracy of the results obtained by employing the CR-EVE approach for a set of numerical examples in plane elasticity. These are properly selected from Table 1 Polynomial degree qto be adopted for the single element with nVvertices and nDOF degrees of freedom, and effective mode number of ˇ Nfor the divergence-free polynomial representation. (Table adapted from [45]) nVnDOF qˇ Nmodes 36 03 48 17 510 17 612 212 714 212 816 318 918 318 10 20 3 18 the literature for the sake of comparison with established formulations. Four classical examples of structures exhibiting geometric nonlinear response are considered, namely, thin and thick cantilevers, an L-bracket beam, and a circular shallow arch. For the EVE, the first-order (k=1) displacement approximation is adopted, while the degree qof the strain interpolation varies according to the number of nodes of the element. For comparison purposes, the VEM stabilized according to the procedure proposed by Artioli et al. [29]is also considered. Table 2summarizes the relevant information useful for identifying, for each example, the operational data concerning the CR-EVE and the adopted reference solutions. The degrees of freedom (DOFs), the average diameter havof the element in the mesh are reported. Furthermore, in the case of the CR-EVE formulation, the polynomial degree q of the strain approximation is detailed. Acronyms identifying the various formulations are used for brevity. In particular, CR-VE-Q4 denotes a stabilized four-noded VE while incorporating the CR approach. On the other hand, CR-EVE-Q4 indicates the homologous enhanced VE. Similarly, the formulations based on Voronoi meshes, generated using Lloyd’s iteration-based algorithm implemented in PolyMesher [48], are designated with the acronyms CR-VE-V and CR-EVE-V for the stabilized and enhanced self-stabilized VE formulation, respectively. From Table 2it clearly emerges that VEs denoted with CR-VE-Q4 and CR-VE-V are based on the use of polynomial approximation of the strain field with order q=0, while the CR-EVE-Q4 and CR-EVE-V elements are associated to q=1 and q=2, respectively. As for the VE simulations, the solution of the nonlinear problem is determined by means of an incremental-iterative Newton Raphson method using the Modified Generalized Displacement Control (MGDC) method [50] with a convergence tolerance of 10−6. The load–displacement curves obtained with the proposed approach are compared with both reference solutions provided by previously proposed CRFE formulations [11,22,49,51] and FE results obtained by Abaqus using overkill meshes made of 4-node CPS4 123 Computational Mechanics (2024) 74:379–392 385 Table 2 Data concerning the adopted VE and FE formulations in terms of: used acronyms, DOFs, average diameter of the element, stabilization instances, and polynomial degree qof the strain approximation Example ID Element DOFs havStabilization q Thin cantilever CR-VE-Q4 1280 0.0870 Yes 0 CR-EVE-Q4 1280 0.0870 No 1 CR-VE-V 1282 0.0920 Yes 0 CR-EVE-V 1282 0.0920 No 2 Hybrid Stress FEM [49] 126 n.a n.a n.a Abaqus-CSP4 21794 0.0182 n.a n.a Thick cantilever CR-VE-Q4 1298 0.2641 Yes 0 CR-EVE-Q4 1298 0.2641 No 1 CR-VE-V 1280 0.3842 Yes 0 CR-EVE-V 1280 0.3842 No 2 Meshfree [22] 738 n.a n.a n.a Abaqus-CPS4 41730 0.0442 n.a n.a L-Bracket CR-VE-Q4 772 0.3536 Yes 0 CR-EVE-Q4 772 0.3536 No 1 CR-VE-V 1218 0.3288 Yes 0 CR-EVE-V 1218 0.3288 No 2 Enhanced FEM [11] 772 n.a n.a n.a Abaqus-CPS4 10370 0.0884 n.a n.a Circular shallow arch CR-VE-Q4 1390 27.0931 Yes 0 CR-EVE-Q4 1390 27.0931 No 1 CR-VE-V 2204 25.1400 Yes 0 CR-EVE-V 2204 25.1400 No 2 Meshfree [22] 5522 n.a n.a n.a Abaqus-CPS4 67650 3.5083 n.a n.a elements [52], based on bi-linear shape functions for the displacement interpolation and plane stress assumption in finite elasticity. 5.1 Cantilever beams The first set of analyses investigates the structural response of cantilevers with two different slenderness ratios, corresponding to a thin and a thick cantilever subjected to uniform vertical load at the free end, illustrated in Sects.5.1.1 and 5.1.2, respectively. 5.1.1 Thin cantilever The response of a thin cantilever subjected to a shear load at the free end is analyzed assuming length l=10 mm, height h=0.1478 mm, and thickness t=0.1 mm. Elastic material with Young’s modulus E=108MPa and null Poisson’s ratio is employed. The same problem was studied by Karkon and Rezaiee-Pajand [49] through a hybrid-stress FE formulation. Figure2shows the undeformed beam configuration and the boundary and loading conditions. The adopted CR-EVE meshes consist of the same DOFs for quadrilateral and Voronoi elements, respectively. An Fig. 2 Thin cantilever: geometry (dimensions are in millimeters), boundary and loading conditions initial load factor λ =10.0 is set. Figure 3ashows the deformed geometry of the cantilever at the last load increment. The load–displacement curves, relating the vertical displacement of point A in Fig. 2to the total applied force, for the different formulations adopted are displayed in Fig.3b. The superior accuracy of the CR-EVE with respect to the classic stabilized VEM clearly emerges. It can be observed that the load–displacement curve obtained with the higher-order enhanced CR Voronoi elements (CR-EVEV) overlaps with the reference solution obtained using the Abaqus overkill mesh (CPS4) and the hybrid-stress formulation proposed by Karkon et al. [49]. Conversely, the load–displacement curve obtained with the CR-EVE-Q4 is less accurate with the same number of DOFs. Finally, to demonstrate the robustness of CR-EVE formulation, Table 3 123 386 Computational Mechanics (2024) 74:379–392 Fig. 3 Thin cantilever: adeformed configuration at the final step of the analysis and bload–displacement curves obtained with the standard and the CR-EVE compared with those derived from the hybrid stress FE formulation [49] and the Abaqus overkill mesh made of CPS4 FEs shows the vertical displacement δof the beam end section for variable DOFs with both quadrilateral and Voronoi meshes. The reference overkill-mesh result is also reported for the reader’s convenience. 5.1.2 Thick cantilever A thick cantilever, borrowed from Yaw et al. [22], is subjected to a shear load applied at the free end. The cantilever has a length l=10 in, height h=2 in, and thickness t=2 in and is made of a material with Young’s modulus E=100 ksi and null Poisson ratio. Figure4shows the cantilever geometry, boundary and loading conditions. The analysis is carried out with meshes of quadrilateral and polygonal virtual elements sharing the same DOFs. It has been verified that the use of 20000 quadrilateral CR-EVETable 3 Thin cantilever: vertical displacement δat P=200 N for variable DOFs ID Element DOFs havδ[mm] CR-VE-Q4 810 0.1303 5.794 CR-EVE-Q4 810 0.1303 7.243 CR-VE-V 812 0.1163 6.661 CR-EVE-V 812 0.1163 7.529 CR-VE-Q4 1280 0.0870 6.807 CR-EVE-Q4 1280 0.0870 7.534 CR-VE-V 1282 0.0920 7.137 CR-EVE-V 1282 0.0920 7.646 Abaqus-CPS4 21794 0.0182 7.699 Fig. 4 Thick cantilever: geometry (dimensions are in inches), boundary and loading conditions Q4 leads to an almost identical load–displacement profile. The computation of load–displacement curves, evaluated by monitoring the displacement of point A in Fig. 4, is performed with the MDGC with prescribed initial load factor λ = 0.2. Figure5a illustrates the deformed configuration of the thick cantilever, while Fig.5b displays the load–displacement curves for different VE meshes and enhancement types. The current results match well both the analytical solution obtained by Yaw [51] for an Eulero-Bernoulli-beam assuming large displacements and axial deformations and that computed through a meshfree CR-formulation [22]. For completeness, the Abaqus solution obtained with an overkill mesh of CPS4 elements incorporating a finite elasticity formulation is also reported. It can be observed that, at least in the present example, the incidence of non-infinitesimal strains emerges at the final stages of the loading history, as soon as the deformed cantilever tends to take an almost vertical configuration so that any further change of configuration should be ascribed to axial finite elastic deformations rather than to large displacements and rotations, a circumstance where Total or Updated Lagrangian finite elasticity formulations are more accurate [15]. The example confirms nevertheless that the presented CR-EVE approach leads to satisfactory results in terms of accuracy with respect to the homologous reference numerical solutions even when coarse meshes are employed. The influence of non-convex and distorted meshes on the robustness of the proposed CR-EVE formulation has been also tested. For this purpose, the meshes shown in Fig.6, made of distorted quadrilateral elements (CR-EVE-Q4D), 123 Computational Mechanics (2024) 74:379–392 387 Fig. 5 Thick cantilever beam: adeformed configuration at the end of the analysis, bload–displacement curves obtained with the standard and the enhanced CR-VEM compared with the available analytical Euler– Bernoulli beam solution [51] and results derived from the meshfree CR-formulation [22] and the Abaqus overkill mesh made of CPS4 FEs Voronoi Lloyd-iteration-based elements (CR-EVE-V), and non-convex elements (CR-EVE-NC) have been used. In Fig.6, the regular quadrilateral mesh (CR-EVE-Q4) is also displayed for the sake of comparison. Further details, in terms of number of DOFs, polynomial degree qand average diameter havof the elements of the meshes, are reported in Table 4. The resulting load–displacement curves in Fig.7confirm that the proposed formulation performs well even with nonconvex and distorted elements, as the results are not affected by the mesh type. Analogously to the thin beam example, the beam free end vertical displacement is shown in Table 5for all enhanced and classical formulations with both quadrilateral and Voronoi meshes with different DOFs number. The reference overkillmesh result is also shown. 5.2 L-bracket The third example is drawn from the paper of Battini [11] and involves the modeling of the bending behavior of the Lbracket shown in Fig. 8. The L-bracket is fixed at the top and it is subjected to a uniform shear load at the free end. Each branch has a length of 10 mm and a square cross-section. The elastic properties of the material are Young’s modulus E=3×107MPa, and Poisson ratio ν=0.3. Meshes of 304 square and Voronoi elements are used for the analysis. Figure9a reports the deformed shape at the last load increment, while the load–displacement curves for the various types of elements are reported in Fig. 9a. These are derived by relating the vertical displacement of point A at the free end (Fig.8) to the applied load. The obtained solutions are computed with a prescribed initial load factor λ =0.15. The results are consistent with the Abaqus overkill mesh and the numerical solution obtained by Battini [11]. Notably, in general, the self-stabilized CR-EVE formulation exhibits greater accuracy than the standard stabilized one, providing global load–displacement curves overlapped to those derived from the other reference numerical approaches. 5.3 Shallow arch The last studied structural example is illustrated in Fig.10 and concerns the thin circular shallow arch previously studied by Yaw et al. [22] by means of a CR-meshless formulation. The arch is simply supported with hinges located at both ends along the axis line. A concentrated force is applied on the symmetry axis at the upper surface. The arch geometry is featured by radius r=10581.6 mm, thickness t=79.2mm, depth w=25.4 mm, and rise f=76.48 mm, being the horizontal span length l=2540 mm. The material is modeled assuming Young’s modulus E=68.948 kN/mm2and null Poisson ratio. Meshes made of quadrilateral and Voronoi elements are assumed, adopting the initial load factor λ equal to 1.5. Figure11b shows the load–displacement curve of the considered shallow arch, which relates the force to the vertical displacement of its application point A. The snap-through behavior is satisfactorily captured, as shown by the comparison with the reference solutions derived from Abaqus with an overkill mesh of CPS4 elements and the CR meshless approach developed by Yaw et al. [22]. In particular, the results confirm that accuracy greatly improves when the enhanced VEM is employed instead of the classical VEM with stabilization. It is suggested that the slight discrepancies between the comparison and the present curves 123