本發(fā)明涉及計算流體力學應用領域中的一種數(shù)值方法，具體是一種模擬橋墩周期受力的數(shù)值方法。

背景技術：

當河流經(jīng)過大型橋墩物時，因為橋墩的迎風面和背風面的壓力差的存在，在橋墩的后方會產(chǎn)生水流的漩渦。漩渦隨著流水向下方繼續(xù)流動，形成一系列脫落渦。這一系列脫落渦的流動也被稱為卡門渦街，是一個充滿了旋渦的旋流場，河流水流流場可以被視為是一個不可壓縮流動。實際上，脫落渦會對橋墩周期地產(chǎn)生誘導作用力，嚴重的會使建筑物晃動，如果渦脫落的周期和橋墩建筑的自激頻率相接近，建筑物的晃動會被放大，甚至發(fā)生毀壞。

對于橋墩周期受力有多種研究手段。其中，計算機的數(shù)值模擬技術有著重要地位，是計算流體力學在該領域的一個擴展應用。計算流體力學綜合了流體力學、應用數(shù)學、計算機科學，是一門應用性極強的學科。流體力學問題的數(shù)值模擬以其低成本、直觀性強的優(yōu)勢，在流體流動的機理探索、工業(yè)產(chǎn)品設計等各個相關領域占據(jù)重要地位。橋墩后方脫落渦的數(shù)值模擬面臨的最大的問題即是如何提高數(shù)值模擬的精度、降低誤差，忠實地表現(xiàn)橋墩后方脫落渦流動的特性。

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

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

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

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

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

2.加入的渦量的空間離散精度是被限定的，無法進一步提高；

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

為了更加精確地模擬橋墩后方脫落渦流動—不可壓縮無粘流的旋渦運動，需要進一步改進對于流場中接觸不連續(xù)的捕捉機制。一種途徑是根據(jù)不可壓縮流的特點，改進在流場中通過加入渦量來抵消數(shù)值耗散的內(nèi)在工作機理，通過保持渦量的精度，更精確地用模擬流場中的旋渦運動，形成一種新的，針對橋墩后方的渦脫落長生的周期受力的數(shù)值模擬技術。該技術可以使渦量的空間離散具有高階精度的格式，同時該可以利用源項增進收斂進程的穩(wěn)定性。

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

本發(fā)明涉及計算流體力學的應用領域中一種模擬橋墩周期受力(不可壓縮無粘旋流流動)的數(shù)值方法，具體是根據(jù)不可壓縮流的特點，在流場中通過加入兩種不同形式的渦量力來抵消數(shù)值耗散的數(shù)值方法。該方法可以使加入的渦量的空間離散具有高階精度的格式，同時還可以利用源項增進數(shù)值解的收斂進程的穩(wěn)定性。通過保持渦量的精度，更精確地用模擬流場中的旋渦運動，是一種新的不可壓縮流的旋渦運動的數(shù)值模擬技術。

首先寫出不可壓縮、無粘流流動的控制方程，包括連續(xù)方程和動量方程，分別為

式中，v是速度矢量，包含直角坐標系下的三維i，j，k分量u，v，w；ρ、p、t分別是密度、壓力和時間。式中算子表示為符號·代表內(nèi)積計算。動量方程(2)的右邊是體積力項形式的源項按照渦量限制法中的定義

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

式中符號×代表叉乘運算，代表渦量梯度變化的方向，即

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

動量方程(2)中的源項需要公式(3)中的乘以一個比例因子ε，則

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

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

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

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

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

其中，

因為是螺旋度的定以，所以被稱作渦量在變化梯度方向的螺旋力；

因為拉普拉斯算子的耗散特性，被稱作在渦量在變化梯度方向的粘性耗散力。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說明

圖1橋墩后方的流動模式。圖中，4橋墩、5周期性渦脫落、6來流方向。

圖2渦量在變化梯度方向的螺旋力和渦量變化梯度方向的粘性耗散力形成原理和關系圖圖中，1渦量作用平面、2螺旋角、3渦量限制平面。

圖3模擬橋墩周期受理的流程圖

具體實施方式

以下以一個具體實施方案進一步說明本發(fā)明提出的模擬橋墩周期手里的數(shù)值方法。這種數(shù)值模擬技術可以用fortran90計算機高級程序語言實現(xiàn)，并通過計算機運行來模擬三維流動的周期性渦脫落流場。

具體實施方案中，橋墩其周圍的流場是以旋渦運動為主的流場。計算網(wǎng)格采用多塊結構化適體六面體網(wǎng)格，橋墩周圍采用o-型網(wǎng)格，圍繞橋墩剖面方向192個網(wǎng)格；法線方向68個網(wǎng)格；翼展方向96個網(wǎng)格。數(shù)值模擬的空間離散采用有限體積法(finitevolumemethod)，每個計算網(wǎng)格成為控制單元，流動變量在控制單元的中心。

數(shù)值模擬的條件是：

來流速度：57m/s；

流體密度：1000kg/m^3；

不可壓縮無粘流的動量方程，即方程(2)中不計粘性項。本發(fā)明中要求加上在計算網(wǎng)格邊界上產(chǎn)生了一個由于渦量在變化梯度方向的螺旋力產(chǎn)生的向量意味著在動量方程左邊的對流項后面加上這個向量，在等號右邊，即源項的位置加上渦量變化梯度方向的粘性耗散力按照守恒定律，將動量方程(2)簡化寫成適合河流流動模擬的粘性不可壓縮淺水方程的形式，并寫出積分形式，

其中，守恒變量對流項向量和渦量的螺旋力向量(見公式(19))分別為

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

可以將對流項向量寫成展開形式，

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

其中，

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

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

圖3給出了本實施例模擬橋墩后方脫落渦的流程圖。公式(20)至(29)體現(xiàn)了這個流程中直至生成模擬橋墩后方脫落渦的不可壓縮流控制方程，即包含渦量在變化梯度方向的螺旋力向量和渦量變化梯度方向的粘性耗散力的不可壓縮無粘流的動量方程。下面的步驟是求解這個方程。

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

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

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

而使用一個假設的初始速度計算一個不計壓力項的右手項即

顯然有，

不計壓力項時候，可以從公式(31)獲得這個初始速度場

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

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

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

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

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

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

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

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

以上是本發(fā)明提出的一種模擬橋墩周期性受力流動的數(shù)值方法的論述。