Finite-strain three-dimensional solids with rotational degrees of freedom: non-linear statics and dynamics

In this paper we seek refined yet efficient computational models of large overall motion in statics and dynamics. The efficiency is achieved by the proposed model of 8-node brick element with rotational degrees of freedom which allows to separate large displacements and large rotations. The independent rotation field leads to an intrinsic representation of the rotation tensor, ensuring a smooth interaction between 3D solids and beam elements. The element is based on a sound variational formulation and the incompatible mode method which allows to construct enhanced strain representation. Several numerical examples are presented to show an excellent performance of this element in the whole range of large overall motion in the statics and dynamics problems.


Introduction
Solid elements with low-order interpolations are generally recommended in structural mechanics because they can be used efficiently in nonlinear applications since they have a more robust performance in the distorted configurations. However, in many cases, especially in bending dominated problems, standard brick elements show severe stiffening effects known as locking problems [1,2]. Several methods have been proposed over years in order to overcome these locking phenomena. For example, the enhanced assumed strain (EAS) have been used in geometrically nonlinear version [3][4][5] where the strain field can be enriched in order to improve the element's performance under certain conditions. 3D solid elements with rotational degrees of freedom are also proposed with different discrete approximations like the method of the mixed interpolations of tensorial components [6] for the brick element of Wils on and Ibrahimbegovic [2]. The Space Fiber Rotation concept (SFR) [7] and the shift of the mid-side displacement DOFs of the classical 20-node hexahedral element into corner nodal translations and rotations are proposed by Yunu and al. [8].
In this paper, we propose the method of incompatible modes for constructing enhanced strain approximation as the most viable approach for improving the performance of low order elements. A very important choice pertains to large strain measures, which allows to separate large displacements and large rotations. We illustrate the proposed approach on an 8-node displacement-based hexahedron, including both usual displacement degrees of freedom and rotation degrees of freedom that are independently interpolated. We recapitulate the variational formulation including the incompatible displacement modes for 3D finite displacement elasticity problems, along with the fine details of the numerical implementation.
Given the goal of using this element in dynamics, we introduce the inertial effects in the variational formulation in order to handle dynamic analysis. It is important to state that the proposed formulation is shown to be more efficient than others formulations using floating frames. In fact, it is set in a fixed frame that allows to eliminate the Coriolis effects and lead to a simple quadratic form of the kinetic energy. This results with a mass matrix with constant components and consequently simplifies the time-integration computational phase.
The proposed enhanced solid element is compatible with shell finite elements [9] as well as beam elements [10]. This offers an efficient performance for modeling linear and nonlinear dynamic behavior of complex structures. The rotational degrees of freedom are presented by orthogonal tensor, as an intrinsic representation of large rotations.
The outline of the paper is as follows. In "Standard variational formulations for statics" section we present the variational formulations for the continuum with independent rotation fields in geometrically nonlinear theory. Then, we discuss the regularized form of the variational formulations, with some details of numerical implementation in "Modified variational formulations for statics and its discrete approximation" section. We extend to dynamics analysis in "Variational formulations for dynamics" section. The solution procedures are developed to obtain the corresponding solutions in "Newmark implicit time-stepping scheme" section. Several numerical simulations dealing with static and dynamic problems are presented in "Numerical examples" section. Some closing remarks are given in "Closing remarks" section.

Standard variational formulations for statics
In this section we discuss the variational formulations for the continuum with independent rotation field in geometrically nonlinear theory. The starting point in our considerations can be provided by the classical potential energy principle. It is a function of the finite displacement field u, writing as a sum of the strain energy and the potential energy of external forces where W (H(u)) is a stored energy given as a function of the chosen Biot's finite strain measure and u is the finite displacement field. The second integral represents the external work for the Dirichlet boundary value problem. The finite strain measure H, often called Biot strain [11], can be the most explicitly defined via the polar decomposition theorem (e.g., see [12], p. 14). Namely, if the deformation is a vector field ϕ, which is, without loss of generality for our purposes, specified with respect to the Euclidean coordinate system with the base vectors e i , i.e. if then, the deformation gradient can be written as where ⊗ denotes the tensor product. Introducing the displacement vector field: u = ϕ(x) − x, the deformation gradient can be rewritten in terms of the displacement gradient ∇u as The polar decomposition theorem [13] states that the deformation gradient can be factored in a unique way into F = RU where U is the right stretch tensor describing deformation, while R is the orthogonal rotation tensor. Then, the Biot strain tensor, defined with: H = U − I, is used to rewrite the polar decomposition theorem By eliminating the rotation field in (5) via orthogonality of R, we get a functional relationship between H and u H + which is needed in (1). The relation in (6) is quite complex given geometrically nonlinear nature of the problem. Thus, we can choose much simpler formulation if the Biot strain is included in the functional (1).
Hence, we want that the polar decomposition in (5) be recovered as the Euler Lagrange equation of a new variational formulation rather than having it as a subsidiary condition of the variational formulation in (1). The Lagrange multiplier procedure can be used to impose that condition in the form P is the non-symmetric Piola-Kirchhoff stress, which plays the role of the Lagrange multiplier [14]. The weak form of the polar decomposition in (7) can be added to the variational formulation in (1) The formulation we propose, considers Biot strain and stress measures. The Biot stress tensor is obtained by the pull-back of the first Piola-Kirchhoff stress P, by using the rotation part R of the deformation gradient: T = R T P. The variational formulation in (8) can be written as The Euler-Lagrange equations associated with the principle in (9) can be obtained by taking the directional derivative in the direction of virtual displacements δu, virtual rotations δR, virtual stresses δT and virtual strains δH. We thus obtain: linear momentum balance, angular momentum balance, definition of strains H and rotations R and constitutive equations We use the Legendre transform to introduce the complementary energy Σ(symmT) By introducing the result in (11) into the variational formulation in (9), we can eliminate the strain field H to get the three-field variational formulation The Euler-Lagrange equations in (10) remain preserved, apart from constitutive equations, which connects directly the stresses with the displacements and rotations

Enhanced displacement gradient
In seeking to improve performance of the discrete approximation, we consider the incompatible modes method to reparameterize the assumed displacement gradient field in terms of additive decomposition of 'compatible' ∇u and enhanced ('incompatible') part d. In particular, we assume the enhanced displacement gradient given as The condition above can be included in the variational principle in (12) via Lagrange multiplier procedure to get where P is the first Piola-Kirchoff stress tensor, conjugate with the enhanced displacement gradient d.
Note that in the variational principle in (15) it is sufficient that the enhanced displacement gradient belongs only to the space of square-integrable functions in the domain V denoted L 2 (V). Thus, the enhanced displacement gradient in the finite element approximation can be discontinuous over an element boundaries ('incompatible'). We note in passing that the strong form and the Euler-Lagrange equation will not change, since

Regularized variational formulations
We assume that the complementary energy potential can be constructed as a quadratic form in terms of stress where C is the elasticity tensor defined by its standard form [13]. Hence, the variational principle in (15) becomes The variational principle needs to be regularized in order to preserve stability. In geometrically linear theory, Hughes and Brezzi proposed a regularized form of the variational principle [14] in order to be able to use any convenient discrete approximation. The corresponding generalization of their proposal for the present geometrically nonlinear case can be written as skewT · γ −1 skewT dV (19) where γ is a regularization parameter. An optimal value of γ , γ = μ, was identified in geometrically linear cases [14]. The regularized form of the variational principle preserves the Euler-Lagrange equations (10) and (13), while producing an additional Euler-Lagrange equation which is given as By means of Eq. (20), the regularized variational principle can be obtained featuring only the kinematics variables being able to reduce to a minimum the number of unknown fields.

Variational equations in compact notation
For computational efficiency, we will further switch to the matrix notation. The variational principle in (21) can be restated as where the symmetric part of the Biot strain measure is mapped into a 6-dimensional vector while the shew-symmetric part is mapped into a 3-dimensional vector Here we used the following compact notation The perturbed configuration is defined by u = u + δu and R = exp( δW)R with δu as the virtual displacement and δW as the infinitesimal skew-symmetric tensor of virtual rotation where We note the key difference in dealing with large displacements (vector) and large rotations (tensor) in this kind of computation (e.g. see [15] for more elaborate presentations). This allows to write the variation of strain measures where Y(u) and D are defined as follows We further simplify the writing with virtual displacement and virtual rotation vectors grouped together in δa T = δu T , δw T . The variational equations (principle of virtual work) that follow from (22) are obtained by the directional derivative in the direction of virtual displacement and virtual rotation δa, virtual enhanced displacement gradient δd and virtual stresses δp.

Finite element incompatible mode interpolations
We deal here with the discrete problem, defined with the finite element approximations. The displacement and rotation fields are approximated by using the standard isoparametric interpolations for three-dimensional solid elements with 8 nodes (30) where N I are the standard shape functions while u I and w I are the corresponding nodal values. The virtual displacements δu and virtual rotations δw are approximated in the same manner With these interpolations, the displacement gradient of the virtual displacement field can be written as The interpolation of the enhanced displacement gradient is constructed by using derivatives of the incompatible displacement α [2] are chosen as quadratic polynomials. This choice is made to enhance the bending dominated performance [2].
The proposed interpolation should be modified to ensure that the incompatible modes are not active for constant stress, which guarantees the convergence of the incompatible mode method in the spirit of the patch test [16]; we impose: which will in turn eliminate the variable p from the variational equations. For a given interpolationĜ, we can always construct the modified interpolation enforcing (34) as shown in [17] .
Such modified incompatible mode interpolation will in turn eliminate the variable p from the variational equations. We should thus compute the stress tensor as follows which further reduces the variational equations in (29) to a set of equilibrium equations, which can be written as One possibility to solve the non-linear system in (37) is to linearize the complete system, as presented in [18] for a two-dimensional case. In present 3D case, the linearization leads to more complexity and involves the use of the secondary storage. Note however that the system size is not increased thanks to using the operator split method, which is presented in detail in Appendix I.

Variational formulations for dynamics
The equations of motion for the transient problem is based upon the variational formulations, by appealing to the Hamilton principle. In fact, the only new term with respect to the variational formulations in statics concerns the directional derivative of the kinetic energy, written as The kinetic energy in (39) is chosen as quadratic form only in terms of displacement fields. The independent rotation field is not involved, for it only serves to introduce the Biot strain measures. Thus, the variational equations in dynamics are written as modified form of these in (29) assuming for inertial effects In dynamics, displacement and rotation fields are functions of both space and time. We use the separation of variables approach in order to construct the finite element approximations of displacements and rotations The variational equations in (40) can now be rewritten as with the accelerations fieldü interpolated in the same way as the compatible displacements field. It is easy to see that the mass matrix will have constant entries and take the following form We note that the zero masses associated to angular accelerationsẅ can prove troublesome and affect computation stability in dynamic problems. In order to illustrate that clearly, we consider the linearized form of equations of motion in (42) (1) for free vibration case: Uexp(i ωt) and Wexp(i ωt), which leads to This clearly shows by choosing: U = 0, W = 0 that the presence of zero terms in the mass matrix would require infinite frequencies, which is the ultimate case of stiff equations [19] with all difficulties that will impose.
The simplest way to overcome this deficiency is by using the penalty method and introducing the mass matrix contribution of the rotational degrees of freedom through a kind of regularization. This contribution is assigned to be equal to the ones coming from translational degrees of freedom, multiplied by a regularization parameter η which varies from 0 and 1. Finally, the regularized mass matrix can be written as

Newmark implicit time-stepping scheme
For nonlinear problems, the evolution of state variables is obtained by step-by-step integration schemes. To that end, the time interval of interest is partitioned into a number of time steps (0 ≺ t 1 ≺ t 2 ≺ · · · ≺ t n ≺ t n+1 ≺ · · · ≺ T ). At a typical time t n , the values of u n and R n are known. The corresponding values at time t n+1 are computed in each step independently for single-step schemes, such as Newmark. For displacement vector we have sample additive updates where Δu (i) n+1 is the incremental displacement at each iteration (i). The rotation update is somewhat more involved in that we have to choose between various possibilities of parameters for rotation representation (e.g. see [15]). If the spatial representation is used, we can carry out the rotation update as where Δθ (i) n+1 is the incremental rotation vector at each iteration (i) and Λ(•) corresponds to a direct representation of the orthogonal rotation tensor via the rotation vector.
The intrinsic representation of the finite rotations by an orthogonal tensor can be reduced to a set of four quaternion parameters. The rotational update can be written as The quaternion parameters of the iterative rotation parameter of the orthogonal tensor {q where Δw n+1 is the axial vector of incremental rotation. Note that Δw n+1 and Δθ n+1 are interconnected The scalar θ is the euclidean norm of θ and the tensor Θ is the shew-symmetric matrix associated to θ.
In dynamics, besides computation of displacements and rotations, one also needs to provide the values of velocities and accelerations at each time step. In what follows, the Newmark time integration scheme is used. The linear velocity and acceleration are advanced from time t n to t n+1 by using the standard Newmark approximationṡ where Δt is the time step while β and γ are the Newmark coefficients. The Newmark implementation for finite rotations is a bit more laborious [20]. Here, we follow previous works on 3D beams [15] with an extension of the standard Newmark algorithm valid only for so-called spatial incremental rotation vector θ n+1 . We assume that the angular velocity and acceleration, respectivelyẇ n+1 =ẇ(t n+1 ) andẅ n+1 =ẅ(t n+1 ), can be provided by the Newmark approximations for finite rotationṡ whereω n+1 andα n+1 are given bỹ From the proposed form of the Newmark approximations for both displacements and rotations, we can compute the linearization of velocities and accelerations, providing corresponding iterative updates. First, we can write for displacementṡ Similarly, angular velocity and acceleration iterative updates are obtained using the iterative incremental rotation vectoṙ where Λ(•) is given by the exponential mapping formula of Rodrigues [15], written as By exploiting the relationship presented in (51), angular velocity and acceleration updates becomė where we use the result in [15], It is interesting to note that the rotation updates in terms of incremental rotation vector can be easily computed thanks to the additive procedure of updates, θ Moreover, velocity and acceleration updates in (55) and (57) have formally the same structure for both linear and angular configurations which is the main advantage of the proposed Newmark algorithm for finite rotation.
The linearized form of the equations of motion in dynamics is given as

Numerical examples
Several numerical examples are presented in order to demonstrate a very satisfying performance of the enhanced 3D solid element proposed herein. The presented simulation results concern both static and dynamic problems. All the computations are performed with a research version of the computer program FEAP, written by Prof. R.L. Taylor at UC Berkeley [16].

Short cantilever beam in large compressive deformation: static problem
The first example presents a comparison between the proposed solid element and a standard solid element, which employs the Green-Lagrange strain measure. A single finite element cantilever beam is subjected to large compression, by imposing the displacement at one end while keeping the other end fixed (see Fig. 1). The material properties are selected as: E = 2000 N/m 2 and μ = 0. The imposed displacement (u x = −0.99) is increased in 10 equal increments until practically reaching maximum possible compressive strains. The plots of the xx-component of the PiolaKirchhoff stress as function of stretch are presented in Fig. 1 for both models. We can easily see that the zero stress value accompanying maximum compressive strain for the Saint-Venant-Kirchhoff material cannot be justified for any real material. In fact, a large compression must logically be accompanied by a large value of stress which is shown, moreover, by the results of the proposed model.
The deficiency in representing very large compressive strains by the standard solid element is due to the loss of convexity of the strain energy. However, this drawback of violating poly-convexity conditions is easily repaired by the proposed element. For this, it is useless to compare the proposed element to others solid element, which employ Green Lagrange deformations measures.

Large deflection of a cantilever: static problem
In this second validation example, we consider the cantilever beam of Fig. 2 clamped at one end and loaded, at the another, with four equal concentrated moments which add to the total value of m = 0.01π Nm. The beam is modeled with 10 × 1 × 1 enhanced 3D elements.  The exact rotation angle θ can be computed with the classical Euler formula based upon the classical beam theory [10]: θ = m/EI. For the given value of the moment, geometric and material proprieties the reference solution of the deformed shape is a semi circle.
The deformed shape is presented in Fig. 2, showing an excellent agreement with the reference shape. The loading is applied in a single load step, and the solution is obtained after 7 iterations, with a satisfying rate of convergence (see Table 1).

Circular arch under single load: static instability problem
The next example shows the efficient use of the proposed enhanced 3D solid element to solve a geometric instability problem, first solved by 3D beam. We consider a circular arch, hinged at one end and clamped at another, under a vertical point load in the middle. The selected properties of the arch are given in Fig. 3a. The numerical solution is obtained by a mesh having thirty 8-node enhanced elements.
The computation is performed with a time step Δt = 0.01 s and needs only 5 iterations to converge. The corresponding result is also computed with the same number of 2-node beam elements (see [10]). In Fig. 4, the load/displacement curves for both approximations are given, indicating a very good agreement for two curves obtained with the different finite element models.

Beam distortion by a solid plate: static problem
In order to show the compatibility of the proposed element with beam element, we propose the following example. We consider a straight beam of length l = 4 m, aligned with the z-axis. The finite element model of the beam consists of five 3D geometrically exact beam elements. The mechanical properties of the beam are chosen as: EA = 100 N, GA = 50 N

Free vibration analysis: dynamic problem
For the dynamics part, the first example is proposed for free vibration analysis. For a cantilever beam fixed at one end, we are concerned with the identification of natural frequencies and mode shapes. This beam is modeled with 10×1×1 enhanced 3D elements (Fig. 6). Table 2 gives the natural frequencies of the first ten modes, sorted in increasing order for both enhanced solid element and standard beam element. The results reveal that solid element detects mode shapes of distortion which are unseen by the beam element. We can also see that the natural frequency values of both elements correspond well to each other. Duplicated frequencies are justified by the fact that each bending mode in y direction has a similar bending mode in z direction.
The mode shapes obtained with mesh of elements are presented in Fig. 7. Compared to the standard beam element, the proposed element is capable of taking into account sectional changes as showing in the mode shapes 3 and 6. In order to captive that in the case of beam, one needs higher order geometrically-exact beam models [21,22], incorporating in-plane cross sectional changes and out-plane warping by adding supplementary degrees

Cantilever beam subjected to end shear force: dynamic problem
In the second version of the same test, we keep all the data the same as in the first version, except that we change the loading proprieties. First, the beam is subjected to four vertical forces at the free end. Each force has a weak amplitude (f = 0.001 N) in order to assimilate linear behavior and a sinusoidal time variation with the frequency ω = 3.12413236 × 10 −1 rad/s, which is the first natural frequency in the initial configuration. The vertical displacement under force applied at the free end is presented in Fig. 8. The resonance condition is satisfied and the dynamic response presents rapidly increasing as expected.
We have proposed two alternatives to define the mass matrix. In order to compare these possibilities, we examine the bending problem of a cantilever beam under four concentrated forces applied at its free ends. The proprieties of the beam are maintained as before (Fig. 6). At each node, the applied force is a triangular pulse, linearly increasing for 0 ≤ t ≤ 1 to attain the maximum value F = 0.5 N, then decreasing for 1 ≤ t ≤ 3 to −F and finally increasing again for 3 ≤ t ≤ 4 to return to 0 at t = 6. The computations are carried out with the time step Δt = 0.1 s. Time histories for the free end displacement Fig. 8 Vertical displacement in the free end of the cantilever beam under force: resonance components for the different configurations are given in Fig. 9. We note that the amplitude of vibration is of the same order as the total sum of the applied forces.The obtained results show a good agreement between the geometrically exact beam and 3D enhanced solid elements. The beam reference response is very close to the response of the 3D solid, with somewhat greater accuracy in the case of the incomplete mass matrix (neglecting the rotation inertia) regarding the dominant response frequencies.
Despite better accuracy of computed response, the presence of a zero matrix block associated with rotation degrees of freedom is the source of major difficulty in the computational treatment of such problem over a very long interval of time. In order to illustrate this, we increase the amplitude of the loading to F = 2 N and we perform the computation of the fields of displacement over a longer interval in time. The time step is selected as Δt = 0.1 for both configurations. The dynamic response for displacement at the free end is plotted in Fig. 10. First, we note high oscillations in the computed displacement which are obtained by the Newmark scheme. Moreover, the configuration using the standard mass matrix can no longer converge for time exceeding T = 81.3. However, with a regularized mass matrix, the residual remains sufficiently small to ensure the convergence for substantially longer period.  Adding terms in the mass matrix associated to the rotational degrees of freedom ensures a satisfying performance of the proposed element in computations over very long time interval. However, this implies additional computation operations corresponding to the angular accelerations.

Large scale wind turbine: dynamic problem
For the last example, we consider a dynamic analysis for a model of wind turbine with flexible blades. We choose to make simpler blade design, since our first objective is to check the efficiency of the proposed approach to handle numerically the structure behavior under large overall motion without getting into too much details of shape optimal design.
The turbine consists of three blades with a length of 50 m. Each blade is modeled with the proposed enhanced 3D solid. The revolute joint, providing connectivity between blades, is defined with two rigid bodies [23]. The interaction between flexible and rigid components doesn't need any special requirements given that the rotational degrees of freedom appear in the rigid body formulation for the large displacement in a non-linear form [23]. The tower, which is 85 meters tall and with circular cross section, is modeled with geometrically exact beam [24] (see Fig. 11). The wind turbine blades are subjected to a uniform distributed load using follower pressure approach (see Appendix II). The blades facets where the loading is applied are defined with 4 nodes two-dimensional elements that belong to the boundaries of the three-dimensional 8 nodes brick elements of the blades. We choose an amplitude of the follower pressure p = 360 N/m 2 , equivalent to an average wind speed 60 km/h (see Appendix II for special treatment of follower pressure loads). The pressure function p is presented in Fig. 11.
The chosen material and geometric characteristics are presented in Fig. 11. The computations are carried out with the constant time step Δt = 0.2. The number of iterations per step doesn't exceed 20, showing rather robust convergence with regard to multi-finite elements used in the same structure.
The proposed element describes properly the rotational motion of the wind turbine. For a long time period, we start observing a large deformation related to a combined bending and twisting motion which causes premature damage of the turbine. This demonstrates the usefulness of the proposed solid element with respect to the geometrically exact beam element for this kind of modeling. Namely, beam formulation is based on the rigid crosssection assumption, which can not captive such damage phenomena. The tower is submitted to large bending, involving harmful effect on the wind turbine stability which may be even more pronounced for large scale structures in offshore environments. In order to overcome this problem, we introduce pre-stressed cables connecting the wind tower to the support (see Fig. 12). The finite element approach for computing cable structure undergoing large displacement is coded for non linear analysis, based on the previous work [25].
The cross-section area of the cable used here is A = 0.018 m 2 and the mass density is ρ = 0.1 kg/m 3 . The undeformed cable configuration is straight line, with an initial pre-tension stress S 0 = 800 N/m 2 . The Saint Venant material model is adopted E = 2 × 10 7 N/m 2 . We increase the pressure amplitude to 900.
The time histories of the horizontal displacement at the top of the tower for both configurations, with or without cables, are plotted in Fig. 13. We note that tower's deflection is reduced, on average, by half when cables are added. Hence, the proposed solution combines robustness and low weight to improve the stability of the wind turbine without increasing costs.

Closing remarks
In this paper, we have discussed the regularized variational formulation for 3D solid element with the incompatible mode method to analyze geometrically nonlinear problems in statics and dynamics. The statics formulation is follow up of our previous works limited only to 2D statics [18,26]. The main novelty is the generalization of the formulation as well as the finite element implementation to framework of 3D nonlinear dynamic problems.
In particular, the mass matrix is defined from a simple quadratic form of the kinetic energy as opposed to the co-rotational formulation and the floating frame formulation used by the others authors. Moreover, a regularized form of the mass matrix is proposed in order to ensure high computational performance without adding complexity. It is important that the regularized mass matrix ensures a good performance by being able to avoid infinite frequencies, and thus facilitating the convergence of the Newmark implicit scheme. Further gain in computational efficiency, which is worth to note, is brought about by the operator split methodology reducing the computation of the the final value of the incompatible mode parameters to a single iteration.
It is interesting to note that the proposed elements with rotational degrees of freedom can be easily combined with geometrically exact beam models and can also be efficient in the nonlinear dynamic analysis of thick plates and shells [9]. This enables a smooth transition between solid and structural elements. Moreover, they are constructed to allow in plane cross-sectional changes as well as out of plane cross-sectional warping.
In order to reduce the impact of zero masses associated to the rotational degrees of freedom, besides of a simple regularization of the mass matrix, another alternative seems more relevant by controlling the dissipation of high frequency modes contribution as proposed for beam case in [19,27] .