專利名稱：可壓縮旋流場的數(shù)值模擬方法
技術(shù)領(lǐng)域：
本發(fā)明涉及計算流體力學(xué)(CFD Computational Fluid Dynamics)領(lǐng)域中的一種數(shù)值模擬方法，具體是一種模擬 可壓縮旋流場的數(shù)值方法。
背景技術(shù)：
計算流體力學(xué)綜合了流體力學(xué)、應(yīng)用數(shù)學(xué)、計算機科學(xué)，是一門應(yīng)用性極強的學(xué)科。流體力學(xué)問題的數(shù)值模擬以其低成本、直觀性強的優(yōu)勢，在流體流動的機理探索、工業(yè)產(chǎn)品設(shè)計等各個相關(guān)領(lǐng)域占據(jù)重要地位。計算流體力學(xué)面臨的最大的問題和挑戰(zhàn)之一即是如何提高數(shù)值模擬的精度，忠實地表現(xiàn)流體流動的特性。影響流體力學(xué)問題的數(shù)值模擬的精度的重要因素之一是當使用數(shù)值方法求解流體控制方程，即歐拉(Euler)方程或者納維爾-斯托克斯(Navier-Stokes)方程時,會產(chǎn)生數(shù)值耗散(numerical diffusion),造成數(shù)值解的誤差。例如數(shù)值方法中對控制方程中的對流項的各種空間離散方法(如中心差分、迎風(fēng)差分)、時間離散方法(如顯示時間積分、隱式時間積分)、湍流模型(如雙方程模型、大渦模擬)的使用，以及計算網(wǎng)格的正交性都會產(chǎn)生不同程度的數(shù)值耗散。此外，數(shù)值模擬中還經(jīng)常要使用一種人工數(shù)值耗散(artificialdiffusion)技術(shù),其目的是通過適當降低計算精度而獲得穩(wěn)定的數(shù)值解。例如流場中的不連續(xù)(discontinuity)界面的捕捉,如激波(shock)的捕捉，即是依靠加入適量的人工耗散項，以避免二階以上精度的數(shù)值解在流動變量梯度較大的地方出現(xiàn)數(shù)值振蕩現(xiàn)象。數(shù)值耗散可以理解為是流場中的一種能量損失，這種能量損失在某種程度上使得數(shù)值模擬結(jié)果不能忠實的體現(xiàn)流體的流動特性，降低了計算精度。先進的數(shù)值方法應(yīng)該是在保證獲得穩(wěn)定性的數(shù)值解的前提下，將數(shù)值耗散減至最小。數(shù)值耗散對流場的數(shù)值模擬結(jié)果最明顯的影響體現(xiàn)在對流動變量的間斷界面的捕捉。如前所述，激波是強間斷界面，跨過激波，速度、密度，壓力均不連續(xù)，對其捕捉必須加入一定的人工數(shù)值耗散，以抵消數(shù)值解在不連續(xù)處的振蕩。因為激波前后存在熵增，即能量的損失，所以，通過加入人工數(shù)值耗散捕捉激波具有合理的物理意義。但是，過大的數(shù)值耗散會使數(shù)值解獲得的激波界面變得模糊，降低了數(shù)值解對激波強度和空間位置的預(yù)測精度。流場中存在另一類流動不連續(xù)現(xiàn)象，即接觸不連續(xù)(contact discontinuity)。相對激波而言，接觸不連續(xù)的特點是弱不連續(xù)，不連續(xù)界面是滑移線(slip-line)，跨過滑移線，沒有質(zhì)量傳遞、密度和切向速度不連續(xù)，但徑向速度和壓力連續(xù)。一類旋流場(Vortex-dominated Flows)中具有最典型的接觸不連續(xù)界面的流動現(xiàn)象。例如，具有三角機翼的飛行器，具有較大展翼比，在大攻角飛行時周圍的流場是一明顯的旋渦運動為主。其他例如直升機旋翼后方的流場、渦輪發(fā)動機轉(zhuǎn)子周圍的流場，等等，同樣會遇到類似的現(xiàn)象。計算機的數(shù)值模擬中對于這種接觸不連續(xù)的捕捉更加困難，因為數(shù)值方法中的數(shù)值耗散即使很小也會使弱不連續(xù)界面變得模糊，降低了數(shù)值解對流場的預(yù)測精度，這也是旋流場的數(shù)值模擬技術(shù)成為CFD領(lǐng)域的重大挑戰(zhàn)的原因。描述可壓縮粘性流體運動的控制方程包括連續(xù)方程和動量方程，即Navier-Stokes方程,分別由以下兩式給出，^+ν·(ρν)=0,Cl)
dt^+(ν·ν)ν = --Vp + i^v(v · v)+ ^V2V + δ,(2)
式中，V是速度矢量，V = i + i^ + wfc ,包含直角坐標系下的三維i, j, k分量u,V,W ; P、p、t分別是密度、壓力和時間。式中▽算子表示為7 = + + ;符號 代
表內(nèi)積計算。動量方程(2)中的代表壓縮性引起的剪切動力粘性項、代
PP
表正壓運動粘性項(，是拉普拉斯算子V2)、·《代表體積力項。對于
V2dx dy2 dz7· B
一類以旋渦運動為主的流場，常引入描述旋流流動的變量-渦量(vorticity)對流場進行
描述。渦量用凾表示，在直角坐標系下有+ + 按照渦量的定義有
ij k
Fd d d (dw av\ f du(dv 5 Vco =s Vx V = — — — = --— I +---7 + -----— \k 9\o)
dx dy dzdz) ydz dx J (St dy J
uvw式中符號X代表叉乘運算。對方程(2)等號兩邊作叉乘運算，并利用方程(I)和渦量的定義，可以推導(dǎo)出描述
可壓縮粘性旋渦流場的渦量方程，
營+( ·ν)δ=φ·ν)^-*(ν·ν)+ι^2ω + ^^+νχ(|^7(ν·々)) + νχδ。 (4)上式使用了可壓縮流的動力粘度V，定義為V= μ/P。公式⑷中等號左邊表示渦量的變化率；右邊第一項表示渦量的三維變形產(chǎn)生的渦量；第二項表示可壓縮性對渦量的影響；第三項和第五項表示粘性產(chǎn)生的渦量耗散；第四項表示壓力梯度和密度梯度不平行時產(chǎn)生的渦量；第六項表示體積產(chǎn)生渦量的作用。公式(4)對于二維情況下，( ·ν)^ = 0:
對于正壓流場，^^=0:對于無粘流，且▽xgwiv.vjl =。，該式可以被簡化。對于渦量場的影響最顯著的是渦量的耗散項VV2O,主要是因為流體粘性作用引起的，同樣的道理，數(shù)值方法中需要加入的人工耗散是以粘性形式起作用，直接導(dǎo)致了流體粘性引起的渦量耗散以外的渦量耗散。流體的粘性作用，包括人工粘性的作用具有橢圓型特征，即粘性作用在空間是各向同性的。為了提高旋流場的數(shù)值模擬的精度，一類方法是加密計算網(wǎng)格，在更加細小的空間尺度內(nèi)求解流體控制方程。加密計算網(wǎng)格首先會使計算量加大，增加計算成本。此外，數(shù)值計算的誤差隨著計算網(wǎng)格的增加會不斷積累，在一定程度上造成相反的效果。另一類方法是在流場中采用物理模型來增加流場中渦量的強度。例如在流場中加入點渦模型，可以人為地增加渦量，或者在流場局部直接求解渦量方程(公式4)，以減小渦量的輸運過程中的數(shù)值耗散。但是，這些方法在應(yīng)用上仍受到一定限制。點渦模型是在預(yù)先明確旋渦發(fā)生位置的前提下才能使用，僅適合一些簡單的流動現(xiàn)象。除了二維不可壓縮正壓流場，渦量方程的求解比與欲求解的Euler、Navier-Stokes方程更為復(fù)雜。最后一類方法是開發(fā)具有最小數(shù)值耗散的空間離散算法，例如，對流項離散方法ENO、WENO方法，主要是針對求解直角坐標系下控制方程(I)和(2)的求解過程中使用的。但是，這些數(shù)值方法的數(shù)值耗散依然過大，會降低弱不連續(xù)的捕捉精度，使不連續(xù)界面變得模糊?？傊瑸榱烁泳_地模擬可壓縮流的旋渦運動，需要進一步改進對于流場中諸如旋渦運動一類的接觸不連續(xù)的捕捉機制。一種途徑是按照公式(4)給出的渦量變化的原理，在流場中補償人工耗散項所產(chǎn)生的渦量耗散，同時在數(shù)值解求解構(gòu)成中將補償?shù)臏u量耗散進行高精度空間離散，更精確地模擬流場中的旋渦運動，形成一種新的可壓縮流的旋渦運動的數(shù)值模擬技術(shù)。

發(fā)明內(nèi)容
本發(fā)明涉及計算流體力學(xué)領(lǐng)域中一種模擬可壓縮旋流場的數(shù)值方法，具體是一種在流場中補償人工耗散項產(chǎn)生的渦量耗散，同時在數(shù)值解求解構(gòu)成中將補償?shù)臏u量耗散進行高精度空間離散，更精確地模擬可壓縮流的旋渦運動的數(shù)值模擬方法-渦量耗散的補償技術(shù)。本發(fā)明采用的技術(shù)方案可壓縮流體流動的控制方程是一組偏微分方程，包括連續(xù)、動量和能量方程。無粘流的動量方程又被稱為Euler方程；粘性流動量方程又被稱為Navier-Stokes方程。工程應(yīng)用上，流體流動的控制方程的數(shù)值解的求解常采用有限體積法(Finite Volume Method)將控制方程在預(yù)先生成計算網(wǎng)格上進行離散化。有限體積法需要將控制方程寫成守恒積分形式，每一個計算網(wǎng)格單元被是為控制體，控制體內(nèi)流場守恒變量隨時間的變化是通過在每一個特定時間間隔內(nèi)，通過控制體表面(即計算網(wǎng)格邊界)上的守恒變量的通量值的和來確定。這也是的流體控制方程數(shù)值解的基本原理。上述通過計算網(wǎng)格界面的通量值包括對流通量項、人工耗散項和粘性項。各種不同的空間離散方法，如中心差分的代表JST方法、考慮偏微分方程組的特征的代表VanLeer、⑶SP、AUSM、Roe等方法，都是用來逼近在計算網(wǎng)格邊界上的對流項和人工耗散項通量值。粘性項的空間離散，因為粘性項的橢圓形特性，均采用中心差分。如前所述，空間離散的數(shù)值方法中需要加入一定的人工數(shù)值耗散。其加入的形式可以是直接的，如JST、CTSP、Roe或者是隱含式的，如VanLeer、AUSM方法。人工耗散在使偏微分方程組獲得數(shù)值解穩(wěn)定的作用也降低了旋渦結(jié)構(gòu)的捕捉精度，構(gòu)成了數(shù)值模擬誤差的一個來源。對于旋流場，加入的人工耗散使得流場中的渦量過度耗散。如前所述，對于渦量場的影響最顯著的是渦量的耗散項，所以要在旋流場中減小因加入人工耗散而造成的渦量耗散。一個使得旋流場的渦量耗散保持最小化的原則是;渦量的耗散方向和渦量的模的最大變化梯度方向保持一致。而人工耗散的加入使得這兩個方向分離。因此，本發(fā)明采取的解決方法是，在計算網(wǎng)格邊界上加入一個渦量耗散補償項通量，補償因為渦量的耗散方向和渦量的模的最大變化梯度方向的分離而造成的渦量過度耗散。這個渦量耗散補償項應(yīng)該與人工耗散項等價、作用相反。該方法被稱為渦量耗散的補償技術(shù)。其原理圖見圖I。圖中表示，通過計算網(wǎng)格的界面，通常是由兩個通量，即對流項通量(包含粘性項通量)和加入的人工耗散項，通過加入與人工耗散項相反的作用、等價的渦量耗散補償項通量，可以抵消部分數(shù)值耗散。渦量耗散的補償技術(shù)有以下三個特點I.這個加入渦量耗散補償項是在對流項加入人工耗散項以后加入的，所以不影響原來數(shù)值解的穩(wěn)定性策略；2.這個被加入的渦量耗散補償項不是由有流體粘性引起的，而由數(shù)值粘性引起的；3.渦量的耗散項與渦量的的模的最大梯度變化的方向一致，則不需要渦量耗散補 償項。按照量綱分析，渦量耗散補償項I；是密度和單位渦量耗散補償項ξ的乘積， ·定義為fo:式中，S 代表渦量的模的梯度最大變化的方向；νω代表數(shù)值粘性系數(shù)，被稱為渦運動的補償粘性系數(shù)，與流體的運動粘度有相同的量綱，是一個標量；R。代表補償渦量的特征半徑。式中符號X代表叉乘運算。當?shù)販u量的模的梯度最大變化的方向S 表示為(6)其中，Φ是渦量的模；是渦量的模的梯度；是渦量的模的梯度的模，分別定義為盧=間=+ω; + ω] >(7)= Ujjlk，(8)
ax dy dz
_ NI =摘 + ⑵ +{~]。 (9)渦量的模的梯度的模表示了最大渦量的變化程度。數(shù)值粘性系數(shù)的確定因為人工耗散項是以流體的粘性耗散的形式體現(xiàn)的，對于動量方程具有速度的拉普拉斯算子(V2V)的量綱，對具有速度的拉普拉斯算子的量綱人工耗散項取旋度
則有
(作 Vx 運算)，VX(V2^)=V2W(V(V2))Xt(10)如果舍去式中右邊的三階導(dǎo)數(shù)項，速度的拉普拉斯算子的旋度完全變成渦量的拉普拉斯算子，既是渦量方程(4)中的耗散項，說明了公式(2)中的人工耗散項與公式(3)中的渦量耗散項的等價性。所以渦量補償也就是以渦量的形式補償了數(shù)值方法中的人工耗散，消除了這部分誤差，提高數(shù)值模擬的精度。各種空間離散的數(shù)值方法中的人工耗散項均需要用一個具有速度量綱的參數(shù)作為人工耗散的尺度放大系數(shù)。例如，在JST方法中使用譜半徑(spectral radii)作為人工耗散的尺度放大系數(shù)。譜半徑是指偏微分方程在計算網(wǎng)格邊界處的最大特征值。對流項的譜半徑Λ是A = I V|+C，(11)式中C是當?shù)氐穆曀?；V是抗變速度，定義為=
B = In1, Hyt zj是計算網(wǎng)格界面的法線方向向量。如果對具有速度量綱的人工耗散的放大尺度系數(shù)取旋度，可得到具有渦量的量綱的參數(shù)(見公式3)。這個放大系數(shù)可以用渦量的形式表示，再用來放大渦量耗散項V2S,而且這個系數(shù)具有流體運動粘度的量綱，因此，該系數(shù)可以被認為是一個數(shù)值粘性引起的，被稱為渦運動補償粘性系數(shù)。因為渦運動補償是在空間不同方向進行的，渦運動補償粘性系數(shù)向量定義為

Ve*(12) 其中，補償渦量的半徑向量Ret -Rii+RjaJ+ Mk，定義成Rtf=Ato(14)上式和公式(5)中，補償渦量的特征半徑R。在二維和三維空間分別定義為Rc = ,(13)/Jc=IqX0(14)
2補償渦量的特征半徑的幾何意義解釋在計算網(wǎng)格上的補償渦量一定有一個與之等效的環(huán)量，而環(huán)量是圍繞該網(wǎng)格內(nèi)某一封閉區(qū)域上的速度的線積分。該封閉區(qū)域的大小用公式(13)或(14)表示，計算網(wǎng)格單元的空間尺度的二分之一作路徑形成的區(qū)域的半徑，是渦量作用在該計算網(wǎng)格內(nèi)的平均空間范圍。公式(14)中，渦量在最大作用的方向系數(shù)向量，+ + 。其作用是為了在流體的粘性層內(nèi)減小渦量補償，以及考慮到計算網(wǎng)格的大變形率，需要根據(jù)渦量向量的各向異性程度進行補償渦量的半徑向量的坐標方向的各向異性化，則k1 =l + max —, — (15)
ωχkJ =l + max —, —，(16)
COy G>yAri= I+ max —o(17)
A ωζ最終，在計算網(wǎng)格界面上的抗變渦運動補償粘性系數(shù)νω定義為Va SB ν ·η =nX + [Ri Jny+ (r^ J nz] 。(18)按照以上推導(dǎo)，圖2給出計算渦量耗散補償項的步驟流程。公式(5)給出的渦量耗散補償項可以在三維直角坐標系下展開，其形式表示為
權(quán)利要求
1.一種模擬可壓縮旋流場的數(shù)值方法，其特征在于在計算網(wǎng)格的界面加入渦量耗散補償項，其等于流體密度P乘以單位質(zhì)量渦量耗散補償項向量t，t的表達形式是其中，代表渦量的模的梯度最大變化的方向、νω代表渦運動補償粘性系數(shù)，與流體的運動粘度有相同的量綱，是ー個標量、R。代表補償渦量的特征半徑；渦量用A表示，在直角坐標系下有 符號X代表叉乘運算；v2是拉普拉斯算子。
2.根據(jù)權(quán)利要求I所述的ー種可壓縮旋流場的數(shù)值模擬方法，其特征在于所述的渦量的模的梯度最大變化的方向··表示為 唭中，Φ是渦量的摸、是渦量的模的梯度、|V#|是渦量的模的梯度的摸，分別定義為
3.根據(jù)權(quán)利要求I所述的ー種可壓縮旋流場的數(shù)值模擬方法，其特征在于所述的補償渦量的特征半徑R。在ニ維和三維空間分別定義為^和& =^Qk，其中Ω在ニ維和三維空間分別代表計算網(wǎng)格面積和體積。
4.根據(jù)權(quán)利要求I所述的ー種可壓縮旋流場的數(shù)值模擬方法，其特征在于所述的抗變渦運動補償粘性系數(shù)νω定義為渦運動補償粘性系數(shù)向量與計算網(wǎng)格界面的法線方向向Mn = [η,, rty, wzj的點乘，即νβ ■辛β · ,其中，符號·代表點乘計算。
5.根據(jù)權(quán)利要求4所述的渦運動補償粘性系數(shù)向量Ve，其特征在干，定義為補償渦 量的半徑向量與渦量的模值比，即
6.根據(jù)權(quán)利要求5所述的補償渦量的半徑向量I，其特征在干，^定義成補償渦量的特征半徑R。與渦量在最大作用的方向系數(shù)向量的乘積，即^ =At，其中，
7.一種模擬可壓縮旋流場的數(shù)值方法，其特征在于在計算網(wǎng)格的界面加入渦量耗散補償項，其在三維直角坐標系下展開形式表示為
全文摘要
本發(fā)明是一種模擬可壓縮旋流場的數(shù)值方法。其主要特征是在計算網(wǎng)格邊界上加入一個渦量耗散補償項通量，補償因為渦量的耗散方向和渦量的模的最大變化梯度方向的分離而造成的渦量過度耗散。這個渦量耗散補償項應(yīng)該與人工耗散項等價、作用相反。該方法被稱為渦量耗散的補償技術(shù)。這個被加入渦量耗散補償項是在對流項加入人工耗散項以后加入的，不影響原來數(shù)值解的穩(wěn)定性策略、是由數(shù)值粘性引起的、方向是沿著渦量的耗散項與渦量的模的最大梯度變化方向的叉乘方向。
文檔編號G06F17/50GK102682146SQ20111038871
公開日2012年9月19日 申請日期2011年11月30日 優(yōu)先權(quán)日2011年11月30日
發(fā)明者路明 申請人:天津空中代碼工程應(yīng)用軟件開發(fā)有限公司