復(fù)合材料層合結(jié)構(gòu)應(yīng)力分析

復(fù)合材料層合結(jié)構(gòu)應(yīng)力分析

隨著工業(yè)和技術(shù)的進步，一體化板的外殼結(jié)構(gòu)在航空航天等領(lǐng)域得到了越來越多的應(yīng)用。關(guān)于復(fù)合材料層合板殼的有限元分析一直是一個令人關(guān)注的問題。為了克服傳統(tǒng)的退化殼元難以處理大的轉(zhuǎn)動增量的問題,最近幾年人們把研究目光轉(zhuǎn)向了固體殼元,但隨之而來的是固體殼元的厚度自鎖問題?，F(xiàn)有的解決方法有采用平面應(yīng)力假設(shè)的本構(gòu)方程、加強假定應(yīng)變列式和雜交應(yīng)力元列式等,但這些方法都存在著各種各樣的問題。采用平面應(yīng)力假設(shè)的本構(gòu)方程在處理厚度變化較大的層合結(jié)構(gòu)時會導(dǎo)致厚度方向正應(yīng)力的明顯不連續(xù)分布,加強假定應(yīng)變列式同樣不能得到連續(xù)的厚度方向的應(yīng)力分布,且它要在單元一級上消除內(nèi)部參數(shù),降低了計算效率,雜交應(yīng)力元列式雖可得到連續(xù)的厚度方向的應(yīng)力分布,但有時會影響位移的計算精度。本文中則引入廣義應(yīng)力的概念,建議了一種新的改進本構(gòu)陣的方法,克服了厚度自鎖問題,尤其適用于復(fù)合材料層合板殼結(jié)構(gòu)的數(shù)值分析。正交化方法是近年來多變量有限元研究領(lǐng)域取得的一個主要成果,它能明顯地提高低階單元的計算效率。在高階單元的研究方面,S.W.Lee等建議了一種雜交穩(wěn)定列式,將應(yīng)變插值函數(shù)分為低階和高階兩部分,與低階項對應(yīng)的剛度陣和相應(yīng)的減縮積分單元等價,與高階項對應(yīng)的部分用來克服單元的零能模式,但其建議的列式中對應(yīng)高階項的部分還要采用高階的數(shù)值積分,且要求一滿陣的逆。本文列式中借鑒了低階雜交元列式中的正交化方法,消除了低階和高階項的交叉項,并提出一個精化措施,只用低階積分即可求出高階穩(wěn)定陣,且避免了矩陣求逆運算,顯著提高了計算效率。1根據(jù)ddlagrangaan的位移法，有限列方程1.1增量應(yīng)變變分的計算九節(jié)點固體殼元(每個節(jié)點六個自由度)的坐標(biāo)、位移采用下列的標(biāo)準(zhǔn)形式: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分別為變形前、后的中面法線)由兩個轉(zhuǎn)動矢量來表示,目前的列式中拋棄了這一假設(shè),因此v′i-vi為一有三個分量的矢量。ΔU(ξ,η,ζ)=9∑i=1Νi(ξ,η)(ΔUi+ζΔVi)=Ν(ξ,η)ΔˉU+ζΝ(ξ,η)ΔˉV=ΔU0(ξ,η)+ζΔUn(ξ,η)(3)其中:xi、Ui、ΔUi分別是節(jié)點坐標(biāo)、節(jié)點位移和增量節(jié)點位移矢量,Ni是二維拉格朗日插值函數(shù)?！={x1?x9}={X1+U1?X9+U9},ˉV={V1?V9},ΔˉU={ΔU1?ΔU9},ΔˉV={ΔV1?ΔV9},Ν=[Ν1,?,Ν9]由應(yīng)變定義,t時刻的Almansi逆變應(yīng)變和增量Almansi逆變應(yīng)變?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)在增量應(yīng)變中忽略關(guān)于ζ的二次項,得Δε==Δε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)顯然增量應(yīng)變的變分可寫成δΔελ=(BLλ+BΝλ)δΔˉUδΔεz=(BLz+BΝz)δΔˉU(9)其中:λ=m,b,o,n。當(dāng)前時刻的Cauchy應(yīng)力和增量Kirchhoff應(yīng)力可寫成{σ=σ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是復(fù)合材料的總鋪層數(shù),h是所有層的總厚度,hl是各層的厚度,J是當(dāng)前構(gòu)形與等參構(gòu)形間的雅可比轉(zhuǎn)換矩陣,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η此外,我們定義當(dāng)前時刻的廣義應(yīng)力如下(增量廣義應(yīng)力類似):Ν⊥={ΝΤΜ}=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ΤεΤ將廣義應(yīng)力和增量廣義應(yīng)力代入變分泛函(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改進的本構(gòu)陣將式(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)得改進的本構(gòu)陣為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)廣義應(yīng)力分為低、高階兩部分,即:Δσue001φ=Pue001φLΔβue001φL+Pue001φHΔβue001φH,其中應(yīng)力插值函數(shù)滿足下列正交性條件:〈PTue001φLC-1ue001φPue001φH〉=0,引入記號φ,φ代表L和H(分別代表低、高階項,以下各式中同),則應(yīng)力插值函數(shù)可簡化為:Δσφ=∑φΡφφΔβφφ,將應(yīng)力插值函數(shù)和式(9)、(20)代入式(21)有δΠΗRe=∑φφ?-δ(ΔβφφΤ)[ΗφφΔβφφΤ-GφφΔUˉ-Ζφφ]+δ(ΔUˉΤ)GφφΔβφφΤ?+δ(ΔUˉ)Τ(ΚσΔUˉ+Q-F)(22)其中:Ηφφ=?ΡφφΤCφ-1Ρφφ?,Gφφ=?ΡφφΤBφ?,Ζφφ=?ΡφφΤ(εφ-Cφ-1σφ)?(23)由式(22),得獨立假定的增量應(yīng)力參數(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)為克服減縮積分元的零能模式,選取下列高階應(yīng)力插值函數(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材料及結(jié)構(gòu)幾何四邊固支圓柱殼受如圖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。由于結(jié)構(gòu)和載荷的對稱性,僅將四分之一殼體離散為2×2網(wǎng)格。圖2反映了使用十個加載步時,變形隨載荷變化的情況,并與文獻中的結(jié)果進行了比較。3.2流固體的結(jié)構(gòu)和載荷一端固支圓柱殼,如圖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,由于結(jié)構(gòu)和載荷的對稱性,僅將二分之一半圓柱殼體離散為8×8網(wǎng)格。圖4反映了A點的撓度隨載