Modeling laminated composite structures as assemblage of sublaminates

oo_xL768388 13.00+ .w IV88 Pngwnon Pms plc

MODELING LAMINATED STRUCTURES AS ASSEMBLAGE

COMPOSITE OF SUBLAMiNATES

S. N. CHATTERJEEand V. RAMNATH MaterialsSciences Corporution. Gwynedd Plaza II. Bethlehem Pike. Spring House. PA lYJ77. U.S.A.

Abstract-A mixed variational principle is presented for stress analyses of laminated composite structures with or without delamin;ctions and discontinuitics. Of particular interest is the use of a functional involving interlaminar stresses and in-plane strains. instead of strain or complementary energy. A laminate is then treated as an assemblage of sublaminates with assumed variations of displacements and interlaminar stresses through the thickness of each sublaminate. Application of the variational principle. at this stage. reduces the problem to IL set of dilferenGal equations and constitutive relations for each sublaminate similar to those encountered in higher order plate theories. as well as continuity conditions of displacements and equilibrium of surface tractions at pcrfeclly handed intcrf&s. It is dcmonstratcd with an example of ;L dishonded luminutc modeled as an ;issembl;lge of two subl;imin;l(cs that Ihe mixed boundary conditions at P partly bonded interface yield conccntratcd forces aI the disband lips. This is the exact form of stress singularity in this model. Finite clement analysts arc utili/cd for obtaining solutions with more sublaminates. Stress tieIds and energy release rates at dcktmination tips arc compared with rlnsticity and standard Cnite elcnicnt solutions in anofhcr illuslrativc cxamplc.

IN-l-ROIXJCTION

Growth

of disbands

created due to impcrfcct bonding or foreign object impact can scvcrcly

inllucncc the load bearing capacity of Inminatcd composites. For this reason. delamination fracture has been thcsubjcct ofmany investigations. and singular

Mcthodscmploying

integral equations discussed in scvcral studics(l-31

Fourier

transforms

are suitable for stress

analygis when the disbands arc located away from cdgcs or other discontinuitics.

In many

casts. such disbonds can. however. occur near edges or other geometric discontinuities, may in fact originate at such locations. joints

arc CxamplCs. Rigorous

elasticity

and

Free cdgcs in Iaminatcs. ply drops, and bonded solutions

and numerical

problems arc impossible from practical considerations.

Standard

calculations

for such

finite clement methods

(or those employing c‘lcmcnts with singular liclds) have been utilized in such problcms[4,5]. In such currently available mod&

with assumed displacrment ticlds, it is nczcssary to model

each layer ofdifl’erent

scparittcly. In many GISL‘S. howcvcr (SW results for elliptic

orientation

disbonds and those for two-dimensional an approximate

solution

based on

(2-D)

problcms[Z,

the assumption

sublaminatcs which deform according to shear deformation rclcasc rates close to elasticity solutions. thcsc approximate

solutions

This

3. 6. 71 for large disbonds),

that disbonds

arc located between

plate theory yields strain energy

is true even though the stress singularities

appear ;IS conccntratrd

forces at disbond tips[6-81.

This

in is

especially true when the disbond dimensions are large compared to the thickness and the plate models yield a very good approximation. Since disbonds of comparatively small sizes arc usually not critical. it appears that for all practical purposes it may bc unncccssary to USC standard tinitc elcmcnt tcchniqucs modcling each layer. For more accurate solutions, Iaminatcd plate thcorics of higher order. or a large number of sublaminates through the thickness of the luminatc with each sublaminatc still containing many layers, can possibly bc utilized. In higher order plate thcorics. howcvcr. calculation of plate properties based on assumed strain t&Ads dots not accurately rcprcscnt ctl?cts of some stiffnesses. For example. shear stiffnesscs arc more accurately rcprcsentcd by the assumption of constant stress fields. and correction factors having values close to those of isotropic or orthotropic plates[9, IO]. than by the assumption of constant strain fields. For these reasons. in this study a mixed variational principle with assumed linear variation of displacement fields and transverse shear stresses. and constant transverse normal stress through the thickness, is

S. N.

440

CHAITERJEE and V. 4

Lammae

’

R~sfN~rn

2. x3

I

mterfaces Lammae

Sub(ammates

Delammotw7s

An dement of sublaminate n Boundary curve C

Fig. I. A hminxtc

formulated

(as dcse&ed

divitluil inlo N suhlamin~~lcs

later) to obtain a higher order plate theory for caeh sublaminate,

which is applicably to the general three-din~~nsion~ll (3-D) di~ontinuities

(Fig.

1). Displacement

imposed in terms of Lagrangian

contin~lity

multipliers,

although

similar

in form

at bonded interfaces are

which are the tractions

Existence of concentrated line forces at disbond formulation,

conditions

problem with d&bonds and at those locations.

boundaries is also demonstrated.

to theories developcui by other investigators

Refs [I I, I?] for example) is more straightforward

This (see

for practical applications.

With the imposed condition that displacements of and tractions on adjacent layers on the

two sides of a bonded interface bc the same (opposite in sign for traction) the

Euler

equations for generalized plane strain (quasi-3-D or 2-D) problems can be reduced to a set of simultaneous ordinary differential

equations as in Refs [I I, I?]. Counting

forces at disbond or crack tips as unknowns,

the concentrated

there exists the required number of unknowns

to satisfy the boundary conditions or other required relations across disbond or crack tips. It may be noted that such conditions on displacements (not displacement derivatives) and moment or force equilibrium (not stresses) can only be satisfied. For this reason, concentrated forces exist at such locations as in similar

approximate

formulations[&-8,

131.

Such a system of equations for a finite (but not small) number of sublaminates is, however, not very simple

and the procedure for development of a gcncral algorithm

involved because of the possibility

can bc quite

that some of the characteristic roots of the system of

equations can be repeated[ I I, 121. Also, such an approach is useful only in linear problems. For these reasons, a finite element model is developed for quasi-3-D problems based on assumed linearly varying displacement fields in the direction parallel to laminations. String stiffnesses of rectangular elements to consider effects of large deformations, as well as change in shear stiffness to model non-linear shear response, which are important in bonded joint or other problems, can then be easily incorporated into the analysis.

Modeling laminated composite structures as assemblage of sublaminates MIXED VARIATIONAt

441

PRINCIPLE

Consider a laminated plate with K laminae as shown in Fig. 1. Let $, .$ and d, denote the displacements, strains and stresses in lamina k referred to the Cartesian coordinate system shown in the figure. Let A.,+’and I;’ (i = 1. 2, 3; k = 2, 3,. . . , K) be the surface tractions on laminae k- 1 and k, respectively, on surface At, parts of the interface k which are not delaminated and 8 be the displacement on the corresponding surface. For prescribed surface tractions ?? and/or displacements u,” on the boundary of such a laminated plate as well as body forces bi one can pose the following mixed variational principle : Variation (I) = 0

(1)

where functional I is given by

In the above equation, a/I takes values I I, 22. 12, 21. yB varies over 33. 13. 31.23, 32 and i and j vary from I to 3 ; repeated indices indicate summation. U, and ?p in the last two surface integrals are equal to the displacements of laminae the surfaces of which coincide with the bounding surface and unknown tractions on S., respectively. S, (or S,) denotes part of the bounding surface, where stresses (or displacements) are prescribed. u, are the direction cosines of the outward normal to the bounding surface. The functional Wk also depends on the material propcrtics of lamina k (elastic. non-linear elastic, or elastic-plastic total deformation theory) and is chosen to satisfy the constitutive relations such that dW’

---_ = 4

dWk

and

ac;,,

- -a&

= 0.

(3a,

b)

It may be noted that if a stress component is expressed by eqn (3a) the corresponding strain does not appear in eqn (3b). If one seeks condition (I), i.e. I to be stationary permitting the functions u;“.c$, (a/l = I I, 12.22). t$,(k = I,. . . , K), 1:‘. t?: (k = 2.3,. . . , K) and i;” to vary over appropriate domains, then the Euler equations for the variational problem are given by the relationships described by eqn (3a) and the following set of equations:

awk -q=

1/2(4,J+&)*

--a$- ’ = i.l+k, $

J,

=

i.,-‘,

O,,V,

a;k,J+b;L

=

f,*

=

0, 4/I = W&

on S,.

bi,Vj

=

f,‘,

+u”,.,) (k = 1.2,. . . , K) U, =

U,’ Ofl

S, (k = 2,3,...,K) (4a-h)

d-’

= $ = d,

l,‘k+A;k

= 0

(k = 2.3 ,...,

K on A,).

(4i.j)

It is clear, therefore, that if all variables are chosen to fulfill eqn (I) all equations of the boundary value problem under consideration are identically satisfied provided

.s$ = 1/2(&-t-r&);

yS = 33,13,31,23,

and 32.

(5)

However, in the next section it is shown how approximate solutions may be obtained by treating a laminated structure as an assemblage of sublaminates, choosing linear variations of displacements and interlaminar stresses through the thickness of each sublaminate which may contain laminae of different orientation or materials, and directly using the variational

s. N.

442

tiNAll

CHAlTEnJEEand V.

principle, eqns (I) and (2). The authors are not aware of any available mixed variational principle which allows such a choice of the field variables in an assemblage of sublaminates. each of which may be inhomogeneous. APPROXIMATE

FORMULATION

AS AN

ASSEMBLAGE

We subdivide the laminate into N sublaminates,

OF SUBLAMINATES

nth sublaminate

(n = 1,2,. . . , N)

note that % P n = K all perfectly bonded to one ( > n= I another and assume the following displacement fields. and interlaminar (Go,. o3 i, a,*) stress fields in each sublaminate. zn being the :-coordinate measured from the midplane of that laminate containing

Pm numbers of laminae

u:= fP(.lc,y) + GJ4y.L y)

(64

We also assume that the strain fields can be derived from displacement fields (6a) by eqns (4c) and (5) and that at all perfectly bonded interfaces inside a sublaminate A!” and I# are chosen to satisfy eqns (4d). (4~) (so that eqn (4j) is also satisfied because of eqns (6b)) and (4i). rcspcctivcly. WC denote the surface tractions on the surfaces of sublaminatc n at Z” = k/r,,/2 (cithcr prcscribcd in ?‘,” or unknown in i., or f0) as t:“. Further, the unknown surface tractions on other bounding surfaces (parallel to the c-axis) of sublaminate n where displacements are prescribed are assumed in the form given below. It may be noted that N, M, Q, R (unknown on the bounding contour) are the stress resultants used in higher order plate theorics[ 141. ”

Q

i= 1,2and

On

?!=7;.+!$k”“. n

n

(6~)

Contracted notations for interlaminar stresses are used in eqns (6b). In what follows, we will use similar contracted notions for other stresses and strains given below and omit superscript n denoting the sublaminate u

II

=

61.

El1 = &I*

2C12

=

El,

622

&22

2EI3

=

=

62. E2,

=

U,? E,)

= =

E59 2.523

CTJ

&J

=

Eb.

(7)

If the relationships between stresses and strain in the lamina k are known in terms of the (6 x 6) stillness matrix [Cl, e.g. (3, =

C,,E/+ C,&

6, =

C.j&,+C*Bql

(8)

where i,j vary over I, 2, and 3 and tx, fi vary over 4, 5, and 6, then the functional W (the superscript is omitted) in eqns (2) and (3) can be expressed as Hf = 1/2[C;,&,E,+ 2C::Ei0, - S.#qu,u,].

(9)

It should be noted that use of this functional permits the choice of interlaminar stresses

Modeling Iaminatcd composite structures as assemblage of sublaminates

443

(6b) irrespective of assumed displacement fields @a), which yields a better representation of interlaminar stresses than that obtained by assuming displacement field only, especially in the case of sublaminates containing laminae of different orientations or materials. The elements of the three (3 x 3) matrices C’, C” and S can be calculated from the elements of C by expressing Ui and E, in terms of si and o, (i = 1,2,3 ; Q = 4,5,6) using eqns (8). In many laminates the following relations hold and will be assumed here although one can keep the generality if desired :

G = O(i= 1,2,3;r

= 5,6)

Sror=0(x = 4;j? = 56).

and

(10)

Substitution of eqns (9) and (IO) and the assumed displacement and interlaminar stress fields (6a) and (6b), strain fields (4c) and (5), interface tractions 2:” (~tisfying eqns (4d), (4e),k=2.3 ,..., P,, and expressed as t:” at other locations) as well as I$” (6c) in eqn (2) and integration with respect to zn over the thickness h, yields a modified expression for I in terms of the new variables introduced above and q, the unknown displacements over perfectly bonded parts of the interfaces. If one now seeks I to be stationary, permitting 4. JI:, r,f”, Q to vary over appropriate domains, one obtains a set of Euler equations of the variational problem which are similar to those encountered in higher order plate theories[l4]. There are, however, some differences because of the assumption of displacements as well as intcrtaminnr stresses and the mixed boundary conditions considered here. For these reasons, we give below the final form of these equations which is suitable for practical applications. The intermediate steps are omitted here for brevity. but can bc found in Ref. [Is]. In the following equations N, M, Q, R are the stress resultants used in higher order theories. so, X, y’, y’, arc extensional strains, curvatures, and midplane and linearly varying shear strain components, respcctivciy (tho superscript n, denoting sublaminate n is omit&d). A is the plan area of the laminated structure, A, is the bonded part of the intcrfacc n and A,* is a part of the unbondcd portion where known displacements are prcscribcd. C: (or Cz) are parts of the boundary contour of sublaminate n where displaccmcnts (or stresses) are prescribed. Diflrentiul

equutions for suhlutninate n

N,.,+N3,*+1; +r; +f, = 0

+fz= 0

N,.,+N2,2+1: +r; Qt.,

M,q,+Ml,2+

Ri*l+&.l

+Qz,z+G

+t; 4-J; = 0

$(r;-ri)+g,-Q, -Na+

=0

$t; -ri)+y, f;.9*=

Conslilutive relations for suhlunriflate n

= 0

I

h,( I, :,,) dz,,,.

(11)

S. N. CHAITERIEE and V.

444

RAMNATH

where 4, = A,,+E,.&,~iL,

A ,J = EJEu

B,, = B,‘,+E,~F,JIEu,

B,d = F,,IE,,

D,, = D:/+FlaFjJE4dr

A,, = lIEdo

A:,, B:,, D;, =

C;,( I,:,, z;) d:, I

(13)

Strain component displacement E” I

=

K,

ufi.

~I=$LZ+$Z.~. Continuity of displacements 2.3.. . . , N)

relations for sublaminate =

Jl,,jr

8:

b=k.

=

UC,z

+&I

n

(no sum) i = 1.2

;

Yt=$,+u!L

und quilihrium

yh =$,.,;

oj’ tractions

h ;;-JIy-’ = q, u:‘“+ ;;“*y = I,,“” ‘- 11”

a=

1,2.

over bonded interfaces

t,-(n-‘)+t,+n = 0.

L

(14) (A,, n =

(15)

Prescribed di.~plucemm t corrditions u:“+j;*;

= u:+”

n = 1,2,...,Novcr

u*-” ,

u: = l,,‘”

n=2,3,...,N+I

N = 1,2....,N

A,+

overA,*

(l&l)

over C’:

*y = I(ly

(16b) Prescribed

traction conditions on C:

N,v, +N,v,-(t:‘+t;O) N,~,+N~v~--(t;~+t~~) M,v,+M,v~-$;O-t;O)

M,v, +hfzvz-

?qr;Lt;O) 2

QIV, +QIVI -(t:“+t;O) R,v, + R2vZ- ;(t+t;“)

(17)

Modelinglaminated

composite structures

b*,i;c* =

A set of equations similar unknown

stress resultants

(I 7) r,?“(C) and t,‘*(C)

as assemblage of sublaminates

i-:(k)

us

(18)

dz,.

to eqns (17) with superscript

* replaced by 0 relate the

on C: where eqns (16b) hold. In these equations, and in eqns

are line forces which act on laminate surfaces I,, = k/t,,/2 at contour

C. These are either prescribed or unknown.

Interactive

forces between two sublaminates

fall in the latter category. Such forces are expected at the contours which are the free edges of a laminate as discussed later. Similar contours

which are the delamination

line forces are contained in tp in eqns (I 1) at the

boundary curves (shown by dotted lines in Fig.

across which the displacement components Uiq tii are continuous, not necessarily so. These line forces are caused by the discontinuities gradients. and can be called the stress singularities[b&

1)

but their gradients are in these displacement

131 in the mixed boundary value

problem of bonded sublaminates (or plates) described by eqns (I I )-( IS). In the next section we will demonstrate

the existence of these forces from

the exact solution

of a sample

problem. In parts of area ,4. whcrc all the sublaminatcs cqns (I I) can bc rcduccd to a set 3(N+ terms of 3(N+

I) displaccmcnts at (Nf

exist. similar

equations arc obtained,

interface displaccmcnts. as wcil as equations. will bc more.

When no displaccmcnt conditions

arc prcscribcd on surfaces pcrpcndicuiar to the z-axis.

thcsc equations can bc solved with enough unknowns.

which togcthcr with the unknown

conccntratcd forces. can bc cvaluatcd to satisfy cqns (i6b)

and (17). as well as continuity

of u:, $: and force and moment ccluilihrium across dciamination l:or cxamplc, let us consider the problem to bc quasi-3-D the plant ,4.4 which I&.

passes

through the tip of a delamination

I. On the right-hand side 0fclA.

thcrc arc 3(N+

hand side which contains one delamination, the solution 12X+

of the second-or&r

I8 unknowns,

continuity

satisfying conditions conditions

ilt

boundaries.

and restrict our attention to as shown on the top half of

I) displaccmcnt variables. On the left-

thcrc are 3(N-t2) dilfercntial

dispiaccmcnt variables. In

equations,

there will be a total of

arc available for satisfying

displacement

conditions across thr plane AA and the rest are to be used for

ends .\:= constant away from AA. At AA. 3(N+2)

(three at cJch of N+2

of N sublaminatcs) unknowns

ordinary

half of which (i.e. 6N+9)

and equilibrium

equations in

I) surfaces by making use of eqns (IS) as well as

relations (12). In other parts, whcrc dslaminations but the number of unknown

arc perfectly bonded, the differential

I) second-order partial differential

interfaces) and 6N equilibrium

have to br satislicd.

which come from

conditions

In general, one therefore

three concentrated

forces at AA

displacement (six for each

needs 3(N-

I) extra

in s, y, : direction

at

each intcrfacc other than the top and bottom surfaces of the whole laminate. An exact mathematical proof of the cxistcncc of such forces in the gcncral probicm is possible but complicated. Wc will give this proof for a simple problem considered in the next section involving two sublaminates. It is clear, however, that these forces that displaccmcnts arc continuous

WC

caused due to the fact

across A/I, but their gradients arc not. The cxistcnce of

such forces at intcrfJccs other than those containing delamination tips may seem impossible at first glance but is a natural conscqucncc of such discrctization

schemes and assumed

ficlds (similar phcnomcna arc obscrvcd in linitc clcmcnt solutions).

However.

the forces at

intcrfaccs other than the ones containing dclaminations. will be small and will rcducc with increasing N. It is shown in the double cantilcvcr beam problem considered later that even the force at the delamination tip rcduccs as the number of sublaminates, N. is increased. and the stress licid approaches the exact stress solution. As shown in the examples. energy rcleasc rates at delamination tips can. however. be evaluated with sufficient accuracy even with a small number of sublaminates. Pagano[ I I, 121 obtained a similar approximate formulation for homogeneous sublaminatcs. and solved the 2-D free edge problem ofa laminated composite. In his formulation,

S. N. CHAT-I-WFE and V. RAHNA~

4%

+

E!astlclty (I ) Constant stroln flnute etements (2) Higher order theory for four sublaminates shown by dotted lines (41 Closed formsotutlon (31

x 0

0

Method

(1)

(2)

(3)

(41

(~~

2 48

2.39

2.42

234

Eo =6 895 GPO

-100 -075-050

-025

025

0

050

075

loo

x/a

LocatIon,

Fig. 1. Results for unit normal stress n on delamination surfaxs in ;I laminate.

however, the ~~?n~~r~tr~~t~d forces art’ absent. and the intCrfiIcc strcsscs iit the edges (which art‘ assumed finitc) appcx to diverge IS the number (~l.subt~lrnini~t~s is incrcascd. Even for such 2-D problems. the solution

algorithm

bccomcs quite cumbersome as the number of

subktmin;ttcs is incrcasod. Morcovcr, ;I major objcctivcof this invcstigittion the influcncc of large deformation bonded joints.

For thcsc reasons. iI linitc clcmcnt solution of such 2-D problems (or quasi-

3-D) since three dispkcmcnts of y is formulntcd pondcnt ofy.

was todctcrminc

;~nd non-linear shear rcsponsc of adhcsivcs in single klp

arc retained, but all quantities arc i\sstmlcd to bc indcpcndcnt

its discusscc! later.

Other quasi-3-D

problems whcrc strcsscs arc indc-

but the displnccmcnts i\rc not ncccssarily so

it laminntc which is long and loudcd in the _4rcction)

STRESS

SlN~iULA~l’rY

IN

I’fIE

(as

in the free cdgc problem of

may br obtained by superposition.

FORM

Of- f:ORCES

Existcncc of conccntratcd forces at disbond tips hits been dcmonstrntcd plate models employing shcur deformation dcmonstratc

this with

the exilct solution

for bonded

thcory[X. 131. For the prcscnt model. WL‘ will for iI silmpk

problem.

Consider

a midplanc

symmetric laminate which is intinitcly long in the .‘c-and,,,-directions, but contains a disbond of length 2~ at the midplnne (Fig. prtxsurc distribution

2). The disbond surfaces arc subjected to il normal

(cquill and opposite on tho two surfilccs) which is indcpcndcnt of the

jq-coordinate so that the problem is quasi-3-D.

For simplicity. wc will consider the laminate

as an asscmblagc of two sublaminatcs (one above and the other below the disbond) which arc themselves midpl~n~ symmetric and balanced. We also assume th3t II’,, good ussilmption

= 0 which is 3

for the c;tsc of it rtpcuting layup pattern. The n~~l~~v~lnishin~displ~i~~rnents

and stress resuitnnts itrt I#:‘, rb,. tr:‘, J/+ N,. rCf,. Q,. N,, and RI. Bccitusc of symmetry conditions about the disbondod intcrfucc wc need to consick only the top sublaminate which is stress free iit the top surface. WC. therefore. hrrvc no body forces and

It is clear that one has to solve the first. third. fourth and sixth equation of eqns (I I) with second terms and body force terms (,fZg) on the teft-hand side of each equation equal to zero. The rest of the equations in eqns (I i) are identically satisfied.

Modelinglaminatedcompositestructura PI usemblagcof sublaminatu Writing

where h is the thickness of the sublaminate, we note that

u, = 0,

1x1> a.

(21)

Making use of constitutive relations (12) and eqns (19a)-( 19c) one obtains from eqns (1 I)

+(++$)KII+*(G-*;).(22a-d) One approach often used in crack problems in elasticity is to employ Fourier transforms and reduce the system of equations to a system of integral equations for determining the displacement discontinuities. We will follow the same procedure and introduce the Fourier transform of a quantity as

t(d)

eM dx’

(23a-c)

whcrei=J-Iandx’= x/h, the non-dimensionalized x-coordinate. Let a’ = a/h. the non-dimensionaiized half disbond length. With the help of known relations between transforms of derivatives of quantities and their transforms and qns (22a)-(22c) one can express 0:. 0;. and 0: in terms of 0;. Performing the required algebra and substituting in eqn (22d) the following relation is obtained :

t

’=

--

I

f,+I,

sys*+c:) (s*+c:)(s*+c:) vJ

where ?, and P, are Fourier transforms of the quantities

(244

S. N. CHAITERJEE and

v. t ;

k4MNATH

(x’) d-r’

(24b. c)

and

(25a-f) Because of eqn (21) one must have ‘I’

V,(Y)

dx’ = 0.

(26)

u’ Equations

(24) imply that T3 and C’, each satisfy fourth-order

lx’1 > a’ and lx’1 < u’, respectively. which may bc obtained of the sublaminate

separately.

imposing appropriate conditions

continuity

This yicids an altcrnativc (of displaccmcnts)

at lx’1 = u’ with due consideration

diffcrcntial

by considcriny

equations

way for solving the problem

and equilibrium

to conccntratcd

for

the two regions by

(of stress resultants)

forces at disbond tips. In

many casts. such a solution is easier to obtain, and WC will follow this procedure in another example considered later. However. to show the existence of these forces, we take the inverse transform

of eqn

(24a). Noting

that l’,(Y)

= 0; Ix’1 > (I’ because of (21) one obtains

(27a)

where/(s)

is the rational function of s’ in the right-hand

side of eqn (24a). In what follows,

we will assume that c: are real, and denote [k as their positive square roots. Writingf(s) in the form of partial

fractions and making

use of some integral formulae[l6]

eqn (27a)

reduces to

where

b, = h? =

ICf(rI-1,)+121/{lo(l, +MG-i31 -[Tf(1~-1,)+1*1/{10(1, +MG-r;:)}.

(27b-d)

Since T,(Y) is known (see eqns (24b) and (19d)) for 1.~~1 c a’. this is a Fredholm integral equation of the second kind for determining V, in that domain. We now note that V, must be a continuous function for I.Y’I < a’, but there is no other requirement except for condition (26). It may be noted that discontinuous derivatives of displacements are possible if the pressure p(x) is in the form of a concentrated force for Ix’1 < a’ as in shear deformable beam theory.

Modeling laminated composite structures as assemblage of sublaminates

A

solution of the integral equation (which is given later for uniform

indicate that V&r’) of eqn (?I).

pressure) will

usually approaches a non-zero value as 1x1’-* a’ and therefore, because

T, has a discontinuity

at x’ = a’ given by

T,(d’)--

r&7’-)

V,(a’-)

=

and F(a’+)

(Ba)

I

z

I

where F(a’-)

4-W

indicate the values of the function

left and from the right, respectively. A similar discontinuity

as u’ is approached from the

exists at I = -a’. Therefore.

because of eqns (24)

- V,(-a’+)ii(.K’+a’)]+F(x’)

(28b)

where d denotes the delta-dirac function and F(x’) is a function which is equal to p(.~‘) (the prescribed pressure) for lx’1 < a’ and for Is’1 > a’ it is the sum of two exponentially decaying functions, majority

since it

(k = 1.2) are positive and real numbers.

of laminated composites this is true. However,

It may be noted that for a

the expressions given above are

also valid for complex values of ii (note that Cl = c) provided & are chosen as their roots with positive real parts. so that i2 = r, and h2 = 6,. WC have shown the existence of concentrated forces at disbond tips for a particular problem with only two sublaminatcs.

In a gcncral quasi-3-D

problem for an infinitely

long

laminate with many sublaminatcs and self-equilibrating tractions on prescribed delamination surfaces or on top and bottom surfaces of the laminates a similar Fourier

transform

solution

is applicable, and a set of integral equations similar to eqn (27b) may bc obtained. It may also be shown that the concentrated forces in (r,, f2, IJ at the disbond tips arc proportional to the corresponding

displnccmcnt discontinuity

of proportionality

constants

gradients at the same tips, although the

may change depending upon how many sublaminatcs arc used

in the model. The general nature of the stress (force) singularity,

however, remains the

same and will also hold for laminates of finite length. Such forces will also exist at other discontinuities dilfcrcnt Fourier

like free edges. In such cases, however, it is better to obtain solutions

manner as dcscribcd earlier. since solutions transforms

in addition

in a

to those in the form of

must bc assumed to satisfy conditions at the ends. Once the concentrated

forces I” and displacement discontinuity

gradients (U”)

at the disbond tips are evaluated,

the energy rcleasc rates may be obtained by the application of Irwin’s method considering another crack tip an infinitesimal

virtual crack closure

distance ahead of the present one. It

is easy to show that for the particular singular nature of the stress field, the result for each mode is G = ;(=I/”

(29)

whcrc forces and displacements for a particular mode (I-III) also valid for cracks other than dclaminations

are considered. The result is

(e.g. transverse cracks) and is now being

widely used in conjunction with standard finite element techniques where t’ is taken as the force transfcrrcd

at the tip node.

For complctcncss, we give the solution of eqn (27b) for uniform prcssurc, and compnrc the results for a particular lnminatc with other approximate and rigorous solutions. When

PW = T,(Y)

PO

= p’x’,

Ix’1 < a’

where P ’ = poft’lD,,

(30a-d)

S. N. CHA~EE

450

and V. Rutrum

and an exact solution of eqn (27b) is

C r,+ Cz sinh co-r’ ,, cash {@a

where C, and Cz are constants chosen to make the coefficients of sinh ckx’ (k = 1,2) equal to zero in the expression obtained after substituting eqn (30d) on the left-hand side of eqn (27b) and integrating the second term. For tanh &,a 8 1 (although it is not necessary to make this assumption) the results can be expressed in the following simple form :

(31)

where

From eqns (30d) and (28b) the displacement forces at x’ = a’ are

u;

discontinuity

gradients and concentrated

= 2Y,(a’) z - 2(C, - C2 +u”/6)

D,,

I; =

VI h2UI

(4

(32)

+I*)

and the energy release rate G, is obtained from eqn (29). The same result will be obtained from the consideration of change in crack or complementary energy with respect to a. The term crack energy means the work done by the pressure on crack opening displacement and it is often used in the literature on fracture mechanics. Figure 2 shows the displacement gradients and energy release rates obtained by the foIlowing four methods for a 25.4 mm disbond in a 64 ply ((O,/ + 45?/ T 45,/O,),), ASl-350 I plate. (I) Numerical solution of integral equations to which the exact elasticity problem is reduced[2,3]. (2) Constant strain finite elements. (3) Present theory employing four sublaminates (as shown in Fig. 2) with the two sublaminates above and below the delamination containing two 0” layers only. A finite element method discussed later is used for this purpose. (4) Closed form solution for two sublaminates based on present theory discussed above.

Modeling laminated composite structures as assemblage of sublaminates Table ASI-3501-6

451

I. Properties used for calculation

Material, ply thickness = 0.1398 mm

0“ Layers: Er = 125 GPa.

“w = Y, = 0.28,

E, = E; = 10 GPa. Y”_= 0.35

+ 45’ Lugers (smeoreJ) E, = E. = 19.93 GPa. v,, = 0.718,

G,; J 5.8 GPa

E, = II.06 GPa, Y, = wVz= 0.104.

G,, = 4.52 GPa

A shear correction factor of 5/6 is used in computation of effective shear stiffnessesin the case of one sublaminate model since a repeated lay-up pattern exists(9, IO]. In case more sublaminates are used. the correction factors are omitted.

Properties used for the calculation are given in Table 1. To reduce the number of degrees of freedom for method 2 we have replaced the +45” layers by a homogeneous material with smeared properties (Table I) based on constant in-plane strain and constant transverse stress assumptions. The same properties are also used in all other methods for the purpose of comparison, although it is not necessary to do so. As can be seen from the figure, the crack opening displacement gradients as well as energy release rates from the first three methods are extremely close to one another. It should be noted, however, that the displacement gradients at the tips are unbounded for the exact elasticity solution, but they are bounded for all other methods. The displacement gradients near the tips obtained from the closed form solution are quite different from those from other methods. Similar differences will also exist in the stress field adhead of the disbond tip (see results for another problem given later). It is clear from the results of method 3, that by choosing thinner sublaminates near delaminations (or increasing the number of sublaminates) more accurate representations are obtained. This phenomena is analogous to that observed in standard finite element solutions with finer mesh size near crack tips. It should be noted, however, that the energy release rate from the closed form solution is extremely close to other solutions. This is not surprising since strength of materials type solutions are often found to yield accurate estimates of energy release rates in many problems. A shear deformation theory solution of this problem is given in Ref. [3]. It should be pointed out that for method 2, it is necessary to model the smeared 245” material and 0” layers separately, whereas for method 3, they are lumped in a sublaminate (except for the Crisublaminate near the disbonded interface). For this reason the degrees of freedom and the band width for method 2 are about three times those in method 3. It is clear, therefore, that the present method should yield substantial reduction (60%, or more) in computational times. In the next section we give the outline of finite element modeling using the present approach, and later we will compare the stress fields ahead of a crack tip in another sample problem obtained by increasing the number of sublaminates with those calculated from exact numerical solutions.

FINITE

ELEMENT

MODEL

FOR 2-D PROBLEMS

In finite element modeling, the laminate (or laminated structural component) is divided into a finite number of parts. For problems where some of the sublaminates are truncated before others (as in a bonded joint) each truncation location is made to coincide with one of the boundaries of the partition. Let the nth sublaminate be divided into M parts. The mth element shown in Fig. I has length I,, height h, and four nodes numbered l-4. To obtain the local stiffness matrix it is referred to a local coordinate system x’, z’ (the primes are omitted in the expressions which follow), and it is assumed that the three displacements ui in the element are given by (the subscripts I, 2, 3 of u as well as the subscripts M and n to I and h, respectively, are also omitted)

S. N. CHAITERJEE

and V. RAMNAIH

u'-$+u3-~4

II/=

+ 3h

u’_uL_((L(‘_u~)

ih

+

11’ = displacement

x

at node j, i = I.. . . ,4.

(33)

The nodal forces can be easily calculated from values of N. M. Q and R at .K= +I/2 and r’ at: = + It/2 (which are obtained by the use ofeqns (I 1) so as to satisfy static equilibrium). They could also be obtained by using the virtual work principle. This results in the following relationship between nodal forces and nodal displacements [ _N] = [K’] [ ~1 where

[FIT = [WI, N:, N:). (N:, . . .), (N:. . . .). (N:.N;,N:)]

The superscript denotes the location (node) and the subscript indicates the direction of the forces and displacements and [K’] is the local ( I2 x 12) symmetric stiffness matrix which is better written as an assemblage of I6 (3 x 3) submatrices [k”“] given in the Appendix. Note that

[K”‘]

= [Py.

should bc pointed out that if local coordinates of some of the elements arc inclined with rcspcct to global s.: coordinates (as in the problem of ply drops) appropriate transformation laws have to bc utilized to obtain rotated local stitl’ncss matrices bcforc the global matrix is ilSSCfllblCd. In some CilSCS. use of triangular clcmcnts arc ncccssary. local stitTncss matrices of which arc discussed in Ref. [ 151. String stilTncsscs due to large deformations which should bc added to the terms of the stitTncss matrix of rectangular elements and modilications ncccssary to consider the elastoplustic shear stress-strain relation of adhesive clcmcnts illX! ill!40 given in Kcf. [ I5]. It

TIIE

DOUBLE

CANTILEVER

IIEAM

PROBLEM

In this section. we consider a wide double cantilever beam shown in Fig. 3 subjected to two equal and opposite forces P on the two arms. We will assume that these forces are thr resultants of parabolic shear stress distributions, although any general distribution may bc considered without any difficulty. We first give a closed form solution for two sublaminatcs (one above and the other below the crack) each of which is midplane symmetric and balanced. As in the previous illustrative example, we also assume D’, , = 0 and consider the top sublaminate only utilizing the symmetry conditions about : = 0. Howcvcr. we will follow a ditferent routr by choosing solutions of the differential equations, eqns (I 1), in the cracked ( -N < s c 0) and untracked parts (X > 0) as discussed earlier. We will also assume thilt the uncrackcd length is large compared to u and it can be considcrcd to bc infinitely long. For brevity, WComit the details and give the expression of non-zero quantities as a function of .Y’= s,‘/I. /t is the thickness of the sublaminatcs. a’ = o//z and h is the width of the beam. i,,. i,. iz. 1,. I?. I, are the same as in eqns (25). Utilizing the end conditions at I = -u’ (N, = M, = R, = 0 and Q, = -P/h). and the fact that 11;-II/ ,/J/2 = 0 at x’ = 0. the solutions of the equ.ations for --a’ < X’ -X0 arc as follows : liI:/il

=

-Ps’~hL,,+PJ1’.r”/6hD,,+Plt~a’.r’~/2hD,,-C,.u’+C,~2

II/>= C,(cosh co-v’+ tanh coa’ sinh Co.\-‘)+ Cz

$, = - Pi~‘.r”/ZhD,, - Ph%i.r'/hD, u:)/h = C:-

, + C,

C,(sinh i,,s’+ tanh ino’ cash I&).4 ,4/A, Jo

Modeling laminated composite structures as assemblage of sublaminatcs

453

Shepr.QrfPrmotig?!

1 Presentsc4utloM one&Jte (cloud form) (31 One PI& (closed form It2

60

- --_

Two ptotcs (flnlte element I(41 Four phates (finite element 1(5)

I

Elostlc~ty

so . -

(1)

NumerIcal solution [I81

40

Force

crack top Dasd. qrod.

~lutlon

(;b/P

V;bD,, /

131

449

5.13

23.0

(4)

4*t

4.92f

20.9

3.72~

568f

21. I

Ouonthes

Energy release rate

As

a

G~b*&/P*h*

30

\

2o

P c I

(5) IO

t Indlcotes

values at tip node.

0

-101

-20 I

I OS

(

I

I5 x/h

IO

I

I

20

25

I 30

Fig. 3. Strcssdistribution ahcad ofcrack tip in an isotropic DCB spccimcn-from diGrent solutions.

N, -0;

Q,

=

M, = - Ph(x’+u’)/h

-P/h,

R , = L’:,C,C,,(sinh

co-r’+ tanh <,,u’ cash C,_r’)/h

N, = C,(cosh cOx’+ tanh C,,u’ sinh Co-r’)(Ad., - A ;,/A, ,) where C,. Cl. C, are arbitrary

(35)

constants.

All quantities should decay as x is increased in the solutions for x’ > 0, which can be written in the following form (B,, k = I, 2 being arbitrary constants; when Cz are complex (; are their roots with positive real parts such that Cz = [, and Bz = B,)

k-l

IL1=

24/h;

$1 =

i

& exP

(-c&x’)

k-l

u:/h = -2A,4

2

f&(1,--I/j,?)

II

A,,

exp (-LX’)

k-l

Nl =O;

Q,

= L,,l,

M, = -D,,

2 &i,2eXp(-h’)

t-1

t k-l

BA

=P

(--L-r’)

II

h

G = DII

1h3.

t B& eXp (-i&i) k-l

1

(36)

1

Note that U! -Jl&2 = 0. Denoting the value of the quantity f(Y) as .Y‘-+ 0 from the left and from the right as f(O-) and f(O’), respectively, the continuity of displacement components at and equilibrium conditions across x’ = 0 (with the consideration of a concentrated force t: in the z-direction at the crack tip and taking an infinitesimal element extending on both sides of x = 0) will be satisfied if

J/,(0-) = J/,(0’); Mt(O-)

= M,(O+);

u,(o+) = u,(O-);

Q,(O+)-Q,(O-)+r;

=0:

u,(O”) = u1(0-) R,(O+)-R,(O-)-t’lhj2

(37a-c) = 0. (38a-c)

Note that there are six unknowns (Cr. C1. CJ, B,, Bz. f?) to bc evaluated from six equ~~ti~~ns. Whitney[~7] obtained a solution to the problem (wiih minor diffcrcncos) for the cast of ir homogeneous beam without any consideration to the conccnttatcd force r;. Thr equilibrium equation involving RI, eqn (38~). was left unsatisfied, since thcrc arc not enough unknowns without r:. One should note that even if ordinary shear doformation theory is used to solve the problem the concentrated force r; must be considered, othcrwisc one will bc unable to satisfy the equilibrium equation, eqn (38b), involving Q,. Based on this fact we may conjecture that even with the use of theories of still higher order. concentrated forces will exist in such sublaminate assemblage models. Hopefully, this will bc investigated in future studies. For completeness, the constants obtained from eqns (37a). (37~) and (3&t) (38~) which will be needed later are given below. C2 may be obtained from eqn (37b) B, = -[P*(l

+I+‘)-+l;*]/I,(C,

Bz = [P*(I +T,a’)+r;*]/i:(i,

-iz) -iz)

C, = B,+Bz

(SO)

where f,, II are given in eqns (3 I) and P+ = Ph’jhn,, p* = t;h’/D,, 3 /, = tanh Coo’ IS = f,io+l,

+I2

I, = (i,+12)(~*+~:).

(40)

Other relevant quantities like stress distribution for X’ > 0 (rj) can be obtained from eqns

Modeling laminated composite structures as asscmbiagc of sublaminates

455

Tablt 2 Value of C,

Theory Elasticity[ IS] Prescat theory Shear deformation theory

2.34 2.77

1.85

(35) and (36). The gradient of the displacement di~ontinuity energy release rate are

6’: at the crack tip and the

u: = 2 * tS+(i, +I,) G, = gu:.

(41)

For ordinary shear deformation theory the concentrated force and the displacement discontinuity gradient (denoted with a subscript S) at the tip are

For values of n’ such that fi z I (see eqns (40)) the energy release rate in eqns (41) can bc exprcsscd in the form G, =

P2h2(u’+C~)2/b2D,,

(43)

und a similar expression holds for shear dcformation theory C., being a constant dcpcnding on the theory used. For an isotropic plate with Poisson’s ratio Y= 0.3 and correction factors of 516 (with l., , and L; ,) the stress itlt~nsity factor ~al~ulilt~d to give the same Gr can h written in the same form as that for elasticity theory/ IS] for cf‘ > 4

whcrc the values of C6 are given in Table 2. Although the contribution of C, to Gr (or !-I,) is small for large a’, this term is absent in Whitney’s solution[l’l]. It should be noted that the cncrgy release rate can also be calculated from

G=

P du!( -a) 7;

&

and the result is the same as in eqns (41) or eqns (42) although the algebra is quite complex to show this equivalence in the case of eqns (41). Figure 3 gives the comparison of the stress distribution ahcad of the crack tip obtained for the isotropic case calculated from five solutions. (I) Elasticity solution(l8). (2) Closed form solution (shear deformation theory), (3) Closed form solution (present formulation), one sublaminate for z > 0. (4) Finite element solution based on present formulation with two sublaminates used to model the half: > 0. (5) Same as above with four sublaminates. For methods 4 and 5 we have not used any correction factors, although 5/6 is utilized for 2 and 3. Also given in this figure are the values of concentrated forces and displacement

S. N. CHAI-I-EIUEE

156

and V.

RAMNATH

gradients at the tip obtained from methods 2-5 which show that the concentrated forces reduce and the gradients increase gradually as more sublaminates (or higher order theory) are used. This

is consistent

with the fact that in the case of the exact solution,

the dis-

placement gradients become unbounded. it should be noted, however. that the products of the force and the gradients for all cases are roughly the same, since energy release rates are usually not sensitive to the model used. It is seen, however. that the stress distribution

ahead

of the tip obtained from the present formulation

as the

approach the elasticity solution

number of sublaminates is increased.

CONCLUDING

The

mixed variational

assemblage of sublaminates computational

principle

REMARKS

and the approximate

formulation

are shown to provide a powerful.

in terms of an

yet relatively inexpensive

algorithm for approximate stress analysis in composite laminates in presence

of delaminations.

Non-zero

values of displacement gradients and concentrated forces at

disbond tips are shown to be the exact mathematical forms of singular fields in the solution of mixed boundary value problems using the sublaminate assemblage model. Obviously, they differ from the singular fields in exact elasticity solutions,

but better representations

of such fields are obtained by increasing the number of sublaminates.

It is also shown that

energy release rates can be obtained from the local field at crack tips, which is very useful in complicated problems. The finite element model described is for quasi-3-D problems and can be cxtcnded to 3-D problems substantial

with a littlc cKort and such extensions will result in

reduction in computer time rcquircd for 3-D analysts as compared to standard

linite clcmcnt methods whcrc layers of ditTcrcnt orientations Results

for other delamination

have to bc modeled separately.

problems computed with this model have been compared

cxtcnsivcly in Ref. [IS] with those from othsr linitc clcmcnt methods modcling each layer. as well as from numerical solutions

based on rigorous elasticity theory. Applications

problem of stress analysis near a ply-termination joints have also been dcmonstratcd rcquircd in transition

to the

and fracture mechanics analysis of bonded

in Ref. [ 151.The stilrncss matrix of triangular clcmcnts

regions for ply drop problems can bc found in Ref. [I 51. It should bc

pointed out that in 3-D or 2-D problems of dcluminations

bctwccn two dissimilar

isotropic

or anisotropic layers, oscillatory type stress singularitics are often observed which sometimes result in part closure of disbonds as discussed in Rcfs [I, 3, 5. 81. In the prcscnt formulation, such oscillatory

singularities

expcctcd if the intrructivc

arc absent, but part closures arc still possible. and should bc

concentrated force at a tip is comprcssivc in nature. The extent

of such part closures may be detcrmincd if an appropriate purpose. Further

algorithm

is devised for this

studies in this direction will be helpful for obtaining a better understanding

of the effect of interlaminar

stress fields on fatigue and fracture behavior on laminated

composite components. .~c.k,to,~~~,c~~c,,,r~*,tr.v-Thcauthors sxprcss lhcir sinccrc ;~pprcciaGon lo lhc Naval sponsoring this work under Conlract No. N-62269-82-C-0705.

Air Dcvclopmcnl

Ccnhx

for

REFERENCES I. f:. Erdogan and G. D. Gupta. Laycrcd composircs wilh an inlcrfxc flaw. /nr. J. Solic.!~ Sfrucfures 7. 1080 (1971). 2. S. N. Chat~crjcc. On inlcriaminar dcfecls in taminlwd composilrs. In ~\l&w~ Deeelopmenf.s in Composiv bfureridv uml S~rucrures (Edilcd by J. R. Vinson). p. I. ASME (1979). 3. S. N. Chaltwrjce. Three and IWO dimcnsionat sircssliclds near dclaminaGons in laminated composire plntcs. brr. J. So/icLySwucrurer 23. I535 ( l9U7). 4. 0. J. Wilkins, J. R. Eiscnman, R. A. Cnmin. W. S. Mqolis and R. A. Rcnson. Chnraclcrizing dclaminxtion growth in graphite epoxy. In Dmtuqr in Composire ~\fu/erirrfs : f3~sic,~fcrfr[r,fisnt.~. A ccumufurion. Tokercmw cd Churacreri~utbm. ASTM STP 775. p. 162 ( 1982). 5. S. S. Wang and F. G. Yuan. A hyhrid finilc element approach to composite k~minale ctaslicity problems with singularilics. J. rlppl. .tfedr. SO. X35 ( t 9R3). 6. S. N. Chattcrjce and R. B. Pipcs. Composite dcfcct signilicancc. Contract Report. NADC ROO4R-60. August (1981). 7. S. N. Chatbxjcc and R. B. Pipcs. Study of graphilc epoxy composites for matcrial flaw criticali(y. Contract Reporl. NADC 7HZJI-60. November (1980).

Modeling laminated composite structures as assemblage of sublaminates

-07

8. S. N. Chatterjee. R. B. Pipes and W. A. Dick. Mixed mode delamination fracture in laminated composites. CumposiresSri. Technol. Ls, 49 ( 1986). 9. P. C. Chou and J. Carleone. Transverse shear in laminated plate theories. A/AA J. Il. 1333 (1973). IO. S. N. Chatterjee and S. V. Kulkami. Shear correction factors for laminated plates. A/AA J. 17.498 (1979). II. N. J. Papano. Stress fields in composite laminates. fnr. J. Soli& Srructures 14, 385 (1978). 12. N. J. Pa8ano. Free edge stress fields in composite laminates. Inc. 1. Soliu!v Srrucrures 144OI (1978). 13. K. Hellan. Debond dynamics of an elastic strip, I. Timoshenko beam properties and steady motion. fnr. J. Fracture I4 91 (1978). I-t. R. M. Christensen. Mechanics uf Cumposire Muteri&. Wiley. New York (1979). IS. S. N. Chatterjee. V. Ramnath. R. B. Pipes and W. A. Dick. Composite defmt significance. Contract Report. NADC Contract No. N-62269-82-C-0705. February (1985). 16. I. S. Gadshteyn and I. M. Ryzhtk. Tuhle o//nregrols, Series, Products. Academic Press. New York (1965). 17. J. M. Whitney. Stress analysis of the double cantilever beam specimen. Composites Sci. Technol. 23. 201 (1985). 18. D. J. Chang. R. Muki and R. A. Westmann. Double cantilever beam models in adhesive mechanics. Int. J. S&b Srnrcfures 12, I3 ( 1976).

APPENDlX:

ELEMENTS

OF SUBMATRICES

UII +2Pll

[/f-“1 =

b,, +?PI*

-ulr-o?l

J,,+?p:,

-4,-d,,

Sym.

1

El, +2p.!r

-u12

-PII

-hll-PlI

U,,+~~l,

-b,,

-p,:

-d,,-p,,

d,,+dlz

(I,, +or,

Ih”‘l =

OF LOCAL

cl,, fd?,

-e,:-P,,

-~~,:-P,,

-h

-b,,-p,z

-c/,1-p:.!

-d,,

-o11-o11

-d,:-cl,,

-e,:-PI,

L

21

-PI?

_LI,# -u:, -d:J

-bx+p,,

u:1 -UC

-d,:+p:,

d,,-d,l

d12 -d?,

-e,z

+p,,

1 1

u,,+zp,, [KJ’J =

[

Sym. b ,:-b,t

o,,--a,1

d ,:-2P:1

d,2-d,,

d :, -d,,

e,,-2p,,

h,,+2p,,

-(I:,-a,:’

d:,+2p:z

-d:.,-d,, e:: CZp,,

1

STIFFNESS

MATRIX

S. N. CHAITEIUEE

P.8

=

L,p1J6h2:

a. p = I. 2.

and V.

RAMNATH