Dynamic Nonlinear Co-rotational Formulation for Two-dimensional Continua

This paper presents an extension of BATTINI ́s formulation to two-dimensional analysis of nonlinear dynamical problems. The main interest of the co-rotational approach is to separate rigid body motions from pure deformations at the local element level through co-rotated framework. Employing co-rotated framework to derive both internal and inertial terms is the objective of the formulation presented in this study. Numerical examples are presented to illustrate the ability of the proposed co-rotational formulation.


Latin American Journal of Solids and
Decomposing the motion of the element into rigid-body and pure deformational parts is the main idea behind the co-rotational approach.This goal is attained through a local reference system that continuously rotates and translates with the element.The deformational response is expressed at the level of a local system, whereas the geometric non-linearity persuaded by the large rigid-body motion is incorporated in transformation matrices relating the local and global internal force vectors and tangent stiffness matrix (Battini, 2008;Eriksson and Faroughi, 2014).The main advantage of the co-rotational approach is the fact that different material models and different geometrical linear theory can be easily used in the local deformation when the pure deformation part is assumed to be small (Battini, 2008;Eriksson and Faroughi, 2014).Another interesting feature of the CR approach is that the transformation matrices are identical for elements which have the same number of nodes and degrees of freedom.
In the literature, many applications of 4-node plane in different fields of analysis have been investigated.Seki and Atluri (1994) developed a 2D plane element using an assumed stress hybrid element.Fores (2006) employed an assumed strain approach for triangular element based on a total Lagrangian formulation.Akasha and Mohmed (2012) developed the geometrical nonlinear formulations for a 4-node plane element based on total Lagrangian formulation.On the other hand, many researchers have applied the co-rotational (CR) technique as an alternative approach to study the geometrical nonlinear analysis of structures.This approach has been extensively used in construction of beam (Crisfiled, 1990;Battini, 2008;Le et al., 2013), plate (Izzuddien) and shells elements (Eriksson et al., 2002;Battini, 2007).The CR technique has been rarely utilized in plane structures (Mohammed andMajied, 2014), andBattani (2008) was the first who developed a nonlinear static co-rotational 4-node plane element.Regarding the nonlinear dynamic formulation of 2-dimensional 4-node plane element, to the best of author's knowledge, there is no complete numerical analyses developed using the co-rotational approach.
The main target of this paper is to develop a 2D CR nonlinear dynamic formulation for flexible plane structures.The concept of CR framework (decompose the motion of element into pure deformation and rigid-body) is adopted to obtain the internal force vector, tangent stiffness matrix and inertia force vector and tangent dynamic matrix.With a CR formulation, a linear 4-node plane element is used to model the geometrical nonlinear behaviour of plane structures in local coordinates.
Regarding the solution of nonlinear dynamic equation of motion, the Newmark method is implemented in this work.A predictor-corrector procedure is used to solve the nonlinear dynamic equation of motion.Indeed, to satisfy the nonlinear dynamic equation of motion at each time step, a modified Newton-Raphson method is employed (Bathe, 1996).

CO-ROTATIONAL FORMULATION
This section presents the dynamics of the two-dimensional co-rotational formulation for 4-noded plane element.To obtain the dynamic nonlinear co-rotational formulation four terms must be defined.The first one determines angles of rotation between co-rotating frame and global coordinate; the second one expresses the relation between local and global variables.The third term is defined to attain the internal force vector and tangent stiffness matrix, and lastly, the fourth term is considered to obtain the inertia force vector and dynamic tangent stiffness matrix.In following sections, Latin American Journal of Solids and Structures 12 (2015) 477-491 Sections 2-1 and 2-2, the internal force vector and tangent stiffness matrix developed by Battini (2008) are reviewed and reproduced for simplicity and clarity.Then, in Section 2-3, the extension of the formulation is presented for two-dimensional analysis of nonlinear dynamical problems.

Coordinate systems, 4-node plane kinematics
A plane finite element with four corner nodes is shown in Fig. 1.It is described by node coordinates in its referenced positions, and their displacements.The global coordinate system is defined by
The global coordinates and global displacements of node i are denoted by respectively.The origin of local coordinates, ( ) C, x, y is taken at center point, C, which is calculated by The axes of the local and global coordinate systems are considered parallel to the undeformed configuration.Therefore, the local coordinates of the node i is defined by ( ) r , ; ; The element motion initially up to deformed configuration is divided in two steps.The first one is the rigid translation and rotation of the initial element.The rigid translation is described by the movement of the center point of the element, C, expressed in the global coordinate, and denoted by u c that is calculated as The rigid rotation, denoted by the angle θ , describes the orientation of the local coordinate system in the current configuration.After some algebraic manipulations, the rigid rotation is obtained as follows (Battani, 2008): Latin American Journal of Solids and Structures 12 (2015) 477-491 The second step is a local deformation concerning to local coordinate which is defined by the local deformation displacements u i , i=1,...,4.These displacements are expressed as cos sin sin cos Finally, the local and global displacement vectors are

internal force vector and tangent stiffness matrix
Since the internal virtual work in both local and global coordinate system is equal, the transformation matrix between global and local coordinate systems can be obtained using the variational relation between local and global displacements.The relation between the local internal force vector and the global one f g is The global internal force vector can be specified by In the literature, there are many linear formulations to calculate the local stiffness matrix K l and internal force vector f l associated with P l .The global tangent stiffness matrix is expressed as: Taking the variation of Eq. 8, (Le et al., 2013), it leads to the global stiffness matrix as (Battani, 2008): where Latin American Journal of Solids and Structures 12 (2015) 477-491 The matrices A and G show the effects on the deformational displacements from a rotation of the reference system and they can be expressed as (Battani, 2008): ( ) where ( ) , ; 1,..., 8 i i x y i = denotes the initial coordinates of node i in the undeformed element system, and ( ) x y i = denotes the deformed coordinates of node i which is calculated as di x = ;

Inertia force vector, mass matrix and tangent dynamic matrix
In this section, the inertia force vector, mass and tangent dynamic matrices are derived.In order to formulate internal and inertia terms, the same kinematic explanation of co-rotational is implemented.
The kinetic energy of an element, KE , is obtained as where, ρ stands for the material density [kg/m 3 ], t is the thickness and G u ɺ and G vɺ denote the global velocity components.The inertial force vector is computed from kinetic energy by employing the Lagrangian equation of motion, in which the kinetic energy, KE , reads as: where M is the global mass matrix and it is given by

M T T
T l m = (16) In Eq (16), T and l m denote the rotation matrix and local mass matrix of four-node plane element.In the literature, several formulations are usually taken for the local mass matrix l m , such as lumped mass matrix, consistent mass matrix obtained from bi-linear shape functions or superconvergent mass matrix (Fried and Chavez, 2004).
The rotation matrix, T is defined as Latin American Journal of Solids and Structures 12 (2015) 477-491 The first term of the inertial force vector is computed by differentiations of the kinetic energy: In Eq (19), the term M θ denote the derivative of M with respect to θ and is calculated as follows: where: ( ) Substituting Eq.21 into Eq.20 leads to another expression for the term M θ as ( ) Further, the second term of inertial force vector can be computed by differentiation of the kinetic energy with respect to global displacement.Therefore it leads to: By substituting Eqs 18-23 into Eq.14, the inertial force vector is expressed as (24)

Nonlinear equation of motion
As shown in the above equations, the inertia force vector depends on where K f , g f and ext F are the inertial, internal and external force vectors, respectively.The Newmark method can be extended in order to solve Eq. 25.It is worth mentioning that here all matrices are not fixed (Crisfield, 1997). where is the equivalent dynamic out-of-balance forces.If all necessary information is available at time step n, one can solve Eq. 25 by using a predictor-corrector method.The term , 1 i n h + can be written as follow using Taylor series: ( ) in which K , d n denotes the dynamic tangent matrix and is given by Substituting Eq. 27 into Eq.26, an expected incremental predictor step can be computed as ( ) The use of Newmark time integration leads to the corrective updates displacement, velocity and acceleration as described in below In Eq. 30, t ∆ shows the time step and the β and γ are the parameters of the Newmark method.Using the updated displacement, velocity and acceleration (Eq.30), Eq. 26 can be again computed.If the equivalent dynamic out-of-balance forces are not zero, the Newton-Raphson corrector must be used.After some algebraic manipulations, one can obtain the following for the corrector This procedure is repeated until the equivalent dynamic out-of-balance forces are smaller than a predictor tolerance.
In order to obtain the tangent dynamic matrix, the differentiation of each term must be calculated.The tangent stiffness matrix, mass and gyroscopic matrix are given by ; The derivation of each term of inertia force vector respect to global displacement is computed in order to obtain the expression for f P K g ∂ ∂ .This leads to ( ) ) Therefore, the expression of f P where the term M θ θ ∂ ∂ is obtained using Eqs (20-22) as: with ( ) , , , ; 1 0 Latin American Journal of Solids and Structures 12 (2015) 477-491 Substituting Eqs 10, 33 and 34 into Eq.28, and after some algebraic work, the tangent dynamic matrix is found as

NUMERICAL EXAMPLES
Two numerical examples are used to show the performance of the proposed formulation.MATLAB Version 7.4 (R2007a) was used for the calculations.In presented formulation, the linear Qm6 (Taylor et al., 1976) has been used for the local formulation.The results obtained with the new corotational formulation are evaluated with results obtained by nonlinear formulation of the Qm6 element corresponding to the updated Lagrangian considered as a reference solution.
Note that all considered structures are without damping, and β and γ are fixed at ¼ and ½, respectively.

Example 1: angle frame
For the first test-case, an angle frame depicted in Fig. 2 is considered.The structure is clamped at one end, and its module of elasticity, Poisson ratio and material density are considered to be 40 GPa, 0.3 and 1740 kg/m 3 , respectively.The dimension of the structure loaded by a harmonically distributed vertical force 0 sin( ) is also shown in Fig. 2. The amplitude of the external load, 0 F , and its frequency, ω , are assumed to be 4000 N and 25 π rad/s, respectively.The size of the time step is taken 1.2E-3 s, and a fine mesh consists of 304 square elements is used.In Figs. 3 and 4, the time history response of node A in X and Y directions are illustrated, respectively.According to Figs. 3 and 4, it can be concluded that the results obtained by the proposed formulation agree very well with the selected reference solution obtained from nonlinear formulation based on the updated Lagrangian.According to Table 1, it can be concluded that the presented co-rotational formulation is faster than the reference solution (which in this study is the nonlinear formulation based on the update Lagrangian).As mentioned earlier, one of the limitations of the proposed CR method is that it is only valid for small pure deformation.To show the applicability of these formulation to example 1, where the structure is meshed by 240 elements, the element deformation (calculated from Eq. ( 5)) for a sample element (shown in Fig. 7) at three different time steps are calculated and tabulated in Table 2.
One can conclude from Table 2 that the pure deformations for the considered element are small in all time steps.This argument is also corrects for all elements of the structure, and thus the CR formulation is a robust approach to model these kind of problems.

Example 2: angle frame
As a second test-case example, a half ring depicted in Fig. 8 is considered.Due to the symmetry, only half of the structure is modeled.This structure is made of steel with the density of 7800 kg/m 3 , Poisson ratio of 0.3, and the module of elasticity is considered to be 200 GPa.The inner and outer radial of the structure are 200 and 210 mm, respectively.The structure is clamped at one end, and the other end is fixed in X direction.The structure is loaded by a harmonic concentrated vertical force 0 sin( ) 25 ω π = rad/s.Here, the size of time step is assumed to be 3.2E-3 s, and a fine mesh consists of 240 elements is considered.Similar to the first example, the time history response of node A is considered and its variation along Y direction is plotted in Fig. 9.
From Fig. 9, one can find a very good agreement between results obtained by the proposed corotational formulation and the one calculated by the reference classical method.
The same comparison as above, but using a coarse mesh for proposed model, is depicted in Fig. 10.In this example the coarse mesh consists of 160 square elements.This qualitative comparison along with the quantitative one reported in Table 3 clearly highlights the robustness of the proposed formulation.It should be noted that in this example all elements' pure deformation are small (at the same order for that of example 1).

CONCLUSION
In this study, an extension of BATTINI´s formulation for nonlinear dynamical analysis of fournode plane element was presented.In CR approach, the small deformation quantities are extracted from the given displacement field, named CR-framework.Here, the CR-framework is adopted to develop internal forces vector and tangent stiffness matrix along with the inertia forces vector and tangent dynamic matrix.The proposed formulation was successfully tested against two numerical examples that show its applicability to investigate the plane structures possessing large displacement, but small strains.Furthermore, it should be emphasized that using the described CR framework, competent linear elements are automatically transformed to nonlinear formulations.

Figure 2 :
Figure 2: Angle frame with dimension; left fine mesh, right coarse mesh.

Figure 3 :
Figure 3: The X displacement history of point A.

Figure 4 :
Figure 4: The Y displacement history of point A.

Figures 5
Figures5 and 6show the time history response of node A in X and Y directions using new corotational formulation using a coarse mesh.Here, the coarse mesh consists of 240 square elements.As shown in Figs.5 and 6, the results obtained from new co-rotational formulation with coarse mesh are much closer to the ones obtained from the reference solution (fine mesh).

Figure 5 :
Figure 5: The X displacement history of point A.

Figure 6 :
Figure 6: The Y displacement history of point A.

Figure 7 :Table 2 :
Figure 7: Angle frame with dimension; coarse mesh along with sample element.

Figure 8 :
Figure 8: quarter thin ring with dimension; left fine mesh, right coarse mesh.

Figure 9 :
Figure 9: The Y displacement history of point A.

Figure 10 :
Figure 10: The Y displacement history of point A.

Table 1
is presented the CPU time and total number of iterations necessary for two formulations.

Table 1 :
Angle frame-CPU time (total number of iterations).

Table 3 :
quarter thin ring-CPU time (total number of iterations).