基于廣義位移法的薄壁結(jié)構(gòu)超單元建模方法

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

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

0 引 言

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

在超單元法中，整體結(jié)構(gòu)被分為若干個(gè)超單元，每個(gè)超單元可以視為一系列有限元素的集合。超單元法通過對系統(tǒng)自由度進(jìn)行減縮，從而大幅減少計(jì)算量，適用于大型復(fù)雜結(jié)構(gòu)的分析[3]。Guyan[4]將結(jié)構(gòu)自由度分為主自由度和副自由度，通過強(qiáng)制忽略副自由度的影響，對剛度矩陣和質(zhì)量矩陣進(jìn)行縮聚，適用于靜力分析和低階模態(tài)分析；陳國平[5,6]將界面自由度和內(nèi)部自由度均進(jìn)行減縮，提出了適用于結(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é)點(diǎn)位移縮聚為廣義位移，對于3種簡化模型分別推導(dǎo)出超單元?jiǎng)偠染仃?，通過Matlab編寫相應(yīng)求解程序。通過3個(gè)算例結(jié)果與有限元結(jié)果的對比，驗(yàn)證了該方法的準(zhǔn)確性和有效性。

1 薄壁梁中的超單元法

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

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

續(xù)圖1

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

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

2 超單元?jiǎng)偠染仃?/h2>

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

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

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

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

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

對第i根縱向筋條：

對第j塊板：

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

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

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

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

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

得到單個(gè)超單元的剛度矩陣后，按照圖3所示的方式進(jìn)行組集，即可以得到整個(gè)結(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)力板

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

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

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

在局部坐標(biāo)系下桿元的剛度矩陣[9]為

在局部坐標(biāo)系下，對 4個(gè)三角形單元?jiǎng)偠染仃嚱M集，就得到四邊形單元的剛度矩陣，它有5個(gè)節(jié)點(diǎn)：

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

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

對于桿：

對于板：

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

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

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

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

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

且有：

最終，對每個(gè)超單元，有以下表達(dá)式：

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

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

2.3 梯形受剪板

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

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

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

文獻(xiàn)[9]中給出了梯形受剪板在整體坐標(biāo)系下的剛度矩陣：

其中，

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

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

3 不同模型的計(jì)算比較

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)程序得到的計(jì)算結(jié)果。為了便于比較，也利用通用有限元軟件ABAQUS對該算例進(jìn)行了計(jì)算。式（28）為本算例采用虛力原理得到的解析解[10]。

式中 第1項(xiàng)為彎曲引起的變形，第2項(xiàng)為剪切引起的變形。

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

圖7 各計(jì)算方法計(jì)算時(shí)間比較Fig.7 Comparison of Calculation Time of Each Calculation Method

圖8 錐度角為0時(shí)計(jì)算誤差與單個(gè)超單元長細(xì)比的關(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 計(jì)算誤差與錐度角的關(guān)系Fig.9 Ralationship Curve between Calculation Error and Taper Angle

3.2 小展弦比后掠翼

該算例為帶4根梁和3組桁條的后掠機(jī)翼。結(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

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

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

3.3 小錐度機(jī)身盒段

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

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

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

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

4 結(jié) 論

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