基于顯式動力學的軟式飛艇流固耦合計算框架

張 宇, 王曉亮

(上海交通大學 航空航天學院，上海 200240)

飛艇能完成諸如通信、偵察、運輸等任務，在軍事和民事應用中具有巨大的前景.2009年3月，上海交通大學“致遠一號”飛艇驗證艇在蘇州成功試飛，驗證了導航控制、數傳、測控等技術方案的可行性，隨后于同年11月完成飛行試驗[1].2014年，美國國家航空航天局(NASA)開展了“平流層飛艇設計”競賽，目標是以平流層飛艇代替衛星[2].2015年，北京航空航天大學在內蒙古成功開展了平流層飛艇的長時留空試驗[3].飛艇從結構上分為硬式飛艇、軟式飛艇和半硬式飛艇[4].其中，軟式飛艇不含內部骨架，艇身由柔性蒙皮材料組成，通過在氣囊內注入氦氣形成內外壓差從而維持外形.與同體積的硬式飛艇相比，軟式飛艇基于內外壓差成型，易于制造和保養，且由于不含內部骨架、自重較輕、能攜帶更多的任務載荷，所以小型飛艇多為軟式飛艇.但由于軟式飛艇的結構特點，內外壓差不易維持容易導致其流固耦合現象明顯.

對軟式飛艇而言，其流固耦合現象主要來源于蒙皮薄膜結構，針對薄膜結構的流固耦合分析方法已經獲得了較好地發展.Lee等[5]通過求解三維Navier-Stokes(N-S)方程獲得氣動力，結構分析采用共旋有限元和共旋殼單元方法，分析了流固耦合效應對微型機器人薄膜撲翼推進性能的影響.De Nayer等[6]對一充氣柔性膜結構在不同來流雷諾數下的瞬時響應進行了數值模擬，流場由大渦模擬獲得.在飛艇的應用方面，Bessert等[7]通過流體求解器VSAERO和固體力學求解軟件Abaqus，研究CL160飛艇在不同俯仰角下的升力系數變化.王曉亮等[8]通過Fluent和Abaqus軟件，形成了非線性動態流體-結構交錯積分耦合法，使用隱式方法對某平流層飛艇在突風環境下的結構響應進行了分析.Liu等[9]通過SIMPLE算法和非線性有限元方法對比了彈性與剛性飛艇模型之間的壓力分布與流場參數差異.吳小翠等[10]研究了在不同剛度構型下，飛艇的定常流固耦合特性，并利用量綱分析證明了模型與實物間的相似律.由此可見，目前已經初步形成了關于飛艇的流固耦合分析方法，但上述方法主要針對飛艇的單向流固耦合特性進行研究，且由于隱式方法在求解固體結構響應時需要迭代計算，對于高度非線性問題可能難以收斂.

本文以Fluent與Abaqus作為求解器，使用雙向流固耦合方法，通過求解非定常雷諾平均NS(URANS)方程獲取艇身表面氣動力，由顯式動力學方法獲得飛艇的結構響應，由徑向基函數(RBF)方法進行氣動力和位移數據交換，由Delaunay映射方法完成流場網格更新，由Fortran子程序完成氣動載荷加載.最后，應用上述框架分析某軟式飛艇在不同壓差下的結構響應規律.

1 幾何及環境描述

一般的硬式飛艇結構如圖1所示.由于研究對象為軟式飛艇，所以在硬式飛艇模型的基礎上去除了吊艙和內部骨架，艇身外形由參數化截面(PARSEC)方法獲得[11].為體現艇身的旋成體特點，在原有的11個參數上只需保留8個即可確定飛艇外形.PARSEC方法示意圖如圖2所示.其中：rh為頭部半徑；xd為最大橫截面位置；rd為最大橫截面半徑；kd為最大半徑處曲率；αt為尾部偏移角；βt為尾部張開角；dt為尾部厚度；ht為尾部高度.艇身外形所對應的各參數取值如表1所示.飛艇尾翼翼型為NACA0012.y坐標由如下的6次多項式方程確定[12]：

表1 艇身外形參數Tab.1 Shape parameters of hull

圖1 硬式飛艇外形示意圖Fig.1 Schematic diagram of rigid airship

圖2 PARSEC方法Fig.2 PARSEC method

(1)

式中的未知量ai可由如下線性方程組獲得：

(2)

式中：cte為無量綱翼型弦長，此處取為1.

軟式飛艇所處的環境示意圖如圖3所示.其中：p為艇身內部壓力；p∞為外界遠場壓力；L為艇身總長；D為最大直徑；denv為蒙皮厚度；Eenv、ρenv和νenv分別為蒙皮彈性模量、密度和泊松比；g為重力加速度，方向垂直向下；“×”表示垂直紙面向里的橫風.飛艇內外壓差為Δp=p-p∞，飛行高度h=20 km，空氣密度與動力黏度分別為ρair=0.089 kg/m3和μair=1.422 Pa·s.實際艇身通過PARSEC模型放大45倍而來，尾翼前緣距離原點0.75L，尾翼后掠角為30°，翼尖長度為0.08L，兩十字尾翼翼尖距離為1.5D.艇身內部壓力和外界遠場壓力的大小隨高度的變化而變化，來流速度為v∞.上述有關變量的取值如表2所示.

圖3 飛艇幾何及其環境簡圖Fig.3 Diagram of airship geometry and its environment

表2 飛艇幾何及蒙皮材料參數Tab.2 Parameters of airship geometry and envelope material

2 數值方法

2.1 計算流體力學(CFD)數值方法

使用連續性方程和N-S方程即可求解三維不可壓縮流動問題，微分形式的連續性方程為

(3)

(4)

湍流模型采用Spalart-Allmaras模型，該模型適合求解航空外流場問題，使用壁面增強函數以便獲得較好的近壁面流場信息.求解器采取不可壓縮形式的壓力基求解器，壓力-速度耦合格式為“Coupled”，空間梯度離散方法為“Least Squares Cell Based”，壓力項采用二階格式離散，動量和修正湍流黏度采用二階迎風格式離散.CFD時間增量參數為dtf，每一增量步內迭代步數設置為50以保證計算收斂.

流體計算域如圖4(a)所示，速度入口距離飛艇中心20L=900 m，壓力出口距離飛艇中心40L=1 800 m，遠場邊界距離飛艇中心20L=900 m.網格劃分工具為Pointwise，在飛艇表面布置三角形網格，采取T-Rex方法生成邊界層網格，遠場流域采用非結構四面體網格填充.

圖4 流體計算域及其網格Fig.4 Fluid computational domain and its meshes

2.2 可靠性驗證

為降低流場網格因素帶來的誤差，有必要進行網格無關性檢驗.通過對艇身表面網格疏密程度的調整，達到既能保證結果精度又能節約計算資源的目的.為此，劃分了4種疏密程度的網格，分別統計在v∞=12 m/s橫風作用下飛艇z方向的阻力系數CDz，相關飛艇流場的網格信息如表3所示.由表3可知，中等疏密程度的網格已經能較好地滿足計算結果的精度要求.為節約計算資源，選取中等網格劃分策略進行后續分析，最終的流場網格如圖4(b)所示.

表3 飛艇流場網格信息Tab.3 Mesh information of airship fluid

在上述數值方法及網格單元設置下，飛艇受到v∞=12 m/s橫風作用時艇身表面與不同x剖面交線上的無量綱網格高度y+分布如圖5所示.其中：

圖5 不同交線上的y+分布Fig.5 Distribution of y+ on different intersection lines

徑向坐標為艇身表面的y+值；圓周坐標θ為交線繞x軸的夾角.由圖5可知，在不同艇身交線上的y+最大值均約為3，大部分處于0.5～2范圍內.故所選取的數值方法及網格劃分策略滿足增強壁面函數及湍流模型對邊界層質量的要求，能較好地捕捉艇身外流流場信息.

2.3 動網格更新方法

流場網格更新通過Delaunay映射方法實現，該方法針對中等變形問題具有很高的穩健性，避免了求解大型矩陣的高耗時缺點.關于該算法的具體原理參考文獻[13].借助Fluent中的用戶自定義函數(UDF)功能，將Delaunay方法集成在DEFINE_GRID_MOTION宏中，即可實現網格更新效果.

2.4 計算結構動力學(CSD)數值方法

對于結構動態響應，一般可通過隱式或顯式動力學算法獲得.隱式算法需在每個時間步內迭代求解平衡方程，當前后兩次迭代值滿足誤差要求時才會進入下一時間步，而顯式算法在任一時間步內無需迭代求解大型矩陣.相較于隱式算法，在結構設置合理的前提下，顯式算法不存在收斂問題.因此，隱式算法一般只用于求解常規的動力學問題，而顯式算法更適合高度非線性問題且能更好地捕捉動態響應細節.

顯式動力學關于運動方程的中心差分形式可表示為

(5)

(6)

(7)

固體網格采用三角形殼單元劃分，經過網格無關性分析后確定網格單元總數為 3 180.采取基于中心差分的顯式動力學對運動方程進行積分求解.但中心差分算法是條件穩定的，其時間增量與應力波傳播速度vd與網格特征尺寸Lc密切相關，在此設定CSD的時間增量參數為dts，其值一般要滿足dts≤Lc/vd.計算時將十字尾翼固定，在飛艇內部施加均勻分布的法向壓力作為內壓，由CFD計算得到的氣動力通過Fortran語言二次開發子程序施加到軟式飛艇蒙皮外表面.

3 數據傳遞方法及驗證

采取RBF實現氣動力與位移的傳遞.RBF插值方法不依賴于網格拓撲結構，只須獲取CFD和CSD的插值坐標即可進行插值.需要說明的是，進行數據傳遞時CFD與CSD的插值點分別為物面網格節點和殼單元中心點.

3.1 位移傳遞

假設已知S個結構網格單元中心的坐標?s和某一方向上的位移u，則可建立關于位移的插值格式為

u=Φω

(8)

式中：Φ為徑向基函數矩陣；ω為系數向量.式(8)的展開形式為

(9)

式中：ξji為點?sj與點?si的距離與徑向基函數作用半徑ref的比值，即

(10)

徑向基函數選取 Wendland二階光順格式[14]，可表示為

(11)

由于矩陣Φ是對稱正定的，其逆可通過Cholesky分解快速獲得，所以可以獲得系數向量ω的分量值，于是CFD網格節點的位移可由下式獲得：

(12)

k=1,2,…,Sf

式中：m為結構殼單元編號；k為流體網格節點編號；Sf為流體物面網格節點個數.

3.2 氣動力傳遞

氣動力傳遞過程與上節類似，在此不再贅述.如果所面臨的流體插值點數目遠大于固體插值點數目的情況，可使用局部徑向基函數插值方法來提高效率，具體過程參考文獻[15].

3.3 流固耦合過程

關于軟式飛艇的非定常流固耦合流程如圖6所示.一個子循環步的耦合過程如下：① 給軟式飛艇內部施加均勻分布的法向壓力，并在外部施加幅值為0的氣動力，此步驟用于形成因內部壓力而變形的飛艇初始外形；② CSD模塊在當前子循環步中以增量dts向前推進Δt，判斷CSD模塊向前推進的時間是否達到dtf，即是否滿足|Δt-dtf|≤ε，若滿足則進行下一步，否則繼續進行CSD模塊計算，其中ε為時間允許誤差，在數值上等于dts；③ 將CSD模塊的單元中心變形量通過RBF方法傳遞給CFD模塊網格節點，并由Delaunay映射方法更新網格節點位置，開啟CFD模塊計算進程；④ 當CFD模塊完成當前時間步的計算后，將氣動力傳遞到CSD網格單元.重復上述過程即可完成雙向流固耦合的交替計算，且從上述過程可知實際的耦合步長為dtf.

圖6 流固耦合計算流程Fig.6 Calculation process of fluid-structure interaction

3.4 流固耦合方法驗證

為驗證所提方法的可靠性，以如圖7所示的基于網格的并行化代碼耦合接口(MpCCI)案例為研究對象，將密度為 1 000 kg/m3、彈性模量為0.1 GPa、泊松比為0.49的彈性板置于長方體導管中.入口速度為8 m/s，出口為壓力出口，表壓為0，管壁為無滑移壁面，彈性板的上端固定在管壁上，取彈性板下端頂點為控制點，其余參數及設置參考文獻[16].分別利用MpCCI和本文所提方法(FSI-ED)耦合Fluent和Abaqus進行計算，最終獲得控制點的振動幅值隨時間的變化情況如圖8所示.其中：ux為控制點在來流方向上的位移.由圖8可知，MpCCI與FSI-ED計算所得的振動頻率分別為44.8 Hz與45.0 Hz，其相對誤差為0.45%，所得最大振幅的平均相對誤差為6.1%.

圖7 平板沖擊案例計算域(m)Fig.7 Computational domain of plate impact case (m)

圖8 FSI-ED與MpCCI計算結果對比Fig.8 Comparison of calculation results of FSI-ED and MpCCI

選取NACA0014機翼進行氣動力傳遞的精度驗證，其機翼表面網格如圖9所示.來流馬赫數為 0.839 5，攻角為3.06°，當地大氣壓為 100 311.75 Pa，當地溫度為300.9 K.表4統計了通過積分獲得的作用在流體和固體網格表面單元上3個方向的合力和力矩大小，其中：Fx、Fy和Fz為機翼在x、y和z方向受到的合力；Mx、My和Mz為機翼在x、y和z方向受到的合力矩.由表4可知，在氣動力的傳遞過程中幾乎沒有能量損失，所有對比項的相對誤差均在0.1%以下.

表4 傳遞前后的力與力矩對比Tab.4 Comparison of forces and moments before and after transmission

圖9 NACA0014機翼表面網格Fig.9 Surface mesh of NACA0014 wing

綜上，所提數值方法能較精確地對振動過程進行描述且能以較高的精度保證氣動力的傳遞，可滿足流固耦合問題的數值計算.

4 軟式飛艇的動態結構響應

4.1 流固耦合時間步長

對于流固耦合分析而言，時間步長的選取至關重要，不當的時間步長甚至可能導致錯誤的結果.CSD模塊的時間增量步dts與應力波傳遞速度vd和網格特征尺寸Lc相關，這兩個參數直接反映了飛艇模型的固有頻率.Abaqus在模型數據檢查階段會自動計算上述兩個參數，并對顯式時間步長進行估計.根據其估計量選取dts=5×10-5s.

對CFD模型而言，一般采取隱式算法，因此理論上其時間步長不受限制，但過小的時間步長會導致計算資源有所增加，過大的時間步長可能導致失去真實的結構響應細節.對軟式飛艇而言，其振動頻率大多低于5 Hz[17].出于對計算資源及計算精度的綜合考量，選取dtf=1×10-2s.

4.2 結構振動響應

設定軟式飛艇的內外壓差Δp的范圍為[208,852] Pa，增量為92 Pa.應用所提計算框架對不同壓差下整體艇身的結構響應進行分析.

對能維持一定整體剛度的軟式飛艇而言，在尾翼固定的情況下，艇身頭部的振動幅值最大，能夠較好地反映軟式飛艇在橫風作用下的振動特性[8]，故取飛艇頭部頂點作為控制點進行監測.在不同壓差情況下，飛艇控制點在時域上的響應曲線如圖10所示，其中uz為控制點在z方向上的位移.從時域曲線可見，當Δp低于576 Pa時，艇身振動幅值包線隨時間單調遞減；當Δp高于668 Pa后，振動幅值包線不再單調變化，而是呈現出較為復雜的動態響應.通過快速Fourier變換對時域曲線進行轉換得到如圖11所示的頻域響應曲線，其中：f為頻率；A為幅值.Δp與fm的關系曲線如圖12所示，通過對不同壓差下的振動主頻fm進行統計，可見振動主頻與飛艇的內外壓差之間呈近似線性關系：

fm=5.8×10-5Δp+3.667 9

(13)

圖10 控制點的時域響應曲線Fig.10 Time domain response curves of control points

圖11 控制點的頻域響應曲線Fig.11 Frequency domain response curves of control points

圖12 Δp與fm的關系曲線Fig.12 Δp versus fm

由式(13)可見，斜率僅為5.8×10-5，因此可認為在文中設定的壓差范圍內，振動主頻fm幾乎不受壓差Δp的影響，即當軟式飛艇受到的壓差能基本維持艇身剛度之后，繼續增加內壓對振動頻率的影響較小，這與相關試驗結論吻合[17].

5 結語

本文基于URANS和顯式動力學方程，通過Delaunay映射方法完成流場網格更新，采用徑向基函數完成氣動力和結構位移在流體單元和固體單元之間的傳遞，進而構建軟式飛艇在非定常氣動力激勵下的顯示動力學結構響應計算框架.該框架適合雙向流固耦合計算，能較好地反映軟式飛艇的結構響應規律，體現顯式計算格式效率高的優點.飛艇艇身由PARSEC方法獲得，應用上述框架對該軟式飛艇的計算結果表明，當軟式飛艇能基本維持艇身剛度后，繼續增加內壓對振動頻率的影響較小，計算結果與相關試驗結論吻合.本文形成的計算框架對于高空氣球以及其他膜結構高空飛行器在瞬態氣動力作用下的結構響應分析均適用.