基于共旋有限元格式的非線性連接翼氣動(dòng)彈性響應(yīng)分析

竇怡彬，孫文釗，樊 浩，梅星磊，曾清香

(上海機(jī)電工程研究所，上海 201109)

0 引 言

連接翼具有較大的展弦比，在飛行過程中變形較大，具有典型的幾何非線性特點(diǎn)。同時(shí)，上翼受到下翼的擠壓作用，承受較大的面內(nèi)壓力導(dǎo)致結(jié)構(gòu)呈現(xiàn)非線性屈曲特性，BLAIR等通過試驗(yàn)驗(yàn)證了這種現(xiàn)象[1]。機(jī)翼在飛行過程中的非線性屈曲現(xiàn)象是連接翼特有的，正是由于這種特點(diǎn)，在機(jī)翼結(jié)構(gòu)的氣動(dòng)彈性問題上，連接翼有別于傳統(tǒng)的單翼布局機(jī)翼。LIVNE總結(jié)了2001年之前關(guān)于連接翼結(jié)構(gòu)氣動(dòng)彈性研究的現(xiàn)狀和面臨的挑戰(zhàn)，指出在進(jìn)行連接翼氣動(dòng)彈性分析時(shí)必須考慮機(jī)翼的屈曲和氣動(dòng)彈性之間的耦合情況，并且屈曲是比顫振影響更大的失穩(wěn)模式[2]。WEISSHAAR等研究了連接翼結(jié)構(gòu)的體自由度顫振問題，認(rèn)為體自由度顫振受機(jī)翼后掠角和飛行器質(zhì)心位置的影響較大，應(yīng)該在設(shè)計(jì)中予以考慮[3]。LEE等采用分段線性方法研究了時(shí)域下考慮屈曲影響的連接翼非線性氣動(dòng)彈性響應(yīng)問題，認(rèn)為響應(yīng)的形式與初始配平條件及陣風(fēng)響應(yīng)載荷有關(guān)[4]。DEMASI等研究了連接翼結(jié)構(gòu)的非線性屈曲行為、非線性氣動(dòng)彈性穩(wěn)定性問題、體自由度顫振和操縱面間隙引起的連接翼極限環(huán)振動(dòng)問題，認(rèn)為研究連接翼結(jié)構(gòu)的非線性屈曲有助于設(shè)計(jì)出氣動(dòng)彈性靜/動(dòng)穩(wěn)定的飛行器，并提出了非線性結(jié)構(gòu)靜氣動(dòng)彈性發(fā)散的突跳-發(fā)散概念[5-7]。SOTOUDEH等提出了連接翼結(jié)構(gòu)幾何非線性分析的增量法，該方法基于完全本征方程，通過數(shù)值算例驗(yàn)證可以用于分析連接翼結(jié)構(gòu)的非線性靜力學(xué)和穩(wěn)定性問題[8]。國(guó)內(nèi)對(duì)于連接翼的研究起步較晚，且大多集中在氣動(dòng)特性分析、線性氣動(dòng)彈性分析、結(jié)構(gòu)優(yōu)化以及可靠性等方面[9-21]，竇怡彬等采用有限元方法研究了不同幾何參數(shù)和載荷分布形式對(duì)Prandtl平板連接翼結(jié)構(gòu)非線性屈曲特性的影響[22]。對(duì)于Prandtl平板連接翼結(jié)構(gòu)的非線性氣動(dòng)彈性瞬態(tài)響應(yīng)研究目前未見公開報(bào)道。

已有的文獻(xiàn)中關(guān)于連接翼結(jié)構(gòu)非線性氣動(dòng)彈性分析大多是在一定的靜態(tài)載荷作用下進(jìn)行穩(wěn)定性分析[23]，或者氣動(dòng)力部分使用簡(jiǎn)單的工程方法進(jìn)行簡(jiǎn)化，導(dǎo)致分析結(jié)果的精度不高、通用性不強(qiáng)[24]。在楊勁松等[25]提出的基于三角形平板殼單元(由離散Kirchhoff板彎單元[26]和優(yōu)化膜單元[27]組合而成)的薄殼結(jié)構(gòu)非線性動(dòng)力學(xué)計(jì)算方法基礎(chǔ)上，本文結(jié)合計(jì)算流體力學(xué)方法和松耦合求解策略，提出了一種薄殼結(jié)構(gòu)非線性氣動(dòng)彈性計(jì)算方法，并以Prandtl平板連接翼為分析對(duì)象，通過數(shù)值仿真方法研究該結(jié)構(gòu)在不同工況條件下的非線性瞬態(tài)響應(yīng)。

1 近似能量守恒算法

由Generalized-α?xí)r間積分算法可知，在tn+1-α?xí)r刻(全文用下標(biāo)n，n+1和n+1-α分別表示tn時(shí)刻，tn+1時(shí)刻和tn+1-α?xí)r刻。文中出現(xiàn)不同時(shí)刻下的同一變量時(shí)，為了簡(jiǎn)單起見只對(duì)其中一個(gè)進(jìn)行變量說明)，結(jié)構(gòu)單元非線性動(dòng)力學(xué)平衡方程可表達(dá)為

ge,n+1-α=fmas,n+1-α+fint,n+1-α-fext,n+1-α=0

(1)

式中：fmas,n+1-α、fint,n+1-α和fext,n+1-α分別表示結(jié)構(gòu)的慣性力、內(nèi)力和外力矢量；ge,n+1-α表示單元的不平衡力矢量；α為積分參數(shù)，本文計(jì)算時(shí)取α=0.5。式(1)中慣性力、內(nèi)力和外力矢量的具體表達(dá)式分別為

(2)

fint,n+1-α=

(3)

fext,n+1-α=(1-α)fext,n+1+αfext,n

(4)

Generalized-α算法和傳統(tǒng)的只利用起始點(diǎn)或者終點(diǎn)時(shí)刻數(shù)據(jù)算法的區(qū)別在于：慣性力、內(nèi)力和外力的平衡不是建立在n+1時(shí)刻或者n時(shí)刻，而是建立在n+1-α?xí)r刻。根據(jù)能量守恒原理，每個(gè)時(shí)間步內(nèi)的結(jié)構(gòu)內(nèi)能和動(dòng)能的變化等于外力在時(shí)間步內(nèi)做的功，該關(guān)系可以通過式(5)所述的平衡方程得到滿足。

(5)

式中：ΔK、ΔS和ΔW分別表示內(nèi)能、動(dòng)能和外力做功的增量；ΔE是時(shí)間步內(nèi)的總能量變化；Δd是全局坐標(biāo)系下迭代計(jì)算過程中的單元位移增量。只要單元的平衡力矢量ge在n+1-α?xí)r刻等于0，在任意位移增量下式(5)就能得到滿足，也就意味著總能量守恒。

2 基于預(yù)估 - 校正過程的響應(yīng)求解

一般可以采用預(yù)估 - 校準(zhǔn)方法對(duì)式(1)進(jìn)行隱式求解，即在時(shí)間段[tn,tn+1]內(nèi)，以tn的結(jié)果為初始值預(yù)估出tn+1時(shí)刻的位移。一般來說這個(gè)預(yù)估結(jié)果并不滿足式(1)，需要對(duì)其進(jìn)行重復(fù)迭代校正，直到其解滿足tn+1時(shí)刻動(dòng)力學(xué)平衡條件。

2.1 預(yù)估步

(6)

(7)

式中：Δua和Δθa分別表示單元節(jié)點(diǎn)a上的平動(dòng)位移增量和轉(zhuǎn)動(dòng)偽矢量增量。將式(6)代入式(2)，并利用關(guān)系Re,n+1=Re,n，整理后得到fmas,n+1-α的預(yù)估值為

(8)

式中：上標(biāo)pred表示預(yù)估步(后文相同)。

令

(9)

(10)

(11)

式中:Kint,n-α和fint,n-α分別是tn-α?xí)r刻的單元內(nèi)力切線剛度矩陣和內(nèi)力，沒有上標(biāo)pred是因?yàn)閠n-α為已知時(shí)刻。將式(8)～(11)代入式(1)，并進(jìn)行單元組裝,得到預(yù)估步動(dòng)力學(xué)方程為

(12)

(13)

(14)

式中：Σ表示組裝單元矩陣。

2.2 校正步

通常來說預(yù)估步得到的位移結(jié)果是無法滿足式(1)的，因此，可將結(jié)構(gòu)動(dòng)態(tài)不平衡力矢量gall,n+1-α在第l次迭代處做一階Taylor展開，得到校正步的第l+1次迭代近似結(jié)果為

(15)

(16)

對(duì)于式(16)右邊第一項(xiàng)和第三項(xiàng)，利用δRa,n+1=spin(δωa)Ra,n+1以及spin(·)函數(shù)(矢量的旋量矩陣運(yùn)算符[25])的性質(zhì)(其中δωa為節(jié)點(diǎn)a的三維無窮小旋轉(zhuǎn)變分)，進(jìn)行簡(jiǎn)化后可得

(17)

(18)

(19)

對(duì)式(6)求變分，然后將其結(jié)果代入式(16)右邊第二項(xiàng)和第四項(xiàng)可得

(20)

式中：Kmas2是式(15)中慣性力切線剛度矩陣的第二個(gè)組成部分；矩陣HΔθ是三維無窮小旋轉(zhuǎn)變分δω和轉(zhuǎn)動(dòng)偽矢量變分δθ間的轉(zhuǎn)換矩陣，其表達(dá)式為

ΗΔθ=diag[I3H(Δθ1)I3H(Δθ2)I3H(Δθ3)]

(21)

將式(17)和式(20)代入式(16)，可得式(15)中由慣性力切線剛度矩陣，即Kmas=Kmas1+Kmas2。同理對(duì)式(3)進(jìn)行變分可得

(22)

根據(jù)單元無關(guān)共旋格式理論可知，式(22)右邊第一項(xiàng)等于幾何剛度矩陣和位移變分的乘積，即

(23)

(24)

(25)

對(duì)式(25)求解后得到迭代步的位移增量，結(jié)構(gòu)的坐標(biāo)直接采用平動(dòng)位移增量更新，采用式(26)和(27)對(duì)節(jié)點(diǎn)的轉(zhuǎn)動(dòng)進(jìn)行更新

(26)

(27)

式(21)中的轉(zhuǎn)動(dòng)偽矢量增量可以采用四元數(shù)法從式(28)中提取。

(28)

3 氣動(dòng)彈性模型

3.1 氣動(dòng)彈性的松耦合方法

圖1 松耦合策略流程Fig.1 Loose coupling procedure

3.2 位移-載荷插值方法

基于松耦合策略的流固耦合分析，結(jié)構(gòu)域和流體域采用各自的網(wǎng)格進(jìn)行空間離散，一般而言這兩套網(wǎng)格之間并不完全重合，因此需要通過插值算法實(shí)現(xiàn)兩者在交接面上的位移-載荷傳遞。常用的方法有徑向基函數(shù)法(RBF)、常體積轉(zhuǎn)換法(CVT)、有限元插值法等。本文提出一種改進(jìn)的常距離徑向基函數(shù)法(CDRBF)，其具體計(jì)算步驟在后續(xù)章節(jié)說明。

3.2.1位移傳遞

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

3.2.2載荷傳遞

(37)

式中：faero和fstru分別表示作用在氣動(dòng)網(wǎng)格和結(jié)構(gòu)網(wǎng)格上的載荷向量；下標(biāo)all表示所有節(jié)點(diǎn)。將式(31)代入式(37)，可得到載荷傳遞關(guān)系為

fstru=HTfaero

(38)

3.3 網(wǎng)格運(yùn)動(dòng)

對(duì)于空間分布的氣動(dòng)網(wǎng)格節(jié)點(diǎn)其位移傳遞也采用RBF方法，和公式(31)相似，氣動(dòng)網(wǎng)格點(diǎn)的位移傳遞公式為

(39)

(40)

(41)

dy=dyRBF(1-e-σy)

(42)

式中：y為氣動(dòng)網(wǎng)格節(jié)點(diǎn)到對(duì)稱面的距離；σ為用戶指定參數(shù)，本文計(jì)算時(shí)取σ=100。

4 數(shù)值仿真

圖2 Prandtl平板連接翼結(jié)構(gòu)Fig.2 Prandtl plane joined-wing configuration

結(jié)構(gòu)模型采用三角形平板殼單元離散，網(wǎng)格數(shù)為2×2×42個(gè)，氣動(dòng)力載荷采用CFD方法計(jì)算得到。流體域采用非結(jié)構(gòu)網(wǎng)格進(jìn)行空間離散，網(wǎng)格總數(shù)為322 481個(gè)單元。圖3給出了連接翼和對(duì)稱面的網(wǎng)格情況；圖4是采用CDRBF方法得到的變形后連接翼氣動(dòng)網(wǎng)格，其中黑色點(diǎn)表示變形后的結(jié)構(gòu)節(jié)點(diǎn)，兩者吻合較好。圖5對(duì)比了RBF方法和CDRBF方法對(duì)連接翼結(jié)構(gòu)位移插值的效果，可以發(fā)現(xiàn)RBF方法插值得到的氣動(dòng)網(wǎng)格物面有較大失真，這主要是由于結(jié)構(gòu)和氣動(dòng)網(wǎng)格幾何外形之間不能完全貼合。采用RBF方法時(shí)位移傳遞引起的誤差會(huì)導(dǎo)致氣動(dòng)力計(jì)算的誤差，進(jìn)而影響整個(gè)計(jì)算結(jié)果的精度，而CDRBF方法利用變形前后氣動(dòng)網(wǎng)格點(diǎn)到結(jié)構(gòu)單元所在平面的距離不變這個(gè)約束條件解決了這個(gè)問題。

圖3 連接翼和對(duì)稱面網(wǎng)格Fig.3 Mesh of joined-wing and symmetrical boundary

圖4 連接翼網(wǎng)格變形(CDRBF)Fig.4 Deformed aerodynamic mesh of joined-wing (CDRBF)

圖6和圖7分別給出了工況1狀態(tài)下P1點(diǎn)和P2點(diǎn)各個(gè)方向上位移隨時(shí)間變化的曲線，可以發(fā)現(xiàn)響應(yīng)是收斂的，說明該動(dòng)壓下結(jié)構(gòu)是穩(wěn)定的，沒有出現(xiàn)發(fā)散現(xiàn)象。圖6中z向位移隨時(shí)間變化可以分為三個(gè)階段，每經(jīng)歷一個(gè)階段z向振動(dòng)的平衡位置都會(huì)向下移動(dòng)，直到最終振動(dòng)到達(dá)平衡位置。圖8給出了幾個(gè)特征時(shí)間點(diǎn)上連接翼的變形構(gòu)型，可以看出，在擾動(dòng)的初始階段，連接翼上下翼都發(fā)生了振動(dòng)，當(dāng)振動(dòng)從第一階段進(jìn)入第二階段后，下翼的振動(dòng)減小，以上翼振動(dòng)為主。圖9給出了結(jié)構(gòu)動(dòng)能、內(nèi)能和總能量的變化曲線，可以看出，在每個(gè)時(shí)間步上結(jié)構(gòu)的總能量守恒，算法穩(wěn)定。

圖5 RBF和CDRBF方法對(duì)連接翼結(jié)構(gòu)位移插值的結(jié)果Fig.5 The RBF and CDRBF method for joined-wing displacement interpolation results

圖6 P1點(diǎn)的位移Fig.6 Displacement components of P1

圖7 P2點(diǎn)的位移Fig.7 Displacement components of P2

圖8 不同時(shí)刻連接翼變形構(gòu)型Fig.8 Deformed configuration of joined-wing at different time

圖9 動(dòng)能、內(nèi)能和總能量曲線Fig.9 Kinetic, internal and total energies

圖10和圖11分別給出了工況2狀態(tài)下P1點(diǎn)和P2點(diǎn)各個(gè)方向的位移，比較圖6和圖10可以發(fā)現(xiàn)，圖10中連接翼結(jié)構(gòu)經(jīng)歷了一般的靜氣動(dòng)彈性過程，即在來流有正攻角的情況下，機(jī)翼發(fā)生向z正方向的變形并達(dá)到平衡位置。圖12比較了工況1(黑色)和工況2(紅色)平衡時(shí)刻的變形構(gòu)型，在不同的動(dòng)壓下，連接翼有著不同的靜氣動(dòng)彈性平衡構(gòu)型：工況1情況下，連接翼弦向低頭同時(shí)展向向下彎曲；工況2情況下，連接翼弦向低頭同時(shí)展向向上彎曲，和普通的后掠大展弦比單翼的靜氣動(dòng)彈性平衡構(gòu)型相同。圖13是工況2情況下的結(jié)構(gòu)動(dòng)能、內(nèi)能和總能量變化曲線，同樣地，在每個(gè)時(shí)間步上能量都能保持守恒且算法穩(wěn)定。

圖10 P1點(diǎn)的位移Fig.10 Displacement components of P1

圖11 P2點(diǎn)的位移Fig.11 Displacement components of P2

圖12 連接翼結(jié)構(gòu)的平衡構(gòu)型Fig.12 Balance configuration of joined-wing

5 結(jié)束語(yǔ)

本文基于Generalized-α?xí)r間積分算法建立了基于共旋有限元格式的薄殼結(jié)構(gòu)非線性動(dòng)力學(xué)響應(yīng)的近似能量守恒算法，并以Prandtl平板連接翼為分析對(duì)象，研究了不同飛行條件下連接翼的非線性氣動(dòng)彈性瞬態(tài)響應(yīng)過程。結(jié)果表明，由于結(jié)構(gòu)的特殊性，Prandtl平板連接翼的氣動(dòng)彈性響應(yīng)和普通單翼結(jié)構(gòu)的氣動(dòng)彈性響應(yīng)有所不同。本文研究了兩種動(dòng)壓條件下連接翼的響應(yīng)，有以下發(fā)現(xiàn)。

1) 當(dāng)攻角為正且來流動(dòng)壓較大時(shí)，連接翼的氣動(dòng)彈性響應(yīng)中出現(xiàn)了繞多個(gè)平衡位置振動(dòng)的現(xiàn)象。在氣動(dòng)阻尼的作用下連接翼結(jié)構(gòu)收斂到靜氣動(dòng)彈性平衡位置，此時(shí)的平衡構(gòu)型弦向低頭且展向向下彎曲，和傳統(tǒng)的后掠單翼變形構(gòu)型不同(弦向低頭展向向上彎曲)。

2) 當(dāng)來流動(dòng)壓較小時(shí)，連接翼的氣動(dòng)彈性響應(yīng)中只出現(xiàn)一個(gè)平衡位置，且其靜氣動(dòng)彈性變形構(gòu)型和后掠單翼的變形構(gòu)型相同。這種現(xiàn)象和連接翼載荷分布以及結(jié)構(gòu)非線性有關(guān)。因此，在結(jié)構(gòu)設(shè)計(jì)過程中需要采用更為精細(xì)的分析手段對(duì)連接翼結(jié)構(gòu)的氣動(dòng)彈性特性進(jìn)行分析。