boundaryFoam解析

引言

大量的LES与DNS被应用于求解槽道流(案例可参考icoFoam解析)。在模拟槽道流的过程中,要获得充分发展的湍流,通常需要此管道足够的长,并且模拟足够长的周期。这样会大大增加网格数量并增加计算时间。另一方面,充分发展的槽道流,可用近似一维的方程进行描述。boundaryFoam既为这样一个求解器,可用于求解充分发展的槽道湍流。其主要用于求解充分发展的槽道流流中的径向湍流变量分布。在boundaryFoam中可调用不同的RANS湍流模型,因此求解的湍流场可用于为其他算例提供一个精准的初始场或进口边界条件。同时,对于周期发展的槽道流,为了保证所需的流量不被粘性耗散,或需要添加固定的压力梯度驱动,或提供固定的流量参数。boundaryFoam采用第二种方式保证流速。

控制方程与算法

现考虑\(x\)方向的附加湍流粘度的动量方程:

(1)\[\begin{split} \frac{\partial u_1}{\partial t} + \frac{\partial u_1 u_1}{\partial x}+\frac{\partial u_2 u_1}{\partial y}+\frac{\partial u_3 u_1}{\partial z} \hspace{15em} \\ = -\frac{\partial p}{ \partial x} + \frac{\partial}{\partial x}\left(\nu_t\frac{\partial u_1}{\partial x}\right)+\frac{\partial}{\partial y}\left(\nu_t\frac{\partial u_1}{\partial y}\right)+\frac{\partial}{\partial z}\left(\nu_t\frac{\partial u_1}{\partial z}\right) \end{split}\]

其中\(u_1\)\(u_2\)\(u_3\)为速度的分量,\(p\)为除以密度之后的压力,\(\nu_t\)为有效粘度。在对其进行简化的时候,下面的情况是适用的:

  • 稳态情况,方程(1)中的第一项为\(0\)

  • 一维问题,\(u_2=u_3=0\)同时\(\frac{\p}{\p z}=0\)

  • 对于充分发展的流动,对于除了压力外的变量有\(\frac{\p}{\p x}=0\)

因此,方程(1)可以化简为:

(2)\[ -\frac{\partial}{\partial y}\left(\nu_t\frac{\partial u_1}{\partial y}\right)=-\frac{\partial p}{ \partial x} \]

方程(2)即为控制一维槽道流的控制方程。

连续性方程在哪里?

一维槽道流的连续性方程隐含在\(\frac{\p}{\p x}=0\)。若施加周期性边界条件,\(x\)方向一个网格,因此必定满足连续性方程。

同时,因为\(x\)方向只有一层网格,其中的\(-\frac{\partial p}{ \partial x}\)为一个单一值。

该方程具有如下数学特征:

  • 其类似壁面温度为零的固体导热方程。压力梯度用于提供一个平均的温度,最终温度分布呈现抛物线结构。这和槽道流中的速度型线是一致的。

  • 若压力梯度为零,方程(2)的解为零,这也解释了为什么周期的槽道流为什么需要提供压力梯度驱动力。

方程(2)中有两个未知数\(u_1\)\(p\)并不封闭。在这里可以用两种方法对其封闭:

  • 压力梯度法:若已知压力梯度,则可以直接代入压力梯度值,求解\(\frac{\partial}{\partial y}\left(\nu_t\frac{\partial u_1}{\partial y}\right)\)即可;

  • 速度修正法:若已知平均流速,可以通过附加流速限定性条件,对方程(2)进行封闭;

通常情况下,平均流速比较容易获取,压力梯度的值很难得到。但是对于充分发展的湍流,可以获得压力梯度与壁面剪切力的平衡关系。对方程(2)\(y\)方向的积分:

(3)\[ \int_0^h\left(-\frac{\partial \tau_{xy}}{\partial y}\right)\rd y=-\int_0^h\left(\frac{\partial p}{ \partial x}\right)\rd y \]
(4)\[ \tau_{xy}^h -\tau_{w}=-\tau_{w}=-\frac{\partial p}{ \partial x}h \]

因此,如果已知壁面剪切力,就可以获得压力梯度。壁面剪切力参数在某些直接模拟、大涡模拟的工作中为已知量。在这种情况下,可以给定压力梯度来进行求解。

boundaryFoam采用的是速度修正法。类似的方法也被植入到OpenFOAM中的meanVelocityForce源项中。这是一种更加普适性的做法。速度修正法的主要思想是

  1. 给定一个猜测的压力梯度,求解传输速度。

  2. 计算传输速度和预期速度的差值。

  3. 更新压力梯度修正量,速度修正量,并返回第一步重新计算传输速度,使得结果一步步的逼近真实值。

由于方程(2)右侧为一个单一值,因此现在将方程(2)改写为

(5)\[ -\frac{\partial}{\partial y}\left(\nu_t\frac{\partial u_1}{\partial y}\right)=-\mathrm{gP} \]

首先,对方程(5)做离散有:

(6)\[ A_{\rP}u_{\rP}^*+\sum A_{\rN}u_{\rN}^*=-\frac{1}{\Delta V_{\rP}}\int \mathrm{gP} ^{n}\rd V=-\mathrm{gP}^{n} \]

其中方程右侧的压力梯度取已知量。对其求解可以获得预测速度\(u^*\)。再次强调,方程(6)右侧\(\mathrm{gP}^{n} \)针对所有网格点都是单一值。在收敛的情况下应该有:

(7)\[ A_{\rP}u_{\rP}^{n+1}+\sum A_{\rN}u_{\rN}^{n+1}=-\mathrm{gP}^{n+1} \]

方程(6)与方程(7)相减有修正速度以及修正压力的关系式:

(8)\[ A_{\rP}u_{\rP}^{'}+\sum A_{\rN}u_{\rN}^{'}=-\mathrm{gP}^{'} \]

在这个方程中,略去临点的影响:

(9)\[ A_{\rP}u_{\rP}^{'}=-\mathrm{gP}^{'} \]

看起来像SIMPLE算法?

如果继续回代到方程(6)中有:

(10)\[ A_{\rP}u_{\rP}^{**}=-\mathrm{gP}^{*} \]

后续可以用来组建压力泊松方程。但是因为一维槽道流无法调用连续性方程。因此不存在压力泊松方程,需要使用另外一种策略。

方程(9)有下述特点:

  • 右侧\(-\mathrm{gP}^{'} \)为一个单一值;左侧\(A_\rP,u_{\rP}^{'}\)则定义在每个网格点上;

  • 若每个网格的\(A_\rP\)不相同,则\(u_{\rP}^{'}\)必然不同;

现在要获取方程(9)\(-\mathrm{gP}^{'} \)修正量的值。因为平均速度是用户给定的,因此可以首先通过获取速度修正量的值,然后通过方程(9)中来获取压力修正。首先定义下面体平均速度:

(11)\[ \bar{u}=\frac{\sum \bar{u_\rP} V_\rP}{\sum V_\rP} \]

其中\(\bar{u}\)为用户给定的值。另一方面,方程(6)获得预测速度\(u^*\)的体平均值为\(\bar{u^*}=\sum u^*_PV_P/\sum V_P\),因此我们近似的将速度修正量理解为:

(12)\[ u^{'}\approx \bar{u}-\bar{u^*}=\bar{u}-\frac{\sum u_\rP^* V_\rP}{\sum V_\rP} \]

注意这里的\(u^{'}\)为一个单一值而非场。在这里继续引入假设:每个网格点的速度修正\(u_{\rP}^{'}\)均是相同的,即\(u^{'}\)。将方程(12)带入到方程(9)左侧有:

(13)\[ A_{\rP}\left(\bar{u}-\frac{\sum u_\rP^* V_\rP}{\sum V_\rP}\right) \]

上述方程若成立,\(A_{\rP}\)只能为一个值而不能为一个场(参考前文讨论:若每个网格的\(A_\rP\)不相同,则\(u_{\rP}^{'}\)必然不同。若\(u_{\rP}^{'}\)相同,则\(A_\rP\)需要相同)。因此可以将\(A_{\rP}\)体平均化,且有方程:

(14)\[ \frac{\sum V_\rP A_{\rP}}{\sum V_\rP} \left(\bar{u}-\frac{\sum u_\rP^* V_\rP}{\sum V_\rP}\right)=-\mathrm{gP}^{'} \]

方程(14)可以用来求解\(\mathrm{gP}^{'}\)。同时有更新后的\(\mathrm{gP}^{*}\)为:

(15)\[ -\mathrm{gP}^{*}=-\mathrm{gP}^{n}-\mathrm{gP}^{'}=-\mathrm{gP}^{n}+\frac{\sum V_\rP A_{\rP}}{\sum V_\rP} \left(\bar{u}-\frac{\sum u_\rP^* V_\rP}{\sum V_\rP}\right) \]

且速度可以通过下式更新:

(16)\[ u_{\rP}^{**}=u_{\rP}^{*}+u^{'} \]

随后可以继续进行下一次迭代。

关键代码

boundaryFoam中的速度方程为:

fvVectorMatrix UEqn
(
    divR == gradP
);

其对应方程(6)

dimensionedVector UbarStar = flowMask & U.weightedAverage(mesh.V());

上述代码对应符号\(\bar{u^*}\)

U += (Ubar - UbarStar);

上述代码即方程(16)

gradP += (Ubar - UbarStar)/(1.0/UEqn.A())().weightedAverage(mesh.V());

上述代码即方程(15)