國家數(shù)值風洞工程結構網格變形程序SGDP V1.0開發(fā)與應用

孫 巖,江 盟,孟德虹,黃 勇

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所,綿陽 621000)

0 引 言

物面變形或運動下的非定常流動數(shù)值模擬、氣動外形優(yōu)化設計、靜/動氣動彈性計算中外形隨時間發(fā)生變化,需要對流場計算網格進行更新以獲取新外形下的流場,動態(tài)網格生成技術是解決這一問題的有效途徑[1]。動態(tài)網格生成主要有兩種實現(xiàn)方式:網格變形與網格重構。相比網格重構,網格變形具有算法簡單、計算量小、自動化程度高和不改變網格拓撲及不引入新的離散誤差等優(yōu)點,能夠減小網格變化對數(shù)值模擬結果的影響,在氣動優(yōu)化和流/固耦合模擬中被廣泛采用[2]。但網格變形技術也存在一個明顯的缺點,即在變形后網格質量會降低,尤其是大變形情況下,可能會出現(xiàn)網格單元交叉或負體積單元的現(xiàn)象,導致流場模擬失敗。另一方面,在非定常計算中,網格變形方法會被反復調用,網格更新的時間消耗是影響流動模擬效率的關鍵因素。因此,發(fā)展高效/魯棒/高質量的網格變形算法一直是動態(tài)網格生成領域的研究熱點[3]。

近十年,一些優(yōu)秀的網格變形算法相繼被提出,如Delaunay背景映射[4]、徑向基插值(Radial Basis Functions,RBF)[5-9]、反向距離加權(Inverse Distance Weight,IDW)[10]等代數(shù)方法,逐步取代早期采用基于偏微分方程求解的方法,成為目前網格變形采用的主要思路。而一些綜合兩種或多種方法的復合型網格變形策略也展現(xiàn)出更好的網格變形能力和更新效率,例如RBF-Delaunay[11]、RBFs-MSA[12]、RBF-TFI[13]等,其中RBF-TFI方法因其優(yōu)異網格變形能力在多種問題中得到應用[14]。

盡管網格變形方法逐步趨向成熟和多樣化,已經能夠解決十分復雜的工程問題[15],但目前的算法實現(xiàn)主要以研究型代碼為主,尚未有高效/魯棒/易用的模塊化網格變形程序供學術研究和工程問題模擬使用。

國家數(shù)值風洞(National Numerical Windtunnel,NNW)項目是中國政府在2018年底啟動的大型自主軟件研發(fā)項目,致力于發(fā)展以計算流體動力學(Computational Fluid Dynamics,CFD)為核心的流體數(shù)值模擬軟件群,構建CFD 產業(yè)應用生態(tài)圈,通過以點帶面,推動中國計算機輔助設計行業(yè)的發(fā)展,助力未來中國制造業(yè)向著智能化、數(shù)字化方向升級。作為NNW 項目共性基礎庫的重要組成部分,網格變形技術模塊在吸收前期課題組研究代碼的基礎上,通過邏輯梳理、代碼重構、接口規(guī)范、功能擴展和算例測試,首先實現(xiàn)了多塊結構網格變形程序(Structured Grid Deformation Program,SGDP)V1.0 版 本 的 開發(fā)。目前,SGDP V1.0正在NNW 項目內部開展測試,待成熟完善后,將在全國范圍內推廣使用,作為學術研究或工程問題模擬的有力工具逐步應用于氣動外形優(yōu)化、流/固耦合模擬等領域。

本文的目的是通過對SGDP V1.0程序的理論方法、組成架構、主要功能、運行方法、運行性能和工程應用等進行簡要介紹,讓使用者對SGDP 程序有更深入的認識,促進用戶之間進行交流和發(fā)現(xiàn)問題,推動SGDP的持續(xù)良性發(fā)展,為整個CFD 行業(yè)提供性能更加優(yōu)異的基礎程序庫。

1 SGDP V1.0簡介

1.1 理論方法

結構網格變形程序SGDP V1.0采用了徑向基函數(shù)(Radial Basis Functions,RBF)與超限插值(Transfinite Interpolation,TFI)結合的復合型動態(tài)網格變形方法來對物面變形后的網格進行更新。RBF方法能夠獲得高質量的變形網格,但標準形式的RBF為了保證物面的精確變形,需要將所有物面網格點選做徑向基基點,導致插值矩陣龐大、計算效率低和魯棒性差等問題。TFI方法基于棱線的網格點分布特性即可重構出新的面和體網格,具有良好的計算效率和網格質量,但TFI無法對物面變形后的網格棱線進行更新。復合方法RBF-TFI利用RBF對物面變形后的空間棱線坐標進行更新,然后利用TFI方法重構得到變形后的空間面網格和體網格。RBF-TFI方法中,變形后的物面網格直接利用物面網格點位移進行更新得到,物面的精準變形得到保證,從而RBF可以選擇少量的徑向基基點,能夠極大改善網格變形的效率和魯棒性。RBF-TFI方法具體的網格更新流程如圖1所示,其中為了保證物面和空間棱線交叉位置的一致性,選擇基點時需要將空間棱線與物面連接的點選做徑向基基點?；c選擇算法首先選擇所有物面網格窗口的角點,然后調用重復點查找模塊去除重復的角點,并通過“貪婪”算法選擇其他物面網格點作為剩余的徑向基基點。

圖1 RBF-TFI方法網格更新流程Fig.1 Grid updating procedure of RBF-TFI method

SGDP V1.0中采用了不帶多項式的RBF,可實現(xiàn)遠場網格保持靜止,并對網格塊的棱線進行坐標更新。

徑向基函數(shù)插值可以表示為:

式中:d 為網格點的位移,x 為網格點的坐標矢量,N為徑向基插值基點的數(shù)量,αi(i＝1,2,…,N)為插值系數(shù),φ 為基函數(shù),‖x－xi‖為兩點之間的距離。

兩點之間距離‖x－xi‖可以表示為:

已知N 個插值基點位置xi(i＝1,2,…,N)的位移di,利用式(1)可以構建關于N 個關于插值系數(shù)ai的線性方程組

式中

求解方程組(3)得到插值系數(shù)ai,然后利用式(1)可以插值得到任意位置x 的位移d(x)。

變形網格質量與基函數(shù)的類型有關,根據大量網格算例測試結果,本文選擇網格質量調控能力較好的Wendland’s C2函數(shù),具體表示為:

式中,η＝‖x－xi‖/R,其中R 為基函數(shù)的影響半徑,更多有關RBF的理論可以參考文獻[5]。

SGDP V1.0中采用了Soni等發(fā)展的標準形式TFI[16],能夠更好地處理復雜外形和極端分布下的網格構造,有關TFI方法的詳細可以參考文獻[16]。

1.2 組成架構

圖2給出了SGDP V1.0程序的文件架構,整個程序采用Fortran95語言編寫,其中main.f90定義了程序的主入口,然后通過launch.f90文件可以選擇調用不同的網格變形方法,目前僅定義了RBF-TFI方法,通過launch.f90 可以快速擴展其他網格變形方法。文件rbf_tfi.f90定義了網格變形的流程和具體實現(xiàn),包括物面網格提取、物面網格位移計算、徑向基基點選擇、插值與系數(shù)矩陣計算、面/體網格更新等,grid.f90定義了網格數(shù)據結構、網格數(shù)據、網格質量分析和網格輸入/輸出等模塊,tfi.f90定義了二維/三維TFI的實現(xiàn)函數(shù),logo.f90定義了一些程序的logo和版權信息,distance.f90定義了不同距離計算的函數(shù),base.f90定義了矩陣運算、字符串處理、Tecplot文件IO 和重復點查找等基礎模塊。

SGDP V1.0程序的功能架構主要由網格文件輸入/輸出、網格變形、變形結果輸出三個大模塊組成。網格文件輸入/輸出模塊用于讀入和寫出結構網格數(shù)據與邊界信息定義文件,目前能夠處理2種常用結構網格格式:Plot3d 和Gridgen,3 種文件存儲格式:Binary、Unformatted、Ascii,2 種 數(shù) 據 精 度:Single Precision和Double Precision,組合共計12種不同的結構網格數(shù)據文件格式,邊界信息文件格式采用Gridgen中Generic Solver定義的方式,使用Ascii格式進行存儲。網格變形模塊根據不同的變形類型,計算物面網格點的位移運動,然后插值得到空間棱線網格點的位移,并利用TFI算法對空間網格數(shù)據進行構造更新。變形結果輸出模塊將變形前后網格質量分布、徑向基基點選點結果、基點插值誤差分布、變形外形等寫出到指定的文件,便于程序運行后對網格變形效果進行分析。

圖2 SGDP V1.0程序文件架構Fig.2 File architecture of SGDP V1.0

1.3 主要功能

SGDP V1.0 程序目前支持4 種類型的網格變形:平動、轉動、彎曲和用戶自定義。平動通過輸入參數(shù)文件中X、Y、Z 方向的位移進行定義;轉動為繞固定軸轉動,目前僅支持繞X、Y、Z 三個軸轉動,通過輸入參數(shù)文件中的轉動軸、轉動中心和轉動角度定義;彎曲主要用于描述沿特定方向的彎曲變形,通過輸入參數(shù)文件中的最大彎曲變形、彎曲變形起始位置、彎曲變形最大位置進行定義;用戶自定義通過讀入文件的方式定義物面網格點的位移,可以描述任意種類的變形,輸入文件的路徑、文件名在參數(shù)文件中指定。

在實際工程應用中,單一的運動或變形方式非常少見,因此平動、轉動、彎曲3種變形方式主要用來對程序進行測試和演示算例制作,用戶自定義方式用來實現(xiàn)復雜變形或運動的輸入。例如,氣動彈性計算或氣動外形優(yōu)化設計中,流場物面網格點的位移文件可以通過結構變形插值或者幾何參數(shù)化模塊輸出獲得,然后調用SGDP V1.0,即可得到物面變形后的流場計算網格。

此外,SGDP V1.0通過對更新前后的網格單元數(shù)據進行分析,可以獲得每個網格塊的所有網格單元質量,并統(tǒng)計得到每個網格塊的最差單元質量(worst)、最優(yōu)單元質量(best)和平均單元質量(average),然后輸出變形前后網格單元質量在不同塊中的分布特性。

1.4 運行流程

SGDP V1.0程序目前可以在Windows和Linux平臺下運行,均采用命令行和配置文件的方式啟動。

在Window 下啟動CMD 或在Linux 下啟動Terminal,定位至工作目錄后,SGDP V1.0可以通過下面3種方式運行:

1)SGDP

2)SGDP project_config

3)SGDP method project_config

其中,第1種運行方式,默認采用RBF-TFI作為網格變形方法,config.txt文件作為默認參數(shù)配置文件;第2中運行方式,RBF-TFI作為默認網格變形方法,配置參數(shù)文件的名稱project_config可以用戶自己定義,例如SGDP naca0012.txt;第3種運行方式,網格變形方法method和配置參數(shù)文件名稱project_config均可以用戶自己定義,例如SGDP RBF_TFI naca0012.txt。

SGDP V1.0程序運行前,需要準備好原始網格和邊界信息文件,并配置好輸入參數(shù)文件。SGDP V1.0借鑒了斯坦福非結構流場解算軟件SU2的參數(shù)輸入模式,采用扁平化參數(shù)組織結構,參數(shù)讀入沒有先后順序,從而可以方便地添加或刪除任意的參數(shù),參數(shù)文件中沒有定義的參數(shù)將使用程序中的默認定義值。圖3給出了某個網格變形算例的輸入參數(shù)文件示例,這里采用了用戶自定義的變形輸入方式,只對相關的參數(shù)進行了定義。

圖3 SGDP V1.0輸入參數(shù)文件定義示例Fig.3 Input parameter file example of SGDP V1.0

SGDP V1.0程序啟動后,會對讀入的部分特定參數(shù)進行解析,并進行參數(shù)安全檢查,然后將全部的參數(shù)信息輸出到屏幕,方便用戶了解程序運行情況,及時發(fā)現(xiàn)參數(shù)輸入錯誤的情況并進行相應修正。圖4給出了SGDP V1.0的全部運行參數(shù)打印信息。有關每個參數(shù)的具體含義,可以參考SGDP V1.0的使用手冊。

2 性能測試

計算時間消耗是網格變形程序的關鍵性能指標,在重復調用程序的仿真工程作業(yè)中,對整個仿真的時間和計算成本有著重要的影響。本節(jié)利用DLR-F6翼身組合體構型[17],采用不同規(guī)模的多塊結構網格,對SGDP V1.0的運行時間性能進行測試。

測試中采用的不同規(guī)模網格,基于1套基礎網格數(shù)據,利用課題組自主開發(fā)的多塊結構網格自適應稀疏/加密程序,成比例的加密得到,不同規(guī)模網格的拓撲結構和網格塊數(shù)量是不變的。基礎網格數(shù)據含有72 個網格塊,有3 868 996個網格點、3 566 336個六面體單元,圖5給出了DLR-F6翼身組合體外形和基礎網格的拓撲結構。

圖4 SGDP V1.0全部運行參數(shù)輸出信息Fig.4 Printing information of all running parameters of SGDP V1.0

本次測試,網格加密系數(shù)在0.5～4.0之間變化,網格點數(shù)量從50萬增長到2.3億,選擇DLR-F6 模型繞質心旋轉20°的算例,徑向基插值最大基點數(shù)目為300。程序在Windows 7系統(tǒng)上運行,計算機硬件采用Intel i7-8700 CPU 處理器,含有12個處理器單元,主頻為3.2GHz,內存為32Gb。SGDP V1.0目前只有串行版本,采用單CPU 方式運行。

圖6給出了SGDP V1.0程序總的CPU 運行時間隨網格點數(shù)目Nnode的變化規(guī)律。隨著Nnode的增加,網格變形消耗的CPU 時間在逐步增大,但增加的速率是減小的,即每百萬網格點消耗的時間隨著網格點數(shù)的增加逐漸減小。

SGDP V1.0程序總的運行時間主要包括4個部分:網格文件讀取、徑向基基點選擇、體網格更新和網格文件寫出。圖7給出了SGDP V1.0程序運行時間分布隨網格節(jié)點數(shù)目的變化圖,基于“貪婪”法的徑向基基點選擇占據程序全部運行時間的絕大部分,比例隨著網格點數(shù)據的增加逐步減小。這是由于貪婪選點算法是將物面網格點位移誤差最大的點逐步選擇出來,算法的計算量正比于物面網格點數(shù)量。而隨著網格規(guī)模增加,物面網格點數(shù)量是按照加密比例平方增加,而體網格點數(shù)量是按照加密比例立方增加,所以物面網格點數(shù)量的增加速度要遠小于體網格點的增加速度。網格更新時間、網格數(shù)據讀入/寫出時間占總運行時間的比例隨著網格點數(shù)目逐步增加的,但網格讀入/寫出時間占總運行時間的比例非常小,基本在1.0%以下。

圖5 DLR-F6翼身組合體外形及基礎網格拓撲結構Fig.5 Shape and basis grid topology of DLR-F6 wing-body configuration

圖6 SGDP V1.0總運行時間隨網格節(jié)點數(shù)目的變化Fig.6 Total CPU time cost of SGDP V1.0 with different grid nodes

圖7 不同網格節(jié)點數(shù)目下SGDP V1.0運行時間分布Fig.7 CPU time cost percent of SGDP V1.0 with different grid nodes

對于網格變形程序多次調用的算例,由于物面精確變形不是通過徑向基函數(shù)插值RBF來保證,所以基點貪婪選擇算法可以利用最大變形只執(zhí)行一次,然后將基點在物面網格點上的索引編號通過數(shù)組或文件記錄下來,下一次程序調用直接使用此編號去獲取徑向基基點的坐標,然后更新系數(shù)和插值矩陣,將極大地改善程序的運行效率。對于這類情況,例如氣動彈性多次迭代或優(yōu)化中不同參數(shù)下氣動外形變化,SGDP V1.0程序的運行時間主要依賴于網格讀寫時間和體網格更新時間。而對于利用內存進行網格數(shù)據傳遞的例子,網格讀取時間也將是一次性的。

圖8給出了網格讀寫時間、網格變形更新時間隨總的網格點數(shù)目的變化規(guī)律。可以看出,網格讀寫時間、網格更新時間隨網格點數(shù)目基本是線性變化的趨勢,每百萬網格點的讀取時間約為0.01 s,每百萬網格點的寫出時間約為0.008 s,每百萬網格點的更新時間約為0.1 s。因此,對于千萬量級的網格,SGDP V1.0程序可以在數(shù)秒內完成一次變形網格生成,綜合來說是十分高效的。

圖8 網格讀寫時間、更新時間隨網格節(jié)點數(shù)目的變化Fig.8 CPU time cost of grid reading,writing and updating with different grid nodes

3 變形算例

本節(jié)提供一些典型情況下的網格變形算例,用來驗證SGDP V1.0網格變形能力和網格質量維持特性。

圖9給出了NACA0012翼型繞軸旋轉90°前后的變形網格,其中旋轉中心位于弦線1/4位置。旋轉情況下,變形向遠場的衰減過程會使網格單元產生明顯的剪切變形,對網格單元質量有著顯著的影響,如圖9(b)。

變形前,NACA0012網格最差單元質量為0.800,變形后,NACA0012 網格最差單元質量為0.597。90°旋轉后,盡管網格單元質量有了顯著的降低,但依然維持較高的水平(＞0.30),說明SGDP V1.0在顯著旋轉變形下仍然具有很好的網格質量維持能力。從圖9(b)還可以看出,物面附近網格隨物面一起運動,能夠保證物面附近黏性流場計算區(qū)域網格單元保持良好的連續(xù)性和正交性。

圖9 NACA0012旋轉90°的變形網格Fig.9 Deformed grid of NACA0012 airfoil with 90°rotation

圖10 為某大展弦比柔性機翼模型在氣動載荷作用下變形前后的流場計算網格。該大展弦比機翼半展長為1542,翼梢最大彎曲變形達到402,約為半展長的26%。在如此大幅度變形下,機翼附近網格單元仍然能夠隨著機翼一起運動,保證了變形后的網格質量,說明SGDP V1.0具有處理大變形情況的能力。

圖10 大展弦比柔性機翼模型變形前后網格剖面Fig.10 Grid section of high aspect-ratio flexible wing model before and after deformation

圖11 為NASA 共同研究模型翼/身/掛架/短艙組合體構型CRM-WBPN 彎曲變形前后的流場計算網格。CRM-WBPN 模型計算網格有153 個網格塊和約1200萬網格點(半模),網格拓撲在機翼/掛架/短艙連接位置十分復雜。從圖11 可以看出,SGDP V1.0能夠很好地處理大彎曲變形下CRM 的網格變形,表明SGDP V1.0具有非常好的復雜外形網格變形處理能力。

圖11 NASA CRM-WBPN 模型變形前后網格剖面Fig.11 Grid section of NASA CRM-WBPN model before and after deformation

圖12 為NASA CRM-WBPN 大彎曲變形前后單個網格塊最差和最優(yōu)單元質量隨網格塊的分布。變形前后,最差、最優(yōu)單元質量基本沒有變化,依然具有較高的網格質量,甚至某些塊變形后最差單元質量還有所改善(第88塊),說明SGDP V1.0在大變形情況具有十分優(yōu)異的網格質量維持能力。

4 工程應用

目前,結構網格變形程序SGDP V1.0已經被完整地集成到了亞跨超流場模擬平臺(Trisonic Platform,TRIP)中,用于開展不同類型的靜/動氣動彈性耦合數(shù)值模擬,如機翼載荷分配、副翼效率、風洞模型結構靜變形影響修正[18-19]、戰(zhàn)斗機垂尾抖振、機翼顫振[15]、空調管道流/固耦合振動等。

圖12 NASA CRM-WBPN 模型變形前后最差、最優(yōu)網格單元質量分布Fig.12 Worst and best grid cell quality distribution of NASA CRM-WBPN model before and after deformation

圖13 HIRENASD機翼變形計算結果Fig.13 Deformation calculation results of HIRENASD wing

圖13 給出了HIRENASD 機翼靜氣動彈性計算的變形結果。HIRENASD 機翼是由德國亞琛工業(yè)大學設計,用于研究跨聲速飛行雷諾數(shù)條件下氣動彈性現(xiàn)象與機理的風洞試驗模型,曾被選為第一屆氣動彈性預測研討會(AePW1)的標準算例,詳細信息見文獻[20]。本文算例的計算條件取自文獻[21],用于驗證氣動彈性計算程序的精準度,具體計算參數(shù)為馬赫數(shù)Ma＝0.80、雷諾數(shù)Rec＝1.4×107、載荷因子q/E＝4.7×10－7、迎角α＝3°,其中載荷因子q/E 為風洞來流速壓q 與風洞模型材料彈性模量E 的比值。圖13中LE、TE分別表示機翼前后緣,TRIP表示采用SGDP V1.0的氣動彈性計算模塊,TAU 為德國宇航院的CFD 計算軟件,EXP表示風洞試驗得到的結果?？梢钥闯?采用SGDP V1.0的TRIP 獲得和TAU、EXP比較一致的彎曲和扭轉變形結果,在靠近翼梢位置TRIP 與TAU 計算的扭轉變形存在略微的差異,可能是TAU 僅采用了前4階模態(tài)導致扭轉變形計算截斷誤差增大引起的。

圖14 HIRENASD機翼壓力系數(shù)計算結果Fig.14 Pressure coefficients calculation results of HIRENASD wing

圖14 給出了HIRENASD 機翼3個展向剖面位置的壓力系數(shù)計算結果,圖中Rigid表示機翼剛性條件下獲得的壓力分布,Elastic表示采用SGDP V1.0計算的彈性狀態(tài)下的壓力分布。可以看出,彈性狀態(tài)下的壓力分布與試驗值十分吻合,在靠近翼梢位置與剛性結果差異明顯,側面驗證了SGDP V1.0程序的可靠性,能夠提供高質量的變形網格。

5 結 論

對國家數(shù)值風洞工程NNW 標準共性基礎庫的第1款結構網格變形程序SGDP V1.0的特性進行了簡要的概述,并通過DLR-F6不同規(guī)模網格對SGDP V1.0的性能進行了測試,以及開展了SGDP V1.0網格變形能力和網格質量維持特性驗證、靜氣動彈性數(shù)值模擬等工作。從測試與計算結果可以得到:

1)SGDP V1.0擁有非常高的網格變形效率,可以在數(shù)秒內完成千萬量級結構網格更新,能夠滿足目前絕大多數(shù)工程實際中的網格變形需求。

2)針對不同類型結構網格變形問題,SGDP V1.0均能獲得較高質量的變形網格,具有良好的網格變形處理和質量維持能力,程序展示了優(yōu)良的魯棒性。

3)SGDP V1.0采用了扁平化和模塊化的參數(shù)輸入方式,可以通過少量幾個輸入參數(shù)的修改應用于不同類型的網格變形算例,展示了良好的通用性。

下一步SGDP V1.0將結合使用經驗和反饋意見,對程序不斷地進行迭代完善,并結合未來超大規(guī)模網格使用需求,開展并行版本的開發(fā)。SGDP V1.0將根據單位下一步的知識產權政策,在逐步完善后確定是否開源,如有使用需求,可以通過郵箱y.sun＠cardc.cn與作者聯(lián)系,研究團隊將提供執(zhí)行程序版本。

致謝:感謝中國空氣動力研究與發(fā)展中心洪俊武研究員在多塊網格自適應稀疏/加密方面提供的建議和幫助。