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

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

0 引 言

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

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

1 近似能量守恒算法

由Generalized-α時間積分算法可知，在tn+1-α時刻(全文用下標n，n+1和n+1-α分別表示tn時刻，tn+1時刻和tn+1-α時刻。文中出現(xiàn)不同時刻下的同一變量時，為了簡單起見只對其中一個進行變量說明)，結構單元非線性動力學平衡方程可表達為

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

(1)

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

(2)

fint,n+1-α=

(3)

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

(4)

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

(5)

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

2 基于預估 - 校正過程的響應求解

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

2.1 預估步

(6)

(7)

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

(8)

式中：上標pred表示預估步(后文相同)。

令

(9)

(10)

(11)

式中:Kint,n-α和fint,n-α分別是tn-α時刻的單元內(nèi)力切線剛度矩陣和內(nèi)力，沒有上標pred是因為tn-α為已知時刻。將式(8)～(11)代入式(1)，并進行單元組裝,得到預估步動力學方程為

(12)

(13)

(14)

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

2.2 校正步

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

(15)

(16)

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

(17)

(18)

(19)

對式(6)求變分，然后將其結果代入式(16)右邊第二項和第四項可得

(20)

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

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

(21)

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

(22)

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

(23)

(24)

(25)

對式(25)求解后得到迭代步的位移增量，結構的坐標直接采用平動位移增量更新，采用式(26)和(27)對節(jié)點的轉動進行更新

(26)

(27)

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

(28)

3 氣動彈性模型

3.1 氣動彈性的松耦合方法

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

3.2 位移-載荷插值方法

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

3.2.1位移傳遞

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

3.2.2載荷傳遞

(37)

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

fstru=HTfaero

(38)

3.3 網(wǎng)格運動

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

(39)

(40)

(41)

dy=dyRBF(1-e-σy)

(42)

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

4 數(shù)值仿真

圖2 Prandtl平板連接翼結構Fig.2 Prandtl plane joined-wing configuration

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

圖3 連接翼和對稱面網(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點和P2點各個方向上位移隨時間變化的曲線，可以發(fā)現(xiàn)響應是收斂的，說明該動壓下結構是穩(wěn)定的，沒有出現(xiàn)發(fā)散現(xiàn)象。圖6中z向位移隨時間變化可以分為三個階段，每經(jīng)歷一個階段z向振動的平衡位置都會向下移動，直到最終振動到達平衡位置。圖8給出了幾個特征時間點上連接翼的變形構型，可以看出，在擾動的初始階段，連接翼上下翼都發(fā)生了振動，當振動從第一階段進入第二階段后，下翼的振動減小，以上翼振動為主。圖9給出了結構動能、內(nèi)能和總能量的變化曲線，可以看出，在每個時間步上結構的總能量守恒，算法穩(wěn)定。

圖5 RBF和CDRBF方法對連接翼結構位移插值的結果Fig.5 The RBF and CDRBF method for joined-wing displacement interpolation results

圖6 P1點的位移Fig.6 Displacement components of P1

圖7 P2點的位移Fig.7 Displacement components of P2

圖8 不同時刻連接翼變形構型Fig.8 Deformed configuration of joined-wing at different time

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

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

圖10 P1點的位移Fig.10 Displacement components of P1

圖11 P2點的位移Fig.11 Displacement components of P2

圖12 連接翼結構的平衡構型Fig.12 Balance configuration of joined-wing

5 結束語

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

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

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