OpenFOAM中重力的植入方式
考虑重力的NS方程可以写为:
ρ ∂ u ⃗ ∂ t + ρ ∇ ( u ⃗ u ⃗ ) = ∇ ( μ ∇ u ⃗ ) − ∇ P + ρ g ⃗ (1) \rho \frac{\partial \vec u}{\partial t}+\rho \nabla(\vec u \vec u)=\nabla(\mu\nabla \vec u)-\nabla P +\rho \vec g \tag {1} ρ∂t∂u+ρ∇(uu)=∇(μ∇u)−∇P+ρg(1)在openFOAM中,许多求解器中都耦合了重力,以interFoam为例,其中的重力定义处的代码块为:
fvc::reconstruct
(
(
mixture.surfaceTensionForce()
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
将其中的重力抽离出来,其表达式为fvc::reconstruct( (ghf*fvc::snGrad(rho) )* mesh.magSf() )
。通过溯源代码,我们可以进一步得到其中ghf的表达式:
new surfaceScalarField
(
"ghf",
(gFluid[i] & fluidRegions[i].Cf()) - ghRef
)
其中gFluid[i]
返回了通过字典(constant/g)读取得到的重力加速度项,fluidRegions[i].Cf()
返回了当前网格表面的面心矢量,ghRef是给定的相对重力加速度,一般为0。综上所述,openFOAM中的重力的表达式为:
G
⃗
=
−
(
g
⃗
⋅
s
⃗
)
∇
ρ
(2)
\vec G=-(\vec g \cdot \vec s)\nabla \rho\tag {2}
G=−(g⋅s)∇ρ(2)其中
s
⃗
\vec s
s是待求点的位置矢量,这显然与方程1使用的重力表达式相差很大,接下来我们就一步一步推导这个结果。
首先明确,在涉及到重力的求解器中,其求解的压力p并不是总压P,而是动压,即:
P
=
p
+
ρ
g
⃗
⋅
s
⃗
(3)
P=p+\rho \vec g \cdot \vec s\tag {3}
P=p+ρg⋅s(3)之所以要这样处理是为了方便用户设置初场,大家可以设想最简单的溃坝算例,如果直接对总压设置初场,那我们要设置的则是一个梯形的压力分布,计算域的y坐标越高其净水压强就越大,而如果我们拆分成动压和静压的形式,就可以直接设置初场为0了。可能有小伙伴对
g
⃗
⋅
s
⃗
\vec g \cdot \vec s
g⋅s项不太理解,其实我们只要对其展开就很好说明问题了:
g
⃗
⋅
s
⃗
=
(
0
,
g
,
0
)
⋅
(
x
,
y
,
z
)
=
g
y
(4)
\vec g \cdot \vec s=(0,g,0)\cdot(x,y,z)=gy\tag {4}
g⋅s=(0,g,0)⋅(x,y,z)=gy(4)在动量方程中,我们存在对总压的压力梯度项,将总压拆成静压和动压即可得到:
∇
P
=
∇
(
p
+
ρ
g
⃗
⋅
s
⃗
)
=
∇
p
+
∇
(
ρ
g
⃗
⋅
s
⃗
)
=
∇
p
+
(
g
⃗
⋅
s
⃗
)
∇
ρ
+
ρ
∇
(
g
⃗
⋅
s
⃗
)
(5)
\nabla P=\nabla (p+\rho \vec g \cdot \vec s)=\nabla p+\nabla(\rho \vec g \cdot \vec s)=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho\nabla( \vec g \cdot \vec s)\tag {5}
∇P=∇(p+ρg⋅s)=∇p+∇(ρg⋅s)=∇p+(g⋅s)∇ρ+ρ∇(g⋅s)(5)其中关键点在于处理
∇
(
g
⃗
⋅
s
⃗
)
\nabla(\vec g \cdot \vec s)
∇(g⋅s)项,根据式4有:
∇
(
g
⃗
⋅
s
⃗
)
=
∇
(
g
y
)
=
(
∂
(
g
y
)
∂
x
,
∂
(
g
y
)
∂
y
,
∂
(
g
y
)
∂
z
)
=
(
0
,
g
,
0
)
=
g
⃗
(6)
\nabla(\vec g \cdot \vec s)=\nabla (gy)=(\frac{\partial (gy)}{\partial x},\frac{\partial (gy)}{\partial y},\frac{\partial (gy)}{\partial z})=(0,g,0)=\vec g\tag {6}
∇(g⋅s)=∇(gy)=(∂x∂(gy),∂y∂(gy),∂z∂(gy))=(0,g,0)=g(6)因此,式5最终写为:
∇
P
=
∇
p
+
(
g
⃗
⋅
s
⃗
)
∇
ρ
+
ρ
∇
(
g
⃗
⋅
s
⃗
)
=
∇
p
+
(
g
⃗
⋅
s
⃗
)
∇
ρ
+
ρ
g
⃗
(7)
\nabla P=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho\nabla( \vec g \cdot \vec s)=\nabla p+ ( \vec g \cdot \vec s)\nabla\rho+\rho \vec g\tag {7}
∇P=∇p+(g⋅s)∇ρ+ρ∇(g⋅s)=∇p+(g⋅s)∇ρ+ρg(7)将式7带回方程1,即可得到OpenFOAM代码中使用的含重力的动量方程了:
ρ
∂
u
⃗
∂
t
+
ρ
∇
(
u
⃗
u
⃗
)
=
∇
(
μ
∇
u
⃗
)
−
∇
p
−
(
g
⃗
⋅
s
⃗
)
∇
ρ
(8)
\rho \frac{\partial \vec u}{\partial t}+\rho \nabla(\vec u \vec u)=\nabla(\mu\nabla \vec u)-\nabla p - ( \vec g \cdot \vec s)\nabla\rho \tag {8}
ρ∂t∂u+ρ∇(uu)=∇(μ∇u)−∇p−(g⋅s)∇ρ(8)