復(fù)合材料層合殼體熱彈耦合簡(jiǎn)化模型

鐘軼峰，陳 磊，余文斌，張亮亮

（1．重慶大學(xué)a．土木工程學(xué)院；b．山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，重慶400045；2．美國(guó)猶他州立大學(xué) 機(jī)械與航空航天工程系，洛根84322）

近20年來(lái)，先進(jìn)復(fù)合材料結(jié)構(gòu)因其高強(qiáng)度、高模量、可設(shè)計(jì)性等優(yōu)點(diǎn)已廣泛應(yīng)用于航天航空、機(jī)械、土木等領(lǐng)域。由于復(fù)合材料不同成分、不同纖維方向的熱膨脹系數(shù)相差很大，相比于各向同性材料對(duì)溫度的改變更加敏感，復(fù)合材料板／殼結(jié)構(gòu)的熱響應(yīng)已成為近年來(lái)復(fù)合材料力學(xué)熱點(diǎn)課題之一［1－2］，出現(xiàn)了各種復(fù)合材料板／殼理論。Ng等［3］應(yīng)用古典層合理論（CLT）分析了復(fù)合材料圓柱殼的熱殘余應(yīng)力和自由振動(dòng)，由于忽略了橫向剪切效應(yīng)，該理論僅對(duì)薄殼有效；為克服CLT的缺陷，Pradyumna等［4］基于一階剪切變形理論（FOSDT）分析了復(fù)合材料層合殼體在熱環(huán)境下非線性動(dòng)態(tài)穩(wěn)定性，但模擬彎曲時(shí)數(shù)值計(jì)算結(jié)果嚴(yán)重失真，產(chǎn)生“橫剪自鎖”現(xiàn)象，其原因在于總應(yīng)變能中包含橫向剪切應(yīng)變能的項(xiàng)在量級(jí)上不正確，沿厚度方向變化的剪切應(yīng)變線性假設(shè)與上下表面剪應(yīng)力為零自相矛盾，需人為地引入與幾何形狀和材料相關(guān)的剪切修正因子來(lái)解決這一矛盾；為提高沿厚向的應(yīng)力／應(yīng)變預(yù)測(cè)精度，各國(guó)學(xué)者紛紛把研究重點(diǎn)轉(zhuǎn)向了高階剪切變形理論，Li等［5］用高階殼理論分析了復(fù)合材料圓柱殼在熱環(huán)境下的非線性屈曲和后屈曲性能，Nosier等［6］基于層合理論確定濕熱環(huán)境下正交鋪層／角鋪層殼體的局部位移函數(shù)和層間應(yīng)力，但仍無(wú)法準(zhǔn)確預(yù)測(cè)面內(nèi)的不連續(xù)性和沿厚度方向的橫向位移分量。

由上述分析可知，現(xiàn)有板殼熱彈性分析理論大多基于位移場(chǎng)假設(shè)，無(wú)法較好地反映層合板殼的三維效應(yīng)和鋪層之間的相互作用。筆者嘗試采用變分漸近法［7－8］分析層合圓柱殼的單向耦合熱彈性問(wèn)題。通過(guò)對(duì)降維模型近似能量中變分項(xiàng)的漸近修正，得到與原三維能量盡可能接近的近似能量，從而構(gòu)建一種無(wú)需先驗(yàn)性假設(shè)且便于工程應(yīng)用的簡(jiǎn)化殼體模型，為復(fù)合材料層合殼體在工程中的實(shí)際應(yīng)用提供有價(jià)值的參考。

1 復(fù)合材料層合殼體三維能量表達(dá)式

基于Hamilton擴(kuò)展原理，三維結(jié)構(gòu)的彈性動(dòng)力性能可表述為

式中：t1和t2是任意兩固定時(shí)間；K和U分別是動(dòng)能和Helmholtz自由能；δˉW是外力所做虛功，上劃線表明不需要對(duì)該項(xiàng)精確變分。

圖1 層合殼體幾何構(gòu)型及坐標(biāo)系

由于僅考慮單向耦合熱彈性問(wèn)題，因殼體變形引起的溫度改變忽略不計(jì)，得到無(wú)二次項(xiàng)的Helmholtz自由能泛函［9］：

式中：D為6×6階三維材料矩陣；α為6×1階三維熱膨脹系數(shù)列陣；T為相對(duì)于零應(yīng)力狀態(tài)的溫差；為殼體幾何修正系數(shù)，其中kαβ為殼體面外曲率，因選用坐標(biāo)系為圖1所示曲率線，k12＝k21＝0。

Γ為三維應(yīng)變場(chǎng)［10］，可定義為

式中：為相應(yīng)的算子；

其中：εαβ，καβ統(tǒng)稱二維廣義應(yīng)變，其階數(shù)分別用n，n／h表示；w i為未知翹曲函數(shù)，若將參考面定義為殼體中面，則w i須滿足如下3個(gè)約束

殼體參考面上質(zhì)點(diǎn)的動(dòng)能和外力所做虛功可分別表示為

式中：μ為質(zhì)量密度；U i為沿坐標(biāo)軸x i的位移，其上的點(diǎn)表示對(duì)時(shí)間求導(dǎo)；＝－δU3；α為虛擬旋轉(zhuǎn)；廣義力f i和力矩mα可定義為

式中：p i，τi，βi分別為體力、殼體頂／底面表面力。

由式（1）、（6）和式（7），殼體參考面的動(dòng)力性能可改寫為

由于式（9）中存在未知翹曲函數(shù)w i，若直接求解該式，將會(huì)遇到與求解原三維問(wèn)題同樣的困難。常規(guī)做法是對(duì)翹曲場(chǎng)作先驗(yàn)性假設(shè)，由于復(fù)合材料的各向異性和非均質(zhì)性，這種假設(shè)可能會(huì)產(chǎn)生較大的誤差。變分漸近法可利用殼體固有的小參數(shù)對(duì)式（9）的變分項(xiàng)進(jìn)行漸近分析得到w i，從而構(gòu)建與原三維模型盡可能接近的二維殼體模型。

2 降維方法及近似能量推導(dǎo)

對(duì)殼體結(jié)構(gòu)，可選擇h／l?1，h／R?1作為變分漸近計(jì)算所需的小參數(shù)，由此可估計(jì)各荷載階數(shù)為

式中μ為彈性材料常量的階數(shù)。

利用h／l，h／R可將三維能量漸近擴(kuò)展為如下形式的系列二維近似能量泛函

為能處理體多層殼結(jié)構(gòu)，并與二維有限元求解器相銜接，可將三維翹曲場(chǎng)離散為一維有限單元形式，

式中：S為形函數(shù)；V為沿橫法線方向的翹曲場(chǎng)節(jié)點(diǎn)值。

將式（11）代入式（9），得到離散形式能量泛函為

式中：L為荷載相關(guān)項(xiàng)；新引入的與幾何形狀和材料屬性有關(guān)的變量為

式（5）翹曲約束的離散形式可表示為

式中：H＝ ［STS］；ψ為零初始曲率E0的正交化核心矩陣，ψTHψ＝I。這樣，未知翹曲函數(shù)的求解問(wèn)題轉(zhuǎn)化為式（14）約束下式（12）最小化問(wèn)題。

2．1 零階近似

應(yīng)用變分漸近法，需根據(jù)不同階數(shù)找到泛函的主導(dǎo)項(xiàng)。式（12）中與未知翹曲函數(shù)有關(guān)的零階近似主導(dǎo)項(xiàng)為

式中：E0，D hε0，Dεε0，αh0，αε0分別由式（14）中φ＝1（無(wú)幾何修正）定義的相關(guān)矩陣。

相應(yīng)的零階翹曲函數(shù)為

將式（16）代入式（15）得到漸近修正到O（1）階的能量泛函為

這與如下形式的熱彈性古典層合模型相一致，

式中：A，NT分別為二維剛度矩陣和溫度產(chǎn)生的應(yīng)力合力，其計(jì)算式為

2．2 一階近似

零階近似可用于分析薄殼的全局性能和面內(nèi)分量，更高階近似可用來(lái)分析對(duì)中厚殼失效十分重要的面外應(yīng)力和應(yīng)變（σi3，Γi3），從而構(gòu)建更精確殼體模型。

首先，將能量泛函漸近修正到O（h／R）階以考慮初始曲率效應(yīng)，O（h／R）階翹曲因?qū)δ芰繜o(wú)貢獻(xiàn)，可不必計(jì)算。得到

式中：

其中帶星矩陣由φ－1代替式（14）的φ定義。

其次，將能量漸近修正到O（h2／l2）以考慮橫向剪切變形。為此，將零階翹曲函數(shù)攝動(dòng)為

將式（23）代回式（17），可得到一階近似的總能量泛函主導(dǎo)項(xiàng)為

式中：

相應(yīng)的一階翹曲可求解為

最后得到漸近修正到O（h2／l2），O（h／R）的總能量泛函為

式中：

2．3 近似能量轉(zhuǎn)換及三維場(chǎng)重構(gòu)

盡管式（27）能量泛函漸近修正到O（h2／l2），O（h／R）階，但因含有廣義應(yīng)變的導(dǎo)數(shù)ε；α，難以直接應(yīng)用。為得到實(shí)用的能量范函，可將式（27）轉(zhuǎn)換為工程中常用的Reissner－Mindlin模型形式。

在Reissner－Mindlin模型中有兩個(gè)附加橫向剪切應(yīng)變Reissner－Mindlin模型應(yīng)變量R與ε的關(guān)系可表示為

式中：

將式 （29）代 入 式 （27），可得到 由 Reissner－Mindlin應(yīng)變量表示的能量泛函為

而實(shí)際應(yīng)用的廣義Reissner－Mindlin模型形式為

為得到與式（32）等效的 Reissner－Mindlin模型，可通過(guò)以下2個(gè)彎矩平衡方程消除所有二維應(yīng)變量的偏導(dǎo)數(shù)

由式（33），可將式（31）改寫為

式中：

若對(duì)任何R，U＊都趨于零，則可得到漸近修正Reissner－Mindlin殼體模型。對(duì)于一般各向異性殼體，該項(xiàng)往往并不為零，可通過(guò)最小二乘法等優(yōu)化技術(shù)最小化U＊，等效Reissner－Mindlin模型的精確性取決于U＊趨近于零的程度。

由于降維模型的可靠性最終取決于其對(duì)三維場(chǎng)預(yù)測(cè)的精確度，因此還需提供重構(gòu)關(guān)系以完善降維模型。由式（3），重構(gòu)的三維應(yīng)變場(chǎng)可表示為

三維應(yīng)力場(chǎng)可使用材料本構(gòu)關(guān)系得到，

3 算例

柱形彎曲問(wèn)題已作為評(píng)估新提出的3D／2D分析模型精確性的基準(zhǔn)。本章將所述理論和方法編制成變分漸近板／殼分析程序VAPAS，對(duì)圖2所示4層簡(jiǎn)支復(fù)合材料層合殼體在熱／載荷下的柱形彎曲問(wèn)題進(jìn)行分析。

圖2 層合殼體結(jié)構(gòu)示意圖

3．1 模型參數(shù)

殼體 各 層 傾 角 為 ［90°／0°／90°／0°］；殼 厚h＝1 mm，半徑R＝10 cm，夾角φ＝π／3；采用的坐標(biāo)系為x1∈ ［0，φ］，x2∈ ［0，∞］，x3∈ ［－h(huán)／2，h／2］；材料為石墨／環(huán)氧復(fù)合材料，材料屬性為

殼體承受的溫度變化和正弦面荷載分別為

3．2 數(shù)值分析與討論

圖3（a）～（f）繪出了重構(gòu)的沿厚度方向應(yīng)力分布，并與一階剪切變形理論（FOSDT）、古典層合理論（CLT）和精確解［11－12］進(jìn)行對(duì)比。應(yīng)變和位移的變化趨勢(shì)和精度與應(yīng)力相同，限于篇幅未在此繪出。由于應(yīng)力分量σαβ，σ33分別是正弦和余弦函數(shù)，應(yīng)力分布分別繪于x1＝π／6和x1＝π／3處。為方便比較，圖中縱橫坐標(biāo)分別正則化為

由圖3可看出：VAPAS、FOSD和CLT都能較準(zhǔn)確地預(yù)測(cè)面內(nèi)應(yīng)力分量σαβ分布，由于VAPAS能得到優(yōu)化剪切剛度矩陣G，其結(jié)果比FOSDT和CLT更精確；對(duì)于橫向應(yīng)力分量σα3，因特殊的鋪層設(shè)計(jì)（正交鋪設(shè)），CLT無(wú)法計(jì)算該值，而FOSDT因?qū)ξ灰茍?chǎng)所做先驗(yàn)性假設(shè)，誤差較大，VAPAS與精確解相一致，且重構(gòu)的應(yīng)力分布在不同層分界面處連續(xù)，這與大多數(shù)位移基殼理論完全不同。

圖3 4層簡(jiǎn)支復(fù)合材料圓柱殼體在熱／載荷下沿厚度方向的應(yīng)力分布

4 結(jié)論

1）基于變分漸近法構(gòu)建了復(fù)合材料圓柱殼熱彈性簡(jiǎn)化Reissner－Mindlin模型和重構(gòu)關(guān)系，可考慮熱／載荷單向耦合效應(yīng)，并使用4次多項(xiàng)式表示沿厚度方向的任意溫度分布，相較于假設(shè)溫度沿厚度方向線性分布（單層板理論）或沿層線性分布（層合理論）更符合實(shí)際情況。

2）在推導(dǎo)降維模型過(guò)程中，使用變分漸近法求解未知翹曲函數(shù)，不需任何動(dòng)力學(xué)假設(shè)和剪切修正因子，也不同于傳統(tǒng)的板／殼理論假設(shè)翹曲場(chǎng)的一般形式，用高階翹曲作參數(shù)來(lái)求解假設(shè)函數(shù)中的未知參數(shù)的方法。

3）通過(guò)算例驗(yàn)證：重構(gòu)的沿厚度方向應(yīng)力分量與三維精確解吻合很好，且簡(jiǎn)化模型為等效單層殼模型，計(jì)算量與一階剪切變形理論相當(dāng)，若殼體層數(shù)增加，其高效性更加顯著。

［1］Shen H S．Postbuckling of nanotube－reinforced composite cylindrical shells in thermal environments， Part II：pressure－loaded shells［J］．Composite Structures，2011，93（10）：2496－2503．

［2］Li Z M，Yang D Q．Thermal postbucking analysis of 3D braided composite cylindrical shells［J］．Journal of Mechanics，2010，26（2）：113－122．

［3］Ng T Y，Lam K Y，Reddy J N．Dynamic stability of cross－ply laminated composite cylindrical shells ［J］．International Journal of Mechanical Sciences，1998，40（8）：805－823．

［4］Pradyumna S，Gupta A．Nonlinear dynamic stability of laminated composite shells integrated with piezoelectric layers in thermal environment［J］．Acta Mechanica，2011，218（3／4）：295－308．

［5］Li Z M，Lin Z Q，Chen G L．Nonlinear buckling and postbuckling behavior of 3D braided composite cylindrical shells under external pressure loads in thermal environments［J］．Journal of Pressure Vessel Technology，2009，131（4）：258－269．

［6］Nosier A，Miri A K．Boundary－layer hygrothermal stresses in laminated，composite，circular，cylindrical shell panels［J］．Archive of Applied Mechanics，2010，80（4）：413－440．

［7］Pinhas B Y．New variational－asymptotic formulations for interlaminar stress analysis in laminated plates［J］．Journal of Applied Mathematics and Physics，1986，37（8）：305－321．

［8］Lee C Y．Dynamic variational asymptotic procedure for laminated composite shells：part I：low－frequency vibration analysis［J］．Journal of Applied Mechanics，2009，76（9）：110－122．

［9］Reddy J N．Mechanics of laminated composite plates and shells，theory and analysis［M］．2 ed．CRC Press，2004．

［10］鐘軼峰，余文斌．用變分漸近法進(jìn)行復(fù)合材料層合板仿真及三維場(chǎng)重構(gòu)［J］．復(fù)合材料學(xué)報(bào)，2010，27（4）：174－179．ZHONG Yifeng，YU Wenbin．Simulation and 3D field recovery of composite laminated plates by use of variational asymptotic method ［J］． Acta Materiae Compositae Sinica，2010，27（4）：174－179．

［11］Ren J G．Exact solutions for laminated cylindrical shells in cylindrical bending ［J］．Composite Science and Technology，1987，29（3）：169－187．

［12］Jing H S，Tzeng K G．Elasticity solution for laminated anisotropic cylindrical panels in cylindrical bending［J］．Composite Structures，1995，30（4）：307－317．