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

上一篇博文解读了 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}}}}
$$

g0Primeg0alpha 的导数

$$
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}
$$
同前面一样,frictionalPressurePrimefrictionalPressure 对固相体积分率的导数。

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 ]
$$