Dynamic Nonlinear Co-rotational Formulationfor Two-dimensional Continua

Abstract

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.

Keywords:
Co-rotational approach; nonlinear dynamic; two-dimensional continua

1 INTRODUCTION

The 4-node plane element is considered as two-dimensional continua. Many contributions are conducted in the literature for developing accurate and competent 4-node plane elements. Most of these studies only consider linear formulations (Przemieniecki et al., 1964; Felippa, 2006; Ahmadian and Faroughi, 2011).

Co-rotational (CR) approach is a simple way to derive non-linear elements at least for problems possessing large displacements and rotations, but small strain in local coordinates. The CR approach was originally introduced by Wempner (1969); Belytschko et al. (1973). Several authors used CR approach to develop efficient beam and plate elements in nonlinear static and dynamic analysis of structures (Alsafadie et al., 2010; Battini, 2007; Behdinan et al., 1998; Le et al., 2011; Le et al., 2013). Recently, Eriksson and Faroughi (2013) developed a triangular space membrane element based on CR approach for quasi-static inflation simulation. Felippa and Haugen (2005) provided a comprehensive review of the state of the art of CR formulations.

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 and Majied, 2014), and Battani (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).

2 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, 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.

2.1 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 (O, X, Y).

Figure 1:
Kinematic of element (Battini, 2008).

The global coordinates and global displacements of node i are denoted by ri = (X_i, Y_i ) and ui = (U_i, V_i ) 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

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 ucthat 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):

The second step is a local deformation concerning to local coordinate which is defined by the local deformation displacements , i = 1, ..., 4. These displacements are expressed as

Finally, the local and global displacement vectors are

2.2 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 fg is

The global internal force vector can be specified by

In the literature, there are many linear formulations to calculate the local stiffness matrix Kl and internal force vector fl associated with Pl. 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

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 (xi, yi); i = 1, ..., 8 denotes the initial coordinates of node i in the undeformed element system, and (xdi, ydi); i = 1, ..., 8 denotes the deformed coordinates of node i which is calculated as xdi = + x_i ; y_di = + y_i .

2.3 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³], t is the thickness and and 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

In Eq (16), T and m_l 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 m_l , 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

The first term of the inertial force vector is computed by differentiations of the kinetic energy:

and the global mass matrix is only function of θ that varies with time:

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

2.3.1 Nonlinear equation of motion

As shown in the above equations, the inertia force vector depends on Pg, , and the elastic force vector depends on Pg. Therefore, one can express the non-linear governing equation of motion as

where fK, fg and Fext 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 h_{n + 1} 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 h_{i, n + 1} can be written as follow using Taylor series:

in which Kd, 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. 26can 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 δPg

At this stage, one needs to again calculate the update displacement, velocity and acceleration using

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 ∂fK/∂Pg. This leads to

Therefore, the expression of ∂fK/∂Pg is given by

where the term ∂M_θ/∂_θ is obtained using (Eqs 20-22) as:

with

Substituting Eqs 10, 33 and 34 into Eq. 28, and after some algebraic work, the tangent dynamic matrix is found as

3 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 co-rotational 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 1/4 and 1/2, respectively.

3.1 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³, respectively. The dimension of the structure loaded by a harmonically distributed vertical force F = F ₀sin(w_t ) is also shown in Fig.2. The amplitude of the external load,F ₀, and its frequency, w, 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.

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

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.

Figure 3:
The X displacement history of point A.

Figure 4:
The Y displacement history of point A.

Figures 5 and 6 show the time history response of node A in X and Y directions using new co-rotational 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:
The X displacement history of point A.

Figure 6:
The Y displacement history of point A.

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.

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

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

Number of plane elements	New formulation	Reference solution
Coarse mesh: 240	14.32 (180)	17.43 (197)
Fine mesh: 304	20.23 (185)	24.25 (199)

Table 2:
Element pure deformation for a sample element shown in Fig.7.

node	308	309	313	314
t = 0.0055 s
x	0.0866E-5	0.0031E-5	-0.0136E-5	-0.0101E-5
y	0.3084E-5	0.019E-5	-0.1215E-5	-0.0390E-5
t = 0.055 s
x	0.2058E-4	-0.2121E-4	0.0508E-4	-0.0445E-4
y	0.2173E-4	-0.0989E-4	-0.2249R-4	0.1065E-4
t = 0.3 s
x	-0.0939E-4	0.082E-4	-0.0083E-4	0.014E-4
y	-0.1976E-4	-0.0007E-4	0.1908E-4	0.0075E-4

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.

3.2 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³, 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 F= F ₀sin(w_t ) with F ₀ = 4000 N, ω = 25π rad/s.

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

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.

Figure 9:
The Y displacement history of point A.

From Fig. 9, one can find a very good agreement between results obtained by the proposed co-rotational 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.

Figure 10:
The Y displacement history of point A.

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

Number of plane elements	New formulation	Reference solution
Coarse mesh: 160	8.51 (83)	12.23 (105)
Fine mesh: 240	12.16 (86)	16.57 (107)

It should be noted that in this example all elements' pure deformation are small (at the same order for that of example 1).

4 CONCLUSION

In this study, an extension of BATTINI´s formulation for nonlinear dynamical analysis of four- node 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.

References

• Ahmadian, H., Faroughi, S., (2011). Development of super-convergent plane stress element formulation using an inverse approach. Finite Element Analysis Design 196: 173-181.
• Akasha, N. M., Mohmed A. E., (2012). Geometrically Nonlinear analysis using Plane Stress / strain Elements based on Alternative Strain Measures. Journal of Science and Technology 13: 1-12.
• Alsafadie, R., Hjiaj, M., Battini, J.M., (2010). Co-rotational mixed finite element formulation for thin-walled beams with generic cross-section. Computer Methods Applied Mechanic Engineering 199: 3197-3212.
• Bathe, K.J., (1996). Finite element procedures, Prentice-Hall Inc.
• Battini, J.M., (2007). A modified co-rotational framework for triangular shell elements. Computer Methods Applied Mechanic Engineering 196: 1905-1914.
• Battini, J.M., (2008). A non-linear co-rotational 4-node plane element. Mech. Resea. Commun. 35: 408-413.
• Battini, J. M, (2008). A rotation- free co-rotational plane beam element for nonlinear analysis. International Journal for Numerical methods in Engineering 75: 672-689.
• Behdinan, K., Stylianou, M.C, Tabarrok, B., (1998). Co-rotational dynamic analysis of flexible beams. Computer Methods Applied Mechanic Engineering 154: 151-161.
• Belytschko, T., Hsieh, B. J., (1973). Non-linear transient finite element analysis with convected co-ordinates. International Journal Numeric Method Engineering. 7: 255-271.
• Crisfiled, M.A., (1990). A consistent co-rotational formulation for non-linear, three dimensional, beam elements. Computer Methods Applied Mechanic Engineering 18: 131-150.
• Crisfield, M.A., (1997). Non-Linear Finite Element Analysis of Solids and Structures. In: Advanced topics. Vol 2. Wiley, Chischester, 455-456.
• Eriksson, A., Faroughi, S., (2013). Quasi-static inflation simulations based on co-rotational triangular space membrane elements. International Journal of Strutural Stability and Dynamics 13: 1250067 (26 pages).
• Eriksson, A., Pacoste, C., (2002). Element formulation and numerical techniques for stability problems in shells. Computer Methods Applied Mechanic Engineering 19: 3775-3810.
• Felippa, C.A., (2006). Supernatural Quad 4: a template formulation. Computer Methods Applied Mechanic Engineering 195: 5316-5342.
• Felippa, C.A., Haugen, B., (2005). A unified formulation of small-strain co-rotational finite elements: I. Theory. Computer Methods Applied Mechanic Engineering 194: 2285-2335.
• Flores F.G., (2006). A two dimensional linear assumed strain triangular element for finite deformation analysis. American Society of Mechanical Engineers, 73: 966-970.
• Fried, I., Chavez, M., (2004). Super-accurate finite element eigenvalue computation. Journal of Sound and Vibration 275: 415-22.
• Izzuddin, B.A., (2005). An enhanced co-rotational approach for large displacement analysis of plates. International Journal for Numerical methods in Engineering 64: 1350-1374.
• Le, T.N., Battini, J.M., Hjiaj, M., (2011). Efficient formulation for dynamics of co-rotational 2D beams. Computer Mechanic 48: 153-161.
• Le, T.N., Battini, J.M., Hjiaj, M., (2013). A consistent 3D corotational beam element for nonlinear dynamic analysis of flexible structures, Computer Methods in Applied Mechanics and Engineering 269: 538-565.
• Przemieniecki, J.S., Berke, L., (1964). Digital computer program for the analysis of aerospace structures by the matrix displacement method. Flight Dynamic Lab Rep. FDL 18-64.
• Rezaiee-Pajand, M., Maghoobi, M., (2014). An efficient formulation for linear and geometric non-linear membrane elements. Latin American Journal of Solids and Structures 11: 1012-1035.
• Seki, W., Atluri, S.N., (1994), Analysis of strain localization in strain softening hyper-elastic materials, using assumed stress hybrid element. Computational Mechanics Centre, Atlanta, Georgia, No.14, 549-585.
• Taylor, R., Beresford, P., Wilson, E., (1976). A non-conforming element for stress analysis. International Journal for Numerical Methods in Engineering 10: 1211-1219.
• Wempner, G., (1969). Finite elements, finite rotations and small strains of flexible shells. International Journal of Solids and Structures 5: 117-153.

Publication Dates

• Publication in this collection
  Mar 2015

History

• Received
  29 Jan 2014
• Reviewed
  08 Sept 2014
• Accepted
  10 Sept 2014