王昭力，曾濤，周志宏,*，陳宇

1. 四川大學 建筑與環(huán)境學院，成都 610065

2. 宜賓四川大學產業(yè)技術研究院，宜賓 644000

飛機穿過含有過冷水滴的云層時，水滴撞擊到飛機表面可能凝結成冰，結冰會改變飛機升力部件的外形并影響其氣動性能，導致飛機的操縱性和穩(wěn)定性惡化；其他部位結冰也可能導致傳感器失靈、操縱面卡死、發(fā)動機失火等不良影響，嚴重時直接導致空難事故[1-2]。飛機結冰研究主要有飛行試驗、風洞試驗和數值模擬3種研究手段，其中數值模擬具有成本低、效率高等優(yōu)點，在飛行器防除冰系統設計、適航取證、空難事故分析等階段被廣泛應用。飛機結冰是水滴在空氣繞流作用下，撞擊到飛機表面發(fā)生凝結的現象，這一物理過程的數值模擬主要包含水滴軌跡及撞擊到物體表面的收集率、傳熱相變及水膜溢流、冰形增長3個主要的計算模塊，其中，水滴收集率計算是整個結冰計算的基礎，是其他模塊的輸入依據。同時，水滴收集率也是防冰系統設計的重要輸入依據。

水滴收集率計算主要有基于水滴粒子的拉格朗日方法和基于兩相流觀點的歐拉方法。最先發(fā)展起來的水滴撞擊特性的計算方法是拉格朗日軌跡法，在遠場釋放一定數量水滴粒子，通過求解水滴受力方程，獲得這些水滴的運動軌跡，通過統計撞在物面上的水滴數量或撞擊面積比等方式，獲得飛機表面的水滴收集率，國內楊倩[3]、易賢[4]等在這些方面做了許多工作。

水滴軌跡計算的歐拉兩相流法是將空氣和水滴視為兩相流，水滴相分布在空氣相中，通過在全流域內求解水滴相的控制方程得到水含量的分布以及水滴的速度，再依據物面處的水含量及法向速度計算水滴收集率，相比拉格朗日法，歐拉法在處理復雜外形水滴撞擊特性計算時更有優(yōu)勢。歐拉兩相流計算方法的問題在于，控制方程的高度非線性容易帶來數值振蕩。Bourgault等[5]建立了兩相流動的控制方程并使用有限單元法進行求解，后又在此基礎上添加了間斷捕捉項來解決歐拉法數值不穩(wěn)定的問題，并將此算法植入到FENSAP中[6]。Iuliano等[7]對多段翼的水滴流場進行了計算，并對比了歐拉法與拉格朗日法之間的區(qū)別以及空氣黏性對水滴撞擊特性的影響。Tong等[8]分析了歐拉法數值不穩(wěn)定的原因，在原始的水滴控制方程基礎上添加了擴散項用來消除迭代過程中的奇異性，該擴散項能根據網格粗細和水含量梯度自動調整。國內張大林等[9]在歐拉法中對流項使用MUSCL格式計算了二維NACA0012翼型的水滴收集率，林貴平等[10]在水滴的質量守恒方程中添加自適應的擴散項來提高使用對流項高階離散格式時的穩(wěn)定性，易賢等[11]使用歐拉法發(fā)展了適合于三維、復雜構型的水滴收集率計算程序，其對流項的離散是采用迎風插值和線性插值相結合的方法。招啟軍等[12]為解決尾流等區(qū)域的密度脈沖現象引起的穩(wěn)定性和收斂性問題，提出并建立了遮蔽區(qū)擴散模型。陳海昕等基于歐拉法，發(fā)展了一套結冰模擬軟件[13]。

由于一階格式的精度不夠，而單純的高階格式又容易引起數值振蕩和密度脈沖，目前針對空氣流場對流離散格式的研究，已經有較多的研究[14-17]，但絕大部分未應用于水滴流場的數值求解中。水滴流場和空氣流場的求解有所不同，水滴流場的控制方程相比常規(guī)的空氣流場沒有黏性項，因此其對數值振蕩或者奇異值更加敏感；此外，針對空氣流場的對流項離散格式的研究大多是集中于超高速工況，由于超高速會使空氣流場產生激波，但激波前后流動參數形成的間斷性質與水滴流場中遮蔽區(qū)的界面間斷性質不同。因此，研究水滴流場的對流項離散格式是有意義的。

本文在不添加擴散項的前提下，通過在二階精度格式中引入系列通量限制器來解決收斂性的問題，以發(fā)展飛機結冰計算中水滴流場的高階精度計算格式。

1 數值方法

1.1 水滴流場控制方程

針對飛機結冰過程中水滴-空氣兩相流歐拉法的建立，由于云霧中水滴含量很小，通常有如下假設：

1) 將水滴視為是分布在空氣中的連續(xù)相，采用液態(tài)水含量(Liquid Water Content, LWC)表示單位體積中的水質量，單位為g/m3。引入水滴體積分數α，定義為控制單元內水滴體積占總體積的比例，自由來流的α可采用如下公式計算：

(1)

式中：ρ為水的密度，運動過程中水滴密度不發(fā)生改變。

2) 常規(guī)水滴條件(水滴平均粒徑<40 μm)下，水滴為球形，單個水滴的尺寸用水滴當量直徑deq描述。

3) 水滴在運動過程中不會發(fā)生變形、融合和破碎，且撞擊物面后不發(fā)生飛濺。

4) 水滴的動量方程中不考慮黏性和壓力項。

5) 由于水滴體積占比很小，往往在10-6級別，因此，可近似假設水滴和空氣是單向耦合的，即只考慮水滴受到的空氣阻力，水滴對空氣流場的影響可以忽略，且水滴與空氣不發(fā)生任何能量或質量的交換。

6) 作用在水滴上的力只有阻力和重力，且阻力是穩(wěn)態(tài)的。

基于上述假設，可以獲得如下微分形式的水滴流場的控制方程：

(2)

(3)

使用中心有限體積方法，對積分形式的水滴控制方程在控制體內進行積分，并使用高斯定理，寫成向量的形式，可得到：

(4)

式中：

(5)

式中：ua和u分別為空氣和水滴的速度向量；F為水滴受到的體積力(只有重力)；u、v、w和ua、va、wa分別為水滴速度向量u和空氣速度向量ua的分量；Fx、Fy、Fz為體積力F的分量；τm為動量反應時間：

(6)

f為阻力函數：

(7)

由于常見的飛機飛行時的雷諾數超出了斯托克斯公式計算水滴黏性阻力的理論允許范圍，因此采用修正的阻力函數[11]：

(8)

Rer為相對雷諾數：

(9)

式中：μa為空氣的動力黏性系數。

1.2 空氣流場計算

歐拉法的求解是將空氣和水滴視為兩相流，且只考慮空氣相對水滴相的單向耦合作用。因此在計算水滴流場前需先計算空氣流場。

本文的空氣流場采用中國空氣動力研究與發(fā)展中心計算空氣動力研究所開發(fā)的NNW-PHengLEI軟件進行計算。計算采用SSTk-ω模型，通量分裂采用Roe格式，對流項通量限制器使用的是Minmod限制器。

1.3 水滴流場控制方程的TVD格式

水滴流場的控制方程是一個帶有源項的雙曲型非線性對流方程，因為不具有壓力梯度和擴散項來及時抹平水含量的高梯度，導致了數值不穩(wěn)定性，因此其求解的關鍵在對流項的離散，需要同時考慮數值解的準確性和穩(wěn)定性。

將水滴流場控制方程的積分形式半離散化：

(10)

(FC·S)i,j+1/2,k-(FC·S)i,j-1/2,k+

(FC·S)i,j,k+1/2-(FC·S)i,j,k-1/2

(11)

式中：S為控制體表面的面積矢量；FC為控制體表面的對流矢量。

FC先使用通量矢量分裂(Flux Vector Splitting, FVS)方法分裂成非負和非正兩組：

FC(Q)=FC+(QL)+FC-(QR)

(12)

式中：QL和QR分別為控制體表面左右兩側的守恒變量[ααuαvαw]T。

由于水滴流場控制方程的Jacobian系數矩陣有4重特征值Vc，Vc為逆變速度，是流場速度與網格線切向矢量的點積，ξ、η、ζ分別為曲線坐標系下的3個方向：

Vc=ukx+vky+wkzk=ξ,η,ζ

(13)

(14)

控制體表面的值是根據控制體上的均值插值得到的，不同的插值方式決定了對流項的離散格式的空間精度和收斂性。

1.3.1 控制體表面的對流通量插值

對于任意的物理量q，一階迎風格式(First-Order Upwinding, FOU)可以寫為

qi+1/2=qi

(15)

普通的低階格式，包括一階迎風格式、混合格式、指數格式等擁有良好的穩(wěn)定性和有界性，但精度較低，往往導致過度的數值擴散。而單純的二階格式會使流場的間斷附近產生非物理振蕩，有可能導致迭代的發(fā)散。為了提高精度，研究者們提出了一系列線性對流格式，這些線性對流格式可以用van Leer提出的“k-格式”來概括：

qi+1/2=qi+

(16)

當k取不同的值時能夠獲得不同的高階對流格式，這些格式能在一個維度上使用5點模塊來計算，且在計算線性對流方程時都是線性的。但是“k-格式”的所有成員都不具備有界性，Godunov定理闡明了任何截斷誤差高于一階的線性對流格式都不會是單調的[18]，因此在計算有間斷或梯度很大的情況時會出現嚴重的非物理數值振蕩，從而導致計算的發(fā)散。

為了解決低階格式的低精度和高階格式的不穩(wěn)定的問題，需要構造非線性的對流格式。最常用和有效的一種構造非線性對流格式的方法就是“通量限制器”的使用。

重構過程的插值可以寫成如下“標準半離散通量限制形式”，在均勻網格下：

(17)

式中：ri+1/2為梯度率，表示中心梯度和迎風梯度的比值，本質上是針對變量平滑度的局部測量：

(18)

以ξ方向為例，在非均勻結構網格中水滴流場控制體表面左右兩側的守恒變量可以寫為

(19)

(20)

(21)

式中：ψ(r)即為通量限制器函數，不同的通量限制器函數決定了對流項離散格式的精度和收斂性。

1.3.2 水滴流場發(fā)散原因

對于不使用通量限制器的高階格式計算水滴流場時，導致發(fā)散的原因一般有兩個：

1) 水滴流場中，雖然不存在空氣流場中的激波現象，但是由于物面遮擋效應所導致的無水區(qū)邊緣的α的梯度很大，使用高階精度格式計算時也會出現數值振蕩。

2) 由于歐拉法中水滴速度的單值性引發(fā)水滴流動交匯處出現奇異值,如圖1所示。在物體后部一般會有水滴交匯的現象出現，對于拉格朗日法，每個水滴是獨立運動的，其軌跡可以發(fā)生交叉和重疊，此處α的值為兩股交匯水流的α之和；但對于歐拉法，由于水滴被視為是連續(xù)相，流線不能出現交叉或重疊，兩股水滴流動分別從上至下和從下至上匯聚成一股向同一方向流動，導致了此處的α異常增大，若離散格式缺乏耗散項及時將此處的異常流動抹平，則極其容易導致計算的發(fā)散[8]。

圖1 水流交匯示意圖

因此需要引入限制準則，并在對流項的離散格式中添加通量限制器，以此來抑制水滴流場間斷處的波動和抹平歐拉法中速度單值性帶來的奇異值。

1.3.3 TVD限制準則

總變差減小(Total-Variation Diminishing，TVD)的概念最早是由Harten提出的[19]?？傋儾钍嵌攘亢瘮挡▌訝顩r的一種指標，對于離散的數值解來說，第n個時間步的總變差的定義為

(22)

而TVD的特征為

TV(qn+1)≤TV(qn)

(23)

總變差能夠衡量數值解的振蕩情況，滿足TVD的格式能夠保證數值解的單調和有界。Sweby[20]在顯式離散的背景下推導了用通量限制器函數表達的Sweby’s TVD限制準則：

(24)

Sweby’s TVD限制準則不依賴于CFL數，對于求解半離散對流方程的恒定狀態(tài)解來說，不依賴CFL數的限制準則更加適用。而在空氣流場恒定的情況下，水滴流場就是定常的，其求解過程是使用“偽時間步”推進的?！発-格式”系列所包括的區(qū)域與Sweby’s TVD限制區(qū)域重合的部分，既可以保證總殘差減小又能保證空間二階精度。

國內外有許多的完全或部分滿足上述限制準則的格式被相繼開發(fā)出來，可以統稱為TVD格式，這些格式大致可以分為兩類：逐段線性(Piecewise-Linear)格式和平滑(Gradually-switching Smooth)格式，本文將5種經典的通量限制器應用于水滴流場控制方程的對流項離散，這些通量限制器的ψ-r圖像如圖2所示。

圖2 5種通量限制器的ψ-r圖像

1) Minmod

ψ(r)=max(0,min(r,1))

(25)

Minmod是一種最簡單的限制器，使用時，實際上是只用迎風梯度或中心梯度來修正一階迎風量。當在極值處時r≤0，離散格式直接降為一階迎風格式，而遠離極值處時，離散格式會在中心差分和二階迎風格式中切換。張涵信[16]提出的NND格式就是使用了此限制器，并闡述了中心差分與二階迎風格式的組合使用何以能限制間斷(梯度陡峭)處上下游的數值振蕩。

2) MUSCL

(26)

MUSCL[20-21]限制器是在Sweby’s TVD限制區(qū)域內使用Fromm格式，因此在光滑區(qū)域具有Fromm格式的精度和收斂性。

3) Superbee

ψ(r)=max(0,min(2r,1),min(r,2))

(27)

Superbee[22]限制器在ψ-r圖中是二階Sweby’s TVD限制區(qū)域的上界，此格式是被設計用來詳細刻畫流場間斷的，但是由于在高曲率區(qū)域大量使用一階順風格式導致了光滑區(qū)域的過度壓縮[20]。和Minmod限制器一樣，Superbee限制器在光滑區(qū)域也是使用中心差分與二階迎風格式的組合，不過判斷條件與其相反。

4) Harmonic

(28)

Harmonic[23-24]格式實質上是用中心梯度和迎風梯度的調和平均值來修正一階迎風項，其也是基于Fromm格式的斜率限制器，但只在r=1處具有和Fromm格式相同的斜率。

5) Van Albada

(29)

Van Albada限制器在整個定義域上都是光滑、連續(xù)且可微的，這為離散格式帶來了良好的收斂性。與其他限制器不同的是，此限制器在r≤0時并不為0，即完全不使用一階迎風格式。

將上述通量限制器代入到式(19)、式(20)中即可得到對應的對流項離散格式。

1.4 時間推進格式

對于水滴流場的半離散式(10)，令

(30)

略去下標i,j,k，則式(10)可以寫為

(31)

式中：R稱為方程的殘差。

采用顯示的四步Runge-Kutta法來進行時間推進，具體格式如下：

(32)

式中：α1=1/4;α2=1/3;α3=1/2;α4=1;Δt為時間步長。

1.5 邊界條件和初始條件

水滴流場的求解采用有限體積法，流場的邊界條件通過設置虛網格格心值的方式給出，邊界面上的通量基于虛網格格心值插值得到，水滴流場的邊界條件分為自由來流邊界和物面邊界。

自由來流邊界上α如1.1節(jié)中假設1)所述，速度向量根據攻角設置；物面邊界的設置需要判斷撞擊區(qū)和無水區(qū)，u為水滴速度向量，n為物面的外法向量。當u·n≤0時，即為撞擊區(qū)，此時水滴撞上物面，并被全部捕捉后結冰，計算時將水滴視為穿過此邊界，該區(qū)域虛網格的α和u均與貼近物面的第1層網格的值保持相同；當u·n>0時,即為無水區(qū)，此時相當于水滴從物面上流出，實際上此區(qū)域是被物體遮擋，沒有水滴的流入，理論上此邊界的α應設置為0，為了避免計算過程中出現數值錯誤，該區(qū)域虛網格的α取物面上第1層網格的α的10-9倍，u與其保持一致。

針對物體之后在迭代過程中會存在奇異值的問題，將初始流場的α設置成一個接近于0的小值可以有效的解決此問題，相當于水滴從自由來流邊界漸漸傳輸到整個計算域，物體之后的無水區(qū)一直是保持沒有水滴的狀態(tài)，從根本上避免了水滴流動交匯產生奇異值的問題。

但上述初始條件的計算較慢，本文的初始條件設置如下：水滴流場全場的初始α用工況設置中的LWC計算，初始速度沿攻角方向，速度大小由來流馬赫數給出。TVD格式在這種初始條件下物體之后雖然存在密度脈沖，但不會引起計算的發(fā)散。

1.6 水滴收集率計算

水滴流場的求解主要是為了獲得物面水滴收集率，為結冰或防冰計算提供輸入，水滴收集率的定義是物面表面微元的實際水收集量與最大可能的水收集量間的比值，用于衡量物面收集水滴能力的大小，水滴收集率的計算采用如下公式：

(33)

式中：αn為物面附近的水滴體積分數；α∞為無窮遠處的水滴體積分數；u∞為無窮遠處水滴的速度矢量。

2 結果與分析

2.1 NACA0012翼型

對于NACA0012翼型，計算了如表1所示的4種工況和如表2所示的3種粗細的網格。

圖3為Case1工況下一階迎風格式和5種帶不同限制器的TVD格式計算得到的水滴收集率β與實驗數據的對比。圖中：s為到駐點的距離；c為機翼的弦長。圖4為Case1工況下計算過程中α的收斂曲線，可以看出，相比于一階迎風格式，具有限制器的TVD格式均能更好的計算水滴收集率。一階格式的誤差主要出現在駐點附近，駐點附近會出現計算的水滴收集率偏大的現象，結合圖5中水滴體積分數的分布，可看出，一階格式的計算結果會在駐點附近產生水含量的集中。

表1 NACA0012翼型計算工況

表2 NACA0012翼型計算所用網格參數

TVD格式的計算結果與實驗數據[25]相比，在駐點附近偏小，撞擊極限偏小，且有輕微的不對稱現象。收集率和撞擊極限偏小的原因可能在于實驗過程無法保證所有水滴的直徑均相同，而是在一個區(qū)間內分布；而不對稱的來源可能是網格的不對稱和重力的作用效應，但小水滴的重力和空氣對水滴的作用力相比是一個小量，小水滴的重力幾乎可以忽略不計，因為有限體積方法對流項的離散需要用到每個單元網格的面法向量，所以網格的不對稱可能會造成機翼上下區(qū)域水含量

圖3 Case1的水滴收集率

圖4 Case1中α的迭代收斂曲線

圖5 Case1的水滴體積分數云圖

輕微出入，因此導致水滴收集率的不對稱。

圖6為Case2工況下網格適應性分析的結果。對于文中所使用的5種TVD格式，中網格和細網格的計算結果基本一致，只在駐點處輕微有一點出入。實際上粗網格的計算結果趨勢也與中網格和細網格的一致，只是因為駐點處的點不夠密，無法精細的描述出物面駐點處的水滴收集率分布?？傮w來說TVD格式網格適應性較好，在普通密度的網格中就能計算出良好的水滴收集率分布結果。

圖7、圖8分別為Case3、Case4下，使用粗細不同的兩套網格計算得到的水滴收集率，網格的參數如表2所示，計算結果與Lewice和DROP3D軟件計算的[26]進行對比。不論是粗網格還是細網格，一階迎風格式均在駐點處有水含量過度集中導致的水滴收集率偏大。當水滴直徑較小時(20 μm)，Minmod限制器也出現了駐點處水滴收集率偏大的情況。對比粗網格的計算結果，使用細網格時，一階格式與TVD格式的計算結果更加接近，幾乎只在駐點附近不同。說明當網格比較稀疏時，TVD格式能夠較好的提升計算水滴收集系數的精度。使用細網格計算時，Case4狀態(tài)下，Superbee格式計算水滴流場的迭代過程發(fā)生了發(fā)散，主要原因在是翼型遮擋部位的無水區(qū)邊緣產生了無法抹平的數值波動。

圖6 Case2中用TVD格式在不同尺寸網格下計算的水滴收集系數

圖7 Case3的水滴收集系數

圖8 Case4的水滴收集系數

2.2 NACA23012翼型

NACA23012翼型的計算工況如表3所示。

表3 NACA23012翼型計算工況

從圖9水滴收集率的計算結果來看，一階格式的計算結果在駐點處還是偏大，所采用的5種TVD格式的計算結果相差無幾。TVD格式間的差別主要體現在對間斷的分辨率的差別上，由于水滴收集率只涉及物面前緣的流場，這一部分并沒有流動的間斷，因此所采用的5種TVD格式的計算結果非常接近。

對于Case5(圖9和圖10)，除了撞擊極限外，一階格式的計算結果偏大，而TVD格式的計算結果與實驗數據基本一致。撞擊極限偏小的原因可能是由于本文程序是按照單一的水滴尺寸計算的，而實驗和Lewice中的水滴并不是單一尺寸，而是呈現一定的分布的，其中包含有一定量大尺度粒徑的水滴，這部分水滴撞擊會導致撞擊極限更大；對于Case8(圖11和圖12)，一階格式在駐點處收集率偏大，其余部分與TVD格式基本一致。而與實驗結果相比，不論是一階格式還是TVD格式，均在駐點處偏大，且撞擊極限偏小，但總體來說，本文程序的計算結果較Lewice[26]更接近實驗數據[27]。

圖9 Case5的水滴收集系數

圖10 Case5中α的迭代收斂曲線

圖11 Case8的水滴收集系數

圖12 Case8中α的迭代收斂曲線

從α的收斂曲線上看，殘差變化平穩(wěn)后，不同格式的殘差從大到小的順序為Superbee、MUSCL、Van Albada、Minmod、Harmonic、一階迎風格式。對于Case8，使用Superbee格式計算時發(fā)生了發(fā)散，發(fā)散的原因是在無水區(qū)邊緣產生了數值波動且無法及時抹平。

圖13為使用Minmod限制器時在Case7下計算得到水滴體積分數云圖，可以看到尾部有明顯的高密度脈沖。圖14是在其中取y/c=0.036截面得到的α隨橫坐標x/c的變化趨勢，可以看出，單值性導致的物體后方的高密度區(qū)域會隨著deq的增加而愈發(fā)趨于明顯和陡峭。

圖15是使用不同的對流項離散格式后取截面后得到的α隨橫坐標x/c的變化，使用TVD格

圖13 Case7下使用Minmod格式的水滴體積分數云圖

圖14 Minmod格式下不同粒徑尾部大密度區(qū)域對比

圖15 Case7下不同限制器尾部大密度區(qū)域對比

式時在物體后方產生的高密度區(qū)域要比使用一階格式時更加的明顯，陡峭程度從大至小依次為：Superbee、MUSCL、Van Albada、Minmod、Harmonic、一階迎風格式，這與殘差從大至小的順序一致。由此可知，在殘差變化平穩(wěn)后，影響水滴流場迭代求解的因素就是物體之后的密度脈沖。

圖16為使用另一種初始條件，即將初始α設置成一個接近于0的小值時，在Case7下使用Minmod限制器計算得到的水滴體積分數云圖?？梢钥吹酱朔N初始條件的設置可以有效地解決物體之后密度脈沖的問題。

2.3 30P30N三段翼

30P30N三段翼的計算工況如表4所示。圖17為30P30N三段翼的外形示意圖，機翼由前緣縫翼、翼身和后緣襟翼3部分組成。機翼弦長c的定義為當前緣縫翼和后緣襟翼均未偏轉時整體的長度，在此算例中，前緣縫翼和后緣襟翼的偏轉角度均為30°。前緣縫翼的上下兩端容易發(fā)生邊界層分離和轉捩。

圖16 在初始水滴體積分數接近于0時Case7下使用Minmod格式的水滴體積分數云圖

Table 4 Calculation conditions for 30P30N three-element airfoil

圖17 三段翼30P30N的外形示意圖

圖18～圖20分別為前緣縫翼、翼身和后緣襟翼上的水滴收集系數分布曲線。使用Van Albada格式計算時，在前緣縫翼的上端出現了不可抹平的奇異值。對于其他幾種計算格式，在前緣縫翼上，TVD格式在駐點附近和上端撞擊極限處跟一階迎風格式有所出入；在翼身上，5種格式的計算結果基本一致；在后緣襟翼上，4種TVD格式的計算結果幾乎沒有出入，而一階迎風格式相比下水滴收集率的分布整體向左側偏移。

圖21為Case9中的收斂曲線，平均殘差從大到小依次為Superbee、MUSCL、Minmod、Harmonic、一階迎風格式，順序與之前計算的幾種工況一致。TVD格式中的Van Albada格式無法處理此算例中邊界層分離的情況，而其他的格式對復雜構型的適應性較好。

圖18 Case9中前緣縫翼的水滴收集系數

圖19 Case9中翼身的水滴收集系數

圖20 Case9中后緣襟翼的水滴收集系數

圖21 Case9中α的迭代收斂曲線

3 結 論

將TVD格式應用到水滴流場求解的對流項離散中，通過計算分析，有如下結論：

1) 基于水滴收集率計算和實驗數據的對比，TVD格式較一階迎風格式的精度有較好的提升，特別是在駐點附近，但所采用的5種TVD格式計算的精度無明顯差別。

2) TVD格式有較好的網格適應性，在較稀網格條件下，兩種格式的計算精度差距更大。

3) 在歐拉法水滴流場迭代計算，殘差的主要來源之一是物體之后的密度脈沖，隨著水滴粒徑的增大，高密度區(qū)域范圍和密度都隨之增加。將初始條件設為接近于0的小值時可以避免密度脈沖。若使用本文中的的初始條件設置方式，不同的TVD格式對此區(qū)域的抹平能力不同，從強至弱依次為：Harmonic、Minmod、Van Albada、MUSCL、Superbee。

4) Superbee格式的穩(wěn)定性較差，不僅體現在物體尾部高密度區(qū)域抹平能力較差，而且在部分工況下，無法處理水滴流場間斷導致的數值振蕩。

5) 針對復雜構型，Van Albada格式在計算過程中發(fā)生了發(fā)散，而其他格式都其表現出良好的適應性。

6) 在不修改水滴流場控制方程的前提下，探索用通量限制器構造的TVD格式求解水滴流場的可能性，其計算結果也證明此類格式是有效的。