本文介绍了标量运输求解器scalarTransPortFoam
基本介绍
标量输运方程
$$
\frac{\partial}{\partial t}\int_V\rho\phi dV+\int_S\rho\phi\mathbf v\cdot\mathbf ndS=\sum f_\phi\tag{1}
$$
其中
$$
f_\phi^d=\int_S\Gamma\nabla\phi\cdot\mathbf ndS\tag{2}
$$
$$
\frac{\partial \rho\phi}{\partial t}+\nabla\cdot(\rho\phi\mathbf v)=\nabla\cdot(\Gamma \nabla\phi)\tag{3}
$$
库朗数的计算
- Courant–Friedrichs–Lewy condition：简称CFL条件，在许多显式的时间行进模拟中，时间步长必须小于特定时间，否则模拟会产生错误的结果
- 表达式
$$
C=\frac{u\Delta t}{\Delta x}\leqslant C_{max}\tag{4}
$$
$$
C=\Delta t(\sum_{i=1}^n\frac{u_{x_i}}{\Delta x_i})\leqslant C_{max}\tag{5}
$$
- 在OpenFOAM中，定义为：
$$
\Co=0.5\frac{\Delta t\sum_f|\phi_f|}{\Delta V}\tag{6}
$$
其中$\phi_f$表示网格单元面f的通量，$\Delta V$表示网格单元体积。方程(6)也被称之为面库朗数。由于通量守恒，进入网格单元的通量等于流出网格单元的通量，因此在计算面库朗数的时候，要乘以0.5。
程序详读
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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
| \*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fvIOoptionList.H"
#include "simpleControl.H"
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
-----------------------------------createFields.H----------------------------
Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Info<< "Reading diffusivity DT\n" << endl;
dimensionedScalar DT
(
transportProperties.lookup("DT")
);
-----------------------------------createFields.H----------------------------
-----------------------------------createPhi.H-------------------------------
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(U) & mesh.Sf()
);
-----------------------------------createPhi.H-------------------------------
-----------------------------------createFvOptions.H-------------------------
fv::IOoptionList fvOptions(mesh);
-----------------------------------createFvOptions.H-------------------------
simpleControl simple(mesh);
Info<< "\nCalculating scalar transport\n" << endl;
-----------------------------------CourantNo.H-------------------------------
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
if (mesh.nInternalFaces())
{
scalarField sumPhi
(
fvc::surfaceSum(mag(phi))().internalField()
);
CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum << endl;
-----------------------------------CourantNo.H-------------------------------
while (simple.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
while (simple.correctNonOrthogonal())
{
solve
(
fvm::ddt(T)
+ fvm::div(phi, T)
- fvm::laplacian(DT, T)
==
fvOptions(T)
);
}
runTime.write();
}
Info<< "End\n" << endl;
return 0;
}
|