本文总结OpenFOAM中代码语句与模型或公式的对应关系
icoFoam solver
连续性方程
∇
⋅
U
=
0
\nabla \cdot U=0
∇⋅U=0
动量方程
∂
U
∂
t
+
∇
⋅
(
U
U
)
=
−
∇
p
+
∇
⋅
(
ν
∇
U
)
\frac{\partial U}{\partial t}+\nabla\cdot (UU)=-\nabla p + \nabla \cdot (\nu \nabla U)
∂t∂U+∇⋅(UU)=−∇p+∇⋅(ν∇U)
fvm::ddt(U)
返回fvVectorMatrix，返回N*N的矩阵，矩阵中每个元素都是括号内的U，因此是Vector。
- 与动量方程对应关系
∂ U ∂ t \frac{\partial U}{\partial t} ∂t∂U - 与离散方程对应关系
∂ ∂ t ∫ V P U d V ≈ ∂ ∂ t ( U P V P ) = a ⋅ U P n + 1 + b ⋅ U P n \frac{\partial }{\partial t}\int_{V_P} U dV \approx \frac{\partial }{\partial t} (U_P V_P)=a\cdot U_P^{n+1} +b\cdot U_P^n ∂t∂∫VPUdV≈∂t∂(UPVP)=a⋅UPn+1+b⋅UPn
最后一个等号所得的形式取决于定义在fvSchemes中的ddtSchemes离散格式
Field和Matrix的关系是什么？Field为(N,1)的向量，其中N为网格数。Matrix顾名思义为(N,N)的矩阵。Field为待求量，我们想按照AX=B的计算格式来求解Field X，由于B同样为Field，这就要求A为Matrix or Scalar
fvm::div(phi, U)
返回fvVectorMatrix
- 与动量方程对应关系
∇ ⋅ ( U U ) \nabla\cdot (UU) ∇⋅(UU) - 与离散方程对应关系
∫ V P ∇ ⋅ ( U U ) d V = ∫ ∂ V P ( U U ) d S ≈ ∑ ( U U ) f S ≈ ∑ F f n U f n + 1 = a ⋅ U P n + 1 + ∑ b N U N n + 1 \int_{V_P}\nabla\cdot (UU) dV=\int_{\partial V_P} (UU) dS\approx\sum(UU)_f S \approx \sum F^n_f U^{n+1}_f = a\cdot U_P^{n+1} + \sum b_N U_N^{n+1} ∫VP∇⋅(UU)dV=∫∂VP(UU)dS≈∑(UU)fS≈∑FfnUfn+1=a⋅UPn+1+∑bNUNn+1
因此， d i v ( a , b ) \bf{div(a,b)} div(a,b)操作对应的运算符为 ∇ ⋅ ( a b ) \nabla \cdot (ab) ∇⋅(ab)，对应的离散为 ∑ ( a b ) \sum (a b) ∑(ab)
上述离散结果返回Matrix，该矩阵对角线元素为上式中最右端项中的a
第四项到第五项的过渡取决于如何将 U f U_f Uf转化为体心值 U P U_P UP，其离散方法定义在fvSchemes/divSchemes
fvm::laplacian(nu, U)
返回fvVectorMatrix
- 与动量方程对应关系
∇ ⋅ ( ν ∇ U ) \nabla \cdot (\nu \nabla U) ∇⋅(ν∇U) - 与离散方程对应关系
∫ V P ∇ ⋅ ( ν ∇ U ) d V = ∫ ∂ V P ( ν ∇ U ) d S ≈ ∑ ν ( ∇ U n + 1 ) f S = ∑ ν [ ( ∇ U n + 1 ) f ⋅ S S f ] S f = a ⋅ U P n + 1 + ∑ b N U N n + 1 \int_{V_P}\nabla \cdot (\nu \nabla U) dV = \int_{\partial V_P} (\nu \nabla U) dS \approx \sum \nu (\nabla U^{n+1})_f S = \sum \nu \left[ (\nabla U^{n+1})_f \cdot \frac{S}{S_f} \right]S_f = a\cdot U_P^{n+1} + \sum b_N U_N^{n+1} ∫VP∇⋅(ν∇U)dV=∫∂VP(ν∇U)dS≈∑ν(∇Un+1)fS=∑ν[(∇Un+1)f⋅SfS]Sf=a⋅UPn+1+∑bNUNn+1
其中 n = S / S f n=S/S_f n=S/Sf为单位矢量， S f S_f Sf为标量，值为面积，上式第四项对应代码 f v m : : s n G r a d ( U ) \bf{fvm::snGrad(U)} fvm::snGrad(U)
上式第四项到第五项的过渡取决于面上梯度 ( ∇ U ) f (\nabla U)_f (∇U)f的离散方式，即如何将面上梯度转化为体心值，定义在
fvSchemes/laplacianSchemes
fvc::grad(p)
返回scalarField，该返回类型是fvVectorMatrix的子集，执行solve(UEqn == -fvc::grad(p))时将后者做了类型强制转换
- 与动量方程对应关系
∇ p \nabla p ∇p - 与离散方程对应关系
∫ V P ∇ p d V = ∫ ∂ V P p d S ≈ ∑ p f n S = p P n + ∑ p N n \int_{V_P}\nabla p dV = \int_{\partial V_P} p dS \approx \sum p^n_f S = p_P^n+\sum p_N^n ∫VP∇pdV=∫∂VPpdS≈∑pfnS=pPn+∑pNnfvSchemes/gradSchemes决定了如何由面上值 p f p_f pf插值得到体心值 p P p_P pP
solve(UEqn == -fvc::grad(p))
求解括号内以fvMatrix作为系数矩阵的AX=B，很显然UEqn为A，双等号右边为B，两者相等解出X，即U。另外，该语句没有返回值，但是会将所求解的场量更新，即更新U
fvc::flux(HbyA)
返回surfaceScalarField，其大小为内部面的个数，等于constant/polyMesh/neighbour文件中元素的个数。HbyA为volVectorField，其大小为网格数，可见该操作是将体心矢量场转变为面心标量场，由矢量转变为标量，一定是将该矢量做了点乘，因此不难分析出该语句对应的操作是
f
l
u
x
[
(
H
b
y
A
)
P
]
=
(
H
b
y
A
)
f
⋅
S
\bf{flux[(HbyA)_P]}=(\bf{HbyA})_f\cdot S
flux[(HbyA)P]=(HbyA)f⋅S
注意，上式没有多余的符号（如求和号）。该过程涉及将体心值插值为面心值
fvc::div(phiHbyA)
由上文关于
d
i
v
\bf{div}
div的分析可知，该代码对应的离散过程如下，因此可分析出其返回值类型应为scalarField
f
v
c
:
:
d
i
v
(
p
h
i
H
b
y
A
)
=
∑
(
H
b
y
A
)
f
⋅
S
\bf{fvc::div(phiHbyA)}=\sum (\bf{HbyA})_f\cdot S
fvc::div(phiHbyA)=∑(HbyA)f⋅S
fvm::laplacian(rAU, p)
由上文分析可知，fvm的返回值类型为fvMatrix，具体为标量还是矢量取决于拉普拉斯项中的未知数类型。该代码对应的过程为
∫
V
P
∇
⋅
(
1
A
∇
p
)
d
V
=
∫
∂
V
P
(
1
A
∇
p
)
d
S
≈
∑
(
1
A
∇
p
n
+
1
)
f
S
=
∑
[
(
1
A
∇
p
n
+
1
)
f
⋅
S
S
f
]
S
f
=
a
⋅
p
P
n
+
1
+
∑
b
N
p
N
n
+
1
\int_{V_P}\nabla \cdot \left( \frac{1}{A} \nabla p \right) dV = \int_{\partial V_P} \left( \frac{1}{A} \nabla p \right) dS \approx \sum \left( \frac{1}{A} \nabla p^{n+1} \right)_f S = \sum \left[ \left( \frac{1}{A} \nabla p^{n+1} \right)_f \cdot \frac{S}{S_f} \right]S_f = a\cdot p_P^{n+1} + \sum b_N p_N^{n+1}
∫VP∇⋅(A1∇p)dV=∫∂VP(A1∇p)dS≈∑(A1∇pn+1)fS=∑[(A1∇pn+1)f⋅SfS]Sf=a⋅pPn+1+∑bNpNn+1
上式右数第二项为laplacian项对应的形式，最后一项为返回值
原文地址：https://www.jb51.cc/wenti/3282288.html
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
