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} \]

注意：

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

方程(2)即为控制一维槽道流的控制方程。该方程具有如下数学特征：

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

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

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

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

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

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

(3)\[ \int\left(-\frac{\partial \tau_{xy}}{\partial y}\right)\rd y=-\int\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)做离散有：

(5)\[ A_Pu_P^*+\sum A_Nu_N^*=-gradP^{n} \]

对其求解可以获得预测速度\(u^*\)。再次强调，方程(5)右侧针对所有网格点都是单一值。在收敛的情况下应该有：

(6)\[ A_Pu_P^{n+1}+\sum A_Nu_N^{n+1}=-gradP^{n+1} \]

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

(7)\[ A_Pu_P^{'}+\sum A_Nu_N^{'}=-gradP^{'} \]

在这个方程中，略去临点的影响：

(8)\[ A_Pu_P^{'}=-gradP^{'} \]

现在要获取方程(8)中的修正量的值。在这里，我们有下列已知平均速度关系：

(9)\[ \bar{u}=\bar{u}\frac{\sum V_P}{\sum V_P}=\frac{\sum \bar{u} V_P}{\sum V_P} \]

方程(5)获得预测速度\(u^*\)的体平均值为\(\sum u^*_PV_P/\sum V_P\)，因此我们近似的将速度修正量理解为：

(10)\[ u^{'}\approx \bar{u}-\frac{\sum u_P^* V_P}{\sum V_P} \]

注意这里的\(u^{'}\)为一个单一值而非场。在有了方程(10)之后，可以近似获得修正速度：

(11)\[ u_P^{**}=u_P^{*}+u_P^{'}\approx u_P^{*}+u^{'}=u_P^{*}+\left(\bar{u}-\frac{\sum u_P^* V_P}{\sum V_P}\right) \]

注意其中\(u_P^{'}\)与\(u^{'}\)的区别，前者表示每个网格的修正速度，后者表示平均修正速度。现在要获取压力修正量，如果将方程(10)代入到(8)，会发现方程两边存在问题。例如对于下面的方程：

(12)\[ A_P\left(\bar{u}-\frac{\sum u_P^* V_P}{\sum V_P}\right)=?-gradP^{'} \]

方程左边的\(A_P\)针对每个网格点有不同的值，但方程右侧的\(-gradP^{'}\)却为单一值。因此二者并不能作为一个等式成立。因此，类似的，我们可以把这个方程改写为：

(13)\[ \frac{\sum V_PA_P}{\sum V_P}\left(\bar{u}-\frac{\sum u_P^* V_P}{\sum V_P}\right)=-gradP^{'} \]

因此，修正压力可以表示为：

(14)\[ -gradP^{*}=-gradP^{n}-gradP^{'}=-gradP^{n}+\frac{\sum V_PA_P}{\sum V_P}\left(\bar{u}-\frac{\sum u_P^* V_P}{\sum V_P}\right) \]

至此，我们有了修正压力项\(-gradP^{*}\)，修正速度项\(u_P^{**}\)。由于在此过程中引入了若干假定。因此，需要在修正压力、修正速度的基础上，继续进行迭代。直至每一次的修正压力、修正速度趋向于不变。