復合材料層合結構應力分析

復合材料層合結構應力分析

隨著工業(yè)和技術的進步，一體化板的外殼結構在航空航天等領域得到了越來越多的應用。關于復合材料層合板殼的有限元分析一直是一個令人關注的問題。為了克服傳統(tǒng)的退化殼元難以處理大的轉動增量的問題,最近幾年人們把研究目光轉向了固體殼元,但隨之而來的是固體殼元的厚度自鎖問題?，F(xiàn)有的解決方法有采用平面應力假設的本構方程、加強假定應變列式和雜交應力元列式等,但這些方法都存在著各種各樣的問題。采用平面應力假設的本構方程在處理厚度變化較大的層合結構時會導致厚度方向正應力的明顯不連續(xù)分布,加強假定應變列式同樣不能得到連續(xù)的厚度方向的應力分布,且它要在單元一級上消除內(nèi)部參數(shù),降低了計算效率,雜交應力元列式雖可得到連續(xù)的厚度方向的應力分布,但有時會影響位移的計算精度。本文中則引入廣義應力的概念,建議了一種新的改進本構陣的方法,克服了厚度自鎖問題,尤其適用于復合材料層合板殼結構的數(shù)值分析。正交化方法是近年來多變量有限元研究領域取得的一個主要成果,它能明顯地提高低階單元的計算效率。在高階單元的研究方面,S.W.Lee等建議了一種雜交穩(wěn)定列式,將應變插值函數(shù)分為低階和高階兩部分,與低階項對應的剛度陣和相應的減縮積分單元等價,與高階項對應的部分用來克服單元的零能模式,但其建議的列式中對應高階項的部分還要采用高階的數(shù)值積分,且要求一滿陣的逆。本文列式中借鑒了低階雜交元列式中的正交化方法,消除了低階和高階項的交叉項,并提出一個精化措施,只用低階積分即可求出高階穩(wěn)定陣,且避免了矩陣求逆運算,顯著提高了計算效率。1根據(jù)ddlagrangaan的位移法，有限列方程1.1增量應變變分的計算九節(jié)點固體殼元(每個節(jié)點六個自由度)的坐標、位移采用下列的標準形式:x(ξ,η,ζ)=9∑i=1Νi(ξ,η)(xi+ζvi)=Ν(ξ,η)ˉx+ζΝ(ξ,η)ˉv=x0(ξ,η)+ζxn(ξ,η)(1)U(ξ,η,ζ)=9∑i=1Νi(ξ,η)[Ui+ζ(v′i-vi)]=U0(ξ,η)+ζUn(ξ,η)(2)在傳統(tǒng)的退化殼元中,假定變形前后殼體的厚度保持不變,故v′i-vi(v′i、vi分別為變形前、后的中面法線)由兩個轉動矢量來表示,目前的列式中拋棄了這一假設,因此v′i-vi為一有三個分量的矢量。ΔU(ξ,η,ζ)=9∑i=1Νi(ξ,η)(ΔUi+ζΔVi)=Ν(ξ,η)ΔˉU+ζΝ(ξ,η)ΔˉV=ΔU0(ξ,η)+ζΔUn(ξ,η)(3)其中:xi、Ui、ΔUi分別是節(jié)點坐標、節(jié)點位移和增量節(jié)點位移矢量,Ni是二維拉格朗日插值函數(shù)?！={x1?x9}={X1+U1?X9+U9},ˉV={V1?V9},ΔˉU={ΔU1?ΔU9},ΔˉV={ΔV1?ΔV9},Ν=[Ν1,?,Ν9]由應變定義,t時刻的Almansi逆變應變和增量Almansi逆變應變?yōu)棣舏=εii=xT′iU′i-UT′iU′i/2,2εij=xΤ′iU′j+xΤ′jU′i-(UΤ′iU′j+UΤ′jU′i)/2(4)Δεi=Δεii=XΤ′iΔU′i+ΔUΤ′iΔU′i/2,2Δεij=XΤ′iΔU′j+XΤ′jΔU′i-(ΔUΤ′iΔU′j+ΔUΤ′jΔU′i)/2(5)其中:i,j=ξ,η,ζ,且i≠j。ε=={εx′εy′2εx′y′}=Τ={εξεη2εξη},εt={2εz′x′2εz′y′}=Τt{2εζξ2εζη},εz=Τzεζ(6)Δε=={Δεx′Δεy′2Δεx′y′}=Τ={ΔεξΔεη2Δεξη},Δεt={2Δεz′x′2Δεz′y′}=Τt{2Δεζξ2Δεζη},Δεz=ΤzΔεζ(7)在增量應變中忽略關于ζ的二次項,得Δε==Δεm+ζΔεb=(BLm+BΝm/2)ΔˉU+ζ(BLb+BΝb/2)ΔˉU,Δεz=(BLz+BΝz/2)ΔˉU,Δεt=Δεo+ζΔεn=(BLo+BΝo/2)ΔˉU+ζ(BLn+BΝn/2)ΔˉU(8)顯然增量應變的變分可寫成δΔελ=(BLλ+BΝλ)δΔˉUδΔεz=(BLz+BΝz)δΔˉU(9)其中:λ=m,b,o,n。當前時刻的Cauchy應力和增量Kirchhoff應力可寫成{σ=σz}=[C=C×CΤ×Cz]{ε=εz},σt=Ctεt,{ΔσΔσz}=[C=C×CΤ×Cz]{Δε=Δεz},Δσt=CtΔεt(10)每個單元的最小勢能原理可寫為ΠeΡ=12n∑l=11∫-11∫-11∫-1[{ε=+Δε=εz+Δεz}Τ[C=C×CΤ×Cz]?{ε=+Δε=εz+Δεz}+(εt+Δεt)ΤCt(εt+Δεt)]?Jh1hdξdηdζl-(ˉU+ΔˉU)ΤF(11)其中:n是復合材料的總鋪層數(shù),h是所有層的總厚度,hl是各層的厚度,J是當前構形與等參構形間的雅可比轉換矩陣,F是外部載荷形成的等效節(jié)點力。為提高計算效率,近似取ˉJ=J|ζl=0,同時代入式(10),有ΠΡe=12?(ε⊥+Δε⊥)ΤC?⊥(ε⊥+Δε⊥)+(εΤ+ΔεΤ)ΤCΤ(εΤ+ΔεΤ)?-(Uˉ+ΔUˉ)ΤF(12)其中:ε⊥={εmεzεb},Δε⊥={ΔεmΔεzΔεb},C?⊥=12∑l=1n∫-11[C=C×ζlC=C×ΤCzζlC×ΤζlC=ζlC×ζl2Cz]hlhdζl,εt={εoεn},Δεt={ΔεoΔεn},CΤ=12∑l=1n∫-11[CtζlCtζlCtζl2Ct]hlhdζl,?o?=2∫-11∫-11(o)Jˉdξdη此外,我們定義當前時刻的廣義應力如下(增量廣義應力類似):Ν⊥={ΝΤΜ}=12∑l=1n∫-11{σ=σzζlσ=}hlhdζl=12∫-11{C=(εm+ζlεb)+C×εzC×Τ(εm+ζlεb)+CzεzζlC=(εm+ζlεb)+ζlC×εz}hlhdζl=C?⊥ε⊥QΤ={QoQn}=12∑l=1n∫-11{σtζlσt}hlhdζl=12∫-11{Ct(εo+ζlεn)ζlCt(εo+ζlεn)}hlhdζl=CΤεΤ將廣義應力和增量廣義應力代入變分泛函(12),有δΠΡe=?δ(Δε⊥)Τ(Ν⊥+C?⊥Δε⊥)+δ(ΔεΤ)Τ(ΝΤ+CΤΔεΤ)?-δ(ΔUˉ)ΤF(13)經(jīng)線性化處理,有δΠΡe=δ(ΔUˉ)Τ[(ΚL+ΚΝ)ΔUˉ+R-F](14)其中單元剛度矩陣為ΚL=?B⊥LΤC?⊥B⊥L?,ΚΝΔUˉ=?B⊥ΝΤΝ⊥+BΤΝΤQΤ?(15)等效節(jié)點力為R=?B⊥LΤΝ⊥+BΤLΤQΤ?(16)其中:B⊥L[BmLBzLBbL],B⊥L=[BmΝBzΝBbΝ],BΤL=[BoLBnL],BΤΝ=[BoΝBnΝ]1.2改進的本構陣將式(10)改寫為{σ=εz}=[AB-BΤD]{ε=σz}=[S=-1-S=-1S×-S×ΤS=-1Sz-S×ΤS=-1S×]{ε=σz}(17)為了克服厚度自鎖,定義εz=12∑l=1n∫-11([-BΤD]{ε=σz})hlhdζl=12∑l=1n∫-11([-BΤD]{εm+ζlεbσz})hlhdζl=[-B0ΤD0-B1Τ]{εmσzεb}(18){ΝΜ}=12∑l=1n∫-11{σ=ζlσ=}hlhdζl=12∑l=1n∫-11[ABζlAζlB]{εm+ζlεbσz}hlhdζl=[A0B0A1A1B1A2]{εmσzεb}(19)其中:[A0A1A2B0B1D0]=12∑l=1n∫-11[AζlAζl2ABζlBD]hlhdζl由式(10),(18),(19)得改進的本構陣為C⊥=[A0+B0B0Τ/D0B0/D0A1+B0B1Τ/D0B0Τ/D01/D0B0Τ/D0A1+B1B0Τ/D0B1/D0A2+B1B1Τ/D0](20)2單剛單干單階項與減縮積分元等價的計算為了表述的簡潔性,引入記號ue001φ,ue001φ代表⊥和T(⊥和T所代表的意義與1.1中所述相同)。Hellinger-Reissner變分原理可寫成ΠΗRe=∑φ=⊥,Τ?-(σφΤ+ΔσφΤ)Cφ-1(σφ+Δσφ)/2+(σφΤ+ΔσφΤ)(εφ+Δεφ)?-(Uˉ+ΔUˉ)ΤFδΠΗRe=∑φ=⊥,Τ?δ(ΔσφΤ)[Cφ-1(σφ+Δσφ)-(εφ+Δεφ)]+(σφΤ+ΔσφΤ)δ(Δεφ)?-δ(ΔUˉ)ΤF(21)廣義應力分為低、高階兩部分,即:Δσue001φ=Pue001φLΔβue001φL+Pue001φHΔβue001φH,其中應力插值函數(shù)滿足下列正交性條件:〈PTue001φLC-1ue001φPue001φH〉=0,引入記號φ,φ代表L和H(分別代表低、高階項,以下各式中同),則應力插值函數(shù)可簡化為:Δσφ=∑φΡφφΔβφφ,將應力插值函數(shù)和式(9)、(20)代入式(21)有δΠΗRe=∑φφ?-δ(ΔβφφΤ)[ΗφφΔβφφΤ-GφφΔUˉ-Ζφφ]+δ(ΔUˉΤ)GφφΔβφφΤ?+δ(ΔUˉ)Τ(ΚσΔUˉ+Q-F)(22)其中:Ηφφ=?ΡφφΤCφ-1Ρφφ?,Gφφ=?ΡφφΤBφ?,Ζφφ=?ΡφφΤ(εφ-Cφ-1σφ)?(23)由式(22),得獨立假定的增量應力參數(shù)為Δβφφ=Ηφφ-1(GφφΔUˉ+Ζφφ)(24)將式(24)代入式(22),得δΠΗRe=δ(ΔUˉΤ){(∑φ,φGφφΤΗφφ-1Gφφ+Κσ)ΔUˉ-F+Q+∑φ,φGφφΤΗφφ-1Ζφφ}(25)為了提高計算效率,低階項部分選擇如下的插值函數(shù),使本雜交元列式的單剛低階項部分與減縮積分元等價Ρ⊥L=[Ρ1Ι7Ρ2Ι7Ρ3Ι7Ρ4Ι7],ΡΤL=[Ρ1Ι4Ρ2Ι4Ρ3Ι4Ρ4Ι4]其中:Ρ1=(1-ξ3)(1-η3)/4,Ρ2=(1+ξ3)(1-η3)/4,Ρ3=(1+ξ3)(1+η3)/4,Ρ4=(1-ξ3)(1+η3)/4則式(23)可簡化為GφL=?ΡφLΤBφ?L=[Jˉ1Bφ1?Jˉ4Bφ4]?ΖφL=[Jˉ1(εφ-Cφ-1σφ)1?Jˉ4(εφ-Cφ-1σφ)4]?ΗφL=?ΡφLΤCφ-1ΡφL?L=diag[Jˉ1Cφ1-1,?,Jˉ4Cφ4-1]其中∶?o?L=∑i=14Jˉi(o)i?(o)i為矩陣或變量‘o’在第i個2×2高斯積分點的值。所以有GφLΤΗφL-1GφL=?BφΤCφBφ?L,GφLΤΗφL-1ΖφL=?BφΤ(Cφεφ-σφ)?L故式(25)可變?yōu)棣摩唉e=δ(ΔUˉΤ){[∑φ?BφΤCφBφ?L+GφΗΤΗφΗ-1GφΗΤ)+Κσ]ΔUˉ-F+Q+∑φ[?BφΤ(Cφεφ-σφ)?L+GφΗΤΗφΗ-1ΖφΗ]}(26)為克服減縮積分元的零能模式,選取下列高階應力插值函數(shù):Ρ⊥Η=1J[Τ=-1Ρ⊥03×201×201×203×2Τ=-1Ρ⊥],ΡΤΗ=1J[Τt-1Ρ⊥02×102×1Τt-1ΡΤ],Ρ⊥=[Ρξ00Ρη00],ΡΤ=[ΡξΡη],Ρξ=ξ(η2-1/3)?Ρη=η(ξ2-1/3)則G⊥Η=?Ρ⊥ΗΤB⊥?L=[GmξGmηGbξGbη]=[xˉΜmξxˉΜmηxˉΜbξxˉΜbη]?Η⊥Η=?Ρ⊥ΗΤC⊥-1Ρ⊥Η?=?1J02[Ρ⊥03×201×201×203×2Ρ⊥]ΤS⊥[Ρ⊥03×201×201×203×2Ρ⊥]?GΤΗ=?ΡΤΗΤBΤ?=[GtoGtn]=xˉΤ[ΜtoΜtn],ΗΤΗ=?ΡΤΗΤCΤ-1ΡΤΗ?=?1J02[Ρ⊥02×102×1ΡΤ]ΤSΤ[ΡΤ02×102×1ΡΤ]?其中:S⊥=diag{Τ=-1,0,Τ=-1}?C⊥-1?diag{Τ=-Τ,0,Τ=-Τ},SΤ=diag{Τt-1,Τt-1}?CΤ-1?diag{Τt-Τ,Τt-Τ}其中:Mmξ、Mmη、Mbξ、Mbη、Mto、Mtn可利用Maple(WaterlooMapleInc.)等軟件推得顯式。近似取:Jˉ≈J0=Jˉ|ξ=η=0,S⊥≈Sˉ⊥=Sˉ⊥|ξ=η=0?SΤ≈SˉΤ=SˉΤ|ξ=η=0,則H⊥H,HTH均可由2×2的高斯積分求出。代入式Q=∑φ?BφΤσφ?L?ΖφΗ=0,式(26)變?yōu)棣摩唉e=δ(ΔUˉΤ){[∑φ?BφΤCφBφ?L+GφΗΤΗφΗ-1GφΗ+Κσ]ΔUˉ-F+∑φ?BφΤCφεφ?L}(27)3計算值的示例3.1材料及結構幾何四邊固支圓柱殼受如圖1所示均勻分布徑向力作用,其幾何尺寸為:R=2540mm,L=508mm,θ=0.1°,殼體總厚度為h=2.54mm,考慮[0°/90°]鋪層,所使用的材料常數(shù)為:E1=25×106N/mm2,E2=E3=106N/mm2,G12=G31=0.5×106N/mm2,G23=0.2×106N/mm2,ν12=ν23=ν13=0.25。由于結構和載荷的對稱性,僅將四分之一殼體離散為2×2網(wǎng)格。圖2反映了使用十個加載步時,變形隨載荷變化的情況,并與文獻中的結果進行了比較。3.2流固體的結構和載荷一端固支圓柱殼,如圖3所示自由端曲邊中點受一集中力作用,其幾何尺寸為:R=101.6mm,h=3mm,L=304.8mm,材料為橫向各向同性,彈性常數(shù)為:EL=2068.5N/mm2,ET=517.125N/mm2,νLT=νTL=0.3,GLT=795.6N/mm2,考慮90°/0°/90°和0°/90°/0°兩種鋪層。殼體所受外力為Pmax=36N,由于結構和載荷的對稱性,僅將二分之一半圓柱殼體離散為8×8網(wǎng)格。圖4反映了A點的撓度隨載