原 瀟，王志瑾，夏 津，Kretov Anatolii

（1. 南京航空航天大學(xué)，飛行器先進設(shè)計技術(shù)國防重點學(xué)科實驗室，南京，210016；2. 上海機電工程研究所，上海，201109）

0 引 言

飛行器結(jié)構(gòu)中廣泛采用大長細比的薄壁梁結(jié)構(gòu)型式，如機身、機翼。對類似結(jié)構(gòu)進行準確的強度計算分析是非常復(fù)雜并且耗時耗力的事情。在計算機技術(shù)出現(xiàn)之前，要對這樣的結(jié)構(gòu)進行分析計算，經(jīng)常采用一系列假設(shè)，通過工程模型得到近似解。隨著計算技術(shù)和數(shù)值計算方法的快速發(fā)展，解析法和近似法逐漸被有限元方法取代。采用有限元方法對結(jié)構(gòu)進行應(yīng)力應(yīng)變狀態(tài)分析前，一般先對結(jié)構(gòu)進行必要的簡化，建立等效的有限元數(shù)學(xué)模型[1,2]。然而，有限元模型中包含大量的結(jié)構(gòu)信息數(shù)據(jù)，在初始設(shè)計階段，無法給出完整的結(jié)構(gòu)信息。另外，在結(jié)構(gòu)選型階段，還需要通過初始的計算分析對結(jié)構(gòu)布置提供參考依據(jù)。因此，建立簡單、高效、可靠并有足夠計算精度的數(shù)學(xué)模型仍然是非常必要。

在超單元法中，整體結(jié)構(gòu)被分為若干個超單元，每個超單元可以視為一系列有限元素的集合。超單元法通過對系統(tǒng)自由度進行減縮，從而大幅減少計算量，適用于大型復(fù)雜結(jié)構(gòu)的分析[3]。Guyan[4]將結(jié)構(gòu)自由度分為主自由度和副自由度，通過強制忽略副自由度的影響，對剛度矩陣和質(zhì)量矩陣進行縮聚，適用于靜力分析和低階模態(tài)分析；陳國平[5,6]將界面自由度和內(nèi)部自由度均進行減縮，提出了適用于結(jié)構(gòu)動特性分析的超單元法；在靜力學(xué)領(lǐng)域，Ireneusz Kreja[7]提出一種梁超單元，適用于具有對稱開剖面的薄壁梁結(jié)構(gòu)；康龍輝[8]提出了適用于恒定截面薄壁結(jié)構(gòu)靜力分析的超單元法，但只適用于無錐度的薄壁結(jié)構(gòu)。

本文針對具有錐度的薄壁梁結(jié)構(gòu)，基于橫截面內(nèi)剛性假設(shè)，將節(jié)點位移縮聚為廣義位移，對于3種簡化模型分別推導(dǎo)出超單元剛度矩陣，通過Matlab編寫相應(yīng)求解程序。通過3個算例結(jié)果與有限元結(jié)果的對比，驗證了該方法的準確性和有效性。

1 薄壁梁中的超單元法

薄壁梁結(jié)構(gòu)一般由板、桿組成，可以有任意橫剖面形狀，其典型外形特征是大長細比的圓柱或錐形殼體。結(jié)構(gòu)設(shè)計條件是滿足強度和剛度要求下結(jié)構(gòu)重量最輕。飛機機翼、機身，飛行器機體都是具有這種特征的薄壁梁結(jié)構(gòu)，它們主要由橫向構(gòu)件（翼肋或框）、縱向構(gòu)件（梁緣條、桁條）、腹板及蒙皮組成，如圖1所示。為了建立計算模型，在橫向構(gòu)件處引入與xy面平行的剖面，共（c+1）個。這樣，長度為 l的薄壁梁就由順序相連的c個單元組成，每個單元的長度為le。每個橫剖面同 z軸的交點就是廣義節(jié)點，兩個橫剖面之間的一段結(jié)構(gòu)實際上就是在兩個剖面上具有廣義位移的超單元。在對超單元進行分析時，忽略了橫向構(gòu)件，只考慮桁條、蒙皮、腹板的力學(xué)行為，并將桁條（緣條）等效為n個桿單元，蒙皮、腹板等效為m個板單元。

圖1 典型薄壁梁結(jié)構(gòu)模型Fig.1 Typical Thin-walled Beam Structure Model

續(xù)圖1

根據(jù)計算精度從低到高，假設(shè)每個橫剖面從絕對剛性到自由變形。假定每個剖面上的面內(nèi)廣義位移函數(shù)為：?1(z)，?2(z)，…，?M(z)，面外廣義位移函數(shù)為：ψ1(z)，ψ2(z)，…，ψN(z)。其中 M、N的值根據(jù)計算精度要求選取，且Mmin=3，Mmax=m；Nmin=3，Nmax=n。剖面上任意節(jié)點k的位移分量可以表示為

式中 z為與橫剖面垂直的坐標軸；s為與橫剖面輪廓相切的方向；函數(shù)?(z)，ψ(z)為尋找的廣義位移列向量，僅與坐標z相關(guān)。橫向和縱向分布函數(shù)J(s)，r(s)為廣義坐標s的函數(shù)。

2 超單元剛度矩陣

2.1 無錐度結(jié)構(gòu)

直翼盒模型如圖2所示。蒙皮和腹板簡化為受剪板，縱向筋條受軸向拉壓。

圖2 直翼盒模型示意Fig.2 Diagram of Straight-wing Box Model

模型中，M=3，N=n，即每個截面具有平面內(nèi)分別沿x軸和y軸的平動自由度?1，?2，繞z軸的轉(zhuǎn)動自由度?3以及 n個節(jié)點分別沿 z軸的平動自由度ψ1，ψ2，ψ3，…，ψn。相應(yīng)地，?1，?2，?3分別對應(yīng)于截面上的廣義集中力 Qx，Qy，Mz，ψ1，ψ2，ψ3，…，ψn對應(yīng)的廣義力為作用在相應(yīng)節(jié)點的沿z軸方向的集中力Z1，Z2，Z3，…，ZN。

由于結(jié)構(gòu)沒有錐度，所以每根桿的長度相同，且桿軸向與整體坐標系 z軸平行，故沿桿軸向的積分可以轉(zhuǎn)換為沿 z軸的積分，同時可以交換積分與求和的順序。

對第i根縱向筋條：

對第j塊板：

式中 Sj為第j個板單元沿著s方向的長度；vje和vjb分別表示板單元末端節(jié)點和起始端節(jié)點在 z軸方向上的位移。

寫出由n個桿單元和m個板單元組成的第e個超單元的應(yīng)變能和外力功并進行變分，以矩陣形式表示如下：

在一個超單元范圍內(nèi)，認為可能的位移及虛位移?(z)、ψ(z)、δ?(z)、δψ(z)符合線性變化規(guī)律，即：

將式（9）代入式（7），根據(jù)線性位移假設(shè)進行積分。根據(jù)能量原理：δAe=δUe，經(jīng)過整理得到：

對于小錐度結(jié)構(gòu)，需要對式（11）中的廣義力向量進行修正：

得到單個超單元的剛度矩陣后，按照圖3所示的方式進行組集，即可以得到整個結(jié)構(gòu)的剛度矩陣。

圖3 面內(nèi)剛性假設(shè)下小錐度結(jié)構(gòu)剛度矩陣組裝Fig.3 Assembly of Stiffness Matrix of Small Taper Structure under the Assumption of Cross-section Rigidity

2.2 梯形平面正應(yīng)力板

實際計算時，當(dāng)無錐度或者結(jié)構(gòu)錐度很小時，符合2.1節(jié)中對應(yīng)的情況，當(dāng)錐度較大時，就必須有另外的方法。

與此相關(guān)的主要問題就是要求解四邊形板的剛度矩陣。借助有限元中的經(jīng)典方法，用4個三角形代替四邊形，如圖4所示。

圖4 大錐度盒段模型示意Fig.4 Diagram of Large Taper Box Segment Model

在局部坐標系下桿元的剛度矩陣[9]為

在局部坐標系下，對 4個三角形單元剛度矩陣組集，就得到四邊形單元的剛度矩陣，它有5個節(jié)點：

由于中心5號節(jié)點沒有節(jié)點力，根據(jù)此條件，基于式（15），可以得到5號節(jié)點位移與其余4個節(jié)點位移之間的關(guān)系，進一步得到4節(jié)點四邊形平面應(yīng)力板的剛度矩陣：

對節(jié)點力和位移進行坐標矩陣變換，就將剛度矩陣從局部坐標變換成了總體坐標。

對于桿：

對于板：

根據(jù)上一步結(jié)果，基于節(jié)點力平衡原則，將 n根桿和m塊板的剛度矩陣在整體坐標系下進行組集，得到第e個超單元以有限元法為基礎(chǔ)的剛度矩陣，即：

式中 ue，ue+1分別對應(yīng)于 e和 е+1剖面的位移向量；Fe，F(xiàn)e+1為相應(yīng)的節(jié)點力向量。

將超單元兩端面的真實位移轉(zhuǎn)換成廣義位移：?e、ψe和?e+1、ψe+1分別為剖面 e和 e+1的廣義位移向量；Re、Ze和Re+1、Ze+1分別為剖面e和e+1的廣義力向量。

在超單元剖面上，真實位移和廣義位移之間、廣義力和真實力之間的關(guān)系分別為

式中 Т為位移轉(zhuǎn)置矩陣；λ為力轉(zhuǎn)置矩陣。

且有：

最終，對每個超單元，有以下表達式：

式中 KF_p為對應(yīng)于節(jié)點力、節(jié)點位移的剛度矩陣；KS_p為對應(yīng)于超單元廣義力、廣義位移的剛度矩陣。

由于每一個剖面只有 3個面內(nèi)自由度，N個面外自由度，則通過真實位移向廣義位移轉(zhuǎn)換，超單元的剛度矩陣階數(shù)由 6n×6n降為 2(3+N)×2(3+N)。

2.3 梯形受剪板

2.2 節(jié)中超單元剛度矩陣的計算工作量很大，對于需要多次迭代的設(shè)計和優(yōu)化過程不方便也不可取?；诖?，設(shè)想用更簡單的方法形成板的剛度矩陣，假設(shè)此時板只受剪切。

梯形受剪板如圖5所示，其板厚為t，剪切彈性系數(shù)為G。取梯形板4個角點為節(jié)點，分別記為1、2、3、4。圖中表示節(jié)點i的節(jié)點力。

圖5 受剪梯形板剛度矩陣簡化方式Fig.5 Simplified Method of Stiffness Matrix of Shear Trapezoidal Plate

文獻[9]中給出了梯形受剪板在整體坐標系下的剛度矩陣：

其中，

類似于2.2節(jié)，每一超單元的剛度矩陣由n根桿和m塊板的剛度矩陣按照節(jié)點力平衡原則進行組集，然后通過廣義坐標變換得到。

式中 KF_sp為對應(yīng)于節(jié)點力、節(jié)點位移的剛度矩陣；KS_sp為對應(yīng)于超單元廣義力、廣義位移的剛度矩陣。

3 不同模型的計算比較

3.1 四梁盒段

盒段的形狀、尺寸、載荷如圖6所示，材料為鋁合金，彈性模量E=72 GPa，泊松比μ=0.33。

圖6 算例1模型示意Fig.6 Diagram of Example1 Model

模型1、模型2和模型3是分別對應(yīng)2.1，2.2，2.3節(jié)所提出的建模方法，采用Matlab語言編寫相應(yīng)程序得到的計算結(jié)果。為了便于比較，也利用通用有限元軟件ABAQUS對該算例進行了計算。式（28）為本算例采用虛力原理得到的解析解[10]。

式中 第1項為彎曲引起的變形，第2項為剪切引起的變形。

圖7給出了采用各種計算模型和方法所需的計算時間（以ABAQUS為基準）。采用模型1計算最快，模型3與ABAQUS相當(dāng)，模型2計算時間顯著高于其它模型和方法。圖8給出采用模型1、模型2、模型3、錐度角β=0°時單個超單元長細比對計算誤差的影響。當(dāng)結(jié)構(gòu)的總長度一定時，從圖8可以看出，隨著超單元個數(shù)的增加，即單個超單元長細比的減小，結(jié)構(gòu)最大位移的計算精度逐漸提高，即可通過控制單個超單元長細比（超單元個數(shù)）來獲得不同的計算精度。同時，模型1與模型3在錐度角β=0時計算結(jié)果的精度相當(dāng)，且均低于模型2的計算精度。圖9給出了采用模型2、模型3時，結(jié)構(gòu)錐度角對計算誤差的影響。從圖9中可以看出，超單元法的計算精度隨著結(jié)構(gòu)錐度角的增加而下降。

圖7 各計算方法計算時間比較Fig.7 Comparison of Calculation Time of Each Calculation Method

圖8 錐度角為0時計算誤差與單個超單元長細比的關(guān)系Fig.8 Relationship Curve between the Calculation Error and Long Fineness Ratio of a Single Superelement at a Taper Angle of 0°

圖9 計算誤差與錐度角的關(guān)系Fig.9 Ralationship Curve between Calculation Error and Taper Angle

3.2 小展弦比后掠翼

該算例為帶4根梁和3組桁條的后掠機翼。結(jié)構(gòu)的外形、尺寸及載荷見圖10，材料為鋁合金，彈性模量E=72 GPa，泊松比μ=0.3；梁緣條面積為173 mm2，外面兩組桁條的面積為50 mm2，中間一組桁條面積為76.4 mm2，蒙皮和腹板厚度為3 mm。在結(jié)構(gòu)端部施加集中力Qy，大小為15 000 N，與y軸平行。

圖10 后掠翼示意Fig.10 Rear Wing Diagram

文獻[11]中給出了這個結(jié)構(gòu)的試驗結(jié)果以及計算結(jié)果；模型1中，計算板剛度矩陣時采用2.2節(jié)方法，即將蒙皮和腹板當(dāng)作平面應(yīng)力板處理；模型2中，把板當(dāng)作受剪板處理，此時必須把板的面積（20～30δ2）等效到兩邊的筋條上。圖11給出了沿結(jié)構(gòu)展長方向各方法得到的位移對比。模型 1的計算結(jié)果與文獻[11]值大致相當(dāng)，而最大位移值較文獻值更接近試驗值。與試驗值相比，模型2計算結(jié)果最大誤差為7.7%，計算精度低于模型1的最大誤差5.8%。

圖11 沿展長方向y向位移計算結(jié)果對比Fig.11 Comparison of the Calculation Results of y-direction Displacement Along the Length of the Exhibition

3.3 小錐度機身盒段

該小錐度機身段長度為8 m，半徑由0.6 m過渡到0.75 m，其示意見圖12。沿機身軸向分布24根長桁，其橫截面積為50 mm2；蒙皮厚度為1.6 mm，材料為鋁合金，彈性模量 E=72 GPa，泊松比 μ=0.3。從機身固支端起沿機身長度方向，每隔0.5 m，在橫截面上施加大小為4000 N的集中力Qy，總載荷大小為64 000 N。將機身沿長度方向劃分為16個超單元。

圖12 小錐度機身盒段示意Fig.12 Diagram of Small Taper Fuselage Box Segment

模型1、模型2分別對應(yīng)于應(yīng)力板和受剪板兩種簡化方式。將 Matlab程序求解結(jié)果與有限元軟件ABAQUS分析結(jié)果進行比較，沿機身長度方向橫截面y向位移結(jié)果對比見圖 13。由計算結(jié)果可知，兩種模型的計算結(jié)果與有限元分析結(jié)果相接近，模型1橫截面y向最大位移的計算誤差大于模型2，為3.2%。則超單元法亦適用于機身類薄壁結(jié)構(gòu)。沿著結(jié)構(gòu)的長度方向，超單元法的計算誤差逐漸增大。對于大長細比的薄壁結(jié)構(gòu)，可以通過增加超單元的劃分個數(shù)（即減小單個超單元的長細比）來提高超單元法的計算精度。

圖13 沿機身長度方向y向位移計算結(jié)果對比Fig.13 Comparison of the Calculation Results of y-direction Displacement Along the Length of the Fuselage

4 結(jié) 論

本文提出的針對大長細比薄壁結(jié)構(gòu)的超單元建模方法，單元數(shù)量少，剛度矩陣階數(shù)低，可以大大減少數(shù)據(jù)準備工作量，減少人力消耗。同時超單元方法的計算精度可以根據(jù)不同位移假設(shè)依次提高，適合不同設(shè)計階段的精度要求。對于3種設(shè)計模型，2.1節(jié)中的模型和2.3節(jié)中模型計算效率高，2.2節(jié)中的模型給出的結(jié)果精度最高，但計算效率相對較低。3種模型都可以有效地用于實際結(jié)構(gòu)設(shè)計中。