劉 璟 王 玲，2 胡東飛 鐵 鳴，2 吳旭生

1.北京臨近空間飛行器系統(tǒng)工程研究所，北京100076

2.空間物理重點(diǎn)實(shí)驗(yàn)室，北京100076

臨近空間通常是指20～100km的空域［1］。臨近空間高超聲速飛行器具有快速、機(jī)動(dòng)、靈活的特點(diǎn)，具有極大的發(fā)展?jié)摿Γ?］，對于維護(hù)國家安全和利益具有極為重要的戰(zhàn)略意義。

臨近空間飛行器的研制面臨許多復(fù)雜的問題和挑戰(zhàn)，其功能復(fù)雜、機(jī)理復(fù)雜、工藝復(fù)雜，具有多自由度、多變量、非線性、強(qiáng)耦合、參數(shù)時(shí)變等綜合復(fù)雜性［3］，而且涉及總體、氣動(dòng)、結(jié)構(gòu)、控制、防熱、動(dòng)力等多個(gè)學(xué)科和專業(yè)領(lǐng)域的相互作用、高度耦合的子系統(tǒng)。傳統(tǒng)的各學(xué)科獨(dú)立設(shè)計(jì)與仿真驗(yàn)證的模式難以體現(xiàn)各學(xué)科之間的耦合關(guān)系。需要開展多學(xué)科、全系統(tǒng)、全流程、多物理場耦合過程的協(xié)同仿真研究，驗(yàn)證各學(xué)科耦合關(guān)系下的更接近于真實(shí)情況的飛行性能。學(xué)科間存在著大量的耦合與交互關(guān)系，其中一些涉及學(xué)科間交互的復(fù)雜仿真問題需要多個(gè)學(xué)科領(lǐng)域的仿真模型、軟件相互協(xié)作共同完成［4］。

本文研究并給出了臨近空間高超聲速飛行器的各學(xué)科之間的一種耦合建模方法。由于計(jì)算能力和研制習(xí)慣等方面的制約，氣動(dòng)、結(jié)構(gòu)、防熱等專業(yè)的仿真建模和模型解算工具不同。流體力學(xué)領(lǐng)域側(cè)重于使用有限體積法，固體力學(xué)領(lǐng)域廣泛采用有限元方法，且會(huì)按需求選用不同的計(jì)算工具，包括ABAQUS，ANSYS，NASTRAN和FORTRAN等，這就帶來了不同的網(wǎng)格劃分方式。在多物理場耦合仿真計(jì)算過程中，學(xué)科間復(fù)雜的交互關(guān)系和異構(gòu)網(wǎng)格之間的信息傳遞就成為一個(gè)難點(diǎn)。為解決該問題，研究并實(shí)現(xiàn)了異構(gòu)網(wǎng)格耦合界面之間以及耦合域之間的信息傳遞方法，并應(yīng)用于臨近空間高超聲速飛行器的多物理場耦合仿真。

1 臨近空間高超聲速飛行器多物理場耦合建模

臨近空間高超聲速飛行器各學(xué)科之間的耦合關(guān)系復(fù)雜，其中包括流固耦合問題［5－6］、熱固耦合問題和熱流固耦合問題等。

針對計(jì)算流體力學(xué)(CFD)/計(jì)算固體力學(xué)(CSD)耦合已經(jīng)發(fā)展了全耦合、緊耦合和松耦合［7］3種耦合方式。飛行器氣動(dòng)與結(jié)構(gòu)之間的耦合關(guān)系主要是:在高超聲速條件下，氣動(dòng)載荷引起的飛行器結(jié)構(gòu)表面位移不可忽略，而表面位移反過來會(huì)影響流場和氣動(dòng)載荷。

飛行器熱固領(lǐng)域之間和熱流領(lǐng)域之間的耦合關(guān)系主要是:在高超聲速條件下，氣動(dòng)加熱引起的高溫使得飛行器結(jié)構(gòu)中產(chǎn)生熱應(yīng)力，進(jìn)而導(dǎo)致結(jié)構(gòu)應(yīng)變，溫度的變化也使得材料強(qiáng)度和剛度相應(yīng)改變;飛行器的結(jié)構(gòu)應(yīng)變以及端頭和翼前緣等處的熱燒蝕改變了氣動(dòng)外形，進(jìn)而對飛行器所受的氣動(dòng)力產(chǎn)生影響。熱流固的相互耦合關(guān)系還會(huì)影響飛行器的控制系統(tǒng)精度和可控性。

由于臨近空間高超聲速飛行器的高度復(fù)雜性，目前難以進(jìn)行全系統(tǒng)的多場緊耦合的仿真建模與驗(yàn)證，因此，采用松耦合方法給出了一種多物理場耦合的協(xié)同仿真流程。氣動(dòng)模型基于CFD網(wǎng)格計(jì)算飛行器外表面所受載荷，防熱模型和結(jié)構(gòu)強(qiáng)度模型分別采用ABAQUS和ANSYS作為網(wǎng)格生成工具和模型解算工具。如圖1，在一個(gè)仿真步長內(nèi)，氣動(dòng)、防熱和結(jié)構(gòu)強(qiáng)度仿真模型按其耦合關(guān)系依次解算，步驟序號(hào)說明了解算和數(shù)據(jù)交互的順序，雙箭頭標(biāo)明了數(shù)據(jù)傳遞方向。

圖1 耦合仿真流程圖

2 異構(gòu)網(wǎng)格映射方法

在物理場耦合的各學(xué)科模型之間的信息交互中，存在以下幾類問題:

1)異構(gòu)網(wǎng)格耦合界面位移映射;

2)異構(gòu)網(wǎng)格耦合界面載荷映射;

3)異構(gòu)網(wǎng)格耦合域位移映射;

4)異構(gòu)網(wǎng)格耦合域載荷映射。

其中，耦合界面上映射方法作為流固耦合等問題的一個(gè)重要方面，是一個(gè)研究熱點(diǎn)。常用的方法有常體積轉(zhuǎn)換法［8］、加權(quán)余量法［9］、平板樣條法、反距離加權(quán)法、徑向基函數(shù)法等很多種。

2.1 映射方法的原理和思想

異構(gòu)網(wǎng)格之間信息傳遞方法應(yīng)當(dāng)滿足連續(xù)條件和守恒條件，即網(wǎng)格信息所體現(xiàn)的外形等物理信息是連續(xù)的，網(wǎng)格的劃分只是一種離散方法;映射完成后，兩套網(wǎng)格的整體或局部的能量都必須是守恒的。連續(xù)條件和守恒條件可以作為映射算法的推導(dǎo)依據(jù)或驗(yàn)證依據(jù)。

可將各種異構(gòu)網(wǎng)格之間信息傳遞方法的核心思想總結(jié)為以下幾類:

1)擬合:由網(wǎng)格離散點(diǎn)的信息擬合出連續(xù)的形函數(shù)或場函數(shù)等，并在另一套網(wǎng)格中應(yīng)用該擬合函數(shù)求得各點(diǎn)的值;

2)插值或加權(quán):利用周圍各點(diǎn)對目標(biāo)點(diǎn)影響程度的不同，且相鄰網(wǎng)格點(diǎn)的信息是連續(xù)變化的特點(diǎn)，求得目標(biāo)點(diǎn)的信息;

3)守恒:從方法上保證整體和局部各處的信息在映射前后是嚴(yán)格守恒的，例如基于靜力等效、能量守恒等原理推導(dǎo)出的映射公式。

2.2 本文研究和應(yīng)用的映射方法

比較幾類映射方法，函數(shù)擬合的方法可以在一定程度上還原離散網(wǎng)格點(diǎn)的連續(xù)函數(shù)，但對于復(fù)雜外形或復(fù)雜結(jié)構(gòu)的情況，擬合函數(shù)的精度會(huì)嚴(yán)重下降;守恒類方法推導(dǎo)過程十分嚴(yán)謹(jǐn)，往往存在矩陣求逆過程，難以避免奇異矩陣，在實(shí)際工程應(yīng)用中處理起來比較復(fù)雜;插值或加權(quán)類方法最為常見，因此采用局部信息的插值或加權(quán)方法處理異構(gòu)網(wǎng)格耦合界面和耦合域的信息傳遞問題。

信息從一套網(wǎng)格傳遞到另一套網(wǎng)格，為便于區(qū)分，下文把這兩套網(wǎng)格分別稱作“源網(wǎng)格”和“目標(biāo)網(wǎng)格”。

2.2.1 表面網(wǎng)格載荷和位移映射

常見的表面網(wǎng)格單元的形狀只有三角形和四邊形2種，將四邊形也作為2個(gè)三角形處理。對于耦合界面上的目標(biāo)網(wǎng)格的每個(gè)網(wǎng)格點(diǎn)Ps0都可以唯一確定1個(gè)對應(yīng)的源網(wǎng)格三角單元。將該點(diǎn)Ps0投影在三角單元所在平面上，得到如圖2所示的點(diǎn)Ps。

［8］的面積坐標(biāo)定義，如圖2，得到點(diǎn)Ps的面積坐標(biāo)Ni為:

圖2 面積坐標(biāo)的定義

式中，A為ΔP1P2P3的面積，A1為ΔPsP2P3的面積，A2為ΔPsP1P3的面積，A3為ΔPsP2P1的面積。

處理氣動(dòng)載荷(單位為帕斯卡)映射時(shí)，可以將CFD網(wǎng)格三角單元面積坐標(biāo)視為三角單元各頂點(diǎn)載荷fi對ANSYS網(wǎng)格點(diǎn)Ps的權(quán)重，得到Ps點(diǎn)的載荷為:

對于每個(gè)ANSYS網(wǎng)格點(diǎn)按上述面積權(quán)重法進(jìn)行處理，即可完成異構(gòu)網(wǎng)格耦合界面上的載荷映射過程。

處理位移映射時(shí)，以常體積轉(zhuǎn)換法為基礎(chǔ)，對上述搜索到的三角單元逐個(gè)進(jìn)行計(jì)算。按式(3)和式(4)即可得到耦合界面上網(wǎng)格點(diǎn)Ps0進(jìn)行位移映射后的坐標(biāo)

式中，加撇的點(diǎn)表示該點(diǎn)位移后的坐標(biāo)，如A'表示位移后的三角單元面積。

2.2.2 體網(wǎng)格載荷映射

常見的體網(wǎng)格單元形狀包括四面體、六面體、三棱柱等多種，需要采用一種與網(wǎng)格形狀無關(guān)的方法處理異構(gòu)網(wǎng)格耦合域上的載荷映射問題。

在反距離加權(quán)法的基礎(chǔ)上，以從ABAQUS網(wǎng)格向ANSYS網(wǎng)格進(jìn)行熱載荷映射為例，得出如下映射方法:

1)將兩套網(wǎng)格重疊在一起，在耦合域上搜索每個(gè)ANSYS網(wǎng)格點(diǎn)的坐標(biāo)位于哪個(gè)ABAQUS網(wǎng)格單元中。

已知P1，P2，P3，Pn4個(gè)點(diǎn)的坐標(biāo)，其中，點(diǎn)P1，P2，P3構(gòu)成一個(gè)平面。按行列式的計(jì)算法則，有

判斷點(diǎn)是否在原網(wǎng)格單元內(nèi)的原理是:

將某個(gè)ABAQUS網(wǎng)格單元中心點(diǎn)P0坐標(biāo)和ANSYS網(wǎng)格點(diǎn)Ps分別代入式(5)，計(jì)算結(jié)果同號(hào)，則可判斷點(diǎn)P0與點(diǎn)Ps在平面S123的同側(cè)。

對于一個(gè)均質(zhì)凸n面體(n＞3)，可知其重心(中心)一定在n面體內(nèi)。因此只要按上述方法判斷點(diǎn)Ps與某體網(wǎng)格單元的中心相對于該單元所有的面同側(cè)，即可知該點(diǎn)在該單元內(nèi)部。

2)體網(wǎng)格單元各頂點(diǎn)對點(diǎn)Ps的權(quán)重Ni為:

式中，li表示第i個(gè)頂點(diǎn)與Ps之間的距離。

對于在由n個(gè)頂點(diǎn)構(gòu)成的網(wǎng)格單元內(nèi)的點(diǎn)Ps，其熱載荷映射公式為:

需要指出的是反距離加權(quán)法使用了除法，且除數(shù)有可能為0或接近0，因此在編程實(shí)現(xiàn)時(shí)需要對兩點(diǎn)重合或近似重合的情況進(jìn)行判斷。

2.2.3 體網(wǎng)格位移映射

飛行器結(jié)構(gòu)受力發(fā)生變形時(shí)，在變形區(qū)域，每個(gè)微小的網(wǎng)格單元都會(huì)隨之發(fā)生位移、壓縮或拉伸。由于網(wǎng)格足夠微小，發(fā)生扭曲變形的量所造成的誤差可以忽略。在此假設(shè)基礎(chǔ)上，本文提出一種解決體網(wǎng)格位移映射的“中值位移法”，即由每個(gè)ABAQUS網(wǎng)格節(jié)點(diǎn)所在的由n個(gè)節(jié)點(diǎn)構(gòu)成的ANSYS體網(wǎng)格單元各頂點(diǎn)的平均位移求得該ABAQUS網(wǎng)格節(jié)點(diǎn)的位移，如式(9):

式中，Ai表示ANSYS體網(wǎng)格單元各節(jié)點(diǎn)，Q表示ABAQUS網(wǎng)格節(jié)點(diǎn)，加撇表示位移后的節(jié)點(diǎn)。

式(9)為中值位移法的數(shù)學(xué)表述，該方法的前處理過程，即確定節(jié)點(diǎn)與單元的位置關(guān)系的方法與式(5)所述相同。該方法的優(yōu)點(diǎn)是步驟簡單、易于實(shí)現(xiàn)，計(jì)算效率高。

2.3 異構(gòu)網(wǎng)格信息傳遞算法庫

將異構(gòu)網(wǎng)格之間信息傳遞的各種方法整合到一起，形成異構(gòu)網(wǎng)格信息傳遞算法庫，如圖3所示。算法庫由搜索算法庫、映射算法庫、前處理模塊和后處理模塊組成。前處理模塊按用戶選擇的網(wǎng)格類型動(dòng)態(tài)生成信息傳遞所需的所有數(shù)據(jù)結(jié)構(gòu)。后處理模塊將映射結(jié)果按相應(yīng)的網(wǎng)格文件格式輸出，并釋放內(nèi)存。搜索算法庫內(nèi)包含:節(jié)點(diǎn)與表面網(wǎng)格單元位置關(guān)系搜索算法、節(jié)點(diǎn)與體網(wǎng)格單元位置關(guān)系搜索算法。映射算法庫可以處理各種異構(gòu)網(wǎng)格之間的信息傳遞問題，有包括上文所介紹的常體積轉(zhuǎn)換法、面積權(quán)重法、中值位移法、反距離加權(quán)法在內(nèi)的多種映射算法，并且可以不斷豐富。

2.4 算法庫精度分析

為了驗(yàn)證算法庫的算法精度，選擇了半徑為0.1m的半球體和長0.2m、半徑0.05m的圓柱體模型，各建立2套網(wǎng)格，如圖4和圖5所示。

給半球模型建立溫度場分布函數(shù)，如式(10)所示:

其耦合界面與耦合域上的載荷映射平均誤差見表1。

給圓柱模型表面加載位移分布函數(shù)，如式(11)所示:

δx沿徑向加載，表面以下各層網(wǎng)格等比例壓縮。其耦合界面與耦合域上的載荷映射平均誤差見表1。

可見以上算例的精度符合工程計(jì)算要求，可應(yīng)用于多物理場耦合建模與仿真試驗(yàn)。

3 結(jié)論

研究了臨近空間高飛行器多物理場耦合建模方法，為解決領(lǐng)域間因不同設(shè)計(jì)習(xí)慣、使用不同的工具軟件等造成復(fù)雜數(shù)據(jù)交互問題，研究并實(shí)現(xiàn)了適用于復(fù)雜外形和結(jié)構(gòu)的各種異構(gòu)網(wǎng)格耦合界面及耦合域之間信息傳遞的方法，封裝成算法庫，并可應(yīng)用于臨近空間飛行器多物理場耦合仿真。

圖3 異構(gòu)網(wǎng)格信息傳遞算法庫

表1 算法庫驗(yàn)證結(jié)果

參考文獻(xiàn)

［1］ Stephens H.Near-space［J］.Air Force Magazine，2005，88(7):31.

［2］ 李怡勇，李智，沈懷榮.臨近空間飛行器發(fā)展與應(yīng)用分析［J］.裝備指揮技術(shù)學(xué)院學(xué)報(bào)，2008，19(2):61-65.(LI Yi-yong，LI Zhi，SHEN Huai-rong.Analysis on Development and Application of Near Space Vehicle［J］.Journal of the Academy of Equipment Command＆Technology，2008，19(2):61-65.)

［3］ 鐵鳴，王玲，孫學(xué)功.臨近空間飛行器一體化仿真分析平臺(tái)體系研究［C］.臨近空間飛行器技術(shù)論壇，2009.(TIE Ming，WANG Ling，SUN Xue-gong.Study on Simulation Analysis Platform System for Near Space Vehicle［C］.Nearspace Vehicle’s Technology Forums，2009.)

［4］ 鐵鳴，胡東飛，王玲，吳旭生.基于HLA的臨近空間飛行器虛擬飛行試驗(yàn)方法［J］.信息與電腦，2010，7:51-52.(TIE Ming，HU Dong-fei，WANG Ling，WU Xu-Sheng.Virtual Flight Testing Method of Near Space Vehicle Based on HLA［J］.China Computer＆Communication，2010，7:51-52.)

［5］ 錢若軍，董石麟，袁行飛.流固耦合理論研究進(jìn)展［J］.空間結(jié)構(gòu)，2008，14(1):3-15.(QIAN Ruo-jun，DONG Shi-lin，YUAN Xing-fei.Advances in Research on Fluidstructure Interaction Theory［J］.Spatial Structures，2008，14(1):3-15.)

［6］ 安效民，徐敏，陳士櫓.多場耦合求解非線性氣動(dòng)彈性的研究綜述［J］.力學(xué)進(jìn)展，2009，39(3):284-298.(AN Xiao-min，XU Min，CHEN Shi-lu.An Overview of CFD/ CSD Coupled Solution for Nonlinear Aeroelasticity［J］.Advances in Mechanics，2009，39(3):284-298.)

［7］ Hurka J，Ballmann J.Elastic Panels in Transonic Flow［J］.AIAA Paper 2001-2722，2001.

［8］ Goura G S，Badcock K J，et al.A Data Exchange Method for Fluid-structure Interaction Problems［J］.The Aeronautical Journal，2001，105:215-221.

［9］ De Boer C，Van Zuijlen，Van A H，Bijl H.Review of Coupling Methods for Non-matching Meshes［J］.Computer Methods in Applied Mechanics and Engineering，2007，196:1515-1525.