openfoam-解析流体技术二icoFoam详解

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"//通过fvOption来修正源项，主要用于在运行时添加多孔介质、MRF多重参考系、显性源项、热"//*************************************intmain(intargc,char{ #include"createMesh.H"//创建网格对象，必备头文件。请忽略#include"createFields.H"//从求解器根 #include"createFvOptions.H"//创建源项，无需源项可删除。请忽略//***********************************InfonCalculatingscalartransport\nendl;在OpenFOAM中采用Info来在终端输出信息，内为需要输出的信息#include"CourantNo.H"//计算库郎数，详见《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"//*************************************intmain(intargc,char{"#include"createTimeH创建时间#include"createMesh.H创建网格#include"createFields.H"//创建场//***********************************//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"""//#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""//"//*************************************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"//*************************************//intmain(intargc,char*argv[]){#include"setRootCase.H"#include"createTimeH"#include"createMesh.H"#include"readControls.H"#include"createFields.H"//***********************************//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"//***********************************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