本發(fā)明涉及計(jì)算流體力學(xué)應(yīng)用領(lǐng)域中的一種數(shù)值方法，具體是一種模擬風(fēng)力葉片尾跡的數(shù)值方法。

背景技術(shù)：

大型風(fēng)力發(fā)電機(jī)目前已經(jīng)有著廣泛的應(yīng)用。對(duì)于風(fēng)力葉片周圍的流場(chǎng)運(yùn)動(dòng)規(guī)律性的掌握，對(duì)于風(fēng)力葉片的設(shè)計(jì)起著重要作用。一個(gè)重要的流動(dòng)特征即是葉片周圍的旋渦流動(dòng)。如圖1所示，風(fēng)力葉片的尾跡流場(chǎng)包括兩類以旋渦運(yùn)動(dòng)為主(vortex-dominated)的流動(dòng)。一類是葉片的翼展方向的尖端處，即翼尖處產(chǎn)生的翼尖渦；另一類是葉片表面產(chǎn)生的脫落渦。翼尖渦的形成是由于旋翼上下表面的壓力不同，會(huì)引起旋翼下方的氣流繞經(jīng)翼尖向機(jī)翼上方流動(dòng)，因而形成了旋渦，該旋渦從翼尖脫落向下游延伸，在旋翼尾跡形成了螺旋形狀的流動(dòng)模式，是無(wú)粘流動(dòng)。旋翼表面的脫落渦的產(chǎn)生是由于旋翼表面的粘性邊界層中的逆壓梯度引起粘性邊界層的分離所致。翼尖渦和脫落渦形成了風(fēng)力葉片尾跡的旋渦體系。

對(duì)于風(fēng)力葉片的尾跡的流動(dòng)有多種研究手段。其中，計(jì)算機(jī)的數(shù)值模擬技術(shù)有著重要地位，是計(jì)算流體力學(xué)在該領(lǐng)域的一個(gè)擴(kuò)展應(yīng)用。計(jì)算流體力學(xué)綜合了流體力學(xué)、應(yīng)用數(shù)學(xué)、計(jì)算機(jī)科學(xué)，是一門應(yīng)用性極強(qiáng)的學(xué)科。流體力學(xué)問(wèn)題的數(shù)值模擬以其低成本、直觀性強(qiáng)的優(yōu)勢(shì)，在流體流動(dòng)的機(jī)理探索、工業(yè)產(chǎn)品設(shè)計(jì)等各個(gè)相關(guān)領(lǐng)域占據(jù)重要地位。直升機(jī)旋翼尾跡的數(shù)值模擬面臨的最大的問(wèn)題和挑戰(zhàn)之一即是如何提高數(shù)值模擬的精度，降低誤差，忠實(shí)地表現(xiàn)翼尖渦流動(dòng)的特性。

影響流體力學(xué)問(wèn)題的數(shù)值模擬的精度的重要因素之一是：當(dāng)使用數(shù)值方法求解流體控制方程，即歐拉(euler)方程或者納維爾-斯托克斯(navier-stokes)方程時(shí)，會(huì)產(chǎn)生數(shù)值耗散(numericaldiffusion)，造成數(shù)值解的誤差。例如數(shù)值方法中對(duì)控制方程的對(duì)流項(xiàng)的空間離散方法(如中心差分、迎風(fēng)差分)、時(shí)間離散方法(如顯示時(shí)間積分、隱式時(shí)間積分)、湍流模型(如雙方程模型、大渦模擬)的使用，以及計(jì)算網(wǎng)格正交性都會(huì)產(chǎn)生不同程度的數(shù)值耗散。此外，數(shù)值模擬中還經(jīng)常要使用一種人工數(shù)值耗散(artificialdiffusion)技術(shù)，其目的是通過(guò)適當(dāng)降低計(jì)算精度而獲得穩(wěn)定的數(shù)值解。

數(shù)值耗散對(duì)流場(chǎng)的數(shù)值模擬結(jié)果的最明顯的影響體現(xiàn)在對(duì)流動(dòng)變量的不連續(xù)界面(discontinuity)的捕捉。流場(chǎng)中一種強(qiáng)不連續(xù)的界面的捕捉，例如激波(shock)的捕捉， 即是依靠加入適量的人工耗散項(xiàng)，以避免數(shù)值解在流動(dòng)變量梯度變化較大的地方出現(xiàn)數(shù)值振蕩現(xiàn)象。數(shù)值耗散可以理解為是流場(chǎng)中的一種能量損失，這種能量損失在某種程度上使得數(shù)值模擬結(jié)果不能忠實(shí)的體現(xiàn)流體的流動(dòng)特性，降低了計(jì)算精度。先進(jìn)的數(shù)值方法應(yīng)該是在保證獲得穩(wěn)定性的數(shù)值解的前提下，將數(shù)值耗散減至最小。激波是強(qiáng)間斷界面，對(duì)其捕捉必須加入一定的人工數(shù)值耗散。因?yàn)榧げㄇ昂蟠嬖陟卦?，即能量的損失，所以，通過(guò)加入適當(dāng)?shù)娜斯?shù)值耗散捕捉激波具有合理的物理意義。但是，諸如直升機(jī)旋翼尾跡一類的流動(dòng)，存在著流場(chǎng)中另一類流動(dòng)不連續(xù)現(xiàn)象，即接觸不連續(xù)(contactdiscontinuity)。旋渦的產(chǎn)生自然和其周圍的流體產(chǎn)生一個(gè)不連續(xù)面。這個(gè)接觸不連續(xù)面相對(duì)激波而言是弱不連續(xù)，跨過(guò)不連續(xù)界面，壓力和法向速度是連續(xù)的。數(shù)值模擬中對(duì)于這種接觸不連續(xù)的捕捉更加困難，因?yàn)閿?shù)值方法中的數(shù)值耗散即使很小也會(huì)使弱不連續(xù)界面變得模糊，降低數(shù)值解對(duì)流場(chǎng)的預(yù)測(cè)精度，這也是直升機(jī)旋翼尾跡一類的旋流場(chǎng)的數(shù)值模擬技術(shù)成為cfd領(lǐng)域的重大挑戰(zhàn)的原因。

為了提高風(fēng)力葉片尾跡流場(chǎng)的數(shù)值模擬的精度，一種方法是加密計(jì)算網(wǎng)格，在更加細(xì)小的空間尺度內(nèi)求解流體控制方程。加密計(jì)算網(wǎng)格首先會(huì)使計(jì)算量加大，增加計(jì)算成本。此外，數(shù)值計(jì)算的誤差隨著計(jì)算網(wǎng)格的增加會(huì)不斷積累，在一定程度上造成相反的效果。另一種方法是在流場(chǎng)中采用物理模型來(lái)增加流場(chǎng)中描述旋流流動(dòng)的變量—渦量(vorticity)的強(qiáng)度。例如在流場(chǎng)中加入點(diǎn)渦模型，可以人為地增加渦量；或者在流場(chǎng)局部直接求解渦量方程，以減小渦量的輸運(yùn)過(guò)程中的耗散。但是，這些方法在應(yīng)用上仍受到一定限制。點(diǎn)渦模型是在預(yù)先明確旋渦發(fā)生位置的前提下才能使用，僅適合一些簡(jiǎn)單的流動(dòng)現(xiàn)象。除了二維不可壓縮正壓流場(chǎng)，渦量方程比與欲求解的euler、navier-stokes方程更為復(fù)雜。

二十世紀(jì)初期，美國(guó)科學(xué)家johnsteinhoff提出了一種提高不可壓縮(incompressible)旋流場(chǎng)的求解精度的數(shù)值方法，渦量限制法(vorticityconfinement)。該方法的原理是依靠在流場(chǎng)中加入限定的渦量以抵消數(shù)值耗散來(lái)模擬旋流場(chǎng)的流動(dòng)狀態(tài)。具體表現(xiàn)形式是在流體控制方程的動(dòng)量方程中的源項(xiàng)位置，加入一個(gè)渦量形式的體積力項(xiàng)，從數(shù)值耗散中將渦量減去，克服數(shù)值耗散造成的旋渦場(chǎng)的接觸不連續(xù)界面的模糊，從而更精確地捕捉旋渦結(jié)構(gòu)，實(shí)現(xiàn)提高旋渦場(chǎng)的計(jì)算精度的目的。原始的渦量限制法是公知的，這里不再敘述。盡管該方法在捕捉不可壓縮流場(chǎng)中的接觸不連續(xù)的方面已經(jīng)取得了明顯的改進(jìn)效果，但存在以下缺陷：

1.對(duì)加入的渦量調(diào)整完全靠系數(shù)ε；

2.加入的渦量的空間離散精度是被限定的，無(wú)法進(jìn)一步提高；

3.源項(xiàng)對(duì)動(dòng)量方程的數(shù)值解的收斂的穩(wěn)定性影響是不確定的。

風(fēng)力葉片的轉(zhuǎn)速較低，馬赫數(shù)小于0.3，其周圍的流場(chǎng)可以被視為是一個(gè)不可壓縮(incompressible)流動(dòng)，同時(shí)是一種粘性流流動(dòng)。為了更加精確地模擬風(fēng)力葉片尾跡的流 動(dòng)—不可壓縮粘性流的旋渦運(yùn)動(dòng)，需要進(jìn)一步改進(jìn)對(duì)于流場(chǎng)中接觸不連續(xù)的捕捉機(jī)制。一種途徑是根據(jù)不可壓縮流的特點(diǎn)，改進(jìn)在流場(chǎng)中通過(guò)加入渦量來(lái)抵消數(shù)值耗散的內(nèi)在工作機(jī)理，通過(guò)保持渦量的精度，更精確地用模擬流場(chǎng)中的旋渦運(yùn)動(dòng)，形成一種新的，針對(duì)風(fēng)力葉片尾跡的數(shù)值模擬技術(shù)。該技術(shù)可以使渦量的空間離散具有高階精度的格式，同時(shí)該可以利用源項(xiàng)增進(jìn)收斂進(jìn)程的穩(wěn)定性。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明涉及計(jì)算流體力學(xué)的應(yīng)用領(lǐng)域中一種模擬風(fēng)力葉片尾跡的流動(dòng)(不可壓縮粘性旋流流動(dòng))的數(shù)值方法，具體是根據(jù)不可壓縮流的特點(diǎn)，在流場(chǎng)中通過(guò)加入兩種不同形式的渦量力來(lái)抵消數(shù)值耗散的數(shù)值方法。該方法可以使加入的渦量的空間離散具有高階精度的格式，同時(shí)還可以利用源項(xiàng)增進(jìn)數(shù)值解的收斂進(jìn)程的穩(wěn)定性。通過(guò)保持渦量的精度，更精確地模擬流場(chǎng)中的旋渦運(yùn)動(dòng)，是一種新的模擬風(fēng)力葉片尾跡的數(shù)值模擬技術(shù)。

首先寫出風(fēng)力葉片尾跡流動(dòng)的控制方程，即不可壓縮、粘性流流動(dòng)的控制方程，包括連續(xù)方程和動(dòng)量方程，分別為

式中，v是速度矢量，包含笛卡爾坐標(biāo)系下的三維i，j，k分量u,v,w；ρ、p、t、ν分別是密度、壓力，時(shí)間和粘性系數(shù)。式中算子表示為符號(hào)·代表內(nèi)積計(jì)算。動(dòng)量方程(2)的右邊是體積力項(xiàng)形式的源項(xiàng)按照渦量限制法的定義

式中，代表渦量，在直角坐標(biāo)系下有按照渦量的定義，

式中符號(hào)×代表叉乘運(yùn)算，代表渦量梯度變化的方向，即

其中，φ是渦量的模；是渦量的模的梯度的模，即

動(dòng)量方程(2)中的源項(xiàng)需要公式(3)中的乘以一個(gè)比例因子ε，則

其物理意義是指向渦量中心的力。的原始的渦量限制法中給出ε為常數(shù)，為0.01-0.05，其作用是用來(lái)調(diào)整渦量限制的大小。該方法形式簡(jiǎn)單，不需要加密計(jì)算網(wǎng)格，在不可壓縮旋流場(chǎng)的數(shù)值模擬中得到廣泛應(yīng)用。對(duì)該方法在隨后的改進(jìn)主要集中在參數(shù)ε的表達(dá)。

如公式(2)表示的，渦量限制法中在不可壓縮流動(dòng)的動(dòng)量方程的等號(hào)右邊加入了一個(gè)體積力項(xiàng)，它代表渦量在其變化梯度相反的方向的受力，如公式(7)所示。實(shí)際上，如果將公式(4)帶入公式示(7)并運(yùn)用算子的運(yùn)算規(guī)律可得

在上式中再次運(yùn)用算子的運(yùn)算規(guī)律可得

式中是拉普拉斯算子，定義為又根據(jù)公式(1)不可壓縮流的連續(xù)方程，將公式(9)代入公式(8)，可得

對(duì)于不可壓縮流，公式(10)與公式(7)完全等價(jià)?，F(xiàn)將公式(10)中的重新寫成如下形式，即

其中，

因?yàn)?imgfile="bda0000899215750000051.gif"wi="124"he="61"img-content="drawing"img-format="gif"orientation="portrait"inline="no"/>是螺旋度的定以，所以被稱作渦量在變化梯度方向的螺旋力；

因?yàn)槔绽顾阕拥暮纳⑻匦裕?imgfile="bda0000899215750000054.gif"wi="70"he="78"img-content="drawing"img-format="gif"orientation="portrait"inline="no"/>被稱作在渦量在變化梯度方向的粘性耗散力。

至此公式(2)中以動(dòng)量方程的源項(xiàng)形式表示的渦量在其變化梯度相反的方向的受力(渦量限制)限制被分解為兩部分，形成兩種力的形式。其中渦量在變化梯度方向的螺旋力與渦量的變化方向和渦量大小有關(guān)；而渦量在變化梯度方向的粘性耗散力僅與渦量的變化梯度方向有關(guān)。圖2給出兩種形式的力的形成原理和關(guān)系圖。圖中表明，由渦量和渦量梯度變化的方向向量構(gòu)成平面1，被稱作渦量作用平面。原始的渦量限制與該平面垂直，并與它的兩個(gè)分量，渦量在變化梯度方向的螺旋力和渦量變化梯度方向的粘性耗散力共同構(gòu)成平面2，被稱作渦量限制平面。平面1垂直與平面2。和的夾角被稱作螺旋角，該夾角的余弦即為渦量限制被分解后形成的兩種形式的力，和將在數(shù)值方法中起到不同的作用。

公式(12)和(13)中是用了兩個(gè)參數(shù)：ε1和ε2，分別作為不同的兩個(gè)力的放大系數(shù)。如果ε1和ε2相等，且等于公式(7)中的ε，則成為同向渦量限制(isotropicvorticityconfinement)，其效果等同于公式(7)表示的原始的渦量限制；而ε1和ε2不相等時(shí)，成為異向渦量限制(anisotropicvorticityconfinement)。而不可壓縮旋流場(chǎng)的數(shù)值模擬精度的提高可以通過(guò)異性渦量限制獲得。

作為源項(xiàng)形式的體積力對(duì)偏微分方程(2)的數(shù)值解的收斂過(guò)程有重要影響。其作用體現(xiàn)在加入源項(xiàng)后，方程(2)有可能變成具有剛性的(stiff)的性質(zhì)，致使一類依靠顯式時(shí)間迭代獲得數(shù)值解的收斂性降低，同時(shí)也降低數(shù)值解的穩(wěn)定性。因而，此類方程的求解為了保證時(shí)間精度，必須采用較小的時(shí)間步長(zhǎng)，結(jié)果造成計(jì)算速度變慢。衡量方程的剛性的指標(biāo)是源項(xiàng)的雅各比矩陣的特征值。如果雅各比矩陣的最大特征值和最小特征值的實(shí)部的比例小，則說(shuō)明偏微分方程的剛性小，數(shù)值解容易收斂，而且穩(wěn)定性強(qiáng)。

考慮兩種情況：第一，源項(xiàng)分別為即原始的渦量限制法中的體積力，源項(xiàng)的雅各比矩陣為特征值為i等于2或3；第二，源項(xiàng)為即渦量變化梯度方向的粘性耗散力，源項(xiàng)的雅各比矩陣為特征值為在兩種形式的體積力的情況下，源項(xiàng)的雅各比矩陣分別表示為

按照雅各比矩陣的特征值的求解方法，求解下式

和

其中i是單位矩陣，即可獲得兩種源項(xiàng)的特征值，而且有

上式中re表示取實(shí)部運(yùn)算。表明了以為源項(xiàng)時(shí)，偏微分方程的剛性比以為源項(xiàng)時(shí)小。因而，在使用一類顯式時(shí)間積分方法求解不可壓縮流的動(dòng)量方程(2)時(shí)會(huì)增強(qiáng)收斂性和數(shù)值穩(wěn)定性。所以，利用這個(gè)特性，可以將(渦量在變化梯度方向的螺旋力)從等號(hào)右邊移動(dòng)到等號(hào)左邊，等號(hào)右邊的源項(xiàng)僅保留(渦量變化梯度方向的粘性耗散力)。

如公式(12)所示，渦量在變化梯度方向的螺旋力與渦量的變化方向和渦量大小有關(guān)，數(shù)值解中需要高精度的空間離散，以提高模擬精度。源項(xiàng)中的被移動(dòng)等號(hào)左邊后，可以利用高斯定理，將控制體單元內(nèi)的體積力形式轉(zhuǎn)變?yōu)楸砻嫱康姆e分形式。有限體積法(finitevolumemethod)中，某一控制體單元，即某一計(jì)算網(wǎng)格的體積為ω，而其表面積為面積元為ds。根據(jù)高斯定理，可將體積積分變換為面積積分，則渦量在變化梯度方向的螺旋力從體積積分變?yōu)槊娣e積分的過(guò)程可表示為

其中計(jì)算網(wǎng)格邊界的單位向量。于是，在計(jì)算網(wǎng)格邊界上產(chǎn)生了一個(gè)由于渦量在變化梯度方向的螺旋力產(chǎn)生的向量，被稱作渦量的螺旋力向量

在數(shù)值解的求解過(guò)程中，可以對(duì)該向量在計(jì)算網(wǎng)格邊界上進(jìn)行高階精度的空間離散。由于該向量與渦量的變化方向和渦量大小有關(guān)，為進(jìn)一步提高了旋渦流場(chǎng)的空間模擬精度提供了可能性。例如可以同過(guò)流動(dòng)變量的高階插值獲得計(jì)算網(wǎng)格界面上用來(lái)計(jì)算向量的值。

此外，如果適當(dāng)增大ε1、減小ε2，意味著增大渦量在變化梯度方向的螺旋力、降低渦量變化梯度方向的粘性耗散力，即采用異向渦量限制，利用調(diào)節(jié)這兩個(gè)力的內(nèi)部關(guān)系的機(jī)制，可以進(jìn)一步調(diào)整旋渦運(yùn)動(dòng)中的接觸不連續(xù)界面的捕捉效果。

總之，本發(fā)明提出的提高不可壓縮旋流場(chǎng)的數(shù)值模擬精度的數(shù)值方法包括四個(gè)技術(shù)元素，分別為：

1.將渦量限制法中的源項(xiàng)分解稱為兩個(gè)力，渦量在變化梯度方向的螺旋力和渦量變化梯度方向的粘性耗散力，并將渦量在變化梯度方向的螺旋力移動(dòng)到不可壓縮流的動(dòng)量方程的等號(hào)左變；

2.通過(guò)高斯定理將計(jì)算網(wǎng)格內(nèi)渦量在變化梯度方向的螺旋力轉(zhuǎn)化為計(jì)算網(wǎng)格邊界上的上的力，進(jìn)行高精度的空間離散，按照通量進(jìn)行計(jì)算；

3.源項(xiàng)保留渦量變化梯度方向的粘性耗散力用來(lái)提高數(shù)值解的收斂性和穩(wěn)定性；

4.渦量在變化梯度方向的螺旋力和渦量變化梯度方向的粘性耗散力采用不同的放大系數(shù)，即采用異向的渦量限制。

本發(fā)明提出的提高直升機(jī)旋翼尾跡不可壓縮旋流場(chǎng)的數(shù)值模擬精度的數(shù)值方法是根據(jù)不可壓縮流的特點(diǎn)，通過(guò)加入兩種不同形式的力來(lái)改進(jìn)在流場(chǎng)中抵消數(shù)值耗散的內(nèi)在工作機(jī)制，使計(jì)算網(wǎng)格內(nèi)的渦量在變化梯度方向的螺旋力轉(zhuǎn)化為計(jì)算網(wǎng)格邊界上的上的力，可以使其空間離散具有高階精度的格式；同時(shí)源項(xiàng)保留渦量變化梯度方向的粘性耗散力，用來(lái)提高數(shù)值解的收斂性和穩(wěn)定性。這兩個(gè)力采用不同的放大系數(shù)，可以進(jìn)一步保持渦量的精度，更精確地模擬風(fēng)力葉片尾跡中的旋渦運(yùn)動(dòng)。

附圖說(shuō)明

圖1風(fēng)力葉片翼尾跡中的旋渦體系。圖中，10無(wú)粘翼尖渦、11粘性脫落渦。

圖2渦量在變化梯度方向的螺旋力和渦量變化梯度方向的粘性耗散力形成原理和關(guān)系圖

圖中，7渦量作用平面、8螺旋角、9渦量限制平面。

圖3葉片的外形和計(jì)算域。

圖中，1遠(yuǎn)場(chǎng)、2尾跡、3近場(chǎng)、4尾跡、5機(jī)翼、6近場(chǎng)。

圖4模擬風(fēng)力葉片尾跡的流程圖。

具體實(shí)施方式

以下以一個(gè)具體實(shí)施方案進(jìn)一步說(shuō)明本發(fā)明提出的這種數(shù)值模擬技術(shù)?？梢杂胒ortran90計(jì)算機(jī)高級(jí)程序語(yǔ)言實(shí)現(xiàn)，并通過(guò)計(jì)算機(jī)運(yùn)行來(lái)模擬橫截面為naca0012形狀的三維旋轉(zhuǎn)葉片在高攻角飛行狀態(tài)下的尾跡流場(chǎng)。

具體實(shí)施方案中，風(fēng)力發(fā)電機(jī)采用四個(gè)旋轉(zhuǎn)的葉片，葉片展翼比為7。如圖3所示，剖面為naca0012的三維葉片和計(jì)算域。計(jì)算域取四分之一象限，包含一個(gè)葉片，z-軸為旋轉(zhuǎn)軸。計(jì)算網(wǎng)格采用多塊結(jié)構(gòu)化六面體網(wǎng)格。根據(jù)風(fēng)力葉片流場(chǎng)的特點(diǎn)將計(jì)算域分為遠(yuǎn)場(chǎng)(1)、尾跡(2)、近場(chǎng)(3)，五個(gè)塊。機(jī)翼周圍采用o-型網(wǎng)格，圍繞機(jī)翼剖面方向192個(gè)；法線方向48個(gè)；翼展方向32個(gè)。數(shù)值模擬的空間離散采用有限體積法(finitevolumemethod)，每個(gè)計(jì)算網(wǎng)格成為控制單元，流動(dòng)變量在控制單元的中心處。數(shù)值模擬的飛行條件是：

來(lái)流速度：57m/s；

大氣壓力：100kpa；

大氣密度：1.2kg/m³；

攻角：20^o。

本發(fā)明中要求加上在計(jì)算網(wǎng)格邊界上產(chǎn)生了一個(gè)由于渦量在變化梯度方向的螺旋力產(chǎn)生的向量意味著在動(dòng)量方程左邊的對(duì)流項(xiàng)后面加上這個(gè)向量，在等號(hào)右邊，即源項(xiàng)的位置加上渦量變化梯度方向的粘性耗散力按照守恒定律重新寫出動(dòng)量方程(2)的三維積分形式，

其中，守恒變量對(duì)流項(xiàng)向量粘性項(xiàng)向量和渦量的螺旋力向量分別為

上式中，密度ρ為一固定值；v代表速度不變量，是網(wǎng)格邊界的外法線向量與速度向量的點(diǎn)乘，即

正粘性應(yīng)力項(xiàng)和剪切粘性應(yīng)力項(xiàng)非別為

可以將對(duì)流項(xiàng)向量和渦量的螺旋力向量寫成展開形式，

對(duì)流項(xiàng)經(jīng)壓力分離后，可以寫成不計(jì)壓力的對(duì)流項(xiàng)向量和壓力向量之和，顯然有，

其中，

進(jìn)一步可以寫出本發(fā)明提出的渦量在變化梯度方向的螺旋力向量和源項(xiàng)形式的渦量變化梯度方向的粘性耗散力的具體形式分別為，

其中，渦量的模φ和渦量的模的梯度的模由公式(5)和(6)給出。公式(28)和(29)中涉及的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)項(xiàng)均可利用格林公式求取或者直接按照泰勒展開求取，其過(guò)程是公知的，不再敘述。

圖4給出了具體實(shí)施例中采用本發(fā)明提出的模擬直升機(jī)旋翼尾跡的數(shù)值方法的流程圖。以上公式(20)至(29)體現(xiàn)了這個(gè)流程中直至生成模擬直升機(jī)旋翼尾跡的不可壓縮流控制方程，即包含渦量在變化梯度方向的螺旋力向量和渦量變化梯度方向的粘性耗散力的不可壓縮粘性流的動(dòng)量方程。下面的步驟是求解這個(gè)方程。

以下是在同位網(wǎng)格(collocatedgrid)中運(yùn)用標(biāo)準(zhǔn)的基于壓力的方法(pressure-basedmethod)求解包含渦量在變化梯度方向的螺旋力向量和渦量變化梯度方向的粘性耗散力的不可壓縮粘性流控制方程的過(guò)程。

首先將公式(20)寫成離散形式，

其中，ωi,j,k代表空間任一離散控制體體積；△t代表積分時(shí)間步長(zhǎng)；代表流場(chǎng)更新后速度；△sm代表控制體任一表面面積；nf代表控制體表面?zhèn)€數(shù)，三維情況下為6。重新整理右手項(xiàng)，有

而使用一個(gè)假設(shè)的初始速度計(jì)算一個(gè)不計(jì)壓力項(xiàng)的右手項(xiàng)即

顯然有，

不計(jì)壓力項(xiàng)時(shí)候，可以從公式(31)獲得這個(gè)個(gè)初始速度場(chǎng)

其中需要在計(jì)算網(wǎng)格邊界m上的不計(jì)壓力的對(duì)流項(xiàng)向量渦量在變化梯度方向的螺旋力向量以及在計(jì)算網(wǎng)格內(nèi)的渦量變化梯度方向的粘性耗散力再帶入公式(33)以獲得不計(jì)壓力的右手項(xiàng)

在按照公式(28)計(jì)算網(wǎng)格邊界m處的渦量在變化梯度方向的螺旋力向量的時(shí)候，對(duì)于其中的一階導(dǎo)數(shù)項(xiàng)采用高精度離散方法，例如在i-方向的三階空間離散表示為，

網(wǎng)格中心點(diǎn)的初始速度在計(jì)算網(wǎng)格邊界上利用公知的muscl方法進(jìn)行二階精度插值，可獲得在網(wǎng)格邊界m上的速度因而，根據(jù)公式(35)不計(jì)壓力項(xiàng)的情況下，在邊界m處的速度表示為

其中，是包含網(wǎng)格邊界m的輔助網(wǎng)格的體積。在邊界m處，考慮壓力項(xiàng)的速度為

其中，是公式(32)在網(wǎng)格邊界m處的值，任然用muscl方法進(jìn)行二階插值獲得。整理公式(37)和(38)，有

考慮壓力項(xiàng)的速度應(yīng)該滿足不可壓縮流的連續(xù)方程(1)，即將連續(xù)方程(1) 寫成離散形式，

將公式(39)帶入(40)，有

求解方程(41)，即可獲得t+△t時(shí)刻，流場(chǎng)的計(jì)算網(wǎng)格邊界m處的速度和壓力值pm，這兩個(gè)值可以通過(guò)集合平均處理，轉(zhuǎn)化為計(jì)算網(wǎng)格中心處的值。

以上是本發(fā)明提出的一種模擬直升機(jī)旋翼尾跡的數(shù)值方法的論述。