Efficient 3D analysis of laminate structures using ABD-equivalent material models

15 Pages·2015·6.16 MB·English

by  Goldy Kumar

FiniteElementsinAnalysisandDesign106(2015)41–55
ContentslistsavailableatScienceDirect
Finite Elements in Analysis and Design
journal homepage: www.elsevier.com/locate/finel
Efficient 3D analysis of laminate structures using ABD-equivalent
material models
Goldy Kumarn, Vadim Shapiro
SpatialAutomationLab,MechanicalEngineering,UniversityofWisconsin,Madison,USA
a r t i c l e i n f o a b s t r a c t
Articlehistory: Laminatecompositesarewidelyusedinautomotive,aerospace,andincreasinglyinconsumerindustries,
Received28August2014 duetotheirreducedweightandsuperiorstructuralproperties.However,structuralanalysisofcomplex
Receivedinrevisedform laminatestructuresremainschallenging.2Dfiniteelementmethodsbasedonplate/shelltheoriesmay
17July2015 be accurate and efficient, but they generally do not apply to the whole structure and require
Accepted27July2015 identificationandpreprocessingoftheregionswheretheunderlyingassumptionshold.Fullyautomated
structural analysis using solid 3D elements with sufficiently high order basis functions is possible in
Keywords: principle,butisrarelypracticedduetothesignificantincreaseinthecostofcomputationalintegration
Composite overalargenumberoflaminateplies.
Laminate
We propose a procedure to replace the original laminate by much simpler new virtual material
Finiteelementmethod
models.Thesevirtualmaterialmodels,undertheusualassumptionsmadeinlaminationtheory,havethe
Materialmodel
sameconstitutiverelationshipasthecorresponding2Dplatemodeloftheoriginallaminate,butrequire
onlyasmallfractionofcomputationalintegrationcostsin3DFEA.Wedescribeimplementationof3D
FEA using these material models in a meshfree system using second order B-spline basis functions.
Finally,wedemonstratetheirvaliditybyshowingagreementbetweencomputedandknownresultsfor
standardproblems.
&2015PublishedbyElsevierB.V.
1. Introduction commonpracticeistoassumethatthelayersarepermanentlyfused
together and ignore any fluctuation in stress–strain fields at the
1.1. Motivation interfaces of layers [3–6]. These assumptions allow to approximate
laminate'sglobalbehaviorasthatofaplateorashell.Interlaminar
Laminatecompositesarenowusedwidelyinvarietyofindustries, stresses and strains may be significant in boundary regions and
including aerospace, automobile, medical and sports [1–3]. Lami- regions of discontinuities [7] where full three-dimensional and/or
natesarelightweightandstiffwithcustomizablematerialproperties, layered methods should be used, but the plate/shell assumptions
resulting in structures superior to those made of homogeneous givesufficientlyaccuratestressandstrainestimatesforregionsaway
materials[2–4].Highstiffness-to-weightratioisachievedusingfiber from those regions [8]. In this paper, we will show that the same
reinforced plies. These plies, when fused together under high plate/shell assumptions, when applicable, may be used within a
temperatureandpressure,formcomplexmonolithiclaminateparts. general 3D finite element analysis to dramatically speed up the
Thefiberreinforcements,laidusingtechniquesrangingfrommanual analysisprocedure.
to fully automatic, are generally parallel and unidirectional and, Structural analysis of laminates can be carried out using
therefore, result in plies which are anisotropic in nature. Material differentfiniteelementmethods,andsomeofthemareillustrated
properties are customized by varying fiber angle within each ply, foratypicallaminatepartinFig.1.Duringfiniteelementanalysis
controlling the number of plies, and adding additional materials (FEA), stiffness matrix K for each element must be computed,
e
between plies such as cores and fillers. The presence of numerous whichingeneralformisgivenas[9]
Z
plies,however,leadstocomplexgeometryandmaterialdistribution
inlaminatestructures,and,therefore,structuralanalysisoflaminates Ke¼ BT(cid:2)Q(cid:2)BdΩ; ð1Þ
Ω
bytreatingeachplylayerindividuallyisprohibitivelyexpensive.The e
where B is the strain–displacement matrix, Q is the material
Ω
constitutiverelationmatrix,and istheelement'sdomainover
e
whichintegrationisdone.Sincetherearenumerousplies,mesh-
nCorrespondingauthor.Tel.:þ16083349646.
ing each ply independently requires a large number of elements
E-mailaddresses:[email protected](G.Kumar),
[email protected](V.Shapiro). andis,therefore,prohibitivelyexpensive.Amuchsmallernumber
http://dx.doi.org/10.1016/j.finel.2015.07.009
0168-874X/&2015PublishedbyElsevierB.V.
42 G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55
Fig.1. Figure shows a structure made of 3 laminates analyzed using different finite element methods. (a) Plate/shell element. (b) Element of 3D conforming mesh.
(c)Elementofa3Dnon-conformingmesh.(d)Elementofanon-conformingmeshwithcurvedlaminateinside.
of elements is needed if layered element method [10–12] is used composite laminate structures that is as general as 3D FEA and
(Fig. 1b and c), wherein multiple plies can cross through an asefficientas2DFEAwhendimensionalreductionmakessense.
element. Such finite element methods (also called 3D FEA since Specifically,weproposeamethodtoreducetheexcessivecostof
integration domain Ω is 3D) may be further classified into integrationforlayeredelementsbytakingadvantageofplate/shell
e
conforming and non-conforming (or meshfree), depending on natureoflaminates,wheneversuchassumptionsarereasonable.To
whether the finite element mesh conforms to the geometry of this end, we have devised a procedure to obtain material models
thelaminate.Threedimensionalelementsmayexhibitlockingand whicharesimplerbutareequivalenttotheoriginallaminate,under
ill-conditioning of stiffness matrix when used for laminates that theassumptionmadeinlaminationtheories.Werefertothesenew
are thin in the plylayup direction [9,13]. These problems can be materialmodelsasABD-equivalent material models, asthey result
alleviated or eliminated by using higher order hierarchical1 basis in the same ABD matrices as the original laminate and, therefore,
functions[14,15]. can replace the original laminate during integration if plate/shell
All3DFEAmethodsarecomputationallyexpensive,astheabove assumptionsapply.We demonstrate the effectiveness oftwosuch
integrationhastobecarriedoverlargenumberofplies(tensoreven material models—a 3-ply and a graded material model—in a non-
hundreds). Integration is performed using quadrature rules that conforming FEA system using layered solid elements. We validate
depend on the geometry of the element as well as the degree of thetwoABD-equivalentmaterialmodelsbyusingthemtoanalyze
theintegrand,andamountstosamplingtheintegrandatanumberof several benchmark problems, and compare obtained results from
quadraturepoints[9].Togetanideaofthehighcostofintegration knownresults.Thefullyimplementednon-conformingFEAsystem
forlaminates,weconsiderthelayeredelementusedinreference[10] uses layered solid elements with second-degree B-spline basis
to analyze a laminate made of 100 plies. The element used is an functionsthatarehierarchicalinnature.
eight-nodebrickelementwithtri-linearbasisfunctions,which,fora Abriefoutlineofthepaperisasfollows.InSection2,wesurvey
homogeneousmaterial,isfullyintegratedusing2integrationpoints therelatedwork.Section3developstheconceptofABD-equivalent
ineachdirection,or8integrationpointsintotal[10].However,ina material models and proposes two specific examples of such
laminate,8integrationpointsareneededforeachply,whichresults models.Implementationoftheproposedapproachincombination
ina100foldincreaseforour100-plylaminate.Sinceintegrationcost with a non-conforming finite element method is described in
represents a significant portion of the overall solution procedure, Section 4. Its effectiveness is demonstrated using a number of
analysis of composite laminates using layered elements is an benchmark problems in Section 5, followed by conclusions and
expensiveproposition. futureworkinSection6.
Plate or shell assumptions reduce the computation cost and
increase accuracy of FEA for laminates owing to their thinwalled
nature.Theseassumptionsmayleadtodifferentlaminationtheories, 2. Relatedwork
where the material matrices Q of allthe plies are replacedbythe
so-called ABD matrices [4,16]. The structure and the integration The finite element methods for simulating global behavior of
domainΩ effectivelyreducestoa surface (Fig.1a),whichiswhy laminatestructures[16–18]maybebroadlyclassifiedasa2Dora
e
thismethodisalsocalled2DFEA.However,2DFEAisnotvalidin 3D FEA. For the purposes of this paper, we only consider those
regions near boundaries and discontinuities (Fig. 1), which have methodswhichignoreinter-plyphenomena,butwenotethatthe
significant 3D stresses and, therefore, plate/shell assumptions are inter-layer stress–strain can be partially predicted from global
invalid.Inthissense,2DFEAmethodsarenotgeneral,becausesuch deformation[4,19].
regionsarecommoninlaminatestructures.
2.1. Two-dimensionalFEA
1.2. Contributionsandoutline Laminates usually behave as plates or shells, and are analyzed
using 2D FEA. Depending on the strain field assumed in the
Basedontheabovediscussion,thechoicebetween2Dand3D laminate's thickness direction, different lamination theories exist
FEAamountstoatrade-offbetweengeneralityandcomputational [5,17,16], and can be classified as one of the following: Classical
efficiency. We seek to develop an approach to analysis of Lamination Plate Theory (CLPT), First Order Shear Deformation
Theory(FSDT),orHigherOrderShearDeformationTheories(HSDTs).
CLPT assumes that laminates undergo only stretching and pure
1A basis function is called hierarchical when a higher order basis function
bending (Fig. 2C): in-plane strains vary linearly in the thickness
containsallthelowerorderbasisfunctions;forexample,B-splinesarehierarchical
basisfunctions. direction, and out-of-plane strains are absent. On the other hand,
G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55 43
Fig.2. (A)Plywithfiberangle1201.(B)xzcrosssectionofalaminatewithpliesatangles45/(cid:4)60/120/0.(C)Linearvariationofstrainϵxalongzwhilelaminateisin
stretchingandpurebending.
Althoughtheseelementsdonotrequiremid-surfacesexplicitly,they
mustconformtothegeometryofthelaminate,withtheirzdirection
aligned to the transversal direction of structure's offset thickness,
because the in-plane and out-of-plane behaviors are assumed or
enhanceddifferently.Werefertosuchelementsasconforminglayered
elements(Fig.1b).Theseelements,justlikelayeredelements,arestill
expensiveforFEAoflaminates.
Fig.3. (A)Alapjointbendingunderin-planeloading,whichleadstohighstress
Aligning element's z direction to the laminate's transverse
concentrationnearthejoint[20,21].(B)Whenanalyzedasa2Dstructure,bending
thicknessdirectionalsosimplifiesvolumeintegration.Ifnotaligned,
inlapjointisnotcapturedatall.
pliescanintersectelementsatarbitraryanglesandrequirecompu-
tation of intersections between individual plies and elements,
FSTDandHSDT,asthenamessuggest,assumenon-zeroout-of-plane which is both non-trivial and computationally expensive. For this
shearstrains;FSDTassumesconstantvariationwhileHSDTsassume reason,3DFEAoflaminatesusingnon-conformingmeshbecomes
linear or even higher order [19] variation for out-of-plane shear lessattractive,despiteitsadvantagesoverusingconformingmesh
strains.In-planestrainsareassumedtobeidenticaltoCLPT,andout- [26]. For completeness, we mention that several other methods
of-plane normal strain is absent. Since our goal in this paper is to havebeenproposed[5,16,17,27],buttheyeitherlackgenerality,or
showthatABD-Equivalentmaterialmodelscanbeusedtosimplify sufferfromoneormoreofthechallengesmentionedabove.
theactualmaterialand,therefore,reducethecostof3DFEA,weonly
focus on CLPT and FSDT and not higher order theories. Fig. 2
illustrates a laminate and the variation of strain field in the x 3. Materialmodelsforcompositelaminates
directionalongthethicknessobtainedusingCLPT.
Asdiscussedabove,although2DFEAisefficient,itisnotgeneral. Inthissection,webrieflysummarizetheclassicaltheoriesused
Inaddition,2DFEAsuffersfromadditionaldrawbacksthatlimitits toestablishtheconstitutiverelationshipsinalaminate,andderive
applicability in complex laminate structures. First, it requires the ABD-equivalent material models which remain valid under iden-
structuretobedimensionallyreducedtoasurface,whichcouldbea ticalassumptions.
complex task in itself [22]. Even if the structure is successfully
reduced,modelinganassemblyofmultipleplatesandshellscanbe 3.1. Constitutiverelationsfororthotropicplies
problematic [6]. In addition, due to dimensional reduction, 2D FEA
cansometimescompletelymissa3Dphenomenon.Forexample,in Plies made of uni-directional fibers are generally modelled as
thelapjointproblemshowninFig.3,2DFEAmissesthemoments transversely isotropic materials, but we are considering a more
due to eccentric forces when the lap joint is reduced to a surface. general orthotropic model for the plies [4]. To characterize a
Finally,differenttheoriesforthinandthickplatesandshellsfurther material in linear elasticity, stiffness matrix C is used, which
complicate2DFEA[23]. requires 9 independent elastic constants in case of orthotropic
materials.Theconstitutiverelationbetweenstressandstraintakes
2.2. Three-dimensionalFEA ageneralformgivenby[4]
σ ¼C (cid:2)ϵ; i;j¼1;2;‥;6; ð2Þ
In theory, three-dimensional FEA using layered elements will i ij j
σ ϵ
accurately simulate deformation in laminate structures. In practice, where, as usual, is stress and is strain. The equation is in
however,usingsolidelementsisexpensive,andtheyalsomightgive contractednotationwherei;j¼1;2,and3arex,y,andzcoordinate
inaccurate results due to locking. As a compromise, several hybrid axes, while i;j¼4;5, and 6 are yz, zx, and xy planes respectively.
methodsthatincorporatetwo-dimensionalplateandshellbehaviors Unlikeisotropicmaterials,Cisdirectiondependentfororthotropic
in3DFEAhavebeenproposed.Forexample,solid-shellelementsare3D materials and needs to be transformed from its principal material
elements that use Assumed Natural Strain method to deform like directionstoelement'scoordinatedirections.FullmatrixformofEq.
platesandshells[6,24].Theirthreedimensionalnatureiswellsuited (2),includingthetransformationofCfromitsprincipaltoageneral
for interfacing with other solid elements in assemblies. These ele- coordinatesystem,isgiveninAppendixA.1.
ments, however, still require mid-surface extraction and alsocannot The plane–stress constitutive relationship for dimensionally
simulatebehaviorsotherthanplateandshellbehaviors. reduced laminates is characterized by a 3(cid:3)3 stiffness matrix Q
Continuum solid-shell elements, unlike solid-shell elements, are [4]andisgivenas
standard displacement based elements, but use advanced finite σ ¼Q (cid:2)ϵ; i;j¼1;2;3; ð3Þ
elementtechniqueslikeAssumedStrainMethodandEnhancedStrain i ij j
Methodtoimprovetheirperformanceforthinstructures[10,11].The where i;j¼1 and 2 stand for x and y axes respectively, while
higherthenumberofassumedandenhancedparameters,thebetteris 3 stands for the xy plane. Full matrix form of the equation is in
the element's performance, but at the expense of generality [10,25]. AppendixA.2.
44 G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55
ϵ ϵ
In thick plates, out-of-plane shear strains and are ABDmatrixmodelforlaminates:
yz xz
significant,andtherefore,out-of-planeshearstiffness,inaddition (cid:2)N(cid:3) (cid:2)A B(cid:3) (cid:2)ϵo(cid:3)
to in-plane stiffness, are needed to characterize a ply. In an ¼ (cid:2) κ : ð8Þ
σ σ M B D
arbitrarycoordinatesystem,shearstresses and arerelated
yz xz
toshearstrainsϵ andϵ as[5]: The individual coefficients of A, B and D matrices for indices
yz xz
"σ # " # "ϵ # i¼1,2,and3aregivenas
Q Q
σ4 ¼ 44 45 (cid:2) ϵ4 ; ð4Þ Z t=2 Z t=2 Z t=2
5 Q45 Q55 5 A ¼ Q dz; B ¼ Q zdz; D ¼ Q z2dz: ð9Þ
ij ij ij ij ij ij
(cid:4)t=2 (cid:4)t=2 (cid:4)t=2
whereindices4and5standforplanesyzandxzrespectively.
ThedetailsofthederivationcanbefoundinAppendixB.Intuitively,
Let us now assume that, in a laminate, z direction is the
A and D are extensional and bending components of the stiffness
thickness direction for both the laminate and its plies, which is
respectively, while B couples stiffness between bending and
alsothethirdprincipaldirectionof theorthotropicplymaterials.
stretching that occurs in a laminate if its material properties are
Recall that plate theory assumes that the thickness of a plate in
asymmetrical about its mid-plane. If B is a non-zero matrix, a
stretchingandpurebendingremainsconstant,orinotherwords,
ν ν normalpullinxorydirectioncanleadtobendingandviceversa.
Poisson's ratios and are zero[6]. This assumption reduces
thegeneralstressx–zstrainryezlationinEq.(2)to The out-of-planΓe shearΓstresses can also be reduced to mid-
2 3 2 3 2 3 planeshearforces 4and 5usingaprocedureverysimilartothe
σ ϵ
66σx77 66QQ11 QQ12 QQ13 00 00 00 77 66ϵx 77 Eoqnse.u(7se)dantodr(e8d),ulceeadthineginto-planestressestomid-planeforcesNiin
6666666τττxyyyz7777777¼6666666Q01123 Q02223 Q01333 Q044 Q045 00 7777777(cid:2)6666666γγγxyyyz7777777; ð5Þ "ΓΓ45#¼K(cid:2)"AA4445 AA4555#(cid:2)"ϵϵ45#: ð10Þ
4 xz5 4 0 0 0 Q45 Q55 0 5 4 xz5 ϵ ϵ
σ ϵ Here,strains and areassumedconstantinthezdirection,and
z 0 0 0 0 0 E3 z anydeviation4fromth5eactualfieldiscorrectedusingacorrection
wheretheindicesCijarereorderedtomatchthoseofQij.Also,note factorK[5].TheextensionalshearstiffnesscoefficientsA44,A45,and
that C33 reduces to E3, Young's Modulus in the third principal A55 aredefinedasAijinEq.(9).
materialdirection.
Setting transverse Poisson's ratios to zero is invalid when 3D 3.3. ABD-equivalentmaterialmodelsoflaminate
stresses,includinginterlaminarstresses,inlaminatestructuresare
being computed. However, the assumption is valid for ABD- IfweacceptthatABDmatricesareanaccurateapproximation
equivalentmaterialmodelsbecausetheyarebasedonplate/shell ofalaminate'sbehavior,itstandstoreasonthatanytwomaterial
assumptions for laminates which ignore interlaminar effects. As models that result in identical ABD matrices should be deemed
mentioned in Section 1.1, plate/shell assumptions are reasonable equivalent. In fact, multiple material models with identical ABD
assumptionsforregionsawayfromboundariesanddiscontinuities matrices exist, forming an equivalence class of material models.
inlaminates,whereinterlaminarstressesarenotsignificant. Fromthisclass,wecanseeksomeofthesimplestmodels,which
canbeusedtoreplacetheoriginallaminateinthe3Dintegration
procedureinFEA(Fig.4).
3.2. Constitutivemodelforlaminates,orlaminationtheories
Any new ABD-equivalent model must satisfy the following
equivalencerelationship:
The classical lamination plate theory (CLPT) assumes that
Z Z
slfaotrmraeii,nnsatϵrteoiasianncϵadinactuoarnvnlyaytuuprnoeidneκtriiganotttshhtereelatmcmhidiinn-pagtlaeanncedanapnbudereirseblgaeitnvededninliagns.eTa[4hrl]eyrteo- Aoij¼Z(cid:4)tt=t==222Qoijdz¼ Z(cid:4)t=tt=2=22Qnijdz;
Bo¼ Qozdz¼ Qnzdz;
ϵi¼ϵoiþz(cid:2)κi; i¼1;2;3; ð6Þ ij Z(cid:4)t=2 ij Z(cid:4)t=2 ij
t=2 t=2
where z is the distance of the point from the mid-plane. Fig. 2C Do¼ Qoz2dz¼ Qnz2dz; ð11Þ
showstheplotofϵ1atatypicalcrosssectionofalaminate. ij (cid:4)t=2 ij (cid:4)t=2 ij
Sincetheamountsofstretchingandbendinginaplatedepend wheretheoriginallaminateQo definesthematricesAo,BoandDo,
onlyonthenetforcesandmoments,thestressesinthelaminate's andthenewmaterialmodelQijn resultsinthesameABijDmij atriceisj.
cross-section can be reduced to just mid-plane forces N and ij
i Since the above integral equations are a system of three
σ
moments M. This is done by integrating stresses over the
i i equations for each entry of ABD matrices, theycan be completely
thicknessofthelaminatefrom (cid:4)t=2tot=2andisgivenas determined by a material model Qn that varies in the thickness
Z t=2 Z t=2 direction,andthevariationisfullyspijecifiedby3ormoreindepen-
Ni¼ σidz and Mi¼ σizdz: ð7Þ dent coefficients. There are infinitely many such models, and any
(cid:4)t=2 (cid:4)t=2
two of them are interchangeable if the assumptions made in
Combiningthisequationwithplane–stressconstitutiverelation- lamination theory hold. In other words, for the purpose of finite
shipfromEq.(3)andstrainfieldfromEq.(6)leadstotheso-called elementanalysis,anyarbitrarilycomplexlaminatewithnumerous
Fig.4. OriginallaminatewithmaterialpropertiesQoijcanbereplacedbyanewmaterialmodelQnijthatisABDequivalenttoQoij.
G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55 45
plies can be replaced by a much simpler material model yielding given laminate. Again, there is a unique quadratically varying
identicalresults.Furthermore,sincethenewmodelisvirtual,itis gradedmaterialforthegivenABDmatrices.
not subject to manufacturingconstraintsand does not necessarily
havetobeply-based.
3.3.3. TransversematerialpropertiesoftheABD-equivalentmaterial
Fordemonstration,wechoosetwosuchmaterialmodels:a3-
models
ply laminate and a quadratically graded material model. In the
In addition to in-plane, we also require out-of-plane material
nexttwosubsections,thein-planematerialpropertiesofthetwo
propertiestocompletelycharacterizetheABD-equivalentmodels.
ABD-equivalent models respectively are derived, followed in
These out-of-plane, or transverse, material properties can be
Section3.3.3bythederivationofout-of-planematerialproperties
derived using approaches similar to above, and are common for
thatarecommonforboththemodels.
boththeABD-equivalentmodels.
From Eqs. (9) and (10), the transverse shear properties of a
3.3.1. Three-plymaterialmodel laminatearegivenbytheextensionalstiffnessmatrixA fori;j¼4
ij
In the equivalence class of material models with given ABD and 5. Since we are only interested in the equivalent material
matrices,a3-plymaterialmodelisthesimplestply-basedmodel. behavior, transverse shear stiffness Qnk can be assumed constant
ij
Fig.5Ashowsa3-plymodelwithQnkasthematerialpropertiesof along the laminate's thickness. When substituted in Eq. (11), Qn
ij ij
thekthply,andforsimplicity,eachplyisassumedtobeofequal are obtained as the average values of Ao over the laminate's
ij
thickness. With these assumptions, Eq. (11) can be solved for thicknesst:
materialpropertiesQnk,andintermsofABDmatricesofthegiven
ij Ao
laminate,aregivenas Qn¼ ij fori;j¼4;5: ð15Þ
ij t
36Do(cid:4)18tBo(cid:4)t2Ao 13t2Ao(cid:4)36Do
Qn1¼ ij ij ij; Qn2¼ ij ij; Transverse normal stiffness of laminates, which is not at all
ij 8t3 ij 4t3
requiredfor2DFEA,isneededfortheABD-equivalentmodels.From
Qn3¼36Doijþ18tBoij(cid:4)t2Aoij; ð12Þ Eq. (5), transverse normal stiffness for plate structures reduces to
ij 8t3 Young's Modulus E . It is well known that for thin layered
3
wheretisthelaminate'stotalthicknessandindicesi;j¼1;2,and structures, the resultant out-of-plane Young's Modulus E3o can be
approximatedastheharmonicaverageoftheYoung'sModulusEokof
3.Clearly,therealwaysexistsaunique3-plylaminatethatisABD- 3
theindividualpliesoftheoriginallaminate[28,4].Again,sincewe
equivalent to the original laminate. Note that only the top and
are only interested in the equivalent material behavior of the
bottom plies depend on the B matrix and capture any material
proposedABD-equivalentmodels,theirYoung'sModulusEn canbe
asymmetry about the mid-plane. If the original laminate is 3
symmetrical, theBmatrixiszero,andthetwopliesQn1 andQn3 assumedconstantthroughoutthelaminatethickness.Thisassump-
ij ij tionmakesEn identicaltoEo:
areidentical. 3 3
1 1 kX¼nhk
¼ ¼ ; ð16Þ
3.3.2. Quadraticallygradedanisotropicmaterial En3 Eo3 k¼1Eo3k
Insteadofaply-basedmodel,wecanalsoreplacetheoriginal
laminate by a continuously varying, or graded, material. Since wherehkisthethicknessofthekthplyandnisthetotalnumberof
there are 3 equations to be satisfied, a quadratic variation with pliesintheoriginallaminate.
3 independent coefficients Λk with k¼1;2;3 is sufficient. A To summarize, we can efficiently construct various virtual
ij
quadratically varying material model is shown in Fig. 5B and is material models that are ABD-equivalent to the original laminate.
givenas If the usual lamination theory assumptions hold for the original
laminate, these virtual models will result in identical stiffness
Qnij¼z2Λ2ijþzΛ1ijþΛ0ij: ð13Þ matrices duringany finiteelementanalysisprocedure.Anyoneof
Λk canbefoundfromtheequivalencerelationsinEq.(11),andare thesemodelscanbeusedduringanalysis;however,somemodels
ij couldbeeasiertoimplementthanothersinaparticularsystem.For
givenas
example,the3-plylaminatemodelisstraightforwardtoimplement
15ð12Do(cid:4)tAoÞ 12Bo 3ð3tAo(cid:4)20DoÞ insystems which alreadysupport representationoflaminates. On
Λ0¼ ij ij ; Λ1¼ ij; Λ2¼ ij ij : ð14Þ
ij t5 ij t3 ij 4t3 theotherhand,thegradedmaterialmodelcanbeusedtoanalyze
Interestingly,coefficientsΛk takedifferentrolesinthematerial laminates in systems that are meant for graded materials, but do
model: together, the quadratiijc coefficient Λ2 and the constant notsupportlaminates.Inthenextsection,wediscussimplementa-
coefficient Λ0 capture the bending and in-pliajne stiffness, while tionofthesetwoABD-equivalentmaterialmodelsanddemonstrate
the linear coijefficient Λ1 captures the coupling stiffness of the theirvalidity.
ij
4. Implementation
Laminates,asproposedbythecurrentandemergingstandards,
are commonly represented as a base surface and an associated
layup table with an entry for each ply [29]. Base surfaces are
generallythetoolingsurfacesonwhichpliesarelaid,andthetable
specifiestheorder,materials,andfiberdirectionsoftheindividual
plies.Insystemssupportinglaminatesbasedontheabovestandard,
implementingthe3-plymodelisstraightforward:theoriginallayup
table with any numberof entries is replaced bya new table with
only 3 entries. However, for implementing the graded material
Fig.5. (A)A3-plyABD-equivalentlaminatemodelwithstiffnessmatrixQinjkforthe model, instead of a table for plies, coefficients Λ defining the
kth ply. (B) A graded ABD-equivalent material with its properties varying ij
quadratically. quadraticallygradedmaterialinEq.(13)arerequired.
46 G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55
We implemented the ABD-equivalent models in a meshfree
systemcalledScanandSolve(SnS)[26].InSnS,displacementsand
stresses are approximated using multi-variate B-spline functions
thatareconstructedoverauniformCartesiangrid.Thischoiceof
thebasisfunctionaddressestheconcernsofshearlockingaswell
asnumericalill-conditioningofstiffnessmatrixforawiderangeof
laminate thicknesses. More details about Scan and Solve can be
foundintheReference[26].
Duringfiniteelementanalysis,theABD-equivalentmaterialmodels
comeintopicturewhilecomputingelementstiffnessmatricesgivenin
Eq.(1).Ifthemeshisconformingandtheelement'sz-axisisalignedto
the laminate's thickness direction (Fig.1b), computing volume inte-
gration for the 3-ply model is straightforward. This is because, in a Fig.7. (A)Anxycross-sectionoftheCartesiangrid,andanarbitraryfiberinthat
conformingmesh,anindividualply'sexactlocationcanbecompletely cross-section.(B)Figurezoomsinoneofthegridelementsandshowsatriangle
determined by its position in the z direction. However volume thatisbeingintegrated.Fromthefiberorientation,thematerialprincipaldirections
integration is more involved for a non-conforming mesh. Plies can 1and2arefound,whicharenotalignedtoelementdirectionsxandyingeneral.
intersectagridelementatarbitraryangles(Fig.1candd),andalsoin
structures made of multiple laminates, more than one laminate can AppendixA.1andinmatrixformisgivenas
intersectanelement.Asaresult,computingintersectionofeachply Qn0¼GT(cid:2)Qn(cid:2)G; ð17Þ
with an element can be both complicated as well as expensive.
Therefore, for an ease of implementation, we approximate volume whereGisthetransformationmatrix.
integrationbyintegrationoversurfaces:eachlaminateisreplacedbya Finally, we assume plate behavior everywhere in our simple
setof surfaces paralleltothe toolingsurface (Fig.6).These surfaces, benchmark laminate problems, with no distinction between plate
whichwecallintegrationsurfaces,canbeeasilygeneratedastooling and non-plate regions. This assumption is justified because we
surface's offsets, a standard geometric operation in a CAD software. compareourresultswiththoseobtainedusingdimensionallyreduced
Thelocationoftheseintegrationsurfacesinthelaminate'sthickness models of laminates. We skip rest of the implementation details as
directioncanbeobtainedusingoneofthequadraturerules,andinthe theyarenotdirectlyrelevanttothecurrentwork.
currentimplementation,weusetheLobattoquadraturerules[30].In
addition to simplifying volume integration, this integration scheme
also makes implementation of the ABD-equivalent graded model 5. Numericaltestresults
much easier: since an integration surface is an offset at a constant
distancefromthelaminate'smid-plane,coefficientsΛ ofthegraded In this section, we compare results of linear static analysis
ij
material (Eq. (13)) are also constant within that integration surface computed using ABD-equivalent materials to results in the stan-
and, therefore, need to be computed once. Integration over each dard reference text [5], as well as to results computed using
surface is done by first triangulating it, and then integrating the commercial software SolidWorks [31] and ANSYS [32]. Since we
obtained triangles using Gauss quadrature rules. The triangles are areusing3DmeshfreeFEAimplementation,wedonotexpectitto
constrainedtoconformtotheCartesiangrid,orinotherwords,each be more efficient than highly optimized commercial 2D FEA
trianglecompletelylieswithinanelementoftheCartesiangrid. solutions. Our only goal is to demonstrate the accuracy of ABD-
Whileintegratingoverthetriangles,wealsoneedtotransform equivalent material models in comparison with other analytical
thematerialmatrixQn fromitsprincipalcoordinatesystemtothe andcomputationalmethodsthatrelyontheABDformulation.This
elementcoordinatesystem.AsexplainedinFig.7,foreverytriangle, alsoexplainswhywedonotcomparetheaccuracyofourresults
we transform Qn once for the triangle's centroid, and use the to more accurate models based on higher-order lamination the-
transformed properties Qn0 for all the quadrature points of that oriesorfull3Delasticityformulations.
triangle. The transformation relation is explained in detail in We are considering four benchmark problems: a rectangular
plate,acylindricalshell,acylindricalroof,andalapjoint.Thefirst
threestructuresaresinglemulti-plylaminates,whilethelapjoint
is a bonded assembly of two laminates. For each of the test
structures,weconsiderthefollowingthreeconfigurations:
1. Cross-plylaminates½0=90(cid:5) ,
n
2. Angle-plylaminates½(cid:4)45=45(cid:5) ,
n
3. 50-plylaminates(AppendixD).
Cross-ply and angle-ply laminate configurations are often standard
configurations.The50-plylaminatewaschosentotestiftheproposed
method scales for a large number of plies. The plies in the 50-ply
laminatewereselectedrandomly,andtheirorientationsaregivenin
Appendix Appendix D. We restricted the number of plies to 50
because that was the maximum number of plies supported by
SolidWorks.Thethreelaminateconfigurationsareasymmetricalabout
theirmid-plane,andthereforeshowstretching–bendingcoupling,as
predictedbytheclassicallaminationtheory.Allpliesineachlaminate
areassumedtobemadeofoneofthetwomaterials:
Material 1: E ¼25:0e6psi, E ¼E ¼1:0e6psi, ν ¼0:25, ν ¼
1 2 3 12 23
Fig.6. LaminatestructureinFig.1replacedbyintegrationsurfaces. ν13¼0:0,G12¼G13¼5:0e5psiandG23¼2:0e5psi,
G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55 47
Fig.8. (A)Platewithgeometryparametersa¼10inandh¼1;0:1;0:01in,clampedfromallfoursideswithasurfacepressure,q¼100psi.(B)Platewitha¼10inand
h¼0:1in,fixedononeendwithaforce,F¼1:0e4lbf ontheoppositeend.
Table1
MaximumdisplacementvalueininchesforplateprobleminFig.8Afordifferentcases.
Laminate Thickness a=h ANSYS SW SnS–3-Ply SnS–Graded
(Numberofelements) 10k 1k 1k 3k 1k 3k
1 10 4.748e(cid:4)3 4.058e(cid:4)3 3.855e(cid:4)3 3.762e(cid:4)3 4.567e(cid:4)3 3.755e(cid:4)3
½0=90(cid:5)5 0.1 100 1:543eþ0 1:550eþ0 1:532eþ0 1:543eþ0 1:532eþ0 1:543eþ0
0.01 1000 1:510eþ3 1:552eþ3 1:145eþ3 1:661eþ3 1:145eþ3 1:662eþ3
1 10 5.094e(cid:4)3 4.360e(cid:4)3 4.057e(cid:4)3 4.152e(cid:4)3 4.880e(cid:4)3 4.050e(cid:4)3
½(cid:4)45=45(cid:5)5 0.1 100 1:629eþ0 1:620eþ0 1:597eþ0 1:611eþ0 1:598eþ0 1:611eþ0
0.01 1000 1:581eþ3 1:578eþ3 1:163eþ3 1:684eþ3 1:163eþ3 1:685eþ3
Material 2: E1¼7:5e6psi, E2¼E3¼2:0e6psi, ν12¼0:25, ν23¼ Table2
ν ¼0:0,G ¼1:25e6psi,andG ¼G ¼0:625e6psi, VonMisesstressatthemid-point(x¼0,y¼0)ofthetopfaceoftheplateinFig.8A.
13 12 13 23
ν
whereE, andGareYoung'sModulus,Poisson'sRatioandShear
Laminate Thickness a=h ANSYS SW SnS–3-Ply SnS–Graded
Modulus respectively, and indices 1, 2 and 3 represent the three
principalmaterialdirections.Thesematerialsareidenticaltothose (Numberofelements) 10k 1k 1k 3k 1k 3k
ν ν
usedin[5]thatweuseforcomparison,exceptweset and
23 13
tobezero. 1 10 2.42e3 2.43e3 2.50e3 2.40e3 3.20e3 2.35e3
½0=90(cid:5) 0.1 100 2.54e5 2.55e5 2.04e5 2.06e5 2.05e5 2.06e5
Finally we note that the elements available for analyzing 5
0.01 1000 2.50e7 2.60e7 1.60e7 1.94e7 1.60e7 1.93e7
laminatesinthethreesystemsaredifferent:
1 10 2.70e3 2.66e3 2.67e3 2.60e3 3.20e3 2.60e3
½(cid:4)45=45(cid:5) 0.1 100 2.25e5 2.26e5 1.80e5 1.82e5 1.80e5 1.82e5
5
Solidworks: We used two-dimensional parabolic triangular shell 0.01 1000 2.25e7 2.22e7 1.52e7 1.70e7 1.52e7 1.69e7
elements.
ANSYS: Weusedtwo-dimensionalShell181elementsforanalyz-
ingsinglelaminatestructures,andSolid186elementsfor and Solve can successfully capture bending in thin structures,
analyzingthelapjoint.Shel181are4nodeelementswith since conventional 3D basis functions tend to underestimate
3 displacement and 3 rotational degrees of freedom at bending deformations due to locking. The plate was made of
each node. A penalty method is used to relate the Material1,andconsistedof10plieslaidincross-plyandangle-ply
independent rotational degrees of freedom about the configurations.
normal (to the shell surface) with the in-plane compo- Table 1 compares the maximum displacements using ANSYS,
nents of displacements. Solid186 elements are 20-node SolidWorks (SW), and proposed method (SnS) for different lami-
layered solid elements that exhibit quadratic displace- nates. Tests were done for three different aspect ratios: thin (a/
mentbehavior[32]. h¼1000), moderately thick (a/h¼100), and thick (a/h¼10). For
Scan&Solve: Each benchmark problem is solved using 1000 and both cross-ply and angle-ply laminates with moderate thickness,
3000 second-order tri-variate B-spline functions on a SnS accurately predicts the maximum displacement values, and
uniform Cartesian non-conforming grid. The Lobatto theresultsfromallthesystemsareincloseagreement.Thereare
quadrature rule implies that 3 integration surfaces per morenoticeabledifferencesinthedisplacementscomputedbythe
ply for the 3-ply laminate model, and 4 integration threesystemsforthinandthicklaminates,e.g.,SnSandSWtend
surfaces for the quadratically graded laminate model, todifferby5%.Importantly,SnSisnotunder-predictingdisplace-
aresufficientforfullintegration. mentsforthinplates,suggestingthatlockingisnotanissue.
Wealsocompare,inTable2,theVonMisesstressesatthemid-
point ðx¼0;y¼0Þ of the top face for different plates. For both
5.1. Clampedrectangularplate cross-plyandangle-plylaminates,stressesfromdifferentmethods
areinagreementforthethickplatewithaspectratio10,butthereis
Our first benchmark problem is a plate clamped on all four somedeviationforotheraspectratios,whichincreases withplate
sides,andanormal pressureonthetopsurface (Fig.8A).Aplate aspectratio.Acomparisonoftheentirestressfieldforthemoderate
under these boundary conditions shows pure bending, with thicknessplate(aspectratio100)isshowninFig.9forcross-plyand
maximumdisplacementatthecenteroftheplate.Thisparticular angle ply laminates. The stress patterns match for the two lami-
problemwaschosentotestifsecondorderB-splinebasisinScan nates, and the highest stress due to stress concentration are also
48 G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55
Fig.9. FigurecomparesVonMisesstressfortopfaceoftheplateinFig.8withaspectratio100for(A1)½0;90(cid:5) analyzedinANSYS(Maxstress-5:83e5),(A2)½0;90(cid:5) analyzed
5 5
usingtheproposedmethod(Maxstress-4:91e5),(B1)½(cid:4)45;45(cid:5) analyzedinANSYS(Maxstress-4:08e5),(B2)½(cid:4)45;45(cid:5) analyzedusingtheproposedmethod(Maxstress-
5 5
3:44e5).
Fig.10. Figureshowsdeformedplateandcolormapofout-of-planedisplacementobtainedusing(A)ANSYS(B)Scan&SolvewithABD-equivalent3-plymaterialmodel,fora
laminateplatemadeof50plies.Asexpected,thein-planeloadleadstoout-of-planebending.(Forinterpretationofthereferencestocolorinthisfigurecaption,thereaderis
referredtothewebversionofthispaper.)
G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55 49
Fig.13. Figureshowingbarrelvaultproblem.Verticalpressureisq¼0:625psiand
thecurvedendsarefixed.Otherdimension:β¼801,R¼300in,a¼600in,h¼3;6
and15in.
Fig.11. ClampedcylinderwithPo¼2:04ksi,R¼20in,a¼20in,h¼1in.
z displacement are in agreement with 0:3730in in ANSYS and
0:3734ininSnS.
Wealsocarriedoutatimeanalysisinordertoestimatethenet
Table3 efficiency achieved using proposed method when using 3D FEA.
Maximumradialdeflection(e(cid:4)1in)forcylindricalshellinFig.11.Q4andQ9results
Complete analysis of the 50-ply laminate plate using the ABD-
from[5].SW-SolidWorks,SnS-ProposedMethod.
equivalent3-plymodeltook14.9s,outofwhich12.8swerespent
Laminate Q4 Q9 ANSYS SW SnS–3-Ply SnS–Graded integrating 9 integration surfaces (3 surfaces per ply). Therefore,
an average of 1.42s were spent integrating each surface. This
16 4 15k 1.2k 1k 3k 1k 3k
implies that integrating over 150 surfaces in the original 50-ply
½0=90(cid:5) 1.870 1.803 1.706 1.848 1.820 1.773 1.820 1.773 modelwouldrequireroughly215sforthesameanalysis.Thegain
½(cid:4)45=45(cid:5) – – 2.204 2.350 2.356 2.291 2.355 2.290 inefficiencyisevenhigherwhenusinggradedmaterialmodel,as
½0=90(cid:5)5 – – 1.719 1.830 1.814 1.776 1.815 1.776 itneedsonly4integrationsurfacesincomparisonto9forthe3-
½(cid:4)45=45(cid:5) – – 2.210 2.340 2.334 2.282 2.334 2.281
5 ply laminate. The total time taken for analysis was only 6.8s,
50-ply – – 2.350 2.542 2.465 2.424 2.464 2.424
decreasingthetotalcomputationcostofanalysisbymorethan30
times. Similar speedup should be expected in any 3D FEA of
within 15% of each other. Solid elements generally better capture laminatedstructuresrelyingonlayeredelements.
thestressesnearedgesastheyexplicitlymodelthefinitethickness
ofplates. 5.2. Clampedcylinderwithinternalpressure
Next,wetestthesamerectangularplate,butmadeof50random
plies.Theboundaryconditionsaredifferentfromtheprevioustests; Next, we consider a cylinder, made of Material 2, that is
the plate is under in-plane load of 1:0e4lbf on one end and clamped at the two ends and has an internal pressure of P .
o
clamped at the opposite end (Fig. 8B). This particular problem DetailsofgeometryandboundaryconditionsareshowninFig.11.
waschosentovalidatetwoimportantclaims:(a)proposedmethod ThemaximumradialdeflectionsfromReddy[5],ANSYS,Solid-
can handle large number of plies, (b) proposed ABD-equivalent Works, and proposed method are comparedin Table 3.Elements
material models successfully capture coupling behavior in lami- used by Reddy [5] are Q4 elements, a four-node (linear) quad-
nates that are asymmetrical about the mid-plane. Due to this rilateralelement,andQ9elements,anine-node(quadratic)quad-
coupling, the in-plane load F will lead to bending and produces rilateral element. The total number of elements used in different
out-of-plane deformation. The proposed method did capture this methodsarespecifiedinthesecondrowofTable3.Thetestswere
coupling phenomena accurately as shown in Fig.10, which com- performed forone aspect ratio,but 5 differentlaminates.Results
paresthezdisplacementfieldsfromSnSandANSYS.Themaximum fromall the methods are clearly in agreement. Fig.12 shows the
Fig.12. Figureshowsdeformedcylinderandthecolormapofradialdisplacementobtainedusing(A)ANSYS(B)Presentmethod,foralaminatecylindermadeof50plies.
(Forinterpretationofthereferencestocolorinthisfigurecaption,thereaderisreferredtothewebversionofthispaper.)
50 G.Kumar,V.Shapiro/FiniteElementsinAnalysisandDesign106(2015)41–55
Table4
Barrelvaultproblem-maximumdeflection(in.)forcross-plyandangleplylaminates.
Laminate R/h Reddy ANSYS SW SnS–3-Ply SnS–Graded
(No.elements) 16 10k 1.2k 1 k 3k 1k 3k
100 2:339eþ0 2:407eþ0 2:460eþ0 2:307eþ0 2:415eþ0 2:307eþ0 2:416eþ0
½0=90(cid:5) 50 5.082e(cid:4)1 5.291e(cid:4)1 5.659e(cid:4)1 5.480e(cid:4)1 5.810e(cid:4)1 5.496e(cid:4)1 5.503e(cid:4)1
20 7.292e(cid:4)2 7.449e(cid:4)2 7.560e(cid:4)2 7.877e(cid:4)2 8.067e(cid:4)2 7.956e(cid:4)2 7.016e(cid:4)2
100 3:597eþ0 3:871eþ0 3:866eþ0 3:411eþ0 3:743eþ0 3:413eþ0 3:699eþ0
½(cid:4)45=45(cid:5) 50 6.760e(cid:4)1 7.652e(cid:4)1 7.170e(cid:4)1 6.675e(cid:4)1 7.157e(cid:4)1 6.671e(cid:4)1 7.780e(cid:4)1
20 1.205e(cid:4)1 1.397e(cid:4)2 1.130e(cid:4)1 1.061e(cid:4)1 1.127e(cid:4)1 1.131e(cid:4)1 1.386e(cid:4)1
100 1:415eþ0 1:434eþ0 1:564eþ0 1:542eþ0 1:593eþ0 1:540eþ0 1:593eþ0
½0=90(cid:5) 50 2.940e(cid:4)1 2.979e(cid:4)1 3.270e(cid:4)1 3.335e(cid:4)1 3.412e(cid:4)1 3.337e(cid:4)1 3.412e(cid:4)1
5
20 5.234e(cid:4)2 5.246e(cid:4)2 5.370e(cid:4)2 5.361e(cid:4)2 5.398e(cid:4)2 5.361e(cid:4)2 5.399e(cid:4)2
100 1:818eþ0 1:836eþ0 1:955eþ0 1:821eþ0 1:912eþ0 1:821eþ0 1:912eþ0
½(cid:4)45=45(cid:5) 50 4.096e(cid:4)1 4.082e(cid:4)1 4.089e(cid:4)1 3.796e(cid:4)1 3.940e(cid:4)1 3.799e(cid:4)1 3.940e(cid:4)1
5
20 1.004e(cid:4)1 9.727e(cid:4)2 8.959e(cid:4)2 7.856e(cid:4)2 8.009e(cid:4)2 7.857e(cid:4)2 8.088e(cid:4)2
50-ply 100 – 1.4106e0 1.478e0 1.391e0 1.445e0 1.391e0 1.445e0
Fig.14. Figureshowsdeformedshellandcolormapofout-of-planedisplacementobtainedusing(A)ANSYS(B)Presentmethod,foralaminateplatemadeof50plies.(For
interpretationofthereferencestocolorinthisfigurecaption,thereaderisreferredtothewebversionofthispaper.)
Table5
MaximumdisplacementvalueininchesforlapjointinFig.15.Secondrowshows
thenumberofelementsused.
Laminate ANSYS SnS–3-Ply SnS–Graded
760 1k 3k 1k 3k
½0;90(cid:5) 1.251 1.147 1.212 1.147 1.212
5
½(cid:4)45;45(cid:5) 7.975 6.540 7.641 6.540 7.641
5
50Plies 3.304 2.937 3.194 2.938 3.195
deformed shape, aswell as, the color map of displacement value
for the cylinder made of 50 plies. Maximum radial displacement
usingANSYSis0:2350in,andusing3-plymaterialmodelinScan
Fig.15. Figureshowsalapjointmadeoftwolaminatesidenticalingeometry,with
andSolveis0:2419in,whichislessthan3%difference. a¼10in,b¼2in,c¼2:5inandh¼0:1in.Theleftendisfullyfixed,andtheright
endisallowedtoslideinthexdirection.AforceF¼10e4lbisalsoappliedonthe
5.3. BarrelVaultproblem rightend.
In this section, we consider another popular benchmark shell laminate (see Fig. 14). The maximum displacement values from
problemknownastheBarrelvaultproblem[5]:acylindricalroof Scan&Solve using 3000 elements and ANSYS are within 2.5% of
underitsownweight.ThestructureismadeofMaterial1,andthe eachother(Table5).
detailedboundaryconditionsareshowninFig.13. As before, when analysing 50-ply laminates using layered 3D
Inreference[5],maximumverticaldisplacementsaregivenfor elements,the3-plylaminateimprovesefficiencyroughly15times,
cross-plyandangle-plylaminateswithdifferentaspectratios.We whereas the graded material was 30 times more efficient than
usethesevaluesforcomparingvaluesobtainedusingSW,ANSYS, usingtheactuallaminatelayup.
andproposedmethodinTable4.ReddyusedQ81elements,which
areeighthorderelement(p¼8)with405degreesoffreedom.For 5.4. Multi-laminatestructure—lapjoint
othermethods,elementsaresameasbefore.Again,thereisafairly
closeagreementbetweentheresultsfromReddy[5],SolidWorks, Finally, to show that proposed method can be extended to
ANSYS and proposed method in all cases, including the 50-ply structures made of multiple laminates, we analyze a lap joint