# 加载中...

WayneWangY
• 博客等级：
• 博客积分：0
• 博客访问：91,777
• 关注人气：155
• 获赠金笔：0支
• 赠出金笔：0支
• 荣誉徽章：

## OpenFOAM求解器（二）：通量与求解器算法

(2015-10-30 17:14:58)

### peqn.flux

1.

phi = phiHbyA - pEqn.flux();

fvc::div(phiHbyA)的问题是由于还要做一个积分，可参考

template
tmp >
div
(
const GeometricField& ssf
)
{
return tmp >
(
new GeometricField
(
"div("+ssf.name()+')',
fvc::surfaceIntegrate(ssf)
)
);
}

2.
For cases which do no have a pressure boundary adjust the balance of
fluxes to obey continuity.  Return true if the domain is closed.

(
surfaceScalarField& phi,
const volVectorField& U,
volScalarField& p
);

```
```
```        if
```
```        (
```
```            magAdjustableMassOut > VSMALL
```
```         && magAdjustableMassOut/totalFlux > SMALL
```
```        )
```
```        {
```
```            massCorr = (massIn - fixedMassOut)/adjustableMassOut;
```
```        }
```
```        else if (mag(fixedMassOut - massIn)/totalFlux > 1e-8)
```
```        {
```
```            FatalErrorIn （这是这个函数唯一有报错的地方）
```
```
```
```   通俗解释上面的代码：啥时候报错呢？就是可调节流量很小，而流量的不平衡较大时。
```
```如果没有报错的话，就按照比例系数massCorr对所有的可调节流量进行修正。
```
```参考 http://www.cfd-china.com/topic/305/adjustphi的作用是检查边界条件
```
```
```
```http://www.cfd-online.com/Forums/openfoam-programming-development/83876-pref-adjustphi.html
```

3.correctBoundaryConditions();
```template<<span style=" color:#808000;">class Type, template<<span style=" color:#808000;">class> class PatchField, class GeoMesh>
```
```void Foam::GeometricField<<span style=" color:#800080;">Type, PatchField, GeoMesh>::
```
```correctBoundaryConditions()
```
```{
```
```    this->setUpToDate();
```
```    storeOldTimes();
```
```    boundaryField_.evaluate();
```
```}
```
```  这里涉及evaluate()函数在每个边界条件中的重载和边界条件的调用，参考http://www.cfd-online.com/Forums/openfoam-programming-development/129271-how-boundary-conditions-called-openfoam-solvers.html
```

4.ddtCorr(U, phi)

```template<<span style=" color:#808000;">class Type>
```
```tmptypename flux<<span style=" color:#800080;">Type>::type, fvsPatchField, surfaceMesh> >
```
```ddtCorr
```
```(
```
```    const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
```
```    const GeometricField
```
```    <</pre>
```
```        typename flux<<span style=" color:#800080;">Type>::type,
```
```        fvsPatchField,
```
```        surfaceMesh
```
```    >& phi
```
```)
```
```{
```
```    return fv::ddtScheme<<span style=" color:#800080;">Type>::New
```
```    (
```
```        U.mesh(),
```
```        U.mesh().ddtScheme("ddt(" + U.name() + ')')
```
```    )().fvcDdtPhiCorr(U, phi); 这个函数位于ddtScheme.H
```
```}
```
```
```
```  virtual tmp<<span style="color: rgb(128, 0, 128);">fluxFieldType> fvcDdtPhiCorr
```
```        (
```
```            const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
```
```            const fluxFieldType& phi
```
```        ) = 0;
```
```
```
```纯虚函数。其中一个实现如下：
```

```template<<span style=" color:#808000;">class Type>
```
```tmp<<span style=" color:#808000;">typename EulerDdtScheme<<span style=" color:#800080;">Type>::fluxFieldType>
```
```EulerDdtScheme<<span style=" color:#800080;">Type>::fvcDdtPhiCorr
```
```(
```
```    const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
```
```    const fluxFieldType& phi
```
```)
```
```{
```
```    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();
```
```
```
```    fluxFieldType phiCorr
```
```    (
```
```        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
```
```    );
```
```
```
```    return tmp<<span style=" color:#800080;">fluxFieldType>
```
```    (
```
```        new fluxFieldType
```
```        (
```
```            IOobject
```
```            (
```
```                "ddtCorr(" + U.name() + ',' + phi.name() + ')',
```
```                mesh().time().timeName(),
```
```                mesh()
```
```            ),
```
```            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
```
```           *rDeltaT*phiCorr
```
```        )
```
```    );
```
```}
```

5.pEqn.flux()的代码实现

http://www.cfd-online.com/Forums/openfoam/114933-description-flux-method.html

===========================================================

phiP = pEqn.flux();

phi = phiHbyA - pEqn.flux(); //如icoFoam, simpleFoam, pimpleFoam

phi == pEqn.flux(); //如reactingFoam, rhoReactingFoam
phi = phiHbyA + pEqn.flux();

phic = phiHbyA - pEqn.flux()/alphacf; //如DPMFoam
+rAUc*fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf);

phi += (phiGradp + pEqn.flux())/rhof; //如cavitatingFoam

phi = pEqn.flux(); //如sonicFoam

phiP = pEqn.flux(); //如porousMultiphaseFoam

（*）

（3.2）

（3.17）

（3.136）

（3.139）

（3.140）

（3.144）

（3.137）

（3.138）

（3.141）

（3.142；3.143）

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)  //层流
);
simpleFoam中的动量预测UEqn.H与上面的相同，使用turbulence->divDevReff(U)考虑了湍流，使用fvOptions考虑了源项，并引入了方程的松弛。（这块比较熟了，有任何问题可留言讨论）

// --- PISO loop
然后是压力求解，
volScalarField rAU(1.0/UEqn.A()); //A() 返回fvVectorMatrix的中心系数，这里即a_P。This is the famous "operator splitting".

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H(); //H() 返回fvVectorMatrix的H源项，这里即H(U)。HbyA见[7]中（3.1）
surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //simpleFoam中这里是 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); ddtCorr最后分析。
);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA) //用不满足无散度的通量求解压力方程3.143，[7]中的（3.9）
);

pEqn.setReference(pRefCell, pRefValue); //压力参考值
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi = phiHbyA - pEqn.flux();  //得到无散度的通量[6]，最后分析
}
}

U = HbyA - rAU*fvc::grad(p); //3.139，此速度满足连续性方程
U.correctBoundaryConditions();

pisoFoam.c与上述过程基本一致。

"phiHbyA",(fvc::interpolate(HbyA) & mesh.Sf()) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)phi = phiHbyA - pEqn.flux();

fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)的作用与PISO的非正交修正有关，大意可能是利用前一次的速度场插值通量场和通量场的差别来对本次的通量场进行修正，对于结构化，完全正交的网格可能不起作用，可参考[9-10]，代码在finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C，其中函数及之间的调用关系见[12]
pEqn.flux()是由fvScalarMatrix返回通量phi = phiHbyA - pEqn.flux(); 是为了得到无散度的通量[5-6]。具体来说，flux()函数作用于fvScalarMatrix，返回的应该是(fvc::interpolate(rUA*fvc::grad(p)) & mesh.Sf（见3.144），这与[1]中的用法是相符合的，因为在porousMultiphaseFoam中速度本身没有输运（达西定律），也就没有H(U)相关的项（见3.137。在实际实现中，flux()返回的值可能还包含非正交修正相关的项[5]。flux()具体的返回值与方程中的fvm项有关，更具体的讨论见[13][14]，可以多看几个求解器来build confidence。

reference:
1. porousMultiphaseFoam详解（一）& OpenFOAM求解器开发入门（一）：求解策略
2. 东岳流体技术文档：何为通量
3. OpenFOAM guide/The SIMPLE algorithm in OpenFOAM
4. OpenFOAM guide/The PISO algorithm in OpenFOAM
5. http://www.cfd-online.com/Forums/openfoam/114933-description-flux-method.html
6. http://www.cfd-online.com/Forums/openfoam-programming-development/83638-phi-peqn-flux-vs-linearinterpolate-u-mesh-sf.html（关于PISO的讨论）
7. 东岳流体技术文档：icoFoam 详解
8. 东岳流体技术文档：potentialFoam 详解
9. http://www.cfd-online.com/Forums/openfoam-solving/60096-ddtphicorr.html
10. http://www.cfd-online.com/Forums/openfoam/84210-fvc-ddtphicorr-rua-u-phi.html
11. http://www.cfd-online.com/Forums/openfoam-programming-development/89144-wrong-operator-transonic-pressure-equation.html
12. http://www.cfd-online.com/Forums/openfoam-programming-development/70444-programming-questions.html
13. http://www.cfd-online.com/Forums/openfoam-solving/109440-peqn-flux.html
14. http://www.cfd-online.com/Forums/openfoam-programming-development/68027-different-implementation-solving-peqn.html
15. http://www.cfd-online.com/Forums/openfoam-solving/90435-interpolation-pressure.html
16. http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture11.pdf

0

• 评论加载中，请稍候...

发评论

以上网友发言只代表其个人观点，不代表新浪网的观点或立场。

新浪BLOG意见反馈留言板　电话：4000520066 提示音后按1键（按当地市话标准计费）　欢迎批评指正

新浪公司 版权所有