openfoam-解析流體技術(shù)二icoFoam詳解

icoFoamicoFoam采用的是標(biāo)準(zhǔn)瞬態(tài)PISO算法，本文檔從最基本的NS方程離散開始推導(dǎo)，并和OpenFOAM植入的代碼相對應(yīng)以供理解瞬態(tài)PISO算法如何在OpenFOAM中實現(xiàn)以及NS方程在OpenFOAM中的離散。對于穩(wěn)態(tài)的標(biāo)準(zhǔn)SIMPLE算法以及PIMPLE算法的對比請參考《簡單易懂的PISO，SIMPLE算法對比》NS方程離我們有動量方程

u

u n我們令當(dāng)前時間步為n，下一個時間步為n+1，預(yù)測的速度為urRhieandChow插值精神，首先忽略壓力梯度項uN-S方程進(jìn) ndVdt tt 其中

urSur

unu

方程1.3代表一個向量，速度的下標(biāo)n表示臨點的速度，上標(biāo)n表示當(dāng)前的時間步。進(jìn)其可以 f ur f pt

FnururS -S（bA詳見下述。批注[W用3]:HbyA因此需要下文的對HbyA批注[W用3]:HbyA因此需要下文的對HbyAAu Au prpn Au Ap n- prpAu+Au-E=-ppn r nA Au- p n ru其中由于為面預(yù)測速度，我們需要進(jìn)行插值，在此步我們可以引入各種格式。假設(shè)uf有中心差分uruur

代入

uru urun unurVn pt Ff nVn Fn

2

Su

Sunpun t d 我們可以把上式簡化為Aur+Aun-En=

d

批注批注[W用1]:+fvm::div(phi,-fvm::laplacian(nu,)p值得注意的是，在整個時間步內(nèi)

n UEqn.A()即????保持不[W用2]:solve(UEqnUEqn.H()中的系數(shù)[W用2]:solve(UEqn至此，如果我們加入壓力梯度的顯性項。即其為速度預(yù)測方程。求解此方程后（solve…），我們有預(yù)測速度??????第一次修正步（由????計算????再計算

首先我們通過預(yù)測的????組建HbyAr1Aur+Er n p??????為第一次修正速度（應(yīng)符合連續(xù)性方程）????為第一次修正后的壓力。那么有：

purrHbyArp

我們令（依據(jù)連續(xù)性方程

uprrurrS p依據(jù)公式3.3，我們對面值進(jìn)行插值A(chǔ)urrHbyAr1 Apf p,f

r 整理有

HbyAf fpS p rf fp fA參考icoFoam中phiHbyA的源代碼

p iHArHbyAr pr hiyAsiiHUrrpp求解后我們有壓力場????。這一步我們稱之為第一步修正。利用求解的????回代到公式3.3。我們有滿足連續(xù)性的速度場??????。但我們可以發(fā)現(xiàn)HbyA項并沒有進(jìn)行更新且有所滯后，但實際上??????應(yīng)該由更新的HbyA和更新的壓力場進(jìn)行組建，因此我們進(jìn)行下一個時間步更HbyA第二次修正步（由??????計算??????再計算參考第一次修正步，直至速度方程左右的誤差可以忽評經(jīng)過Issa證明，2步的PISO循環(huán)已經(jīng)足夠。其速度場的誤差為deltaT4階無窮小，壓力場的誤差為deltaT3階無窮小。deltaT越小，相對于方程本身的離散誤差來說，PISO循環(huán)的誤差越小。PISO循環(huán)數(shù)越多，deltaT的誤差階數(shù)越高icoFoamicoFoam采用的是標(biāo)準(zhǔn)瞬態(tài)PISO算法，本文檔從最基本的NS方程離散開始推導(dǎo)，并和OpenFOAM植入的代碼相對應(yīng)以供理解瞬態(tài)PISO算法如何在OpenFOAM中實現(xiàn)以及NS方程在OpenFOAM中的離散。對于穩(wěn)態(tài)的標(biāo)準(zhǔn)SIMPLE算法以及PIMPLE算法的對比請參考《簡單易懂的PISO，SIMPLE算法對比》NS方程離我們有動量方程

u

u n我們令當(dāng)前時間步為n，下一個時間步為n+1，預(yù)測的速度為urRhieandChow插值精神，首先忽略壓力梯度項uN-S方程進(jìn) ndVdt

ff

udVdt udSdtuSt=

S 其中

urSur

unu

方程1.3代表一個向量，速度的下標(biāo)n表示臨點的速度，上標(biāo)n表示當(dāng)前的時間步。進(jìn)其可以 f ur f pt

FnururS -S（bA詳見下述。批注[W用3]:HbyA因此需要下文的對HbyA批注[W用3]:HbyA因此需要下文的對HbyAAu Au prpn Au Ap n- prpAu+Au-E=-ppn r nA Au- p n ru其中由于為面預(yù)測速度，我們需要進(jìn)行插值，在此步我們可以引入各種格式。假設(shè)uf有中心差分uruur

代入

uru urun unurpnptp pFfpnVnV

ndp Fn

p 2

Su

Sunpun t d 我們可以把上式簡化為Aur+Aun-En=

d

批注批注[W用1]:+fvm::div(phi,-fvm::laplacian(nu,)p值得注意的是，在整個時間步

n UEqn.A()即????保持不UEqn.H()中的系數(shù)即????保持不至此，如果我們加入壓力梯度的顯性項。即[W[W用2]:solve(UEqn其為速度預(yù)測方程。求解此方程后（solve…），我們有預(yù)測速度??????第一次修正步（由????計算????再計算首先我們通過預(yù)測的????組建HbyAr1Aur+Er

n p??????為第一次修正速度（應(yīng)符合連續(xù)性方程）????為第一次修正后的壓力。那么有：

purrHbyArp

我們令（依據(jù)連續(xù)性方程

uprrurrS p依據(jù)公式3.3，我們對面值進(jìn)行插值urrHbyAr1 Apf p,fA

r 整理有

HbyAf fpS p rf fp fA參考icoFoam中phiHbyA的源代碼

p iHArHbyAr pr hiyAfsiiHUrr求解后我們有壓力場????。這一步我們稱之為第一步修正。利用求解的????回代到公式3.3。我們有滿足連續(xù)性的速度場??????。但我們可以發(fā)現(xiàn)HbyA項并沒有進(jìn)行更新且有所滯后，但實際上??????應(yīng)該由更新的HbyA和更新的壓力場進(jìn)行組建，因此我們進(jìn)行下一個時間步更HbyA第二次修正步（由??????計算??????再計算參考第一次修正步，直至速度方程左右的誤差可以忽評經(jīng)過Issa證明，2步的PISO循環(huán)已經(jīng)足夠。其速度場的誤差為deltaT4階無窮小，壓力場的誤差為deltaT3階無窮小。deltaT越小，相對于方程本身的離散誤差來說，PISO循環(huán)的誤差越小。PISO循環(huán)數(shù)越多，deltaT的誤差階數(shù)越高,scalarTransportFoamOpenFOAM之中最基本的求解擴(kuò)散率0/U中設(shè)置不隨時間變化的速度場。所以，這個求解器不能用于求解流場#include"fvCFD.H在OpenFOAM的求解器中，涉及到的構(gòu)建時間、組建矩陣、有限體積離散、組建網(wǎng)格、量綱設(shè)置等大量內(nèi)#include"fvIOoptionList.H"http://通過fvOption來修正源項，主要用于在運行時添加多孔介質(zhì)、MRF多重參考系、顯性源項、熱"http://*************************************intmain(intargc,char{ #include"createMesh.H"http://創(chuàng)建網(wǎng)格對象，必備頭文件。請忽略#include"createFields.H"http://從求解器根 #include"createFvOptions.H"http://創(chuàng)建源項，無需源項可刪除。請忽略//***********************************InfonCalculatingscalartransport\nendl;在OpenFOAM中采用Info來在終端輸出信息，內(nèi)為需要輸出的信息#include"CourantNo.H"http://計算庫郎數(shù)，詳見《LTS局部時間步，庫郎數(shù)分析》，請忽略{while(simple.correctNonOrthogonal())//是否進(jìn)行非正交修正。如果再fvSolution字典文件中設(shè)置為0，那么求解下面的《OpenFOAM中的有限{)//-fvm::laplacian(DT,T)//fvOptions(T)//源項，若不需要源項請忽略};//}return}createFields.H。幾乎所有求解器都會有這樣一個.H文件用來創(chuàng)建每volScalarFieldT//volScalarField表示創(chuàng)建標(biāo)量場，本例T為一種典型的定義，自定義求解器的時候可以直接并改(IOobetpnFAM采用IOobectIOobet(runTime.timeName(),//在運行時mesh,//于網(wǎng)//Info<<"ReadingfieldU\n"<<Info<<"ReadingtransportProperties\n"<<IOdictionarytransportProperties//IOdictionary即為創(chuàng)建字典文件，用戶可以在transportProperties中參數(shù)，其他參見上(IOobject::MUST_READ_IF_MODIFIED在字典文件被更改的時候進(jìn)行IOobject::NO_WRITE從不輸出，字典文件不需要輸出) "laplacianFoamOpenFOAM中最簡單的求解器之一，它可以作為用戶了解OpenFOAM求解器的一個起點。本求解器的?T??2(?????)= 碼和scalarTransportFoam相同，簡要說明不再贅述) #include"fvCFD.H"http://*************************************intmain(intargc,char{"#include"createTimeH創(chuàng)建時間#include"createMesh.H創(chuàng)建網(wǎng)格#include"createFields.H"http://創(chuàng)建場//***********************************//Info<<"\nCalculatingtemperaturedistribution\n"<<endl;({{}"<<"ClockTime="<<runTime.elapsedClockTime()<<"}return}if(runTime.outputTime())//如果輸出的時間步和運行時間步相同的時候，為真，進(jìn)行括號內(nèi)的程{《OpenFOAM編程指南ponent(vector::X)//將gradT矢量的x分量賦值給此gradTx標(biāo)量;//}deltaT，這會使質(zhì)量失衡（瞬態(tài)問題每個時間步的質(zhì)量守恒CoNum=sumphi=∑|????? CoNum= Co=

加和了兩次，因此公式1.2需要除2。scalarFieldsumPhi(alphaCoNum=sumPhi=(??1?0.01)(0.99???1)∑|?????????|,0.01<??1< sumPhi≈??1??2∑|?????????|,0.01<??1< scalarmaxDeltaTFact=min(maxCo/(CoNum+SMALL),maxAlphaCo/(alphaCoNum+scalardeltaTFact=min(min(maxDeltaTFact,1.0+0.1*maxDeltaTFact),1.2);()

的庫郎數(shù)為止。然后就是interFoam方程循環(huán)。LTSinterFoam求解deltaT為非均一的場。所需要的各種量1/dimensionedScalar("maxDeltaT",dimTime,maxDeltaT),?t

∑|?????????????|2???????????????alphainterFoam中，我們采用歐拉法進(jìn)行時間步離散。但是在LTSinterFoam中，我們對時間項采用局部時間步離散；正后的alpha場；1fvc::smooth,fvc::spread,fvc::sweepLTS中使用其他場參數(shù)對局部時間步場進(jìn)行光順，參LTS代碼請參考之前的文檔。本文檔涉及一些C++類相關(guān)內(nèi)容，有關(guān)C++類的內(nèi)容推薦《C++名字為fluid，他需要傳入?yún)?shù)U和phi。略去構(gòu)造機(jī)理。singlePhaseTransportModel.C1代碼，有下述成員函數(shù)：=0，表明這是虛基類（抽象基類viscosityModel.H相同的文件下，1 2 至此singlePhaseTransportModel類，名在執(zhí)行fluid.correct()的時候進(jìn)入：singlePhaseTransportModel.correct()，viscosityModel.correct()，由于viscosityModel為虛基類，進(jìn)入具體類：nu_= returnmax(…,…)本語句返回計算粘度值注意到在牛頓流體中我們的速度方程其中拉斯項為laplacian(nu,U)，在肥牛臀流體中我場。方法和fluid.correct()一致。從略。其他內(nèi)容和icoFoam一致。分析，一步一步推導(dǎo)VOF模型。VOF定義一個示蹤函數(shù)11，流20，零到一之間的值則表示自由表面?？梢?，VOF的界面并不是直接算出來的，而是簡介的計算每個網(wǎng)格內(nèi)的相場值而后處和界面壓力降?P均衡。這個界面壓力降主要取決于表面張力和曲面弧度。?生的力為(??1??2)????1????2。表面張力作用在曲面微元的四個邊上，合力向上，考慮C????2在垂2

1

p1p2C

R2nx(xdx)nx(x)nx(xdx)xnx(xdx)sind xdxC=σ，以及11

nd n n=nn(11 ppn= 在VOF中，我們用表示體積分?jǐn)?shù)，有α=

100<??<1n|n=-(|p

c0c

c c pp()( || u

u(uu)pf (1 (1

uu

1

1tu u1D1D12212D Du uuT u

(uu)puuTg-

(1

1(1) u 續(xù)性方程中的對流項中的所有項中，因為我們要保證相分?jǐn)?shù)的有界和相界面的，在過去而這種方法雖然保持了守恒性但是卻會是的。另法是：添加一個對流項以壓縮界u1u interFoamtwoPhaseEulerFoam中植入上述方程進(jìn)行求解。經(jīng)測試mixerVesselAMI算例，和使用MULES求解的結(jié)果一致。u1u 1Dropimpactontoaliquidlayeroffinitethickness:Dynamicsof 2ANewApproachtoVOF-basedInterfaceCapturingMethods pressibleandCompressible

ucu

u|uccmaxu|uu

c

u,maxu|u1mincu,maxu

|u(uu)pfg ： pgprghpghu(uu) gh

uu

1

u

u

|(uu)uuT gh (1 (1 OpenFOAM中存在了大量的類，他們主要位于安裝 下的src文件夾中，對計算機(jī)語言不我們再看OpenFOAM2.3.0中的一段描述：fluidsinglePhaseTransportModel類，我們先不管他fluid在這里被調(diào)用。我們不考慮他是如何實現(xiàn)的。可以這樣來想：在非利用當(dāng)前的速度來求解粘度。上文的fluid.correct();就是實現(xiàn)的這個功能。while{Info<<"Time="<<runTime.timeName()<<nl<<#include"CourantNo.H"+fvm::div(phi,(fvc::grad(U)&定義一個粘度場，powerLaw模型需要的k寫，會大大增加了代碼。這不符合OpenFOAM的代碼規(guī)范Pimplepiso算法進(jìn)行時間步進(jìn)，時間步內(nèi)若有需要對速度矩陣系數(shù)、非最終迭代步的壓力場進(jìn)行低松弛求解提高穩(wěn)定性。pimpleFoam中的大部分peapooam邊界條件。雷同的代碼不再注釋，請參考《simpleFoam解析》以及《pisoFoam解析》另外提醒的是，幾個基本的OpenFOAM求解器（laplacianFoam,potentialFoam,scalarTransportFoam,icoFoam,pisoFoam,simpleFoam,pimpleFoam）已經(jīng)剖析完成。其余的OpenFOAMpimple、simpleLTS局部時間步框架下的求解器。CFD技術(shù)居多。其余的求解器不再代碼進(jìn)行通篇分析，只分析#include#include"turbulenceModel.H"""http://#include"fvIOoptionList.H"#include"IOMRFZoneList.H"#include在在pimpleFoam調(diào)用setSnGrad模板函數(shù)，計算表達(dá)式表達(dá)式計算之后，調(diào)動setSnGrad內(nèi)的函數(shù)updateCoeffs()更新邊界條件通過fixedFluxPressureFvPatchScalarField.C函數(shù)內(nèi)的updateCoeffs()對場的gradient()參考《icoFoam》，我們有邊界處的速度uHbyA1 AuSHbyAS1pA pSpn

AHbyASuS SpnAHbyASuSSffintmain(intargc,char{#include"setRootCase.H"#include"createTimeH"#include"createMesh.H"#include"createFields.H"#include"createFvOptions.H"{#include"CourantNo.H"#includesetDeltaT.H通過庫郎數(shù)設(shè)定Info<<"Time="<<runTime.timeName()<<nl<<{#include{#include}if{}}<<"ClockTime="<<runTime.elapsedClockTime()<<"<<nl<<}return0;}//SolvetheMomentumif(pimple.momentumPredictor())//此處如果開啟動量預(yù)測，求解下述方程，默認(rèn)為關(guān)閉。具體情況需要具體分析，Henry{}1《LTSHbyA=rAU*UEqn().H();if(pimple.nCorrPISO()<={}adjustPhi(phiHbyA,U,p);//UpdatethefixedFluxPressureBCstoensureflux);//此處即為壓力中通過setSnGrad設(shè)定壓力邊界//Non-orthogonalpressurecorrectorwhile{//PressurecorrectorfvScalarMatrixpEqnpEqn.setReference(pRefCell,pRefValue);if(pimple.finalNonOrthogonalIter()){phi=phiHbyA-}}#include//Explicitlyrelaxpressureformomentum();//U=HbyA-rAU*fvc::grad(p);pimpleFoampisoFoamsimpleFoam對于不可壓縮層流我們有NS味著需要用已知時間步的速度來計算Ap和An。這進(jìn)而又分為兩個情況：針對以上兩個問題的穩(wěn)態(tài)瞬態(tài)的不同，于是有了非迭代的瞬態(tài) SIMPLEPISO也可以用于穩(wěn)態(tài)并且SIMPLE也可以用于瞬態(tài)。其區(qū)別就在于舉例：如果我們考慮SIMPLE和PISO的話，會有很大的不同。其在于瞬態(tài)的PISO是iterative量的時間和資源。瞬態(tài)SIMPLE也具有以下特點：SIMPLE本身側(cè)重于對非線性化的處理。在大時間步下，速度的線性化滯后變的時候，瞬態(tài)SIMPLE算在每個時間步內(nèi)進(jìn)行大量迭代來計算。PIMPLE可以使用松弛技術(shù)。值得注意的是在源代碼中pisoFoam依然存在(),（SIMPLE相同PISO的多次壓力修正以對速度壓力耦合問題壓力方程要實現(xiàn)的關(guān)鍵問題就是從連續(xù)性方程衍生出壓力方程以使得每個網(wǎng)格質(zhì)量守恒PISOSIMPLE1PISO的區(qū)別參見《東岳流體技術(shù)文檔（二）——icoFoam2《東岳流體技術(shù)文檔（二）——icoFoam詳解》內(nèi)有從連續(xù)性方程推導(dǎo)壓力泊松方程的詳細(xì)步驟，此處是動量預(yù)測：依據(jù)p0預(yù)測速度進(jìn)入當(dāng)前時間步，我們有壓力p1,注意每個時間步內(nèi)的Ueqn.H()是變化組建速度方程的Ueqn.A和否否收收斂候引入了滯后，因此判定收得HbyA后清理速度方程便于內(nèi)存管依據(jù)p0計算速度u1，并對速度方程低松進(jìn)入當(dāng)前迭代步，我們有是否否本文檔剖析pisoFoam的源代碼，大量了《icoFoam詳解》中NS方程離散的相關(guān)內(nèi)容。需要注意的是icoFoampisoFoamPISONS方程求解，icoFoam只能用來求解層流。pisoFoam可以求解湍流以及層流。更加具有普適性的求解器為概念。因此可以用《icoFoam詳解》作為補(bǔ)充閱讀。本文檔初級代碼分析將省略。請參見《laplacianFoam,scalarTransportFoam流模型。OpenFOAM-2.3.0版本之前尚不支持fvOptions源項。首先我們分析createFields.H：Info<<"Readingfieldp\n"<<endl;volScalarFieldp(Info<<"ReadingfieldU\n"<<endl;volVectorFieldU( labelpRefCell=0;1,);//需要構(gòu)造的類名 構(gòu)造的類的具體名稱傳入的參此行代碼的意思即為將UphisinglePhaseTransportModel類，其名稱為laminarTransport1請參考《potentialFoam 一個成指針來使用，有關(guān)如何使用auto_ptr建議閱讀C++專業(yè)書籍。我們在此進(jìn)行簡單介紹。本行代碼規(guī)范autoPtr智能指

類名類名 類之中的New函個名為turbulence的指針。對于初學(xué)者理解此意義即可。下面我們進(jìn)入pisoFoam主函數(shù)#include""http://"http://*************************************intmain(intargc,char{#include"createTimeH"#include"createMesh.H"#include"createFields.H"為//***********************************Info<<"\nStartingtimeloop\n"<<endl;while(runTime.loop()){Info<<"Time="<<runTime.timeName()<<nl<<"=//constintnOuterCorrpisoDict.lookupOrDefault<int>("nOuterCorrectors1);//PISO并沒有外循環(huán)步，此步用來和PIMPLE算法相統(tǒng)一，對于沒有用到的，我們會在編譯中遇到提示：=//constintnNonOrthCorrpisoDict.lookupOrDefault<int>("nNonOrthogonalCorrectors0);//PISO非矩形修正步3，默認(rèn)為0//constboolmomentumPredictorpisoDict.lookupOrDefault("momentumPredictor"true);//動量預(yù)測開=",23《potentialFoam{//MomentumUEqnif{}kEpsilon湍流模型，我們在kEpsilon會找到下面的函數(shù)：{}

代碼中我們使用了UEqn.relax()。松弛技術(shù)在穩(wěn)態(tài)求解中我們經(jīng)常遇到，在OpenFOAM中有倆種松弛技術(shù)：對于瞬態(tài)求解器pisoFoam動量預(yù)測并不是必須的，Henry表示在諾數(shù)流動的時候動量預(yù)測對收斂不利，在某//---PISOfor(intcorr=0;corr<nCorr;corr++)//此項為PISO的壓力泊松方程修正{surfaceScalarFieldphiHbyA(adjustPhi(phiHbyAUp7//Non-orthogonalpressurecorrector{4《LTS5http://w 6下文提及的所有數(shù)學(xué)公式符號請參考《icoFoam7adjustPhi詳見《potentialFoam8非矩形修正詳見《potentialFoam//Pressurecorr==nCorr-&&nonOrth==){//}{}

if(nonOrth=={}}#includeU}}turbulence->correct();//在壓力速度求解之后，我們來修正湍流項。加入我們求解的為kEpsilon湍流模型，我們在行correct()函數(shù)的時候，會對k、epsilon方程進(jìn)行計算<<"ClockTime="<<runTime.elapsedClockTime()<<"<<nl<<}Info<<"End\n"<<return}9為何設(shè)定參考壓力詳見《potentialFoam??U=?2??=其中φ為速度勢，在OpenFOAM中用p來表Info<<"Readingfieldp\n"<<endl;volScalarFieldp(p=dimensionedScalar("zero",p.dimensions(),0.0);//詳見Info<<"ReadingfieldU\n"<<endl;volVectorFieldU(",surfaceScalarFieldphi//構(gòu)建通量，依據(jù)連續(xù)性方程，我們要確保速度的散度為0。也即通量守恒。OpenFOAM通過守恒的通量場組建速度。potentialFoam在后續(xù)的代碼中需要構(gòu)建壓力的拉斯方程并通過迭代求解(//{Uphi=fvc::interpolate(U)&}labelpRefCell0;//此處設(shè)置參考壓力網(wǎng)格單元為0，對于某些算例并沒有提供壓力固定值邊界條件，這對于求解壓力泊

(

)=其解為p=ax+b其中a和b的值我們通過p(0)=0，p(1)=1來獲得，這即為固定值邊界條件。如果我們制定兩=scalarpRefValue=0.0;//設(shè)定參考壓力為setRfCell//在OeFOM早期的版本，對于沒有提供固定值邊界條件的壓力條件，在求解器中強(qiáng)制執(zhí)行參考壓力單元為0以及參考壓力值為0（或者類似的值）。因此并不需要函數(shù)，在M目前的版本中，用戶具有更大的靈活，可以通過上述代碼已經(jīng)設(shè)置了參考壓力單元以及參考壓力值。此處函數(shù)的作用為加入用戶設(shè)定的參考壓力網(wǎng)格單元碰巧是c、y、r邊界條件，則取最近的非上述邊界1(#include"fvCFD.H"http://*************************************//intmain(intargc,char*argv[]){#include"setRootCase.H"#include"createTimeH"#include"createMesh.H"#include"readControls.H"#include"createFields.H"http://***********************************//Info<<nl<<"Calculatingpotentialflow"<<endl;//Sincesolvercontainsnotimeloopitwouldnever//functionobjectssodoitourselvesadjustPhi(phiUp);//遍歷邊界條件并進(jìn)行通量加和以確保通量守恒。對于不守恒的算例，修改壓力為zeroGradient邊界條Continuityerrorcannotberemovedbyadjusting"outflow.\nPleasecheckthevelocityboundaryconditions""and/orrunpotentialFoamtoinitialisetheoutflow."12(;;{1P//此處dimTime的單位為(0,0,1,0,0),除以壓力的單位再乘以壓力的單位后單位還是(0,0,1,0,0)dimensionSet(02,2,0,0)后單位變?yōu)?0,2,1,0,0)，進(jìn)行梯度操作后單位變?yōu)?0,0,1,0,)fvc::div(phi)//此處div(phi)的單位為(0,0,-1,0,0)，因為在div()函數(shù)執(zhí)行的時候里面除以了體積的單位，由于通,if(nonOrth==nNonOrthCorr){phipEqn.flux();//從最后收斂的壓力場組建守恒通量場，公式為?(???}}Info<<"continuityerror=<<<<=Info<<"InterpolatedUerror=<<(sqrt(sum(sqr((fvc::interpolate(U)&mesh.Sf())-<<if{}<<"ClockTime="<<runTime.elapsedClockTime()<<"<<nl<<Info<<"End\n"<<return}?2??=用SimplePISO進(jìn)行迭代求解。因此進(jìn)行非矩形修正是必要的。5《OpenFOAM用戶指南》1186 78 du d du uu 然滿足連續(xù)性方程。動量方程在1D網(wǎng)格離散化之后有： duudx d dx uu du dudx FuFuuEuPuPuWpe w aPuPaEuEaWuWpwPuaEuEaWuWpw P 先我們看沒有采用Rhie-Chow插值時候的壓力修正方程。

Rhie-ChowaPuP*aEuE*aWuW*pw*peaPuP'aEuE'aWuW'pw'pe

u'1PwPw

'

'uP

*a

pw'pe'P

uu uw'ue' au*1p'p'u*1p'p'aa aP

P pP' pW pE'uw*ue*

uu0 速度稍有不同其也為Rhie-Chow插值的精髓。將面速度表示為參類似形式：dpuu u

pE 其中ue可使用各種格式來計算比如uuuEe2dpdxau

uEu uEu uuEuPpE dpaP p dxdx p u E u P 2 2 uEuPpeepepepwpE u pEEpEpE pEpPpP p P uEuPpEEpPpEpWpE uEuPpEE3pE3pP uE pEE3pE3pPPue P

u*u*pEpP,u'u'pE'pP pE'pPaue'awPw

pE'pP',u

pP'pWu*pE'pP' a

pP'pW'aPa

2p'1p'1p'u*ua a

u*uE*uP*pEE*3pE*3pP*pW uw*uW*uP*pE*3pP*3pW*pWW

aa2pP'aa

1pE'1

pW' u*u *4p*6p*4p*

*d E 2pP'

pE'

pW'

RhieChow插1 1u*u *4p*6p*4p* *d E a P'a

pW'

pE'd非RhieChow1 1duw*ue 在非RC插值的情況下，如果存在壓力波，如圖：將不為0.消除壓力波。H.M.“Anintroductionofcomputationalfluiddynamics”LaDa,“Numericalmethodsforturbulentflow”本文檔分析simpleFoam求解器，從本文檔開始，重復(fù)的代碼段不再進(jìn)行贅述和。未注改為simpleFoam表示這個過程。simpleFoam分 #include#include"RASModel.H"#include//*************************************intmain(intargc,char{#include"setRootCase.H"#include"createTimeH"#include"createMesh.H"#include"createFields.H"#include"createFvOptions.H"http://***********************************boolFoamsimpleControl{…if{…}{}

//Settofinalise} {Info<<"Time="<<runTime.timeName()<<nl<<//---Pressure-velocitySIMPLE{#include}1<<"ClockTime="<<runTime.elapsedClockTime()<<"<<nl<<}Info<<"End\n"<<return}//Momentum,)//中的速度方程使用{volVectorFieldHbyA("HbyA",U);HbyA=rAU*UEqn().H();surfaceScalarFieldphiHbyA("phiHbyA",fvc::interpolate(HbyA)&mesh.Sf());//Non-orthogonalpressurecorrectorwhile{if{phi=phiHbyA-}}#include//Explicitlyrelaxpressureformomentum//MomentumU=HbyA-rAU*fvc::grad(p);}pisoFoam——型我們使用UEqn.A()，但是tmp類型為UEqn().A()：simpleFoamtutorialpitzDaily算例，817步收斂。和原始simpleFoam相同。非迭代求解求解器。simpleFoam為迭代求解器使用了低松弛技術(shù)。具體請查看《PISO，SIMPLE詳解》。此處從略。本文檔所述內(nèi)容大部分CFD教科書均有介紹從最基本的流體單元開始推導(dǎo)至最終的NS方程，本文檔的特點在于每一步均進(jìn)行推導(dǎo)，過連續(xù)性方歐拉觀點X方向面積流體流量:(???? 流出微元凈質(zhì)量流率Y方向面積流體流量????Z方向面積流體流量????

X方向面積流體流量流入微元凈質(zhì)量流率Y方向面積流體流量Z方向面積流體流量 (??????????????) 有

流入微元凈流率?流出微元凈流率=固定微元流體質(zhì)量增加即 )????????????

+???(????)=拉格朗日觀點

=

+

+

+

=

+??? ???? ???? ???? = +??? +???(????)=

+???(??????)=?? +??????]+?? +??(????)]=

??(??????(??????)即歐拉方法上流動微元上每單位體積上φ流體微元??的增加率 流體微元??凈流出率 流動微元??的增加對于上面的公式，由于??(????)??(??????)

可用壓力p和剪應(yīng)力來表示，剪應(yīng)力張量τ????可表示為作用在xyz9個力來表示，其中j表示方向，i表示與作用與i垂直的平面。

????)????????? ????????????+

????)????????? +

????)???????????????????????+???????????????(??+ ??(???+=

+

+

??(???+

+

+

??(???+=

+

+

+ ??(???+ ?????? +????+????+

??(???+

+

+ ??(????)+???(??????)=??(???+??????)+????????+????????+

??(???+

+???(??????) +????+????+??(????)+???(??????)=??(???+??????)+????????+????????+ +??(??????)=????+????+NS方

+

??=2(????+ )=2????+ 2 ????+

+

+

2

??????=2?????2??????????=??(????+??????)?2???????? 3+??(??????)=????+??(??(????+

+

+??(??????)=????+?? ????+ 2 ????+????

2 ???? + + + + ???(??????)=??(???????)+??(??(?????))=??(???????)=

+

+

??：+ +

+

)=

+

??????+??????+??

??

)+ ??( +)+??( + = +

??

)

?? )

??()

??()

??()

??( = + ??????(????)+??????(????)+??????(????)]+= +???(??????)+??=[

??(????)+????(????)+????(????)]=[????(????)+????(????)+????(???? + ++)=+++)=++)=+++)=

+??(?????)++??(?????)++??(?????)++??(?????)++???(??????)=????+??(?????)++???(??????)=????+??(?????)+{ +???(????)=

+

=

+

??(????????)+

=

+???(??????)=????+??(?????)+NSU,P為時均形式。因此在控制方程中便包含了脈動時時均雷諾應(yīng)大渦模直接模

????=

?????=

????(?????????即+????()+(??????′]=???+???令????=有省略源項形式（以下均省略源項）適用于各種流體的RANS+???(??????)=?????+???????? ??????????????????????)即

+???(??????)=????+??(??(???+??????))????

+

=

+

(??

+

))?

???‘????’?? ?? ???為黏性應(yīng)力，????為雷諾應(yīng)力。???????????????為壁面剪切力。由黏性應(yīng)力和雷諾應(yīng)力構(gòu)成。略黏性應(yīng)力。粘性系數(shù)法需要對????進(jìn)行?；?。????的確定是時均法的，各種湍流模型既由湍流粘性系數(shù) ?????=????????+ ?3?????????????3 ???=??(????+??????)?2?????? 其中????為湍流粘度，????為湍流運動粘度。取決于流動狀態(tài)而非物性。k為單位質(zhì)量流體的湍1

??=2 + +?? ??????+???(??????)=?????+???(??(????+??????))+???(??(????+???? 2

23??????)+=???(??令

????)+???[(??+??)(????+??????)]+ ????????=??+2有

=??

3

+???(??????)= +??? [(????+??????)]+

+???(????)=

+??? [(????+??????)] kEpsilon模NS方程乘以????

?????????????????????????? +????????

=?????

+????

+????RANS方程乘以????

???‘????’ +??

+

(??

))?

?? ?????????? ???????? ????????

?? ??′

?

=???

???‘????’??′ +???? (??′ 即

??

))+

??

??

???????

+(??′??

+????

+????

+

??

??

?? ??

???‘????’=??? +??′ (?? ?? ??))+ ?? 4

??????? 5

6即??????′ ?????????? ?????????????????′?) + ????? + ???? 1 ??673

???‘????’=??? +??′ (?? ?? ??))+ ??

4

??????? 5 1??(????′??

6

?? =2 有

????????+

2?2

????

??

?

????????????????????????????????? +?? )??? +?? )??? =?

?????? ????????????????? ′ ′ ′??=2(????+????+????

???????)

+???? =? ????

?????????′ ??+ )1

???2

?2 7

4

3

?????????????5 1???????????)

??

?? ??

=

(??′??′

??′??′??′)=

(

??????

2????

????????′ (??????+有不可壓縮流體k

????

+

????=???((??+????)????)+

(

? ?? ???? ?? ?????? ????+????

=?

((??

)????)+????1??????

+????)

?????2CFD中對于通量的概念特別重要，在各種基

SfU

SffU面單位矢量n，和面矢量同方向。nSf/SbufdSffufdSf。這個即為通過某個網(wǎng)格面的通量ufdS

dSfSdSfufu 本適用于在本地使用paraview對服務(wù)器上的結(jié)果進(jìn)行后處理，省去數(shù)據(jù)的麻 顯示 etingtmp<type><type>usingnamespaceFoam;int{scalarField ,1.0), scalarFieldf3=f1+//f1(new //f2(new //tmp<scalarField>f3=f1+f2;Info<<f3<<endl;cout<<"RunTime="<<finish-start<<"/"<<CLOCKS_PER_SEC<<"s"<<std::endl;return} 內(nèi)存使用如下（2302musingnamespaceFoam;int{clocktstart,finish;//scalarField ,1.0), //scalarFieldf3=f1+f2;f1(new f2(new tmp<scalarField>f3=f1+Info<<f3<<endl;cout<<"RunTime="<<finish-start<<"/"<<CLOCKS_PER_SEC<<"s"<<return} 內(nèi)存使用如下（776mscalarField ,1.0), scalarFieldf3=f1+會調(diào)用重載的操作符“+”，在返回的時候，他會創(chuàng)建一個類型然后并返回。這樣對是很安全和方便的。但如果對于CFD計算涉及的大數(shù)據(jù)處理，這會增加內(nèi)存占用以及的時f1(new f2(new tmp<scalarField>f3=f1+f1f2f3tmptmpf3f1f2f1f2的和賦f3，f1,f2被釋放。程序目前只占用一個f3的內(nèi)存，約77m。tmp會減少內(nèi)存的使用量以及加快運行速度。下面舉例說明一個如何在求解器中使用tmp類，參見simpleFoam，我們有這一段代碼：fvm::div(phi,+turbulence-volVectorFieldHbyA("HbyA",U);HbyA=rAU*UEqn().H();我們構(gòu)造了一個用于求解動量方程的矢量矩陣系數(shù)的tmp+fvm::div(phi,+turbulence-其并沒有使用tmpPISO算法中，在每個內(nèi)循環(huán)中，壓力方程進(jìn)行了多次修HbyA項后，UEqnpeak把icoFoam求解器的思想搞到精通以方便對其他求解器的算法理解。

UUUURpg U UUUUUU

UUU

URpg

RRUUURRpg 注意到,R均為對稱張量，以及對上述方程左右除以

U

URRpg UTU2 3 依據(jù)Boussinesqeddyviscosity

RUTU2IU2 3 k 1此處的壓力為1p12

R UTU IU2eff eff eff eff,U

UU

pg eff eff