迎風格式在水滴歐拉方程數值求解中的應用

寧義君，顧洪宇，朱東宇，喬龍，李艷亮

1.中國航空工業空氣動力研究院 遼寧省飛行器防除冰重點實驗室， 遼寧 沈陽 110034

2.中國航空工業空氣動力研究院 沈陽市計算流體力學重點實驗室， 遼寧 沈陽 110034

數值模擬是研究飛機結冰的重要手段，其中水滴撞擊特性是結冰預測以及防/除冰系統設計的必備輸入條件。水滴撞擊特性計算即通過求解水滴運動方程，模擬水滴運動軌跡，并確定水滴撞擊物面的區域以及撞擊區域內的水滴收集系數。

對于水滴運動軌跡的數值模擬，最先發展起來的是拉格朗日法[1-2]，國內早期也普遍采用這種方法[3-4]。但模擬復雜三維構型時，該方法釋放粒子數、網格搜索和空間插值計算量都大幅增加，且難以精確捕捉幾何表面邊界并計算水滴撞擊特性[5]。而歐拉法將水滴看作空間連續分布的液態相，以水滴經過的空間位置為研究對象建立水滴運動控制方程，引入水滴場的概念，通過設置物面邊界條件獲取水滴收集系數。相比拉格朗日法，歐拉法在計算復雜外形的水滴撞擊特性時更具明顯優勢。有研究者對比了歐拉法與拉格朗日法的計算結果，表明水滴收集系數分布基本一致[6]。通常認為這兩種方法的計算結果是等價的，且都被飛機適航分析人員所接受[7]。

受益于計算流體力學的發展，歐拉法于20世紀80年代從兩相流的思想發展而來[8]，此后越來越多的研究者開始采用歐拉法計算水滴撞擊特性[9]。當采用歐拉法時，若求解格式不穩定，極易導致非物理解出現[10]。為延緩計算崩潰的出現，很多計算代碼都采用一階格式[11]。Y.Bourgault等[12-13]基于有限單元法實現了水滴歐拉方程二階精度求解，并將該方法引入FENSAP-ICE[14-15]軟件當中（現已被ANSYS 公司收購）。而目前對水滴歐拉方程的求解，絕大多數仍然采用有限體積法[16]。S.K.Jung等[10,17-18]提出了一種水滴歐拉方程改造方法，將其由退化雙曲型轉化為完全雙曲型，并基于有限體積法開發了一種二階保正性迎風離散格式。

國內在這方面早期以FLUENT等商業軟件用戶定義函數研究為主[19]。隨著數值開發工作的深入，開始編制水滴歐拉方程的數值求解代碼。對于對流通量的空間離散，最早出現MUSCL[20-21]、QUICK[22-23]等格式，隨后迎風與線性插值相結合的方法[24]、VanLeer 迎風插值方法[25-26]等被提出。近期趙文朝等[27]采用歐拉法對NACA0012 翼型、冰風洞風道、S 形進氣道及三段翼的水滴撞擊特性進行了網格影響分析，表明當水滴運動沒有受到上游部件的導流或者遮擋時，歐拉法計算結果具有網格無關性。

在對流通量離散格式上，一階格式具備較高的穩定性，但數值耗散大，計算精度較低；而單純的高階格式容易引起數值振蕩和密度脈沖[10]。有研究者通過添加擴散項以提高穩定性，但人為增加了數值耗散[28-29]。一階迎風疊加線性重構與限制器函數的二階離散格式已廣泛應用于空氣歐拉方程的求解[30]，但鮮有在水滴歐拉方程上的應用。

本文發展了一種基于迎風離散格式的具備二階精度的水滴歐拉方程數值求解方法，且無須添加擴散項，并已集成至航空工業氣動院自主開發的航空飛行器結冰仿真軟件穿云（ChuanYun）V1.0 中。通過對NACA0012 翼型水滴撞擊特性的模擬，對一階與二階精度計算結果進行了分析，并根據SAE ARP5903[31]規范性要求，采用若干特征參數對穿云V1.0計算結果與FENSAP-ICE軟件進行量化比較，以驗證本文研究方法的有效性與計算精度，最后以相同計算方法應用于NACA23012翼型和NLF0414翼型。

1 數值方法

1.1 控制方程

將水滴看作分布在空氣中的連續相，引入液態水含量（LWC）表示單位體積空氣中水的質量。水滴在空氣中的LWC 通常在10-3量級，因此在建立控制方程時一般忽略水滴對空氣的作用，采用松耦合方式進行求解，即先求解空氣場，然后基于空氣場再求解水滴運動[32]。

為建立水滴運動的數學模型，對其物理過程進行適當的假設[8,32]：（1）水滴尺寸足夠小，一般在微米（μm）量級，在空氣中均勻分布；（2）水滴為球體，運動過程中不變形、不破碎、不碰撞或凝聚；（3）水滴在運動過程中不發生熱交換、不蒸發，水滴溫度、密度等物性參數不變；（4）作用于水滴的外力只有重力、空氣阻力和浮力；（5）水滴與物面撞擊時不反彈、不飛濺。

基于以上假設，根據質量守恒和動量守恒定律，水滴控制方程可表述為[32]

式中，ρ為水滴在空氣中的LWC；ρa為空氣密度；ρw為水密度；u為水滴速度矢量；ua為空氣速度矢量，K為慣性因子；g為重力加速度矢量。ρa和ua都由空氣場計算結果中直接提取。

慣性因子K為

式中，μa為空氣動力黏度；d為水滴直徑；f為阻力函數；CD為水滴阻力系數；Red為水滴雷諾數。

水滴雷諾數Red為

在計算水滴阻力時，采用經典圓球阻力系數計算公式[33]

采用有限體積法時，對于控制體Ω，水滴歐拉方程的積分形式為

式中，W為守恒量；F為對流通量；Q為源項；可表述如下

式中，ua，va，wa分別為空氣速度ua的三個分量；u，v，w分別為水滴速度u的三個分量；gx，gy，gz分別為重力加速度g的三個分量。

V為逆變速度，計算公式為

式中,nx，ny，nz分別為面法向矢量n的三個分量。

對于二階精度有限體積法，物理量在控制體內的分布是均勻的，將式（6）半離散化，表達形式為

式中，m為控制體面的序號；NF為控制體的面數；ΔS為控制體面的面積。

1.2 空間離散格式

基于通量矢量分裂（FVS）方法，將通量矢量F分解為非負和非正兩部分，非負部分采用左側值，用L 來表示；非正部分采用右側值，用R來表示

由于水滴控制方程雅可比（Jacobi）系數矩陣有多重特征值V，而V為逆變速度，因此根據V與0 的大小關系定義分裂逆變速度，非負和非正逆變速度有

則通量矢量F的具體分裂方案為

1.3 重構與梯度

迎風格式需要確定控制體面的左側值和右側值。當假設每個控制體內的物理量為常數，并且左右值分別為相應控制體單元值時，決定了這種空間離散只具備一階精度。一階迎風格式雖然具備良好的穩定性和有界性，但會導致過度的數值擴散[30]。

若假設物理量在控制體內是變化的，則左右值插值方式決定了具體的空間精度。對于二階精度，物理量在控制體內是線性的。為了計算左右值，必須對物理量進行重構。

線性重構是目前最流行的重構方法。假設物理量在控制體內是分段線性分布的，任意物理量?左右值重構如下

式中，下標I和J分別為左右單元編號；r為單元體心到面心距離；Ψ為限制器函數；??為物理量?所在單元的梯度。根據格林-高斯（Green-Gauss）定理，任意物理量?的梯度近似為?與向外單位法向量乘積的控制體面積分，計算公式如下

式中，nIJ和SIJ分別為控制體面的單位法矢量和面積；N為與單元I相鄰的單元數量。

1.4 限制器

二階或高階迎風格式在空間離散時需要使用限制器函數，以防止在大梯度區域（如水滴遮蔽區邊界）產生數值振蕩和非物理解。因此，所使用限制器至少能夠保證單調性，即求解域中物理量的極大值不增加，極小值不減少，并且在時間推進過程中不會產生新的局部極值。

使用兩種限制器函數：MinMod[34]和Venkatakrishnan[35-36]，見表1。

表1 限制器函數列表Table 1 List of limiter functions

各限制器函數中，r計算公式為

其中

式中，minJ和maxJ分別為單元I所有相鄰單元J的最小值和最大值，ε為機器精度近似值。

單元限制器取面限制器的最小值

1.5 時間推進格式

將水滴歐拉方程的半離散式（8）簡化為

式中，R為水滴歐拉方程的殘差。

采用多步Runge-Kutta時間推進方法

式中，α1, …,αm為迭代系數；Δt為時間步長。

為加速收斂，采用當地時間步長，計算公式為

式中，CFL為庫朗數。

1.6 邊界條件

對于初始條件ρ∞= LWC∞，u∞=ua,∞，∞代表自由來流值。對于遠場邊界條件，如果是流入方向，邊界值等于初始條件；如果是流出方向，邊界值等于相鄰單元值。

對于對稱邊界條件，速度無穿透，其值為ub-(ubn)n，ub為相鄰單元速度矢量；其他物理量邊界值等于相鄰單元值。

對于物面邊界條件，根據物面第一層網格水滴速度物面法向分量ubn，若小于0，代表指向物面內，邊界值采用相鄰單元值；若大于0，代表指向物面外，邊界ρ值取接近于0的小量（如LWC∞×10-7），并代入相應的守恒量ρu中。

1.7 水滴收集系數

水滴撞擊特性由水滴收集系數來表征，水滴收集系數計算公式為

式中，U∞為自由來流速度；Un為水滴物面法向速度。

2 計算結果與分析

空氣流場基于航空工業氣動院自主開發的航空飛行器氣動力仿真軟件飛廉（FeiLian）V1.0[37]進行計算，采用MENTER_SST 湍流模型、Roe 對流通量離散格式以及Venkatakrishnan限制器。

水滴撞擊特性基于航空工業氣動院自主開發的航空飛行器結冰仿真軟件穿云（ChuanYun）V1.0進行計算，計算結果將與FENSAP-ICE軟件進行對比。FENSAP-ICE軟件已普遍應用于飛行器防除冰系統設計，具備工程應用可信度，與其對比可以驗證本文研究方法的有效性與準確性。

根據SAE ARP5903[31]規范性要求，對比計算結果時采用的水滴收集系數特征參數有：βmax為水滴收集系數最大值誤差，(βmax,C-βmax,F)/βmax,F；Sβmax為水滴收集系數最大值所在位置誤差，(Sβmax,C-Sβmax,F)/Chord；Sup為上翼面水滴撞擊極限誤差，(Sup,C-Sup,F)/Chord；Sdw為下翼面水滴撞擊極限誤差，(Sdw,C-Sdw,F)/Chord，其中下標C 為穿云軟件，下標F 代表FENSAP-ICE軟件，Chord代表翼型弦長。

2.1 NACA0012翼型

NACA0012 翼型弦長為0.3048m，計算狀態見表2，主要從不同限制器函數和不同狀態參數兩個方面進行對比分析。

表2 NACA0012翼型計算狀態Table 2 Calculation condition for NACA0012 airfoil

采用一階迎風和表1所列兩種限制器計算的水滴收集系數與FENSAP-ICE對比如圖1所示，圖1中橫坐標為量綱一化弧長S/c（c為翼型弦長），以翼型前緣點為原點，特征參數對見表3。從表3中可以看出，一階迎風和兩種限制器計算的水滴撞擊極限與FENSAP-ICE 一致，誤差主要集中在駐點區域。一階迎風在駐點處的水滴收集系數明顯高出其他計算結果，水滴收集系數最大值誤差βmax達到6.38%，且出現明顯的數據波動，如圖1 所示。一階迎風駐點附近LWC 云圖如圖2 所示，LWC 在駐點附近出現集中，根據水滴收集系數計算公式（19），該現象導致了駐點處水滴收集系數突然增大。

圖1 不同限制器水滴收集系數對比Fig.1 Comparison between collection coefficient for different limiters

圖2 一階迎風LWC云圖Fig.2 LWC cloud of first order upwind

表3 不同限制器水滴收集系數特征參數對比Table 3 Comparison between collection coefficient characteristic parameters for different limiters

在兩種限制器中，MinMod和Venkatakrishnan計算結果與FENSAP-ICE 基本一致，水滴收集系數最大值βmax誤差分別為0.33%和0.34%。說明與一階迎風相比，MinMod 和Venkatakrishnan兩種限制器對水滴收集系數的改進作用是一致的。

與一階迎風相比，使用線性重構與限制器的主要作用在于對LWC 空間分布間斷面進行精確捕捉。一階迎風和兩種限制器、FENSAP-ICE計算的翼型尾跡附近LWC云圖如圖3所示。從圖3可以看出，一階迎風存在明顯的過度耗散，其水滴遮蔽區尺寸明顯小于其他限制器和FENSAPICE 計算結果。從翼型尾跡水滴遮蔽區的寬度來看，兩種限制器計算的水滴遮蔽區比FENSAP-ICE 計算結果略??；從水滴遮蔽區邊界附近的LWC強度來看，兩種限制器計算的LWC峰值比FENSAP-ICE計算結果有所減小。

圖3 不同限制器LWC云圖對比Fig.3 Comparison between LWC cloud of different limiters

根據以上分析，MinMod和Venkatakrishnan兩種限制器對水滴遮蔽區的捕捉精度明顯高于一階迎風，但水滴遮蔽區空間分布比FENSAP-ICE 略小，說明在數值耗散上比FENSAP-ICE 略大。并且兩種限制器計算的水滴遮蔽區基本無差別，對水滴遮蔽區的改進作用是一致的。針對MinMod和Venkatakrishnan兩種限制器，以表1狀態為基準，不同水滴直徑計算結果對比如圖4所示，不同風速計算結果對比如圖5所示，特征參數對比分別見表4和表5。從圖4和圖5、表4和表5可以看出，兩種限制器計算的水滴收集系數是一致的，與FENSAP-ICE 軟件相比，水滴收集系數最大值誤差βmax不超過0.73%。這說明MinMod 和Venkatak rishnan兩種限制器對水滴收集系數的計算不受水滴直徑和風速等計算狀態的影響。

圖4 不同水滴直徑水滴收集系數對比Fig.4 Comparison between collection coefficient for different droplet diameters

圖5 不同風速水滴收集系數對比Fig.5 Comparison between collection coefficient for different velocities

表4 不同水滴直徑水滴收集系數特征參數對比Table 4 Comparison between collection coefficient characteristic parameters for different droplet diameters

表5 不同風速水滴收集系數特征參數對比Table 5 Comparison between collection coefficient characteristic parameters for different velocities

2.2 NACA23012和NLF0414翼型

根據NACA0012 翼型計算結果，可以看出MinMod 和Venkatakrishnan兩種限制器計算結果彼此之間幾乎是一致的，與FENSAP-ICE 軟件計算結果也基本一致。將相同計算方法應用于其他翼型以進一步驗證該方法的計算精度。

采用NACA23012翼型，弦長為0.4572m，基于表2計算狀態，MinMod和Venkatakrishnan兩種限制器計算的水滴收集系數與FENSAP-ICE 對比如圖6 所示，特征參數對比見表6，水滴收集系數最大值βmax誤差分別為0.86% 和0.63%。

圖6 NACA23012翼型水滴收集系數對比Fig.6 Comparison between collection coefficient for NACA23012 airfoil

表6 NACA23012翼型水滴收集系數特征參數對比Table 6 Comparison between collection coefficient characteristic parameters for NACA23012 airfoil

采用NLF0414 翼型，弦長為0.9144m，基于表2 計算狀態，MinMod和Venkatakrishnan兩種限制器計算的水滴收集系數與FENSAP-ICE 對比如圖7 所示，特征參數對比見表7，水滴收集系數最大值βmax誤差分別為-2.35%和-2.29%。

圖7 NLF0414翼型水滴收集系數對比Fig.7 Comparison between collection coefficient for NLF0414 airfoil

表7 NLF0414翼型水滴收集系數特征參數對比Table 7 Comparison between collection coefficient characteristic parameters for NLF0414 airfoil

3 結束語

本文將一種基于逆變速度分裂的迎風格式應用于水滴歐拉方程對流通量的空間離散，通過計算與分析，可以得到如下結論：

（1） 該迎風格式只具備一階精度，在求解水滴歐拉方程時，駐點處的水滴收集系數出現較大誤差，通過疊加線性重構與限制器函數，對水滴歐拉方程的求解可以達到二階精度。

（2） MinMod 和Venkatakrishnan 限制器計算結果彼此之間基本沒有差別，與FENSAP-ICE軟件相比，水滴收集系數基本一致。

（3） 由于一階迎風格式的數值耗散較大，水滴遮蔽區空間分布遠小于FENSAP-ICE 軟件，通過上述二階求解方法，水滴遮蔽區計算精度大幅提升，已接近于FENSAP-ICE軟件。