- Email: [email protected]
Large deformation analysis for anisotropic and inhomogeneous beams using exact linear static solutions
Composite Structures 72 (2006) 91–104 www.elsevier.com/locate/compstruct
Large deformation analysis for anisotropic and inhomogeneous beams using exact linear static solutions S. Agarwal, A. Chakraborty, S. Gopalakrishnan
*
Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India Available online 19 December 2004
Abstract The existing statically exact beam finite element (FE) based on the first order shear deformation theory (FSDT) is used to study the geometric nonlinear effects on static and dynamic responses in isotropic, composite and functionally graded material (FGM) beams. In this beam element, the exact solution of the static part of the governing differential equations is used to construct the interpolating polynomials for the element stiffness and mass matrix formulation. The rotary inertia is also taken into account while formulating the mass matrix. Consequently, the stiffness matrix is statically exact and the mass matrix is much more accurate compared to any other existing FEs. These two aspects make the element devoid of shear locking and an efficient instrument for analysing wave propagation problems, wherein refined mesh is essential. A total Lagrangian formulation of the linear beam element is employed for large displacement and large rotation analysis. The superconvergent property of the beam element is shown by comparing the convergence rate of this beam element to the reduced integrated FSDT beam element formed using the linear shape functions. Numerical examples deal with the static and wave propagation problems to highlight the nonlinear effects. In the static case, both Von-Karman strain tensor and Greens strain tensor is used, whereas, for the wave propagation studies only the Von-Karman strains are used. Power-law variations of the material property distribution is used for the FGM beam. Both high and low frequency pulse loading is employed to bring out the effect of nonlinearity on the transient response. 2004 Elsevier Ltd. All rights reserved. Keywords: Geometric nonlinearity; Green strain; Von-Karman strain; Wave propagation; Composite material; FGM
1. Introduction The study of nonlinear effects on static/dynamic analysis of structural elements is essential for engineers/ researchers working in the field of aerospace engineering wherein the weight saving is of paramount importance. Further, the increasing use of beams as structural components in various fields such as marine, civil and aerospace engineering has necessitated the study of their static/dynamic behavior at large amplitudes. Also, knowledge of free vibration characteristics and postbuckling behavior of beams forms a important aspect *
Corresponding author. Tel.: +91 080 22933019; fax: +91 080 23600134. E-mail address: [email protected] (S. Gopalakrishnan). 0263-8223/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2004.10.019
in assessing the structural integrity. In modern engineering design, there is increasing use of composite beams or of beams made of FGM. Analysis of beams made of such materials has to be based on shear deformation theory as its influence is significant in predicting their global characteristics. The understanding of the response of composite laminates to foreign object impact is also of fundamental importance to the design and accurate assessment of aircraft and aerospace structures. Accurate and consistent modeling of the damage caused by the conditions can only be achieved after the basic deformation response of the system has been characterized. The current efforts to extend the performance and flight envelope of high-speed aerospace vehicles have resulted in structures which may respond to the imposed
92
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
dynamic loads in a geometrically nonlinear (large deformation) fashion. This nonlinearity is due to bending-membrane coupling and gives rise to membrane stretching (in plane stresses) when out-of-plane loading is applied. The result is stiffening of the structure and reduction of dynamic response relative to that of a linear system. Also, stresses due to membrane action, usually neglected in plate flexure may cause a considerable decrease of displacements as compared with the linear solution even though the displacements are still quite small. Conversely, it may be found that for a given load the deflections may increase more rapidly than predicted by a linear solution and indeed a state may be attained wherein the load-carrying capacity decreases with continuing deformation. Linear analysis do not account for these effects and consequently may significantly overpredict the response, leading to grossly conservative designs. Without practical design tools capable of capturing the nonlinear dynamics, further improvements in vehicle performance and system design will be hampered. In aerospace, automobile and defense industries there is increasingly use of composite and graded materials. These materials provides the best combination of the properties of metals like the toughness, low density, high strength and high stiffness. In general, composite beams exhibit coupling between bending, shear, extensional and torsional motions. The coupling can arise from anisotropy, and/or the beams geometry (e.g., pretwist, initial curvature, or gross cross-sectional shape). It is well known that layered composite structures have low impact resistance. The impact resistance depends upon the ply layup, its orientation and also on the material properties of the composite layer. Therefore, it is necessary to know the high frequency behavior of the structure for designing the composite structures for impact loading. Due to its wide applications, nonlinear static and dynamic analysis of composite beams has been the subject of many papers [1–7]. The graded materials are formed by varying the percentage content of two or more materials spatially which will have desired property gradation in spatial directions. The gradation in properties reduces residual stresses and stress concentration factors. Graded materials are also used to adhere two different materials in structures subjected to different loading environments (thermal or mechanical). If the graded material is absent at the interface, then there is a high chance that debonding will occur at some extreme loading conditions. Cracks are likely to initiate at interfaces and grow into the weaker material section. These problems can be resolved by gradually varying the volume fraction of the constituents. FGMs are materials or structures in which the material properties vary with location in such a way as to optimize some function of overall FGM. FGM has gained widespread applicability as thermal-barrier struc-
tures, wear and corrosion-resistant coatings other than joining dissimilar materials [8]. The literature on such advanced materials subjected to impact loading are limited in numbers. Analyses of shear deformable plate with through-thickness material property variation in the presence of Von-Karman nonlinearity are carried out by Reddy and Chin [9], Reddy and Praveen [10], Reddy [11] and Reddy and Hsu [12]. In this paper, a Total Lagrangian formulation for beam element [13], is employed for large deformation and large rotation analysis. This formulation must be implemented using appropriate displacement interpolation functions. There exists statically exact FEs (SEFEs), which are based on the exact solution of the static part of the governing differential equations and thus posses exact stiffness matrices. As a result, single element of this kind suffices for static analysis, subjected to concentrated loading. Further, by virtue of their exact interpolating polynomials the elements enjoy refined mass matrices, which are much more accurate compared to any other existing FEs. SEFEs are developed earlier for Herman-Mindlin rod [14], first-order shear deformable (FSD) isotropic beam [15,16], FSD composite beam [17], FSD FGM beam [18] and beam with Poissons contraction [19]. The exact shape functions used in these elements are not only a function of length of beam but also depend upon cross-sectional and material properties. The degree of interpolating polynomials is dictated by the order of governing differential equation, which attributes to the super-convergent property of the elements. In this paper, the exact shape functions of FSD beam are utilized for geometric nonlinear static and dynamic analysis, where investigations are directed towards finding the difference between linear and nonlinear response and computationally more effective than the existing FSDT beam formed using the linear shape functions. The efficiency and accuracy of this super-convergent beam element is compared with FSDT beam element formed using the linear shape functions. This element is also used to study wave propagation phenomena that results due to high frequency impact loading in composite and FGM structures. There are several works available in the literature, where for the large deformation analysis of the beams with Von-Karman strain field is employed. Menzala and Zuazua [20], has proved that the one-dimensional Von-Karman system of equations describing the planar motion of a uniform prismatic beam approaches to a beam equation of Timoshenkos type. The flexural linear and nonlinear analysis of composite beams under transverse loading based on higher order shear deformation theory is studied by Chandrashekhra and Bengera [21], wherein geometric nonlinearity is incorporated in the formulation by considering Von-Karman strain. Large deflection of multilayer Timoshenko beams using Von-Karman strain–displacement relations has been carried out earlier by Sciuva and Icardi
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
[22]. Their results deal with the nonlinear flexure response of unsymmetrically laminated beams under transverse and compressive axial loads. However, there exists no work which deals with all the components of the nonlinear strain field, which is referred here as the Greens strain. Greens strain contains all the nonlinear terms upto the second order whereas Von-Karman strain contains nonlinear term only in the axial direction. Since, there are several terms that are neglected in the Von-Karman strain, it is required to investigate the effect of the Greens strain over the conventional Von-Karman strain for assesing the overall nonlinear behavior of beams. In this work, we have considered both the strain fields and their effect is demonstrated in the static examples. The paper is organized as follows. In the next section, the details of the element formulation is given along with the required matrices of nonlinear analysis. Subsequently, in Section 3, static and wave propagation studies are carried out to demonstrate the effect of nonlinearity and the high rate of convergence of the present element in the analysis of isotropic, composite and FGM beam. In Section 4, conclusions are drawn.
2. Mathematical formulation The objective of nonlinear analysis is to describe various nonlinearities and the form of basic finite element equations that are used to analyze the nonlinear response of a structural system. For this purpose, the equilibrium of the body considered must be established in the current configuration. We consider the motion of a general body in a stationary Cartesian coordinate system, and assume that the body can experience large displacements and large strains. The aim is to establish the equilibrium positions of the complete body at the discrete time points 0, Dt, 2Dt, 3Dt, . . . , where Dt is an increment in time. To develop solution, we assume that the solutions for the static and kinematic variables for all time steps from time 0 to time t, have been obtained. Then the solution process for the next required equilibrium position corresponding to time t + Dt is typical and is applied repetitively until the complete solution path has been solved for.
is a (virwhere Rt+Dt is the external virtual work, tþDt 0 tual) variation in the Cartesian components of the strain tensor in the configuration at time t + Dt referred to initial configuration, and rtþDt are the Cartesian compo0 nents of the stress tensor in the configuration at time t + Dt and measured in configuration at time t = 0. Since the stresses and strains at time t are known, for solution at time t + Dt the following incremental decompositions are used: r0tþDt ¼ rt0 þ r0 ; 0tþDt ¼ t0 þ 0 : The strain increment components can be separated into linear and nonlinear parts 0 ¼ e 0 þ g0 ; where e0 is linear and g0 is nonlinear incremental strains which are referred to the configuration at time t = 0. The constitutive relations with tensor components C0 is used to relate incremental stresses to incremental strains r0 ¼ C 0 0 : Using above equations, the equilibrium equation in T.L. formulation is obtained as Z Z Z C 0 e0 de0 dV þ rt0 dg0 dV ¼ RtþDt rt0 de0 dV : V
V
In the T.L. formulation all static and kinematic variables are referred to the initial configuration at time t = 0. In T.L. formulation, we consider the basic equation from principle of virtual displacements considering equilibrium of the body at time t + Dt Z r0tþDt d0tþDt dV ¼ RtþDt ; V
v
This relation is used to calculate an increment in the displacements, which is then used to evaluate values of displacements, strains and stresses corresponding to time t + Dt. In the current configuration, the equilibrium conditions between internal and external forces have to be satisfied whether the displacements (or strains) are large or small. Thus, if the displacements are prescribed in the usual manner by a finite number of nodal parameters a, we can obtain the necessary equilibrium equations using the virtual work principle [23]. In all the cases we can write: Z T WðaÞ ¼ B r dV f ; ð1Þ V
where W represents the sum of the external (f) and internal generalized forces dependent upon the displacement field a and B is defined from the strain definition as d ¼ B da;
2.1. Total lagrangian(T.L.) formulation
93
ð2Þ
For large displacements, the strains depend nonlinearly on displacement field and the matrix B is dependent on a. In this case B ¼ B0 þ BL ðaÞ;
ð3Þ
in which B0 is the linear infinitesimal strain displacement matrix and BL is the nonlinear strain–displacement matrix, which depends on the displacement field. The general elastic relation is given by
94
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
r ¼ Dð 0 Þ þ r0 :
ð4Þ
In Eq. (1) integrals are in fact carried out element by element and contributions to nodal equilibrium summed in usual manner. The solution of Eq. (1) is approached iteratively and in the present work Newton–Raphson method is adopted for which the relation between da and dW is required. Differentiating w of Eq. (1) with respect to da we have Z Z T T dW ¼ dB r dV þ B dr dV ¼ K T da; V
V
and using Eqs. (2)–(4) we have Z dW ¼ dBTL r dV þ K da;
ð5Þ
where Z T K¼ B DB dV ¼ K 0 þ K L ;
€ A55 ðw;xx /;x Þ ¼ 0; dw : I 0 w € d/ : I 2 / I 1 €u þ B11 u;xx D11 /xx A55 ðw;x /Þ ¼ 0; where ( )x and ð€Þ denote differentiation with respect to x and t, respectively. The associated force boundary conditions are: A11 u;x B11 /;x ¼ N x ; A55 ðw;x /Þ ¼ V x ; where stiffness coefficients are obtained as Z Z 2 ½A11 ; B11 ; D11 ¼ EðzÞ½1; z; z dA; A55 ¼ GðzÞ dA; A
V
in which K0 represents the small displacements stiffness matrix, i.e., Z K0 ¼ BT0 DB0 dV : ð6Þ V
ð7Þ
V
which is known as the large displacement matrix and contains only terms which are linear and quadratic in a. The first term in Eq. (5) can generally be written as Z dBTL r dV ¼ K r da; ð8Þ V
where Kr is a symmetric matrix dependent on stress level. This matrix is known as the initial stress matrix or the geometric matrix. Combining Eqs. (5)–(8), dW ¼ ðK 0 þ K r þ K L Þ da ¼ K T da;
€ A11 u;xx þ B11 / ¼ 0; du : I 0 €u I 1 / ;xx
B11 u;x þ D11 /x ¼ M x
V
The matrix KL is given by Z K L ¼ ðBT0 DBL þ BTL DBL þ BTL DB0 Þ dV ;
Hamiltons principle, the following differential equations of motion are obtained in terms of three degrees of freedom:
ð9Þ
with KT being the total tangential stiffness matrix. 2.2. Development of exact shape functions The statically exact shear deformable beam element is derived in [17]. For the sake of completion, this is introduced here briefly. The strain energy (S) and kinetic energy (T) for the beam are given by Z Z 1 L S¼ ðrxx zz þ rzz xx þ rxz xz Þ dA dx; 2 0 A Z Z 1 L 2 2 T ¼ qðzÞðU_ þ W_ Þ dA dx; 2 0 A where (Æ) denotes the time derivative. Here, q(z) represents density and L and A represents length and crosssectional area of the beam respectively. On applying
A
and the mass moments are Z qðzÞ½1; z; z2 dA; ½I 0 ; I 1 ; I 2 ¼ A
and Nx,Vx and Mx are respectively, the axial force, shear force and bending moment acting on the boundary nodes. The interpolation functions for displacement field are obtained by solving a system of ordinary differential equations which is a static part of the governing partial differential equations. The exact solution has the form u ¼ c 1 þ c 2 x þ c 3 x 2 ; w ¼ c4 þ c5 x þ c6 x2 þ c7 x3 ; / ¼ c8 þ c9 x þ c10 x2 : These exact solutions have 10 constants and only 6 boundary conditions (3 degrees of freedom per node) are available. Therefore, there are only 6 independent constants. The additional 4 constants are obtained in terms of other constants by substituting exact solution into differential equations. c3 ¼ c10 B11 =A11 ; c7 ¼ c10 =3; c6 ¼ c9 =2; c10 ¼ bðc8 c5 Þ=2; or c3 ¼ aðc8 c5 Þ=2;
c7 ¼ bðc8 c5 Þ=6;
where a ¼ B11 A55 =ðA11 D11 B211 Þ;
b ¼ A11 A55 =ðA11 D11 B211 Þ:
From above relations, exact solution is written as 1 u ¼ c1 þ c2 x þ aðc8 c5 Þx2 ; 2 1 1 w ¼ c4 þ c5 x þ c9 x2 þ bðc8 c5 Þx3 ; 2 6 1 / ¼ c8 þ c9 x þ bðc8 c5 Þx2 : 2
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
The above solution can be written in matrix form as b ðxÞfag; fug ¼ fu ; w ; /g ¼ ½ N fag ¼ fc1 ; c2 ; c4 ; c5 ; c8 ; c9 g; T
ð10Þ
b ðxÞ is the matrix containing functions of x. The where ½ N column vector {a} of independent constants can be expressed in terms of nodal displacements by substituting six displacement boundary conditions while evaluating Eq. (10) for x = 0 (node 1) and x = L (node 2). The relation is " # b ð0Þ N 1 1 ½G ¼ ; f^ ug ¼ ½G fag; fag ¼ ½Gf^ ug; ð11Þ b N ðLÞ where f^ ug ¼ f u1 w1 /1 u2 w2 /2 g is the nodal displacement vector for the element. Now the displacement at any point in the element can be expressed in terms of nodal displacements by substituting Eq. (11) in Eq. (10): T b ðxÞfag ¼ ½N ^ ðxÞ½Gf^ fug ¼ fu ; w ; /g ¼ ½ N ug
¼ ½N ðxÞf^ ug :
ð12Þ T
Further, [N(x)] = [Nu(x) Nw(x) N/] , where Nu(x), Nw(x) and N/ are the exact shape functions for axial, transverse and rotational degrees of freedom respectively and are given in [19]. The element formed using these shape functions will be called statiscally exact shear deformable element (SESDE). Traditionally, linear shape functions are used for transverse displacement w and rotation /, for the analysis of Timoshenko beams. To avoid shear locking, the shear stiffness matrix is reduced integrated. The behavior of this element is thoroughly discussed in [13]. This element is referred here as reduced integrated shear deformable element (RISDE). One of the main objectives of this paper is to study the convergence of SESDE as opposed to RISDE. Thus, the shape function matrix N(x) for RISDE is 2 3 1 x=L 0 0 x=L 0 0 6 7 N ðxÞ ¼ 4 0 1 x=L 0 0 x=L 0 5: 0 0 1 x=L 0 0 x=L ð13Þ
2.3. Formulation of the nonlinear stiffness matrices According to the first order shear deformation theory (FSDT), the axial and transverse displacement field is U ðx; y; z; tÞ ¼ u ðx; tÞ z/ðx; tÞ; W ðx; y; z; tÞ ¼ w ðx; tÞ; where u and w are the axial and transverse displacement of the reference plane and / is the cross-sectional rotation (see Fig. 1(a)). The general strain vector can be
95
defined in terms of infinitesimal and large displacement components ¼ þ L ; where is the linear strain part and L is the large strain part. For the static analysis two different cases of stressstrain measures are employed. 2.3.1. Green’s strain: The total Greens strain is taken as 9 8 1 2 > 2 > > U x þ ðU x þ W x Þ > > > > > 2 = < 1 2 ; ¼ Wzþ / > > > > > > 2 > > ; : ðW x /Þ U x / where the linear strain part is given by 9 8 > = < Ux > ¼ Wz > > ; : Uz þ W x and large strain matrix is given by 9 8 2 t 3 hx 0 U 2x þ W 2x > > = < h 16 1 7 x L ¼ 4 0 htz 5 ¼ U 2z þ W 2z > 2 2> hz ; : htz htx 2ðU z U x þ W z W x Þ 9 8 2 >U þ W2> 1< x 2 x = ; ¼ ð14Þ / > 2> ; : 2/ux in which 2 3 2 3 Ux Uz 6 7 6 7 hx ¼ 4 0 5; hz ¼ 4 0 5: Wx Wz From Eq. (14), L can be written as 1 L ¼ Ah; 2
ð15Þ
where A and h matrices are 2 3 u;x z/;x w;x 0 0 6 0 0 / 0 7 A¼4 5 / 0 u;x z/;x wx 2 3 0 0 ½dN u zdN / a dN w a 0 0 N / a 0 5 ; ¼4 N / a 0 ½dN u zdN / a dN w a 34 9 8 9 8 ux > > u;x z/;x > > > > > > > > > > < = < = hx wx w;x ¼ : h¼ ¼ > > uz > > hz / > > > > > > > > ; : ; : 0 wz
ð16Þ
Cauchy stress measure is used in conjunction with the Green strain.
96
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
Fig. 1. (a) The beam element and the associated degrees of freedom. (b) Depth-wise material distribution for FGM beam.
2.3.2. Von-Karman’s strain: The large strain matrix [24], is given by 9 8 u2;x þ W 2x > > = < 1 : L ¼ 0 > 2> ; : 0
Therefore the large strain–displacement matrix is given as BL ðzÞ ¼ AG:
Now, from expression of L in Eq. (15), the A matrix can be written as 2 3 u;x z/;x w;x 0 0 6 7 0 0 0 05 A¼4 2
0
0
½dN u zdN / a 6 0 ¼4 0
0
0
dN w a 0
0 0
3 0 7 0 5;
0
0
0
Second Piola–Kirchoff stress tensor is used in analysis with Von-Karman strain tensor. Now, the nodal parameters a and displacement field in terms of shape functions are given as w1
/1
u2
w2
ð17Þ
The h matrix in terms of the shape function N and nodal parameters a can be written as 2 3 dN u zdN / 6 7 dN w 6 7 h¼6 ð18Þ 7a ¼ Ga. 4 5 N / 0
A
Now, the initial stress matrix Kr has to be formed. By Eq. (8), we have Z Z K r da ¼ dBTL r dV ¼ GT dAT r dV V V Z ¼ ðGT MGÞ dV da; ð19Þ V
where the stress vector r is given as {rxx, rzz, rxz}. The M in case of Greens strain and Cauchy stress is given as 2 3 rxx 0 rxz 0 6 0 r 0 07 xx 6 7 M ¼6 7 4 rxz 0 rzz 0 5 0
/2 g;
u ¼ N u a; w ¼ N w a; / ¼ N / a:
Now we can readily form the matrix Z LZ T K¼ ðB þ BL Þ QðzÞðB þ BL Þ dA dx: 0
and h matrix is same as in Eq. (16). The total VonKarman strain is now 9 8 1 2 2 > > > > ðu U þ þ W Þ < x x = 2 ;x . ¼ 0 > > > > ; : ðW x /Þ
a ¼ f u1
The linear strain–displacement matrix is 2 3 dN u zdN / 6 7 0 B ðzÞ ¼ 4 5: dN w N /
0
rxz
0
and in case of Von-Karman strain and Second Piola– Kirchoff stress is given as 2 3 rxx 0 0 0 6 0 r 0 07 xx 6 7 M ¼6 7. 4 0 0 0 05 0
0
0
0
Therefore, the total tangent stiffness matrix is given by
Taking differential of Eq. (15) we get
K T ¼ K þ K r.
1 1 dL ¼ dAh þ A dh ¼ A dh ¼ AG da: 2 2
The stresses are related to the strains by the Q matrix given as
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
2
Q11
6 6Q 4 13 0
Q13 Q33 0
3
0
3. Numerical examples
7 0 7 5;
ð20Þ
Q55
where the expressions of Qij can be found in [25]. For FGM material they are functions of z and related to the Youngss modulus Y(z) and Poissons ratio m in the same way as defined for isotropic materials. The E(z) is varying over the depth of the beam in a continuous fashion. It is desired to be continuous, simple and should have the ability to exhibit curvature, both ‘‘concave upward’’ and ‘‘concave downward’’ [26]. Here, the power law, having the desired properties and introduced by Wakashima et al. [27] is used and is given by P ðzÞ ¼ ðP t P b Þ
z 1 þ h 2
97
n þ P b;
ð21Þ
P(z) denotes a typical material property (E, G, q). Pt and Pb denote values of the variables at topmost and bottommost layer of the beam, respectively, and n is a variable parameter, the magnitude of which denotes the curvature. 2.4. Steps followed in nonlinear analysis In the iteration based solving method of nonlinear equations, both in static and dynamic cases, following steps are carried out at every load step (static case) or time step (dynamic case): (a) Elastic linear solution is obtained as a first approximation of the displacement field, a0. (b) w0 is found using Eq. (1) (using the definitions of B and the chosen stress measure). (c) Matrix KT is established using Eqs. (6), (7) and (19). (d) Correction to the displacement field is computed as 1
Da0 ¼ ðK 0T Þ W0 and processes (b)–(d) are repeated until the error in displacement or out-of-balance force vectors are within the prescribed limit. For dynamic analysis, the load at any step is further divided into pseudo time steps if they are large. The Newmark method of time integration scheme is adopted because of its unconditional stability. The only difference that the nonlinear analysis incurs is that at each time step the effective stiffness matrix is displacement dependent and thus above procedure needs to be followed, instead of a linear equation solution procedure.
To validate the developed SESDE, a number of examples have been solved for both static and dynamic cases. The main focus is on studying the geometric nonlinear effects in both the cases, and specifically the faster rate of convergence of SESDE as compared to RISDE. 3.1. Static loading on the beam First, an isotropic cantilever beam of length 1.0 m is taken and the nonlinear effects due to transverse load at its tip is investigated. One of the important parameters governing the behavior of RISDE is the slenderness ratio (L/h) of the beam. It is well known that fully integrated shear deformable beam element constructed using linear shape functions tends to lock for higher L/h ratios. The geometrically nonlinear behavior of this element is studied here for various L/h ratios to bring out the difference in the convergence rates in the two beam models, namely the SESDE and RISDE. The nonlinear effects on the beam is investigated by considering first the Von-Karman strain tensor with second Piola– Kirchoff (PKII) stress tensor and secondly using the Greens strain tensor with the Cauchy stress tensor. 3.1.1. Analysis using von-Karman strain and PKII stress: A cantilever beam is considered where different material properties and geometries are taken for different L/h ratios. For L/h = 10, the Youngs modulus Y is taken as 20 MPa and the width is 0.001 m. Similarly, for L/h = 100, E = 2 GPa and b = 0.01 m. Finally, for L/h = 200, E = 2 GPa and b = 0.005 m. The Poissons ratio is assumed zero in all the cases. The nonlinear deflection at the tip and number of iterations required to converge to the exact solution for the SESDE based analysis for different cases of L/h ratio are observed with increasing transverse load and compared with the solution of RISDE based analysis. In the case of RISDE, selective integration is employed, where full integration is used for constructing the contribution of the bending stiffness and reduced integration is used for the contribution of the shear stiffness. The validity of the results obtained from the SESDE and the RISDE is established by comparing them with results obtained from twodimensional plane stress analysis in ANSYS (Version 5.4). The results are tabulated in Tables 1–3. The displacement variation with load for L/h = 200 is shown in Fig. 2. The results of the linear analysis are also presented in the figure to highlight the difference between linear and nonlinear responses with increasing load. The SESDE and RISDE based models beams are modeled in FEM with 100 elements and the tolerance in displacement field is set at 0.001. From the figure it can be observed that at higher loads nonlinear effects are predominant and deflection
98
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
Table 1 Comparison of convergence rate of SESDE and RISDE, L/h = 10 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.198 0.29 0.38 0.46 0.54 0.59 0.63
4 4 6 7 7 7 8
0.198 0.293 0.384 0.467 0.54 0.595 0.634
4 5 6 7 8 9 9
0.195 0.28 0.36 0.423 0.482 0.535 0.59
Table 2 Comparison of convergence rate for SESDE and RISDE, L/h = 100 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
0.5 1.0 2.0 2.5 3.0 3.5 4.0
0.099 0.194 0.35 0.394 0.43 0.46 0.48
4 6 6 8 9 10 11
0.099 0.177 0.354 0.442 0.530 0.620 0.676
4 7 7 10 14 19 35
0.098 0.190 0.348 0.410 0.452 0.510 0.543
Table 3 Comparison of convergence rate for SESDE and RISDE, L/h = 200 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
0.0625 0.0937 0.1250 0.1562 0.1875 0.21875 0.25
0.198 0.29 0.37 0.43 0.47 0.51 0.53
8 9 9 9 9 11 13
0.177 0.265 0.354 0.441 0.532 0.620 0.676
10 11 11 12 15 23 35
0.195 0.28 0.36 0.42 0.474 0.52 0.536
corroborate well with those obtained from the ANSYS two-dimensional plane stress analysis. As seen in the figure, the results obtained from RISDE are more than the actual solution while using the same number of elements in both the cases. From the tables it is observed that for L/h = 10, the number of iterations required to converge to the exact solution are almost the same for both the beam elements but as the ratio increases to 100 and 200, the SESDE based model requires less number of iterations. Also the displacement values for SESDE based model are close to those obtained from the ANSYS than the RISDE solutions.
0.8
0.7
Displacement (m)
0.6
0.5
0.4
0.3 Linear SESDE NLSESDE NLRISDE ANSYS 2D Nonlinear
0.2
0.1
0
0
0.05
0.1
0.15
0.2
0.25
Load (N)
Fig. 2. Large deflection analysis of cantilever beam with L/h = 200.
values are smaller than the linear ones owing to the increased structure stiffness due to contribution from the nonlinear terms in the strain field. The results of SESDE
3.1.2. Analysis using green’s strain and cauchy stress: In this example, a cantilever beam of 1.0 m length and 0.01 m width is taken. The beam material is isotropic with E = 200 GPa and Poissons ratio is assumed zero. The beam is analyzed for the L/h ratios 10, 100 and 200 for increasing load. The displacement and CPU time taken for convergence to solution are observed and tabulated in Tables 4–6.
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
99
Table 4 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 10 Load(N) 1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
w, FEM SESDE (m) 6
2 · 10 6 · 106 1 · 105 2 · 105 4 · 105 6 · 105 8 · 105 1 · 104 1.2 · 104 1.4 · 104 1.6 · 104 1.8 · 104 2.0 · 104
CPU time (s) 0.4 0.41 0.41 0.38 0.39 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
w, MATLAB BVP (m) 6
2 · 10 6 · 106 1 · 105 2 · 105 4 · 105 6 · 105 8 · 105 1 · 104 1.2 · 104 1.4 · 104 1.6 · 104 1.8 · 104 2.0 · 104
CPU time (s) 1.34 1.19 0.99 1.28 2.03 1.53 1.44 1.0 0.82 1.15 1.76 1.17 1.0
Table 5 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 100 Load(N) 1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
w, FEM SESDE (m) 3
1.77 · 10 3.82 · 103 5.09 · 103 7.17 · 103 9.77 · 103 1.16 · 102 1.3 · 102 1.43 · 102 1.53 · 102 1.63 · 102 1.72 · 102 1.79 · 102 1.87 · 102
CPU time (s) 1.2 1.98 2.37 2.8 3.2 3.6 3.58 4.0 3.98 4.0 4.0 4.4 4.4
w, MATLAB BVP (m) 3
1.782 · 10 3.83 · 103 5.09 · 103 7.18 · 103 9.78 · 103 1.16 · 102 1.3 · 102 1.43 · 102 1.53 · 102 1.63 · 102 1.72 · 102 1.80 · 102 1.87 · 102
CPU time (s) 2.18 3.06 2.28 2.5 3.03 2.66 2.55 3.08 2.36 2.2 3.17 2.39 3.19
Table 6 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 200 Load(N)
w, FEM SESDE (m)
CPU time (s)
w, MATLAB BVP (m)
CPU time (s)
1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
4.43 · 103 7.02 · 103 8.58 · 103 1.12 · 102 1.44 · 102 1.67 · 102 1.85 · 102 2.01 · 102 2.14 · 102 2.26 · 102 2.37 · 102 2.47 · 102 2.56 · 102
3.19 4.0 4.0 4.9 5.3 5.5 5.6 6.0 6.06 6.02 6.41 6.4 6.4
4.44 · 103 7.03 · 103 8.59 · 103 1.12 · 102 1.44 · 102 1.67 · 102 1.86 · 102 2.01 · 102 2.15 · 102 2.27 · 102 2.37 · 102 2.48 · 102 2.57 · 102
2.23 1.82 2.46 3.15 4.76 5.3 5.5 5.97 6.28 6.51 8.409 9.26 10.11
In absence of suitable FEs, that deal with the complete nonlinear strain field, the solutions from the beam analysis are compared with the Boundary Value Problem (BVP) solver bvp4c in MATLAB (Version 6.0), which solve two-point boundary value problems for ordinary differential equations (ODE). The bvp4c is based on the collocation method that provides a C1 continuous solution that is fourth order accurate uni-
formly in the interval containing the two boundary points. In this analysis, the domain is subdivided into 10 intervals and a spline function is fitted at each interval, which satisfies the governing ODE and C1 continuity at the interval boundaries. The relative tolerance is taken as 0.001. The beam is modeled in FEM with 100 SESDEs and the residual tolerance is set at 0.001.
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
2.5
x 10 –5
1.2
Normalized Force, P(t)/Pmax
From the tables, it is observed that, for all L/h ratio, the displacement values from the formulation are matching exactly with the values obtained from the MATLAB solution. Also it is observed that the CPU time for L/h = 10 is less for the SESDE based analysis as compared to the MATLAB solutions. For L/ h = 100, the CPU time for the beam element is less for smaller magnitude of loads but as the load increases the CPU time obtained from both the analysis are comparable. For L/h = 200, the CPU time taken by the beam element is more for small magnitude loads but as the load increases to higher magnitude the CPU time taken is less as compared to that taken by MATLAB solver. Thus, the static analysis concludes that the developed SESDE is quite accurate with faster convergence rate to capture the nonlinear displacement field (in comparison to the RISDE) and the element is devoid of shear locking at high L/h ratios.
2
Frequency amplitude
100
1
0.8 0.6 0.4
1.5
0.2 0 0
100 150 Time, (µ s)
1
300
0.5
0 0
10
20
30
40
50
60
70
80
90
100
Frequency, kHz
Fig. 3. Broadband pulse load.
3.2. Wave propagation analysis
3.2.1. Impact on a composite beam A cantilever beam of length 1.0 m and width 0.1 m is taken in this study. There are 10 layers of equal thickness taken as 0.001 m, bottom five plies of 0 and top five plies of 90. The composite material is GFRP, with E1 = 144.5 GPa, E2 = 9.68GPa, density, q = 1389 kg/m3, m23 = 0.3 and m12 = 0.02. The beam is modeled in FEM with 1000 elements and is subjected to various axial and transverse dynamic loads at the free end of the beam. The responses, in terms of axial and transverse velocities are observed at the free end. The residual tolerance is taken as 0.001 for all subsequent simulations. The broadband pulse considered in this study is shown in Fig. 3, which is having a frequency content of 46 kHz. The time domain representation, shown in the inset, reveals that the load has a duration of only 50 ls. Also, the inset shows the variation for normalized load P(t)/Pmax and Pmax is varied for different simulations to bring out the effect of nonlinearity. To capture
1
0.5
Normalized axial velocity
The characteristics of the wave propagation problem is that the frequency content of the forcing function (such as high velocity impact or blast loading) is very high. Hence, high frequency modes will participate in the response unlike conventional structural dynamics problems. At higher frequencies, the wavelengths are smaller requiring also the element size to be small (of the order of the wavelength). Thus, FE model for wave propagation has a large system size. The aim of the present section is to highlight the difference between nonlinear and linear response in axial and transverse velocity pattern for both composite and FGM beams. In all the subsequent plots, the linear and nonlinear responses are normalized with the maximum linear response in each case.
the high frequency content, the time step dt in Newmark integration is taken as 1.0 ls. The load is applied first at the free end of the cantilever beam in axial direction and the normalized axial velocity measured at the impact site is plotted in Fig. 4 for both linear and nonlinear analysis. The Pmax is taken as 3000 kN. In the figure, the waveform at 100 ls is the incident pulse, whereas, the second inverted pulses are the reflections from the fixed end. As the figure suggests, the incident pulses in the nonlinear analyses are of lower magnitude compared to that of linear analysis. Further, the reflected pulse arrives little early in the nonlinear case, which shows relatively high group speed in nonlinear axial mode. There is not much difference in the SESDE and RISDE analyses, as the response is due to axial loading.
0
–0.5
NL-SESDE NL-RISDE L-SESDE
–1
–1.5
0
1
2
3
4
Time, [s]
5
6
7 x 10 –4
Fig. 4. Axial velocity due to axial pulse load (composite beam).
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
0.4
Normalized transverse velocity
0.2
0
–0.2
–0.4
–0.6
–0.8 NLSESDE NLRISDE LSESDE
–1
0
1
2
3
4 Time, [s]
5
6
7 x 10 –4
Fig. 5. Transverse velocity due to axial pulse load (composite beam).
Normalized transverse velocity
1.2 NLSESDE -NL-RISDE L-SESDE
1 0.8
0.6
0.4
0.2
0 –0.2
–0.4
0
0.5
1
1.5
Time, [s]
–3
x 10
Fig. 6. Transverse velocity due to transverse pulse load (composite beam).
1 0.8 0.6
Normalized load, P/Pmax
When the same load history is applied axially at the free end with Pmax = 1000 kN, the normalized transverse velocity history measured at the impact site looks like that shown in Fig. 5. Here, large deviation in the magnitude of the incident pulse is visible, whereas, the period error (in the reflected pulse) is not so large because of relatively low group speed of bending. The analysis by SESDE and RISDE differ among each other only by a small amount and that is also at large time. Although not shown here, at later time, the RISDE based analysis fail to converge within the stipulated maximum number of iterations (here fixed at 75). The SESDE based analysis does not suffer any such complications. Finally, the broadband pulse is applied in the transverse direction at the free end of the beam (Pmax = 250 kN) and the transverse velocity measured at the free end is plotted in Fig. 6. As the figure suggests, the RISDE based analysis incorrectly predicts higher magnitude of the incident pulse than that predicted by the linear analysis, whereas, the SESDE predicts lower magnitude of the incident pulse. Further, the effect of nonlinearity can be seen in the relatively large rate of decay and early arrival of the incident pulse (captured by SESDE and overlooked by RISDE). At large time period mismatch between SESDE and RISDE is apparent. Thus, the RISDE based analysis is not accurate enough to predict the behavior when flexural response is measured. In the previous example, high magnitude of loading is required to bring out the effect of nonlinearity as the duration of the loading is very small. In this case, we take a low frequency harmonic loading which is present throughout the measured time history and expect that low magnitude of the load will be sufficient to have significant nonlinear effect. The time domain representation of the load is shown in Fig. 7. The analysis is carried out in this with a higher time step of 1 ms. The
101
0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time (sec)
Fig. 7. Low frequency loading (1 Hz).
Pmax is set at 200 N and the load is applied in the transverse direction at the free end of the same cantilever beam of previous example. The resulting normalized transverse velocity at the impact point is plotted in Fig. 8. As observed in the previous examples, in this case also the the nonlinear wave magnitude is less (almost half) than the linear wave magnitude. Further, the frequency content of the nonlinear wave is higher than that in the linear wave solution, which is due to the excitation of the higher harmonics in nonlinear response. The SESDE and RISDE analyses differ among each other at large time, indicating the mismatch in the flexural group speed. Thus, effect of nonlinearity is predominant for large duration loading and for pulse loading the effect appears only at unusually large loading magnitude.
102
S. Agarwal et al. / Composite Structures 72 (2006) 91–104 1
1
Normalized axial velocity
Normalized transverse velocity
NL-SESDE NL-RISDE L-SESDE
0.5
0
0.5
0
–0.5
–1 NL-SESDE NL-RISDE LS-ESDE –0.5
0
0.05
0.1
0.15
0.2
–1.5
0.25
0
0.1
0.2
0.3
0.4
0.6
0.7
0.8
0.9
1 –3
x 10
Fig. 9. Axial velocity due to axial impact load (FGM beam).
Fig. 8. Transverse velocity due to the low frequency loading.
Normalized transverse velocity
1
0.5
0
–0.5
NL-SESDE NL-RISDE LS-ESDE
–1
–1.5
0
1
2
3
4
5
6
7 –4
Time, [s]
x 10
Fig. 10. Transverse velocity due to axial impact load (FGM beam). 1 NL-SESDE NL-RISDE L-SESDE 0.8
Normalized transverse velocity
3.2.2. Impact on a FGM beam A cantilever beam of the same geometry of the previous study is taken in this example also with different material distribution (see Fig. 1(b)). There are three layers in the beam in the transverse direction. The top layer is that of steel of thickness 0.026 m and the bottom layer is of alumina of thickness 0.015 m. In between, there is an FGM layer of 0.009 m thickness, which blends the mechanical properties of alumina and steel smoothly according to the power law (Eq. (21)) with exponent n = 2. Properties of steel are assumed as E = 210 GPa, m = 0.3, q = 7800 kg/m3 and that of alumina are E = 390 GPa, m = 0.3 and q = 3950 kg/m3. The formulated beam is modeled in FEM code with 1000 elements and is subjected to both axial and transverse loading at the free end of the beam. The same time increment and residual tolerance of the previous example is also taken in this study. The same high frequency impact load of the previous example is applied at the free end of the beam, first axially and then in the transverse direction. For the axial case, the maximum magnitude of the load is 1000 MN. The resulting axial velocity is plotted in Fig. 9. The figure ideally captures the essential difference between the linear and nonlinear analysis, as there are both amplitude and period mismatch. Subsequently, for the same peak load applied in the axial direction, the measured transverse velocity history at the impact point is shown in Fig. 10, where the differences between the linear and nonlinear responses are apparent. Because of higher stiffness, the waves travel much faster in this case and at large time the deviation between linear and nonlinear analyses becomes significant. Finally, the load is applied in the transverse direction with Pmax = 100 MN. The resulting normalized transverse velocity history measured at the impact site is
0.5
Time, [s]
Time, [s]
0.6
0.4
0.2
0
–0.2
–0.4
0
0.1
0.2
0.3
0.4
0.5
Time, [s]
0.6
0.7
0.8
0.9
1 x 10
–3
Fig. 11. Transverse velocity due to transverse impact load (FGM beam).
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
shown in Fig. 11, which shows very little deviation in the amplitude and the group speed of bending. In all the cases, SESDE and RISDE based analyses agree quite well with each other. This enables us to conclude that for isotropic and anisotropic material with zero coupling, RISDE will be good enough to predict the response accurately, which will not be the case for asymmetrically stacked composite laminate. Further, the general effect of nonlinearity in terms of increasing stiffness, which in turn, introduces decreased amplitude and increased group speed can be seen in FGM, albeit for a large loading. Moreover, the effect of nonlinearity is not so pronounced in FGM beam in comparison to composite materials.
4. Conclusion To develop the necessary tool for large displacement and large rotation analysis of beam structures, a total Lagrangian formulation of a geometric nonlinear onedimensional beam element has been presented. The incremental displacement fields within the two-noded beam element are defined using the exact shape functions which are constructed from the exact solution of the static part of the governing differential equations of the beam element. These exact shape functions, that are function of length, cross-sectional and material properties of the beam element are used to generate exact stiffness and consistent mass matrices. This beam element is based on the FSDT. It has been shown that the exact shape function based element (SESDE) has a faster rate of convergence than the FSDT beam element formed using linear shape functions with reduced integration (RISDE). It is shown that as the length to depth ratio of the beam increases the SESDE has a much higher rate of convergence than RISDE. The SESDE based solutions are much closer to the solutions obtained from the plane stress analysis performed in ANSYS. Thus, the formulated SESDE is more accurate and computationally more efficient. The formulated beam element is used to analyze wave propagation behavior in composite and FGM material and the response is compared with linear wave solutions. It is observed that, due to the increased stiffness of the beam in nonlinear case, the magnitude of the nonlinear wave is less than that of linear wave and group velocity of the wave is also higher than the linear wave. Further, the nonlinear response shows presence of higher frequencies, which are not present in the linear response. Moreover, the effect of nonlinearity is more pronounced in the presence of nonzero loading than that is observed for pulse loading, which is present only for a very small duration. The nonlinear effect is predominant in composite beam in comparison to the FGM response. It is also ob-
103
served that to get the same amount of nonlinear effect, the load required in the case of composite beam is less than that required for FGM beam. Further, for composite materials with nonzero axial-flexural-shear coupling, RISDE based analysis will not be accurate enough and the SESDE becomes indispensable.
References [1] Rosen A, Loewy RG, Mathew MB. Nonlinear analysis of pretwisted rods using principal curvature transformation: Part I: Theoretical derivation, Part II: Numerical results. AIAA J 1987;25(3):470–8. [2] Heyliger PR, Reddy JN. A higher order beam finite element for bending and vibration problems. J Sound Vib 1998;126(2): 309–26. [3] Bauchan OA, Hong CH. Nonlinear composite beam theory. J Appl Mech 1998;55:156–63. [4] Atilgan AR, Hodges CH. Unified nonlinear analysis for nonhomogeneneous anisotropic beams with closed cross-sections. AIAA J 1990;29(11):1990–9. [5] Minguet P, Dugundji J. Experiments and analysis for composite blades under large deflections. Part I: Static behavior, Part II: Dynamic behavior. AIAA J 1990;28(9):1573–88. [6] Laulusa A, Bauchau OA, Theron NU. Theoretical and experimentalal investigation of the nonlinear behavior of composite beams. La Recherche Aerospatiale 1995;4:223–40. [7] Petrov E, Geradin M. Finite element theory for curved and twisted beams based on exact solutions for three-dimensional solids: Part I: Beam concept and geometrically exact nonlinear formulation, Part II: Anisotropic and advanced beam models. Comput Meth Appl Eng 1998;165:43–127. [8] Suresh S, Mortensen A. Fundamentals of functionally graded materials. London: IOM Communications Ltd; 1998. [9] Reddy JN, Chin CD. Thermomechanical analysis of functionally graded cylinders and plates. J Thermal Stresses 1998;26(1): 593–626. [10] Praveen GN, Reddy JN. Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates. Int J Solid Struct 1998;35(33):4457–76. [11] Reddy JN. Analysis of functionally graded plates. Int J Numer Meth Eng 2000;47:663–84. [12] Reddy JN, Hsu YS. Effects of shear deformation and anisotropy on the thermal bending of layered composite plates. J Thermal Stresses 1980;3:475–93. [13] Bathe KJ. Finite element procedures. New York: Prentice Hall; 1997. [14] Gopalakrishnan S. A deep rod finite element for structural dynamics and wave propagation problems. Int J Numer Meth Eng 2000;48:731–44. [15] Khedir AA, Reddy JN. An exact solution for the bending of thin and thick cross-ply beams. Compos Struct 1997;37:195–203. [16] Reddy JN. On locking-free shear deformable beam finite elements. Comput Meth Appl Mech Eng 1997;149:113–32. [17] Chakraborty A, Mahapatra DR, Gopalakrishnan S. Finite element analysis of free vibration and wave propagation in asymmetric composite beams with structural discontinuties. Compos Struct 2002;55(1):23–36. [18] Chakraborty A, Gopalakrishnan S, Reddy JN. A new beam finite element for the analysis of functionally graded materials. Int J Mech Sci 2003;45:519–39. [19] Chakraborty A, Gopalakrishnan S. Poissons contraction effects in a deep laminated composite beam. Mech Adv Mater Struct 2003;10(3):205–25.
104
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
[20] Menzala GP, Zuazua E. The beam equation as a limit of a 1-d nonlinear von-Karman model. Appl Math Lett 1999;12: 47–52. [21] Chandrashekhara K, Bangera KM. Linear and geometrically nonlinear analysis of composite beams under transverse loading. Compos Sci Technol 1993;47:339–47. [22] Sciuva MD, Icardi U. Large deflection of adaptive multilayered timoshenko beams. Compos Struct 1995;31:49–60. [23] Zienkiewicz OC, Taylor RL. Finite element in structural mechanics. McGraw Hill; 2000.
[24] Mukherjee A, Chaudhuri AS. Piezolaminated beams with large deformations. Int J Solids Struct 2002;39:4567–82. [25] Reddy JN. Mechanics of laminated composite plates. USA: CRC Press; 1996. [26] Markworth AJ, Ramesh KS, Parks WPJ. Modeling studies applied to functionally grsded materials. J Mater Sci 1995;30: 2183–93. [27] Wakashima K, Hirano T, Niino M. Space applications of advanced structural materials. USA: CRC Press; 1990. SP303:97.
Large deformation analysis for anisotropic and inhomogeneous beams using exact linear static solutions S. Agarwal, A. Chakraborty, S. Gopalakrishnan
*
Department of Aerospace Engineering, Indian Institute of Science, Bangalore 560012, India Available online 19 December 2004
Abstract The existing statically exact beam finite element (FE) based on the first order shear deformation theory (FSDT) is used to study the geometric nonlinear effects on static and dynamic responses in isotropic, composite and functionally graded material (FGM) beams. In this beam element, the exact solution of the static part of the governing differential equations is used to construct the interpolating polynomials for the element stiffness and mass matrix formulation. The rotary inertia is also taken into account while formulating the mass matrix. Consequently, the stiffness matrix is statically exact and the mass matrix is much more accurate compared to any other existing FEs. These two aspects make the element devoid of shear locking and an efficient instrument for analysing wave propagation problems, wherein refined mesh is essential. A total Lagrangian formulation of the linear beam element is employed for large displacement and large rotation analysis. The superconvergent property of the beam element is shown by comparing the convergence rate of this beam element to the reduced integrated FSDT beam element formed using the linear shape functions. Numerical examples deal with the static and wave propagation problems to highlight the nonlinear effects. In the static case, both Von-Karman strain tensor and Greens strain tensor is used, whereas, for the wave propagation studies only the Von-Karman strains are used. Power-law variations of the material property distribution is used for the FGM beam. Both high and low frequency pulse loading is employed to bring out the effect of nonlinearity on the transient response. 2004 Elsevier Ltd. All rights reserved. Keywords: Geometric nonlinearity; Green strain; Von-Karman strain; Wave propagation; Composite material; FGM
1. Introduction The study of nonlinear effects on static/dynamic analysis of structural elements is essential for engineers/ researchers working in the field of aerospace engineering wherein the weight saving is of paramount importance. Further, the increasing use of beams as structural components in various fields such as marine, civil and aerospace engineering has necessitated the study of their static/dynamic behavior at large amplitudes. Also, knowledge of free vibration characteristics and postbuckling behavior of beams forms a important aspect *
Corresponding author. Tel.: +91 080 22933019; fax: +91 080 23600134. E-mail address: [email protected] (S. Gopalakrishnan). 0263-8223/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2004.10.019
in assessing the structural integrity. In modern engineering design, there is increasing use of composite beams or of beams made of FGM. Analysis of beams made of such materials has to be based on shear deformation theory as its influence is significant in predicting their global characteristics. The understanding of the response of composite laminates to foreign object impact is also of fundamental importance to the design and accurate assessment of aircraft and aerospace structures. Accurate and consistent modeling of the damage caused by the conditions can only be achieved after the basic deformation response of the system has been characterized. The current efforts to extend the performance and flight envelope of high-speed aerospace vehicles have resulted in structures which may respond to the imposed
92
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
dynamic loads in a geometrically nonlinear (large deformation) fashion. This nonlinearity is due to bending-membrane coupling and gives rise to membrane stretching (in plane stresses) when out-of-plane loading is applied. The result is stiffening of the structure and reduction of dynamic response relative to that of a linear system. Also, stresses due to membrane action, usually neglected in plate flexure may cause a considerable decrease of displacements as compared with the linear solution even though the displacements are still quite small. Conversely, it may be found that for a given load the deflections may increase more rapidly than predicted by a linear solution and indeed a state may be attained wherein the load-carrying capacity decreases with continuing deformation. Linear analysis do not account for these effects and consequently may significantly overpredict the response, leading to grossly conservative designs. Without practical design tools capable of capturing the nonlinear dynamics, further improvements in vehicle performance and system design will be hampered. In aerospace, automobile and defense industries there is increasingly use of composite and graded materials. These materials provides the best combination of the properties of metals like the toughness, low density, high strength and high stiffness. In general, composite beams exhibit coupling between bending, shear, extensional and torsional motions. The coupling can arise from anisotropy, and/or the beams geometry (e.g., pretwist, initial curvature, or gross cross-sectional shape). It is well known that layered composite structures have low impact resistance. The impact resistance depends upon the ply layup, its orientation and also on the material properties of the composite layer. Therefore, it is necessary to know the high frequency behavior of the structure for designing the composite structures for impact loading. Due to its wide applications, nonlinear static and dynamic analysis of composite beams has been the subject of many papers [1–7]. The graded materials are formed by varying the percentage content of two or more materials spatially which will have desired property gradation in spatial directions. The gradation in properties reduces residual stresses and stress concentration factors. Graded materials are also used to adhere two different materials in structures subjected to different loading environments (thermal or mechanical). If the graded material is absent at the interface, then there is a high chance that debonding will occur at some extreme loading conditions. Cracks are likely to initiate at interfaces and grow into the weaker material section. These problems can be resolved by gradually varying the volume fraction of the constituents. FGMs are materials or structures in which the material properties vary with location in such a way as to optimize some function of overall FGM. FGM has gained widespread applicability as thermal-barrier struc-
tures, wear and corrosion-resistant coatings other than joining dissimilar materials [8]. The literature on such advanced materials subjected to impact loading are limited in numbers. Analyses of shear deformable plate with through-thickness material property variation in the presence of Von-Karman nonlinearity are carried out by Reddy and Chin [9], Reddy and Praveen [10], Reddy [11] and Reddy and Hsu [12]. In this paper, a Total Lagrangian formulation for beam element [13], is employed for large deformation and large rotation analysis. This formulation must be implemented using appropriate displacement interpolation functions. There exists statically exact FEs (SEFEs), which are based on the exact solution of the static part of the governing differential equations and thus posses exact stiffness matrices. As a result, single element of this kind suffices for static analysis, subjected to concentrated loading. Further, by virtue of their exact interpolating polynomials the elements enjoy refined mass matrices, which are much more accurate compared to any other existing FEs. SEFEs are developed earlier for Herman-Mindlin rod [14], first-order shear deformable (FSD) isotropic beam [15,16], FSD composite beam [17], FSD FGM beam [18] and beam with Poissons contraction [19]. The exact shape functions used in these elements are not only a function of length of beam but also depend upon cross-sectional and material properties. The degree of interpolating polynomials is dictated by the order of governing differential equation, which attributes to the super-convergent property of the elements. In this paper, the exact shape functions of FSD beam are utilized for geometric nonlinear static and dynamic analysis, where investigations are directed towards finding the difference between linear and nonlinear response and computationally more effective than the existing FSDT beam formed using the linear shape functions. The efficiency and accuracy of this super-convergent beam element is compared with FSDT beam element formed using the linear shape functions. This element is also used to study wave propagation phenomena that results due to high frequency impact loading in composite and FGM structures. There are several works available in the literature, where for the large deformation analysis of the beams with Von-Karman strain field is employed. Menzala and Zuazua [20], has proved that the one-dimensional Von-Karman system of equations describing the planar motion of a uniform prismatic beam approaches to a beam equation of Timoshenkos type. The flexural linear and nonlinear analysis of composite beams under transverse loading based on higher order shear deformation theory is studied by Chandrashekhra and Bengera [21], wherein geometric nonlinearity is incorporated in the formulation by considering Von-Karman strain. Large deflection of multilayer Timoshenko beams using Von-Karman strain–displacement relations has been carried out earlier by Sciuva and Icardi
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
[22]. Their results deal with the nonlinear flexure response of unsymmetrically laminated beams under transverse and compressive axial loads. However, there exists no work which deals with all the components of the nonlinear strain field, which is referred here as the Greens strain. Greens strain contains all the nonlinear terms upto the second order whereas Von-Karman strain contains nonlinear term only in the axial direction. Since, there are several terms that are neglected in the Von-Karman strain, it is required to investigate the effect of the Greens strain over the conventional Von-Karman strain for assesing the overall nonlinear behavior of beams. In this work, we have considered both the strain fields and their effect is demonstrated in the static examples. The paper is organized as follows. In the next section, the details of the element formulation is given along with the required matrices of nonlinear analysis. Subsequently, in Section 3, static and wave propagation studies are carried out to demonstrate the effect of nonlinearity and the high rate of convergence of the present element in the analysis of isotropic, composite and FGM beam. In Section 4, conclusions are drawn.
2. Mathematical formulation The objective of nonlinear analysis is to describe various nonlinearities and the form of basic finite element equations that are used to analyze the nonlinear response of a structural system. For this purpose, the equilibrium of the body considered must be established in the current configuration. We consider the motion of a general body in a stationary Cartesian coordinate system, and assume that the body can experience large displacements and large strains. The aim is to establish the equilibrium positions of the complete body at the discrete time points 0, Dt, 2Dt, 3Dt, . . . , where Dt is an increment in time. To develop solution, we assume that the solutions for the static and kinematic variables for all time steps from time 0 to time t, have been obtained. Then the solution process for the next required equilibrium position corresponding to time t + Dt is typical and is applied repetitively until the complete solution path has been solved for.
is a (virwhere Rt+Dt is the external virtual work, tþDt 0 tual) variation in the Cartesian components of the strain tensor in the configuration at time t + Dt referred to initial configuration, and rtþDt are the Cartesian compo0 nents of the stress tensor in the configuration at time t + Dt and measured in configuration at time t = 0. Since the stresses and strains at time t are known, for solution at time t + Dt the following incremental decompositions are used: r0tþDt ¼ rt0 þ r0 ; 0tþDt ¼ t0 þ 0 : The strain increment components can be separated into linear and nonlinear parts 0 ¼ e 0 þ g0 ; where e0 is linear and g0 is nonlinear incremental strains which are referred to the configuration at time t = 0. The constitutive relations with tensor components C0 is used to relate incremental stresses to incremental strains r0 ¼ C 0 0 : Using above equations, the equilibrium equation in T.L. formulation is obtained as Z Z Z C 0 e0 de0 dV þ rt0 dg0 dV ¼ RtþDt rt0 de0 dV : V
V
In the T.L. formulation all static and kinematic variables are referred to the initial configuration at time t = 0. In T.L. formulation, we consider the basic equation from principle of virtual displacements considering equilibrium of the body at time t + Dt Z r0tþDt d0tþDt dV ¼ RtþDt ; V
v
This relation is used to calculate an increment in the displacements, which is then used to evaluate values of displacements, strains and stresses corresponding to time t + Dt. In the current configuration, the equilibrium conditions between internal and external forces have to be satisfied whether the displacements (or strains) are large or small. Thus, if the displacements are prescribed in the usual manner by a finite number of nodal parameters a, we can obtain the necessary equilibrium equations using the virtual work principle [23]. In all the cases we can write: Z T WðaÞ ¼ B r dV f ; ð1Þ V
where W represents the sum of the external (f) and internal generalized forces dependent upon the displacement field a and B is defined from the strain definition as d ¼ B da;
2.1. Total lagrangian(T.L.) formulation
93
ð2Þ
For large displacements, the strains depend nonlinearly on displacement field and the matrix B is dependent on a. In this case B ¼ B0 þ BL ðaÞ;
ð3Þ
in which B0 is the linear infinitesimal strain displacement matrix and BL is the nonlinear strain–displacement matrix, which depends on the displacement field. The general elastic relation is given by
94
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
r ¼ Dð 0 Þ þ r0 :
ð4Þ
In Eq. (1) integrals are in fact carried out element by element and contributions to nodal equilibrium summed in usual manner. The solution of Eq. (1) is approached iteratively and in the present work Newton–Raphson method is adopted for which the relation between da and dW is required. Differentiating w of Eq. (1) with respect to da we have Z Z T T dW ¼ dB r dV þ B dr dV ¼ K T da; V
V
and using Eqs. (2)–(4) we have Z dW ¼ dBTL r dV þ K da;
ð5Þ
where Z T K¼ B DB dV ¼ K 0 þ K L ;
€ A55 ðw;xx /;x Þ ¼ 0; dw : I 0 w € d/ : I 2 / I 1 €u þ B11 u;xx D11 /xx A55 ðw;x /Þ ¼ 0; where ( )x and ð€Þ denote differentiation with respect to x and t, respectively. The associated force boundary conditions are: A11 u;x B11 /;x ¼ N x ; A55 ðw;x /Þ ¼ V x ; where stiffness coefficients are obtained as Z Z 2 ½A11 ; B11 ; D11 ¼ EðzÞ½1; z; z dA; A55 ¼ GðzÞ dA; A
V
in which K0 represents the small displacements stiffness matrix, i.e., Z K0 ¼ BT0 DB0 dV : ð6Þ V
ð7Þ
V
which is known as the large displacement matrix and contains only terms which are linear and quadratic in a. The first term in Eq. (5) can generally be written as Z dBTL r dV ¼ K r da; ð8Þ V
where Kr is a symmetric matrix dependent on stress level. This matrix is known as the initial stress matrix or the geometric matrix. Combining Eqs. (5)–(8), dW ¼ ðK 0 þ K r þ K L Þ da ¼ K T da;
€ A11 u;xx þ B11 / ¼ 0; du : I 0 €u I 1 / ;xx
B11 u;x þ D11 /x ¼ M x
V
The matrix KL is given by Z K L ¼ ðBT0 DBL þ BTL DBL þ BTL DB0 Þ dV ;
Hamiltons principle, the following differential equations of motion are obtained in terms of three degrees of freedom:
ð9Þ
with KT being the total tangential stiffness matrix. 2.2. Development of exact shape functions The statically exact shear deformable beam element is derived in [17]. For the sake of completion, this is introduced here briefly. The strain energy (S) and kinetic energy (T) for the beam are given by Z Z 1 L S¼ ðrxx zz þ rzz xx þ rxz xz Þ dA dx; 2 0 A Z Z 1 L 2 2 T ¼ qðzÞðU_ þ W_ Þ dA dx; 2 0 A where (Æ) denotes the time derivative. Here, q(z) represents density and L and A represents length and crosssectional area of the beam respectively. On applying
A
and the mass moments are Z qðzÞ½1; z; z2 dA; ½I 0 ; I 1 ; I 2 ¼ A
and Nx,Vx and Mx are respectively, the axial force, shear force and bending moment acting on the boundary nodes. The interpolation functions for displacement field are obtained by solving a system of ordinary differential equations which is a static part of the governing partial differential equations. The exact solution has the form u ¼ c 1 þ c 2 x þ c 3 x 2 ; w ¼ c4 þ c5 x þ c6 x2 þ c7 x3 ; / ¼ c8 þ c9 x þ c10 x2 : These exact solutions have 10 constants and only 6 boundary conditions (3 degrees of freedom per node) are available. Therefore, there are only 6 independent constants. The additional 4 constants are obtained in terms of other constants by substituting exact solution into differential equations. c3 ¼ c10 B11 =A11 ; c7 ¼ c10 =3; c6 ¼ c9 =2; c10 ¼ bðc8 c5 Þ=2; or c3 ¼ aðc8 c5 Þ=2;
c7 ¼ bðc8 c5 Þ=6;
where a ¼ B11 A55 =ðA11 D11 B211 Þ;
b ¼ A11 A55 =ðA11 D11 B211 Þ:
From above relations, exact solution is written as 1 u ¼ c1 þ c2 x þ aðc8 c5 Þx2 ; 2 1 1 w ¼ c4 þ c5 x þ c9 x2 þ bðc8 c5 Þx3 ; 2 6 1 / ¼ c8 þ c9 x þ bðc8 c5 Þx2 : 2
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
The above solution can be written in matrix form as b ðxÞfag; fug ¼ fu ; w ; /g ¼ ½ N fag ¼ fc1 ; c2 ; c4 ; c5 ; c8 ; c9 g; T
ð10Þ
b ðxÞ is the matrix containing functions of x. The where ½ N column vector {a} of independent constants can be expressed in terms of nodal displacements by substituting six displacement boundary conditions while evaluating Eq. (10) for x = 0 (node 1) and x = L (node 2). The relation is " # b ð0Þ N 1 1 ½G ¼ ; f^ ug ¼ ½G fag; fag ¼ ½Gf^ ug; ð11Þ b N ðLÞ where f^ ug ¼ f u1 w1 /1 u2 w2 /2 g is the nodal displacement vector for the element. Now the displacement at any point in the element can be expressed in terms of nodal displacements by substituting Eq. (11) in Eq. (10): T b ðxÞfag ¼ ½N ^ ðxÞ½Gf^ fug ¼ fu ; w ; /g ¼ ½ N ug
¼ ½N ðxÞf^ ug :
ð12Þ T
Further, [N(x)] = [Nu(x) Nw(x) N/] , where Nu(x), Nw(x) and N/ are the exact shape functions for axial, transverse and rotational degrees of freedom respectively and are given in [19]. The element formed using these shape functions will be called statiscally exact shear deformable element (SESDE). Traditionally, linear shape functions are used for transverse displacement w and rotation /, for the analysis of Timoshenko beams. To avoid shear locking, the shear stiffness matrix is reduced integrated. The behavior of this element is thoroughly discussed in [13]. This element is referred here as reduced integrated shear deformable element (RISDE). One of the main objectives of this paper is to study the convergence of SESDE as opposed to RISDE. Thus, the shape function matrix N(x) for RISDE is 2 3 1 x=L 0 0 x=L 0 0 6 7 N ðxÞ ¼ 4 0 1 x=L 0 0 x=L 0 5: 0 0 1 x=L 0 0 x=L ð13Þ
2.3. Formulation of the nonlinear stiffness matrices According to the first order shear deformation theory (FSDT), the axial and transverse displacement field is U ðx; y; z; tÞ ¼ u ðx; tÞ z/ðx; tÞ; W ðx; y; z; tÞ ¼ w ðx; tÞ; where u and w are the axial and transverse displacement of the reference plane and / is the cross-sectional rotation (see Fig. 1(a)). The general strain vector can be
95
defined in terms of infinitesimal and large displacement components ¼ þ L ; where is the linear strain part and L is the large strain part. For the static analysis two different cases of stressstrain measures are employed. 2.3.1. Green’s strain: The total Greens strain is taken as 9 8 1 2 > 2 > > U x þ ðU x þ W x Þ > > > > > 2 = < 1 2 ; ¼ Wzþ / > > > > > > 2 > > ; : ðW x /Þ U x / where the linear strain part is given by 9 8 > = < Ux > ¼ Wz > > ; : Uz þ W x and large strain matrix is given by 9 8 2 t 3 hx 0 U 2x þ W 2x > > = < h 16 1 7 x L ¼ 4 0 htz 5 ¼ U 2z þ W 2z > 2 2> hz ; : htz htx 2ðU z U x þ W z W x Þ 9 8 2 >U þ W2> 1< x 2 x = ; ¼ ð14Þ / > 2> ; : 2/ux in which 2 3 2 3 Ux Uz 6 7 6 7 hx ¼ 4 0 5; hz ¼ 4 0 5: Wx Wz From Eq. (14), L can be written as 1 L ¼ Ah; 2
ð15Þ
where A and h matrices are 2 3 u;x z/;x w;x 0 0 6 0 0 / 0 7 A¼4 5 / 0 u;x z/;x wx 2 3 0 0 ½dN u zdN / a dN w a 0 0 N / a 0 5 ; ¼4 N / a 0 ½dN u zdN / a dN w a 34 9 8 9 8 ux > > u;x z/;x > > > > > > > > > > < = < = hx wx w;x ¼ : h¼ ¼ > > uz > > hz / > > > > > > > > ; : ; : 0 wz
ð16Þ
Cauchy stress measure is used in conjunction with the Green strain.
96
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
Fig. 1. (a) The beam element and the associated degrees of freedom. (b) Depth-wise material distribution for FGM beam.
2.3.2. Von-Karman’s strain: The large strain matrix [24], is given by 9 8 u2;x þ W 2x > > = < 1 : L ¼ 0 > 2> ; : 0
Therefore the large strain–displacement matrix is given as BL ðzÞ ¼ AG:
Now, from expression of L in Eq. (15), the A matrix can be written as 2 3 u;x z/;x w;x 0 0 6 7 0 0 0 05 A¼4 2
0
0
½dN u zdN / a 6 0 ¼4 0
0
0
dN w a 0
0 0
3 0 7 0 5;
0
0
0
Second Piola–Kirchoff stress tensor is used in analysis with Von-Karman strain tensor. Now, the nodal parameters a and displacement field in terms of shape functions are given as w1
/1
u2
w2
ð17Þ
The h matrix in terms of the shape function N and nodal parameters a can be written as 2 3 dN u zdN / 6 7 dN w 6 7 h¼6 ð18Þ 7a ¼ Ga. 4 5 N / 0
A
Now, the initial stress matrix Kr has to be formed. By Eq. (8), we have Z Z K r da ¼ dBTL r dV ¼ GT dAT r dV V V Z ¼ ðGT MGÞ dV da; ð19Þ V
where the stress vector r is given as {rxx, rzz, rxz}. The M in case of Greens strain and Cauchy stress is given as 2 3 rxx 0 rxz 0 6 0 r 0 07 xx 6 7 M ¼6 7 4 rxz 0 rzz 0 5 0
/2 g;
u ¼ N u a; w ¼ N w a; / ¼ N / a:
Now we can readily form the matrix Z LZ T K¼ ðB þ BL Þ QðzÞðB þ BL Þ dA dx: 0
and h matrix is same as in Eq. (16). The total VonKarman strain is now 9 8 1 2 2 > > > > ðu U þ þ W Þ < x x = 2 ;x . ¼ 0 > > > > ; : ðW x /Þ
a ¼ f u1
The linear strain–displacement matrix is 2 3 dN u zdN / 6 7 0 B ðzÞ ¼ 4 5: dN w N /
0
rxz
0
and in case of Von-Karman strain and Second Piola– Kirchoff stress is given as 2 3 rxx 0 0 0 6 0 r 0 07 xx 6 7 M ¼6 7. 4 0 0 0 05 0
0
0
0
Therefore, the total tangent stiffness matrix is given by
Taking differential of Eq. (15) we get
K T ¼ K þ K r.
1 1 dL ¼ dAh þ A dh ¼ A dh ¼ AG da: 2 2
The stresses are related to the strains by the Q matrix given as
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
2
Q11
6 6Q 4 13 0
Q13 Q33 0
3
0
3. Numerical examples
7 0 7 5;
ð20Þ
Q55
where the expressions of Qij can be found in [25]. For FGM material they are functions of z and related to the Youngss modulus Y(z) and Poissons ratio m in the same way as defined for isotropic materials. The E(z) is varying over the depth of the beam in a continuous fashion. It is desired to be continuous, simple and should have the ability to exhibit curvature, both ‘‘concave upward’’ and ‘‘concave downward’’ [26]. Here, the power law, having the desired properties and introduced by Wakashima et al. [27] is used and is given by P ðzÞ ¼ ðP t P b Þ
z 1 þ h 2
97
n þ P b;
ð21Þ
P(z) denotes a typical material property (E, G, q). Pt and Pb denote values of the variables at topmost and bottommost layer of the beam, respectively, and n is a variable parameter, the magnitude of which denotes the curvature. 2.4. Steps followed in nonlinear analysis In the iteration based solving method of nonlinear equations, both in static and dynamic cases, following steps are carried out at every load step (static case) or time step (dynamic case): (a) Elastic linear solution is obtained as a first approximation of the displacement field, a0. (b) w0 is found using Eq. (1) (using the definitions of B and the chosen stress measure). (c) Matrix KT is established using Eqs. (6), (7) and (19). (d) Correction to the displacement field is computed as 1
Da0 ¼ ðK 0T Þ W0 and processes (b)–(d) are repeated until the error in displacement or out-of-balance force vectors are within the prescribed limit. For dynamic analysis, the load at any step is further divided into pseudo time steps if they are large. The Newmark method of time integration scheme is adopted because of its unconditional stability. The only difference that the nonlinear analysis incurs is that at each time step the effective stiffness matrix is displacement dependent and thus above procedure needs to be followed, instead of a linear equation solution procedure.
To validate the developed SESDE, a number of examples have been solved for both static and dynamic cases. The main focus is on studying the geometric nonlinear effects in both the cases, and specifically the faster rate of convergence of SESDE as compared to RISDE. 3.1. Static loading on the beam First, an isotropic cantilever beam of length 1.0 m is taken and the nonlinear effects due to transverse load at its tip is investigated. One of the important parameters governing the behavior of RISDE is the slenderness ratio (L/h) of the beam. It is well known that fully integrated shear deformable beam element constructed using linear shape functions tends to lock for higher L/h ratios. The geometrically nonlinear behavior of this element is studied here for various L/h ratios to bring out the difference in the convergence rates in the two beam models, namely the SESDE and RISDE. The nonlinear effects on the beam is investigated by considering first the Von-Karman strain tensor with second Piola– Kirchoff (PKII) stress tensor and secondly using the Greens strain tensor with the Cauchy stress tensor. 3.1.1. Analysis using von-Karman strain and PKII stress: A cantilever beam is considered where different material properties and geometries are taken for different L/h ratios. For L/h = 10, the Youngs modulus Y is taken as 20 MPa and the width is 0.001 m. Similarly, for L/h = 100, E = 2 GPa and b = 0.01 m. Finally, for L/h = 200, E = 2 GPa and b = 0.005 m. The Poissons ratio is assumed zero in all the cases. The nonlinear deflection at the tip and number of iterations required to converge to the exact solution for the SESDE based analysis for different cases of L/h ratio are observed with increasing transverse load and compared with the solution of RISDE based analysis. In the case of RISDE, selective integration is employed, where full integration is used for constructing the contribution of the bending stiffness and reduced integration is used for the contribution of the shear stiffness. The validity of the results obtained from the SESDE and the RISDE is established by comparing them with results obtained from twodimensional plane stress analysis in ANSYS (Version 5.4). The results are tabulated in Tables 1–3. The displacement variation with load for L/h = 200 is shown in Fig. 2. The results of the linear analysis are also presented in the figure to highlight the difference between linear and nonlinear responses with increasing load. The SESDE and RISDE based models beams are modeled in FEM with 100 elements and the tolerance in displacement field is set at 0.001. From the figure it can be observed that at higher loads nonlinear effects are predominant and deflection
98
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
Table 1 Comparison of convergence rate of SESDE and RISDE, L/h = 10 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
1.0 1.5 2.0 2.5 3.0 3.5 4.0
0.198 0.29 0.38 0.46 0.54 0.59 0.63
4 4 6 7 7 7 8
0.198 0.293 0.384 0.467 0.54 0.595 0.634
4 5 6 7 8 9 9
0.195 0.28 0.36 0.423 0.482 0.535 0.59
Table 2 Comparison of convergence rate for SESDE and RISDE, L/h = 100 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
0.5 1.0 2.0 2.5 3.0 3.5 4.0
0.099 0.194 0.35 0.394 0.43 0.46 0.48
4 6 6 8 9 10 11
0.099 0.177 0.354 0.442 0.530 0.620 0.676
4 7 7 10 14 19 35
0.098 0.190 0.348 0.410 0.452 0.510 0.543
Table 3 Comparison of convergence rate for SESDE and RISDE, L/h = 200 Load(N)
w, FEM SESDE (m)
No. of iterations
w, FEM RISDE (m)
No. of iterations
w, ANSYS nonlinear (m)
0.0625 0.0937 0.1250 0.1562 0.1875 0.21875 0.25
0.198 0.29 0.37 0.43 0.47 0.51 0.53
8 9 9 9 9 11 13
0.177 0.265 0.354 0.441 0.532 0.620 0.676
10 11 11 12 15 23 35
0.195 0.28 0.36 0.42 0.474 0.52 0.536
corroborate well with those obtained from the ANSYS two-dimensional plane stress analysis. As seen in the figure, the results obtained from RISDE are more than the actual solution while using the same number of elements in both the cases. From the tables it is observed that for L/h = 10, the number of iterations required to converge to the exact solution are almost the same for both the beam elements but as the ratio increases to 100 and 200, the SESDE based model requires less number of iterations. Also the displacement values for SESDE based model are close to those obtained from the ANSYS than the RISDE solutions.
0.8
0.7
Displacement (m)
0.6
0.5
0.4
0.3 Linear SESDE NLSESDE NLRISDE ANSYS 2D Nonlinear
0.2
0.1
0
0
0.05
0.1
0.15
0.2
0.25
Load (N)
Fig. 2. Large deflection analysis of cantilever beam with L/h = 200.
values are smaller than the linear ones owing to the increased structure stiffness due to contribution from the nonlinear terms in the strain field. The results of SESDE
3.1.2. Analysis using green’s strain and cauchy stress: In this example, a cantilever beam of 1.0 m length and 0.01 m width is taken. The beam material is isotropic with E = 200 GPa and Poissons ratio is assumed zero. The beam is analyzed for the L/h ratios 10, 100 and 200 for increasing load. The displacement and CPU time taken for convergence to solution are observed and tabulated in Tables 4–6.
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
99
Table 4 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 10 Load(N) 1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
w, FEM SESDE (m) 6
2 · 10 6 · 106 1 · 105 2 · 105 4 · 105 6 · 105 8 · 105 1 · 104 1.2 · 104 1.4 · 104 1.6 · 104 1.8 · 104 2.0 · 104
CPU time (s) 0.4 0.41 0.41 0.38 0.39 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
w, MATLAB BVP (m) 6
2 · 10 6 · 106 1 · 105 2 · 105 4 · 105 6 · 105 8 · 105 1 · 104 1.2 · 104 1.4 · 104 1.6 · 104 1.8 · 104 2.0 · 104
CPU time (s) 1.34 1.19 0.99 1.28 2.03 1.53 1.44 1.0 0.82 1.15 1.76 1.17 1.0
Table 5 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 100 Load(N) 1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
w, FEM SESDE (m) 3
1.77 · 10 3.82 · 103 5.09 · 103 7.17 · 103 9.77 · 103 1.16 · 102 1.3 · 102 1.43 · 102 1.53 · 102 1.63 · 102 1.72 · 102 1.79 · 102 1.87 · 102
CPU time (s) 1.2 1.98 2.37 2.8 3.2 3.6 3.58 4.0 3.98 4.0 4.0 4.4 4.4
w, MATLAB BVP (m) 3
1.782 · 10 3.83 · 103 5.09 · 103 7.18 · 103 9.78 · 103 1.16 · 102 1.3 · 102 1.43 · 102 1.53 · 102 1.63 · 102 1.72 · 102 1.80 · 102 1.87 · 102
CPU time (s) 2.18 3.06 2.28 2.5 3.03 2.66 2.55 3.08 2.36 2.2 3.17 2.39 3.19
Table 6 Comparison of convergence rate for SESDE and MATLAB boundary value problem solver, L/h = 200 Load(N)
w, FEM SESDE (m)
CPU time (s)
w, MATLAB BVP (m)
CPU time (s)
1.0 3.0 5.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0
4.43 · 103 7.02 · 103 8.58 · 103 1.12 · 102 1.44 · 102 1.67 · 102 1.85 · 102 2.01 · 102 2.14 · 102 2.26 · 102 2.37 · 102 2.47 · 102 2.56 · 102
3.19 4.0 4.0 4.9 5.3 5.5 5.6 6.0 6.06 6.02 6.41 6.4 6.4
4.44 · 103 7.03 · 103 8.59 · 103 1.12 · 102 1.44 · 102 1.67 · 102 1.86 · 102 2.01 · 102 2.15 · 102 2.27 · 102 2.37 · 102 2.48 · 102 2.57 · 102
2.23 1.82 2.46 3.15 4.76 5.3 5.5 5.97 6.28 6.51 8.409 9.26 10.11
In absence of suitable FEs, that deal with the complete nonlinear strain field, the solutions from the beam analysis are compared with the Boundary Value Problem (BVP) solver bvp4c in MATLAB (Version 6.0), which solve two-point boundary value problems for ordinary differential equations (ODE). The bvp4c is based on the collocation method that provides a C1 continuous solution that is fourth order accurate uni-
formly in the interval containing the two boundary points. In this analysis, the domain is subdivided into 10 intervals and a spline function is fitted at each interval, which satisfies the governing ODE and C1 continuity at the interval boundaries. The relative tolerance is taken as 0.001. The beam is modeled in FEM with 100 SESDEs and the residual tolerance is set at 0.001.
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
2.5
x 10 –5
1.2
Normalized Force, P(t)/Pmax
From the tables, it is observed that, for all L/h ratio, the displacement values from the formulation are matching exactly with the values obtained from the MATLAB solution. Also it is observed that the CPU time for L/h = 10 is less for the SESDE based analysis as compared to the MATLAB solutions. For L/ h = 100, the CPU time for the beam element is less for smaller magnitude of loads but as the load increases the CPU time obtained from both the analysis are comparable. For L/h = 200, the CPU time taken by the beam element is more for small magnitude loads but as the load increases to higher magnitude the CPU time taken is less as compared to that taken by MATLAB solver. Thus, the static analysis concludes that the developed SESDE is quite accurate with faster convergence rate to capture the nonlinear displacement field (in comparison to the RISDE) and the element is devoid of shear locking at high L/h ratios.
2
Frequency amplitude
100
1
0.8 0.6 0.4
1.5
0.2 0 0
100 150 Time, (µ s)
1
300
0.5
0 0
10
20
30
40
50
60
70
80
90
100
Frequency, kHz
Fig. 3. Broadband pulse load.
3.2. Wave propagation analysis
3.2.1. Impact on a composite beam A cantilever beam of length 1.0 m and width 0.1 m is taken in this study. There are 10 layers of equal thickness taken as 0.001 m, bottom five plies of 0 and top five plies of 90. The composite material is GFRP, with E1 = 144.5 GPa, E2 = 9.68GPa, density, q = 1389 kg/m3, m23 = 0.3 and m12 = 0.02. The beam is modeled in FEM with 1000 elements and is subjected to various axial and transverse dynamic loads at the free end of the beam. The responses, in terms of axial and transverse velocities are observed at the free end. The residual tolerance is taken as 0.001 for all subsequent simulations. The broadband pulse considered in this study is shown in Fig. 3, which is having a frequency content of 46 kHz. The time domain representation, shown in the inset, reveals that the load has a duration of only 50 ls. Also, the inset shows the variation for normalized load P(t)/Pmax and Pmax is varied for different simulations to bring out the effect of nonlinearity. To capture
1
0.5
Normalized axial velocity
The characteristics of the wave propagation problem is that the frequency content of the forcing function (such as high velocity impact or blast loading) is very high. Hence, high frequency modes will participate in the response unlike conventional structural dynamics problems. At higher frequencies, the wavelengths are smaller requiring also the element size to be small (of the order of the wavelength). Thus, FE model for wave propagation has a large system size. The aim of the present section is to highlight the difference between nonlinear and linear response in axial and transverse velocity pattern for both composite and FGM beams. In all the subsequent plots, the linear and nonlinear responses are normalized with the maximum linear response in each case.
the high frequency content, the time step dt in Newmark integration is taken as 1.0 ls. The load is applied first at the free end of the cantilever beam in axial direction and the normalized axial velocity measured at the impact site is plotted in Fig. 4 for both linear and nonlinear analysis. The Pmax is taken as 3000 kN. In the figure, the waveform at 100 ls is the incident pulse, whereas, the second inverted pulses are the reflections from the fixed end. As the figure suggests, the incident pulses in the nonlinear analyses are of lower magnitude compared to that of linear analysis. Further, the reflected pulse arrives little early in the nonlinear case, which shows relatively high group speed in nonlinear axial mode. There is not much difference in the SESDE and RISDE analyses, as the response is due to axial loading.
0
–0.5
NL-SESDE NL-RISDE L-SESDE
–1
–1.5
0
1
2
3
4
Time, [s]
5
6
7 x 10 –4
Fig. 4. Axial velocity due to axial pulse load (composite beam).
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
0.4
Normalized transverse velocity
0.2
0
–0.2
–0.4
–0.6
–0.8 NLSESDE NLRISDE LSESDE
–1
0
1
2
3
4 Time, [s]
5
6
7 x 10 –4
Fig. 5. Transverse velocity due to axial pulse load (composite beam).
Normalized transverse velocity
1.2 NLSESDE -NL-RISDE L-SESDE
1 0.8
0.6
0.4
0.2
0 –0.2
–0.4
0
0.5
1
1.5
Time, [s]
–3
x 10
Fig. 6. Transverse velocity due to transverse pulse load (composite beam).
1 0.8 0.6
Normalized load, P/Pmax
When the same load history is applied axially at the free end with Pmax = 1000 kN, the normalized transverse velocity history measured at the impact site looks like that shown in Fig. 5. Here, large deviation in the magnitude of the incident pulse is visible, whereas, the period error (in the reflected pulse) is not so large because of relatively low group speed of bending. The analysis by SESDE and RISDE differ among each other only by a small amount and that is also at large time. Although not shown here, at later time, the RISDE based analysis fail to converge within the stipulated maximum number of iterations (here fixed at 75). The SESDE based analysis does not suffer any such complications. Finally, the broadband pulse is applied in the transverse direction at the free end of the beam (Pmax = 250 kN) and the transverse velocity measured at the free end is plotted in Fig. 6. As the figure suggests, the RISDE based analysis incorrectly predicts higher magnitude of the incident pulse than that predicted by the linear analysis, whereas, the SESDE predicts lower magnitude of the incident pulse. Further, the effect of nonlinearity can be seen in the relatively large rate of decay and early arrival of the incident pulse (captured by SESDE and overlooked by RISDE). At large time period mismatch between SESDE and RISDE is apparent. Thus, the RISDE based analysis is not accurate enough to predict the behavior when flexural response is measured. In the previous example, high magnitude of loading is required to bring out the effect of nonlinearity as the duration of the loading is very small. In this case, we take a low frequency harmonic loading which is present throughout the measured time history and expect that low magnitude of the load will be sufficient to have significant nonlinear effect. The time domain representation of the load is shown in Fig. 7. The analysis is carried out in this with a higher time step of 1 ms. The
101
0.4 0.2 0 –0.2 –0.4 –0.6 –0.8 –1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Time (sec)
Fig. 7. Low frequency loading (1 Hz).
Pmax is set at 200 N and the load is applied in the transverse direction at the free end of the same cantilever beam of previous example. The resulting normalized transverse velocity at the impact point is plotted in Fig. 8. As observed in the previous examples, in this case also the the nonlinear wave magnitude is less (almost half) than the linear wave magnitude. Further, the frequency content of the nonlinear wave is higher than that in the linear wave solution, which is due to the excitation of the higher harmonics in nonlinear response. The SESDE and RISDE analyses differ among each other at large time, indicating the mismatch in the flexural group speed. Thus, effect of nonlinearity is predominant for large duration loading and for pulse loading the effect appears only at unusually large loading magnitude.
102
S. Agarwal et al. / Composite Structures 72 (2006) 91–104 1
1
Normalized axial velocity
Normalized transverse velocity
NL-SESDE NL-RISDE L-SESDE
0.5
0
0.5
0
–0.5
–1 NL-SESDE NL-RISDE LS-ESDE –0.5
0
0.05
0.1
0.15
0.2
–1.5
0.25
0
0.1
0.2
0.3
0.4
0.6
0.7
0.8
0.9
1 –3
x 10
Fig. 9. Axial velocity due to axial impact load (FGM beam).
Fig. 8. Transverse velocity due to the low frequency loading.
Normalized transverse velocity
1
0.5
0
–0.5
NL-SESDE NL-RISDE LS-ESDE
–1
–1.5
0
1
2
3
4
5
6
7 –4
Time, [s]
x 10
Fig. 10. Transverse velocity due to axial impact load (FGM beam). 1 NL-SESDE NL-RISDE L-SESDE 0.8
Normalized transverse velocity
3.2.2. Impact on a FGM beam A cantilever beam of the same geometry of the previous study is taken in this example also with different material distribution (see Fig. 1(b)). There are three layers in the beam in the transverse direction. The top layer is that of steel of thickness 0.026 m and the bottom layer is of alumina of thickness 0.015 m. In between, there is an FGM layer of 0.009 m thickness, which blends the mechanical properties of alumina and steel smoothly according to the power law (Eq. (21)) with exponent n = 2. Properties of steel are assumed as E = 210 GPa, m = 0.3, q = 7800 kg/m3 and that of alumina are E = 390 GPa, m = 0.3 and q = 3950 kg/m3. The formulated beam is modeled in FEM code with 1000 elements and is subjected to both axial and transverse loading at the free end of the beam. The same time increment and residual tolerance of the previous example is also taken in this study. The same high frequency impact load of the previous example is applied at the free end of the beam, first axially and then in the transverse direction. For the axial case, the maximum magnitude of the load is 1000 MN. The resulting axial velocity is plotted in Fig. 9. The figure ideally captures the essential difference between the linear and nonlinear analysis, as there are both amplitude and period mismatch. Subsequently, for the same peak load applied in the axial direction, the measured transverse velocity history at the impact point is shown in Fig. 10, where the differences between the linear and nonlinear responses are apparent. Because of higher stiffness, the waves travel much faster in this case and at large time the deviation between linear and nonlinear analyses becomes significant. Finally, the load is applied in the transverse direction with Pmax = 100 MN. The resulting normalized transverse velocity history measured at the impact site is
0.5
Time, [s]
Time, [s]
0.6
0.4
0.2
0
–0.2
–0.4
0
0.1
0.2
0.3
0.4
0.5
Time, [s]
0.6
0.7
0.8
0.9
1 x 10
–3
Fig. 11. Transverse velocity due to transverse impact load (FGM beam).
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
shown in Fig. 11, which shows very little deviation in the amplitude and the group speed of bending. In all the cases, SESDE and RISDE based analyses agree quite well with each other. This enables us to conclude that for isotropic and anisotropic material with zero coupling, RISDE will be good enough to predict the response accurately, which will not be the case for asymmetrically stacked composite laminate. Further, the general effect of nonlinearity in terms of increasing stiffness, which in turn, introduces decreased amplitude and increased group speed can be seen in FGM, albeit for a large loading. Moreover, the effect of nonlinearity is not so pronounced in FGM beam in comparison to composite materials.
4. Conclusion To develop the necessary tool for large displacement and large rotation analysis of beam structures, a total Lagrangian formulation of a geometric nonlinear onedimensional beam element has been presented. The incremental displacement fields within the two-noded beam element are defined using the exact shape functions which are constructed from the exact solution of the static part of the governing differential equations of the beam element. These exact shape functions, that are function of length, cross-sectional and material properties of the beam element are used to generate exact stiffness and consistent mass matrices. This beam element is based on the FSDT. It has been shown that the exact shape function based element (SESDE) has a faster rate of convergence than the FSDT beam element formed using linear shape functions with reduced integration (RISDE). It is shown that as the length to depth ratio of the beam increases the SESDE has a much higher rate of convergence than RISDE. The SESDE based solutions are much closer to the solutions obtained from the plane stress analysis performed in ANSYS. Thus, the formulated SESDE is more accurate and computationally more efficient. The formulated beam element is used to analyze wave propagation behavior in composite and FGM material and the response is compared with linear wave solutions. It is observed that, due to the increased stiffness of the beam in nonlinear case, the magnitude of the nonlinear wave is less than that of linear wave and group velocity of the wave is also higher than the linear wave. Further, the nonlinear response shows presence of higher frequencies, which are not present in the linear response. Moreover, the effect of nonlinearity is more pronounced in the presence of nonzero loading than that is observed for pulse loading, which is present only for a very small duration. The nonlinear effect is predominant in composite beam in comparison to the FGM response. It is also ob-
103
served that to get the same amount of nonlinear effect, the load required in the case of composite beam is less than that required for FGM beam. Further, for composite materials with nonzero axial-flexural-shear coupling, RISDE based analysis will not be accurate enough and the SESDE becomes indispensable.
References [1] Rosen A, Loewy RG, Mathew MB. Nonlinear analysis of pretwisted rods using principal curvature transformation: Part I: Theoretical derivation, Part II: Numerical results. AIAA J 1987;25(3):470–8. [2] Heyliger PR, Reddy JN. A higher order beam finite element for bending and vibration problems. J Sound Vib 1998;126(2): 309–26. [3] Bauchan OA, Hong CH. Nonlinear composite beam theory. J Appl Mech 1998;55:156–63. [4] Atilgan AR, Hodges CH. Unified nonlinear analysis for nonhomogeneneous anisotropic beams with closed cross-sections. AIAA J 1990;29(11):1990–9. [5] Minguet P, Dugundji J. Experiments and analysis for composite blades under large deflections. Part I: Static behavior, Part II: Dynamic behavior. AIAA J 1990;28(9):1573–88. [6] Laulusa A, Bauchau OA, Theron NU. Theoretical and experimentalal investigation of the nonlinear behavior of composite beams. La Recherche Aerospatiale 1995;4:223–40. [7] Petrov E, Geradin M. Finite element theory for curved and twisted beams based on exact solutions for three-dimensional solids: Part I: Beam concept and geometrically exact nonlinear formulation, Part II: Anisotropic and advanced beam models. Comput Meth Appl Eng 1998;165:43–127. [8] Suresh S, Mortensen A. Fundamentals of functionally graded materials. London: IOM Communications Ltd; 1998. [9] Reddy JN, Chin CD. Thermomechanical analysis of functionally graded cylinders and plates. J Thermal Stresses 1998;26(1): 593–626. [10] Praveen GN, Reddy JN. Nonlinear transient thermoelastic analysis of functionally graded ceramic-metal plates. Int J Solid Struct 1998;35(33):4457–76. [11] Reddy JN. Analysis of functionally graded plates. Int J Numer Meth Eng 2000;47:663–84. [12] Reddy JN, Hsu YS. Effects of shear deformation and anisotropy on the thermal bending of layered composite plates. J Thermal Stresses 1980;3:475–93. [13] Bathe KJ. Finite element procedures. New York: Prentice Hall; 1997. [14] Gopalakrishnan S. A deep rod finite element for structural dynamics and wave propagation problems. Int J Numer Meth Eng 2000;48:731–44. [15] Khedir AA, Reddy JN. An exact solution for the bending of thin and thick cross-ply beams. Compos Struct 1997;37:195–203. [16] Reddy JN. On locking-free shear deformable beam finite elements. Comput Meth Appl Mech Eng 1997;149:113–32. [17] Chakraborty A, Mahapatra DR, Gopalakrishnan S. Finite element analysis of free vibration and wave propagation in asymmetric composite beams with structural discontinuties. Compos Struct 2002;55(1):23–36. [18] Chakraborty A, Gopalakrishnan S, Reddy JN. A new beam finite element for the analysis of functionally graded materials. Int J Mech Sci 2003;45:519–39. [19] Chakraborty A, Gopalakrishnan S. Poissons contraction effects in a deep laminated composite beam. Mech Adv Mater Struct 2003;10(3):205–25.
104
S. Agarwal et al. / Composite Structures 72 (2006) 91–104
[20] Menzala GP, Zuazua E. The beam equation as a limit of a 1-d nonlinear von-Karman model. Appl Math Lett 1999;12: 47–52. [21] Chandrashekhara K, Bangera KM. Linear and geometrically nonlinear analysis of composite beams under transverse loading. Compos Sci Technol 1993;47:339–47. [22] Sciuva MD, Icardi U. Large deflection of adaptive multilayered timoshenko beams. Compos Struct 1995;31:49–60. [23] Zienkiewicz OC, Taylor RL. Finite element in structural mechanics. McGraw Hill; 2000.
[24] Mukherjee A, Chaudhuri AS. Piezolaminated beams with large deformations. Int J Solids Struct 2002;39:4567–82. [25] Reddy JN. Mechanics of laminated composite plates. USA: CRC Press; 1996. [26] Markworth AJ, Ramesh KS, Parks WPJ. Modeling studies applied to functionally grsded materials. J Mater Sci 1995;30: 2183–93. [27] Wakashima K, Hirano T, Niino M. Space applications of advanced structural materials. USA: CRC Press; 1990. SP303:97.