Derivation of Symmetric Secant Stiffness Matrices for Nonlinear Finite Element Analysis

Nonlinear finite element analysis requires the derivation of the secant elastic stiffness matrix of the discretised mechanical system. Most approaches of the literature lead to express the secant stiffness matrix as an asymmetric matrix, which then can be made symmetric by complicated numerical calculations or analytical formulae. The paper illustrates a straightforward method for the derivation of symmetric secant (and tangent) stiffness matrices for isoparametric finite elements made of hyper-elastic materials with both geometric and material nonlinearities. The method is based on the adoption of the nodal coordinates in the current configuration, instead of the nodal displacements, as the main unknowns. The method is illustrated in general for isoparametric elements with translational degrees of freedom. Specialised expressions and numerical examples are presented for truss bar elements and structures.


INTRODUCTION
The increasing availability of light-weight and high-strength structural materials, in combination with advanced manufacturing techniques, has led to the conception of ever more enterprising and slender structures in modern aerospace, mechanical, and civil engineering applications [1][2][3]. As a matter of fact, the design of such advanced structures cannot be undertaken without resorting to numerical methods for structural analysis [4].
When a structural element undergoes displacements and deformations that overcome the limits acceptable for the validity of infinitesimal theory of elasticity, a nonlinear analysis must be conducted to accurately predict its mechanical behaviour [5][6][7]. In the finite element method (FEM), the nonlinear static equilibrium equations for a discretised, elastic mechanical system are generally stated in the following matrix form: ( ) -= K u u p 0, (1) where: K -the secant elastic stiffness matrix of the system, depending on the nodal displacement vector u, and p -the nodal load vector [8].
Equations (1) can also be written as: ( ) = f u p, (2) where: (3) is the nodal elastic force vector. Equations (2) highlight the fact that static equilibrium is a balance between the internal and external forces. Furthermore, from Equations (3), it follows that given two vectors f and u, no unique definition of the secant stiffness matrix exists. Indeed, K(u) can be expressed as either a symmetric or an asymmetric matrix. It is clear that having a symmetric matrix is desirable to reduce the use of computer memory and to improve the computational effectiveness of solution algorithms.
In linear analysis, i.e., for small strains and displacements, the stiffness matrix K is evaluated in the reference configuration of the system independently of the nodal displacements. In this case, the stiff ness matrix is always a symmetric matrix. Instead, in geometrically nonlinear analysis, i.e., when strains and displacements may attain large values, K must be evaluated in the current confi guration and depends on the nodal displacements. In the literature, several strategies have been proposed for the evaluation of the secant stiff ness matrix [9,10]. Most approaches lead to express K as an asymmetric matrix, which then can be made symmetric by suitable manipulations [11]. Also, it should be noted that most approaches of the literature are based on complicated numerical calculations or analytical formulae.
The paper illustrates a straightforward method for the derivation of symmetric secant stiff ness matrices for systems made of hyper-elastic materials characterised by both geometric and material nonlinearities. The method is based on the adoption of the nodal coordinate vector in the current confi guration x, instead of the nodal displacement vector u, as the main unknown. As a result, the nodal internal force vector is expressed as: (4) where: S -a symmetric secant stiff ness matrix depending on the nodal coordinate vector.
The proposed approach turns out to be similar to the absolute nodal coordinate formulation (ANCF) used for fl exible multibody dynamics [12], here adapted to the context of computational solid and structural mechanics. In what follows, fi rst the method is presented for general isoparametric fi nite elements with translational degrees of freedom. Then, specialised expressions are deduced for the truss bar element. Lastly, the eff ectiveness of the method is illustrated by the analysis of some well-known planar and spatial truss structures.

Finite element formulation
Let us consider the nonlinear static equilibrium problem for a continuous body made of a hyper-elastic material and occupying a region Ω of the Euclidean space of dimension d. In the finite element formulation (Figure 1), the region Ω is approximated by the union of m regions, 1 2 , , , m W W W ! , corresponding to the fi nite elements. Elements are connected to each other and to the support only at a discrete number n of points called nodes. Also, the distributed surface and body forces, as well as the kinematic restraints, of the continuous model are substituted by their equivalent, concentrated nodal loads and restraints [8].
Let us fi x a Cartesian reference system Thus, physical vectors can be represented in terms of their real-valued components as column vectors of d R (here and in the following, R denotes the set of real numbers). Besides, let us choose a reference confi guration W as one of the many confi gurations that the body can occupy. Let P ÎW be the point corresponding to a given material particle in the reference confi guration and P Î W the point corresponding to the same particle in the current confi guration. We denote with P P O =x and P P O =x the position vectors of the two points. Besides, (5) denotes the displacement vector of the considered particle.
The displacement vector fi eld defi ned over the reference confi guration represents the main unknown of the continuous problem. In the standard fi nite element formulation, the displacement fi eld within each element is approximated by a suitable interpolation of the displacement vectors of the nodes connected to the element. In particular, for an isoparametric element with only translational degrees of freedom, the displacement vector fi eld can be represented as: The nodal displacement vector of the element can also be expressed as follows: e e = u A u, (7) where: Î u u u u ! R respectively are the assembly matrix of the element and the displacement vector of the system. The assembly matrix e A has all null entries, except for the macro-rows and macro-columns corresponding to the element nodes, where identity matrices The displacement vector u collects the displacement components of all of the nodes of the model and constitutes the main unknown of the discretised problem. It should be noted that while the assembly matrix is useful for the analytical formulation, in numerical implementation it is more computationally effective to code the assembly process via direct extraction of the sub-vectors and sub-matrices corresponding to the involved nodes [13].
For what follows, it is also useful to introduce the nodal position vectors of the element in the reference and current configurations, respectively. In line with the notation used for the displacements, R respectively are the nodal position vectors of the system in the reference and current configurations.

Green-Lagrange strain measure
In the presence of large strains and displacements, the deformation can be conveniently measured by using the Green-Lagrange strain tensor [14]. For an element with dimension e d , the Cartesian components of the strain tensor turn out to be: where: s a and s b are local coordinates for the element with the indices α and β ranging from 1 to e d . By introducing Equations (5) and (6) into (8), after some manipulations, the expression for the strain components can be put in the following form: where: the symmetric matrices For what follows, it will be useful to calculate the derivative of the strains:

Nonlinear static equilibrium equations
The nonlinear static equilibrium equations can be deduced by imposing the stationarity of the total potential energy of the system [15]: where: U is the strain energy stored in the whole system and W is the virtual work done by the loads.
Such contributions can be calculated as explained in the following.
The mechanical response of hyper-elastic materials is characterised by the definition of a potential function ( ) E ab j j = , from which the stress-strain relationships are derived. In particular, the work-conjugate stress measure for the Green-Lagrange strain tensor is the second Piola-Kirchhoff stress tensor [14], whose components turn out to be S E ab ab Besides, j can be interpreted as the strainenergy density, i.e. the strain energy stored in a unit volume in the reference configuration. Thus, the strain energy stored in an element e is: (14) and the strain energy of the whole system is given by: Besides, it is assumed that the load vector can be expressed as: where λ is a scalar load multiplier and p is a reference nodal load vector. As a consequence, the virtual work of the loads turns to be: By differentiating Equation 12, the nonlinear equilibrium equations are obtained: where, by recalling also Equations (14)- (17), the nodal elastic force vector is: and the nodal load vector is: (20)

Symmetric secant elastic stiffness matrices
By applying the chain rule for differentiation and recalling Equation (7), from Equation (19), the nodal elastic force vector is obtained: where: is the symmetric secant elastic stiffness matrix of the element, expressed as a function of the nodal position vector. By substituting Equations (23) into (21) and recalling Equation (11), the nodal elastic force vector for the whole system is obtained as: where: The implementation of incremental-iterative solution methods also requires the introduction of the tangent stiffness matrix of the system, which can be calculated as follows: is the tangent stiffness matrix of element e.
In closure, it is noted that the volume integrals appearing in Equations (24) and (29) can be carried out numerically, if necessary [8].

Explicit expressions for the truss bar element
The general expressions of the secant stiffness matrices are now specialised to an isoparametric truss bar element. The formulation will be valid for both planar (d = 2) and spatial (d = 3) structures.
Let us consider a truss bar element labelled e with is fixed along the bar centreline ( Figure 2).
The displacement field is assumed to be described by linear shape functions: where the following constant matrix is introduced: The axial strain in the truss bar can be calculated from Equation (9): In passing, it is noted that the axial strain, as well as the corresponding stress 11 S , turns out to be constant because of the linearity of the assumed shape functions. This result simplifi es the integration over the volume of the element and avoids the use of numerical integration procedures. In particular, by substituting Equation (33) into (24) and performing the integration over the reference volume of the element e e e V A L = , the element secant stiff ness matrix is obtained: It is worth noting that Equations (36) and (37) are valid not only for linearly elastic materials, but also for any type of nonlinear hyper-elastic material [14].

Material laws
For the sake of illustration, two ideal materials are considered for the truss bars [14]: 1) a linearly elastic material, for which the strainenergy density is: where: E is the Young's modulus; 2) an incompressible neo-Hookean material, for which the strain-energy density is: where: / 3 E µ = is the shear modulus and In the two cases, Equation (13) for the second Piola-Kirchhoff stress specialises as follows: 11 11 S EE =

Von Mises' truss
As a fi rst example, the planar von Mises' truss shown in Figure 4 is analysed. Despite its simplicity, this structure can exhibit a wide range of instability phenomena [16]. Therefore, thanks to the availability of analytical solutions, it has been widely used for the validation of numerical   [19]. Recently, Fonseca and Gonçalves [20] have presented an analytical solution for the von Mises' truss with bars made of nonlinear hyper-elastic materials.
The , respectively corresponding to a shallow truss and a deep truss. The equilibrium paths, i.e. the curves representing the equilibrium confi gurations in the load-displacement space, have been obtained by using the arc-length method [17,18] with the admissible direction cone constraint proposed by Ligarò and Valvo [19]. Figure 5 shows the paths obtained for the shallow and deep trusses with both linearly elastic and neo-Hookean materials. For the shallow truss, the overall structural response changes only marginally because of the diff erent material laws adopted. Instead, for the deep truss, remarkable diff erences in the structural response are observed due to the diff erent material laws. For the deep truss, it is known that also bifurcation points and a secondary branch are present on the equilibrium path [16,19,20]. However, their detection and tracing are not described here, as not being related to the main scope of the present work.

Star-shaped dome
As a second example, the star-shaped dome shown in Figure 6 is considered. This three-dimensional structure has been investigated fi rst by Abatan and Holzer [21] and later used by many Authors as a benchmark test.
The following numerical values are assumed: . As for the previous example, the equilibrium path has been obtained by using the method by Ligarò and Valvo [19]. Figure 7 shows the equilibrium path of the dome traced by assuming linearly elastic and  neo-Hookean material laws. For this structure, the overall response changes only marginally because of the diff erent material laws adopted. Anyway, the eff ectiveness of the proposed formulation is demonstrated. Also, it should be noted that the complete equilibrium path of this structure features some bifurcation points and secondary branches [19]. However, their detection and tracing are omitted here, as not being related to the main scope of the present work. The equilibrium path of the dome is a curve in a 22-dimensional space (1 load multiplier + 21 displacement components of the free nodes). Figures 7a and  7b show projections of such curve onto two selected planes. The apparent self-intersections of the path in Figure 7a do not correspond to bifurcation points. In fact, the same curve projected onto a diff erent plane appears as a simple curve (Figure 7b). Lastly, it is noted that one value of vertical displacement of the central node 1 may occur at several values of the load multiplier (Figure 7a), because these correspond in fact to diff erent deformed confi gurations of the structure.

CONCLUSIONS
General analytic expressions have been deduced for the secant and tangent stiff ness matrices of isoparametric fi nite elements made of hyper-elastic materials. Such symmetric matrices are useful for the fi nite element analysis of elasticity problems characterised by both geometric and material nonlinearities. Specialised expressions of the stiff ness matrices have been presented for the truss bar element. Then, the effectiveness of the proposed formulation has been demonstrated for the nonlinear static analysis of some well-known examples of planar and spatial truss structures.
The same formulation for the stiff ness matrices can be usefully exploited also for nonlinear dynamic analysis. As a fi rst application, the dynamic deployment of a cable net has been recently investigated [22].  Specialisation to more general isoparametric elements for the analysis of two-and three-dimensional elasticity problems, as well as the extension to elements with rotational degrees of freedom, will be the subject of future developments.