OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之曳力模型的调用过程

转自：https://xiaopingqiu.github.io/，侵删

仅仅觉得都挺有用的，所以转发一下。仅仅是想做个备份，希望不要某天淹没在互联网的长河里
前面有三篇博文对 OpenFOAM-2.1.x 中的 twoPhaseEulerFoam 求解器进行了解读，然而 OpenFOAM-2.3.x 中，这个求解器的代码有了很大的变化。本文将以一个曳力模型的调用过程为例，介绍 OpenFOAM-2.3.x 中 twoPhaseEulerFoam 是如何调用相间作用力模型的。后续还将对 OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 的其他方面进行解读。

主程序中(“ UEqn.H “)，dragCoeff 定义为：

1
volScalarField dragCoeff(fluid.dragCoeff());

其中，fluid 为 twoPhaseSystem 类的对象，所以，要去找 twoPhaseSystem 类的成员函数 dragCoeff()。
在源文件 twoPhaseSystem.C 中找到如下定义：

1
2
3
4
Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::dragCoeff() const
{
    return drag_->K();
}

而 drag_ 的定义为

1
autoPtr<BlendedInterfacialModel<dragModel> > drag_;

所以，需要找到类 BlendedInterfacialModel<dragModel> 的成员函数 K()的定义。
在源文件 BlendedInterfacialModel.C 中，找到如下定义：

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
53
template<class modelType>
Foam::tmp<Foam::volScalarField>
Foam::BlendedInterfacialModel<modelType>::K() const
{
    tmp<volScalarField> f1, f2;

    if (model_.valid() || model1In2_.valid())
    {
        f1 = blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed());
    }

    if (model_.valid() || model2In1_.valid())
    {
        f2 = blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed());
    }

    tmp<volScalarField> x
    (
        new volScalarField
        (
            IOobject
            (
                modelType::typeName + "Coeff",
                pair_.phase1().mesh().time().timeName(),
                pair_.phase1().mesh()
            ),
            pair_.phase1().mesh(),
            dimensionedScalar("zero", modelType::dimK, 0)
        )
    );

    if (model_.valid())
    {
        x() += model_->K()*(f1() - f2());
    }

    if (model1In2_.valid())
    {
        x() += model1In2_->K()*(1 - f1);
    }

    if (model2In1_.valid())
    {
        x() += model2In1_->K()*f2;
    }

    if (model_.valid() || model1In2_.valid() || model2In1_.valid())
    {
        correctFixedFluxBCs(x());
    }

    return x;
}

对于曳力模型，上述成员函数的 modelType 可以实例化为 dragModel，要理解该函数的行为，有三点需要清楚： pair_， pair1In2_， pair2In1_ 的定义； model_， model1In2_， model2In1_ 的定义； blendingMethod 类的成员函数 f1 与 f2 的定义。下面一一解释：

1. pair_， pair1In2_， pair2In1_ 的定义

这三个是 BlendedInterfacialModel 类的数据成员，回到 twoPhaseSystem 类中去看 drag_ 的初始化，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
drag_.set
   (
       new BlendedInterfacialModel<dragModel>
       (
           lookup("drag"),
           (
               blendingMethods_.found("drag")
             ? blendingMethods_["drag"]
             : blendingMethods_["default"]
           ),
           pair_,
           pair1In2_,
           pair2In1_
       )
   );

可知，BlendedInterfacialModel 类中的 pair_， pair1In2_， pair2In1_ 是将 twoPhaseSystem 类的数据成员传递过去来实现初始化的，所以，真正要看懂的是twoPhaseSystem类中数据成员pair_， pair1In2_， pair2In1_ 的初始化，见如下代码：

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
pair_.set
    (
        new phasePair
        (
            phase1_,
            phase2_,
            g,
            sigmaTable
        )
    );
    pair1In2_.set
    (
        new orderedPhasePair
        (
            phase1_,
            phase2_,
            g,
            sigmaTable,
            aspectRatioTable
        )
    );
    pair2In1_.set
    (
        new orderedPhasePair
        (
            phase2_,
            phase1_,
            g,
            sigmaTable,
            aspectRatioTable
        )
    );

可见，pair_ 是 phasePair 类的指针， pair1In2_ 与 pair2In1_ 是 orderedPhasePair 类的指针。
其中phase1_ 和 phase2_ 是通过从文件”phaseProperties”里读取内容来初始化的：

1
2
3
4
5
6
7
8
9
10
11
12
phase1_
    (
        *this, 
        *this,
        wordList(lookup("phases"))[0]
    ),
    phase2_
    (
        *this,
        *this,
        wordList(lookup("phases"))[1]
    ),

举例说，假设”phaseProperties” 文件里有以下内容：

1
phases (particles air);

则， phase1_ = “particles”， phase2_ = “air” 。

根据 phasePair 类中的定义，成员函数 dispersed() 总是返回对象的phase1（也就是 phasePair 或者 orderedPhasePair 类的构造函数的第一个参数），所以，对于 “particles air” 体系，pair1In2_.dispersed() = phase1_.name() = "particles"， 而 pair2In1_.dispersed() = phase2_.name() = "air"。

2. model_， model1In2_， model2In1_ 的定义

这三个是 BlendedInterfacialModel 类的数据成员，定义和初始化如下：

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
autoPtr<modelType> model_;
autoPtr<modelType> model1In2_;
autoPtr<modelType> model2In1_;

if (modelTable.found(pair_))
    {
        model_.set
        (
            modelType::New
            (
                modelTable[pair_],
                pair_
            ).ptr()
        );
    }

    if (modelTable.found(pair1In2_))
    {
        model1In2_.set
        (
            modelType::New
            (
                modelTable[pair1In2_],
                pair1In2_
            ).ptr()
        );
    }

    if (modelTable.found(pair2In1_))
    {
        model2In1_.set
        (
            modelType::New
            (
                modelTable[pair2In1_],
                pair2In1_
            ).ptr()
        );
    }
}

注意，这里讨论的是曳力模型的调用，所以，如前所述，modelType 可以实例化为 dragModel。
modelTable 是 phasePair::dictTable 类型的引用，本质上是一个 HashTable（ HashTable<dictionary, phasePairKey, phasePairKey::hash> ），其 key 是 phasePairKey 类型的对象，value 是 dictionary 类的对象。 found 函数通过查找 modelTable 对象中是否存在某个 key 来决定返回值是 true 还是 false。

这里要分头说，一边是 modelTable 的初始化，另一边是 pair_， pair2In1_， pair2In1_ 如何与 phasePairKey 类进行对比。

从 twoPhaseSystem 类中对 drag_ 的初始化可知， modelTable 的初始化是由 lookup("drag") 来完成的。lookup 函数的作用是读取”phaseProperties” 文件的内容来实现对一个 HashTable 的初始化（具体过程将会在后续解读中涉及）。举例说，以下 “phaseProperties” 的内容

1
2
3
4
5
6
7
8
9
10
11
12
13
drag
(
    (particles in air)
    {
        type            GidaspowErgunWenYu;
        residualAlpha   1e-6;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }
);

将利用 (particles in air) 来初始化一个 phasePairKey 对象（利用 phasePairKey类中的空白构造函数和重载的 >> 符号）。成员ordeded_ 的值取决于 “in” 或 “and” ，若形如 “particles in air “，ordeded_ = true，若形如 “particles and air “， 则 ordeded_ = false。 则剩余内容将用于初始化一个 dictionary 对象。

而另一方面，phasePair 是 phasePairKey的派生类， orderedPhasePair 则是phasePair 的派生类，所以，将 pair_ ，pair1In2_ 以及 pair2In1_ 作为 found 函数的参数，隐含了将派生类的引用转换成基类引用。 phasePair类默认 ordered_ = false， 而 orderedPhasePair 类则默认ordered_ = true。pair_ ，pair1In2_ 以及 pair2In1_ 与 modelTable_ 的 key 进行比较，比较的是对应的 phase1，phase2 和 ordered_ 三个成员的值是否相等，只有三者都一样时， found 函数才返回 true 。 所以，对于上面提到的设置，即

1
2
3
4
5
6
phases (particles air);
drag
(
    (particles in air)
    .......
);

只有modelTable.found(pair1In2_)的值为true。同样，也就只有 model1In2_.valid() 为 true （即 model1In2_ 指针不为空。）

3. blendingMethod 类的成员函数 f1 与 f2

这两个函数的实现在不同的 blendingMethods中不一样，以最简单的 noBlending 类型为例：

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
53
54
55
56
57
58
59
Foam::tmp<Foam::volScalarField> Foam::blendingMethods::noBlending::f1
(
    const phaseModel& phase1,
    const phaseModel& phase2
) const
{
    const fvMesh& mesh(phase1.mesh());

    return
        tmp<volScalarField>
        (
            new volScalarField
            (
                IOobject
                (
                    "f",
                    mesh.time().timeName(),
                    mesh
                ),
                mesh,
                dimensionedScalar
                (
                    "f",
                    dimless,
                    phase2.name() != continuousPhase_
                ) // 如果 phase2 就是 continuousPhase，那么 f1 = 0；否则 f1 = 1
            )
        );
}

Foam::tmp<Foam::volScalarField> Foam::blendingMethods::noBlending::f2
(
    const phaseModel& phase1,
    const phaseModel& phase2
) const
{
    const fvMesh& mesh(phase1.mesh());

    return
        tmp<volScalarField>
        (
            new volScalarField
            (
                IOobject
                (
                    "f",
                    mesh.time().timeName(),
                    mesh
                ),
                mesh,
                dimensionedScalar
                (
                    "f",
                    dimless,
                    phase1.name() == continuousPhase_
                ) // 如果 phase1 是 continuousPhase，那么 f2 = 1，否则 f2 = 0。
            )
        );
}

再回头看 BlendedInterfacialModel 类的成员函数 K()，

1
2
3
4
5
6
7
8
9
if (model_.valid() || model1In2_.valid())
   {
       f1 = blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed());
   }

   if (model_.valid() || model2In1_.valid())
   {
       f2 = blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed());
   }

根据上面的 “phaseProperties” 的设置，可知 pair1In2_.dispersed() = "particles"， pair2In1_.dispersed() = "air"，而 continuousPhase_ 是从 “phaseProperties” 的 “blending” 子字典里读取的，这里continuousPhase_ = "air"，于是，可以得到 f1 = 0， f2 = 1。再看 K() 的返回值，可知，最终有效的返回值是

1
2
3
4
if (model1In2_.valid())
   {
       x() += model1In2_->K()*(1 - f1);
   }

即，最终 K() 函数的返回值是 model1In2_->K() 。而model1In2_->K()的值就取决于具体调用的曳力模型了，举例说，假如调用的是 Ergun 曳力模型，则 K() 最终返回的值，也就是 “UEqn.H” 中的 dragCoeff 的值是

1
2
3
4
5
6
7
8
9
10
11
Foam::tmp<Foam::volScalarField> Foam::dragModel::K() const
{
    return
        0.75
       *CdRe()
       *max(pair_.dispersed(), residualAlpha_)
       *swarmCorrection_->Cs()
       *pair_.continuous().rho()
       *pair_.continuous().nu()
       /sqr(pair_.dispersed().d());
}

其中 CdRe() 的定义为

1
2
3
4
5
6
7
8
9
10
11
12
Foam::tmp<Foam::volScalarField> Foam::dragModels::Ergun::CdRe() const
{
    return
        (4/3)
       *(
            150
           *max(scalar(1) - pair_.continuous(), residualAlpha_)
           /max(pair_.continuous(), residualAlpha_)
          + 1.75
           *pair_.Re()
        );
}

注意所有曳力模型的 K()函数形式是一样的，不同曳力模型的区别在于 CdRe() 的实现不一样。
此外，virtualMass， heatTransfer，lift，wallLubrication，turbulentDispersion 这些子模型的调用也都是经过类似的过程进行的。