OpenFOAM-2.3.x 中的 twoPhaseEulerFoam 解析之 kineticTheoryModel 的子模型

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

仅仅觉得都挺有用的，所以转发一下。仅仅是想做个备份，希望不要某天淹没在互联网的长河里
上一篇博文解读了 kineticTheoryModel 其中提到需要调用子模型来完成其功能，这里将 OpenFOAM 中 kineticTheoryModel 模型的子模型罗列如下。

1. viscosityModel

viscosityModel 的作用是根据颗粒温度 Theta 来计算固相粘度。
基类代码如下，核心是那个返回固相粘度的 nu 函数。

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
60
61
62
63
64
65
66
67
68
69
70
71
namespace Foam
{
namespace kineticTheoryModels
{

/*---------------------------------------------------------------------------*\
                           Class viscosityModel Declaration
\*---------------------------------------------------------------------------*/

class viscosityModel
{
    // Private member functions

        //- Disallow default bitwise copy construct
        viscosityModel(const viscosityModel&);

        //- Disallow default bitwise assignment
        void operator=(const viscosityModel&);
protected:

    // Protected data

        const dictionary& dict_;
public:

    //- Runtime type information
    TypeName("viscosityModel");

    // Declare runtime constructor selection table
    declareRunTimeSelectionTable
    (
        autoPtr,
        viscosityModel,
        dictionary,
        (
            const dictionary& dict
        ),
        (dict)
    );
    // Constructors

        //- Construct from components
        viscosityModel(const dictionary& dict);
    // Selectors

        static autoPtr<viscosityModel> New
        (
            const dictionary& dict
        );

    //- Destructor
    virtual ~viscosityModel();


    // Member Functions

        virtual tmp<volScalarField> nu
        (
            const volScalarField& alpha1,
            const volScalarField& Theta,
            const volScalarField& g0,
            const volScalarField& rho1,
            const volScalarField& da,
            const dimensionedScalar& e
        ) const = 0;

        virtual bool read()
        {
            return true;
        }
};

2.3.x 版自带四种固相粘度模型，分别如下：

1.1 none

顾名思义，这个模型计算的固相粘度值为零。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 注意这里的 kineticTheoryModels 不是类名，而是命名空间
Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::nu 
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    return dimensionedScalar
    (
        "0",
        dimensionSet(0, 2, -1, 0, 0, 0, 0),
        0.0
    )*alpha1; // 返回 0
}
1.2 Syamlal 模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::Syamlal::nu
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    return da*sqrt(Theta)*
    (
        (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1.0)*sqr(alpha1)/(3.0 - e)
      + (1.0/6.0)*alpha1*sqrtPi/(3.0 - e)
    );
}

公式为
$$
\nu_s = d_p \sqrt{\Theta}\left[ \frac{4}{5}\varepsilon_s^2 g_0 \frac{(1+e)}{\sqrt{\pi}} + \frac{1}{15} \sqrt{\pi} \cdot g_0 \varepsilon_s^2\frac{(1+e)(3e-1)}{\sqrt{\pi}} + \frac{1}{6}\varepsilon_s \frac{\sqrt{\pi}}{3-e}\right ]
$$
其中$g_0$是由 径向分布模型计算得到的。

1.3 HrenyaSinclair 模型
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
//- Characteristic length of geometry
        dimensionedScalar L_; // 新定义的一个变量，
        
 // 构造函数       
Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::HrenyaSinclair
(
    const dictionary& dict
)
:
    viscosityModel(dict),
    coeffDict_(dict.subDict(typeName + "Coeffs")),
    L_("L", dimensionSet(0, 1, 0, 0, 0), coeffDict_.lookup("L")) // 从外部读取 L_ 的值
{}

        
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::nu
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    volScalarField lamda
    (
        scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
    );

    return da*sqrt(Theta)*
    (
        (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1)*sqr(alpha1)/(3.0-e)
      + (1.0/6.0)*sqrtPi*alpha1*(0.5*lamda + 0.25*(3.0*e - 1.0))
       /(0.5*(3.0 - e)*lamda)
      + (10/96.0)*sqrtPi/((1.0 + e)*0.5*(3.0 - e)*g0*lamda)
    );
}

公式如下
$$
\begin{aligned}
\nu_s = & d_p \sqrt{\Theta} , [ \frac{4}{5} \varepsilon_s \cdot g_0 \frac{1+e}{\sqrt{\pi}} + \frac{1}{15}\sqrt{\pi} \cdot g_0 \varepsilon_s^2 \frac{(1+e)(3e-1)}{3-e} \
+ & \frac{1}{6} \sqrt{\pi} \cdot \frac{0.5\lambda + 0.25(3e-1)}{0.5(3-e)\lambda} \varepsilon_s + \frac{10}{96}\sqrt{\pi}\cdot \frac{1}{0.5(1+e)(3-e)g_0\cdot \lambda} ]
\end{aligned}
$$
其中
$$
\lambda = 1+\frac{d_p}{6\sqrt{2}\cdot \varepsilon_s}\cdot \frac{1}{L}
$$

1.4 Gidaspow 模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::viscosityModels::Gidaspow::nu
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    return da*sqrt(Theta)*
    (
        (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*sqr(alpha1)
      + (1.0/6.0)*sqrtPi*alpha1
      + (10.0/96.0)*sqrtPi/((1.0 + e)*g0)
    );
}

$$
\nu_s = d_p\sqrt{\Theta}\left [ \frac{4}{5} \varepsilon_s^2 g_0\cdot \frac{(1+e)}{\sqrt{\pi}} + \frac{1}{15} \sqrt{\pi}\cdot g_0(1+e)\varepsilon_s^2 + \frac{1}{6} \sqrt{\pi}\cdot \varepsilon_s^2 + \frac{10}{96} \frac{\sqrt{\pi}}{(1+e)g_0}\right ]
$$

2. radialModel

这个类的作用是计算径向分布函数 g0

有三种 radialModel 可以选择：

2.1 SinclairJackson
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
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::SinclairJackson::g0
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return 1.0/(1.0 - cbrt(min(alpha, alphaMinFriction)/alphaMax));
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    volScalarField aByaMax
    (
        cbrt(min(max(alpha, scalar(1e-3)), alphaMinFriction)/alphaMax)
    );

    return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax));
}

$$
g_0 = \frac{1}{1-\sqrt[3]{\frac{min(\varepsilon_s, \varepsilon_{s,min})}{\varepsilon_{s,max}}}}
$$

g0Prime 为 g0 对 alpha 的导数

$$
g_{0Prime} = \frac{\partial g_0}{\partial \varepsilon_s} = \frac{\frac{1}{3\cdot \varepsilon_{s,max}}}{(aByaMax-aByaMax^2)^2}
$$

其中

$$
aByaMax=\sqrt[3]{\frac{min(\varepsilon_s,\varepsilon_{s,min})}{\varepsilon_{s,max}}}
$$

2.2 LunSavage
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::LunSavage::g0
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{

    return pow(1.0 - alpha/alphaMax, -2.5*alphaMax);
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::LunSavage::g0prime
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return 2.5*pow(1.0 - alpha/alphaMax, -2.5*alphaMax - 1);
}

$$
g_0 = \left( 1-\frac{\varepsilon_s}{\varepsilon_{s,max}}\right)^{-2.5,\varepsilon_{s,max}}
$$

2.3 CarnahanStarling 模型
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
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{

    return
        1.0/(1.0 - alpha)
      + 3.0*alpha/(2.0*sqr(1.0 - alpha))
      + sqr(alpha)/(2.0*pow3(1.0 - alpha));
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0prime
(
    const volScalarField& alpha,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        2.5/sqr(1.0 - alpha)
      + 4.0*alpha/pow3(1.0 - alpha)
      + 1.5*sqr(alpha)/pow4(1.0 - alpha);
}

$$
g_0 = \frac{1}{1-\varepsilon_s} + \frac{3\varepsilon_s}{2(1-\varepsilon_s)^2} + \frac{\varepsilon_s^2}{2(1-\varepsilon_s)^3}
$$

3. granularPressureModel

顾名思义，这个类是用来计算固相压力的。

OpenFOAM 内置两种固相压力模型：

3.1 Lun 模型

granularPressureCoeff 返回的是固相压力的系数，这个返回值乘以颗粒温度 Theta 才是固相压力。

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
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::Lun::granularPressureCoeff
(
    const volScalarField& alpha1,
    const volScalarField& g0,
    const volScalarField& rho1,
    const dimensionedScalar& e
) const
{

    return rho1*alpha1*(1.0 + 2.0*(1.0 + e)*alpha1*g0);
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::Lun::
granularPressureCoeffPrime
(
    const volScalarField& alpha1,
    const volScalarField& g0,
    const volScalarField& g0prime,
    const volScalarField& rho1,
    const dimensionedScalar& e
) const
{
    return rho1*(1.0 + alpha1*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha1));
}

$$
P_{s,coeff} = \rho\varepsilon_s[1+2(1+e)\varepsilon_sg_0]
$$

granularPressureCoeffPrime 函数计算的是 $\partial P_{s,coeff}/\partial \varepsilon_s$ 。

3.2 SyamlalRogersOBrien 模型
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
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien::
granularPressureCoeff
(
    const volScalarField& alpha1,
    const volScalarField& g0,
    const volScalarField& rho1,
    const dimensionedScalar& e
) const
{

    return 2.0*rho1*(1.0 + e)*sqr(alpha1)*g0;
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::granularPressureModels::SyamlalRogersOBrien::
granularPressureCoeffPrime
(
    const volScalarField& alpha1,
    const volScalarField& g0,
    const volScalarField& g0prime,
    const volScalarField& rho1,
    const dimensionedScalar& e
) const
{
    return rho1*alpha1*(1.0 + e)*(4.0*g0 + 2.0*g0prime*alpha1);
}

$$
P_{s,coeff} = 2\rho(1+e)\varepsilon_s^2\cdot g_0
$$

4. frictionalStressModel

在稠密气固两相流中，当固相体积分率大于某个值时，单纯考虑跟颗粒温度关联的固相压力和固相粘性还不够，还需要考虑所谓的摩擦应力。这个类就是用来计算摩擦应力的。

有两种模型可选：

4.1 Schaeffer 模型
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
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::Schaeffer
(
    const dictionary& dict
)
:
    frictionalStressModel(dict),
    coeffDict_(dict.subDict(typeName + "Coeffs")),
    phi_("phi", dimless, coeffDict_.lookup("phi"))
{
    phi_ *= constant::mathematical::pi/180.0; 
    // 这个phi_是一个角度，从外部读取，在外部设置的时候，按角度的单位来设置，这里是将角度转换成弧度。
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::~Schaeffer()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressure
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        dimensionedScalar("1e24", dimensionSet(1, -1, -2, 0, 0), 1e24)
       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 10.0);
}

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::
frictionalPressurePrime
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return
        dimensionedScalar("1e25", dimensionSet(1, -1, -2, 0, 0), 1e25)
       *pow(Foam::max(alpha1 - alphaMinFriction, scalar(0)), 9.0);
}

$$
P_{f} = 10^{24}\cdot max(\varepsilon_s-\varepsilon_{s,friMin},0)^{10}
$$
同前面一样，frictionalPressurePrime 是 frictionalPressure 对固相体积分率的导数。

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
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMax,
    const volScalarField& pf,
    const volSymmTensorField& D
) const
{
    const scalar I2Dsmall = 1.0e-15;

    // Creating nu assuming it should be 0 on the boundary which may not be
    // true
    tmp<volScalarField> tnu
    (
        new volScalarField
        (
            IOobject
            (
                "Schaeffer:nu",
                alpha1.mesh().time().timeName(),
                alpha1.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            alpha1.mesh(),
            dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0.0)
        )
    );

    volScalarField& nuf = tnu();

    forAll (D, celli)
    {
        if (alpha1[celli] > alphaMax.value() - 5e-2)
        {
            nuf[celli] =
                0.5*pf[celli]*sin(phi_.value())
               /(
                    sqrt(1.0/6.0*(sqr(D[celli].xx() - D[celli].yy())
                  + sqr(D[celli].yy() - D[celli].zz())
                  + sqr(D[celli].zz() - D[celli].xx()))
                  + sqr(D[celli].xy()) + sqr(D[celli].xz())
                  + sqr(D[celli].yz())) + I2Dsmall
                );
        }
    }

    // Correct coupled BCs
    nuf.correctBoundaryConditions();

    return tnu;
}

$$
\nu_f = \left { \begin{aligned}
& 0 , &\varepsilon_s \le \varepsilon_{s,max} \
& \frac{0.5p_f sin\phi}{\sqrt{I_{2d}}}, & \varepsilon_s > \varepsilon_{s,max}
\end{aligned} \right .
$$

其中，$p_f$ 代表的是上面的 frictionalPressure
$$
D = \frac{1}{2}(\nabla U + \nabla U^T)
$$

$$
\sqrt{I_{2D}} = \frac{1}{6}[(D_{11}-D_{22})^2 + (D_{22}-D_{33})^2 + (D_{33}-D_{11})^2] + D_{12}^2 + D_{13}^2 + D_{23}^2
$$

4.2 JohnsonJackson 模型
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
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
class JohnsonJackson
:
    public frictionalStressModel
{
    // Private data

        dictionary coeffDict_;

        //- Material constant for frictional normal stress
        dimensionedScalar Fr_;

        //- Material constant for frictional normal stress
        dimensionedScalar eta_;

        //- Material constant for frictional normal stress
        dimensionedScalar p_;

        //- Angle of internal friction
        dimensionedScalar phi_;
        
        
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
JohnsonJackson
(
    const dictionary& dict
)
:
    frictionalStressModel(dict),
    coeffDict_(dict.subDict(typeName + "Coeffs")),
    Fr_("Fr", dimensionSet(1, -1, -2, 0, 0), coeffDict_.lookup("Fr")),
    eta_("eta", dimless, coeffDict_.lookup("eta")),
    p_("p", dimless, coeffDict_.lookup("p")),
    phi_("phi", dimless, coeffDict_.lookup("phi"))
{
    phi_ *= constant::mathematical::pi/180.0;
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
frictionalPressure
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{

    return
        Fr_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
       /pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_);
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::
frictionalPressurePrime
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMinFriction,
    const dimensionedScalar& alphaMax
) const
{
    return Fr_*
    (
        eta_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_ - 1.0)
       *(alphaMax-alpha1)
      + p_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
    )/pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_ + 1.0);
}


Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::nu
(
    const volScalarField& alpha1,
    const dimensionedScalar& alphaMax,
    const volScalarField& pf,
    const volSymmTensorField& D
) const
{
    return dimensionedScalar("0.5", dimTime, 0.5)*pf*sin(phi_);
}

这个模型的特点是需要设置一些跟材料有关的参数。
$$
P_{f} = F_r \frac{max(\varepsilon_s-\varepsilon_{s,friMin},0)^{\eta}}{max(\varepsilon_{s,max}-\varepsilon_{s},0)^{p}}
$$

$$
\nu_f = 0.5\cdot p_fsin\phi
$$
一些材料的物性参数建议值可参见”Derivation, Implementation and Validation of Computer Simulation Models for Gas-Solids Fluidized Beds” Table-3.5。

5. conductivityModel

这个类的作用是计算 颗粒温度方程中的颗粒温度传导系数，只有在使用偏微分方程求解颗粒温度是才会用到。

有三种可选：

5.1 Syamlal 模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::Syamlal::kappa
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    return rho1*da*sqrt(Theta)*
    (
        2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (9.0/8.0)*sqrtPi*g0*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
       /(49.0/16.0 - 33.0*e/16.0)
      + (15.0/32.0)*sqrtPi*alpha1/(49.0/16.0 - 33.0*e/16.0)
    );
}

$$
\kappa = \rho d_p \sqrt{\Theta}\left [ 2 \varepsilon_s^2 g_0 \frac{1+e}{\sqrt{\pi}} + \frac{\frac{9}{8} \sqrt{\pi}g_0 \cdot 0.25(1+e)^2(2e-1)\varepsilon_s^2}{49/16-33e/16} + \frac{\frac{15}{32}\sqrt{\pi} \cdot \varepsilon_s }{49/16-33e/16} \right ]
$$

5.2 HrenyaSinclair 模型
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
Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::HrenyaSinclair
(
    const dictionary& dict
)
:
    conductivityModel(dict),
    coeffDict_(dict.subDict(typeName + "Coeffs")),
    L_("L", dimensionSet(0, 1, 0, 0, 0), coeffDict_.lookup("L"))
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::
~HrenyaSinclair()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    volScalarField lamda
    (
        scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
    );

    return rho1*da*sqrt(Theta)*
    (
        2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (9.0/8.0)*sqrtPi*g0*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
       /(49.0/16.0 - 33.0*e/16.0)
      + (15.0/16.0)*sqrtPi*alpha1*(0.5*sqr(e) + 0.25*e - 0.75 + lamda)
       /((49.0/16.0 - 33.0*e/16.0)*lamda)
      + (25.0/64.0)*sqrtPi
       /((1.0 + e)*(49.0/16.0 - 33.0*e/16.0)*lamda*g0)
    );
}

这个模型需要输入一个特征长度 L。
$$
\begin{aligned} \kappa = \ & \rho d_p \sqrt{\Theta} , [ 2\varepsilon_s^2 g_0 \frac{1+e}{\sqrt{\pi}} + \frac{\frac{9}{8} \sqrt{\pi}g_0 \cdot 0.25(1+e)^2(2e-1)\varepsilon_s^2}{49/16-33e/16} \
\ + &\frac{\frac{15}{16}\sqrt{\pi} \varepsilon_s \cdot (0.5e^2+0.25e-0.75+\lambda) }{(49/16-33e/16)\lambda} + \frac{\frac{25}{64}\sqrt{\pi}}{(1+e)(49/16-33e/16)\lambda g_0}
] \end{aligned}
$$

其中
$$
\lambda = 1+\frac{d_p}{6\sqrt{2}\varepsilon_s L}
$$

5.3 Gidaspow 模型
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa
(
    const volScalarField& alpha1,
    const volScalarField& Theta,
    const volScalarField& g0,
    const volScalarField& rho1,
    const volScalarField& da,
    const dimensionedScalar& e
) const
{
    const scalar sqrtPi = sqrt(constant::mathematical::pi);

    return rho1*da*sqrt(Theta)*
    (
        2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
      + (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha1)
      + (15.0/16.0)*sqrtPi*alpha1
      + (25.0/64.0)*sqrtPi/((1.0 + e)*g0)
    );
}

$$
\kappa = \rho d_p \sqrt{\Theta}\left [ 2 \varepsilon_s^2 g_0 \frac{1+e}{\sqrt{\pi}} + \frac{9}{8}\sqrt{\pi}g_0\cdot 0.5(1+e)\varepsilon_s^2 + \frac{15}{16}\sqrt{\pi} \varepsilon_s +\frac{25}{64}\frac{\sqrt{\pi}}{(1+e)g_0}\right ]
$$