曲線纖維壁板屈曲/后屈曲建模與快速分析方法

王澤溪，萬志強，王曉喆，楊超

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

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

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

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

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

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

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

1 基本理論與力學模型建立

1.1 曲線纖維描述方法

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

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

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

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

圖2 絲束被剪切轉向導致的厚度變化[23]Fig.2 Tow thickness variation by shearing [23]

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

1.2 小變形分析基本理論

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

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

式中：上標0 為在板的中性面處。

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

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

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

其中

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

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

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

式中：A、B和D分別為復合材料板的面內、耦合和彎 曲 剛度矩陣； ε0為 中性面的應變； κ為 中 性 面 曲率；N和M分別為面內合應力和合彎矩。

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

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

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

根據(jù)式(10)的本構方程，中性面的應變和合應力N之間的關系可以由如下方程得出：

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

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

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

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

在進行相應的分析和公式推導之前，本文對于大變形的定義如下：

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

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

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

式中：e為非線性狀態(tài)下的應變，以區(qū)分與線性狀態(tài)下的應變ε。

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

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

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

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

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

2 瑞利-里茲求解方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 有限元模型與求解策略

3.1 有限元求解模型

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

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

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

3.2 機翼彈性變形下壁板面內位移提取方法

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

圖7 機翼氣動彈性變形前后上翼面壁板邊界結點位置對比示意圖Fig.7 Comparation of node location on edges of focusing panel before and after aeroelastic deformation

坐標系的旋轉和平移采用如下方式：

式中： α、 β 和 γ分別為局部坐標系3 個坐標軸對應的單位向量在全局分析坐標系中的表達式(列向量形式[)；]為邊界處各點在全局坐標系中的坐標；為對應點在局部分析坐標系中的坐標。

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

圖8 轉換為局部坐標系后各結點位移示意圖Fig.8 Diagram of displacement of each node after conversion to a local coordinate system

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

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

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

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

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

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

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

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

3.4 有限元與瑞利-里茲法對比驗證

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

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

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

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

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

3.4.2 屈曲/后屈曲特性對比

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結 論

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

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

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

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

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