Hybrid equilibrium element with high‐order stress fields for accurate elastic dynamic analysis

In the present article the two‐dimensional hybrid equilibrium element formulation is initially developed, with quadratic, cubic, and quartic stress fields, for static analysis of compressible and quasi‐incompressible elastic solids in the variational framework of the minimum complementary energy principle. Thereafter, the high‐order hybrid equilibrium formulation is developed for dynamic analysis of elastic solids in the variational framework of the Toupin principle, which is the complementary form of the Hamilton principle. The Newmark time integration scheme is introduced for discretization of the stress fields in the time domain and dynamic analysis of both the compressible solid and quasi‐incompressible ones. The hybrid equilibrium element formulation provides very accurate solutions with a high‐order stress field and the results of the static and dynamic analyses are compared with the solution of the classic displacement‐based quadratic formulation, showing the convergence of the two formulations to the exact solution and the very satisfying performance of the proposed formulation, especially for analysis of quasi‐incompressible elastic solids.


INTRODUCTION
The present article investigates the analysis of the elastodynamic problem using the stress-based approach of the hybrid equilibrium element formulation (HEE), which can be defined with high-order stress fields and can give great accuracy in stress computation. The equilibrium-based approach was introduced in the pioneering works of Fraeijs de Veubeke 1,2 for analysis of the elastic static problem and is defined in terms of stresses which implicitly satisfy domain and boundary equilibrium equations. The solution of the elastic static problem is determined in a weak form as a stationary condition of the complementary energy functional.
The equilibrium-based formulation was also developed in the hybrid form and several contributions [3][4][5][6][7][8] are available in the literature. In the hybrid equilibrium element (HEE) formulation, stress fields implicitly satisfy domain equilibrium equations and are independently defined for each finite element by polynomial stress functions of arbitrary order, with a set of so-called generalized stresses as degrees of freedom, which do not represent nodal values. The interelement equilibrium condition and the boundary equilibrium condition are imposed by the classic hybrid approach. The displacement The article is organized as follows: in Section 2 the HEE is developed for analysis of the elastic static problem with quadratic, cubic, and quartic stress fields; in Section 3 the proposed formulation is developed for analysis of the elastic dynamic problem in the time domain; in Section 4 some numerical simulations are presented the analysis of static and dynamic problems of both compressible and quasi-incompressible materials; Section 5 reports the conclusions and includes a discussion of future developments; finally, the appendix reports some details of the numerical formulation.

STATIC HYBRID EQUILIBRIUM FORMULATION
Let us consider a 2D elastic body occupying the closed region Ω. The body is referred to a Cartesian reference system (x, y) and is subjected to body force b (x, ) in Ω , traction t (x, ) on the free boundary Γ T , imposed displacement u ( ) on the constrained boundary Γ U and time ∈ (t 0 , t). The equilibrium formulation belongs to the class of stress-based approaches and the weak form solution of the elastostatic problem is given as the stationary condition of the complementary energy functional, with stress fields satisfying the domain and boundary equilibrium equations. The two-dimensional static equilibrium formulation considered in the present article follows the same reasoning path as proposed in References 10,18. The two-dimensional domain is discretized by a set of N e nonoverlapping triangular subdomains Ω e , with ⋃ N e e=1 Ω e = Ω. The element subdomain boundary Γ e = Ω e is composed of three sides Γ e s with s = 1, 2, 3, each of which can lie at the free boundary Γ e s ⊂ Γ T , or can lie at the constrained boundary Γ e s ⊂ Γ U , or can be an internal side between two subdomains Γ int ≡ Ωe 1 ∩ Ωe 2 , with e 1 ≠ e 2 .
The hybrid approach of the equilibrium formulation is developed with the stress fields e independently defined and satisfying the equilibrium equation in each subdomain (Div e + b e = 0 in Ω e ). The interelement equilibrium condition at all internal sides and the boundary equilibrium condition at all free boundary sides are imposed by the classical hybrid formulation, for which independent displacement fields u e s (x) are defined for each element side Γ e s , and they are assumed as Lagrangian variables in order to mutually connect adjacent elements or to apply traction on the free boundary (see, e.g., References 5,7,10,18,19). For the triangular finite element discretization, the hybrid equilibrium formulation gives the following modified complementary energy functional: (1) with u e s = u on Γ e s ∩ Γ U and n e s the outward unit vector normal to side Γ e s . The stationary condition of the functional Π c , with respect to the Lagrangian variable u s gives the weak form of the interelement equilibrium condition for an internal side and the weak form of the boundary equilibrium condition for a free boundary side. The stationary conditions of the functional Π c , with respect to the stress tensor e provides a weak form of the compatibility condition between elastic strains e = D ∶ e and displacement at the boundary sides u e s=1, 3 . Details of the static formulation are given in References 10,18. In hybrid equilibrium formulations the finite element is defined by the element stress fields satisfying the domain equilibrium equation, which does not interpolate nodal degrees of freedom, but is a function of generalized stresses. In the present article, the hybrid equilibrium element is developed only for two-dimensional membrane problems with polynomial stress fields of order n s = 2, 3, 4.
Let a triangular finite element of domain Ω e be considered and referred to a local Cartesian reference (x, y) centered at vertex 1, as shown in Figure 1A-C. The quadratic stress fields (n s = 2) of a two-dimensional element are defined by the following polynomial functions x = a 1 + a 2 y + a 3 y 2 − a 9 x − a 10 x 2 ∕2 − 2a 12 y = a 4 + a 5 x + a 6 x 2 − a 8 y − a 10 y 2 ∕2 − 2a 11 xy = a 7 + a 8 x + a 9 y + 2a 10 xy + a 11 x 2 + a 12 y 2 , where b x and b y are components of uniform volume force on the element, and terms a 1 , … , a 12 are generalized stress variables. The stress fields of Equations (2)-(4) implicitly satisfy the domain equilibrium equation and for the eth element can be represented in the following Voigt notation (C) (B) (A) F I G U R E 1 Nine-node, 12-node, and 15-node triangular hybrid equilibrium elements (HEE) where S e (x) is the coefficient matrix, a e (of dimension n a = 12) collects all generalized stress variables, and considered, but only with polynomial function of order lower than that of stress formulation. The cubic and quartic stress fields are analogously defined and the dimension of generalized stress vector a e is n a = 18 for the cubic stress and n a = 24 for the quartic one. See Appendix for the respective forms of the matrix S e (x) and vector a e in the quadratic, cubic, and quartic stress field formulations. Displacement is independently defined at each element side and continuity at the vertexes is not imposed. In the proposed stress-based approach, equilibrium conditions are satisfied, whereas displacement continuity can only be weakly imposed. The geometry and displacement components of the element sides Γ e s with s = 1, 2, 3 are modeled by a classic isoparametric mapping where n u is the number of nodes per side; m = n u (s − 1) + n is the node number; N n ( ) is the nth shape function (polynomial function of order n d = n u − 1) in the parental coordinate −1 ≤ ≤ 1; U s m and X s m respectively are the displacement and the Cartesian coordinate vectors of node m; N is the matrix collecting the shape functions and the vectors x e s , u e s collect respectively coordinates and kinematic degrees of freedom of side Γ e s . The complementary energy functional Π c of Equation (1) can be written in the discretized form as where with h e s = ∫ Γ e s S T e n e s NdΓ, with t e s = ∫ Γ e The degrees of freedom of the fixed boundaries Γ e s ⊂ Γ U are constrained, with assigned value u. The stationary condition of function Π c in Equation (8) with respect to the generalized stress vector a e gives the following element equation of the discretized static hybrid equilibrium formulation Π c a e = C e a e − H e u e = C e a e − h e 1 u e 1 − h e 2 u e 2 − h e 3 u e 3 = 0, (14) which states the relationship between nodal displacement and generalized stress variables at the element level. The element nodal force vector can also be written as q e = H T e a e and the equation of the single hybrid equilibrium element is where the compliance matrix C e is symmetric, positive definite and not singular, so that it can be inverted and the generalized stress a e can be condensed out at the element level, that is, mathematically where the matrix K e = H T e C −1 e H e is the element stiffness matrix and the HEE can be implemented in a classic displacement-based finite element code.
The use of the same order for the stress fields and for displacement of side (n s = n d ) allows the proposed formulation accurately to verify the interelement equilibrium condition, with co-diffusive stresses e 1 and e 2 through the interelement side Γ int = Ω e 1 ∩ Ω e 2 . Therefore, in order to impose the interelement equilibrium condition, the nine-node HEE in Figure 1A is employed with a quadratic stress field, while the 12-node HEE in Figure 1B is employed a with cubic stress field, and the 15-node HEE in Figure 1C is employed with a quartic stress field.
The proposed HEE formulation is suitable for analysis of quasi-incompressible materials, whereas it fails to analyze the pure incompressible elastic problem. According to Reference 31 (sec. 4.7), for Poisson ratio = 0.5 the compliance matrix C e is rank deficient by one, due to the uniform hydrostatic stress field, which gives null complementary elastic strain energy. As a consequence, the element stiffness matrix K e obtained by the static condensation in Equation (17) and depending on the inverse of compliance matrix, cannot be evaluated for incompressible elastic solids. In the same Reference 31 (sec. 4.7) an alternative strategy is proposed for the analysis of pure incompressible elastic problem, that is based on elimination of the uniform hydrostatic stress component in the stress basis and reduces the vector of generalized stress a e by one (i.e., n a = 11 for the quadratic stress field, n a = 17 for the cubic one, and n a = 23 for the quartic one). However, such approach can produces rank deficiency in the stiffness matrix and the solution uniqueness cannot be always guaranteed in terms of displacement. Therefore, the static and dynamic analyses of pure incompressible elastic problem are not tackled in the present article, but they represent an interesting topic for further development of the HEE formulation.

THE DYNAMIC EQUILIBRIUM FORMULATION
The HEE static formulation has been developed and implemented with high-order polynomial stress fields and this fine accuracy would also be of great interest for dynamic analysis of elastic bodies. The finite element formulation of the elastic-dynamic problem, in the framework of stress-based approaches, has been proposed by several authors [22][23][24][25][26][27][28][29] and it can be rigorously developed in the variational framework of Toupin's principle, 33 which is a complementary form of Hamilton's principle. 34 In Toupin's principle the velocity, and therefore the kinetic energy density, are defined as functions the time integral of stress, which is the areal density of the impulse J (t) = ∫ t t 0 d or, instead, the stresses are defined as time derivatives of the impulse as =̇J. Therefore the dynamic equilibrium equation can be written as follows: where is the mass density. After time integration and with null velocity at initial conditions (u (t 0 ) =u 0 = 0 ) the following relationship between impulse and velocity can be stated: bd ; the kinetic energy at time t can be defined as whereas the complementary potential energy is defined as a function of the time derivative of the impulse tensor in the following form: where the impulse implicitly verifies the free boundary equilibrium equation , n is the outward unit vector normal to the external surface Γ = Γ T ∪ Γ U and u (t) is the imposed displacement at the constrained surface Γ U .
Finally, the complementary form of the Hamilton functional is defined as and the variational approach of the dynamic response of the elastic solid is given in terms of the impulse tensor J (x, ) defined in the domain Ω × (t 0 , t), which implicitly verifies the free boundary equilibrium condition and for which the first variation of the functional in Equation (22) is null, that is The stationary condition in Equation (23) of the variational approach, after integration by part with respect to the time variable and using the divergence theorem, can be rewritten in the following form: where the symmetric gradient operator (∇ s ) emerges from the symmetry of tensor J and the conditions J (x, t 0 ) = J (x, t) = 0 in Ω and J (x, ) ⋅ n = 0 in Γ T × (t 0 , t) have been considered. Due to Equation (19), the first integral of Equation (24) provides the weak form of the kinematic boundary condition, that is and the second integral of Equation (24) provides the weak form of the strain-displacement kinematic conditions, that is with the strain rate given by the elastic constitutive equatioṅ= D ∶̇= D ∶J. Equations (18) and (24) with the elastic constitutive equation and under the hypothesis of the free boundary equilibrium condition implicitly satisfied, constitute the weak form of the partial differential governing equations of the dynamic formulation of the solid elastic problem.

The dynamic hybrid equilibrium element
The complementary form of the two-dimensional dynamic problem is developed in the present article by a hybrid finite element formulation similar to the static one. The domain discretization is considered with a set of N e nonoverlapping triangular subdomains Ω e , with ⋃ N e e=1 Ω e = Ω. The impulse field (with the relevant stress field) is independently defined in each subdomain as J e (x) in Ω e and it does not implicitly verify the free boundary equilibrium equation. The interelement equilibrium condition at all internal sides and the boundary equilibrium condition at all free boundary sides are imposed in the same way as in the quasi-static formulation. 4,18,35 Therefore, an independent displacement field u e s (x) is defined for each element side Γ e s , and it is assumed as a Lagrangian variable in order to mutually connect adjacent elements or to apply traction on the free boundary.
The hybrid dynamic equilibrium formulation with the assumed domain discretization gives the following hybrid form of the complementary Hamilton functional (or Toupin functional) (see References 22,23,31,33,36) with u e s = u on Γ e s ∩ Γ U and n e s the outward normal to side Γ e s . The stationary condition of this functional with respect to a single domain impulse field J e gives the same weak form in Equation (24), but limited to the subdomain Ω e , whereas the stationary condition with respect to the displacement field u s at the internal side between two adjacent subdomains Γ s = Ω e 1 ∩ Ω e 1 produces the following weak form of the interelement equilibrium condition: where n s = n e 1 s = −n e 2 s . Finally, the stationary condition of the functional in Equation (27) with respect to the displacement field of a free boundary sides Γ e s ∪ Γ T produces the following weak form of the free boundary equilibrium equation The hybrid dynamic equilibrium finite element formulation, based on Toupin's principle, is proposed in Reference 31 with the decomposition of the impulse field into the static and dynamic parts J e = J d e + J s e , where the static field implicitly verifies the domain equilibrium equation (div J s e + b e = 0 in Ω e ). The stress decomposition is also proposed in Reference 32 for analysis of the elastic-dynamic problem in the frequency domain and the estimation of the eigenfrequencies is improved by the combination of two dual formulations, that are the compatible DB finite element model and the equilibrated one. In References 31,32 the elastic-dynamic problem is analyzed in the frequency domain, where decomposition of the impulse field allows to remove the zero eigenvalues. In fact, the static component (J s e ) produces a null partition in the mobility matrix (defined in Equation (35)) with the relevant zero eigenvalues, whereas the dynamic component (J d e ) provides a full-rank mobility matrix.
In the present article the elastic-dynamic problem is analyzed in the time domain where the decomposition of the impulse field is an additional and not necessary computational effort. Therefore for each subdomain Ω e a quadratic polynomial law is assumed for all the components of the impulse tensor, which can be written in the following Voigt's notation: where a e ( ) = [a 1 ( ) , … , a 18 ( )] collects the generalized impulse time-dependent functions a i ( ) and S e collects the polynomial terms. Moreover, the stress field is defined as the time derivative of the impulse, and is simply defied as and the divergence of the impulse law, in Voigt notation, is defined as with T e = See Appendix for the respective forms of the matrices S e (x) and T e (x) in the quadratic, cubic, and quartic stress field formulations. In hybrid dynamic equilibrium formulations the finite element is defined by the element impulse and stress fields which do not interpolate nodal degrees of freedom, but are functions of generalized time-dependent variables a i ( ). In the present article, the proposed formulation is developed for two-dimensional membrane problems with polynomial impulse and stress fields of order n s = 2, 3, 4. The cubic and quartic stress fields are analogously defined and the dimension of generalized impulse vector a e is n a = 30 for the cubic stress and n a = 45 for the quartic one.
The independent displacement fields at the element sides are defined by the same isoparametric formulation as adopted for the static solution in Equations (6) and (7), and the stationary condition of the complementary Hamilton functional in Equation (27), with respect to the vector of impulse field a e ( ), and with integration by part with respect to the time variable , is written in the following discretized form: where is termed mobility matrix in References 22,31 and can be interpreted as an inverse mass matrix; is the compliance element matrix, which can be interpreted as an inverse stiffness matrix; the matrix H e is the stress-displacement coupling matrix, which is defined in Equation (10) but with the matrix of polynomial function S e (x) defined in Equation (30); finally is the vector of the element domain integral of the time-integrated body force b e . The stationary condition of the complementary Hamilton functional in Equation (27), with respect to the vector of element side displacement u e ( ), is written in the following discretized form where q e is the vector of element nodal forces. The dynamic hybrid equilibrium formulation, as well as the static one, maps impulse and stress fields as function of the Cartesian coordinate (x, y) and neither the spatial derivative nor the jacobian transformation are required in the formulation. As a consequence, the HEEs, in both static and dynamic formulations, are not pathologically sensitive to the element distortion.

Integration in the time domain
The time domain (t 0 , t) is discretized in a set of time intervals Δt = t n+1 − t n and for the first interval (n = 0) the initial conditions are known in terms of a generalized impulse vector and its derivatives a n e ,ȧ n e , andä n e , and in terms of nodal displacement and velocity u n e andu n e . The objective of dynamic analysis is to obtain, for a generic time interval, the approximation of the generalized impulse vector and nodal displacement with the relevant derivatives at the end of time step (a n+1 e ,ȧ n+1 e ,ä n+1 e , u n+1 e , andu n+1 e ) given the initial conditions at time t n . For a generic time interval the two integrand functions are assumed to be identically null at the initial step condition, that is The classic Newmark time integration procedure is considered for the generalized impulse vector a n+1 e = a n e +ȧ n e Δt + ( a n+1 e =ȧ n e + (1 − )ä n e Δt +ä n+1 e Δt, where and are the Newmark parameters and a n+1 e = a n e +ȧ n e Δt + where the vector f n+1 e is defined in Equation (37) and the trapezoidal integration rule has been applied for the time integration of the body force After substitution of the Newmark time integration rules defined in Equations (41) where Equation (40) is considered. The last two equations can be rewritten in the following element matrix notation where ( where is the element dynamic stiffness matrix and is the residual vector. Dynamic stiffness matrix and residual vector can be numerically evaluated with a simple static condensation procedure and can be implemented in a standard finite element code for the time-stepping analysis of the dynamic problem of elastic solids. The mobility matrix M e is rank deficient and, as pointed out for the static analysis, the proposed formulation is suitable for analysis of quasi-incompressible materials, whereas it fails to analyze the pure incompressible elastic problem. In fact, for Poisson ratio = 0.5 the compliance matrix C e and the sum of matrices (M e Δt 2 + C e ) are rank deficient, so the dynamic stiffness matrix K e obtained by the static condensation in Equation (52) cannot be evaluated for incompressible elastic solids. The proposed dynamic formulation is also applied for the modeling of time discontinuous load, whose solution in the equilibrium based formulation is governed by time discontinuous stress field. However, in the Newmark discretization of the time domain, the load discontinuity can be applied in a discretized form inside a single (small but finite) time increment Δt > 0. The results of a numerical simulation with time discontinuous load, and with two different values of the time increment, are reported in the article.
Finally, the proposed formulation is based on the assumption of null initial condition, with null displacements and null velocities at time t = t 0 (loading step n = 0), that are u e (t 0 ) = u 0 e = 0 andu e (t 0 ) =u 0 e = 0. In the case of a preloaded solid with initial elastic deformation associated to the nonzero initial conditions, the problem of the elastic static solid has to be preliminary solved and the relevant vectors of generalized stressȧ e (t 0 ) =ȧ (0) e and a 0 e =ȧ 0 e Δt have to be evaluated for each finite element, for resolution of Equation (49). However, this problem is not addressed in the present article.

NUMERICAL SIMULATION
The static and the dynamic formulations were implemented in the finite element code FEAP v8.5 37 using a triangular HEE with the quadratic, cubic, and quartic stress fields. In detail, the nine-node hybrid equilibrium element (HEE2) in Figure 1A is implemented with the quadratic stress field, while the 12-node element (HEE3) in Figure 1B is developed with the cubic stress field, and the 15-node element (HEE4) in Figure 1C with quartic stress field. The performance of the proposed formulation with high-order stress field is illustrated by some static and dynamic analyses and the results are compared with the solutions obtained with the well-known nine-node quadrilateral displacement-based element.
To the authors' best knowledge, the HEE formulation with triangular two-dimensional elements can be affected by SKMs on well-known patches of elements, as described in 35. These SKMs can be restrained by the approach proposed in Reference 18. The numerical simulations proposed in this article are based on meshes which can be affected only by the SKM on elements with two free boundary sides, that are the corners discretized by a single element. The SKM involves the two free boundary sides and can be restrained by a rigid constrain between the degrees of freedom of the two coincident nodes at the corners.

Cook's membrane elastic-static problem
Cook's membrane, represented in Figure 2, is a classic two-dimensional elastic problem, initially proposed in Reference 38, and analyzed by several authors for the performance evaluation of the finite element formulations, as proposed in References 18,28,39,40. The structure is constrained at the left side and loaded at the right side by a unitary force which can be applied by several different traction boundary conditions, although it is generally applied as a uniformly distributed tangential load. 18,28,39 While the trend of the tangential traction field at the free boundary side is not of great importance in a displacement-based formulation, it is pivotal in the equilibrium-based approach. At the two corners C 1 and C 2 represented in Figure 2, continuous and statically admissible stress fields (strictly enforcing the equilibrium equations) cannot be defined with the following boundary conditions: nonzero tangential stress n ≠ 0, null tangential stress m at the upper and lower sides, and null normal stresses m = n = 0. In fact, the two-dimensional stress tensor is defined by three independent components and only the null value n = 0 satisfy the equilibrium equations for a continuous stress field. The load with tangential stress n ≠ 0 can be applied in a strict equilibrium condition if the two corners are discretized with two or more elements, with discontinuous stress fields.
The problem of tangential traction distribution at the right boundary side can be addressed by three different solutions: F I G U R E 2 Sizes and geometry of the Cook membrane elastic problem • a uniform tangential load, for which the exact solution is not regular at the corner, can be applied with discontinuous stress fields at the two corners, by splitting the corners in two or more finite elements; • a quadratic tangential load with null values at the two corners, like the Jourawsky beam solution; • the real distribution of tangential load can be neglected by considering the right side to be rigid in its own plane. This approach also overcomes the problem of the SKM at the single element corner, being a kinematic constrain such as that proposed in Reference 18.
Although the Cook's membrane problem is usually solved with the uniform tangential load, the third approach is preferred and computationally applied by a rigid link between all the vertical degrees of freedom of the right side; the unit tangential force is applied to a node of the vertically rigid side.
The two-dimensional numerical analyses were performed under plane strain conditions with the elastic parameters: E = 1, = 1∕3 for a volumetric compressible test and E = 1, = 0.4999 for a quasi-incompressible test. Several discretizations were tested for both the HEE formulation and for displacement-based one and the details of the analyzed meshes, with the number of nodes, are shown in Table 1.
For the volumetric compressible test, the hybrid equilibrium approach and the displacement-based one converge to the exact solution, as shown in the convergence graphs plotted in Figure 3A,B. The vertical displacement at the right side (rigid) of the Cook membrane computed by the proposed hybrid equilibrium element, with quadratic (HEE2), cubic (HEE3) and quartic stress fields (HEE4) is compared in Figure 3A with the displacement obtained by the classic displacement-based nine-node element (Q9) and is compared with the converged solution (U y = 21.2265) obtained with a very fine mesh with Q9 elements. The relative error of the vertical displacement with respect to the converged solution versus the number of nodes is plotted in Figure 3B, showing the very good performance of the equilibrium-based formulation, especially for the coarse meshes. The maps of horizontal and vertical displacements obtained by the HEE and Q9 elements, with meshes of 8 × 12 elements, are plotted and compared in Figure 4. The nodes of the meshes are marked in the figure. In HEEs the displacement maps are plotted through a set of subdomain for each finite element and displacement is generally discontinuous at the vertices. Note that the HEE solutions are indistinguishable from one another and from the Q9 displacement-based response. The displacement discontinuities between adjacent sides are not appreciable in the maps of displacement, although the proposed formulation is based on a nonconforming element and the mesh is quite coarse.
To assess the accuracy of the present procedure the convergence analysis is also carried out in terms of the maximum principal stress in Figure 5A and in terms of relative error between the numerical solutions and the converge one ( 1 = 0.25439 obtained with a very fine mesh with Q9 elements) in Figure 5B. The results of the HEEs with quadratic, cubic, and quartic stress fields are compared with the Q9 solution and the converged value, obtained with the extremely refined meshes of the HEE4 and Q9 solutions. The percentage relative error shown in Figure 5B highlights the excellent performance of the high-order (cubic and quartic) stress field equilibrium formulation for the stress prediction in the elastic static analysis, showing a much lower error than the Q9 solution for a comparable number of nodes of the discretization. The elastic solution of the Cook's membrane produces a stress singularity at the upper left corner with negative principal stresses, therefore the convergence test cannot be represented in terms of the minimum principal stress or in terms of the maximum tangential stress, with enhancing values at the mesh refinement.  The elastic static response of the cook membrane problem is completed by the map of stresses, which are represented in Figure 6 in terms of the maximum principal stress, minimum principal stress, and maximum tangential stress, obtained with the HEE with quadratic stress field (HEE2), cubic stress field (HEE3), and quartic stress field (HEE4), and compared with the results obtained with the Q9 displacement-based formulation.
The accuracy of the proposed procedure is even more noticeable for the quasi-incompressible materials, for which the classic displacement-based formulation comes up against into the volumetric locking and the mixed formulations, or reduced numerical integration of the volumetric strain energy, are usually employed (see Reference 9). In terms of displacement the hybrid equilibrium approach and the displacement-based one converge to the exact solution, as shown in the convergence graphs plotted in Figure 7A,B, where the vertical displacement at the right side (rigid) of the Cook membrane, computed by the proposed hybrid equilibrium element, with quadratic (HEE2), cubic (HEE3), and quartic stress fields (HEE4) is compared with the displacement obtained by the classic DB nine-node element (Q9) and is compared with the converged solution. The converged displacement is evaluated as the mean value of the results obtained with the very fine meshes of the HEE4 and Q9 solutions. The convergence graphs in Figure 7A,B clearly show the better performance of the equilibrium formulation with respect to the displacement-based one, especially for coarse meshes. The accuracy of the HEE formulation is clear in the convergence diagram plotted in Figure 8 in terms of the maximum principal stress, showing an almost perfect convergence of the equilibrium based solution and the lack of convergence of the DB formulation. The maps of horizontal and vertical displacements obtained using the HEE and Q9 elements, with meshes of 8 × 12 elements, are plotted and compared in Figure 9. The HEE-solutions are almost indistinguishable from one another and from the Q9 displacement-based response and the displacement discontinuities between adjacent elements are not appreciable in Figure 9.
The convergence graphs in the Figures 3, 5, and 7 show also that the HEE formulation is a potential tool for the dual analysis and the error estimation of DB solutions. In fact, as stated in the pioneering work of Fraeijs de Veubeke 1 and more recently in References 41,42, the displacement based formulation and the equilibrium based one respectively provide a lower bound and an upper bound in the elastic solution in terms of strain energy. However, bounds of arbitrary nodal displacements or of stresses cannot be guaranteed.
The elastic static response of the cook membrane problem is completed by the map of stresses, which are represented in Figure 10 in terms of the maximum principal stress, minimum principal stress, and maximum tangential stress, obtained with the HEE with quadratic stress field (HEE2), cubic stress field (HEE3), and quartic stress field (HEE4), and with the Q9 DB formulation. The maps of principal stresses in Figure 10 make it evident the limitation of the classic DB formulation for analysis of nearly incompressible elastic material, with inconsistent stress distributions and with nonconverging solution in terms of stress. Such a problem affects the volumetric stress, therefore and the shear stress is less influenced. The HEE formulation is completely free of any volumetric locking numerical problem for all the considered stress orders.

F I G U R E 6
Maps of principal stresses and maximum tangential stress of compressible Cook membrane test, for the hybrid equilibrium element with quadratic stress field (HEE2), cubic stress field (HEE3), quartic stress field (HEE4) and for the nine-node displacement-based element (Q9), with meshes of 8 × 12 elements

The elastic dynamic problem of a cantilever beam
The accuracy of the high-order stress fields in the proposed formulation is tested under dynamic conditions for a simple cantilever beam subjected to a dynamic pulse. The problem is analyzed in Reference 28 with the left side fixed and with two different load setting: the axial and shear uniform pressures p n = p t = 1 applied at the free end. The uniform shear stress at the right end does not satisfy the local equilibrium condition at the two corners and intrinsically conflict with the proposed equilibrium formulation. Therefore, in the present article the shear load is applied as a quadratic tangential traction with null values at the two corners and maximum value p t = 1.5 at point B represented in Figure 11, according to the Jourawsky beam solution and with unitary tangential force (∫ h∕2 −h∕2 p t dy = 1), with h = 1 mm the beam thickness. Finally, a third load condition is applied in terms of a uniform vertical body force b y with a plane strain quasi-incompressible condition. The sizes of the cantilever beam and the three loading setting are represented in Figure 11. The material properties are the following: elastic modulus E = 200, Poisson ratio = 0.33, and mass density = 0.785. For the quasi-incompressible material the Poisson ratio is = 0.4999. The dynamic pulse load is applied by the time continuous loading law f 1 (t) represented in Figure 11, with two values of the pulse load application time: ΔT 1 = 8 s and ΔT 2 = 0.08 s. The body force f 1 (t)b y is applied only with the first impulse loading law, that is ΔT 1 = 8 s. The dynamic pulse load is also applied by the time discontinuous loading law f 2 (t) with pulse application time: ΔT = 4 s.
The dynamic analyses were performed with the DB nine-node element (Q9) and with the proposed dynamic HEE formulation with quadratic, cubic, and quartic stress fields. Two different meshes were considered: the 2 × 20 coarse mesh and the 10 × 100 fine mesh. The Q9 fine mesh was defined by 16 × 160 nine-node elements. The coarse meshes and the fine meshes are represented together with the displacement map in Figure 15 and with the map of normal stress in Figure 16. The number of nodes of the meshes is shown in Table 2. For the proposed meshes, the HEE formulation produces a spurious kinematic mode at the corner A which is restrained by a rigid link between the horizontal degrees of freedom of the two coincident nodes at the corner A (see Reference 18 for details).

F I G U R E 9
Maps of horizontal and vertical displacement of the quasi-incompressible Cook membrane test, for the hybrid equilibrium element with quadratic stress field (HEE2), cubic stress field (HEE3), quartic stress field (HEE4) and for the nine-node displacement-based element (Q9), with meshes of 8 × 12 elements

4.2.1
Slow axial pulse load The numerical simulation of elastic-dynamic analysis of the cantilever beam subjected to the uniform axial load p n is applied by the time continuous loading law f 1 (t) represented in Figure 11, with the greater value of pulse load application time: ΔT 1 = 8 s, with constant time increment Δt = ΔT 1 ∕160 = 5 ⋅ 10 −2 s and with t max = 50 s. The results of the dynamic analysis obtained with the coarse meshes of the dynamic HEE formulation and with the standard DB one, are compared in Figure 12A in terms of horizontal and vertical displacements and in terms of normal stress x , at point B represented in Figure 12B. The four different solutions (HHE2, HEE3, HEE4, and Q9) are practically coincident in terms of horizontal displacement u x at point B, whereas some small differences can be observed in terms of the vertical component u y , especially in the HEE2 formulation. This error is probably due to the asymmetry of the HEE meshes. Conversely, the HEE formulation computes the exact value of the normal stress x at point B, which is imposed as a boundary condition, and the DB formulation produces non null stress values after the pulse load, when the right end is unloaded and null normal stress is expected. The differences between the HEE and DB formulations are very small and vanish completely with the fine meshes, so that it can be stated that the two formulations converge to the exact elastic dynamic solution.

4.2.2
Fast axial pulse load The differences between the two formulations are more pronounced for the faster pulse load, which is defined by the continuous loading law p n f 1 (t) with ΔT = 0.08 s, with constant time increment Δt = ΔT 2 ∕160 = 5 ⋅ 10 −4 s and with t max = 0.50 s. The relevant numerical results, obtained with coarse and fine meshes, with the quadratic, cubic, and quartic HEE F I G U R E 10 Maps of principal stresses and maximum tangential stress of quasi-incompressible Cook membrane test, for the hybrid equilibrium element with quadratic stress field (HEE2), cubic stress field (HEE3), quartic stress field (HEE4) and for the nine-node displacement-based element (Q9), with meshes of 8 × 12 elements formulations and with the DB one are compared in Figure 13 in terms of horizontal and vertical displacement components at point A in the time domain.
The numerical solutions are also compared in Figure 14A-D in terms of evolution of stress components at point B in the time domain. The stress components x and xy are imposed as a boundary condition in the HEE formulations and coincide to the exact solutions for both coarse meshes and fine meshes. The DB solution provides some small errors in terms of stress which become negligible with the fine mesh. Finally, the evolution in the time domain of normal stress y , which is not imposed as a boundary condition in HEE formulation, is plotted in Figure 14B for the coarse mesh, with some differences between the different solutions, and in Figure 14D for the fine mesh, showing perfectly coincident results and convergence of the two formulations to the exact solution. In this sense, the proposed dynamic equilibrium-based formulation can also be considered as a dual approach for the error estimation and for the convergence analysis of the classic DB approach.

F I G U R E 11
The cantilever beam dynamic test with the three load setting: the axial uniform pressure p n , the quadratic tangential pressure p t , and the vertical body force b y . The loads are applied as dynamic pulse by the continuous loading law f 1 (t) and the discontinuous loading law f 2 (t)  The maps of horizontal displacement u x and the maps of normal stress x , at time t = 16.0 s, computed with the HEE formulation and the DB approach with coarse and fine meshes are plotted respectively in Figures 15 and 16. There it can be observed that displacement and stress computed with the coarse meshes do not coincide with each other and the equilibrium-based solutions do not show a symmetric response, due to the asymmetry of the HEE meshes. Moreover, some displacement discontinuity between coincident nodes is noticeable in the solution with quadratic stress fields (HEE2) but vanishes altogether with higher stress fields (HEE3 and HEE4). The numerical results for the fine meshes respect the symmetry condition of the axial load and are all perfectly coincident in terms of both displacement and stress,

Slow tangential pulse load
Good performances of the hybrid dynamic equilibrium formulation are also found in the numerical simulation of the cantilever beam subjected to a slow and continuous pulse tangential load p t f 1 (t), with the greater value of pulse load application time: ΔT 1 = 8 s, with constant time increment Δt = ΔT 1 ∕160 = 5 ⋅ 10 −2 s and with t max = 50 s. The results computed with coarse meshes are plotted in the two graphs in Figure 17 in terms of vertical displacement u y at point A and in terms of horizontal displacement u x at points A and B, showing negligible differences between the equilibrium-based solutions and the displacement-based one. The numerical results are also plotted in the two graphs in Figure 18 in terms of tangential stress xy at corner A, where the exact solution is null in the whole time domain, and at point B where the quadratic tangential load imposes the stress value xy = 1.5.
displacement u y at point A, and in terms of horizontal displacement u x at points A and B, for the quadratic, cubic, and quartic HEEs and for the nine-node DB element, with coarse 2 × 20 meshes F I G U R E 18 Dynamic response of a cantilever beam subjected to a slow tangential pulse load p t f 1 (t) (ΔT = 8 s), in terms of tangential stress xy at points A and B, for the quadratic, cubic, and quartic HEEs and for the nine-node DB element, with the coarse 2 × 10 meshes The numerical results for a cantilever subjected to a slow pulse tangential load confirms the good results of the proposed formulation, showing some slightly incorrect values of displacement u x at point B computed with the quadratic equilibrium formulation (HEE2). The cubic and quartic stress field equilibrium formulations produce exact values. By contrast, the equilibrium-based approach produces exact values of stresses at the free end, where tangential stress is applied as a boundary condition, and the DB solution computes incorrect values of the tangential stress for all time steps. The distribution of the three stress components at the right end at time step t = 4 s, computed with the coarse meshes of quadratic, cubic, and quartic equilibrium elements and with nine-node DB elements, are plotted in Figure 19. This graph clearly shows the limits of the classical DB models in stress response: linear tangential stress instead of quadratic trend; nonzero normal stress x at the free end; and nonzero normal stress y at the two corners. The differences between the three equilibrium formulation are due to the high order of the exact solution, which is defined by a cubic function for the normal stress component y . So the quadratic formulation can reproduce only an approximate solution, whereas the cubic formulation and the quartic one can reproduce the exact solution. The numerical results are accompanied by maps of the vertical displacement and maps of the tangential stress at time step t = 16 s plotted respectively in Figures 20 and  21, which were computed with the four considered formulations using coarse and fine meshes.

4.2.4
Fast tangential pulse load Similarly to the case of axial force, the differences between the equilibrium-based and displacement-based formulations are more pronounced for the faster tangential pulse load, which is defined by the pulse tangential load p t f 1 (t) , with  Figure 22 in terms of vertical and horizontal displacement components u x and u y at point A in the whole time domain, showing significant differences between the four formulations with coarse meshes, especially for the quadratic equilibrium-based one (HEE2). This differences completely disappear with fine meshes, as shown in the second graph in Figure 22, confirming the convergence of the proposed formulation to the exact solution.
The numerical results computed with coarse meshes are compared in the first graph in Figure 23 in terms of normal stress x in the whole time domain at point C, where such component cannot be imposed as a boundary condition. In the same graph the stress components y and xy computed at point C with the Q9 DB approach are plotted. These stresses should be null, being components of the traction at the free boundary, whereas they are exactly computed by the equilibrium formulation. The same stress components are plotted in the second graph in Figure 23 for fine meshes, showing the convergence of the four formulations to the same solution, although some residual error can be observed in the stress components y and xy computed with the DB nine-node element. The trend of tangential stress xy in the whole time domain at the internal point D, therefore not assumed as boundary condition, is plotted in the graph in Figure 24 for the coarse meshes and for the fine meshes.
The maps of vertical displacement u y and the maps of tangential stress xy , at time t = 0.16 s, computed with the HEE formulations and the DB one, with coarse and fine meshes, are plotted respectively in Figures 25 and 26. There it can be observed that some displacement discontinuity between coincident nodes is noticeable in the solution with coarse mesh and quadratic stress field (HEE2) but totally vanishes with the higher stress fields (HEE3 and HEE4) and also with fine meshes.
In order to analyze the numerical error introduced by the Newmark time integration, the numerical simulation of the cantilever beam subjected to the fast tangential pulse load is also performed with a greater time step value dt = ΔT 2 ∕16 = 5 ⋅ 10 −3 s and with the same time domain with t max = 0.5 s. The differences in the numerical solutions due to the greater time step are negligible with the fine meshes both for the HEE formulation and for the DB one. The greater time step does not produces any significant error neither with the coarse meshes and with comparable effects between the HEE formulation and the DB one. For the coarse meshes, the results of the numerical simulations performed with the two different values of the time step are compared in Figure 27A in terms of tangential stress at point D, and are compared in Figure 27B in terms of vertical displacement at the point A. Figure 27A,B shows that increasing the time step with the Newmark time integration method seems to produce similar effects in the HEE formulation and in the displacement based one, that is a delay in the response with an increase of the period.  Figure 29B, showing very similar results, and with the equilibrium solutions almost coincident to the converged one, obtained by the DB formulation with fine mesh (Q9 element).

4.2.6
Slow tangential pulse load for quasi-incompressible material Finally, the last numerical simulation of a cantilever beam subjected to a dynamic pulse uniform body force b y = 0.1 with a slow loading law (ΔT 1 = 8 s) was performed under a plane strain quasi-incompressible condition with Poisson ratio = 0.4999, for which the DB formulation suffers the volumetric locking numerical problem. The numerical results computed with coarse meshes are compared in the graph in Figure 30A in terms of normal stress x at point C, where such component cannot be imposed as a boundary condition, and the results are compared with the converged solution obtained by the fine mesh. The numerical results are plotted in the graph in Figure 30B in terms of tangential stress xy and normal stress x at the internal point D, and the results are compared with the converged solution obtained by the fine mesh. The normal stress x at point D in the exact solution is null in the whole time domain, and this result is confirmed by the HEE formulation. The two graphs also show the well-known problems of the DB formulation in evaluation of the normal stress components in quasi-incompressible elastic-dynamic problems.
The numerical results computed with coarse meshes are compared in terms of vertical displacement u y and horizontal displacement u x at point A in Figure 31A,B. In the same figure the results are compared with

CONCLUSIONS
The present article develops the two-dimensional hybrid equilibrium element formulation, with quadratic, cubic, and quartic stress fields, for accurate static and dynamic analyses of both compressible solids and quasi-incompressible ones. The formulation is developed in the variational framework of the minimum complementary energy principle for static analysis and in the variational framework of the Toupin principle, which is the complementary form of the Hamilton principle, for dynamic analysis. The proposed formulation is defined by independent stress fields for each finite element. The interelement and free boundary equilibrium conditions are applied using a classic hybrid formulation, and an independent displacement field is defined for each element side. The solution provides stress fields which are co-diffusive between adjacent elements and in equilibrium with traction at the free boundary sides. In the static formulation the stress fields verify the domain equilibrium equations. The dynamic formulation is developed under the hypothesis of null initial condition and is based on the impulse field (time integral of stress) and the dynamic equilibrium equation provides the pointwise velocity as the integral of the inertial term. Both the static formulation and the dynamic one are defined with high-order stress field and provide very accurate solutions in terms of stress. The analysis of a Cook membrane static problem and the analysis of a cantilever beam subjected to a pulse dynamic load were performed with several meshes for both compressible and quasi-incompressible elastic material.
The great accuracy of the stress-based proposed formulation is even more evident for static and dynamic analysis of quasi-incompressible materials, for which the classic displacement based formulation comes up against the volumetric locking. Moreover, the static and dynamic HEE formulations represent powerful numerical tools for analysis of elastic solids and also for dual analysis, as well as error estimation of the solutions performed with the classic displacement-based finite element formulations. The drawback of the proposed formulation is the possible presence of spurious kinematic modes, but they are well known and can be controlled or restrained by means of some different approaches.
The main future development of the proposed formulation could be modeling of the interelement fracture and fragmentation phenomena under dynamic load condition, through the approach proposed by the author in Reference 10, where an extrinsic (initially rigid) cohesive interface is embedded at any element side without any remeshing and without additional degrees of freedom. The extrinsic interface can activate once the stress based damaging condition is attained at the element side and the initially rigid behavior of the embedded interfaces does not affect the dynamic response of the pristine material. Conversely, the classic intrinsic interface with initial elastic behavior introduces additional compliance in the overall elastic behavior of a solid with relevant wave propagation issues.
Finally, the extension to nonzero initial conditions represents a basic requirement for analysis of a general elastic-dynamic problem and the static and dynamic analyses of pure incompressible elastic problem represent an interesting topic for further development of the HEE formulation. for the quartic formulation (n s = 4, n a = 25).

APPENDIX B. HIGH-ORDER DYNAMIC COEFFICIENT MATRICES
The two-dimensional impulse components of the elastic-dynamic problem are defined as functions of the Cartesian coordinates in Equation (30) for the quadratic formulation in the Voigt notation. For the quadratic formulation (n s = 2, n a = 18), the coefficient matrices S (2) e (x) is defined in Equation (30), the matrix of divergence of impulse T (2) e (x) is defined in Equation (33) and the vector of generalized variable is a (2) e = [a 1 a 2 a 3 ... a 16 a 17 a 18 ]. For the cubic formulation (n s = 3, n a = 30) the vector of generalized variable and the coefficient matrices can be written in the following compact