重新看icoFoam

[史浩]

今天重新看OpenFOAM中的icoFoam求解器，有个疑问：在压力修正过程中，感觉每次求解的压力方程中的系数都没有变化。压力修正过程代码如下，压力方程中的rAU和HbyA都是从动量预测方程（UEqn）求出的，UEqn并没有随着U变化而变化，所以压力方程按理并没有发生变化，但是PISO算法中，压力方程中的系数是要随着速度的更新而更新的。那修改压力方程中系数的部分在什么地方？

         while (piso.correct())
         {
             volScalarField rAU(1.0/UEqn.A());
             volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
                 fvc::flux(HbyA)
               + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
             );
 
             adjustPhi(phiHbyA, U, p);
 
             // Update the pressure BCs to ensure flux consistency
             constrainPressure(p, U, phiHbyA, rAU);
 
             // Non-orthogonal pressure corrector loop
             while (piso.correctNonOrthogonal())
             {
                 // Pressure corrector
 
                 fvScalarMatrix pEqn
                 (
                     fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
 
                 pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
                 if (piso.finalNonOrthogonalIter())
                 {
                     phi = phiHbyA - pEqn.flux();
                 }
             }
 
             #include "continuityErrs.H"
 
             U = HbyA - rAU*fvc::grad(p);
             U.correctBoundaryConditions();
         }

[队长别开枪]

@史浩 UEqn.H()是关于U的一个函数，因此每次迭代更新U以后也会更新phiHbyA。

[李东岳]

没有修正压力系数。这是可以理解的因为库郎数的作用瞬态计算时间布长内的变化并没有稳态计算那么大。可以这样想：

1. 如果修正压力系数
2. 如果每次求解一次压力、更新一次速度

你就变成了SIMPLE

[winsway_zero]

volScalarField rAU(1.0/UEqn.A());

按照道理说这个应该是矢量场才对，因为三个方向volVectorField。

[李东岳]

@winsway_zero 是三个方向，但是他们的对角线元素是一样的，

[winsway_zero]

@李东岳 如果只是离散内部面是相同的，但是将边界的影响引入到了对角线上的Ap系数后，对角线系数就不一样了。比如Wall边界条件，里面涉及到了应力的分解：

[winsway_zero]

@李东岳 我下面仔细推导了一下：
首先 rAU(1.0/UEqn.A());，这个公式的计算得到的结果是：

对角系数

这是A()函数：

 template<class Type>
 Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
 {
     tmp<volScalarField> tAphi
     (
         volScalarField::New
         (
             "A("+psi_.name()+')',
             psi_.mesh(),
             dimensions_/psi_.dimensions()/dimVol,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
 
     tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
     tAphi.ref().correctBoundaryConditions();
 
     return tAphi;
 }
 

这是D()对角系数的平均化处理

 template<class Type>
 Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
 {
     tmp<scalarField> tdiag(new scalarField(diag()));
     addCmptAvBoundaryDiag(tdiag.ref());
     return tdiag;
 }

边界对对角系数的影响：

 template<class Type>
 void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
 {
     forAll(internalCoeffs_, patchi)
     {
         addToInternalField
         (
             lduAddr().patchAddr(patchi),
             cmptAv(internalCoeffs_[patchi]),
             diag
         );
     }
 }

从上面的代码可以得到：
$$A_p=\frac{\overline{D}}{\mathrm{\Delta }V}$$
其中：
$$\overline{D}=average(a_p^{\prime })+diag$$
式子中average(Ap')表示的是边界对主对角线系数的影响的平均；diag表示的是内部面离散的主对角线系数。

周围系数作为源项：

H()

 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
 Foam::fvMatrix<Type>::H() const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi
     (
         GeometricField<Type, fvPatchField, volMesh>::New
         (
             "H("+psi_.name()+')',
             psi_.mesh(),
             dimensions_/dimVol,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
     GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref();
 
     // Loop over field components
     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
     {
         scalarField psiCmpt(psi_.primitiveField().component(cmpt));
 
         scalarField boundaryDiagCmpt(psi_.size(), 0.0);
         addBoundaryDiag(boundaryDiagCmpt, cmpt);
         boundaryDiagCmpt.negate();
         addCmptAvBoundaryDiag(boundaryDiagCmpt);
 
         Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
     }
 
     Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
     addBoundarySource(Hphi.primitiveFieldRef());
 
     Hphi.primitiveFieldRef() /= psi_.mesh().V();
     Hphi.correctBoundaryConditions();
 
     typename Type::labelType validComponents
     (
         psi_.mesh().template validComponents<Type>()
     );
 
     for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
     {
         if (validComponents[cmpt] == -1)
         {
             Hphi.replace
             (
                 cmpt,
                 dimensionedScalar(Hphi.dimensions(), 0)
             );
         }
     }
 
     return tHphi;
 }

 template<class Type>
 void Foam::fvMatrix<Type>::addBoundaryDiag
 (
     scalarField& diag,
     const direction solveCmpt
 ) const
 {
     forAll(internalCoeffs_, patchi)
     {
         addToInternalField
         (
             lduAddr().patchAddr(patchi),
             internalCoeffs_[patchi].component(solveCmpt),
             diag
         );
     }
 }

这里面的求解过程包含了：
$$H=\frac{[-\sum a_{\mathbf{N}}{\mathbf{U}}_{\mathbf{N}}+(\mathbf{b}+{\mathbf{b}}^{\prime })+(average(a_p^{\prime })-a_p^{\prime }){\mathbf{U}}_{\mathbf{C}}]}{\mathrm{\Delta }V}$$

最后：

$$HbyA=\frac{H}{A}=\frac{\mathrm{\Delta }V}{average(a_p^{\prime })+diag}\frac{[-\sum a_{\mathbf{N}}{\mathbf{U}}_{\mathbf{N}}+(\mathbf{b}+{\mathbf{b}}^{\prime })+(average(a_p^{\prime })-a_p^{\prime }){\mathbf{U}}_{\mathbf{C}}]}{\mathrm{\Delta }V}$$
$$HbyA=\frac{H}{A}=\frac{[-\sum a_{\mathbf{N}}{\mathbf{U}}_{\mathbf{N}}+(\mathbf{b}+{\mathbf{b}}^{\prime })+(average(a_p^{\prime })-a_p^{\prime }){\mathbf{U}}_{\mathbf{C}}]}{average(a_p^{\prime })+diag}$$