icoFoamicoFoam采用的是標準瞬態PISO算法，本文檔從最基本的NS方程離散開始推導，并和OpenFOAM植入的代碼相對應以供理解瞬態PISO算法如何在OpenFOAM中實現以及NS方程在OpenFOAM中的離散。對于穩態的標準SIMPLE算法以及PIMPLE算法的對比請參考《簡單易懂的PISO，SIMPLE算法對比》NS方程離我們有動量方程

u

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

urSur

unu

方程1.3代表一個向量，速度的下標n表示臨點的速度，上標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其中由于為面預測速度，我們需要進行插值，在此步我們可以引入各種格式。假設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 UEqn.A()即????保持不[W用2]:solve(UEqnUEqn.H()中的系數[W用2]:solve(UEqn至此，如果我們加入壓力梯度的顯性項。即其為速度預測方程。求解此方程后（solve…），我們有預測速度??????第一次修正步（由????計算????再計算

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

purrHbyArp

我們令（依據連續性方程

uprrurrS p依據公式3.3，我們對面值進行插值AurrHbyAr1 Apf p,f

r 整理有

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

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

u

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

ff

udVdt udSdtuSt=

S 其中

urSur

unu

方程1.3代表一個向量，速度的下標n表示臨點的速度，上標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其中由于為面預測速度，我們需要進行插值，在此步我們可以引入各種格式。假設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()中的系數即????保持不至此，如果我們加入壓力梯度的顯性項。即[W[W用2]:solve(UEqn其為速度預測方程。求解此方程后（solve…），我們有預測速度??????第一次修正步（由????計算????再計算首先我們通過預測的????組建HbyAr1Aur+Er

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

purrHbyArp

我們令（依據連續性方程

uprrurrS p依據公式3.3，我們對面值進行插值urrHbyAr1 Apf p,fA

r 整理有

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

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

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

∑|?????????????|2???????????????alphainterFoam中，我們采用歐拉法進行時間步離散。但是在LTSinterFoam中，我們對時間項采用局部時間步離散；正后的alpha場；1fvc::smooth,fvc::spread,fvc::sweepLTS中使用其他場參數對局部時間步場進行光順，參LTS代碼請參考之前的文檔。本文檔涉及一些C++類相關內容，有關C++類的內容推薦《C++名字為fluid，他需要傳入參數U和phi。略去構造機理。singlePhaseTransportModel.C1代碼，有下述成員函數：=0，表明這是虛基類（抽象基類viscosityModel.H相同的文件下，1 2 至此singlePhaseTransportModel類，名在執行fluid.correct()的時候進入：singlePhaseTransportModel.correct()，viscosityModel.correct()，由于viscosityModel為虛基類，進入具體類：nu_= returnmax(…,…)本語句返回計算粘度值注意到在牛頓流體中我們的速度方程其中拉斯項為laplacian(nu,U)，在肥牛臀流體中我場。方法和fluid.correct()一致。從略。其他內容和icoFoam一致。分析，一步一步推導VOF模型。VOF定義一個示蹤函數11，流20，零到一之間的值則表示自由表面。可見，VOF的界面并不是直接算出來的，而是簡介的計算每個網格內的相場值而后處和界面壓力降?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中，我們用表示體積分數，有α=

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 續性方程中的對流項中的所有項中，因為我們要保證相分數的有界和相界面的，在過去而這種方法雖然保持了守恒性但是卻會是的。另法是：添加一個對流項以壓縮界u1u interFoamtwoPhaseEulerFoam中植入上述方程進行求解。經測試mixerVesselAMI算例，和使用MULES求解的結果一致。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文件夾中，對計算機語言不我們再看OpenFOAM2.3.0中的一段描述：fluidsinglePhaseTransportModel類，我們先不管他fluid在這里被調用。我們不考慮他是如何實現的。可以這樣來想：在非利用當前的速度來求解粘度。上文的fluid.correct();就是實現的這個功能。while{Info<<"Time="<<runTime.timeName()<<nl<<#include"CourantNo.H"+fvm::div(phi,(fvc::grad(U)&定義一個粘度場，powerLaw模型需要的k寫，會大大增加了代碼。這不符合OpenFOAM的代碼規范Pimplepiso算法進行時間步進，時間步內若有需要對速度矩陣系數、非最終迭代步的壓力場進行低松弛求解提高穩定性。pimpleFoam中的大部分peapooam邊界條件。雷同的代碼不再注釋，請參考《simpleFoam解析》以及《pisoFoam解析》另外提醒的是，幾個基本的OpenFOAM求解器（laplacianFoam,potentialFoam,scalarTransportFoam,icoFoam,pisoFoam,simpleFoam,pimpleFoam）已經剖析完成。其余的OpenFOAMpimple、simpleLTS局部時間步框架下的求解器。CFD技術居多。其余的求解器不再代碼進行通篇分析，只分析#include#include"turbulenceModel.H"""http://#include"fvIOoptionList.H"#include"IOMRFZoneList.H"#include在在pimpleFoam調用setSnGrad模板函數，計算表達式表達式計算之后，調動setSnGrad內的函數updateCoeffs()更新邊界條件通過fixedFluxPressureFvPatchScalarField.C函數內的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通過庫郎數設定Info<<"Time="<<runTime.timeName()<<nl<<{#include{#include}if{}}<<"ClockTime="<<runTime.elapsedClockTime()<<"<<nl<<}return0;}//SolvetheMomentumif(pimple.momentumPredictor())//此處如果開啟動量預測，求解下述方程，默認為關閉。具體情況需要具體分析，Henry{}1《LTSHbyA=rAU*UEqn().H();if(pimple.nCorrPISO()<={}adjustPhi(phiHbyA,U,p);//UpdatethefixedFluxPressureBCstoensureflux);//此處即為壓力中通過setSnGrad設定壓力邊界//Non-orthogonalpressurecorrectorwhile{//PressurecorrectorfvScalarMatrixpEqnpEqn.setReference(pRefCell,pRefValue);if(pimple.finalNonOrthogonalIter()){phi=phiHbyA-}}#include//Explicitlyrelaxpressureformomentum();//U=HbyA-rAU*fvc::grad(p);pimpleFoampisoFoamsimpleFoam對于不可壓縮層流我們有NS味著需要用已知時間步的速度來計算Ap和An。這進而又分為兩個情況：針對以上兩個問題的穩態瞬態的不同，于是有了非迭代的瞬態 SIMPLEPISO也可以用于穩態并且SIMPLE也可以用于瞬態。其區別就在于舉例：如果我們考慮SIMPLE和PISO的話，會有很大的不同。其在于瞬態的PISO是iterative量的時間和資源。瞬態SIMPLE也具有以下特點：SIMPLE本身側重于對非線性化的處理。在大時間步下，速度的線性化滯后變的時候，瞬態SIMPLE算在每個時間步內進行大量迭代來計算。PIMPLE可以使用松弛技術。值得注意的是在源代碼中pisoFoam依然存在(),（SIMPLE相同PISO的多次壓力修正以對速度壓力耦合問題壓力方程要實現的關鍵問題就是從連續性方程衍生出壓力方程以使得每個網格質量守恒PISOSIMPLE1PISO的區別參見《東岳流體技術文檔（二）——icoFoam2《東岳流體技術文檔（二）——icoFoam詳解》內有從連續性方程推導壓力泊松方程的詳細步驟，此處是動量預測：依據p0預測速度進入當前時間步，我們有壓力p1,注意每個時間步內的Ueqn.H()是變化組建速度方程的Ueqn.A和否否收收斂候引入了滯后，因此判定收得HbyA后清理速度方程便于內存管依據p0計算速度u1，并對速度方程低松進入當前迭代步，我們有是否否本文檔剖析pisoFoam的源代碼，大量了《icoFoam詳解》中NS方程離散的相關內容。需要注意的是icoFoampisoFoamPISONS方程求解，icoFoam只能用來求解層流。pisoFoam可以求解湍流以及層流。更加具有普適性的求解器為概念。因此可以用《icoFoam詳解》作為補充閱讀。本文檔初級代碼分析將省略。請參見《laplacianFoam,scalarTransportFoam流模型。OpenFOAM-2.3.0版本之前尚不支持fvOptions源項。首先我們分析createFields.H：Info<<"Readingfieldp\n"<<endl;volScalarFieldp(Info<<"ReadingfieldU\n"<<endl;volVectorFieldU( labelpRefCell=0;1,);//需要構造的類名 構造的類的具體名稱傳入的參此行代碼的意思即為將UphisinglePhaseTransportModel類，其名稱為laminarTransport1請參考《potentialFoam 一個成指針來使用，有關如何使用auto_ptr建議閱讀C++專業書籍。我們在此進行簡單介紹。本行代碼規范autoPtr智能指

類名類名 類之中的New函個名為turbulence的指針。對于初學者理解此意義即可。下面我們進入pisoFoam主函數#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并沒有外循環步，此步用來和PIMPLE算法相統一，對于沒有用到的，我們會在編譯中遇到提示：=//constintnNonOrthCorrpisoDict.lookupOrDefault<int>("nNonOrthogonalCorrectors0);//PISO非矩形修正步3，默認為0//constboolmomentumPredictorpisoDict.lookupOrDefault("momentumPredictor"true);//動量預測開=",23《potentialFoam{//MomentumUEqnif{}kEpsilon湍流模型，我們在kEpsilon會找到下面的函數：{}

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

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

(

)=其解為p=ax+b其中a和b的值我們通過p(0)=0，p(1)=1來獲得，這即為固定值邊界條件。如果我們制定兩=scalarpRefValue=0.0;//設定參考壓力為setRfCell//在OeFOM早期的版本，對于沒有提供固定值邊界條件的壓力條件，在求解器中強制執行參考壓力單元為0以及參考壓力值為0（或者類似的值）。因此并不需要函數，在M目前的版本中，用戶具有更大的靈活，可以通過上述代碼已經設置了參考壓力單元以及參考壓力值。此處函數的作用為加入用戶設定的參考壓力網格單元碰巧是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);//遍歷邊界條件并進行通量加和以確保通量守恒。對于不守恒的算例，修改壓力為zeroGradient邊界條Continuityerrorcannotberemovedbyadjusting"outflow.\nPleasecheckthevelocityboundaryconditions""and/orrunpotentialFoamtoinitialisetheoutflow."12(;;{1P//此處dimTime的單位為(0,0,1,0,0),除以壓力的單位再乘以壓力的單位后單位還是(0,0,1,0,0)dimensionSet(02,2,0,0)后單位變為(0,2,1,0,0)，進行梯度操作后單位變為(0,0,1,0,)fvc::div(phi)//此處div(phi)的單位為(0,0,-1,0,0)，因為在div()函數執行的時候里面除以了體積的單位，由于通,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進行迭代求解。因此進行非矩形修正是必要的。5《OpenFOAM用戶指南》1186 78 du d du uu 然滿足連續性方程。動量方程在1D網格離散化之后有： 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求解器，從本文檔開始，重復的代碼段不再進行贅述和。未注改為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為迭代求解器使用了低松弛技術。具體請查看《PISO，SIMPLE詳解》。此處從略。本文檔所述內容大部分CFD教科書均有介紹從最基本的流體單元開始推導至最終的NS方程，本文檔的特點在于每一步均進行推導，過連續性方歐拉觀點X方向面積流體流量:(???? 流出微元凈質量流率Y方向面積流體流量????Z方向面積流體流量????

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

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

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

=

+

+

+

=

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

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

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

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

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

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

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

+

+

??(???+

+

+

??(???+=

+

+

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

??(???+

+

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

??(???+

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

+

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

+

+

2

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

+

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

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

+

+

??：+ +

+

)=

+

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

??

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

??

)

?? )

??()

??()

??()

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

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

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

+

=

+

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

=

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

????=

?????=

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

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

+

=

+

(??

+

))?

???‘????’?? ?? ???為黏性應力，????為雷諾應力。???????????????為壁面剪切力。由黏性應力和雷諾應力構成。略黏性應力。粘性系數法需要對????進行模化。????的確定是時均法的，各種湍流模型既由湍流粘性系數 ?????=????????+ ?3?????????????3 ???=??(????+??????)?2?????? 其中????為湍流粘度，????為湍流運動粘度。取決于流動狀態而非物性。k為單位質量流體的湍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。這個即為通過某個網格面的通量ufdS

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

UUUURpg U UUUUUU

UUU

URpg

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

U

URRpg UTU2 3 依據Boussinesqeddyviscosity

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

R UTU IU2eff eff eff eff,U

UU

pg eff eff