OpenFOAM

fixedFluxPressure 边界条件

本篇说说我对 fixedFluxPressure 边界条件的理解。

先来看 OpenFOAM-2.2.x 里的

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
class fixedFluxPressureFvPatchScalarField
:
    public fixedGradientFvPatchScalarField
{
    // Private data

        //- Name of the predicted flux transporting the field
        word phiHbyAName_;
        //- Name of the flux transporting the field
        word phiName_;
        //- Name of the density field used to normalise the mass flux
        //  if neccessary
        word rhoName_;
        //- Name of the pressure diffusivity field
        word DpName_;
        //- Is the pressure adjoint, i.e. has the opposite sign
        Switch adjoint_;

public:

    //- Runtime type information
    TypeName("fixedFluxPressure");
    ......
    ......
    //- Update the coefficients associated with the patch field
    virtual void updateCoeffs();
   
    ......
};

其中 updateCoeffs 函数定义如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
    const surfaceScalarField& phiHbyA =
        db().lookupObject<surfaceScalarField>(phiHbyAName_);
    const surfaceScalarField& phi =
        db().lookupObject<surfaceScalarField>(phiName_);
    fvsPatchField<scalar> phiHbyAp =
        patch().patchField<surfaceScalarField, scalar>(phiHbyA);
    fvsPatchField<scalar> phip =
        patch().patchField<surfaceScalarField, scalar>(phi);
    const scalarField *DppPtr = NULL;
    if (db().foundObject<volScalarField>(DpName_))
    {
        DppPtr =
            &patch().lookupPatchField<volScalarField, scalar>(DpName_);
    }
    else if (db().foundObject<surfaceScalarField>(DpName_))
    {
        const surfaceScalarField& Dp =
            db().lookupObject<surfaceScalarField>(DpName_);
        DppPtr =
            &patch().patchField<surfaceScalarField, scalar>(Dp);
    }
    if (adjoint_)
    {
        gradient() = (phip - phiHbyAp)/patch().magSf()/(*DppPtr);
    }
    else
    {
        gradient() = (phiHbyAp - phip)/patch().magSf()/(*DppPtr);
    }
    fixedGradientFvPatchScalarField::updateCoeffs();
}

代码的说明如下：

这个边界条件，继承自 fixedGradientFvPatchScalarField，可见它是一个用于压力场的第二类边界条件。其用如下公式来计算边界上的压力梯度：
$$
\nabla(p) = \frac{\phi_{H/A} - \phi}{|S_f| D_p}
$$

在设置边界条件的时候，可以指定 phiHbyA， phi 和 Dp 对应的场，如果不指定，则三个量对应的默认场名分别是 phiHbyA， phi 和 Dp，这时如果你的程序中找不到这三个场，那就将出错了。

这个边界条件可以这样理解：
首先来看单相流的情形，比如 icoFoam，半离散化的动量方程为
$$
AU=H-\nabla p
$$
即
$$
\nabla p = H - AU
$$
在入口处，理论上当计算收敛后，压力梯度应当是0。但是在构建压力方程

1
2
3
4
fvScalarMatrix pEqn
(
   fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

的时候，速度 $U$ 还没有经过修正，即这里的 $U$ 是预测值，所以，迭代过程中入口压力梯度不一定为零。将入口压力设置为等于 $H/A-U$，有助于提高计算稳定性。
从代码角度看，对压力方程而言，需要设置的是 inlet 边界上的 $(\nabla p)_{\bot} \cdot rAU$，
$$
(\nabla p)_{\bot} \cdot rAU = (\nabla p)_{\bot}/A_f = \frac{1}{A_f}(\nabla p \cdot \vec{n} \cdot |S_f|) = (H/A) \cdot S_f-U\cdot S_f
$$
设置的边界条件应为
$$
\nabla p \cdot \vec{n} = \frac{(H/A) \cdot S_f-U\cdot S_f}{\big |S_f \big| \cdot rAU}
$$

但实际上，fixedFluxPressureFvPatchScalarField 主要是用于两相流中，以 twoPhaseEulerFoam 为例，
半离散的动量方程为：
$$
U_{a}=\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g
$$

$$
U_{b}=\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g
$$

对 $U_a$ 项两边同时乘以 $\alpha_a$ ，$U_b$ 项两边同时乘以 $\alpha_b$ ，然后合并起来，得到
$$
\begin{align*}
\Big(\frac{\alpha_a}{a_{p,a}\rho_a}+ \frac{\alpha_b}{a_{p,b}\rho_b}\Big)\nabla p & = \alpha_a \Big[\frac{1}{a_{p,a}}H(U_a)+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g\Big] + \alpha_b \Big[ \frac{1}{a_{p,b}}H(U_b)+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g\Big] \\
& - (\alpha_a U_a + \alpha_b U_b)
\end{align*}
$$

转换成界面通量形式，则
$$
D_p \nabla p \cdot S_f = \alpha_a \cdot phiHbyA1 + \alpha_b \cdot phiHbyA2 - (\alpha_a \cdot phi1 + \alpha_b \cdot phi2) = phiHbyA - phi
$$

于是得到压力的边界条件为
$$
\nabla p \cdot \vec{n} = \frac{phiHbyA - phi}{Dp \cdot \big|S_f\big|}
$$
这大概就是 fixedFluxPressureFvPatchScalarField 的物理含义吧。

在 OpenFOAM-2.3.x 以后，fixedFluxPressureFvPatchScalarField 做了修改，不需要指定 Dp 等这些量的值，而是在求解器中直接指定压力梯度。
比如 twoPhaseEulerFoam 中

1
2
3
4
5
6
7
8
9
10
11
12
13
14
setSnGrad<fixedFluxPressureFvPatchScalarField>
   (
       p.boundaryField(),
       (
           phiHbyA.boundaryField()
         - mrfZones.relative
           (
               alpha1f.boundaryField()
              *(mesh.Sf().boundaryField() & U1.boundaryField())
             + alpha2f.boundaryField()
              *(mesh.Sf().boundaryField() & U2.boundaryField())
           )
       )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
   );

setSnGrad 这个函数的定义在 fixedFluxPressureFvPatchScalarField 中定义：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

template<class GradBC>
    inline void setSnGrad
    (
        volScalarField::GeometricBoundaryField& bf,
        const FieldField<fvsPatchField, scalar>& snGrad
    )
    {
        forAll(bf, patchi)
        {
            if (isA<GradBC>(bf[patchi]))
            {
                refCast<GradBC>(bf[patchi]).updateCoeffs(snGrad[patchi]);// 调用带一个参数的 updateCoeffs 函数
            }
        }
    }

// 实际调用的是这个函数，这个函数调用完以后，curTimeIndex_ 将等于当前时间步的标签。
void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs
(
    const scalarField& snGradp
)
{
    if (updated())
    {
        return;
    }

    curTimeIndex_ = this->db().time().timeIndex();

    gradient() = snGradp;
    fixedGradientFvPatchScalarField::updateCoeffs();
}

// 但是，很多其他地方仍然需要不带参数的 updateCoeffs 函数接口。
// 这个函数将不起实质作用，但是，如果在调用这个无参的 updateCoeffs 函数之前，没有先调用带一个参数的 updateCoeffs 函数，
// curTimeIndex_ 也没有更新到当前的时间步标签，那就报错，表明不适合使用这个边界条件。
void Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    if (curTimeIndex_ != this->db().time().timeIndex())
    {
        FatalErrorIn("fixedFluxPressureFvPatchScalarField::updateCoeffs()")
            << "updateCoeffs(const scalarField& snGradp) MUST be called before"
               " updateCoeffs() or evaluate() to set the boundary gradient."
            << exit(FatalError);
    }
}

因此，显然，从 2.3 开始，只有在求解器中显式地调用带一个参数的 updateCoeffs 函数来指定 p 的边界条件时，该求解器的算例才能使用 fixedFluxPressure 这个边界条件。