任靖豪，王強(qiáng)，2，李維浩，劉宇，易賢，2，*

1．中國空氣動力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

2．中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

飛機(jī)結(jié)冰數(shù)值計算中，水滴收集率的計算方法分為兩類：歐拉方法［1-3］和拉格朗日方法。拉格朗日方法通過求解水滴運(yùn)動方程，分析單顆水滴的運(yùn)動軌跡，進(jìn)而獲取液滴的撞擊特性。由于其算法簡單，計算過程表征直觀，且在描述過冷大水滴（SLD）動力學(xué)特性［4-5］方面具有天然優(yōu)勢，被廣泛應(yīng)用于結(jié)冰數(shù)值仿真［6-8］。在二維和簡單三維外形問題上，展現(xiàn)出了高效的計算優(yōu)勢。然而，針對三維復(fù)雜構(gòu)型問題，由于水滴軌跡計算量大、收集率計算復(fù)雜，制約了拉格朗日方法的發(fā)展。

近些年，國內(nèi)外的研究者也在不斷挖掘拉格朗日方法的計算潛力。其中Lee 等［9］利用高斯形函數(shù)發(fā)展了一種基于粒子統(tǒng)計的數(shù)值計算方法，并將該方法應(yīng)用于發(fā)動機(jī)葉片的三維水滴撞擊特性研究中。文獻(xiàn)［10］在此基礎(chǔ)上結(jié)合離散隨機(jī)步模型，實(shí)現(xiàn)了考慮湍流影響的水滴運(yùn)動過程模擬。為了進(jìn)一步適應(yīng)較大規(guī)模的水滴軌跡跟蹤計算，Rendall 等［11］發(fā)展了一種水滴運(yùn)動方程的隱式計算方法，顯著提高了拉格朗日水滴軌跡方程的求解效率。文獻(xiàn)［12-13］基于拉格朗日流管法，利用網(wǎng)格加密的方式，自適應(yīng)地控制關(guān)鍵區(qū)域的水滴釋放密度。相比統(tǒng)計法，該方法極大地降低了水滴軌跡計算規(guī)模。文獻(xiàn)［13］還提出了一種限制徑向基函數(shù)插值方法，用于改善撞擊點(diǎn)向物面網(wǎng)格點(diǎn)傳遞水滴收集率的魯棒性。筆者團(tuán)隊(duì)在此基礎(chǔ)上發(fā)展了一種壁面投影算法［14］，改進(jìn)了拉格朗日軌跡法對迎風(fēng)區(qū)和背風(fēng)區(qū)的辨識能力，實(shí)現(xiàn)了三維復(fù)雜構(gòu)型的水滴收集率快速計算。但是上述拉格朗日流管法需要依賴插值技術(shù)將水滴收集率傳遞到網(wǎng)格角點(diǎn)。其中插值精度很大程度上受撞擊點(diǎn)分布特性的影響。當(dāng)撞擊點(diǎn)分布不均勻或樣本較少時，計算結(jié)果的可靠性會隨之下降。

為了克服插值帶來的不利影響，本文提出了一種基于梯度下降模型［15］的水滴收集率計算方法。文中采用該方法開展了三維圓球、30P30N 平直多段翼、ONERA-M6 機(jī)翼和某型發(fā)動機(jī)短艙等典型外形的水滴收集率計算，并分別與文獻(xiàn)試驗(yàn)數(shù)據(jù)和計算結(jié)果進(jìn)行對比分析，用于驗(yàn)證改進(jìn)方法的可行性。

1 液滴撞擊特性計算方法

1.1 流場計算

水滴撞擊特性計算分為兩步：①首先采用歐拉方法求解空氣流場；②基于空氣流場結(jié)果開展水滴軌跡計算，進(jìn)而獲得水滴收集率。本文采用中國空氣動力研究與發(fā)展中心自研的PHengLei求解器［16］開展流場計算。

1.2 水滴運(yùn)動方程及其求解方法

求解的水滴軌跡方程滿足以下幾個假設(shè)條件：

1） 水滴運(yùn)動過程中僅受自身重力、流場的黏性阻力和浮力，且不考慮液滴對空氣流場的影響。

2） 采用球形阻力公式計算液滴所受阻力，且期間不考慮液滴變形、破裂、凝聚等動力學(xué)特性。

3） 液滴運(yùn)動過程中不考慮質(zhì)量損失。

根據(jù)牛頓第二定律，液滴運(yùn)動的控制方程表達(dá)式為

式中：ua和ud分別表示當(dāng)?shù)乇尘傲鲌鏊俣群鸵旱蔚倪\(yùn)動速度；ρd和ρa(bǔ)分別為液滴和空氣的密度；g為重力加速度；f為液滴的弛豫系數(shù)，其表達(dá)式為

其 中：μ為空 氣黏性系數(shù)；d為液滴 粒徑；CD為 液滴阻力系數(shù)；Re為液滴運(yùn)動的相對雷諾數(shù)；組合參數(shù)CD·Re/24 的表達(dá)式為

式（1）為求解一階常微分方程的初值問題。通過給定水滴初始位置和速度條件，采用四階Runge-Kutta 方法［17］對式（1）進(jìn)行數(shù)值積分，求得水滴的運(yùn)動軌跡。

式中：r為水滴位置坐標(biāo)；下標(biāo)0 表示水滴初始參數(shù)；rp為水滴顆粒空間坐標(biāo)。

1.3 液滴定位及流場參數(shù)插值

求解水滴運(yùn)動方程過程中，每個時間步都需要更新水滴所受的氣動力。因此需要明確各時刻水滴顆粒所處網(wǎng)格單元的編號，通過插值得到液滴當(dāng)?shù)匚恢玫目諝饬鲌鲂畔ⅰ?/p>

采用定向查找算法［14］，以最優(yōu)路徑快速定位出液滴所處單元。明確了水滴所處網(wǎng)格單元后，采用徑向基函數(shù)（RBF）插值方法［18］計算水滴所處位置的流場參數(shù)值。該插值方法的優(yōu)勢在于插值精度高，且外插結(jié)果具有較好的連續(xù)性。

徑向基函數(shù)的基本表達(dá)式為

式中：n為樣本點(diǎn)個數(shù)；rtar和ri分別表示插值目標(biāo)點(diǎn)和第i個樣本點(diǎn)在變量空間中的坐標(biāo)；ωi為第i個樣本點(diǎn)的插值權(quán)重系數(shù)；φ（r）為徑向基函數(shù)；R為基函數(shù)的支撐半徑。選用Wendlend’s C2 函數(shù)作為基函數(shù)，表達(dá)式為

將粒子當(dāng)前所處的網(wǎng)格單元角點(diǎn)作為插值樣本S=｛r1，r2，…，rn｝，Sv=｛q1，q2，…，qn｝，通過式（8）求解各樣本點(diǎn)的插值權(quán)重。

得到權(quán)重系數(shù)矩陣W后，代入式（5）估算出當(dāng)前位置的流場信息。

如圖1 所示，在近壁面長寬比較大的網(wǎng)格內(nèi)進(jìn)行插值時，數(shù)值誤差易造成徑向基函數(shù)矩陣為奇異矩陣。因此，需要將樣本點(diǎn)變量空間坐標(biāo)進(jìn)行無量綱化，消除變量量綱對插值結(jié)果的影響。

圖1 無量綱化示意圖Fig. 1 Dimensionless processing

1.4 水滴收集率計算

液滴收集率的定義：單位時間內(nèi)，壁面與自由來流處單位面積上水滴的質(zhì)量通量的比值，表達(dá)式為

式中：為單位時間流過面元A的液態(tài)水質(zhì)量。當(dāng)液滴軌跡不發(fā)生交叉時，可用多條相鄰軌跡組成的“流管”描述該區(qū)域內(nèi)所有水滴的運(yùn)動特性。式（10）可以改寫成如下形式：

式中：n為流管內(nèi)的水滴個數(shù)；md，i為編號i的液滴質(zhì)量。若水滴發(fā)生交叉，并在物面上表現(xiàn)出不同流管的撞擊域發(fā)生重疊，采用面積比率的計算方式易低估水滴收集率。此時，通常采用時均統(tǒng)計的方法，對撞擊特性進(jìn)行分析。

2 基于梯度下降模型的收集率計算

當(dāng)水滴撞擊點(diǎn)與網(wǎng)格角點(diǎn)一致時，可直接計算物面網(wǎng)格面心上的水滴收集率，計算結(jié)果不受插值精度影響?；谶@一思路，建立水滴撞擊準(zhǔn)確度的評價函數(shù)F。將水滴釋放點(diǎn)坐標(biāo)作為函數(shù)輸入，撞擊點(diǎn)到目標(biāo)角點(diǎn)距離作為函數(shù)輸出。采用梯度下降法沿函數(shù)梯度下降方向求解極小值點(diǎn)，即可得到一條逼近目標(biāo)角點(diǎn)的水滴軌跡。

2.1 水滴軌跡的鄰域模型

引入方形鄰域模型，通過在水滴釋放點(diǎn)P鄰近區(qū)域增加額外的水滴軌跡信息輔助求解函數(shù)的偏導(dǎo)數(shù)。圖2中δ表示鄰域尺度，黃色單元格為輔助水滴的釋放位置，其中輔助水滴粒徑可能大于鄰域尺度。在獲取模型各水滴撞擊點(diǎn)與目標(biāo)點(diǎn)的距離值后，通過差分或最小二乘法計算P點(diǎn)處函數(shù)梯度值。

圖2 水滴軌跡鄰域示意圖Fig. 2 Neighborhood model of droplet trajectory

為保證背景流場參數(shù)插值結(jié)果有較好的連續(xù)性，該鄰域模型內(nèi)的5 顆水滴質(zhì)點(diǎn)需采用相同的插值樣本，即對這5 顆水滴進(jìn)行同步追蹤并默認(rèn)其位于同一網(wǎng)格單元內(nèi)。

2.2 計算方法

已知有一條水滴軌跡trai撞壁，其撞壁信息包含：撞擊點(diǎn)位置、釋放點(diǎn)位置、撞擊點(diǎn)所處的壁面網(wǎng)格單元編號。選取trai撞擊的壁面網(wǎng)格單元上某角點(diǎn)作為目標(biāo)點(diǎn)。將trai釋放點(diǎn)坐標(biāo)xi作為函數(shù)F的輸入，trai撞擊點(diǎn)到目標(biāo)點(diǎn)的距離si作為函數(shù)F的輸出，即

由于在拉格朗日算法中水滴通常釋放于垂直來流的某一空間平面上，即xi是一組平面正交基表示的二維矢量，其矩陣形式可寫成xi［u，v］。

采用2.1 節(jié)介紹的水滴鄰域模型，可以額外獲得4 個撞擊點(diǎn)信息。通過最小二乘法求解方程組（13）：

式中：

兩特征方向的偏導(dǎo)數(shù)可表示為

式（14）中Δulm=ul-um，鄰域模型中心點(diǎn)編號為0，其他點(diǎn)編號分別為1，2，3，4。

求得梯度后，采用遞歸式（16），調(diào)整水滴釋放點(diǎn)位置重新計算水滴軌跡，并帶回函數(shù)（13）重復(fù)上述步驟直至函數(shù)值小于給定閾值，式（16）中g(shù)（xik）表示水滴釋放位置修正的位移量，計算示意圖如圖3 所示。

圖3 逼近目標(biāo)點(diǎn)過程示意圖Fig. 3 Demonstration of process of closing to target impact point

為了加快收斂，消除目標(biāo)點(diǎn)附近迭代振蕩，采用Nesterov 動量梯度下降方法，通過引入下一步的梯度方向修正當(dāng)前時間步的前進(jìn)方向，表達(dá)式為

方程（17）中的λk，i和ηk，i分別通過預(yù)估校正的方式獲得，其中預(yù)估值的表達(dá)式為

若λk，i和ηk，i的 預(yù) 估 值 過 大，可 能 導(dǎo) 致 液 滴 撞擊點(diǎn)出現(xiàn)遠(yuǎn)離目標(biāo)點(diǎn)的趨勢，或是直接飛離物面。因此根據(jù)式（20）和式（21），按先后順序?qū)ζ溥M(jìn)行校正。

在校正過程中，釋放位置更新后的液滴若出現(xiàn)上述情況，則減小λk，i和ηk，i，直至更新后的距離函數(shù)值小于更新前。

在已知撞擊某一網(wǎng)格角點(diǎn)的水滴軌跡后，可以通過類似的技巧逼近得到相鄰角點(diǎn)的撞擊軌跡，如圖4 所示。在逐點(diǎn)逼近過程中，若推進(jìn)步長小于特定值后，水滴仍飛出壁面或是撞擊在非目標(biāo)點(diǎn)鄰近的面網(wǎng)格單元上，將該角點(diǎn)標(biāo)記為撞擊極限點(diǎn)，此時不再查找該點(diǎn)相鄰角點(diǎn)的撞擊軌跡。該算法的詳細(xì)計算流程見圖5。

圖4 鄰點(diǎn)檢索算法示意圖Fig. 4 Diagram of neighbor search method

圖5 計算流程圖Fig. 5 Calculation flow chart

根據(jù)物面網(wǎng)格單元各角點(diǎn)所對應(yīng)水滴軌跡組成的“流管”，通過式（11）計算得到單元中心的水滴收集率。

3 計算結(jié)果及分析

分別采用三維圓球、30P30N 平直多段翼、ONERAM6 機(jī)翼以及某型短艙構(gòu)型對上述方法的可行性進(jìn)行驗(yàn)證。

3.1 三維圓球水滴收集率

計算工況如下：遠(yuǎn)場來流速度為75 m/s，圓球直徑為15.04 cm，遠(yuǎn)場靜壓95 840 Pa，，液滴粒徑為18.6 μm，采用LangmuirD 多尺度粒徑分布模型。迭代逼近過程中的偏移誤差設(shè)為1.0 ×10-5m，時間積分步長5.0×10-5s，水滴釋放平面距模型前緣1.2 m 處，即當(dāng)水滴軌跡撞擊點(diǎn)與目標(biāo)角點(diǎn)的空間距離小于1.0 ×10-5m 時，認(rèn)為該軌跡的撞擊點(diǎn)與目標(biāo)角點(diǎn)擬合。

圖6 中給出了本文方法計算得到的撞擊點(diǎn)分布結(jié)果，并與文獻(xiàn)［14］所述方法在撞擊極限附近加密3 次的撞擊水滴陣列進(jìn)行了比較。圖中顯示本文計算的撞擊點(diǎn)位置與網(wǎng)格角點(diǎn)的重合度良好，且預(yù)測的撞擊極限與3 次加密后的撞擊網(wǎng)格結(jié)果基本一致。表明本文方法能夠準(zhǔn)確預(yù)測得到物面水滴撞擊域以及撞擊極限。

圖6 圓球撞擊點(diǎn)分布結(jié)果Fig. 6 Distribution of impinging point of sphere case

圖7為不同計算方法與試驗(yàn)數(shù)據(jù)［19］的水滴收集率β曲線對比結(jié)果。其中，本文計算得到的收集率曲線與試驗(yàn)數(shù)據(jù)除在駐點(diǎn)處的收集率峰值仍存在一定差異外，兩者分布規(guī)律基本吻合。通過局部放大圖可以觀察到，本文改進(jìn)后的計算結(jié)果與FENSAP 計算結(jié)果貼合度更好，并且在駐點(diǎn)處的水滴收集率更接近試驗(yàn)數(shù)據(jù)。圖8 比較了依賴插值的拉格朗日流管法和本文方法計算得到的水滴收集率云圖。顯而易見，本文方法的計算結(jié)果較傳統(tǒng)方法的插值結(jié)果表現(xiàn)出更好的軸對稱特征，更加符合實(shí)際物理分布規(guī)律。上述結(jié)果說明本文方法能夠克服現(xiàn)有拉格朗日方法中插值對水滴收集率計算結(jié)果帶來的不利影響，改善了方法的魯棒性。

圖7 圓球水滴收集率分布曲線比較Fig. 7 Comparison of computational and experimental collection efficiency of sphere case

圖8 圓球水滴收集率分布云圖Fig. 8 Contour of collection efficiency of sphere case

本文方法需要通過迭代運(yùn)算獲得壁面角點(diǎn)的撞擊軌跡，而每次迭代都需要重新計算水滴軌跡。因此，相比傳統(tǒng)流管方法，本文方法需要耗費(fèi)更多的計算機(jī)時。圖9 統(tǒng)計了不同軌跡運(yùn)算次數(shù)的目標(biāo)點(diǎn)個數(shù)分布情況，橫軸表示不同的軌跡運(yùn)算次數(shù)區(qū)間，縱軸表示與之對應(yīng)的目標(biāo)點(diǎn)個數(shù)。結(jié)果顯示大多數(shù)目標(biāo)點(diǎn)的軌跡運(yùn)算次數(shù)小于20。相比在提升計算精度和改善算法魯棒性上等方面帶來的收益，本文方法增加的計算量是可以接受的。另外，可以采用多起點(diǎn)并行計算策略，進(jìn)一步縮短遍歷撞擊域內(nèi)所有角點(diǎn)所需的時間。本節(jié)算例采用10 線程并行計算耗費(fèi)機(jī)時33.86 s，與單核計算耗時343.82 s相比有較顯著的提升。

圖9 撞擊域內(nèi)水滴軌跡迭代運(yùn)算次數(shù)統(tǒng)計Fig. 9 Statistics of iteration times in impact zone

為了克服因部分角點(diǎn)迭代次數(shù)增加導(dǎo)致計算效率下降的困難，下一步工作將對不同的軌跡定位插值方法、鄰域模型計算精度以及迭代算法進(jìn)行研究分析，制定一套能夠降低迭代運(yùn)算量的控制策略。

3.2 30P30N 平直多段翼水滴收集率

計算工況如下：模型弦長111.44 mm，來流速度78 m/s，遠(yuǎn)場靜壓95 840 Pa，迎角0°，液滴粒徑21 μm，采用LangmuirD 多尺度粒徑分布模型。

圖10 給出了30P30N 平直多段翼表面多連通域的水滴收集率分布云圖。圖中顯示，水滴撞擊域主要位于縫翼前緣和襟翼下翼面，而主翼段沒有出現(xiàn)水滴撞擊。通過截取模型中截面的水滴收集率分布曲線，與Bidwell［20］報告中公開發(fā)布的試驗(yàn)測量結(jié)果進(jìn)行比較，驗(yàn)證計算結(jié)果的準(zhǔn)確性。如圖11 和圖12 所示，計算得到的縫翼處的水滴收集率分布規(guī)律與試驗(yàn)數(shù)據(jù)吻合度較高，襟翼上的水滴收集率峰值雖存在一定差異，但整體趨勢基本相同。上述結(jié)果表明，在復(fù)雜氣動特性影響下，本文方法能夠保證較高的計算精度。

圖10 30P30N 平直多段翼水滴收集率計算結(jié)果Fig. 10 Results of collection efficiency of 30P30N multi-elements airfoil

圖11 縫翼水滴收集率分布曲線Fig. 11 Comparison of computational and experimental collection efficiency on the slat

圖12 襟翼水滴收集率分布曲線Fig. 12 Comparison of computational and experimental collection efficiency on the flap

3.3 ONERA-M6 機(jī)翼水滴收集率

計算工況如下：平均氣動弦長為0.54 m，展長為1.0 m，來流 速 度52 m/s，迎角6°，遠(yuǎn) 場 靜 壓101 325 Pa，液滴粒徑20.0 μm。

圖13 給出了ONERA-M6 機(jī)翼水滴收集率計算云圖。圖14 截取了機(jī)翼沿展向站位在20%，60%和90%處的收集率分布曲線，并與文獻(xiàn)［21-22］的計算結(jié)果進(jìn)行了對比。結(jié)果顯示本文計算結(jié)果與文獻(xiàn)數(shù)據(jù)保持了良好的一致性。

圖13 ONERA-M6 機(jī)翼水滴收集率計算結(jié)果Fig. 13 Results of collection efficiency of ONERA-M6 wing

圖14 沿展向不同站位處S1-S3水滴收集率曲線Fig. 14 Comparison of collection efficiency on Sections S1-S3

圖15給出了不同迎角下的水滴收集率計算結(jié)果?？梢?，隨著迎角的增大，除了上下極限位置發(fā)生改變外，機(jī)翼展向的收集率分布規(guī)律也有明顯變化，翼梢處的水滴收集率變化受迎角影響更加敏感。

圖15 不同迎角下的水滴收集率分布云圖Fig. 15 Contour of collection efficiency with different angles of attack

3.4 某型發(fā)動機(jī)短艙及整流罩水滴收集率

本節(jié)算例是一個典型的旋成體模型，采用非結(jié)構(gòu)類型的流場網(wǎng)格，如圖16 所示。計算工況如下：來 流 速 度67 m/s，迎 角0°，遠(yuǎn) 場 靜 壓101 325 Pa，水滴粒徑44.5 μm。

圖16 短艙流場網(wǎng)格展示Fig. 16 Flow field grid of engine nacelle model

圖17比較了本文方法與FENSAP 計算得到的物面水滴收集率云圖。從圖中可見：①兩者計算結(jié)果一致性較好，預(yù)測得到的撞擊域以及撞擊極限位置基本相同，且水滴撞擊區(qū)域主要位于短艙唇口與帽罩前部段；②本文方法結(jié)果表現(xiàn)出了良好的對稱分布特征，且等高線過渡光滑，撞擊極限刻畫清晰。

圖17 短艙水滴收集率云圖Fig. 17 Contour of collection efficiency of engine nacelle model

圖18和圖19 進(jìn)一步比較了垂直于Y軸和Z軸中截面上的收集率曲線。結(jié)果顯示，兩者預(yù)測的撞擊極限位置以及曲線分布趨勢基本一致，僅在帽罩和短艙駐點(diǎn)區(qū)域的收集率峰值存在一定差異。

圖18 垂直Y 軸中截面水滴收集率曲線對比Fig. 18 Comparison of collection efficiency on the section vertical Y-axis

圖19 垂直Z 軸中截面水滴收集率曲線對比Fig. 19 Comparison of collection efficiency on the section vertical Z-axis

4 結(jié) 論

1） 基于梯度下降法發(fā)展了一種水滴收集率計算方法。該方法通過自適應(yīng)地調(diào)控釋放點(diǎn)位置，以匹配撞擊點(diǎn)和網(wǎng)格節(jié)點(diǎn)，避免了插值誤差對計算結(jié)果帶來的不利影響，使物面水滴收集率表征結(jié)果更加光滑。

2） 計算結(jié)果與文獻(xiàn)數(shù)據(jù)吻合良好，并且該方法在三維平直多段翼、變迎角后掠翼及發(fā)動機(jī)短艙的水滴收集率計算中表現(xiàn)出良好的魯棒性，具備模擬三維復(fù)雜外形水滴撞擊特性的能力。

3） 本文方法需要依賴已知的一條或多條水滴撞擊軌跡開展后續(xù)計算。為了進(jìn)一步完善本方法的計算流程，提高程序的通用性能，下一步工作，將利用水滴軌跡的最小壁面距離信息，實(shí)現(xiàn)初始撞擊軌跡的快速預(yù)測。