王澤溪，萬志強(qiáng)，王曉喆，楊超

(1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191；2.北京航空航天大學(xué)無人系統(tǒng)研究院，北京 100191)

復(fù)合材料層合板具有明顯的質(zhì)量優(yōu)勢，已經(jīng)被現(xiàn)代飛行器結(jié)構(gòu)廣泛使用[1]。現(xiàn)代飛行器對(duì)于結(jié)構(gòu)減重具有迫切需求，復(fù)合材料具有比強(qiáng)度高、比剛度大、可設(shè)計(jì)性強(qiáng)等優(yōu)勢[2]，在新研飛機(jī)上的應(yīng)用比重也越來越高，如大型民用飛機(jī)Boeing787，其復(fù)合材料結(jié)構(gòu)質(zhì)量占到了全機(jī)結(jié)構(gòu)質(zhì)量的50%[3]。傳統(tǒng)直線纖維復(fù)合材料層合板承受面內(nèi)壓縮和剪切載荷時(shí)易發(fā)生失穩(wěn)破壞，在使用過程中復(fù)合材料層合板的性能優(yōu)勢得不到充分發(fā)揮[4]。纖維鋪放技術(shù)的傳統(tǒng)設(shè)計(jì)方法采用直線纖維制備材料層合板的鋪放技術(shù)，受限于絲束鋪設(shè)技術(shù)和設(shè)備的制約及優(yōu)化時(shí)的計(jì)算消耗，大多采用0°、±45°、90°鋪層方向，在承受面內(nèi)壓縮和剪切時(shí)，其穩(wěn)定性無法滿足要求。隨著現(xiàn)代自動(dòng)鋪放制造技術(shù)的不斷發(fā)展和成熟，為了滿足更高的性能需求，采取纖維曲線鋪放制備變剛度復(fù)合材料成了一個(gè)新的重要研究方向。

復(fù)合材料纖維曲線鋪放技術(shù)的迅速發(fā)展，驅(qū)動(dòng)了現(xiàn)代飛行器結(jié)構(gòu)采用新鋪設(shè)形式，實(shí)現(xiàn)更好的超輕、高強(qiáng)、高效結(jié)構(gòu)性能。曲線纖維復(fù)合材料層合板的纖維鋪設(shè)角可以連續(xù)變化[4]，剪裁設(shè)計(jì)空間大[5]、減重優(yōu)勢明顯[6]。不斷變化的纖維方向角使得每個(gè)鋪層的剛度在不同位置各不相同，設(shè)計(jì)者可由此調(diào)整鋪層的內(nèi)在剛度分布，提高結(jié)構(gòu)效率[7]。相比于傳統(tǒng)直線纖維復(fù)合材料，曲線纖維復(fù)合材料具有如下突出優(yōu)勢：板厚度不變的情況下可以按需求改變其剛度分布[4-5]；僅通過改變纖維方向角就可以實(shí)現(xiàn)層合板在厚度上按設(shè)計(jì)需求連續(xù)變化[7-8]。為充分發(fā)揮復(fù)合材料纖維的承載性能，改變傳統(tǒng)復(fù)合材料鋪層結(jié)構(gòu)和方式，采用纖維曲線鋪放技術(shù)制備變剛度復(fù)合材料層合板，在降低制造成本和提高復(fù)合材料性能方面顯示出了極為突出的優(yōu)越性和極大的潛力[9]，將會(huì)成為未來高性能纖維復(fù)合材料發(fā)展的主流趨勢之一，其突出的設(shè)計(jì)性和減重特點(diǎn)，也將使曲線纖維復(fù)合材料結(jié)構(gòu)制備這一技術(shù)在飛行器結(jié)構(gòu)的氣動(dòng)彈性剪裁設(shè)計(jì)中發(fā)揮更重要的作用[10]。

壁板局部屈曲問題是機(jī)翼結(jié)構(gòu)設(shè)計(jì)中需要重點(diǎn)考慮的關(guān)鍵約束，而各向異性板的屈曲問題一直是學(xué)界的研究重點(diǎn)。Nardo[11]、Thielemann[12]和Green[13]較早對(duì)膠合板的屈曲問題進(jìn)行了研究，并提出了傅里葉級(jí)數(shù)形式的封閉解。Ashton和Waddoups[14]將瑞利-里茲法應(yīng)用于分析各向異性板的穩(wěn)定性和動(dòng)力學(xué)響應(yīng)。瑞利-里茲法是求解正交各向異性材料即復(fù)合材料層合板時(shí)的常用半解析求解方法，許多學(xué)者應(yīng)用基于正交多項(xiàng)式的瑞利-里茲法進(jìn)行結(jié)構(gòu)分析，發(fā)現(xiàn)正交多項(xiàng)式系較傅里葉級(jí)數(shù)或梁的模態(tài)函數(shù)有更好的收斂性。Loja 等[15]采用正交多項(xiàng)式，基于瑞利-里茲法分析層合板屈曲、自由振動(dòng)和動(dòng)力穩(wěn)定性等問題。Wu 等[16-17]利用瑞利-里茲法對(duì)曲線纖維層合板的穩(wěn)定性問題進(jìn)行了研究，著重考慮了邊界條件的引入和消去剛度的微分項(xiàng)，解決了變剛度層合板中Levy 型解析解難以求出的問題；還提出了一種改進(jìn)的瑞利-里茲方法用于均勻軸壓下的圓孔變角度纖維板的屈曲分析，證明了利用變剛度提高復(fù)合材料結(jié)構(gòu)的抗屈曲性能是可行的。

國內(nèi)也有多位學(xué)者和機(jī)構(gòu)基于現(xiàn)有商業(yè)軟件對(duì)纖維曲線鋪設(shè)技術(shù)開展了相關(guān)研究。如秦永利等[18]對(duì)纖維曲線鋪放的變剛度復(fù)合材料層合板的研究進(jìn)展進(jìn)行了具體的介紹，馬永前等[19]用ABAQUS有限元軟件對(duì)纖維曲線鋪放的復(fù)合材料層合板進(jìn)行了建模計(jì)算，驗(yàn)證其面內(nèi)受力情況下，屈曲荷載顯著提高，幅度達(dá)14%左右。杜宇等[20]通過復(fù)合材料層合板鋪層角度的設(shè)計(jì)，使曲線纖維層合板的極限載荷和強(qiáng)度載荷較直線時(shí)提高30%。牛雪娟等[21]通過對(duì)變剛度鋪放路徑進(jìn)行數(shù)學(xué)建模得到變剛度復(fù)合材料層合板的拓?fù)錁?gòu)型，利用有限元分析軟件GENESIS 實(shí)現(xiàn)層合板拉伸分析。衛(wèi)宇璇等[22]基于有限元軟件對(duì)曲線纖維層合板進(jìn)行了仿真模擬，研究了曲線纖維軌跡對(duì)復(fù)合材料層合板力學(xué)性能的影響。

綜合來看，有限元法能夠較好處理復(fù)雜的結(jié)構(gòu)形式和邊界條件，使研究對(duì)象從層合板拓展到復(fù)雜翼面結(jié)構(gòu)，但有限元法的求解耗費(fèi)大，單純采用有限元模型會(huì)大幅增加結(jié)構(gòu)優(yōu)化的難度和計(jì)算耗費(fèi)；解析法不能準(zhǔn)確考慮氣動(dòng)彈性變形下的任意邊界條件，無法適用于氣動(dòng)彈性優(yōu)化。

為了解決有限元方法建模困難、求解速度慢，而解析法無法準(zhǔn)確描述復(fù)雜邊界條件的難題，本文通過在總勢能方程中引入拉格朗日乘子的方式將復(fù)雜邊界條件引入到壁板控制方程中，并通過Airy應(yīng)力函數(shù)描述待求解的壁板內(nèi)應(yīng)力分布，以此建立了能考慮復(fù)雜邊界條件的曲線纖維壁板力學(xué)模型；通過瑞利-里茲法進(jìn)行解析求解，并采用勒讓德多項(xiàng)式作為形函數(shù)以加快收斂速度、提高計(jì)算精度；采用牛頓-辛普森迭代方法進(jìn)行后屈曲階段的位移和應(yīng)力求解，并使用自適應(yīng)迭代步長來加快收斂速度。通過與有限元計(jì)算結(jié)果進(jìn)行對(duì)比來驗(yàn)證該快速求解方法的準(zhǔn)確性和精度，并比較二者在相同計(jì)算條件下的求解效率。最后基于該快速求解方法總結(jié)曲線纖維壁板在不同邊界條件和設(shè)計(jì)構(gòu)型下的屈曲/后屈曲力學(xué)特性。

1 基本理論與力學(xué)模型建立

1.1 曲線纖維描述方法

本文的參考纖維路徑是通過沿特定方向和壁板的特定區(qū)域內(nèi)線性改變纖維方向來定義的，如圖1 所示；其余各纖維通過將該纖維沿一定方向平移得到。本文使用的命名法 ＆lt;θ1|θ2|θ3＆gt;表示單層的纖維變化，表示在壁板中設(shè)有3 個(gè)控制點(diǎn)，即在邊界處有2 個(gè)控制點(diǎn)(點(diǎn)1 和點(diǎn)3)，在中心處有1 個(gè)控制點(diǎn)(點(diǎn)2)。纖維角度隨位置變化形式可由式（1）描述。

圖1 曲線纖維壁板纖維路徑描述方法Fig.1 Description method of VAT fiber path

式中： θ(x) 為 坐標(biāo)為x處 的纖維方向角；a為平板長度； θ1、 θ2和 θ3分別為3 個(gè)控制點(diǎn)處的纖維方向角。

本文曲線纖維設(shè)計(jì)采用2 種典型方法進(jìn)行，即自動(dòng)纖維 鋪放技術(shù)(automatic fiberplacement,AFP)和連續(xù)纖維束剪切（continuous t ow shearing，CTS)法[23]。其中，AFP 方法不會(huì)改變纖維厚度，而采用CTS 方法進(jìn)行設(shè)計(jì)時(shí)，由于該生產(chǎn)方式使得纖維束在層合板中以連續(xù)受力擠壓的方式改變方向時(shí)保持單位長度的體積(即單位長度中能鋪放纖維束的質(zhì)量)恒定，因而纖維束的厚度增加、寬度減小，如圖2 所示。采用CTS 方法進(jìn)行曲線纖維設(shè)計(jì)時(shí)，纖維束的厚度分布可以采用恒定體積假設(shè)計(jì)算[23]：

圖2 絲束被剪切轉(zhuǎn)向?qū)е碌暮穸茸兓痆23]Fig.2 Tow thickness variation by shearing [23]

式中：t為x處纖維束厚度；t0為參考點(diǎn)處的纖維束厚度； θ0為參考點(diǎn)處纖維方向角。

1.2 小變形分析基本理論

在前屈曲(pre-buckling)階段，即載荷達(dá)到臨界屈曲載荷之前，薄板處于小變形情況下，此時(shí)位移和應(yīng)變是線性的；其中性面的橫向正應(yīng)力相比其他應(yīng)力小得多，因而往往將其忽略。應(yīng)變-位移關(guān)系也為線性，可表示為

式中：u、v、w分別為板各點(diǎn)在x、y、z這3 個(gè)方向的位移； εx、 εy和 εz為 對(duì)應(yīng)方向的正應(yīng)變； γ為切應(yīng)變。式(3)中的應(yīng)變分量在相容條件下變形可以得到：

式中：上標(biāo)0 為在板的中性面處。

引入Airy 應(yīng)力函數(shù) Φ，其定義為

式中：Nx、Ny為壁板在x和y方向的合應(yīng)力，即各層應(yīng)力沿厚度積分后的結(jié)果；Nxy為切應(yīng)力沿厚度積分結(jié)果。為簡潔起見，本文部分式中的偏微分記號(hào)會(huì)以帶逗號(hào)角標(biāo)的形式出現(xiàn)，即w,x=?w/?x，可以將Airy 應(yīng)力函數(shù)簡寫為

為進(jìn)行屈曲分析，本文有如下假設(shè)：在屈曲出現(xiàn)之前，隨著載荷的增加，薄板處于二維應(yīng)力狀態(tài)下，在本文記為[]，其中系數(shù) λ是隨著載荷增加而成比例單調(diào)增加的系數(shù)；應(yīng)變增量與位移增量仍為線性關(guān)系，且應(yīng)力增量依然遵循線性應(yīng)力-應(yīng)變關(guān)系，則面內(nèi)的合應(yīng)力增量與位移仍為如式(3)描述的線性關(guān)系。將式(3)～式(6)代入虛功原理表達(dá)式中，最終可得到各向同性薄板屈曲問題的最小勢能原理表達(dá)式為

其中

式中： Πb為 勢能的泛函；D為板的彎曲剛度；積分符號(hào)中的Sm為 板的積分區(qū)域； μ為泊松比。

1.3 曲線纖維層合板小變形分析方法

對(duì)直線纖維層合板，根據(jù)經(jīng)典層合板理論得到本構(gòu)方程[24]：

式中：A、B和D分別為復(fù)合材料板的面內(nèi)、耦合和彎 曲 剛度矩陣； ε0為 中性面的應(yīng)變； κ為 中 性 面 曲率；N和M分別為面內(nèi)合應(yīng)力和合彎矩。

而對(duì)于曲線纖維層合板，其A、B、D矩陣中各項(xiàng)不再是常數(shù)，而會(huì)隨位置改變，此時(shí)式(8)可以寫成如下形式[25]：

為便于后續(xù)計(jì)算，將式(9)改寫成部分求逆形式[25]：

對(duì)于本文所研究的對(duì)稱均衡鋪層，彎曲-面內(nèi)耦合矩陣B=0，從而B*=0，D*=D。

根據(jù)式(10)的本構(gòu)方程，中性面的應(yīng)變和合應(yīng)力N之間的關(guān)系可以由如下方程得出：

式中：aij為 面內(nèi)剛度矩陣A中各項(xiàng)對(duì)應(yīng)數(shù)值。將相容方程式(4)和平衡方程式(11)代入應(yīng)變能方程，并引入拉格朗日乘子描述的邊界條件，可以得到曲線纖維壁板在小變形情況下由于面內(nèi)拉伸/壓縮所引起的總泛函方程：

式中：C1為 載荷邊界；C2為位移約束邊界；Mv0為中性面處外力矩；V z0為中性面處集中力；u0和v0分別為位移約束邊界沿x和y方向的位移分布。式(12)中受變分影響的函數(shù)為待求的應(yīng)力函數(shù) Φ和邊界應(yīng)力，可通過瑞利-里茲法進(jìn)行求解。同樣可以得到曲線纖維薄板在小變形情況下由于彎曲引起的總泛函方程：

式中： λ為屈曲因子；等號(hào)右邊第1 個(gè)積分項(xiàng)為板的彎曲所具有的應(yīng)變能；等號(hào)右邊第2 個(gè)積分項(xiàng)為相容方程和位移方程的耦合關(guān)系。式(13)中受變分影響的函數(shù)為位移函數(shù)w。

1.4 曲線纖維層合板大變形分析方法

在進(jìn)行相應(yīng)的分析和公式推導(dǎo)之前，本文對(duì)于大變形的定義如下：

1)變形足夠大(即變形是板厚度的數(shù)倍)，但仍足夠小到適用簡單的形式分析板的曲度；

2)在彎曲時(shí)，中性面產(chǎn)生了應(yīng)變。

此時(shí)任意一點(diǎn)處的應(yīng)力-應(yīng)變關(guān)系仍為線性，但位移-應(yīng)變關(guān)系不再是線性，可通過微分形式的馮卡門大變形方程(von Kármánlargedef lection equations)作為控制方程給出[24]。忽略高階項(xiàng)，板的中性面應(yīng)變和中性面位移間的關(guān)系可以通過下式定義：

式中：e為非線性狀態(tài)下的應(yīng)變，以區(qū)分與線性狀態(tài)下的應(yīng)變?chǔ)拧?/p>

需要注意的是，對(duì)于大變形下的屈曲問題，其非線性部分是由式(14)中非線性項(xiàng)的引入導(dǎo)致的，該部分非線性項(xiàng)導(dǎo)致了面內(nèi)變形與彎曲變形的耦合，而這與材料屬性(各向同性或各向異性；曲線纖維或直線纖維)均無關(guān)。但曲線纖維層合板由于其剛度特性隨位置變化，會(huì)使得分析后的應(yīng)力分布呈現(xiàn)更加非線性的形式。

當(dāng)板處于大變形情況時(shí)，中性面產(chǎn)生了較大傾斜，此時(shí)的面內(nèi)合應(yīng)力Nx、Ny和Nxy會(huì)影響到z方向的平衡方程，曲線纖維壁板在大變形狀態(tài)下的非線性平衡方程為[16,25]

同樣地，可以得到曲線纖維板通過Airy 應(yīng)力函數(shù)表示的非線性相容方程[16,25]：

進(jìn)而可以建立曲線纖維壁板在大變形條件下，含有邊界條件的總勢能變分方程[16,25]：

式中：等號(hào)右邊前2 項(xiàng)積分式為薄板的面內(nèi)拉伸/壓縮和彎曲引起的應(yīng)變能；等號(hào)右邊第3 項(xiàng)積分式為以應(yīng)力函數(shù) Φ的形式表示了相容方程和平衡方程之間的耦合關(guān)系；等號(hào)右邊最后2 項(xiàng)積分式為一般邊界條件；下標(biāo)c1和c2分別為定義了應(yīng)力和位移約束的邊界。積分項(xiàng)s和v分別為沿給定邊界的切線和法線方向。通過引入拉格朗日乘子，非線性的應(yīng)變-位移關(guān)系能夠被包含在勢能方程式(17)中，從而不需要單獨(dú)求解相容方程和平衡方程，這意味著僅需求解靜止?fàn)顟B(tài)下的方程式(17)就能同時(shí)滿足非線性平衡方程、相容方程、給定的彎矩及橫向剪切合應(yīng)力(面外)邊界條件和位移(面內(nèi))邊界條件。

2 瑞利-里茲求解方法

2.1 分析坐標(biāo)系與形函數(shù)選擇

本文將對(duì)y方向無約束(見圖3(a))和y方向有位移約束(見圖3(b))2 種不同的邊界條件展開研究；如未加特殊說明，則表示纖維角度沿x方向改變，載荷也沿x方向施加于薄板上。為簡單起見，將分析坐標(biāo)系進(jìn)行歸一化：

圖3 曲線纖維薄板載荷、位移及坐標(biāo)系示意圖Fig.3 Load, displacement and coordinate of VAT fiber plate

式中：a為平板長度；b為平板寬度。

本文采用瑞利-里茲模型進(jìn)行屈曲、后屈曲分析，將位移函數(shù)w和Airy 應(yīng)力函數(shù) Φ假設(shè)為如下的級(jí)數(shù)形式[25]：

式中：Xm、Yn、Xp、Yq為滿足給定邊界條件的形函數(shù)；Φ0(ξ,η)為沿邊界處的未知法向應(yīng)力分布(Nx0,Ny0)?？梢员硎境扇缦滦问揭苑蠈?shí)際應(yīng)力分布：

式中：cl和dl分別為應(yīng)力函數(shù)對(duì)于給定邊界條件的未知應(yīng)力場的各階形函數(shù)的待定系數(shù)。具體地，在薄板各方向邊界處可以得到邊界處應(yīng)力與形函數(shù)之間的對(duì)應(yīng)關(guān)系：

式中： φ (η)、 ψ (ξ)為應(yīng)力的形函數(shù)。給定的位移載荷形式需要與待求的邊界合應(yīng)力分布結(jié)合，并將其代入式(12)的邊界積分項(xiàng)中進(jìn)行未知量( Φpq、cl、dl)的求解。

對(duì)于式(19)～式(22)中形函數(shù)的選擇，需要考慮面內(nèi)邊界條件和對(duì)應(yīng)的面內(nèi)應(yīng)力狀態(tài)。對(duì)經(jīng)典形函數(shù)的相關(guān)研究表明[17]，三角函數(shù)多項(xiàng)式由于具有周期性，適用于求解具有周期性的力學(xué)問題；勒讓德函數(shù)具有非周期性，適用于非周期性問題的求解。對(duì)于本文所研究的穩(wěn)定性問題，非周期的勒讓德多項(xiàng)式作為形函數(shù)具有更好的收斂性和求解精度。為加快收斂速度，本文選擇勒讓德多項(xiàng)式作為形函數(shù)；為使其形式滿足邊界條件，需對(duì)勒讓德多項(xiàng)式乘以( 1?ξ2)或 ( 1?η2)：

其中勒讓德形函數(shù)Lr(ξ)一般形式為[17]

本文進(jìn)行的數(shù)值實(shí)驗(yàn)結(jié)果如圖4 所示，縱軸為計(jì)算值和有限元結(jié)果的比值；橫軸坐標(biāo)“x-y-z”中的數(shù)字含義如下：x、y分別對(duì)應(yīng)表征結(jié)構(gòu)剛度本身的形函數(shù)的階數(shù)，而z對(duì)應(yīng)邊界處應(yīng)力形函數(shù)階數(shù)。前10 階屈曲模態(tài)在選用5 階形函數(shù)時(shí)即可實(shí)現(xiàn)收斂；對(duì)于第1 階屈曲因子，只需表征結(jié)構(gòu)剛度本身形函數(shù)的階數(shù)(即x)達(dá)到7 階、其余形函數(shù)仍為5 階即可具有足夠的精度，第1 階屈曲因子相對(duì)誤差為1.3%。該實(shí)驗(yàn)結(jié)果證明，勒讓德多項(xiàng)式在分析板的屈曲問題時(shí)具有較好的收斂性和準(zhǔn)確性。

圖4 SSSS邊界條件下不同形函數(shù)階數(shù)屈曲因子計(jì)算結(jié)果對(duì)比Fig.4 Buckling factors comparison of different shape function orders under SSSS boundary condition

2.2 瑞利-里茲屈曲求解模型

小變形狀態(tài)下式(12)和式(13)可以獨(dú)立求解。將式(19)～式(22)代入式(12)，就可以得到以勒讓德多項(xiàng)式作為形函數(shù)表示的面內(nèi)拉伸/壓縮對(duì)應(yīng)邊值問題的總泛函；其中需要進(jìn)行求解的未知量共有3 組，分別是 Φpq、cl、dl。利用總泛函對(duì)未知量偏導(dǎo)數(shù)為0 的特性建立瑞利-里茲模型：

式中：fi分 別為 Φpq,cl,dl。

由式(25)可以得到一組線性代數(shù)方程，并可以表示為式(27)的張量形式，式(26)分別由式(12)對(duì)Φpq、cl和dl求偏導(dǎo)數(shù)得來：

式中：Kmc等項(xiàng)為板在小變形狀態(tài)下的變分剛度矩陣，矩陣下標(biāo)中的m、c、d 分別表示面內(nèi)、x方向的邊界(即邊界a)和y方向的邊界(即邊界b)； Φ、c和d分別為未知量 Φpq、cl和dl的列向量形式；F x、F y分別為x方向和y方向的邊界條件。式(27)可簡寫為式(28)，這與經(jīng)典理論中的線性靜力學(xué)問題具有相同表達(dá)形式：

對(duì)于圖3（a）中y方向邊界無約束時(shí)的邊界條件(SSFF)，y方向無應(yīng)力，直接得到dl=0，方程組退化為

由式(12)可得，對(duì)于坐標(biāo)歸一化后的薄板有：

注意到在本文所求解的壁板穩(wěn)定性問題中，邊界位移的分布形式u0(s)已知，且獨(dú)立于應(yīng)力分布，因此可得

由式(5)和式(22)可得

方程式(27)中并不包含位移函數(shù)w，僅能獲取邊界處的載荷和板的應(yīng)力分布形式 Φ。為了獲得板的屈曲因子及對(duì)應(yīng)的屈曲模態(tài)信息，需要通過求解彎曲對(duì)應(yīng)的總泛函方程式(13)獲得。采用與式(25)相似的方法，可以得到由待求位移函數(shù)的系數(shù)W表示的給定合應(yīng)力狀態(tài)下的屈曲問題，并可改寫成求解矩陣特征值的形式如下：

式中：矩陣下標(biāo)中的b 表示彎曲； λ為對(duì)應(yīng)剛度矩陣的特征值，按從小到大排序后即為各階屈曲模態(tài)對(duì)應(yīng)的屈曲因子；W為各階屈曲模態(tài)的特征向量(即振型)。給定載荷F與屈曲臨界載荷Fcr的關(guān)系為

2.3 大變形瑞利-里茲求解方法

與屈曲求解方法相似，將形函數(shù)表達(dá)式(19)～式(23)代入大變形總泛函方程式(17)，就可以得到以勒讓德多項(xiàng)式作為形函數(shù)表示的對(duì)應(yīng)邊值問題的總泛函；其中需要同時(shí)進(jìn)行求解的未知量共有4 組，分別是 Φpq、cl、dl和Wmn。同樣利用總泛函對(duì)未知量偏導(dǎo)數(shù)為零的特性建立瑞利-里茲模型：

式中：fi分 別為 Φpq,cl,dl,Wmn。

由式(35)可以得到一組非線性代數(shù)方程，并可以表示為如式(36)的形式，分別由式(17)對(duì) ?pq、cl、dl和Wmn求偏導(dǎo)數(shù)得來：

式中：各矩陣的記號(hào)含義與式(27)和式(33)中相同；W為式(17)中形函數(shù)系數(shù)Wmn的列向量形式。實(shí)際上，所有涉及到2 個(gè)未知量的耦合矩陣(如Kbm、Kbc、Kbd)實(shí)際為3 階張量，為了保持表達(dá)統(tǒng)一及便于理解，本文仍寫作矩陣的形式，而在其具體表達(dá)形式中寫成三維數(shù)組的形式。在計(jì)算時(shí)，涉及到此類數(shù)組的乘法運(yùn)算時(shí)僅對(duì)后二維所構(gòu)成的n個(gè)矩陣依次進(jìn)行矩陣乘法，并將結(jié)果對(duì)應(yīng)排列，則計(jì)算結(jié)果仍為列向量形式，具體過程不再贅述。

式(36)實(shí)際為式(27)和式(33)及額外的非線性耦合項(xiàng)組合得到，這些耦合項(xiàng)正來源于大變形狀態(tài)下面內(nèi)合應(yīng)力對(duì)于中性面平衡方程的影響，即面內(nèi)-彎曲耦合關(guān)系；該非線性方程無顯式求解方法。在求解時(shí)，首先通過求解式(27)或式(29)得到前屈曲狀態(tài)時(shí)的合應(yīng)力分布形式(即變量 Φ、c和d)，并代入式(33)得到板的臨界屈曲載荷Ncr和臨界狀態(tài)的位移函數(shù)(即W)；之后進(jìn)行后屈曲求解，本文通過牛頓-辛普森迭代方法來求解非線性代數(shù)方程式(36)并獲得后屈曲的平衡路徑，即加載步長載荷與位移的關(guān)系。下面給出具體的迭代求解格式推導(dǎo)。非線性方程式(36)中的4 個(gè)方程可以寫為

實(shí)際計(jì)算中，僅使用經(jīng)典的牛頓-辛普森迭代法可能會(huì)在求解部分算例時(shí)出現(xiàn)數(shù)值不穩(wěn)定，因此，需加入阻尼因子 α以增強(qiáng)穩(wěn)定性：

式中： α一般取(0,1)之間的數(shù)值，越接近0 則阻尼越大， α=1則表示無阻尼。

非線性矩陣方程F(X) 對(duì)自變量X的導(dǎo)數(shù)即方程組(37)的雅可比矩陣J，其形式為

式中：3 階張量的各耦合矩陣與向量相乘后退化為2 階張量，即矩陣；進(jìn)行此類運(yùn)算時(shí)僅對(duì)三維數(shù)組中對(duì)應(yīng)的某一維展開。

在進(jìn)行求解時(shí)，先將給定的載荷分解成若干級(jí)微小載荷逐次加載，并在每個(gè)加載步中采用迭代法求解得到位移函數(shù)w和應(yīng)力函數(shù) Φ的待定系數(shù)，直至2 次迭代的殘差 ε小于給定值(本文取10?4)后認(rèn)為收斂并計(jì)算下一個(gè)加載步。該屈曲/后屈曲求解流程如圖5 所示，其中殘差 ε的計(jì)算形式為

圖5 屈曲/后屈曲求解流程Fig.5 Analysis flow of buckling and post-buckling

3 有限元模型與求解策略

3.1 有限元求解模型

為驗(yàn)證本節(jié)所建立的曲線纖維壁板屈曲分析方法的正確性及其在機(jī)翼蒙皮局部壁板屈曲分析中的適用性，選取某大型商用客機(jī)機(jī)翼上翼面的部分蒙皮壁板作為研究對(duì)象，如圖6(a)中紅色部分所示。壁板在面內(nèi)主要承受沿展向的壓縮載荷，如圖6(b)中的箭頭所示，該壁板四周均為剛度較大的翼肋或翼梁，在機(jī)翼受到氣動(dòng)載荷而產(chǎn)生變形時(shí)，壁板上下邊界受到翼肋沿展向的擠壓而發(fā)生變形，而左右邊界向外延展的趨勢受到翼梁的限制，此時(shí)可認(rèn)為該壁板具有位移約束的邊界條件。且該壁板沿展向曲率為0，沿弦向的曲率很小，對(duì)其受到沿展向壓縮時(shí)的穩(wěn)定性影響有限，因此，可簡化成矩形平板進(jìn)行研究，其邊界條件和分析坐標(biāo)系如圖6(b)所示，瑞利-里茲法求解通過MATLAB 自編程序?qū)崿F(xiàn)。采用Patran 建立曲線纖維復(fù)合材料矩形層合板的有限元模型，如圖6(c)所示，其中紅色箭頭表示對(duì)應(yīng)單元內(nèi)纖維方向；并采用Nastran 對(duì)其進(jìn)行了屈曲分析。壁板的幾何尺寸為 450 mm×750 mm，材料常數(shù)及鋪層厚度如表1 所示。本文除分析圖3 中所示的側(cè)邊自由(SSFF)和側(cè)邊施加位移約束(SSSS)這2 種不同情況外(其中a=450 mm，b=750 mm，所給定的單側(cè)位移邊界條件 ?x=10?3mm)，還對(duì)機(jī)翼在氣動(dòng)彈性變形下的復(fù)雜邊界條件進(jìn)行屈曲/后屈曲分析，并對(duì)比3 種不同邊界條件對(duì)屈曲特性的影響。

表1 壁板采用材料的材料屬性Table 1 M aterial properties of com posite used in panel

圖6 局部屈曲分析所選取的壁板示意圖Fig.6 Panel description of local buckling analysis

3.2 機(jī)翼彈性變形下壁板面內(nèi)位移提取方法

在實(shí)際飛行狀態(tài)中，機(jī)翼由于氣動(dòng)力作用而產(chǎn)生彈性變形。從宏觀角度觀察，該彈性變形表現(xiàn)為向上彎曲的撓度和各剖面的扭轉(zhuǎn)；對(duì)局部壁板而言，該壁板的局部分析坐標(biāo)系發(fā)生了位移和旋轉(zhuǎn)，且該位移遠(yuǎn)遠(yuǎn)大于其邊界處在面內(nèi)的位移分量，如圖7 所示。其中，藍(lán)色的點(diǎn)為機(jī)翼變形前該壁板邊界上的點(diǎn)，為方便分析將壁板平移其形心與坐標(biāo)原點(diǎn)重合；紅色的點(diǎn)為機(jī)翼變形后該壁板邊界上的點(diǎn)。將變形后壁板的形心平移到坐標(biāo)原點(diǎn)便可以清晰地觀察到壁板局部坐標(biāo)系即壁板本身的旋轉(zhuǎn)(見圖7(b))。需計(jì)算變形后局部坐標(biāo)系在原分析坐標(biāo)系下的表示形式，進(jìn)而得出變形后的壁板邊界各點(diǎn)在對(duì)應(yīng)的局部坐標(biāo)系中的坐標(biāo)。

圖7 機(jī)翼氣動(dòng)彈性變形前后上翼面壁板邊界結(jié)點(diǎn)位置對(duì)比示意圖Fig.7 Comparation of node location on edges of focusing panel before and after aeroelastic deformation

坐標(biāo)系的旋轉(zhuǎn)和平移采用如下方式：

式中： α、 β 和 γ分別為局部坐標(biāo)系3 個(gè)坐標(biāo)軸對(duì)應(yīng)的單位向量在全局分析坐標(biāo)系中的表達(dá)式(列向量形式[)；]為邊界處各點(diǎn)在全局坐標(biāo)系中的坐標(biāo)；為對(duì)應(yīng)點(diǎn)在局部分析坐標(biāo)系中的坐標(biāo)。

首先需要獲得氣動(dòng)彈性分析后各結(jié)點(diǎn)的位移，進(jìn)而可得到邊界各點(diǎn)變形后的坐標(biāo)；分別對(duì)變形前和變形后2 組結(jié)點(diǎn)采用式(41)的方法轉(zhuǎn)換為局部坐標(biāo)系中的坐標(biāo)，此時(shí)便可將變形前后的壁板置入同一坐標(biāo)系中進(jìn)行分析。邊界各點(diǎn)的面內(nèi)位移也能隨之獲得，如圖8 所示；注意此位移僅為示意而非真實(shí)位移，由于邊界面內(nèi)位移相比板的尺寸小2～3 個(gè)數(shù)量級(jí)，圖8 中繪制的紅色位移示意為各點(diǎn)實(shí)際位移放大1 000 倍后的結(jié)果。在本文的分析中，主要考慮由于面內(nèi)壓縮/拉伸位移對(duì)穩(wěn)定性的影響，忽略對(duì)屈曲影響相對(duì)較小的剪切位移。最終提取得到的該壁板受壓縮邊界和橫側(cè)邊界的位移邊界條件隨位置變化情況如圖9 所示。

圖8 轉(zhuǎn)換為局部坐標(biāo)系后各結(jié)點(diǎn)位移示意圖Fig.8 Diagram of displacement of each node after conversion to a local coordinate system

圖9 受壓和橫側(cè)邊界的位移分布示意圖Fig.9 Displacement distribution of loading edges and transverse edges

3.3 后屈曲特性有限元分析策略

本文的非線性有限元分析使用MSC.Nastran 的SOL400 模塊完成。在進(jìn)行幾何非線性分析時(shí)，對(duì)于不同構(gòu)型，采用固定步長加載的方式難以獲取各構(gòu)型下的準(zhǔn)確臨界屈曲載荷；為了更準(zhǔn)確地獲取不同曲線纖維路徑壁板的力學(xué)特性，并盡可能減小有限元軟件計(jì)算耗費(fèi)，本文采用了自適應(yīng)的變步長分段加載策略。在進(jìn)行非線性分析前，首先使用線性求解方法得出當(dāng)前構(gòu)型的臨界屈曲位移 ?xcr及1 階屈曲模態(tài)；而后基于此進(jìn)行幾何非線性分析，在不同階段采用不同加載步長，如圖10(a)所示，其中橫坐標(biāo)為單側(cè)位移，縱坐標(biāo)為壁板加載端約束力。

通過計(jì)算相鄰2 個(gè)加載點(diǎn)的斜率即可得到當(dāng)前狀態(tài)下結(jié)構(gòu)的等效剛度，如圖10(b)所示(壁板鋪層參數(shù)為進(jìn)而準(zhǔn)確獲得結(jié)構(gòu)發(fā)生屈曲的臨界點(diǎn)。幾何非線性分析表明，對(duì)于本文的碳纖維復(fù)合材料曲線纖維壁板，采用線性方法得到的屈曲臨界變形壁板端面合力與非線性方法獲得的實(shí)際合力一致，誤差可忽略。而屈曲實(shí)際發(fā)生時(shí)所對(duì)應(yīng)的位移略小于線性分析結(jié)果，大致為線性結(jié)果的97%～98%；達(dá)到屈曲載荷后，壁板面內(nèi)承載能力迅速下降。這表明，進(jìn)行飛行器結(jié)構(gòu)的壁板設(shè)計(jì)時(shí)，采用線性分析方法及線性臨界屈曲載荷作為壁板穩(wěn)定性約束是合理的，但超出該臨界屈曲載荷后線性方法不能得到準(zhǔn)確結(jié)果。

圖10 鋪層曲線纖維壁板(SSFF)幾何非線性屈曲分析結(jié)果Fig.10 Geometry nonlinear buckling results of VAT panel under SSFF boundary condition

通過改變壁板x、y方向單元?jiǎng)澐置芏?，可研究有限元方法的?jì)算耗費(fèi)與網(wǎng)格規(guī)模間的關(guān)系，如圖11 所示，其對(duì)應(yīng)邊界條件為SSSS 邊界，鋪層參數(shù)采用CTS 設(shè)計(jì)，計(jì)算環(huán)境為W indows10 系統(tǒng)、i5-7200UCPU(2.50GHz)、8GB 內(nèi)存的個(gè)人計(jì)算機(jī)；其中端面載荷表示加載至2 倍屈曲臨界位移時(shí)加載邊界的合外力，Ry表示沿y方向劃分的單元數(shù)量。從中可以看出，由于曲線纖維壁板中纖維角度隨位置變化，非線性屈曲較難收斂；本文中曲線纖維的變化沿x方向，因此y方向單元密度對(duì)結(jié)果影響程度較低，當(dāng)單元?jiǎng)澐謹(jǐn)?shù)量達(dá)到20 時(shí)繼續(xù)加密對(duì)結(jié)果影響十分有限；而x方向單元直接影響壁板剛度變化情況，劃分密度不足則會(huì)引起數(shù)值波動(dòng)，單元?jiǎng)澐謹(jǐn)?shù)量在130 以上時(shí)可避免波動(dòng)，后屈曲載荷分析結(jié)果接近收斂。

圖11 端面載荷及計(jì)算耗時(shí)隨網(wǎng)格規(guī)模變化示意圖Fig.11 Edge end load and calculation cost changes as mesh density

本文建立的快速求解框架雖然也使用了迭代求解，但其是基于一個(gè)含有所有未知量的單獨(dú)方程進(jìn)行的，無需在2 個(gè)方程(即平衡方程和相容方程)之間進(jìn)行相互迭代，因此能夠快速進(jìn)行非線性分析。上述計(jì)算實(shí)驗(yàn)表明，經(jīng)典的MSC.Nastran 軟件達(dá)到數(shù)值穩(wěn)定需耗費(fèi)約200 s，其他如ABAQUS等軟件耗時(shí)大體相當(dāng)；而本文建立的快速求解方法可于10 s 左右完成非線性分析，求解速度遠(yuǎn)快于經(jīng)典有限元軟件，適用于非線性屈曲的大規(guī)模優(yōu)化。

3.4 有限元與瑞利-里茲法對(duì)比驗(yàn)證

3.4.1 第1 階模態(tài)屈曲因子對(duì)比

按3.2 節(jié)中描述的方法提取壁板受壓和橫側(cè)邊界的位移，分別采用有限元軟件Nastran 和瑞利-里茲法進(jìn)行屈曲分析。由于機(jī)翼整體的彎扭耦合效應(yīng)等影響，邊界上的位移值并非均布勻或線性變化，而呈現(xiàn)出較為復(fù)雜的分布形式。此時(shí)應(yīng)力分布相比均布位移邊界條件變化較大，應(yīng)選取前n階勒讓德多項(xiàng)式以獲得更好的收斂性。當(dāng)壁板邊界條件為均布位移約束的SSFF 或SSSS 邊界條件時(shí)，對(duì)于直線纖維壁板、AFP 固定厚度設(shè)計(jì)及CTS 變厚度設(shè)計(jì)的曲線纖維壁板，瑞利-里茲法與Nastran 第1階模態(tài)相對(duì)誤差均小于0.5%，前5 階模態(tài)相對(duì)誤差小于1%。前5 階模態(tài)振型均與Nastran 結(jié)果一致，為簡略未進(jìn)行展示。

為總結(jié)壁板屈曲特性隨纖維角度變化的規(guī)律，以壁板雙側(cè)承載邊界處控制點(diǎn)纖維角度設(shè)置為T1，壁板中心點(diǎn)處纖維角度為T2， 進(jìn)行對(duì)稱均衡纖維路徑設(shè)計(jì)，即將壁板纖維路徑設(shè)置為的對(duì)稱均衡鋪層形式；采用CTS 設(shè)計(jì)時(shí)，參考角度為邊界處的控制點(diǎn)角度T1。分別固定其中一控制點(diǎn)角度為0°，計(jì)算第1 階模態(tài)屈曲因子隨另一控制點(diǎn)角度變化趨勢，并與有限元結(jié)果進(jìn)行對(duì)比，如圖12 所示。

圖12 固定一控制點(diǎn)角度、屈曲因子隨另一控制點(diǎn)角度變化趨勢Fig.12 Buckling factor variation along with angle of one control point (with the other point fixed)

從分析結(jié)果中可以看出壁板中心處纖維角度與邊界處纖維角度對(duì)屈曲性能影響的差別。通過增加邊界處的纖維角度能夠取得更好的屈曲性能，尤其是采用CTS 變厚度設(shè)計(jì)時(shí)。在此非均布位移載荷作用下，瑞利-里茲法與有限元方法的分析精度相當(dāng)，相對(duì)誤差小于0.5%。

3.4.2 屈曲/后屈曲特性對(duì)比

通過本文建立的瑞利-里茲求解方法與有限元軟件MSC.Nastran 對(duì)壁板在不同邊界條件下的屈曲/后屈曲特性進(jìn)行對(duì)比，驗(yàn)證2 種不同分析方法的準(zhǔn)確性，及對(duì)不同邊界條件的適用性，如圖13 所示。其中，圖13(a)和圖13(b)的橫坐標(biāo)為加載邊界的單側(cè)邊界位移量；圖13(c)的橫坐標(biāo)為受壓邊界上一點(diǎn)的位移量，其余各點(diǎn)壓縮/拉伸的位移量按照?qǐng)D9中的位移分布形式相應(yīng)變化。從中可以明顯看出，對(duì)于相同位移條件，采用CTS 設(shè)計(jì)可以有效增加臨界載荷，但在后屈曲階段承載效率與AFP 設(shè)計(jì)相當(dāng)。2 種方法結(jié)果一致，驗(yàn)證了解析法和有限元方法進(jìn)行后屈曲求解的準(zhǔn)確性。

圖13 不同邊界條件曲線纖維壁板屈曲/后屈曲特性對(duì)比Fig.13 Buckling and post-buckling performance of VAT panel under different boundary conditions

4 曲線纖維壁板穩(wěn)定性特性分析

為進(jìn)一步分析曲線纖維壁板的屈曲/后屈曲特性，并更清晰地展示曲線纖維壁板穩(wěn)定性隨纖維路徑參數(shù)變化的規(guī)律，采用瑞利-里茲法對(duì)纖維路徑為的壁板進(jìn)行穩(wěn)定性分析，根據(jù)分析結(jié)果繪制相應(yīng)的載荷-位移曲線。

對(duì)邊界條件為SSSS 的曲線纖維壁板，分別固定一個(gè)控制點(diǎn)的纖維方向角，將另一控制點(diǎn)纖維方向角從0°變化至90°(變化間隔為10°)，其屈曲/后屈曲特性如圖14 所示；其中圓圈處表示該構(gòu)型的屈曲臨界位移和載荷。其載荷-位移曲線隨著一個(gè)控制點(diǎn)纖維方向角增大呈現(xiàn)出明顯規(guī)律，即臨界屈曲位移增加、斜率變小，如圖14(a)所示。

圖14 SSSS邊界條件的曲線纖維壁板屈曲/后屈曲特性Fig.14 Buckling and post-buckling performance of VAT panel under SSSS boundary condition

可以發(fā)現(xiàn)，雖然增加控制點(diǎn)的纖維方向角能夠增加其臨界屈曲位移，但對(duì)應(yīng)載荷并非持續(xù)增加；且當(dāng)固定壁板邊界處控制點(diǎn)T1=30?時(shí)，屈曲載荷隨著中心處控制點(diǎn)T2的角度增加而先減小后增加；反之，當(dāng)固定T2=30?時(shí) ，壁板屈曲載荷隨T1的角度先增加后減小。為更清晰展示這一規(guī)律，將各組控制點(diǎn)參數(shù)組合的臨界屈曲載荷在同一坐標(biāo)系中繪制，如圖15 所示。

圖15 不同邊界條件的屈曲臨界載荷Fig.15 Critical load of VAT panel under different boundary conditions

從圖15 中可以觀察到，當(dāng)采用固定厚度的AFP設(shè)計(jì)時(shí)，對(duì)于SSFF 邊界條件(見圖15(a))的壁板，其最大屈曲臨界載荷出現(xiàn)在鋪層形式為 ±45?的對(duì)稱均衡鋪層附近，與一般經(jīng)驗(yàn)相符合；對(duì)采用SSSS 邊界條件或基于氣彈變形邊界條件(見圖15(b))的壁板，采用的纖維路徑能夠使壁板獲得最大屈曲載荷，采用不同角度的直線纖維可得到的屈曲載荷則落在T1=T2的直線上，經(jīng)比較發(fā)現(xiàn)，曲線纖維最大屈曲載荷約為直線纖維最大屈曲載荷的1.28 倍。采用變厚度的CTS 設(shè)計(jì)時(shí)，3 種不同邊界條件的臨界屈曲載荷均有大幅度提升，且分布形式十分接近。

為進(jìn)一步比較屈曲前后壁板等效剛度變化，選取屈曲臨界點(diǎn)的載荷隨位移變化的斜率作為屈曲前等效剛度，選取 ?x=2?xcr即2 倍臨界位移點(diǎn)與屈曲臨界點(diǎn)連線的斜率作為屈曲后等效剛度；壁板等效剛度隨控制點(diǎn)參數(shù)變化的規(guī)律如圖16 所示。

圖16 屈曲前等效剛度分布示意圖(SSSS/AFP設(shè)計(jì))Fig.16 Distribution of equivalented stiffness before and afterbuckling(SSSS/AFP design)

壁板在3 種不同條件下屈曲前等效剛度隨控制點(diǎn)變化趨勢基本一致；采用CTS 設(shè)計(jì)時(shí)，雖能有效提高屈曲載荷，但其等效剛度并未有明顯變化。對(duì)SSSS 邊界條件和基于氣動(dòng)彈性位移的邊界條件，屈曲前后等效剛度分布形式接近，但其具體數(shù)值有大幅度下降。

為進(jìn)一步比較壁板在后屈曲時(shí)出現(xiàn)的剛度下降，將組控制點(diǎn)參數(shù)對(duì)應(yīng)的屈曲后等效剛度與屈曲前等效剛度做商，即可得到該曲線纖維路徑下屈曲后剛度下降的百分比，如圖17 所示。

圖17 SSSS邊界條件、CTS設(shè)計(jì)壁板屈曲后剛度比Fig.17 Proportion of equivalented stiffness after buckling compared with stiffness before buckling

從圖17 中可以看出，雖然纖維方向角度較小時(shí)其等效剛度較大，但在后屈曲階段下降的幅度也更大，即面內(nèi)穩(wěn)定性較差；相反，纖維角度較大時(shí)，其等效剛度較小，但在后屈曲階段的等效剛度損失更小，部分參數(shù)甚至無剛度損失?？傮w而言，曲線纖維由于具有更大的設(shè)計(jì)空間和更強(qiáng)的載荷重新分配能力，其臨界屈曲載荷更高，能夠提供更高的結(jié)構(gòu)效率和面內(nèi)穩(wěn)定性。

5 結(jié) 論

本文基于馮卡門大變形方程，發(fā)展了能考慮任意邊界條件的曲線纖維壁板屈曲/后屈曲特性快速求解方法，得出以下結(jié)論：

1）本文發(fā)展的瑞利-里茲求解方法適用于曲線纖維壁板的屈曲/后屈曲分析，線性階段和幾何非線性大變形下力學(xué)特性分析精度與有限元方法一致，且求解速度大大高于有限元軟件。

2）曲線纖維復(fù)合材料具有更大的設(shè)計(jì)空間和載荷重新分配能力，在SSSS 邊界條件或提取的氣動(dòng)彈性位移邊界條件下、幾何尺寸在m 量級(jí)、厚度在mm 量級(jí)的壁板，其曲線纖維設(shè)計(jì)得到的最大屈曲載荷約為同厚度直線纖維最大屈曲載荷的1.28 倍。

3）在后屈曲階段，壁板承載能力(即等效剛度)隨纖維角度不同而具有不同下降趨勢；纖維角度較小時(shí)該下降更為明顯。

4）本文的快速分析方法適用于在設(shè)計(jì)初步階段快速掌握曲線纖維壁板面內(nèi)穩(wěn)定性的規(guī)律或大規(guī)模優(yōu)化。