On a virtual element formulation for trusses and beams

The virtual element method (VEM) was developed not too long ago, starting with the paper [2] related to elasticity in solid mechanics. The virtual element method allows to revisit the construction of different elements; however, it has so far not applied to one-dimensional structures like trusses and beams. Here we study several VEM elements suitable for trusses and beams and show that the virtual element methodology produces elements that are equivalent to well-known finite elements but also elements that are different, especially for higher-order ansatz functions. It will be shown that these elements can be easily incorporated in classical finite element codes since they have the same number of unknowns as finite beam elements. Furthermore, the formulation allows to compute nonlinear structural problems undergoing large deflections and rotations.


Introduction
Beams and trusses are components of many engineering structures. Thus, many analytical and numerical simulation schemes were developed over the last century that can predict the displacements, deflections and normal, shear forces and bending moments in these structural members. Classical truss and beam models are described in any textbook on engineering mechanics, see, e.g., [6] and [4]. These models are based on kinematical assumptions which relate for beams to the names of Euler and Bernoulli for beams neglecting shear deformation. Equilibrium, kinematical relations and linear elastic behavior lead in the case of trusses to second-order ordinary differential equations and for beams to fourth-order ordinary differential equations in case of statics.
For the trusses, it is possible to use C 0 -continuous ansatz spaces that classically can be formulated in a finite element environment using the isoparametric concept. The Euler-Bernoulli theory requires for beams C 1 -continuous ansatz functions which, for a finite element formulation, is easily available using Hermite polynomials, see, e.g., [8], [16] and [12]. Thus, a virtual element formulation for trusses and beams does not necessarily provide any advantages. However, based on the virtual element method beam and truss elements of any order can be formulated easily as shown in Sects. 2.3 and 3.2. These higher-order elements provide better approximations of the deflections and stress resultants within the elements. An extension of the Euler-Bernoulli Fig. 1 Simple bar problem, geometry, data and discretization beam elements is provided by the Kirchhoff-Love theory for plate. Different formulations have already been discussed in [3], [5], [11] and [15] in this context .
The developed virtual truss and beam elements only introduce displacements, deflections and rotations as nodal unknowns at the end of the elements, even for higher-order ansatz spaces where additional unknowns have to be introduced, but these are only internal variables that are not attached to any node. Thus, the virtual elements can be easily incorporated in finite element software. Based on simple examples, it can be concluded that these structural virtual element are accurate and provide analytical solutions, for the kinematical variables and the stress resultants depending on the order of the ansatz.
An extension of the formulation to nonlinear problems in structural mechanics can be constructed in a straightforward way using the approaches derived for the linear problems. The related formulation is shown at the end of this contribution for different orders of ansatz functions.

Virtual element formulation of trusses
The truss or bar depicted in Fig. 1 can be modeled for constant cross section A by the ordinary differential equation where u(x) is the longitudinal displacement of the bar, E A the stiffness, f (x) the loading along the bar and x the coordinate. This differential equation can be recast in its weak form leading to Equivalently a potential can be formulated To solve the bar problem (1) with a discretization scheme, the bar of length l will be subdivided in n e elements of length l e such that n e i=1 l e = l, see Fig. 1. Either the weak form (2) or the potential (3) can be starting point for a discretization scheme.
In this simple example, for a linear ansatz, the finite element and the virtual element method yield the same matrices and hence produce the same results.
However, the methodology for deriving virtual elements is different from the finite element method. Mathematical details can be found, e.g., in [2] for solids in elasticity. It is based on the following ingredients: • An ansatz space in which the displacements are known at the nodal points i and i + 1. see Fig. 1, and at the element edges (which actually do not exist in our case of one-dimensional problems).
Here, to shorten notation, the projected part is defined by u π = [u h ]. The weighting function p in the first equation has the same polynomial order as the projection u π .

Finite element method for trusses
Here the ansatz functions are classically defined on the basis of an isoparametric mapping such that for an element e a linear ansatz is defined in a reference space −1 ≤ ξ ≤ 1, see Fig. 2, This ansatz approximates the displacement field in the element e and thus can be inserted, e.g., in (3) for all n e elements leading to where J e is the Jacobian of the isoparametric mapping J e = dx dξ = l e 2 . Both integrals can be evaluated using numerical integration. In the special case of the bar problem, the first integral is with simply given by the constant term for an element e By introducing the unknown vector u e = u 1 , u 2 T , the term in the square bracket can be written as (u 2 − u 1 )/l e = 1 l e −1 , 1 u e = B u e which yields the matrix form of (8) and the element stiffness matrix The second integral in (6) which relates to the potential of the distributed load f can be evaluated at element level using numerical integration, like the Gauss quadrature. The number of integration point depends then on the polynomial degree of f (ξ ), taking into account that u(ξ ) is a linear function.

Virtual element method with linear ansatz
In the standard isoparametric formulation used in the finite element method, an additional coordinate system is employed in the reference configuration, see (5). On the contrary, the virtual element method is formulated directly with respect to the coordinate system in the initial configuration. With a linear ansatz u π = a 1 + a 2 x, one can compute the gradient of the projected part, here the derivative u π = a 2 , from the orthogonality condition (4) 1 within an element where p is a weighting function that has the same polynomial degree as the ansatz u π . Since u h is not known within the element, the integral on the right side of (11) cannot be computed. However, based on the identity ( p u h ) = p u h + p u h (11) can be reformulated For the chosen linear ansatz p is constant as well as u π , it follows that p = 0. With u h (0) = u 1 and u h (l e ) = u 2 , see Fig. 2, we obtain the projected gradient from (12) which is now a function of the nodal displacements u 1 and u 2 . Thus, it is possible to compute the gradient u π without knowing the function u h inside the element e . In this simple case, the result u π matches exactly u h provided by the finite element method and hence the remainder is zero: u h − u π = 0. But in general the results are different, see e.g. [2]. Finally, u π can be inserted in (3) which provides the same result as given in (8). Thus, also the stiffness matrix of the virtual element K V E M is exactly the same as K F E M in (10).
A problem is now to compute the loading term in (3) since u h is not known inside the element. A possibility is to compute an approximation of u h by the projected displacement u π using the linear ansatz. Mathematically, it can be shown that this approximation of u h by u π in the loading term will result in optimal error estimates, see, e.g., [1].
The complete projection u π follows from (4) 2 where the average displacement u π is equal to the average displacement u h . Instead of evaluating the integrals, a sum of the values at the vertices of the element is employed. This yields for the two nodes of the bar element with x 1 = 0 and x 2 = l e 2 n=1 u π (x n ) = 2 n=1 u h (x n ) → a 1 + (a 1 + a 2 l e ) = u 1 + u 2 (14) and together with (13) we obtain a 1 = u 1 . Thus, u π = u 1 + u 2 −u 1 l e x. This function can be used to evaluate the integral associated with the loading potential in (3). As an example, we compute the potential energy U f e of the loading term for f (x) = f c = const.
which yields with the definition u e = u 1 , u 2 T the matrix form This loading term is exactly the same as the one for the finite element formulation in (6). Thus, so far the virtual element formulation produces equivalent results as the FEM. This will, however, change for higher-order ansatz functions.

Quadratic ansatz for a one-dimensional virtual truss element
The simplest possible higher-order ansatz function that can be used to derive a virtual truss element is provided by a quadratic ansatz for the problem E Au (x) = − f (x), see Fig. 2.
Analogous to the case of the linear ansatz, we select a quadratic function (ansatz order n = 2) for the projection This projection function has three unknown parameters and the question is: how will it be possible to derive a virtual element since we have only two points at the edges of the element with nodal displacements u 2 and u 1 ?
The solution follows from a closer look at the Galerkin projection (12) where now p has the same polynomial degree as (17) and thus is a quadratic polynomial with p = const. Hence the last integral u h dx on the right-hand side does not vanish. But it is also not computable since u h is not known inside the element. The way out is to define a new (internal) variable. This variable is not associated with any node and is called moment in the virtual element literature since for higher-order ansatz functions (n ≥ 2) integrals x n−2 u h dx up to the order n − 2 appear. These can be associated with moments in mechanics. Here we introduce which is scaled by the element length such that m 0 has the same dimensions as u 1 and u 2 . With this new variable the projection in (18) is computable, as we will see next.
It is convenient to introduce a matrix formulation in order to shorten notation. This leads to the derivative of u π With the matrix form p = 1 2x T , the left-hand side of (18) yields The polynomials can be integrated exactly leading to The right-hand side in (18) can be evaluated with p = 0 2 T resulting in With the unknown defined in (19), the right-hand side of (18) follows as The projection Eq. (18) can now be solved for the unknowns a by combining (22) and (23) Furthermore, the constant a 1 can be obtained by the condition that the average of the projection u π is equal to the average of the ansatz u h over the element This result can be written in a more compact matrix form by introducing the matrices together with the unknown vector u e = u 1 , u 2 , m 0 T and the projection operators leading to u π = N (2) π (x) P (2) u e and u π = ∇N (2) The derivative u π from (31) can now be inserted into the potential energy and integrated using the result from (22) Now the stiffness matrix of the virtual truss (T) element follows by differentiation with respect to u e as It is interesting to note that the VEM stiffness matrix is not equivalent to the FEM matrix for a quadratic element which can be computed from the potential (6) using the quadratic ansatz When looking at the eigenvalues ω of both stiffness matrices we obtain ω T, which are actually the two eigenvalues of the stiffness matrix for the virtual and the finite truss element with linear ansatz function. The zero eigenvalue is associated with the rigid body translation. Both stiffness matrices have the correct rank. Thus, a stabilization for the virtual element is not necessary.
The loading term follows from the potential, see (15), Interestingly, for the case of f (x) = f c = const. the integral u h dx can be approximated directly by the variable m 0 introduced in (19) leading to Thus the matrix form of the loading term is which is counterintuitive since no load term is assigned to the nodal degrees of freedom u 1 and u 2 . The loading term for finite elements is given for a constant load by which is clearly different to (37).
In order to verify the correctness of the virtual element formulation leading to (33) and (37), a simple example of a bar under constant load f c is considered that is fixed at the left end (u 1 = 0) , see Fig. 1. Using only one virtual element for the discretization yields with (33) and (37) This result can be introduced in (31) to compute the displacement u π and normal force N = E Au π which leads after some simple algebra to These results are equivalent to the exact analytical solution of the differential equation E Au = − f c . Thus the projected displacement u π of the quadratic virtual truss element delivers an exact solution in this special case. Actually, the exact analytical solution is also recovered when using the finite element formulation in (34) with the load vector in (38). We note that exact solutions are even recovered for the linear ansatz in Sect. 2.2, however, only for the nodal values, see, e.g., [13]. Finally, we observe that stiffness matrix and load vector of the virtual element are different from the finite element formulation and that the finite element formulation includes a third node at the center of the element.

Generalization for higher-order ansatz functions
The virtual element formulation for a truss element can be easily generalized to ansatz functions of any order. With the polynomial where n is the highest polynomial degree of the ansatz function, the projection can be computed. As before, a matrix B π (x) can be introduced The derivative of the ansatz u π is then given by for 1 ≤ k ≤ n. Furthermore the internal variable, moment, is introduced Here a scaling of the moment performed by l k−1 e such that m k−2 has the same dimension as the displacement. From the right-hand side of the projection (23), it is clear that moments only appear when the ansatz order is at least 2.
By using the same polynomial ansatz for p as for u π , the derivative yields in matrix form p = B T π = 1 2x · · · kx k−1 T and the left-hand side of (21) follows as All polynomials can be integrated exactly which yields The right-hand side in (23) can be evaluated with p using (42) and With the unknowns, u e = (u 1 , u 2 ), defined in the right part of Fig. 3, u h (0) = u 1 , u h (l e ) = u 2 and the moments m = m 0 · · · m k−2 the right-hand side of (47) follows for an ansatz up to order n with 2 ≤ k ≤ n as The projection equation (69) can now be solved for the unknowns a by combining (46) and (48) Once the parametersâ are known as function of the nodal and internal variables, the parameter a 1 can be obtained from (27). But due to the construction of the virtual element a 1 = u 1 since u π (0) = u 1 .

Virtual element formulations for Euler-Bernoulli beams
We consider a beam of length l with a stiffness of E I where E is the Young's modulus and I the moment of inertia of the cross section. The beam is loaded by a line load of magnitude q, see Fig. 3. The fourth-order ordinary differential equation for the deflection w of the Euler-Bernoulli beam with a constant cross section is given by By introducing a potential, the differential equation (50) can be written as a minimization problem As for the truss element, the approximation of the deflection w h is not known along the beam axis within a virtual element formulation. The unknown function w h is only known at the nodal points of the element. For a fourth-order differential equation not only the deflection but also the rotations have to be approximated. The unknowns of a beam element are depicted in the right-hand side of Fig. 3. The four unknowns w 1 , θ 1 , w 2 and θ 2 lead to a choice of a cubic ansatz which has four unknown coefficients a i As in the truss formulation, the parameters can be determined from a L 2 -projection. In case of the beam, we have to perform this for the curvature w Again p is a polynomial of the same order as w π . Since second derivatives appear on the left-hand side of the projection, only the coefficients a 3 and a 4 can be determined from (53). Differentiation of (52) yields which provides after integration a 2 × 2 matrix for the left side of the projection (53) with the vector containing the coefficientsâ = a 3 a 4 T . The right-hand side of the projection (53) cannot be computed directly since w h is not known along the beam axis. However, by using partial integration twice we are able to shift w h to single values at the boundary of the element which are known, In (52), we use a cubic polynomial hence the fourth derivative p is zero. Noting that the derivative p can be written as ∇B π = 0, 6 the matrix form of (56) follows with the right side of (54) as The rotation w h is known at the nodes of the beam element: w h (x = 0) = θ 1 and at w h (x = l e ) = θ 2 , see Fig. 3. The same is true for the deflection w h : w h (x = 0) = w 1 and w h (x = l e ) = w 2 . By introducing these relations into the above equation and combining the result with (55), the explicit matrix form of (53) is obtained. The solution of this equation yields the coefficients a 3 and a 4 as a function of the nodal unknowns w eâ with w T e = w 1 θ 1 w 2 θ 2 . Based on this result, the curvature w π can be expressed with (54) as w π = B πâ = B πP w e . Now the strain energy of the beam element (51) will be approximated using w π We note that the integral in this equation was already evaluated in (55). Thus, the stiffness matrix K B,V of the virtual Euler-Bernoulli beam element can be computed explicitly by a second derivative of the potential with respect to the nodal unknowns which is exactly the stiffness matrix that is obtained when using the finite element method with a cubic Hermitian ansatz function, see, e.g., [12].
What remains is to compute the loading term in the potential (51). Since w h is not known inside the element, it can be approximated by w π in (52). So far only a 3 and a 4 are known. To compute coefficients a 1 and a 2 , the relation (14) in Sect. 2 has to be extended to include also the derivative of the deflection (the rotations). The idea is to equalize not only the average of the nodal degrees of freedom of w h and the value of the projection w π at the element nodes but also the rotations w h and w π . Hence, the following conditions can be employed to compute the remaining coefficients By evaluating these equations at x = 0 and x = l e , we obtain for the deflections and the rotations By inserting a 3 and a 4 from (59), the result a 2 = θ 1 and a 1 = w 1 follows. This result could also be obtained by just looking at the ansatz w π in (52) since for x = 0 the ansatz function has the value w π (0) = a 1 = w 1 . The same is true for the first derivative w π which has to be for x = 0: w π (0) = a 2 = θ 1 . Now all coefficients are known as functions of the nodal degrees of freedom w e which can be expressed by the projection With (52) the ansatz w π is complete By inserting this ansatz into the loading term in (51), here written for the element e, we derive for q(x) = q 0 = const.
which is exactly the same result as provided by a cubic Hermitian finite element ansatz, see, e.g., [12]. Also in the case of the Euler-Bernoulli beam, the virtual element method leads to an identical result as the finite element method which has full rank. Thus, the remainder (w h − w π ) in w h = w π + (w h − w π ) is zero.

Fourth-order ansatz for a one-dimensional virtual beam element
This simple example will underline the methodology used to derive virtual elements for higher-order ansatz function. For this purpose, we will discuss a fourth-order (quadratic) ansatz for the problem E I w (x) = q(x), see Fig. 2.
Analogous to the case of the cubic projection, we select a fourth-order function with N a 1 , a 2 , a 3 , a 4 , a 5 . This projection function has five unknown parameters and the question is: how will it be possible to derive a virtual element since we have only two points at the edges of the element with nodal displacements w 1 and w 2 and the nodal rotations θ 1 and θ 2 ?
The solution is the same as in Sect. 2.3. We look at the Galerkin projection (53) together with (56) where now p has the same polynomial degree as (68) and thus is a fourth-order polynomial with p = const. Hence the last integral w h dx does not vanish. But it is also not computable since w h is not known inside the element. As in Sect. 2.3, we introduce an internal variable (so called moment) which is scaled by the element length such that m 0 has the same dimensions as w 1 and w 2 . With this new variable, the projection in (69) can be determined. It is convenient to introduce a matrix formulation in order to shorten notation. This leads to the second derivative of w π w π = 2a 3 + 6a 4 x + 12a 5 x 2 = B (4) π (x)â where B (4) π (x) = 2 6x 12x 2 andâ T = a 3 a 4 a 5 . By writing the second derivative of the polynomial p in matrix form p = [B (4) π ] T = 2 6x 12x 2 T the left-hand side of (69) yields The polynomials can be integrated exactly, leading to The right-hand side in (69) can be evaluated with p = 0 6 24x T and With the unknowns, defined in the right part of Fig. 3, w h (0) = w 1 , w h (l e ) = w 2 , w h (0) = θ 1 and w h (l e ) = θ 2 , the right-hand side of (74) follows as The projection equation (69) can now be solved for the unknowns a by combining (73) and (75) Furthermore, the constants a 1 and a 2 can be obtained by the conditions that the average of the projection w π and its derivative w π evaluated at the nodal points is equal to the average of the ansatz w h and its derivative w h at the nodal points for one element The two equations above yield a matrix system from which a 1 and a 2 can be determined By inserting a 3 to a 5 from (76) this equation system leads after some algebra to a 1 = w 1 and a 2 = θ 1 . The result is obvious since for x = 0 the ansatz function should give the nodal values w 1 and θ 1 at this point.
As in the case of the virtual beam element with cubic ansatz in section 3 a projector, see (65), can be defined which maps expresses the constants a T = a 1 , a 2 , a 3 , a 4 , a 5 in terms of the nodal values and the moment w T e = w 1 , θ 1 , w 2 , θ 2 , m 0 as The definitions in (68) and (71) The derivative w π from (83) can now be inserted into the potential energy and integrated using the result from (73) Now the stiffness matrix of the virtual element follows by differentiation with respect to u e as It is interesting to note that there is no adequate FEM formulation that yields the a finite beam element with 5 degrees of freedom. The loading term follows from the potential, see (51), Interestingly, for the case of q(x) = q 0 = const. the integral w h dx can be evaluated directly using the definition (70) of the variable m 0 . This yields Thus, the matrix form of the loading term is simply given by which is similar to the constant loading term of the truss element with quadratic ansatz function, see (37). One can evaluate U q in a different way by inserting the ansatz for w π , see (83) into the integral in (87). This will actually lead after some algebraic manipulations to exactly the same result. For a load that varies linearly in the virtual element given as q(x) = (1 − x l e ) q 1 + x l e q 2 with two constants q 1 and q 2 describing the load at the nodes, the potential of this loading term is obtained by inserting w π from (83) This result yields the load vector for one element For an assessment of the accuracy of the developed element, a beam that is clamped at the left side is investigated. Only one element is considered for the case of the constant load, q 1 = q 2 = q 0 leading to f 0 , and the linear varying load, f 1 . The boundary conditions for the clamped beam are (w 1 = θ 1 = 0) , see Fig. 4. In that case, the element stiffness matrix (135) can be used where the first and second rows and columns are eliminated. Furthermore, the first two rows of the load vector are deleted leading for the constant load to the equation system with the solution ⎡ The nodal values w 2 and θ 2 match the analytical results. Furthermore, the solution can be used to compute the moment distribution along the beam M(x) = −E I w π . For that, the results of (91) are introduced in (83) which leads to Also this result matches the analytical solution which is related to the fact that the analytical solution is a fourth-order polynomial, and hence, the fourth-order ansatz (68) Here w 2 and θ 2 coincide with the analytical solution at the beam end which actually can be shown for all complete polynomial ansatz functions that fulfil the homogeneous differential equation of the beam. The result (93) leads to bending moment distribution in the beam that does not match the distribution computed from the analytical solution

Generalization for higher-order ansatz functions
The virtual element formulation for Euler-Bernoulli beams can be generalized for higher-order ansatz function. For that, an ansatz is introduced where n is the highest polynomial degree of the ansatz function. The matrix B π (x) in (71) can be written as with 2 ≤ k ≤ n and c 1 (k) = k(k − 1). The second derivative of the ansatz is then given by for 2 ≤ k ≤ n. Furthermore, we can introduce the moment, which is an internal variable, where the scaling by l k−3 e is made such that m k−4 has the same dimension as the deflection. From the right-hand side of the projection (69), it is clear that moments only appear when the ansatz order is at least 4. By using the same polynomial ansatz for p as for w π , the second derivative yields in matrix form p = B T π = 2 6x · · · k(k − 1)x k−2 T and the left-hand side of (69) follows as The polynomials can be integrated exactly, leading to The right-hand side in (69) can be evaluated with p using (98), p = 0 6 . . . c 2 (k)x k−3 T and p = 0 0 24 . . . c 3 (k)x k−4 T where p and p only exist for k ≥ 3 and k ≥ 4, respectively. Here the constants , depending on the ansatz order, were introduced. This results in . . .
With the unknowns, w e = (u 1 , θ 1 , u 2 , θ 2 ), defined in the right part of Fig. 3, = θ 2 and the moments m = m 0 · · · m k−4 the right-hand side of (103) follow for an ansatz up to order n with 3 ≤ k ≤ n as The projection equation (69) can now be solved for the unknowns a by combining (102) and (104) Once the parametersâ are known as function of the nodal and internal variables, the two parameters a 1 and a 2 follow from (77). But as discussed in the previous section, due to the construction of the ansatz we arrive always at a 1 = w 1 and a 2 = θ 1 . With this the projection function, w π is completely determined.
For an ansatz of order n ≥ 4, we will obtain an exact solution for a load of q(x) = n k=4 q k−4 x k−4 . In that case, the potential related to the load can be written as Thus, it is straightforward to develop Euler-Bernoulli beam elements of arbitrary order. As an example, we provide the relevant matrices for a fifth-order ansatz in appendix.

Static condensation
Internal variables can be removed by static condensation. Thus, it is possible to eliminate the moments m i at element level. For that, we split the matrix and right had side at element level in and build the Schur complement leading withK I I = K I I − K I M K −1 M M K M I tō It can be shown for any order of the ansatz used to derive the virtual element that the matrixK I I is equivalent to the matrix representing the cubic ansatz given in (61). This is consistent with the result that the third-order ansatz exactly solves the homogeneous part of the beam equation (50) and it is a well known fact that the nodal deflections w i and the nodal rotations θ i are exact for any given right-hand side, see, e.g., [13] and [7].
However, the approximation of the deflection within the cubic beam element are only exact for q(x) = 0. Hence, as shown before, it makes sense to use the virtual element formulation to obtain higher-order approximations and thus better a representation of the beam deflection and with this also of the bending moment and shear force within the element.
Based on this observation, it is possible to analyze a beam problem by computing 1. the nodal deflections and rotations using the cubic beam formulation which is equivalent to the finite element beam. This yields the exact nodal degrees of freedom u I . 2. Next the moments m M follow from ((108)) 2 which complete the unknowns for the virtual element. 3. The higher-order approximation of the deflection w π within the element can then be computed using the projector, e.g., (81) for the fourth or (133) for the fifth-order virtual element ansatz.
This approach is similar to a postprocessing step where first the problem is solved with the classical beam element, and then, a step follows that produces a higher-order solution at local element level using the virtual element formulation.
The same procedure can also be applied for the truss element. Here a linear ansatz solves the homogeneous Eq. (1). Thus, again an elimination of the internal degrees of freedom can be performed by using (108) 2 which yields the stiffness matrix (10) related to the linear ansatz function. Again the same procedure, as mentioned above for the beam, can be used to compute the nodal degrees of freedom u i and moments m k which yields an enhanced approximation of u π within the truss element by the virtual element formulation.

Application of the virtual element discretization to nonlinear beam formulations
The developed virtual beam elements can be applied in linear and nonlinear problems related to structural analysis. Here we use the methodology of Sects. 2 and 3 to construct a formulation for nonlinear response of Bernoulli-Euler beams. We assume that the beams undergo small strains but finite deflections and rotations. In this case, the static condensation approach, see above, cannot be employed since the nonlinear truss and beam equations will not be solved exactly by a linear and cubic polynomial, respectively.

Theoretical background and discretization
The associated strain measures for the axial strains ε and the curvature κ that can be found in textbooks, like [14]. The elongation ε and the curvature κ are given by Thus, the element formulations in Sects. 2 and 3 have to be combined in order to approximate the displacement u and the deflection w of the beam. Attention has to be paid that second derivatives of u appear in the curvature term which is contrary to the linear beam theory. In case or a linear constitutive relation, the potential describing the virtual beam element can be written in matrix form as with the definitions Note that U (u) is highly nonlinear in terms of the deflection and axial displacement. Now the strains have to be approximated by the virtual element ansatz functions where many different possibilities exist for the choice of the order of the ansatz functions. We will vary the ansatz order n u for axial displacements from linear to cubic (n u = 1, 2, 3) and the ansatz order for the deflection n w from 3 rd to 4 th order (n w = 3, 4).
Next we provide the the projection matrices related to the axial displacement. The first-order part is given by For the second-order ansatz, the projection operators follow from (31) For the third-order ansatz, the projection can be obtained from the matrices in Appendix A.1. For that the equation, system with (124) and (125) has to be solved leading to The first and second derivatives needed in (109) follow as The same procedure can now be applied for the approximation of the deflection w using the third-and fourth-order projectors in (65) and (81). Again we can write where now the ansatz polynomial N (n w ) π and the projection matrix P (n w ) are related to the virtual beam elements in Sects. 3.1 and 3.2. The first and second derivatives needed in (109) follow as w π (n w ) = N π (n w ) (x) P (n w ) w (n w ) e and w π (n w ) = N π (n w ) (x) P (n w ) w (n w ) The previous set of equations is sufficient to establish discretizations for beams in one dimension since it is related to the local coordinate system of the straight beam axis. For more complex structures in which the beams are located in different positions, the matrices and vectors have to be transformed to a global coordinate system. This can be performed by a transformation of the local nodal displacements and rotations u loc I at a node I which then can be expressed in terms of the global variables u The internal variable m i is not affected by this transformation. The ansatz functions in (115) and (117) and their derivatives in (116) and (118) can now be inserted in (109). This yields the strains in terms of the local nodal unknowns and internal variables. These can be transformed via (119) to the global variables. Finally, the strain measures are inserted into the potential (110) The approximation of the strain measures for an ansatz of order n u , n w is given by The displacement approximation which is needed in the loading term in (120) is provided in terms of the nodal unknowns by u T ≈ u T = N (n u ) π P (n u ) u (n u ) e , N (n w ) π P (n w ) w (n w ) e . Insertion of (121) into the potential (120) yields a nonlinear function of the nodal degrees of freedom U ( u) from which the residual and tangent matrix follow for the selected ansatz order The differentiations are performed using automatic code generation tool AceGen, developed by [9], see also [10], which then generates the residual vector and tangent matrix.

Numerical example and comparison of different discretizations
With an example, we want to demonstrate the behavior of different formulations in the nonlinear range. A frame as depicted in Fig. 5a is loaded by a point load in vertical direction. It is simply supported at the bottom and the top right. The material data are given as follows: bending stiffness E I = 2 · 10 5 , axial stiffness E A = 10 6 . The width of the frame is equal to its height (l = h = 300); 20 beam elements are applied to discretized the structure, see Fig. 5b. The solution using virtual elements with different ansatz functions is studied in this example. We apply a classical Euler-Bernoulli beam element with linear ansatz for the axial displacement and cubic ansatz for the deflection. This element is named "1st/3rd". Such ansatz is not optimal since it will neglect the terms related  (109). To account for a better approximation of u, an element is generated that uses a cubic ansatz for the axial displacement besides the cubic one for the deflection. This element is labeled "3rd/3rd". In order to see how a higher-order ansatz for the deflection affects the convergence behavior of the virtual beam elements a 4 th -order ansatz is selected for the deflection together with a quadratic and cubic ansatz for the axial displacement. The related elements are named "2nd/4th" and "3rd/4th," respectively. For all formulations, the same discretization is used in order to investigate the effect of different ansatz orders.
All these element formulations are based on the Euler-Bernoulli theory. The results will be compared with with a simple finite element based on a Timoshenko theory, see, e.g., [14], which has a linear ansatz for the axial displacement, the deflection and the rotation. It is denoted by "1st/1st.". Different formulations are compared with respect to the response of the frame due to an increasing load. The load is increased in steps of F = 5 which yields the load deflection curves depicted in Fig. 6 It is interesting to see that the virtual element "1st/3rd" which does not approximate all terms in the strain measure correctly yields a too stiff response. This is also true for the "2nd/3rd" and "2nd/4th" virtual element which has a quadratic ansatz for the axial displacement and only yields a constant term for u . Thus, the axial displacement in the nonlinear strain measure has to be approximated with at least a cubic ansatz to obtain results that are equivalent to the Timoshenko-Reissner finite beam element with linear ansatz functions for all terms for the given coarse mesh resolution. However, by raising the ansatz order for the deflection the virtual element produces better results which are also better than the ones using Timoshenko beam element with the same discretization. We remark that the result of the discretization using the "3rd/4th" element with 20 virtual elements is equivalent to the converged solution of the Timoshenko-Reissner beam formulation with 2000 finite elements, as depicted in Fig. 6.
The deformed configuration at a load level of F = 45, solved by using the "3rd/4th" virtual element, is shown in Fig. 7 which underlines the capability of the virtual elements to model nonlinear behavior leading to large deflections and rotations with high coarse mesh accuracy.

Conclusion
Virtual element formulations were developed for truss and beam elements using different ansatz orders. The elements have the advantage that only the nodal values at the ends of the elements enter the formulation together with internal variables (moments) for higher-order ansatz functions. For the ansatz with lowest possible order (truss linear, beam cubic), the virtual elements coincide with classical finite elements. For higher-order ansatz functions, the virtual elements differ from known finite element formulations. The new virtual elements can be applied to linear and nonlinear problems.
Funding Open Access funding enabled and organized by Projekt DEAL. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Appendix
In the appendix, we summarize two higher-order element formulations for the truss and the beam element.

A.1 Cubic ansatz for a virtual truss element
The necessary equations to obtain a fifth-order beam element are summarized below.
A.2 Fifth-order ansatz for a virtual beam element The necessary equations to obtain a fifth-order beam element are summarized below.