理论手册内容请参考理论一与理论二

2. 有限体积法求解二维方腔流–求解

二维模拟区域尺寸为x∈[0,0.1],y∈[0,0.1]x\in [0,0.1], y\in [0,0.1]x∈[0,0.1],y∈[0,0.1]，每一维各有NNN个网格，网格长度为h=1/Nh=1/Nh=1/N，网格体心间距d=hd=hd=h，物理量均定义在网格体心。流体粘度为ν\nuν。待求物理量为速度矢量U=[u,v]U=[u,v]U=[u,v]以及压力标量ppp，上边界xxx方向速度u=1u=1u=1，yyy方向速度v=0v=0v=0，其他三边u=v=0u=v=0u=v=0，四边的压力边界条件均为第二类边界条件即∂p/∂x=∂p/∂y=0\partial p/ \partial x = \partial p/ \partial y =0∂p/∂x=∂p/∂y=0

离散格式：对流项中心差分，扩散项中心差分，压力项中心差分，非稳态项一阶欧拉

模拟区域左下角坐标为(0,0)(0,0)(0,0)，右上角坐标为(1,1)(1,1)(1,1)，网格全局编号从0开始，令网格[i,j][i,j][i,j]表示该网格为xxx方向的第iii个网格、yyy方向的第jjj个网格。对于二维，标量SfS_fSf为表面积，等于hhh

2.1. 代码链接

python代码链接

2.2. 求解流程

给定t0t_0t0时刻网格值，所有网格u=v=p=0u=v=p=0u=v=p=0
根据t0t_0t0时刻的uvuvuv计算系数AP,AW,AE,AS,ANA_P,A_W,A_E,A_S,A_NAP,AW,AE,AS,AN，对应update_coefficient函数
根据t0t_0t0时刻的ppp计算uvuvuv，记为u∗,v∗u^*,v^*u∗,v∗，对应update_velo_GS函数
进入内循环，循环次数nCorrnCorrnCorr
- 由u∗,v∗u^*,v^*u∗,v∗计算HbyA∣∗HbyA\big|^*HbyA∣∣∗，由t_0时刻速度更新APA_PAP，对应update_HbyA_AP函数
- 组建压力泊松方程，由HbyA∣∗HbyA\big|^*HbyA∣∣∗计算ppp，记为p∗p^*p∗，对应update_pressure_GS函数
- 由p∗p^*p∗更新速度，对应update_velo_final函数
- 计算连续性方程误差，对应update_continuity_error函数。开始下一次内循环
时间推进，重复以上步骤

可以看出，以上步骤对应PISO算法

2.3. 注意事项

由t0t_0t0时刻推进到t1t_1t1时刻的过程中，假设求出了中间速度U∗U^*U∗
- 系数AP,ANA_P,A_NAP,AN均由t0t_0t0时刻的速度求解
- HbyA由最新速度U∗U^*U∗更新
压力泊松方程中，不要忘记VPV_PVP，推导过程中很容易忽视VP∇p=∑pfSV_P \nabla p = \sum p_f SVP∇p=∑pfS，有了VPV_PVP该等式两端的量纲才相等
求解完动量方程后可以计算连续性方程误差err_mom，此时误差应该较大，因为还没有对压力进行修正。在完成压力修正之后，连续性误差应小于err_mom，可以通过这个指标判断程序是否执行正确

2.4. 与icoFoam计算结果对比

两方向网格数N=20N=20N=20，粘度ν=0.01m2/s\nu=0.01 m^2/sν=0.01m2/s，对比物理时间t=0.1st=0.1st=0.1s时的流场演化结果

2.4.1. x方向速度

this work
icoFoam

2.4.2. y方向速度

this work
icoFoam

2.4.3. 压力

this work
icoFoam

2.4.4. x方向第10个网格位置处，u轴向分布

横坐标为位置

2.4.5. 最顶层网格的u分布

横坐标为位置

3. 附录

3.1. 附录–张量运算

并矢量：将矢量张成二阶矢量即矩阵

UU=(uv)(uv)=(uuuvvuvv)UU = \left( \begin{matrix} u \\ v \end{matrix} \right) \left( \begin{matrix} u\ v \end{matrix} \right) = \left( \begin{matrix} uu & uv \\ vu & vv \end{matrix} \right) UU=(uv)(u v)=(uuvuuvvv)

∇U=(∂x∂y)(uv)=(∂u/∂x∂v/∂x∂u/∂y∂v/∂y)\nabla U = \left( \begin{matrix} \partial_x \\ \partial_y \end{matrix} \right) \left( \begin{matrix} u\ v \end{matrix} \right) = \left( \begin{matrix} {\partial u}/{\partial x} & {\partial v}/{\partial x} \\ {\partial u}/{\partial y} & {\partial v}/{\partial y} \end{matrix} \right) ∇U=(∂x∂y)(u v)=(∂u/∂x∂u/∂y∂v/∂x∂v/∂y)

∇⋅(UU)=(∂/∂x∂/∂y)(uuuvvuvv)=(∂uu∂x+∂vu∂y∂uv∂x+∂vv∂y)\nabla \cdot (UU) = \left( \begin{matrix} {\partial}/{\partial x} & {\partial}/{\partial y} \end{matrix} \right) \left( \begin{matrix} uu & uv \\ vu & vv \end{matrix} \right) = \left( \frac{\partial uu}{\partial x}+\frac{\partial vu}{\partial y} \quad \frac{\partial uv}{\partial x}+\frac{\partial vv}{\partial y} \right) ∇⋅(UU)=(∂/∂x∂/∂y)(uuvuuvvv)=(∂x∂uu+∂y∂vu∂x∂uv+∂y∂vv)

∇⋅(∇U)=(∂/∂x∂/∂y)(∂u/∂x∂v/∂x∂u/∂y∂v/∂y)=(∂2u∂x2+∂2u∂y2∂2v∂x2+∂2v∂y2)\nabla \cdot (\nabla U) = \left( \begin{matrix} {\partial}/{\partial x} & {\partial}/{\partial y} \end{matrix} \right) \left( \begin{matrix} {\partial u}/{\partial x} & {\partial v}/{\partial x} \\ {\partial u}/{\partial y} & {\partial v}/{\partial y} \end{matrix} \right) = \left( \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \quad \frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2} \right) ∇⋅(∇U)=(∂/∂x∂/∂y)(∂u/∂x∂u/∂y∂v/∂x∂v/∂y)=(∂x2∂2u+∂y2∂2u∂x2∂2v+∂y2∂2v)

3.2. 附录–关于扩散项和压力的讨论

此节内容摘自《传递过程原理》。三维牛顿粘度定律的张量形式为
τ=μ[(∇u)+(∇u)T]−(23μ−κ)(∇⋅u)δ\tau=\mu \left[ (\nabla u) + (\nabla u)^T \right] - \left( \frac{2}{3}\mu - \kappa \right) (\nabla \cdot u)\delta τ=μ[(∇u)+(∇u)T]−(32μ−κ)(∇⋅u)δ

τ\tauτ为应力张量，μ\muμ为第一粘性系数，κ\kappaκ为第二粘性系数，δ\deltaδ为单位张量。上式中等号右边的第一项为变形速度张量SSS
S=12[(∇u)+(∇u)T]S=\frac{1}{2}\left[ (\nabla u) + (\nabla u)^T \right] S=21[(∇u)+(∇u)T]

该张量主对角线上的元素又称为相对拉伸（压缩）速率，非主对角线元素称为角变形速率，详细可以查看流体力学教材。

对于一般情形下的可压缩流体，κ=0\kappa=0κ=0
τ=μ[(∇u)+(∇u)T]−23μ(∇⋅u)δ\tau=\mu \left[ (\nabla u) + (\nabla u)^T \right] - \frac{2}{3}\mu (\nabla \cdot u)\delta τ=μ[(∇u)+(∇u)T]−32μ(∇⋅u)δ

对于不可压缩流体，∇⋅u=0\nabla \cdot u=0∇⋅u=0
τ=μ[(∇u)+(∇u)T]=2μS\tau=\mu \left[ (\nabla u) + (\nabla u)^T \right]=2\mu S τ=μ[(∇u)+(∇u)T]=2μS

上式与动量方程公式(2)(2)(2)对比可发现，τ(2)=ν∇u\tau(2)=\nu \nabla uτ(2)=ν∇u，与上式推导不同。以二维为例，
∇⋅[(∇u)+(∇u)T]=(∂/∂x∂/∂y)(2∂u/∂x∂v/∂x+∂u/∂y∂v/∂x+∂u/∂y2∂v/∂y)=(2∂2u∂x2+∂2v∂x∂y+∂2u∂y2∂2v∂x2+∂2u∂x∂y+2∂2v∂y2)\nabla \cdot \left[ (\nabla u) + (\nabla u)^T \right] = \left( \begin{matrix} {\partial}/{\partial x} & {\partial}/{\partial y} \end{matrix} \right) \left( \begin{matrix} 2{\partial u}/{\partial x} & {\partial v}/{\partial x}+{\partial u}/{\partial y} \\ {\partial v}/{\partial x}+{\partial u}/{\partial y} & 2{\partial v}/{\partial y} \end{matrix} \right) \\ = \left( 2\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 v}{\partial x \partial y}+\frac{\partial^2 u}{\partial y^2} \quad \frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 u}{\partial x \partial y}+2\frac{\partial^2 v}{\partial y^2} \right) ∇⋅[(∇u)+(∇u)T]=(∂/∂x∂/∂y)(2∂u/∂x∂v/∂x+∂u/∂y∂v/∂x+∂u/∂y2∂v/∂y)=(2∂x2∂2u+∂x∂y∂2v+∂y2∂2u∂x2∂2v+∂x∂y∂2u+2∂y2∂2v)

取出第一项
(∇⋅τ)x=2∂2u∂x2+∂2v∂x∂y+∂2u∂y2=∂2u∂x2+∂2u∂y2+∂∂x(∂u∂x+∂v∂y)=∇2u+0(\nabla \cdot \tau)_x = 2\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 v}{\partial x \partial y}+\frac{\partial^2 u}{\partial y^2} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial}{\partial x}\left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) = \nabla^2 u + 0 (∇⋅τ)x=2∂x2∂2u+∂x∂y∂2v+∂y2∂2u=∂x2∂2u+∂y2∂2u+∂x∂(∂x∂u+∂y∂v)=∇2u+0

因此∇⋅τ=∇2U\nabla \cdot \tau=\nabla^2 U∇⋅τ=∇2U，公式(2)(2)(2)的写法无误，即NS方程中的扩散项为∇⋅τ\nabla \cdot \tau∇⋅τ

三维牛顿粘度定律从本质上讲，是描述应力张量和变形速度张量之间的流体粘性规律，如果将流体压力和粘性力合起来考虑，引进张量σ\sigmaσ
σ=−pδ+τ=−pδ+2μ(S−13(∇⋅u)δ)\sigma = -p\delta + \tau = -p\delta + 2\mu \left( S-\frac{1}{3}(\nabla \cdot u)\delta \right) σ=−pδ+τ=−pδ+2μ(S−31(∇⋅u)δ)

3.2.1. 湍流模型中对应力的处理

在直接数值模拟（DNS）中，无需引入湍流模型，求解的扩散项就是上述方程。然而在粗网格模拟中，需要引入湍流模型。

对于RANS类模型，常用的k−ϵk-\epsilonk−ϵ或者k−ωk-\omegak−ω方程中的Reynolds stress被表示为
−ρuu‾=μt(∇u‾+(∇u‾)T)−23ρkδ-\rho \overline{uu}=\mu_t (\nabla \overline{u} + (\nabla \overline{u})^T)-\frac{2}{3}\rho k \delta −ρuu=μt(∇u+(∇u)T)−32ρkδ

其中u‾\overline{u}u为系综平均的速度，kkk为湍动能。则扩散项为
∇⋅(μ(∇u‾+(∇u‾)T)−23μ(∇⋅u‾)δ−ρuu‾)=∇⋅((μ+μt)(∇u‾+(∇u‾)T)−23μ(∇⋅u‾)δ–23ρkδ)=∇⋅(μeff(∇u‾+(∇u‾)T)−23μ(∇⋅u‾)δ)–∇⋅(23ρkδ)\nabla \cdot \left( \mu \left( \nabla \overline{\boldsymbol{u}} + (\nabla \overline{\boldsymbol{u}})^T \right) -\frac{2}{3}\mu (\nabla \cdot \overline{\boldsymbol{u}}) \delta -\rho\overline{\boldsymbol{u}\boldsymbol{u}} \right) \\ = \nabla \cdot \left( (\mu + \mu_t) \left( \nabla \overline{\boldsymbol{u}} + (\nabla \overline{\boldsymbol{u}})^T \right) -\frac{2}{3}\mu (\nabla \cdot \overline{\boldsymbol{u}}) \delta – \frac{2}{3}\rho k \delta \right) \\ = \nabla \cdot \left( \mu_{eff} \left( \nabla \overline{\boldsymbol{u}} + (\nabla \overline{\boldsymbol{u}})^T \right) -\frac{2}{3}\mu (\nabla \cdot \overline{\boldsymbol{u}}) \delta \right) – \nabla \cdot \left( \frac{2}{3}\rho k \delta \right) ∇⋅(μ(∇u+(∇u)T)−32μ(∇⋅u)δ−ρuu)=∇⋅((μ+μt)(∇u+(∇u)T)−32μ(∇⋅u)δ–32ρkδ)=∇⋅(μeff(∇u+(∇u)T)−32μ(∇⋅u)δ)–∇⋅(32ρkδ)

其中μeff\mu_{eff}μeff为层流粘性和湍流粘性之和。可压缩和不可压缩流体均适用于上式。

3.2.1.1. OpenFOAM中的处理

从前文的分析可知，控制方程中体现湍流模型的地方仅有扩散项，因此OpenFOAM中也是在扩散项中囊括了湍流模型，比如在rhoPimpleFoam求解器的tUEqn中有这样一行

turbulence->divDevRhoReff(U)

若选择linear eddy viscosity model线性涡粘模型，则上述函数定义为

Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(volVectorField& U
) const
{return(- fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))- fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U));
}

OpenFOAM中grad(u)=∇ugrad(u)=\nabla ugrad(u)=∇u，为张量，另外
dev(∇u)=(∇u)−13tr(∇u)dev(\nabla u)=(\nabla u)-\frac{1}{3}tr(\nabla u) dev(∇u)=(∇u)−31tr(∇u)

dev2(∇u)=(∇u)−23tr(∇u)dev2(\nabla u)=(\nabla u)-\frac{2}{3}tr(\nabla u) dev2(∇u)=(∇u)−32tr(∇u)
其中tr(∇u)=∇⋅u=tr((∇u)T)tr(\nabla u)=\nabla\cdot u=tr((\nabla u)^T)tr(∇u)=∇⋅u=tr((∇u)T)，即矩阵∇u\nabla u∇u的对角线。因此返回值第一行对应的是
∇⋅(μeff((∇u)T−23(∇⋅u)δ))\nabla \cdot \left( \mu_{eff} \left((\nabla u)^T-\frac{2}{3}(\nabla\cdot u)\delta \right) \right) ∇⋅(μeff((∇u)T−32(∇⋅u)δ))

返回值第二行对应的是
∇⋅(μeff∇u)\nabla \cdot (\mu_{eff} \nabla u) ∇⋅(μeff∇u)

除此之外，还缺少最后一项即
∇⋅(23ρkδ)\nabla \cdot \left( \frac{2}{3}\rho k \delta \right) ∇⋅(32ρkδ)
网上给出的解释是该项被纳入到压力梯度项中进行计算，详情参考

3.2.2. 压力

流体力学中所研究的压力指的是static pressure。不可压缩的伯努利方程为
P+12ρv2=P0P+\frac{1}{2}\rho v^2 = P_0 P+21ρv2=P0

上式中的PPP为static pressure，1/2ρv21/2\rho v^21/2ρv2为dynamic pressure，P0P_0P0为total pressure或者stagnation pressure，即驻点压力，也就是流体冲击到一个静止物体时，该物体所受的流体的冲击力。另外还有一种压力为hydrostatic pressure即静水压，该压力在interFoam求解器中有所解析，可参考这里

更多文献参考wiki

有限体积法求解二维方腔流（三）——代码以及与icoFoam结果对比相关推荐

有限体积法求解二维方腔流（一）——动量方程和连续性方程的离散
1. 有限体积法求解二维方腔流–理论手册 1.1. 不可压缩流体控制方程连续性方程 ∇⋅U=0(1)\nabla \cdot U=0 \tag{1} ∇⋅U=0(1) 动量方程 ∂U∂t+∇⋅(UU ...
有限体积法求解二维方腔流（二）——边界条件处理及SIMPLE和PISO算法
关于动量方程和连续性方程的离散参考这里 1.5. 边界处理 1.5.1. 左边界,uw=vw=0u_w=v_w=0uw=vw=0,不包括边角对流项 ∑FfnUfn+1=(−uwnh)Uwn+1+ ...
【运筹优化】蚁群算法求解二维矩形装箱问题（java代码实现）
文章目录 1 前言 2 代码迁移 3 蚁群算法 3.1 蚂蚁类 Ant 3.2 蚁群算法类 ACO_Packing 4 运行结果 5 后话 [运筹优化]求解二维矩形装箱问题的算法合辑(Java代码实现 ...
二维方格子Ising模型代码
#include <iostream> #include <cstdlib> #include <cmath> #include <fstream> # ...
求解欧拉方程的c语言,用有限体积方法求解欧拉方程
<用有限体积方法求解欧拉方程>由会员分享,可在线阅读,更多相关<用有限体积方法求解欧拉方程(12页珍藏版)>请在人人文库网上搜索. 1.有限体积法求解二维可压缩Euler方程计 ...
求解欧拉方程的c语言,用有限体积方法求解欧拉方程.doc
PAGE 4 - 实用标准文档有限体积法求解二维可压缩Euler方程 --计算流体力学课程大作业老师: 夏健.刘学强学生: 徐锡虎学号: SQ09018013018 日期: 2010年2月5 ...
有限体积法及其网格简介
有限体积法及其网格简介有限体积法使目前CFD领域广泛使用的离散化方法,其特点不仅表现在对控制方程的离散结果上,还表现在所使用的网格上. 1.有限体积法的基本思想有限体积法(Finite Volum ...
二阶常微分方程的数值解法（中心差分法和有限体积法）
二阶常微分方程的数值解法(中心差分法和有限体积法) 这里我们介绍中心差分法和有限体积法求解方程. 题目: 用差分法的中心差分格式和有限体积法求解两点边值问题 u′′−α(2x−1)u′−2αu=0,0 ...
【运筹优化】求解二维矩形装箱问题的算法合辑 + Java代码实现
文章目录一.算法合辑 1.1 结合自定义策略 1.1.1 结合自定义策略的禁忌搜索算法(98.80%) 1.2.1 结合自定义策略的蚁群算法(96.70%) 1.2 堆优化的天际线启发式算法(98. ...
二维有限体积 matlab,二维有限体积法计算热传导及源码.pdf
二维有限体积法计算热传导及源码 //#include "stdafx.h" #include #include #include #include #include using n ...

有限体积法求解二维方腔流（三）——代码以及与icoFoam结果对比