張 宇, 王曉亮
(上海交通大學(xué) 航空航天學(xué)院,上海 200240)
飛艇能完成諸如通信、偵察、運(yùn)輸?shù)热蝿?wù),在軍事和民事應(yīng)用中具有巨大的前景.2009年3月,上海交通大學(xué)“致遠(yuǎn)一號(hào)”飛艇驗(yàn)證艇在蘇州成功試飛,驗(yàn)證了導(dǎo)航控制、數(shù)傳、測(cè)控等技術(shù)方案的可行性,隨后于同年11月完成飛行試驗(yàn)[1].2014年,美國(guó)國(guó)家航空航天局(NASA)開(kāi)展了“平流層飛艇設(shè)計(jì)”競(jìng)賽,目標(biāo)是以平流層飛艇代替衛(wèi)星[2].2015年,北京航空航天大學(xué)在內(nèi)蒙古成功開(kāi)展了平流層飛艇的長(zhǎng)時(shí)留空試驗(yàn)[3].飛艇從結(jié)構(gòu)上分為硬式飛艇、軟式飛艇和半硬式飛艇[4].其中,軟式飛艇不含內(nèi)部骨架,艇身由柔性蒙皮材料組成,通過(guò)在氣囊內(nèi)注入氦氣形成內(nèi)外壓差從而維持外形.與同體積的硬式飛艇相比,軟式飛艇基于內(nèi)外壓差成型,易于制造和保養(yǎng),且由于不含內(nèi)部骨架、自重較輕、能攜帶更多的任務(wù)載荷,所以小型飛艇多為軟式飛艇.但由于軟式飛艇的結(jié)構(gòu)特點(diǎn),內(nèi)外壓差不易維持容易導(dǎo)致其流固耦合現(xiàn)象明顯.
對(duì)軟式飛艇而言,其流固耦合現(xiàn)象主要來(lái)源于蒙皮薄膜結(jié)構(gòu),針對(duì)薄膜結(jié)構(gòu)的流固耦合分析方法已經(jīng)獲得了較好地發(fā)展.Lee等[5]通過(guò)求解三維Navier-Stokes(N-S)方程獲得氣動(dòng)力,結(jié)構(gòu)分析采用共旋有限元和共旋殼單元方法,分析了流固耦合效應(yīng)對(duì)微型機(jī)器人薄膜撲翼推進(jìn)性能的影響.De Nayer等[6]對(duì)一充氣柔性膜結(jié)構(gòu)在不同來(lái)流雷諾數(shù)下的瞬時(shí)響應(yīng)進(jìn)行了數(shù)值模擬,流場(chǎng)由大渦模擬獲得.在飛艇的應(yīng)用方面,Bessert等[7]通過(guò)流體求解器VSAERO和固體力學(xué)求解軟件Abaqus,研究CL160飛艇在不同俯仰角下的升力系數(shù)變化.王曉亮等[8]通過(guò)Fluent和Abaqus軟件,形成了非線(xiàn)性動(dòng)態(tài)流體-結(jié)構(gòu)交錯(cuò)積分耦合法,使用隱式方法對(duì)某平流層飛艇在突風(fēng)環(huán)境下的結(jié)構(gòu)響應(yīng)進(jìn)行了分析.Liu等[9]通過(guò)SIMPLE算法和非線(xiàn)性有限元方法對(duì)比了彈性與剛性飛艇模型之間的壓力分布與流場(chǎng)參數(shù)差異.吳小翠等[10]研究了在不同剛度構(gòu)型下,飛艇的定常流固耦合特性,并利用量綱分析證明了模型與實(shí)物間的相似律.由此可見(jiàn),目前已經(jīng)初步形成了關(guān)于飛艇的流固耦合分析方法,但上述方法主要針對(duì)飛艇的單向流固耦合特性進(jìn)行研究,且由于隱式方法在求解固體結(jié)構(gòu)響應(yīng)時(shí)需要迭代計(jì)算,對(duì)于高度非線(xiàn)性問(wèn)題可能難以收斂.
本文以Fluent與Abaqus作為求解器,使用雙向流固耦合方法,通過(guò)求解非定常雷諾平均NS(URANS)方程獲取艇身表面氣動(dòng)力,由顯式動(dòng)力學(xué)方法獲得飛艇的結(jié)構(gòu)響應(yīng),由徑向基函數(shù)(RBF)方法進(jìn)行氣動(dòng)力和位移數(shù)據(jù)交換,由Delaunay映射方法完成流場(chǎng)網(wǎng)格更新,由Fortran子程序完成氣動(dòng)載荷加載.最后,應(yīng)用上述框架分析某軟式飛艇在不同壓差下的結(jié)構(gòu)響應(yīng)規(guī)律.
一般的硬式飛艇結(jié)構(gòu)如圖1所示.由于研究對(duì)象為軟式飛艇,所以在硬式飛艇模型的基礎(chǔ)上去除了吊艙和內(nèi)部骨架,艇身外形由參數(shù)化截面(PARSEC)方法獲得[11].為體現(xiàn)艇身的旋成體特點(diǎn),在原有的11個(gè)參數(shù)上只需保留8個(gè)即可確定飛艇外形.PARSEC方法示意圖如圖2所示.其中:rh為頭部半徑;xd為最大橫截面位置;rd為最大橫截面半徑;kd為最大半徑處曲率;αt為尾部偏移角;βt為尾部張開(kāi)角;dt為尾部厚度;ht為尾部高度.艇身外形所對(duì)應(yīng)的各參數(shù)取值如表1所示.飛艇尾翼翼型為NACA0012.y坐標(biāo)由如下的6次多項(xiàng)式方程確定[12]:
表1 艇身外形參數(shù)Tab.1 Shape parameters of hull
圖1 硬式飛艇外形示意圖Fig.1 Schematic diagram of rigid airship
圖2 PARSEC方法Fig.2 PARSEC method
(1)
式中的未知量ai可由如下線(xiàn)性方程組獲得:
(2)
式中:cte為無(wú)量綱翼型弦長(zhǎng),此處取為1.
軟式飛艇所處的環(huán)境示意圖如圖3所示.其中:p為艇身內(nèi)部壓力;p∞為外界遠(yuǎn)場(chǎng)壓力;L為艇身總長(zhǎng);D為最大直徑;denv為蒙皮厚度;Eenv、ρenv和νenv分別為蒙皮彈性模量、密度和泊松比;g為重力加速度,方向垂直向下;“×”表示垂直紙面向里的橫風(fēng).飛艇內(nèi)外壓差為Δp=p-p∞,飛行高度h=20 km,空氣密度與動(dòng)力黏度分別為ρa(bǔ)ir=0.089 kg/m3和μair=1.422 Pa·s.實(shí)際艇身通過(guò)PARSEC模型放大45倍而來(lái),尾翼前緣距離原點(diǎn)0.75L,尾翼后掠角為30°,翼尖長(zhǎng)度為0.08L,兩十字尾翼翼尖距離為1.5D.艇身內(nèi)部壓力和外界遠(yuǎn)場(chǎng)壓力的大小隨高度的變化而變化,來(lái)流速度為v∞.上述有關(guān)變量的取值如表2所示.
圖3 飛艇幾何及其環(huán)境簡(jiǎn)圖Fig.3 Diagram of airship geometry and its environment
表2 飛艇幾何及蒙皮材料參數(shù)Tab.2 Parameters of airship geometry and envelope material
使用連續(xù)性方程和N-S方程即可求解三維不可壓縮流動(dòng)問(wèn)題,微分形式的連續(xù)性方程為
(3)
(4)
湍流模型采用Spalart-Allmaras模型,該模型適合求解航空外流場(chǎng)問(wèn)題,使用壁面增強(qiáng)函數(shù)以便獲得較好的近壁面流場(chǎng)信息.求解器采取不可壓縮形式的壓力基求解器,壓力-速度耦合格式為“Coupled”,空間梯度離散方法為“Least Squares Cell Based”,壓力項(xiàng)采用二階格式離散,動(dòng)量和修正湍流黏度采用二階迎風(fēng)格式離散.CFD時(shí)間增量參數(shù)為dtf,每一增量步內(nèi)迭代步數(shù)設(shè)置為50以保證計(jì)算收斂.
流體計(jì)算域如圖4(a)所示,速度入口距離飛艇中心20L=900 m,壓力出口距離飛艇中心40L=1 800 m,遠(yuǎn)場(chǎng)邊界距離飛艇中心20L=900 m.網(wǎng)格劃分工具為Pointwise,在飛艇表面布置三角形網(wǎng)格,采取T-Rex方法生成邊界層網(wǎng)格,遠(yuǎn)場(chǎng)流域采用非結(jié)構(gòu)四面體網(wǎng)格填充.
圖4 流體計(jì)算域及其網(wǎng)格Fig.4 Fluid computational domain and its meshes
為降低流場(chǎng)網(wǎng)格因素帶來(lái)的誤差,有必要進(jìn)行網(wǎng)格無(wú)關(guān)性檢驗(yàn).通過(guò)對(duì)艇身表面網(wǎng)格疏密程度的調(diào)整,達(dá)到既能保證結(jié)果精度又能節(jié)約計(jì)算資源的目的.為此,劃分了4種疏密程度的網(wǎng)格,分別統(tǒng)計(jì)在v∞=12 m/s橫風(fēng)作用下飛艇z方向的阻力系數(shù)CDz,相關(guān)飛艇流場(chǎng)的網(wǎng)格信息如表3所示.由表3可知,中等疏密程度的網(wǎng)格已經(jīng)能較好地滿(mǎn)足計(jì)算結(jié)果的精度要求.為節(jié)約計(jì)算資源,選取中等網(wǎng)格劃分策略進(jìn)行后續(xù)分析,最終的流場(chǎng)網(wǎng)格如圖4(b)所示.
表3 飛艇流場(chǎng)網(wǎng)格信息Tab.3 Mesh information of airship fluid
在上述數(shù)值方法及網(wǎng)格單元設(shè)置下,飛艇受到v∞=12 m/s橫風(fēng)作用時(shí)艇身表面與不同x剖面交線(xiàn)上的無(wú)量綱網(wǎng)格高度y+分布如圖5所示.其中:
圖5 不同交線(xiàn)上的y+分布Fig.5 Distribution of y+ on different intersection lines
徑向坐標(biāo)為艇身表面的y+值;圓周坐標(biāo)θ為交線(xiàn)繞x軸的夾角.由圖5可知,在不同艇身交線(xiàn)上的y+最大值均約為3,大部分處于0.5~2范圍內(nèi).故所選取的數(shù)值方法及網(wǎng)格劃分策略滿(mǎn)足增強(qiáng)壁面函數(shù)及湍流模型對(duì)邊界層質(zhì)量的要求,能較好地捕捉艇身外流流場(chǎng)信息.
流場(chǎng)網(wǎng)格更新通過(guò)Delaunay映射方法實(shí)現(xiàn),該方法針對(duì)中等變形問(wèn)題具有很高的穩(wěn)健性,避免了求解大型矩陣的高耗時(shí)缺點(diǎn).關(guān)于該算法的具體原理參考文獻(xiàn)[13].借助Fluent中的用戶(hù)自定義函數(shù)(UDF)功能,將Delaunay方法集成在DEFINE_GRID_MOTION宏中,即可實(shí)現(xiàn)網(wǎng)格更新效果.
對(duì)于結(jié)構(gòu)動(dòng)態(tài)響應(yīng),一般可通過(guò)隱式或顯式動(dòng)力學(xué)算法獲得.隱式算法需在每個(gè)時(shí)間步內(nèi)迭代求解平衡方程,當(dāng)前后兩次迭代值滿(mǎn)足誤差要求時(shí)才會(huì)進(jìn)入下一時(shí)間步,而顯式算法在任一時(shí)間步內(nèi)無(wú)需迭代求解大型矩陣.相較于隱式算法,在結(jié)構(gòu)設(shè)置合理的前提下,顯式算法不存在收斂問(wèn)題.因此,隱式算法一般只用于求解常規(guī)的動(dòng)力學(xué)問(wèn)題,而顯式算法更適合高度非線(xiàn)性問(wèn)題且能更好地捕捉動(dòng)態(tài)響應(yīng)細(xì)節(jié).
顯式動(dòng)力學(xué)關(guān)于運(yùn)動(dòng)方程的中心差分形式可表示為
(5)
(6)
(7)
固體網(wǎng)格采用三角形殼單元?jiǎng)澐?,?jīng)過(guò)網(wǎng)格無(wú)關(guān)性分析后確定網(wǎng)格單元總數(shù)為 3 180.采取基于中心差分的顯式動(dòng)力學(xué)對(duì)運(yùn)動(dòng)方程進(jìn)行積分求解.但中心差分算法是條件穩(wěn)定的,其時(shí)間增量與應(yīng)力波傳播速度vd與網(wǎng)格特征尺寸Lc密切相關(guān),在此設(shè)定CSD的時(shí)間增量參數(shù)為dts,其值一般要滿(mǎn)足dts≤Lc/vd.計(jì)算時(shí)將十字尾翼固定,在飛艇內(nèi)部施加均勻分布的法向壓力作為內(nèi)壓,由CFD計(jì)算得到的氣動(dòng)力通過(guò)Fortran語(yǔ)言二次開(kāi)發(fā)子程序施加到軟式飛艇蒙皮外表面.
采取RBF實(shí)現(xiàn)氣動(dòng)力與位移的傳遞.RBF插值方法不依賴(lài)于網(wǎng)格拓?fù)浣Y(jié)構(gòu),只須獲取CFD和CSD的插值坐標(biāo)即可進(jìn)行插值.需要說(shuō)明的是,進(jìn)行數(shù)據(jù)傳遞時(shí)CFD與CSD的插值點(diǎn)分別為物面網(wǎng)格節(jié)點(diǎn)和殼單元中心點(diǎn).
假設(shè)已知S個(gè)結(jié)構(gòu)網(wǎng)格單元中心的坐標(biāo)?s和某一方向上的位移u,則可建立關(guān)于位移的插值格式為
u=Φω
(8)
式中:Φ為徑向基函數(shù)矩陣;ω為系數(shù)向量.式(8)的展開(kāi)形式為
(9)
式中:ξji為點(diǎn)?sj與點(diǎn)?si的距離與徑向基函數(shù)作用半徑ref的比值,即
(10)
徑向基函數(shù)選取 Wendland二階光順格式[14],可表示為
(11)
由于矩陣Φ是對(duì)稱(chēng)正定的,其逆可通過(guò)Cholesky分解快速獲得,所以可以獲得系數(shù)向量ω的分量值,于是CFD網(wǎng)格節(jié)點(diǎn)的位移可由下式獲得:
(12)
k=1,2,…,Sf
式中:m為結(jié)構(gòu)殼單元編號(hào);k為流體網(wǎng)格節(jié)點(diǎn)編號(hào);Sf為流體物面網(wǎng)格節(jié)點(diǎn)個(gè)數(shù).
氣動(dòng)力傳遞過(guò)程與上節(jié)類(lèi)似,在此不再贅述.如果所面臨的流體插值點(diǎn)數(shù)目遠(yuǎn)大于固體插值點(diǎn)數(shù)目的情況,可使用局部徑向基函數(shù)插值方法來(lái)提高效率,具體過(guò)程參考文獻(xiàn)[15].
關(guān)于軟式飛艇的非定常流固耦合流程如圖6所示.一個(gè)子循環(huán)步的耦合過(guò)程如下:① 給軟式飛艇內(nèi)部施加均勻分布的法向壓力,并在外部施加幅值為0的氣動(dòng)力,此步驟用于形成因內(nèi)部壓力而變形的飛艇初始外形;② CSD模塊在當(dāng)前子循環(huán)步中以增量dts向前推進(jìn)Δt,判斷CSD模塊向前推進(jìn)的時(shí)間是否達(dá)到dtf,即是否滿(mǎn)足|Δt-dtf|≤ε,若滿(mǎn)足則進(jìn)行下一步,否則繼續(xù)進(jìn)行CSD模塊計(jì)算,其中ε為時(shí)間允許誤差,在數(shù)值上等于dts;③ 將CSD模塊的單元中心變形量通過(guò)RBF方法傳遞給CFD模塊網(wǎng)格節(jié)點(diǎn),并由Delaunay映射方法更新網(wǎng)格節(jié)點(diǎn)位置,開(kāi)啟CFD模塊計(jì)算進(jìn)程;④ 當(dāng)CFD模塊完成當(dāng)前時(shí)間步的計(jì)算后,將氣動(dòng)力傳遞到CSD網(wǎng)格單元.重復(fù)上述過(guò)程即可完成雙向流固耦合的交替計(jì)算,且從上述過(guò)程可知實(shí)際的耦合步長(zhǎng)為dtf.
圖6 流固耦合計(jì)算流程Fig.6 Calculation process of fluid-structure interaction
為驗(yàn)證所提方法的可靠性,以如圖7所示的基于網(wǎng)格的并行化代碼耦合接口(MpCCI)案例為研究對(duì)象,將密度為 1 000 kg/m3、彈性模量為0.1 GPa、泊松比為0.49的彈性板置于長(zhǎng)方體導(dǎo)管中.入口速度為8 m/s,出口為壓力出口,表壓為0,管壁為無(wú)滑移壁面,彈性板的上端固定在管壁上,取彈性板下端頂點(diǎn)為控制點(diǎn),其余參數(shù)及設(shè)置參考文獻(xiàn)[16].分別利用MpCCI和本文所提方法(FSI-ED)耦合Fluent和Abaqus進(jìn)行計(jì)算,最終獲得控制點(diǎn)的振動(dòng)幅值隨時(shí)間的變化情況如圖8所示.其中:ux為控制點(diǎn)在來(lái)流方向上的位移.由圖8可知,MpCCI與FSI-ED計(jì)算所得的振動(dòng)頻率分別為44.8 Hz與45.0 Hz,其相對(duì)誤差為0.45%,所得最大振幅的平均相對(duì)誤差為6.1%.
圖7 平板沖擊案例計(jì)算域(m)Fig.7 Computational domain of plate impact case (m)
圖8 FSI-ED與MpCCI計(jì)算結(jié)果對(duì)比Fig.8 Comparison of calculation results of FSI-ED and MpCCI
選取NACA0014機(jī)翼進(jìn)行氣動(dòng)力傳遞的精度驗(yàn)證,其機(jī)翼表面網(wǎng)格如圖9所示.來(lái)流馬赫數(shù)為 0.839 5,攻角為3.06°,當(dāng)?shù)卮髿鈮簽?100 311.75 Pa,當(dāng)?shù)販囟葹?00.9 K.表4統(tǒng)計(jì)了通過(guò)積分獲得的作用在流體和固體網(wǎng)格表面單元上3個(gè)方向的合力和力矩大小,其中:Fx、Fy和Fz為機(jī)翼在x、y和z方向受到的合力;Mx、My和Mz為機(jī)翼在x、y和z方向受到的合力矩.由表4可知,在氣動(dòng)力的傳遞過(guò)程中幾乎沒(méi)有能量損失,所有對(duì)比項(xiàng)的相對(duì)誤差均在0.1%以下.
表4 傳遞前后的力與力矩對(duì)比Tab.4 Comparison of forces and moments before and after transmission
圖9 NACA0014機(jī)翼表面網(wǎng)格Fig.9 Surface mesh of NACA0014 wing
綜上,所提數(shù)值方法能較精確地對(duì)振動(dòng)過(guò)程進(jìn)行描述且能以較高的精度保證氣動(dòng)力的傳遞,可滿(mǎn)足流固耦合問(wèn)題的數(shù)值計(jì)算.
對(duì)于流固耦合分析而言,時(shí)間步長(zhǎng)的選取至關(guān)重要,不當(dāng)?shù)臅r(shí)間步長(zhǎng)甚至可能導(dǎo)致錯(cuò)誤的結(jié)果.CSD模塊的時(shí)間增量步dts與應(yīng)力波傳遞速度vd和網(wǎng)格特征尺寸Lc相關(guān),這兩個(gè)參數(shù)直接反映了飛艇模型的固有頻率.Abaqus在模型數(shù)據(jù)檢查階段會(huì)自動(dòng)計(jì)算上述兩個(gè)參數(shù),并對(duì)顯式時(shí)間步長(zhǎng)進(jìn)行估計(jì).根據(jù)其估計(jì)量選取dts=5×10-5s.
對(duì)CFD模型而言,一般采取隱式算法,因此理論上其時(shí)間步長(zhǎng)不受限制,但過(guò)小的時(shí)間步長(zhǎng)會(huì)導(dǎo)致計(jì)算資源有所增加,過(guò)大的時(shí)間步長(zhǎng)可能導(dǎo)致失去真實(shí)的結(jié)構(gòu)響應(yīng)細(xì)節(jié).對(duì)軟式飛艇而言,其振動(dòng)頻率大多低于5 Hz[17].出于對(duì)計(jì)算資源及計(jì)算精度的綜合考量,選取dtf=1×10-2s.
設(shè)定軟式飛艇的內(nèi)外壓差Δp的范圍為[208,852] Pa,增量為92 Pa.應(yīng)用所提計(jì)算框架對(duì)不同壓差下整體艇身的結(jié)構(gòu)響應(yīng)進(jìn)行分析.
對(duì)能維持一定整體剛度的軟式飛艇而言,在尾翼固定的情況下,艇身頭部的振動(dòng)幅值最大,能夠較好地反映軟式飛艇在橫風(fēng)作用下的振動(dòng)特性[8],故取飛艇頭部頂點(diǎn)作為控制點(diǎn)進(jìn)行監(jiān)測(cè).在不同壓差情況下,飛艇控制點(diǎn)在時(shí)域上的響應(yīng)曲線(xiàn)如圖10所示,其中uz為控制點(diǎn)在z方向上的位移.從時(shí)域曲線(xiàn)可見(jiàn),當(dāng)Δp低于576 Pa時(shí),艇身振動(dòng)幅值包線(xiàn)隨時(shí)間單調(diào)遞減;當(dāng)Δp高于668 Pa后,振動(dòng)幅值包線(xiàn)不再單調(diào)變化,而是呈現(xiàn)出較為復(fù)雜的動(dòng)態(tài)響應(yīng).通過(guò)快速Fourier變換對(duì)時(shí)域曲線(xiàn)進(jìn)行轉(zhuǎn)換得到如圖11所示的頻域響應(yīng)曲線(xiàn),其中:f為頻率;A為幅值.Δp與fm的關(guān)系曲線(xiàn)如圖12所示,通過(guò)對(duì)不同壓差下的振動(dòng)主頻fm進(jìn)行統(tǒng)計(jì),可見(jiàn)振動(dòng)主頻與飛艇的內(nèi)外壓差之間呈近似線(xiàn)性關(guān)系:
fm=5.8×10-5Δp+3.667 9
(13)
圖10 控制點(diǎn)的時(shí)域響應(yīng)曲線(xiàn)Fig.10 Time domain response curves of control points
圖11 控制點(diǎn)的頻域響應(yīng)曲線(xiàn)Fig.11 Frequency domain response curves of control points
圖12 Δp與fm的關(guān)系曲線(xiàn)Fig.12 Δp versus fm
由式(13)可見(jiàn),斜率僅為5.8×10-5,因此可認(rèn)為在文中設(shè)定的壓差范圍內(nèi),振動(dòng)主頻fm幾乎不受壓差Δp的影響,即當(dāng)軟式飛艇受到的壓差能基本維持艇身剛度之后,繼續(xù)增加內(nèi)壓對(duì)振動(dòng)頻率的影響較小,這與相關(guān)試驗(yàn)結(jié)論吻合[17].
本文基于URANS和顯式動(dòng)力學(xué)方程,通過(guò)Delaunay映射方法完成流場(chǎng)網(wǎng)格更新,采用徑向基函數(shù)完成氣動(dòng)力和結(jié)構(gòu)位移在流體單元和固體單元之間的傳遞,進(jìn)而構(gòu)建軟式飛艇在非定常氣動(dòng)力激勵(lì)下的顯示動(dòng)力學(xué)結(jié)構(gòu)響應(yīng)計(jì)算框架.該框架適合雙向流固耦合計(jì)算,能較好地反映軟式飛艇的結(jié)構(gòu)響應(yīng)規(guī)律,體現(xiàn)顯式計(jì)算格式效率高的優(yōu)點(diǎn).飛艇艇身由PARSEC方法獲得,應(yīng)用上述框架對(duì)該軟式飛艇的計(jì)算結(jié)果表明,當(dāng)軟式飛艇能基本維持艇身剛度后,繼續(xù)增加內(nèi)壓對(duì)振動(dòng)頻率的影響較小,計(jì)算結(jié)果與相關(guān)試驗(yàn)結(jié)論吻合.本文形成的計(jì)算框架對(duì)于高空氣球以及其他膜結(jié)構(gòu)高空飛行器在瞬態(tài)氣動(dòng)力作用下的結(jié)構(gòu)響應(yīng)分析均適用.