- Email: [email protected]
Nonlinear elastic stress analysis in curved composite beams
Pergamoln
PII: SOO45-7949(96)00301-X
NONLINEAR
Compurers & Smctures Vol. 62, No. 5, pp. 837459, 1997 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 00457949/97 s17.00 + 0.00
ELASTIC STRESS ANALYSIS IN CURVED COMPOSITE BEAMS F. Fraternalit and G. Bilottit
TDipartimento di Ingegneria Civile, Universita’ di Salerno, Fisciano (Sa) 84084, Italy $Dipartimento di Costruzioni per I’Architettura, I.U.A.V., Tolentini 191, 30100 Venezia, Italy (Received 10 July 1995) Abstract-4 one-dimensional theory and a finite element model for the stress analysis of laminated curved beams are presented. The theory accounts for moderately large rotations, moderately large shear strains and a different elastic behaviour of the material in tension and in compression. Use is made of a penalty technique to enforce the perfect bonding condition between the adjacent laminae of the beam. Such a technique offers the possibility to obtain an approximation of the interlaminar stress field as a direct result. Some numerical results are given to show the potential of the proposed finite element model. Copyright 0 1996 Elsevier Science Ltd.
NOTATION Bo
c L X3
(or
S)
R x B {AI,
Xl,
A2r
A,}
x2
b h N Bp C’k’ L’k’ X$) (or .Skl) R’k’ x(k) glk, {Alk’,Aik’,Af’} xy’, .v’2 h’k’
1. INTRODUCI’ION
natural configuration of a curved laminated beam axis curve of B. length of C line coordinate along C curvature radius of C cross-section of B. geometric centroid of E Frknet frame associated with C (A, unit binormal vector, A2 unit normal vector, A3 unit tangent vector) coordinates on E with respect to (9, Al, .42) width of Z thickness of Z number of sublaminates composing Bo g,eneric sublaminate of B. axis curve of Bhk’ length of C”) line coordinate along C(k) curvature radius of Cck) at the generic point cross-section of fit) geometric centroid of W F&net frame associated with cYk’ cloordinates on Zk) with respect to {#a(k),A{k’,#‘} thickness of Elk)
Interlaminar stresses play a crucial role in the mechanics of laminated composite structures. Indeed, they may produce a failure of the structure due to delamination. Besides, their computation is essential to evaluate the stress field of such structures. On the other hand, it is well-known that displacement based one-dimensional or two-dimensional theories for laminated structures do not allow for an accurate computation of these stresses, since they give rise to discontinuities in transverse stresses at an interface of dissimilar material layers. The adoption of a layer-wise approximation of the displacement field through the thickness and of a method for post-computing the interlaminar stresses via the three-dimensional (30) equilibrium equations may overcome this difficulty [l, 21. However, even such a theory does not furnish interlaminar stresses as a direct result. Moreover, the use .of the above-mentioned method for computing interlaminar stresses is convenient in the case of a flat structure (e.g. a laminated plate), but is rather complicated in the case of a curved structure (a curved laminated beam or a laminated shell). A different approach was proposed in Refs [3] and [4] for the linear analysis of laminated curved beams and shells, respectively. It employs a layer-wise description of the displacement field and a penalty technique to enforce the perfect bonding constraint between the adjacent laminae of the composite structure. Such a technique allows one to obtain
837
838
F. Fraternali and G. Bilotti
direct and accurate information about interlaminar stresses, avoiding the use of the 3D equilibrium equations. The above-mentioned approaches to the stress analysis of laminated composite structures are confined to the linear-elastic field. The present work deals with the formulation of a one-dimensional, nonlinear elastic theory for the analysis of curved laminated beams made up of composite material. Such a theory generalizes that proposed in Ref. [S] and accounts for moderately large rotations of the cross-sections, moderately large shear strains and a different elastic behaviour of the material in tension and in compression (bimodular behaviour). Such a behaviour has been experimentally observed for several real composites, as weak matrix composites (see [6, 71). The beam is modelled as a collection of one-dimensional elements (sublaminates), corresponding to packages of one or more individual laminae (or layers) of the beam. Each of these elements is schematized through the theory proposed in Ref. [5]. They are supposed to be perfectly bonded at the interfaces and such a constraint is enforced into the theory via a penalty technique. The resulting mechanical model allows for an accurate prediction of the interlaminar stress field. In addition, it give rise to the possibility to compute the stress distribution within the layers of the beam. The work also contains a finite element formulation of the proposed mechanical model, which generalizes that given in Ref. [8]. It is completed by numerical applications, which show the capacities of the model to evaluate the stress fields of straight and curved laminated beams in the nonlinear elastic range. Particular attention is paid to the phenomenon of free edge stress concentration [9][14] and on the effects on the distribution of the stress
field, which derive from considering behaviour of the material.
a bimodular
2. THE MECHANICAL MODEL
2.1. Kinematical assumptions Let us schematize the natural configuration & of a curved laminated beam as a group of N one-dimensional elements BJ), @,*‘,. . . , I$“‘, each one of which corresponds to a package of one or more individual laminae (sublaminates). The following trivial relations exist between the sublaminate coordinate systems XF’, A”‘),J?$) and the beam coordinate system XI, X1, X, (Fig. 1): X,“) =
x,,
$2k’ = x2 - &’ d$”
3=-
dX,
R’k’
R
= 1 -$.
(1)
rik) being the coordinate of Yk with respect to {Y, A,, A2). We describe a generic deformation of each sublaminate by means of the one-dimensional theory proposed in Ref. [5]. This leads us to assume that the displacement field of @)(k = I, 2, . . , N) is of the kind: uck)= y(k) + (R(k) _ I)(XI”A\k’
+ _J$“A\k’) +
,,,(k)(2)
where v(“)= v’“)(_@) is the displacement field of Ckl, Rtk)= R’k)(XJJ))is the tensor field of the cross-section rotations, I is the identity tensor and, w(‘)is an axial displacement field, which models a cross-section warping.
Fig. 1. Illustration of a laminated curved beam.
Nonlinear elastic stress analysis We give to w(“) the following pression [5]:
polynomial
wtk)=: $@‘j”‘)‘l(~;‘)‘z#’
ex-
839
generic sublaminate vector
into the following
numerical
(3)
where the indices rl and r2 range over {2,3, . . . , &)} and {2,3,. . . , W}, respectively, while the quantities w$, are functions of the coordinate #). It is convenient to write R(k)as [ 151:
{WB, w(,?,w&Y,. . . , !&lo, . . , w&x,>}’ (7) whose dimension is 6 + m!? [m!?: number of w$~ coefficients in eqn (3)]. 2.2. Principle of virtual work
R’M= I + @o:’+ +, @(lr)’+ f @j(r)’+ . . . .
(4)
@’ being a skew tensor. Hereafter d(k) will indicate the axial vector of Ga). Giving t(“) as the symmetric part of the gradient utk), we now introduce the following smallness assumptions: (i) the component & is of an order 6 much smaller than 1 in correspondence to the sublaminate axis (small axial strains); (ii) the components &’ and & are of the c112order in correspondence to the sublaminate axis (moderate shear strains); (iii) the cross-secttion rotations 41”) are of the .?12 order (moderate rotations); (iv) the gradient of the warping displacement field is of the E order (infinitesimal warping). In particular, the assumption (iii) requires that terms of an order higher than 2 can be overlooked in eqn (4). Consequently, the components of the displacement field (2) can be written as:
+
[t)$
+
+
,,,(k) ,yJkY’XkP 2 119
In the derivation of the equilibrium equations and boundary conditions of the laminated beam, we make use of the Lagrangian formulation of the principle of virtual work:
s
S.GEdV-I
80
bo.GudV s 80
-1
c
po.GudA =0
JJBo
where: S is the second Piola-Kirchhoff stress tensor, E is the Green strain tensor, 6 is the variational operator, 1 is the load multiplier, b, is a given body force density per unit volume of the undeformed body and po is a given surface force density per unit area of the boundary of the undeformed body. The kinematical assumption introduced in the previous section allow us to reduce eqn (8) to the following one-dimensional expression:
&iJp’c#Jjk’]x$k’
k-l
.
(3
Finally, we assume that the different laminates of the beam are perfectly adherent to other. This hypothesis obviously implies the lowing mutual constraint equations between displacement fields II”‘:
subeach folthe
Vk~{l,2,3 We collect the generalized
-1 ,f, { Q$'Gvj"(L) + (C’:” s
,,(a = ,,‘k+ I’ on a@’ A J&k0 + 1)
+
,...,
N-l}.
displacements
(6) of the
(8)
ay’)W1”(L) + L~yw$,(L)} (i = 1,2, 3;
=0
rl = 2, 3, . . . , ack);
r2 = 2, 3, . . . , bck)).
F. Fraternali and G. Bilotti
S,, dW’
Mf’
=
_
s,,X,~)dEk’
are the corresponding (lo)
and w$, ti$, are the generalized strains of the different sublaminates (dots denote derivatives with respect to Y/j);
Ivl”’ =
s
S,, dZk’ X(X,
generalized
stresses [CCL’=
(R’“’ - JI”“)/R’~‘];
G’“‘bn,dZ’“’ + G”‘p,, ds”’ “‘) = s zll, s;xM’
(k.11 = Cl
G’“‘bo,Xk’
s
L’ll
dC’“’+
G”‘p&’
ds’“’
Nonlinear elastic stress analysis c,(1.2) _--
(;‘k’b
+
841
j@) dZ”’ 02
G"'po,Xjj’ ds”’
+;
G"'b,,,xi;' dZ”’
+
G'k'p,,X;' dsck’ bjk) >
+
Ck.1) =
G'k'beXjkJ dZk’
C3
s s
G(k)(bo2xk) -
bo,X;‘) dC”’
XC”,
G’“‘@,&,” - po,@‘) ds”’
+
8.e)
G'k'bo,~k'dW
+;
+
s
(k2) =
c3
1 G”‘bo,Xj” dEk’ 2 (s pxI
_
G'"'po,~:" ds’” 41”’
BE”
>
+;
Gdl c2
=
1 2
_
+
G’VI~,Yzk’dZk’
s
G’k’p,,Y;’ dsfk’ ppx’
(s
G”‘bo,Xjj dZ’“’ pW
F. Fraternali and G. Bilotti
842
+ ; [
j-kpo,-Xlk’ d~‘“‘]O:“(Xk4
are the generalized forces per unit length of the sublaminate axes (the superscript 2 marks secondorder forces; s’lr)denotes the line coordinate along ask));
po3(Xjk))‘fXJJ))dEk’
L’k’ = Z,l’Z
(13)
are the generalized forces on the bases of the sublaminates [LX = 0, L; 3:) = 0, L’k); x”) = Z’)(O), W( L”‘)] . It is convenient to collect the above quantities in the following numerical vectors:
(0,0, O},
+ ;
[j-k, po,Xl”dZ’n]sp’(W) (ry’, ry, I-p}, (0, 0, {O,
9”)
+ [
=
0,.
{ {Aqk’,
. . . , O},
. . , O}}’
Aqk’,
(13
Ny’},
{ibqk’,
MJk’,
M$k’},
h,po,X”’ dZ’“)]OI”‘(%“‘)
+ ; [
~,kpo,W) dZ’k;lOP’(%k))
{@, $‘, @, . . . , @, . . . , c#}9 ip2) = {{O, 0, a>, {p,
cp, &J’},
IO, 0,. . * , O>Y 6$“~”= {{Q&t’,Qg’, Qb’}, {C$“, q.‘), C$:“}, {L& L$, L$, . . . , L& . . . , Lg}}‘,
(17)
Nonlinear elastic stress analysis
Q$“’= {{O,0, O}, {q’, (O,O,.
.1
cy’, q’},
OjlT
(18)
Qf” = { {Ql"l, Qz:', Ql:'}, {cf;“, cf;“, cf;“},
843
Starting from the above constitutive assumptions and performing integrations over the cross-sections of the sublaminates, we obtain the constitutive equations between generalized stresses and strains: 9”’ = n(k’(@,” + ;I%“)
k = 1, 2, . .
, N.
(23)
Here nca is the elasticity matrix of the kth sublaminate. Its expression is given in Ref. [5].
Qy’ = {{a, o, o), (cy’, cy, cy.z’),
2.4. Total potential energy
(0,0, . . . , O>Y
(19)
2.3. Constitutive equations We assume that the composite material of each individual lamina of the beam is linearly-elastic and orthotropic. Concerning the material-symmetry axes L, T, N (L in the fiber direction, T transverse to it in the plane of the lamina, N perpendicular to such a plane), we suppose that the L, T system is oriented at an arbitrary angle /I with respect to the X,, XI system. Consequently, we write as follows the local constitutive equations which are relevant to the problem under consideration:
By supposing that the body and surface forces acting on the beam form a dead load, the total potential energy functional of the one-dimensional model under examination is the following:
+ f$U’) dS”’
- 1
k$,(ti(~*)(o))~(QK;." +@b”“>
S,, := 2QaEx S,s := QISEX+ 2QsEu
- A ,t, (s(~'(L'"'))'(Q~' + ;Qf.“,. (20)
Q,j being the material
stiffnesses in the XI, X2, X3 system. The equations rcelating the stiffness coefficients to the engineering elastic constants (Young’s moduli, Poisson’s ratios and shear moduli) in the materialsymmetry system, and their transformation in the structural system .I’,, Xz, X,, will not be given here, but they can be fc’und in Ref. [5]. With the aim of’ including the case of composite bimodular materials in the theory under examination, we follow Pate1 et al. [6] and Bert [7] assuming that the material possesses one set of elastic moduli when the fibers are elongated (EL = EA, . AL >, 0, AL being the unit vector in the fiber direction) and a different set of moduli when the fibers are shortened (EL < 0). We hence write: Qti = U - H(Et)lQ,
+ HWQ;,
(21)
H(EL) being the Heavyside function:
H(EL;I =z
0, if EL < 0 1, ifEL>O‘
(22)
(24)
It is defined over the set of kinematically admissible displacement fields of the laminated beam. The assumption of perfect bond between adjacent sublaminates requires that such displacement fields satisfy the constraint eqn (6), which can be rewritten as it follows: #+ ,
=
uf
+ I’-
on
a@’
,e-, a&k
VkE{1,2,3 where
+ I’
,...,
N-l},
(25)
F. Fraternali and G. Bilotti
844
The total potential
energy of the penalty model
is: N-l
4-I+; c
+I’-
_
Ul”
k=l
@t
I’-
=
Vi”+ ‘1 _ [@ + ‘1 _
+
):I!;,; “X;l
f#,‘l” + “&k
+ “]_y,
(27)
It is evident that the explicitation of the perfect bond conditions (25) in the functional fI is rather complicate from an analytical point of view. It is hence helpful to enforce such conditions via a penalty technique [16]-[18]. This leads us to introduce an auxiliary model (penalty model), which differs from the original one because the sublaminates are connected by continuous distributions of elastic springs instead of being perfectly bonded (Fig. 2).
where A”) is the interface between @) and B/2+‘) [kk’ = aBik”’ n ~3%~ + “1, l/q I”)is the stiffness (per unit area) of the springs acting at Nk) in the direction perpendicular to the beam plane, l/d) is the stiffness of the springs in the radial direction and I/# is the stiffness of the springs in the axial direction (Fig. 2). It is clear the sublaminate displacement fields a(‘), fi’2’ . , W, are completely independent in the ftuktional I&. When the quantities t~{~‘, nlk) and qsk’ (penalty parameters) approach 0, i.e. when the spring stiffnesses increase indefinitely, it is easy to conjecture that the penalty model approximates the original one [16H18]. Furthermore, as ~1~)-+ 0, ~4~)+ 0, 91”)--t 0, the spring reactions:
(i = 1,2,3; i not summed)
(29)
approach the values of the sublaminate interactions in the original model (interlaminar stresses). The equilibrium equations of the penalty model can be derived by means of the principle of stationary
Fig. 2. Penalty model.
Nonlinear elastic stress analysis potential
energy [19]. Such a principle
3. FINITE ELEMENT APPROXIMATION
ensures that
the couple (ii, A), where
Let us introduce a finite element discretization for the axis curve C of the laminated beam with n, elements C,, C2, . . . , Cme. We suppose that such elements belong to the Lagrange’s family [ 16, 181and are of the isoparametric type. We hence adopt the same shape functions to approximate both the geometry and the displacement field within the generic elements C, (see Fig. 3)
(yl’
1.1 @2’
it=
(30)
&v
characterizes an equilibrium admissible and it results:
state if and only if, fi is
a
aIIp=gI’(Q+crM,1)
&=o
(31)
In eqn (33) 5 is the normalized coordinate on the undistorted (or master) element (- 1 < r < 1); n is the number of nodes connected to the element; Z is the position vector of the generic point of C,, Z, and 6, are the nodal values of the fields Z and B, respectively; N,(c) (I = 1,2, . . . , N) are the Lagrange’s interpolating polynomials of the order n- 1. The Jacobians .P = aS”)jag can be computed from the equations
for all admissible variations 611. 2.5. Stability criterion We adopt the energetic stability criterion of Trefftz [19] to compute the stability points of the laminated beam. Referring to the penalty model, the Trefftz’s criterion states that the couple (il,, &) characterizes a critical equilibrium state (that is a stability point) if, and only if, there exists some V # 0 (buckling mode) such that it results:
6(62rlP)
=
&
l-I,@ I
+
a,?
+
2
845
=
a2M,n,)
pk’=
(, -$)($)
0
1,=z,=o
(32) for all admissible variations 6B.
1
I
n
g-1 Fig.
3. Finite element mesh and detail of the generic element.
846
F. Fraternali and G. Bilotti
where Z,, and Z,, are the components of Z, with respect to a fixed Cartesian coordinate frame (see Fig. 3), while N,., and NJ.{
denotes the contribution to the stiffness due to the springs acting at the sublaminates interfaces. Following a standard procedure, we build up the matrices Ku and Kp by assembling the contributions due to the individual elements of the finite element mesh: KU<,K, (e = 1, 2, . . , n,). In particular we compute the matrices K, by peforming an ulterior
(35)
3.1. Path-following procedure Utilizing eqn (33), it is possible to reduce the equilibrium equations of the penalty model [eqn (31)] into the following discretized form: SU’{K(U)U - n[Q’” + Q”‘(U)]} = 0 where U is the (secant) nodal force respectively. follows:
(36)
the nodal displacement vector; K is stiffness matrix; Q(l) and Q”) are the vectors of the first-and second-order, It is useful to represent U as it
assembly of the sublaminate contributions K’:j, K$;, . . . , KIN. Such contributions will not be given here, since they can be readily obtained by applying to each sublaminate the finite element model presented in Ref. [6]. In a similar way we build up the global force vectors Q(l) and Q”). Concerning instead the element penalty matrices K,, we refer to the Appendix. Coming back now to eqn (36), we remark that it must hold for each admissible variation HUE@’ (M is the total number of equations). Hence, it is equivalent to the following system of nonlinear algebraic equations: R(U, 1) = K(U)U - L[Q”’ + Q”‘(U)] = 0.
(40)
It is helpful to solve eqn (40) in the context of a displacement control procedure. This can be done by adding to eqn (40) an extra equation, which forces a particular component of U (say the pth one) to assume a prescribed value p (see [21], [6]). We hence introduce the following extended system of equations:
(41)
n, being the total number element mesh
of nodes of the finite
nti=(n,+l)xn+l.
(38)
The stiffness matrix K can be expressed as K=Ku+Kp
(39)
where Ku denotes the stiffness matrix for the unconstrained problem (i.e. for the disconnected sublaminates), while Kp (penalty stiffness matrix)
where e, is the unit vector which has the pth entry equal to 1 and all the other entries equal to 0. An efficient solution scheme for the matrix eqn (41) can be obtained by employing the well known Newton-Raphson method. It consists of a partitioning method also called bordering algorithm in the literature (see Refs [22,23]). For the description of such an algorithm we refer to the work [6]. 3.2. Procedure for the computation of stability points The variational eqn (3l), which rules the problem of computing the stability points of the penalty model
Nonlinear elastic stress analysis and the corresponding buckling modes, admits the following discretized formulation: &F&&f,
n)v =
v
(42)
where V is the nodlal displacement vector corresponding with the buckling mode, while KT is the tangent stiffness matrix: K1V=~R(U+aV,l)i,=c’fV~4.
(431
We decompose KT as follows KT = I& + K,
(44)
where K,, is the tangent stiffness matrix unconstrained problem i.e. KT,V = $
where V0 denotes a prescribed value of the pth component of V. It is convenient to let the position p coincide with the equation number where the smallest diagona! term of the (factored) tangent stiffness matrix appears. Concerning the scaling factor Vo, we set
(49) where V. is the initial approximation to V. Since the matrix KT becomes progressively ill-conditioned as (U, ,I) approaches a stability point, it is numerically convenient to transform system (48) into the following form [24]
of the
{Ku(U + aV)W + aV]
- .IQi2’(U + uV>]
i1=0
vv EP
(45)
while KP” is the tangent stiffness matrix due to the springs acting at the sublaminate interfaces:
KspV = $ (Kp(U + aV)[U + aV]> Ia-r! tlv E SF.
(46)
We build up KS,>by assembling the contributions Such contributions can be obtained by applying to the generic sublaminate the finite element model presented in Ref. [6]. Concerning instead KTO,we give in the Appendix the expression of the contribution due to the generic element of the fide element mesh. Awaiting the arbitrariness of SU, eqn (42) is equivalent to the following system of algebraic equations due to the differenl sublaminates.
;Kf(U, a)v = 0.
(47)
The direct computation of stability points can be performed by adding to (47) the equilibrium eqns (41) and an ulterior equation preventing that V becomes the null vector [24, 61. We hence introduce the following extended system
R**(u,
841
v, A, p)
=
=o
(48)
y being a real number greater than zero. A solution scheme for eqn (50) can be obtained by means of a partitioning method similar to that employed for eqn (41) [24,6]. The numerical procedures for eqns (41), (50) can be combined to formulate an overall algorithm for stability analysis [6]. Within such an algorithm, the path-following method is used to follow a solution branch (e.g. the fundamental equilibrium path). Near a stability point the extended system (50) can be instead applied to compute this point. To estimate how to switch between the two procedures, it is helpful to check for the sign of det Kr along the solution branch. Indeed a change in the sign of det KT indicates that a singular point of KT (that is a stability point) has been passed. Once such a point (U,, A,.)has been computed, it is possible to follow the post-critical branch by employing again the path-following method. In particular, if the stability point is a bifurcation point (VTQ = 0), the switching of the path-following procedure on the bifurcated path can be performed by adding to U,. a vector which is proportional to the buckling mode V:
(51) Here { denotes a scaling factor to be determined by imposing that the couple (U, A<) statisfies the equilibrium eqn (41). 4. NUMERICAL RESULTS
We present in this section the results of a finite element study on stress distribution in straight and
F. Fraternali and G. Bilotti
848
I
L/z :
rigid dars
//
h=O.OSL
_ w--v
unimodulor bimodulor
v/h Fig. 4. Loadxenter
deflection curves of a straight laminated beam (D& = Jx Q&g dZ: bending stiffness of the unimodular beam).
curved laminated beams. The adoption of the finite element mode1 formerly described, allowed us to extend such a study beyond the linear-elastic field, taking into account moderately large rotations and shear strains (geometrical nonlinearity) and a different elastic behaviour of the material in tension and in compression (mechanical nonlinearity). Two examples were studied: a three-layer straight beam (example 1) and a three-layer shallow arch (example 2). Both of them consider a cross-ply lamination scheme ([O/90/0]), a central point force and hinged supports. To obtain accurate stress distributions, we let the number of sublaminates coincide with the number of individual layers (N = 3). In each sublaminate we modelled the cross-section warping through the following function:
iv(k)= w&:‘(#‘)2 + w$(%;))’
(k = 1, 2, 3).
us to avoid the numerical incovenient, which is known as “locking” of the finite solution in the literature [ 16, 181. Concerning the calculation of the stress field, we employed different procedures for tangential and normal stresses. In detail, on a cross-section corresponding to a Gauss point, we computed the through-the-thickness distribution of shear stress S,, by interpolating the interlaminar values (A$‘),Al’)) through a layer-wise quadratic shape function. Such
(52)
In both of the examples, we used four-node isoparametric finite elements and a Gauss formula with 4 points to compute, on each element, either the “unconstrained” stiffness matrix Ku. and the penalty stiffness matrix Kp,. In this last one we set
$k’
=
#’
=
g’
=
IO-3
1
maxWu,, 1
(k = 193,
(53)
where max{&,,) is the largest of the (translational) diagonal elements of Ku. The above positions allowed
Fig. 5. Deformation of the terminal part of the beam analysed in the example 1 (unimodular material behaviour, FL'/48D&h = 0.203, displacement magnification ratio = 5).
Nonlinear elastic stress analysis
a function was calibrated librium conditions:
by imposing
s
&3 dEk’ = Nj) X(k)
849
the equi-
(k = 1,2,3).
3j@’ +bh’k’-
(54)
We hence obtained the following distribution:
,$k)+ $ - I) 4
for
2015-
lo-
5-
upper interface (O/90) unimodular behoviour free edges
I
-1
I,
x-0.001
O-
-5.
PI I31
-10.
h
h-0.05L
-15.
I
I
I
I
I
I
I
0.002
0.000
0.004 x3
/
I
I 0.008
0.006
OA3’IO
L
lower interface
(90/O) unimodulor behaviour free edges
X=0.466
IA HI
lJ2 :
k
; 0.000
-20
I
I
0.002
I
1
I
0.004 x3
I
0.006 /
L/2
1
I
h
h-0.05L
I
0.008
I 0.t3’IO
L
Fig. 6. Interlaminar distributions of the transverse shear stress for the beam analysed in the example I.
850
F. Fratemali and G. Bilotti
We instead computed the distribution stress Sa from the constitutive equations,
of normal obtaining:
The elastic properties
of the material
were taken
as
& = Q$)[(_$$" + $_@*2') + (@J) + f@$".")~/
8, =rnl?t,b;=
IA \
x s
-lO-15-
t
:
I-/*
L/2
h=O.OSL
I
upper interface (O/90) unimodulor behoviour free edges
-2o-
1
-25
I
I
0.000
0.002
I
I 0.004 x3
I
I
0.006 /
I
1
0.008
O.( 10
L
(b)
30
I
251
lower interface (90/O) unimodulor behoviour free edges
I4 IEII
h
h=O.OSL
5-j 0
-5 -10
-15
\
x=0.001
Ii_
_ .
_
_ _ _
_ _ _
_
x-o.203 X-O.468
!
0.000
Fig. 7. Interlaminar
I
I
0.002
distributions
I
I
0.004
I
I
0.006
I
I
0.008
I
O.( IO
of the transverse normal stress for the beam analysed in the example 1.
Nonlinear elastic stress analysis
Here bL, br, gN are Young’s moduii in the materialsymmetry system (L, T, N); YLT, YLN,YTN are shear moduli; vLr, vLN,i’rN are Poisson’s ratios; m is the bimodularity ratio (m = 6, I&); the apex + (-) marks values in tension (compression). The bimodular - behaviour that we considered corresponds to m = 0.5 (material weaker in compression than in tension).
0.8
Example 1. A straight, [O/90/0] beam (/z/L = O-05), simply-supported at the ends and subjected to a central point force, was analysed. Owing to symmetry, only one half of the beam was modelled, with a 25 element mesh. In detail, eight elements were used to discretize the part of the beam which starts from one base and spreads out for 0.025 for the span L; five, four and eight elements were
(=I
0.6
‘0
,*’
,’
0.,2
-0.50
8.0
I
-0.33
lb)
6.0 -
\
4.0
I
-
A-0.001 0.205 0.466
_____ ___
I
/ I
I
I
I
I
I
I fme Bdgt3B l unimodular
I
I
brhoviour
-0.17
0.17
-0.00 /
l I
I
x2
IL
I
’ A@ 0.
,I=’
0.4
0.0
851
---x-i
IO
t-l
I
I
I I
I I
X,=O.O42L
I‘,
-
1:
I free edges
I unimodulor behoviour i0
Fig. 8. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analysed in the example I, when the edges are free to warp and the material is unimodular.
F. Fraternali and G. Bilotti
852
instead used to discretize the successive parts of sizes 0.050 L, 0.10 L and 0.325 L, respectively. Two different kinematical behaviours of the bases were considered: bases free to warp (free edges) and bases constrained to rotate rigidly (rigid ends). In addition, either a unimodular and a bimodular behaviour of the material was taken into account. The
load-deflection curves corresponding to the examined cases are shown in Fig. 4. Figure 5 shows the deformation of the terminal part of the beam when the material is unimodular and the bases are free to warp. A marked warping of the central layer can be observed. It is due to the strong difference existing between the axial Young’s moduli
(a)
0.8
0.6 LL \ $: v,
0.4
5 0.2 1 unimodulor behoviour 1 0.0
I -0.33
-0.50
I -0.00
-0.17 x2
15
(b)
‘O-
~-0.001 _____ 0.205 0.466 __-
/
0.50
h
I
I
I
I
I
I
h=O.OSZL
*--------I *s----1___
I I
5-
I 0.33
0.17
---*_ --mm,
t
I _____---
-----
I
-5-
I -10
-0.50
I -0.33
I
I
I
-0.17
0.17
-0.00 x2
/
rig-d edge8 unimodular
behoviour
I
0.33
0. 10
h
Fig. 9. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analyzed in the example I, when the edges are riBid and the material is unimodular.
Nonlinear elastic stress analysis of the edge (0) tayers (Q, = gL) and of thecentral (90) layer (8, = c?~= E’J25). The above-menlioned through-the-thickness nonhomogeneity of material properties is also responsible for a marked stress concentration near the bases of the beam, when the same are not stiffened (free edges). Here the interlaminar stresses appear
853
to be tmbonded, as it is shown in Figs 6 and 7. These figures contain plots of the axial distributions of the transverse shear stress & and of the transverse normal stress & at the interfaces between the edge (90) layers and the central (0) layer, for several values of the load parameter 1 = FL3/(48D&h).
(a) 0.66
-
-
X-O.001
IL
\
0.4.6
x 2
0.26
0.06
-0.14
-0.50
I
I -0.33
I -0.00
-0.17 x2
8.0
6.0
IL
-
-
_____ ‘=;.;i: __c 01460
I
2.0 -
tr \ I ‘. \
,I
2 n
I
fl_V’ (_=-=-~-%_W~_ --_---
-2.0
-
-4.0 -0.50
I -0.33
X,-0.042L
‘\ I\
, I I I / ,’1 I ,’ I I #) I I’ I 1 ,’
4.0 -
0. i0
h
I
\ 2
/
I 0.33
I 0.17
‘. \ ‘. t ‘. \ l.\
----
I
I
I
I
I
I
free edges bimodulor
khoviour
&-=O.**
I
I 0.17
I -0.00
-0.17 x2
/
I 0.33
0.50
h
Fig. 10. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analysed in the example I, when the edges are free to warp and the material is bimodular.
CAS 62/5-C
F. Fraternali and G. Bilotti
854
It can be observed that the given shear stress distributions are generally unsymmetric. Such result, which is due to the above-mentioned edge stress perturbation, is more evident as the load increases, the edges are free to warp and the material is bimodular. In the same circumstances, the relative ti$X$)‘) “weight” of warping terms (ti$‘@‘, increases in eqn (56) and the axial stress distributions appear to be strongly nonlinear. Moreover, in the bimodular case, such distributions exhibit a marked slope discontinuity in correspondence to the zero (Fig. 10).
Figure 6 shows that the transverse shear stress exhibits edge singularities in correspondence to the interlaminar surfaces (O/90) and (90/O). Such singularities are opposite in sign and their powers increase with the load. Even the transverse normal stress edge singularities are in correspondence to the interlaminar surfaces (see Fig. 7). However, the power of the singularity at the upper surface decreases with the load. The power of the other singularity instead increases. The marked oscillation about 0 of the numerical results presented in Figs 6 and 7 is due to the peculiar nature of the penalty model, in which the layers of the beam are connected by means of very stiff springs. Free edge stress concentration is common in composite laminates and well recognized in the technical literature [9-141. It can produce a failure of the composite structure due to delamination. For this reason, methods of reducing the edge stress concentration have been proposed. They essentially consist of reinforcing the stiffness at the free edges, for instance by means of “cap” (a lamina or a collection of laminae zrapped around the edge) [ 13-141. We now pass to analyse the stress distribution within the layers of the beam. Through-the-thickness distributions of the shear stress & and of the axial normal stress S,, were computed at the Gauss point X, = 0.042 L (Figs 8-10). Distinction was made between the free-edge (Fig. 8) and the rigid-edge conditions (Fig. 9). In the former case, the effect of a bimodular behaviour of the material was also considered (Fig. 10).
2.0
Example 2. The second example that we considered refers to a shallow, [O/90/0] arch, which is simply supported at the ends and subjected to a central point force. The axis-line was supposd to be parabolic and its aspect ratio f/l (J central raise, I: base length) was taken as l/20. The thickness-tolength ratio h/l was instead taken as l/200. A uniform finite element mesh composed of 25 elements was employed to discretize the axis-line. Figure 11 shows the load-deflection curves of the considered arch, for a unimodular (8; = S,t) and a bimodular (8; = 0.5 8:) material behaviour (the deflection is referred to the point P shown in the figure). A bifurcation of the fundamental equilibrium path (with a skew buckling mode), a marked influence of the material bimodularity on the critical load and a slightly unstable asymptotic behaviour of the arch along the post-critical path can be observed. Due to the relevant slenderness of the arch, we
I
1.5
1 .o
h=O.O05t f30.051
0.5 - - - - bimodulor -
unimodular
0.0
Fig. I I. Load-displacement
curves of a shallow, [O/90/0] arch, for a unimodular (8, = Bt ) and a bimodular (8, = 0.5 8:) material behaviour.
Nonlinear elastic stress analysis
practically obtained the same results considering the bases as rigid or free to warp. This led us to omit a specific study on ledge effects in the present case. We instead computed the through-the-thickness distributions of the transverse shear stress SX3and of the axial normal stress &, at the Gauss point X, = 0.3665 L, for the unimodular and the bimohular
855
interlaminar surfaces is apparent when the material behaviour is of the bimodular-type. As in the previous example, a slope discontinuity of the normal stress distribution in correspondence to the zero is either apparent in such a case (Fig. 13). 5. CONCLUSIONS
case (kegs 12, 13).
A
more
pronounced
stress
gradient
near
the
A
theory
and
finite
element
model
for
0.6 IL \
w 2
cn
0.4
0.2 unimodular
khoviour
0.0 -0.60
-0.33
-0.17
-0.00 X2
20.0
10.0
Jb)
-
OS.0__
-
-
_lC~.O _.___------_ *e--
-
-;7
Nd
P’
-0.50
I -0.33
I
I
I
I 1
I I
X,=O.J665L
P’
P’
I 4
I I
I
-
I
I
-__
I
khoviour
-0.17
/
I
A-0.001
--mm_
@OOW
0.472 0.959 1.732
t
I -0.00 X2
Fig. 12. Through-the-thickness
0.
h
I
I unimodulor -30.0
0.33
I
_~_-~-,-=. _-
9’ -20.0
-
/
0.17
I
0.17
0.33
0,
h
distributions of the transverse shear stress and of the axial normal stress for the arch analysed in the example 2, when the material is unimodular.
the
F. Fraternali and G. Bilotti
856
non-linear elastic analysis of laminated curved beams has been proposed. Even remaining in a one-dimensional context, such a theory allows for an accurate description of the interlaminar stress field and gives rise to the possibility to compute the stress field within the laminae of the structure. The considered nonlinear effects are both geometrical and mechanical. The first refer to moderately large rotations
0.6
;“’
and shear strains; the second to a different elastic behaviour of the material in tension and in compression. The proposed finite-element model allows one to follow the fundamental equilibrium path of the beam, to compute its stability points and to analyse the post-critical behaviour. It has been applied to a numerical study on the stress distribution in straight
I
I XS=C.3665L
I 0.6 \
\
0.4
bimodular
0.0 -0.50
I -0.33
I -0.17
behoviour
I -0.00
1 I
0.17
0.33
0.
Xz / h 20.0
(b)
10.0
0.0
t.z 5
I
I
I
I
I
X,=O.3665L
, / /
I
LL \
I
-10.0
I
--______
__-- ___--__--- _.--4 .s( /@ /cc I
I
I
I
I I
cc
-20.0
A=O.OOl _____ 0.472 _--
0.956
I bimodular behoviour I c-=0.5&
-30.0
I
-0.50
-0.33
I
I
I -0.17
-0.00 X2
/
0.17
I 0.33
0
h
Fig. 13. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the arch analysed in the example 2, when the material is bimodular.
Nonlinear elastic stress analysis and curved laminated beams. The presented results show the potentia:l of the adopted theory in analysing local effects, such as free edge stress concentration
and a bimodular behaviour of the material. In future works, an extension of the present theory will be formulated,
in order
to study
delamination
and
failure of composite beams.
REFERENCES 1. J. N. Reddy,
A generalization of two-dimensional theories of laminated composite plates. Commun. Appl.
Numer. Meth. 38, 173-180 (1987).
2. J. N. Reddy, E. J. Barber0 and J. L. Teply, A plate bending element based on a generalized laminate plate theory. ht. J. Numer. Meth. Engng 28, 2275-2292
Nomizu and S. Kobayashi, Foundations of Differential Geometry. Wiley, New York (1963). 21. I. L. Bathoz and G. Dhatt, Incremental displacement algorithms for nonlinear problems. Inr. J. Numer. Meth. Engng 14, 1262-1267 (1979).
20. K.
22. E. Riks, The application of Newton’s method to the problem of elastic stability. J. Appl. Me&. 39, 1060-1066 (1972). 23. H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems. In: Applications of Bijiircation - Theory _ (Edited by . _P. H. Rabinowitz), pp. 359-384. Academic Press, New York (1977). 24. J. C. Simo and P. Wriggers, A general procedure for the direct computation of turning and bifurcation points. Int. J. Numer. Meth. Engng 30, 155-176 (1990).
(1989).
3. L. Ascione and F. Fraternali, A penalty model for the analysis of laminated curved beams. Cornput. Strucr. 45(5/6) 985-999 (1992).
4. F. Fraternali and J. N. Reddy, A penalty model for the analysis of laminated composite shells. Int. J. Solids Struct. 30(24). 3337-3355
(1993).
5. L. Ascione &d F. Fraternali,’ A moderate rotation theory of laminated composite curved beams. Znr. J. Engng Analys. Des. I(2), 161-176 (1994).
6. H. P. Pate], J. L. Turner and J. D. Walter, Radial tire cord-rubber composites. Rubber Chem. Technol. 49, 1095-I 110 ( 1976). 7. C. W. Bert, Models for fibrous composites with different properties in tension and in compression. J. Emma Mater. 344-349 $977).
Technol.
Trans.
ASME
857
APPENDIX Penalty stiffness matrices
In Section 3 we introduced the penalty stiffness matrix KP to distinguish the contribution to the global stiffness matrix, which is due to the springs connecting the sublaminates of the beam. We now consider the contributions to KP due to the generic element of the finite element mesh (say the eth one). It is convenient to order the nodal displacement vector of such an element as it follows:
994(4). ~ ,_
(Al)
8. L. Ascione and F. Fraternali, A finite element analysis of the stability of bimodular composite curved beams. Znt. J. Engng Analys. Des. l(3), 315-334 (1994).
9. N. J. Pagano, On the calculation of interlaminar stresses in composite laminates. J. Comput. Mater. 8, 65-82
where
(1974).
10. N. J. Pagano, #Stress fields in composite laminates. Int. J. Solicis Struct. 14, 385400
(1978).
11. A. S. D. Wang and F. W. Crossman, Some new results on edge effects in symmetric composite laminates. J. Compur. Mater. 11, 92-106 (1977). 12. J. D. Whitcomb, I. S. Raju and J. G. Goree, Reliability of the finite-element method for calculating free-edge stress in compor,ite laminates. Comput. Struct. 15(l), 23-37 (1982). 13. J. N. Reddy, Reduction of free-edge stress concentration. J. Appl. Mech. 52, 801-805 (1985). 14. R. Y. Kim, A technique for prevention of delamination. Mech. Comput. Reo. Oct. 28-30, Dayton, Ohio (1981). 15. T. Argyris, An excursion into large rotations. Comput. (1982).
Methods
Appl.
Mech.
Engng
32,
642)
The matrix Kp
85-155
N. Reddy, An Introduction to the Finite Element Method. 2nd edn. McGraw-Hill. New York (1992). 17. J. N. Reddy, Applied Functional Analysis and Variational Methods in Enaineerina. .- Krieaer, Melbourne.
16. J.
FL (1992). 18. G. F. Carey and J. T. Oden, Finite Elements: A Second Course. Prentice-Hall, Englewood Cliffs, NJ (1983). 19. B. Budiansky, Theory of buckling and post-buckling behaviour ot’ elastic structures. In-Advances in Applie; Mechanics (Edited by Chia-Shun Yih), pp. l-65. Academic Press, New York (1974).
where K$h’ =
(wp’ + wp’+)%“)(wp+
(A4)
F. Fraternali and G. Bilotti
858
We now pass to consider the tangent, penalty stiffness matrix KT?,i.e. the matrix defined by eqn (46) of Section 3. Then the contribution to KT~due to the generic element of the finite element mesh is given by the equation:
(AW
x L”‘(W$‘+ + fW(,“‘+); d< NY)P + ‘)d[
(A6)
where
_
(A7)
In eqns (A4)-(A7) we denoted by NP’ and Nf + ‘) the diagonal matrices which have non-zero terms equal to N,(t). In the same equation we put:
KY;:+“=
b & NF+“J’A.” x (Wb”+ I,- + Wp + I,-) _ +I, & 2
(A14)
W#k’+.‘k+ I)-, =
a
=
&I,
&
+ 1,;
,j
1000
0
0100 0 0
0 -X,
=
b’“’
I
X,
b” + I,
-x2
0
0
X,0 0 0 x;2 x,x*
“’
0
‘.’
0
o...o,..o Xi .”
0
x;
“’
x;
x
II1,
1
(Wb”‘W
+wz”+“-)7p
(A9)
Nonlinear elastic stress analysis
The matrices I”‘), PC+“, which appear in eqns (A13) and (Al6), are given
859
by:
i 00 0 0
pw.cr + Ill
=
0
.’
0
”
0
”
0
‘.
0
.’
0 0 0
0
”
0 0
0
10 o0 0o / o
. .
.
0 . .
0 0 0
‘. .
0
0
fj
(Al7)
,.
.
0
,I\“‘,Ask)and ,I\“’being the reactions of the springs interposed between the kth and the (k + I)th sublaminate [see eqn (29) of Section 21.
PII: SOO45-7949(96)00301-X
NONLINEAR
Compurers & Smctures Vol. 62, No. 5, pp. 837459, 1997 Copyright 0 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 00457949/97 s17.00 + 0.00
ELASTIC STRESS ANALYSIS IN CURVED COMPOSITE BEAMS F. Fraternalit and G. Bilottit
TDipartimento di Ingegneria Civile, Universita’ di Salerno, Fisciano (Sa) 84084, Italy $Dipartimento di Costruzioni per I’Architettura, I.U.A.V., Tolentini 191, 30100 Venezia, Italy (Received 10 July 1995) Abstract-4 one-dimensional theory and a finite element model for the stress analysis of laminated curved beams are presented. The theory accounts for moderately large rotations, moderately large shear strains and a different elastic behaviour of the material in tension and in compression. Use is made of a penalty technique to enforce the perfect bonding condition between the adjacent laminae of the beam. Such a technique offers the possibility to obtain an approximation of the interlaminar stress field as a direct result. Some numerical results are given to show the potential of the proposed finite element model. Copyright 0 1996 Elsevier Science Ltd.
NOTATION Bo
c L X3
(or
S)
R x B {AI,
Xl,
A2r
A,}
x2
b h N Bp C’k’ L’k’ X$) (or .Skl) R’k’ x(k) glk, {Alk’,Aik’,Af’} xy’, .v’2 h’k’
1. INTRODUCI’ION
natural configuration of a curved laminated beam axis curve of B. length of C line coordinate along C curvature radius of C cross-section of B. geometric centroid of E Frknet frame associated with C (A, unit binormal vector, A2 unit normal vector, A3 unit tangent vector) coordinates on E with respect to (9, Al, .42) width of Z thickness of Z number of sublaminates composing Bo g,eneric sublaminate of B. axis curve of Bhk’ length of C”) line coordinate along C(k) curvature radius of Cck) at the generic point cross-section of fit) geometric centroid of W F&net frame associated with cYk’ cloordinates on Zk) with respect to {#a(k),A{k’,#‘} thickness of Elk)
Interlaminar stresses play a crucial role in the mechanics of laminated composite structures. Indeed, they may produce a failure of the structure due to delamination. Besides, their computation is essential to evaluate the stress field of such structures. On the other hand, it is well-known that displacement based one-dimensional or two-dimensional theories for laminated structures do not allow for an accurate computation of these stresses, since they give rise to discontinuities in transverse stresses at an interface of dissimilar material layers. The adoption of a layer-wise approximation of the displacement field through the thickness and of a method for post-computing the interlaminar stresses via the three-dimensional (30) equilibrium equations may overcome this difficulty [l, 21. However, even such a theory does not furnish interlaminar stresses as a direct result. Moreover, the use .of the above-mentioned method for computing interlaminar stresses is convenient in the case of a flat structure (e.g. a laminated plate), but is rather complicated in the case of a curved structure (a curved laminated beam or a laminated shell). A different approach was proposed in Refs [3] and [4] for the linear analysis of laminated curved beams and shells, respectively. It employs a layer-wise description of the displacement field and a penalty technique to enforce the perfect bonding constraint between the adjacent laminae of the composite structure. Such a technique allows one to obtain
837
838
F. Fraternali and G. Bilotti
direct and accurate information about interlaminar stresses, avoiding the use of the 3D equilibrium equations. The above-mentioned approaches to the stress analysis of laminated composite structures are confined to the linear-elastic field. The present work deals with the formulation of a one-dimensional, nonlinear elastic theory for the analysis of curved laminated beams made up of composite material. Such a theory generalizes that proposed in Ref. [S] and accounts for moderately large rotations of the cross-sections, moderately large shear strains and a different elastic behaviour of the material in tension and in compression (bimodular behaviour). Such a behaviour has been experimentally observed for several real composites, as weak matrix composites (see [6, 71). The beam is modelled as a collection of one-dimensional elements (sublaminates), corresponding to packages of one or more individual laminae (or layers) of the beam. Each of these elements is schematized through the theory proposed in Ref. [5]. They are supposed to be perfectly bonded at the interfaces and such a constraint is enforced into the theory via a penalty technique. The resulting mechanical model allows for an accurate prediction of the interlaminar stress field. In addition, it give rise to the possibility to compute the stress distribution within the layers of the beam. The work also contains a finite element formulation of the proposed mechanical model, which generalizes that given in Ref. [8]. It is completed by numerical applications, which show the capacities of the model to evaluate the stress fields of straight and curved laminated beams in the nonlinear elastic range. Particular attention is paid to the phenomenon of free edge stress concentration [9][14] and on the effects on the distribution of the stress
field, which derive from considering behaviour of the material.
a bimodular
2. THE MECHANICAL MODEL
2.1. Kinematical assumptions Let us schematize the natural configuration & of a curved laminated beam as a group of N one-dimensional elements BJ), @,*‘,. . . , I$“‘, each one of which corresponds to a package of one or more individual laminae (sublaminates). The following trivial relations exist between the sublaminate coordinate systems XF’, A”‘),J?$) and the beam coordinate system XI, X1, X, (Fig. 1): X,“) =
x,,
$2k’ = x2 - &’ d$”
3=-
dX,
R’k’
R
= 1 -$.
(1)
rik) being the coordinate of Yk with respect to {Y, A,, A2). We describe a generic deformation of each sublaminate by means of the one-dimensional theory proposed in Ref. [5]. This leads us to assume that the displacement field of @)(k = I, 2, . . , N) is of the kind: uck)= y(k) + (R(k) _ I)(XI”A\k’
+ _J$“A\k’) +
,,,(k)(2)
where v(“)= v’“)(_@) is the displacement field of Ckl, Rtk)= R’k)(XJJ))is the tensor field of the cross-section rotations, I is the identity tensor and, w(‘)is an axial displacement field, which models a cross-section warping.
Fig. 1. Illustration of a laminated curved beam.
Nonlinear elastic stress analysis We give to w(“) the following pression [5]:
polynomial
wtk)=: $@‘j”‘)‘l(~;‘)‘z#’
ex-
839
generic sublaminate vector
into the following
numerical
(3)
where the indices rl and r2 range over {2,3, . . . , &)} and {2,3,. . . , W}, respectively, while the quantities w$, are functions of the coordinate #). It is convenient to write R(k)as [ 151:
{WB, w(,?,w&Y,. . . , !&lo, . . , w&x,>}’ (7) whose dimension is 6 + m!? [m!?: number of w$~ coefficients in eqn (3)]. 2.2. Principle of virtual work
R’M= I + @o:’+ +, @(lr)’+ f @j(r)’+ . . . .
(4)
@’ being a skew tensor. Hereafter d(k) will indicate the axial vector of Ga). Giving t(“) as the symmetric part of the gradient utk), we now introduce the following smallness assumptions: (i) the component & is of an order 6 much smaller than 1 in correspondence to the sublaminate axis (small axial strains); (ii) the components &’ and & are of the c112order in correspondence to the sublaminate axis (moderate shear strains); (iii) the cross-secttion rotations 41”) are of the .?12 order (moderate rotations); (iv) the gradient of the warping displacement field is of the E order (infinitesimal warping). In particular, the assumption (iii) requires that terms of an order higher than 2 can be overlooked in eqn (4). Consequently, the components of the displacement field (2) can be written as:
+
[t)$
+
+
,,,(k) ,yJkY’XkP 2 119
In the derivation of the equilibrium equations and boundary conditions of the laminated beam, we make use of the Lagrangian formulation of the principle of virtual work:
s
S.GEdV-I
80
bo.GudV s 80
-1
c
po.GudA =0
JJBo
where: S is the second Piola-Kirchhoff stress tensor, E is the Green strain tensor, 6 is the variational operator, 1 is the load multiplier, b, is a given body force density per unit volume of the undeformed body and po is a given surface force density per unit area of the boundary of the undeformed body. The kinematical assumption introduced in the previous section allow us to reduce eqn (8) to the following one-dimensional expression:
&iJp’c#Jjk’]x$k’
k-l
.
(3
Finally, we assume that the different laminates of the beam are perfectly adherent to other. This hypothesis obviously implies the lowing mutual constraint equations between displacement fields II”‘:
subeach folthe
Vk~{l,2,3 We collect the generalized
-1 ,f, { Q$'Gvj"(L) + (C’:” s
,,(a = ,,‘k+ I’ on a@’ A J&k0 + 1)
+
,...,
N-l}.
displacements
(6) of the
(8)
ay’)W1”(L) + L~yw$,(L)} (i = 1,2, 3;
=0
rl = 2, 3, . . . , ack);
r2 = 2, 3, . . . , bck)).
F. Fraternali and G. Bilotti
S,, dW’
Mf’
=
_
s,,X,~)dEk’
are the corresponding (lo)
and w$, ti$, are the generalized strains of the different sublaminates (dots denote derivatives with respect to Y/j);
Ivl”’ =
s
S,, dZk’ X(X,
generalized
stresses [CCL’=
(R’“’ - JI”“)/R’~‘];
G’“‘bn,dZ’“’ + G”‘p,, ds”’ “‘) = s zll, s;xM’
(k.11 = Cl
G’“‘bo,Xk’
s
L’ll
dC’“’+
G”‘p&’
ds’“’
Nonlinear elastic stress analysis c,(1.2) _--
(;‘k’b
+
841
j@) dZ”’ 02
G"'po,Xjj’ ds”’
+;
G"'b,,,xi;' dZ”’
+
G'k'p,,X;' dsck’ bjk) >
+
Ck.1) =
G'k'beXjkJ dZk’
C3
s s
G(k)(bo2xk) -
bo,X;‘) dC”’
XC”,
G’“‘@,&,” - po,@‘) ds”’
+
8.e)
G'k'bo,~k'dW
+;
+
s
(k2) =
c3
1 G”‘bo,Xj” dEk’ 2 (s pxI
_
G'"'po,~:" ds’” 41”’
BE”
>
+;
Gdl c2
=
1 2
_
+
G’VI~,Yzk’dZk’
s
G’k’p,,Y;’ dsfk’ ppx’
(s
G”‘bo,Xjj dZ’“’ pW
F. Fraternali and G. Bilotti
842
+ ; [
j-kpo,-Xlk’ d~‘“‘]O:“(Xk4
are the generalized forces per unit length of the sublaminate axes (the superscript 2 marks secondorder forces; s’lr)denotes the line coordinate along ask));
po3(Xjk))‘fXJJ))dEk’
L’k’ = Z,l’Z
(13)
are the generalized forces on the bases of the sublaminates [LX = 0, L; 3:) = 0, L’k); x”) = Z’)(O), W( L”‘)] . It is convenient to collect the above quantities in the following numerical vectors:
(0,0, O},
+ ;
[j-k, po,Xl”dZ’n]sp’(W) (ry’, ry, I-p}, (0, 0, {O,
9”)
+ [
=
0,.
{ {Aqk’,
. . . , O},
. . , O}}’
Aqk’,
(13
Ny’},
{ibqk’,
MJk’,
M$k’},
h,po,X”’ dZ’“)]OI”‘(%“‘)
+ ; [
~,kpo,W) dZ’k;lOP’(%k))
{@, $‘, @, . . . , @, . . . , c#}9 ip2) = {{O, 0, a>, {p,
cp, &J’},
IO, 0,. . * , O>Y 6$“~”= {{Q&t’,Qg’, Qb’}, {C$“, q.‘), C$:“}, {L& L$, L$, . . . , L& . . . , Lg}}‘,
(17)
Nonlinear elastic stress analysis
Q$“’= {{O,0, O}, {q’, (O,O,.
.1
cy’, q’},
OjlT
(18)
Qf” = { {Ql"l, Qz:', Ql:'}, {cf;“, cf;“, cf;“},
843
Starting from the above constitutive assumptions and performing integrations over the cross-sections of the sublaminates, we obtain the constitutive equations between generalized stresses and strains: 9”’ = n(k’(@,” + ;I%“)
k = 1, 2, . .
, N.
(23)
Here nca is the elasticity matrix of the kth sublaminate. Its expression is given in Ref. [5].
Qy’ = {{a, o, o), (cy’, cy, cy.z’),
2.4. Total potential energy
(0,0, . . . , O>Y
(19)
2.3. Constitutive equations We assume that the composite material of each individual lamina of the beam is linearly-elastic and orthotropic. Concerning the material-symmetry axes L, T, N (L in the fiber direction, T transverse to it in the plane of the lamina, N perpendicular to such a plane), we suppose that the L, T system is oriented at an arbitrary angle /I with respect to the X,, XI system. Consequently, we write as follows the local constitutive equations which are relevant to the problem under consideration:
By supposing that the body and surface forces acting on the beam form a dead load, the total potential energy functional of the one-dimensional model under examination is the following:
+ f$U’) dS”’
- 1
k$,(ti(~*)(o))~(QK;." +@b”“>
S,, := 2QaEx S,s := QISEX+ 2QsEu
- A ,t, (s(~'(L'"'))'(Q~' + ;Qf.“,. (20)
Q,j being the material
stiffnesses in the XI, X2, X3 system. The equations rcelating the stiffness coefficients to the engineering elastic constants (Young’s moduli, Poisson’s ratios and shear moduli) in the materialsymmetry system, and their transformation in the structural system .I’,, Xz, X,, will not be given here, but they can be fc’und in Ref. [5]. With the aim of’ including the case of composite bimodular materials in the theory under examination, we follow Pate1 et al. [6] and Bert [7] assuming that the material possesses one set of elastic moduli when the fibers are elongated (EL = EA, . AL >, 0, AL being the unit vector in the fiber direction) and a different set of moduli when the fibers are shortened (EL < 0). We hence write: Qti = U - H(Et)lQ,
+ HWQ;,
(21)
H(EL) being the Heavyside function:
H(EL;I =z
0, if EL < 0 1, ifEL>O‘
(22)
(24)
It is defined over the set of kinematically admissible displacement fields of the laminated beam. The assumption of perfect bond between adjacent sublaminates requires that such displacement fields satisfy the constraint eqn (6), which can be rewritten as it follows: #+ ,
=
uf
+ I’-
on
a@’
,e-, a&k
VkE{1,2,3 where
+ I’
,...,
N-l},
(25)
F. Fraternali and G. Bilotti
844
The total potential
energy of the penalty model
is: N-l
4-I+; c
+I’-
_
Ul”
k=l
@t
I’-
=
Vi”+ ‘1 _ [@ + ‘1 _
+
):I!;,; “X;l
f#,‘l” + “&k
+ “]_y,
(27)
It is evident that the explicitation of the perfect bond conditions (25) in the functional fI is rather complicate from an analytical point of view. It is hence helpful to enforce such conditions via a penalty technique [16]-[18]. This leads us to introduce an auxiliary model (penalty model), which differs from the original one because the sublaminates are connected by continuous distributions of elastic springs instead of being perfectly bonded (Fig. 2).
where A”) is the interface between @) and B/2+‘) [kk’ = aBik”’ n ~3%~ + “1, l/q I”)is the stiffness (per unit area) of the springs acting at Nk) in the direction perpendicular to the beam plane, l/d) is the stiffness of the springs in the radial direction and I/# is the stiffness of the springs in the axial direction (Fig. 2). It is clear the sublaminate displacement fields a(‘), fi’2’ . , W, are completely independent in the ftuktional I&. When the quantities t~{~‘, nlk) and qsk’ (penalty parameters) approach 0, i.e. when the spring stiffnesses increase indefinitely, it is easy to conjecture that the penalty model approximates the original one [16H18]. Furthermore, as ~1~)-+ 0, ~4~)+ 0, 91”)--t 0, the spring reactions:
(i = 1,2,3; i not summed)
(29)
approach the values of the sublaminate interactions in the original model (interlaminar stresses). The equilibrium equations of the penalty model can be derived by means of the principle of stationary
Fig. 2. Penalty model.
Nonlinear elastic stress analysis potential
energy [19]. Such a principle
3. FINITE ELEMENT APPROXIMATION
ensures that
the couple (ii, A), where
Let us introduce a finite element discretization for the axis curve C of the laminated beam with n, elements C,, C2, . . . , Cme. We suppose that such elements belong to the Lagrange’s family [ 16, 181and are of the isoparametric type. We hence adopt the same shape functions to approximate both the geometry and the displacement field within the generic elements C, (see Fig. 3)
(yl’
1.1 @2’
it=
(30)
&v
characterizes an equilibrium admissible and it results:
state if and only if, fi is
a
aIIp=gI’(Q+crM,1)
&=o
(31)
In eqn (33) 5 is the normalized coordinate on the undistorted (or master) element (- 1 < r < 1); n is the number of nodes connected to the element; Z is the position vector of the generic point of C,, Z, and 6, are the nodal values of the fields Z and B, respectively; N,(c) (I = 1,2, . . . , N) are the Lagrange’s interpolating polynomials of the order n- 1. The Jacobians .P = aS”)jag can be computed from the equations
for all admissible variations 611. 2.5. Stability criterion We adopt the energetic stability criterion of Trefftz [19] to compute the stability points of the laminated beam. Referring to the penalty model, the Trefftz’s criterion states that the couple (il,, &) characterizes a critical equilibrium state (that is a stability point) if, and only if, there exists some V # 0 (buckling mode) such that it results:
6(62rlP)
=
&
l-I,@ I
+
a,?
+
2
845
=
a2M,n,)
pk’=
(, -$)($)
0
1,=z,=o
(32) for all admissible variations 6B.
1
I
n
g-1 Fig.
3. Finite element mesh and detail of the generic element.
846
F. Fraternali and G. Bilotti
where Z,, and Z,, are the components of Z, with respect to a fixed Cartesian coordinate frame (see Fig. 3), while N,., and NJ.{
denotes the contribution to the stiffness due to the springs acting at the sublaminates interfaces. Following a standard procedure, we build up the matrices Ku and Kp by assembling the contributions due to the individual elements of the finite element mesh: KU<,K, (e = 1, 2, . . , n,). In particular we compute the matrices K, by peforming an ulterior
(35)
3.1. Path-following procedure Utilizing eqn (33), it is possible to reduce the equilibrium equations of the penalty model [eqn (31)] into the following discretized form: SU’{K(U)U - n[Q’” + Q”‘(U)]} = 0 where U is the (secant) nodal force respectively. follows:
(36)
the nodal displacement vector; K is stiffness matrix; Q(l) and Q”) are the vectors of the first-and second-order, It is useful to represent U as it
assembly of the sublaminate contributions K’:j, K$;, . . . , KIN. Such contributions will not be given here, since they can be readily obtained by applying to each sublaminate the finite element model presented in Ref. [6]. In a similar way we build up the global force vectors Q(l) and Q”). Concerning instead the element penalty matrices K,, we refer to the Appendix. Coming back now to eqn (36), we remark that it must hold for each admissible variation HUE@’ (M is the total number of equations). Hence, it is equivalent to the following system of nonlinear algebraic equations: R(U, 1) = K(U)U - L[Q”’ + Q”‘(U)] = 0.
(40)
It is helpful to solve eqn (40) in the context of a displacement control procedure. This can be done by adding to eqn (40) an extra equation, which forces a particular component of U (say the pth one) to assume a prescribed value p (see [21], [6]). We hence introduce the following extended system of equations:
(41)
n, being the total number element mesh
of nodes of the finite
nti=(n,+l)xn+l.
(38)
The stiffness matrix K can be expressed as K=Ku+Kp
(39)
where Ku denotes the stiffness matrix for the unconstrained problem (i.e. for the disconnected sublaminates), while Kp (penalty stiffness matrix)
where e, is the unit vector which has the pth entry equal to 1 and all the other entries equal to 0. An efficient solution scheme for the matrix eqn (41) can be obtained by employing the well known Newton-Raphson method. It consists of a partitioning method also called bordering algorithm in the literature (see Refs [22,23]). For the description of such an algorithm we refer to the work [6]. 3.2. Procedure for the computation of stability points The variational eqn (3l), which rules the problem of computing the stability points of the penalty model
Nonlinear elastic stress analysis and the corresponding buckling modes, admits the following discretized formulation: &F&&f,
n)v =
v
(42)
where V is the nodlal displacement vector corresponding with the buckling mode, while KT is the tangent stiffness matrix: K1V=~R(U+aV,l)i,=c’fV~4.
(431
We decompose KT as follows KT = I& + K,
(44)
where K,, is the tangent stiffness matrix unconstrained problem i.e. KT,V = $
where V0 denotes a prescribed value of the pth component of V. It is convenient to let the position p coincide with the equation number where the smallest diagona! term of the (factored) tangent stiffness matrix appears. Concerning the scaling factor Vo, we set
(49) where V. is the initial approximation to V. Since the matrix KT becomes progressively ill-conditioned as (U, ,I) approaches a stability point, it is numerically convenient to transform system (48) into the following form [24]
of the
{Ku(U + aV)W + aV]
- .IQi2’(U + uV>]
i1=0
vv EP
(45)
while KP” is the tangent stiffness matrix due to the springs acting at the sublaminate interfaces:
KspV = $ (Kp(U + aV)[U + aV]> Ia-r! tlv E SF.
(46)
We build up KS,>by assembling the contributions Such contributions can be obtained by applying to the generic sublaminate the finite element model presented in Ref. [6]. Concerning instead KTO,we give in the Appendix the expression of the contribution due to the generic element of the fide element mesh. Awaiting the arbitrariness of SU, eqn (42) is equivalent to the following system of algebraic equations due to the differenl sublaminates.
;Kf(U, a)v = 0.
(47)
The direct computation of stability points can be performed by adding to (47) the equilibrium eqns (41) and an ulterior equation preventing that V becomes the null vector [24, 61. We hence introduce the following extended system
R**(u,
841
v, A, p)
=
=o
(48)
y being a real number greater than zero. A solution scheme for eqn (50) can be obtained by means of a partitioning method similar to that employed for eqn (41) [24,6]. The numerical procedures for eqns (41), (50) can be combined to formulate an overall algorithm for stability analysis [6]. Within such an algorithm, the path-following method is used to follow a solution branch (e.g. the fundamental equilibrium path). Near a stability point the extended system (50) can be instead applied to compute this point. To estimate how to switch between the two procedures, it is helpful to check for the sign of det Kr along the solution branch. Indeed a change in the sign of det KT indicates that a singular point of KT (that is a stability point) has been passed. Once such a point (U,, A,.)has been computed, it is possible to follow the post-critical branch by employing again the path-following method. In particular, if the stability point is a bifurcation point (VTQ = 0), the switching of the path-following procedure on the bifurcated path can be performed by adding to U,. a vector which is proportional to the buckling mode V:
(51) Here { denotes a scaling factor to be determined by imposing that the couple (U, A<) statisfies the equilibrium eqn (41). 4. NUMERICAL RESULTS
We present in this section the results of a finite element study on stress distribution in straight and
F. Fraternali and G. Bilotti
848
I
L/z :
rigid dars
//
h=O.OSL
_ w--v
unimodulor bimodulor
v/h Fig. 4. Loadxenter
deflection curves of a straight laminated beam (D& = Jx Q&g dZ: bending stiffness of the unimodular beam).
curved laminated beams. The adoption of the finite element mode1 formerly described, allowed us to extend such a study beyond the linear-elastic field, taking into account moderately large rotations and shear strains (geometrical nonlinearity) and a different elastic behaviour of the material in tension and in compression (mechanical nonlinearity). Two examples were studied: a three-layer straight beam (example 1) and a three-layer shallow arch (example 2). Both of them consider a cross-ply lamination scheme ([O/90/0]), a central point force and hinged supports. To obtain accurate stress distributions, we let the number of sublaminates coincide with the number of individual layers (N = 3). In each sublaminate we modelled the cross-section warping through the following function:
iv(k)= w&:‘(#‘)2 + w$(%;))’
(k = 1, 2, 3).
us to avoid the numerical incovenient, which is known as “locking” of the finite solution in the literature [ 16, 181. Concerning the calculation of the stress field, we employed different procedures for tangential and normal stresses. In detail, on a cross-section corresponding to a Gauss point, we computed the through-the-thickness distribution of shear stress S,, by interpolating the interlaminar values (A$‘),Al’)) through a layer-wise quadratic shape function. Such
(52)
In both of the examples, we used four-node isoparametric finite elements and a Gauss formula with 4 points to compute, on each element, either the “unconstrained” stiffness matrix Ku. and the penalty stiffness matrix Kp,. In this last one we set
$k’
=
#’
=
g’
=
IO-3
1
maxWu,, 1
(k = 193,
(53)
where max{&,,) is the largest of the (translational) diagonal elements of Ku. The above positions allowed
Fig. 5. Deformation of the terminal part of the beam analysed in the example 1 (unimodular material behaviour, FL'/48D&h = 0.203, displacement magnification ratio = 5).
Nonlinear elastic stress analysis
a function was calibrated librium conditions:
by imposing
s
&3 dEk’ = Nj) X(k)
849
the equi-
(k = 1,2,3).
3j@’ +bh’k’-
(54)
We hence obtained the following distribution:
,$k)+ $ - I) 4
for
2015-
lo-
5-
upper interface (O/90) unimodular behoviour free edges
I
-1
I,
x-0.001
O-
-5.
PI I31
-10.
h
h-0.05L
-15.
I
I
I
I
I
I
I
0.002
0.000
0.004 x3
/
I
I 0.008
0.006
OA3’IO
L
lower interface
(90/O) unimodulor behaviour free edges
X=0.466
IA HI
lJ2 :
k
; 0.000
-20
I
I
0.002
I
1
I
0.004 x3
I
0.006 /
L/2
1
I
h
h-0.05L
I
0.008
I 0.t3’IO
L
Fig. 6. Interlaminar distributions of the transverse shear stress for the beam analysed in the example I.
850
F. Fratemali and G. Bilotti
We instead computed the distribution stress Sa from the constitutive equations,
of normal obtaining:
The elastic properties
of the material
were taken
as
& = Q$)[(_$$" + $_@*2') + (@J) + f@$".")~/
8, =rnl?t,b;=
IA \
x s
-lO-15-
t
:
I-/*
L/2
h=O.OSL
I
upper interface (O/90) unimodulor behoviour free edges
-2o-
1
-25
I
I
0.000
0.002
I
I 0.004 x3
I
I
0.006 /
I
1
0.008
O.( 10
L
(b)
30
I
251
lower interface (90/O) unimodulor behoviour free edges
I4 IEII
h
h=O.OSL
5-j 0
-5 -10
-15
\
x=0.001
Ii_
_ .
_
_ _ _
_ _ _
_
x-o.203 X-O.468
!
0.000
Fig. 7. Interlaminar
I
I
0.002
distributions
I
I
0.004
I
I
0.006
I
I
0.008
I
O.( IO
of the transverse normal stress for the beam analysed in the example 1.
Nonlinear elastic stress analysis
Here bL, br, gN are Young’s moduii in the materialsymmetry system (L, T, N); YLT, YLN,YTN are shear moduli; vLr, vLN,i’rN are Poisson’s ratios; m is the bimodularity ratio (m = 6, I&); the apex + (-) marks values in tension (compression). The bimodular - behaviour that we considered corresponds to m = 0.5 (material weaker in compression than in tension).
0.8
Example 1. A straight, [O/90/0] beam (/z/L = O-05), simply-supported at the ends and subjected to a central point force, was analysed. Owing to symmetry, only one half of the beam was modelled, with a 25 element mesh. In detail, eight elements were used to discretize the part of the beam which starts from one base and spreads out for 0.025 for the span L; five, four and eight elements were
(=I
0.6
‘0
,*’
,’
0.,2
-0.50
8.0
I
-0.33
lb)
6.0 -
\
4.0
I
-
A-0.001 0.205 0.466
_____ ___
I
/ I
I
I
I
I
I
I fme Bdgt3B l unimodular
I
I
brhoviour
-0.17
0.17
-0.00 /
l I
I
x2
IL
I
’ A@ 0.
,I=’
0.4
0.0
851
---x-i
IO
t-l
I
I
I I
I I
X,=O.O42L
I‘,
-
1:
I free edges
I unimodulor behoviour i0
Fig. 8. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analysed in the example I, when the edges are free to warp and the material is unimodular.
F. Fraternali and G. Bilotti
852
instead used to discretize the successive parts of sizes 0.050 L, 0.10 L and 0.325 L, respectively. Two different kinematical behaviours of the bases were considered: bases free to warp (free edges) and bases constrained to rotate rigidly (rigid ends). In addition, either a unimodular and a bimodular behaviour of the material was taken into account. The
load-deflection curves corresponding to the examined cases are shown in Fig. 4. Figure 5 shows the deformation of the terminal part of the beam when the material is unimodular and the bases are free to warp. A marked warping of the central layer can be observed. It is due to the strong difference existing between the axial Young’s moduli
(a)
0.8
0.6 LL \ $: v,
0.4
5 0.2 1 unimodulor behoviour 1 0.0
I -0.33
-0.50
I -0.00
-0.17 x2
15
(b)
‘O-
~-0.001 _____ 0.205 0.466 __-
/
0.50
h
I
I
I
I
I
I
h=O.OSZL
*--------I *s----1___
I I
5-
I 0.33
0.17
---*_ --mm,
t
I _____---
-----
I
-5-
I -10
-0.50
I -0.33
I
I
I
-0.17
0.17
-0.00 x2
/
rig-d edge8 unimodular
behoviour
I
0.33
0. 10
h
Fig. 9. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analyzed in the example I, when the edges are riBid and the material is unimodular.
Nonlinear elastic stress analysis of the edge (0) tayers (Q, = gL) and of thecentral (90) layer (8, = c?~= E’J25). The above-menlioned through-the-thickness nonhomogeneity of material properties is also responsible for a marked stress concentration near the bases of the beam, when the same are not stiffened (free edges). Here the interlaminar stresses appear
853
to be tmbonded, as it is shown in Figs 6 and 7. These figures contain plots of the axial distributions of the transverse shear stress & and of the transverse normal stress & at the interfaces between the edge (90) layers and the central (0) layer, for several values of the load parameter 1 = FL3/(48D&h).
(a) 0.66
-
-
X-O.001
IL
\
0.4.6
x 2
0.26
0.06
-0.14
-0.50
I
I -0.33
I -0.00
-0.17 x2
8.0
6.0
IL
-
-
_____ ‘=;.;i: __c 01460
I
2.0 -
tr \ I ‘. \
,I
2 n
I
fl_V’ (_=-=-~-%_W~_ --_---
-2.0
-
-4.0 -0.50
I -0.33
X,-0.042L
‘\ I\
, I I I / ,’1 I ,’ I I #) I I’ I 1 ,’
4.0 -
0. i0
h
I
\ 2
/
I 0.33
I 0.17
‘. \ ‘. t ‘. \ l.\
----
I
I
I
I
I
I
free edges bimodulor
khoviour
&-=O.**
I
I 0.17
I -0.00
-0.17 x2
/
I 0.33
0.50
h
Fig. 10. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the laminated beam analysed in the example I, when the edges are free to warp and the material is bimodular.
CAS 62/5-C
F. Fraternali and G. Bilotti
854
It can be observed that the given shear stress distributions are generally unsymmetric. Such result, which is due to the above-mentioned edge stress perturbation, is more evident as the load increases, the edges are free to warp and the material is bimodular. In the same circumstances, the relative ti$X$)‘) “weight” of warping terms (ti$‘@‘, increases in eqn (56) and the axial stress distributions appear to be strongly nonlinear. Moreover, in the bimodular case, such distributions exhibit a marked slope discontinuity in correspondence to the zero (Fig. 10).
Figure 6 shows that the transverse shear stress exhibits edge singularities in correspondence to the interlaminar surfaces (O/90) and (90/O). Such singularities are opposite in sign and their powers increase with the load. Even the transverse normal stress edge singularities are in correspondence to the interlaminar surfaces (see Fig. 7). However, the power of the singularity at the upper surface decreases with the load. The power of the other singularity instead increases. The marked oscillation about 0 of the numerical results presented in Figs 6 and 7 is due to the peculiar nature of the penalty model, in which the layers of the beam are connected by means of very stiff springs. Free edge stress concentration is common in composite laminates and well recognized in the technical literature [9-141. It can produce a failure of the composite structure due to delamination. For this reason, methods of reducing the edge stress concentration have been proposed. They essentially consist of reinforcing the stiffness at the free edges, for instance by means of “cap” (a lamina or a collection of laminae zrapped around the edge) [ 13-141. We now pass to analyse the stress distribution within the layers of the beam. Through-the-thickness distributions of the shear stress & and of the axial normal stress S,, were computed at the Gauss point X, = 0.042 L (Figs 8-10). Distinction was made between the free-edge (Fig. 8) and the rigid-edge conditions (Fig. 9). In the former case, the effect of a bimodular behaviour of the material was also considered (Fig. 10).
2.0
Example 2. The second example that we considered refers to a shallow, [O/90/0] arch, which is simply supported at the ends and subjected to a central point force. The axis-line was supposd to be parabolic and its aspect ratio f/l (J central raise, I: base length) was taken as l/20. The thickness-tolength ratio h/l was instead taken as l/200. A uniform finite element mesh composed of 25 elements was employed to discretize the axis-line. Figure 11 shows the load-deflection curves of the considered arch, for a unimodular (8; = S,t) and a bimodular (8; = 0.5 8:) material behaviour (the deflection is referred to the point P shown in the figure). A bifurcation of the fundamental equilibrium path (with a skew buckling mode), a marked influence of the material bimodularity on the critical load and a slightly unstable asymptotic behaviour of the arch along the post-critical path can be observed. Due to the relevant slenderness of the arch, we
I
1.5
1 .o
h=O.O05t f30.051
0.5 - - - - bimodulor -
unimodular
0.0
Fig. I I. Load-displacement
curves of a shallow, [O/90/0] arch, for a unimodular (8, = Bt ) and a bimodular (8, = 0.5 8:) material behaviour.
Nonlinear elastic stress analysis
practically obtained the same results considering the bases as rigid or free to warp. This led us to omit a specific study on ledge effects in the present case. We instead computed the through-the-thickness distributions of the transverse shear stress SX3and of the axial normal stress &, at the Gauss point X, = 0.3665 L, for the unimodular and the bimohular
855
interlaminar surfaces is apparent when the material behaviour is of the bimodular-type. As in the previous example, a slope discontinuity of the normal stress distribution in correspondence to the zero is either apparent in such a case (Fig. 13). 5. CONCLUSIONS
case (kegs 12, 13).
A
more
pronounced
stress
gradient
near
the
A
theory
and
finite
element
model
for
0.6 IL \
w 2
cn
0.4
0.2 unimodular
khoviour
0.0 -0.60
-0.33
-0.17
-0.00 X2
20.0
10.0
Jb)
-
OS.0__
-
-
_lC~.O _.___------_ *e--
-
-;7
Nd
P’
-0.50
I -0.33
I
I
I
I 1
I I
X,=O.J665L
P’
P’
I 4
I I
I
-
I
I
-__
I
khoviour
-0.17
/
I
A-0.001
--mm_
@OOW
0.472 0.959 1.732
t
I -0.00 X2
Fig. 12. Through-the-thickness
0.
h
I
I unimodulor -30.0
0.33
I
_~_-~-,-=. _-
9’ -20.0
-
/
0.17
I
0.17
0.33
0,
h
distributions of the transverse shear stress and of the axial normal stress for the arch analysed in the example 2, when the material is unimodular.
the
F. Fraternali and G. Bilotti
856
non-linear elastic analysis of laminated curved beams has been proposed. Even remaining in a one-dimensional context, such a theory allows for an accurate description of the interlaminar stress field and gives rise to the possibility to compute the stress field within the laminae of the structure. The considered nonlinear effects are both geometrical and mechanical. The first refer to moderately large rotations
0.6
;“’
and shear strains; the second to a different elastic behaviour of the material in tension and in compression. The proposed finite-element model allows one to follow the fundamental equilibrium path of the beam, to compute its stability points and to analyse the post-critical behaviour. It has been applied to a numerical study on the stress distribution in straight
I
I XS=C.3665L
I 0.6 \
\
0.4
bimodular
0.0 -0.50
I -0.33
I -0.17
behoviour
I -0.00
1 I
0.17
0.33
0.
Xz / h 20.0
(b)
10.0
0.0
t.z 5
I
I
I
I
I
X,=O.3665L
, / /
I
LL \
I
-10.0
I
--______
__-- ___--__--- _.--4 .s( /@ /cc I
I
I
I
I I
cc
-20.0
A=O.OOl _____ 0.472 _--
0.956
I bimodular behoviour I c-=0.5&
-30.0
I
-0.50
-0.33
I
I
I -0.17
-0.00 X2
/
0.17
I 0.33
0
h
Fig. 13. Through-the-thickness distributions of the transverse shear stress and of the axial normal stress for the arch analysed in the example 2, when the material is bimodular.
Nonlinear elastic stress analysis and curved laminated beams. The presented results show the potentia:l of the adopted theory in analysing local effects, such as free edge stress concentration
and a bimodular behaviour of the material. In future works, an extension of the present theory will be formulated,
in order
to study
delamination
and
failure of composite beams.
REFERENCES 1. J. N. Reddy,
A generalization of two-dimensional theories of laminated composite plates. Commun. Appl.
Numer. Meth. 38, 173-180 (1987).
2. J. N. Reddy, E. J. Barber0 and J. L. Teply, A plate bending element based on a generalized laminate plate theory. ht. J. Numer. Meth. Engng 28, 2275-2292
Nomizu and S. Kobayashi, Foundations of Differential Geometry. Wiley, New York (1963). 21. I. L. Bathoz and G. Dhatt, Incremental displacement algorithms for nonlinear problems. Inr. J. Numer. Meth. Engng 14, 1262-1267 (1979).
20. K.
22. E. Riks, The application of Newton’s method to the problem of elastic stability. J. Appl. Me&. 39, 1060-1066 (1972). 23. H. B. Keller, Numerical solution of bifurcation and nonlinear eigenvalue problems. In: Applications of Bijiircation - Theory _ (Edited by . _P. H. Rabinowitz), pp. 359-384. Academic Press, New York (1977). 24. J. C. Simo and P. Wriggers, A general procedure for the direct computation of turning and bifurcation points. Int. J. Numer. Meth. Engng 30, 155-176 (1990).
(1989).
3. L. Ascione and F. Fraternali, A penalty model for the analysis of laminated curved beams. Cornput. Strucr. 45(5/6) 985-999 (1992).
4. F. Fraternali and J. N. Reddy, A penalty model for the analysis of laminated composite shells. Int. J. Solids Struct. 30(24). 3337-3355
(1993).
5. L. Ascione &d F. Fraternali,’ A moderate rotation theory of laminated composite curved beams. Znr. J. Engng Analys. Des. I(2), 161-176 (1994).
6. H. P. Pate], J. L. Turner and J. D. Walter, Radial tire cord-rubber composites. Rubber Chem. Technol. 49, 1095-I 110 ( 1976). 7. C. W. Bert, Models for fibrous composites with different properties in tension and in compression. J. Emma Mater. 344-349 $977).
Technol.
Trans.
ASME
857
APPENDIX Penalty stiffness matrices
In Section 3 we introduced the penalty stiffness matrix KP to distinguish the contribution to the global stiffness matrix, which is due to the springs connecting the sublaminates of the beam. We now consider the contributions to KP due to the generic element of the finite element mesh (say the eth one). It is convenient to order the nodal displacement vector of such an element as it follows:
994(4). ~ ,_
(Al)
8. L. Ascione and F. Fraternali, A finite element analysis of the stability of bimodular composite curved beams. Znt. J. Engng Analys. Des. l(3), 315-334 (1994).
9. N. J. Pagano, On the calculation of interlaminar stresses in composite laminates. J. Comput. Mater. 8, 65-82
where
(1974).
10. N. J. Pagano, #Stress fields in composite laminates. Int. J. Solicis Struct. 14, 385400
(1978).
11. A. S. D. Wang and F. W. Crossman, Some new results on edge effects in symmetric composite laminates. J. Compur. Mater. 11, 92-106 (1977). 12. J. D. Whitcomb, I. S. Raju and J. G. Goree, Reliability of the finite-element method for calculating free-edge stress in compor,ite laminates. Comput. Struct. 15(l), 23-37 (1982). 13. J. N. Reddy, Reduction of free-edge stress concentration. J. Appl. Mech. 52, 801-805 (1985). 14. R. Y. Kim, A technique for prevention of delamination. Mech. Comput. Reo. Oct. 28-30, Dayton, Ohio (1981). 15. T. Argyris, An excursion into large rotations. Comput. (1982).
Methods
Appl.
Mech.
Engng
32,
642)
The matrix Kp
85-155
N. Reddy, An Introduction to the Finite Element Method. 2nd edn. McGraw-Hill. New York (1992). 17. J. N. Reddy, Applied Functional Analysis and Variational Methods in Enaineerina. .- Krieaer, Melbourne.
16. J.
FL (1992). 18. G. F. Carey and J. T. Oden, Finite Elements: A Second Course. Prentice-Hall, Englewood Cliffs, NJ (1983). 19. B. Budiansky, Theory of buckling and post-buckling behaviour ot’ elastic structures. In-Advances in Applie; Mechanics (Edited by Chia-Shun Yih), pp. l-65. Academic Press, New York (1974).
where K$h’ =
(wp’ + wp’+)%“)(wp+
(A4)
F. Fraternali and G. Bilotti
858
We now pass to consider the tangent, penalty stiffness matrix KT?,i.e. the matrix defined by eqn (46) of Section 3. Then the contribution to KT~due to the generic element of the finite element mesh is given by the equation:
(AW
x L”‘(W$‘+ + fW(,“‘+); d< NY)P + ‘)d[
(A6)
where
_
(A7)
In eqns (A4)-(A7) we denoted by NP’ and Nf + ‘) the diagonal matrices which have non-zero terms equal to N,(t). In the same equation we put:
KY;:+“=
b & NF+“J’A.” x (Wb”+ I,- + Wp + I,-) _ +I, & 2
(A14)
W#k’+.‘k+ I)-, =
a
=
&I,
&
+ 1,;
,j
1000
0
0100 0 0
0 -X,
=
b’“’
I
X,
b” + I,
-x2
0
0
X,0 0 0 x;2 x,x*
“’
0
‘.’
0
o...o,..o Xi .”
0
x;
“’
x;
x
II1,
1
(Wb”‘W
+wz”+“-)7p
(A9)
Nonlinear elastic stress analysis
The matrices I”‘), PC+“, which appear in eqns (A13) and (Al6), are given
859
by:
i 00 0 0
pw.cr + Ill
=
0
.’
0
”
0
”
0
‘.
0
.’
0 0 0
0
”
0 0
0
10 o0 0o / o
. .
.
0 . .
0 0 0
‘. .
0
0
fj
(Al7)
,.
.
0
,I\“‘,Ask)and ,I\“’being the reactions of the springs interposed between the kth and the (k + I)th sublaminate [see eqn (29) of Section 21.