版本对应：

该算法对应所有OpenFOAM中的boundaryFoam。

CFD: 一维充分发展槽道流算法

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

需要注意的是，在真实的三维槽道流中，不仅存在\(x\)方向的速度，也存在\(y,z\)方向的速度，且其时间平均后依然存在。如果需要真实的三维槽道流的径向分布数据（存在3个分量的速度以及压力分布），需要进行纯三维的计算。

控制方程与算法

首先有连续性方程

(1)\[ \nabla\cdot\overline{\bfU}=0 \]

其中\(\overline{\bfU}\)表示时间平均后的速度。首先假定不存在\(z\)方向，上式演变为：

(2)\[ \frac{\p \overline{u}}{\p x}+\frac{\p \overline{v}}{\p y}=0 \]

在充分发展的情况下，\(x\)流向方向的所有变量梯度为0，因此上式演变为

(3)\[ \frac{\p \overline{v}}{\p y}=0 \]

注意到\(\overline{v}\)在固体壁面的时候为0，因此方程(13)可以演变为

(4)\[ \overline{v}=0 \]

因此，充分发展的槽道流的时间平均后的连续性方程，就是\(\overline{v}=0\)。

现考虑动量方程：

(5)\[ \frac{\p\overline{\bfU}}{\p t}+\nabla \cdot (\overline{\bfU}\overline{\bfU})-\nabla \cdot(\nu \nabla \overline{\bfU})=-\nabla \overline{p}-\nabla\cdot\boldsymbol{\tau_t}, \]

同样的，假定不存在\(z\)方向，因此\(x\)方向的动量方程可以写为：

(6)\[ \frac{\partial \overline{u}}{\partial t} + \frac{\partial \overline{u} \overline{u}}{\partial x} + \frac{\partial \overline{v} \overline{u}}{\partial y} - \frac{\partial}{\partial x}\left(\nu\frac{\partial \overline{u}}{\partial x}\right) - \frac{\partial}{\partial y}\left(\nu\frac{\partial \overline{u}}{\partial y}\right) = -\frac{\partial \overline{p}}{ \partial x} -\frac{\p \overline{u'u'}}{\p x} -\frac{\p \overline{v'u'}}{\p y} \]

在稳态情况下不存在时间项，同时，\(x\)流向方向的所有变量梯度为0（除了压力）以及\(\overline{v}=0\)，方程(6)可以写为

(7)\[ - \frac{\partial}{\partial y}\left(\nu\frac{\partial \overline{u}}{\partial y}\right) = -\frac{\partial \overline{p}}{ \partial x} -\frac{\p \overline{v'u'}}{\p y} \]

同理，\(y\)方向的动量方程可以写为：

(8)\[ \frac{\partial \overline{v}}{\partial t} + \frac{\partial \overline{u} \overline{v}}{\partial x} + \frac{\partial \overline{v} \overline{v}}{\partial y} - \frac{\partial}{\partial x}\left(\nu\frac{\partial \overline{v}}{\partial x}\right) - \frac{\partial}{\partial y}\left(\nu\frac{\partial \overline{v}}{\partial y}\right) = -\frac{\partial \overline{p}}{ \partial y} -\frac{\p \overline{u'v'}}{\p x} -\frac{\p \overline{v'v'}}{\p y} \]

经过简化后有

(9)\[ 0 = -\frac{\partial \overline{p}}{ \partial y} -\frac{\p \overline{v'v'}}{\p y} \]

方程(9)、(7)、(4)表示在一个槽道中，充分发展的忽略\(z\)方向的时间平均后的NS方程。也就是说对于一个直接模拟时间平均后的数据，应该满足方程(9)、(7)、(4)。其具有以下特点：

\(\overline{v}=0\)；
\(\overline{p}\)在\(y\)方向上存在梯度；

Warning

\(\overline{p}\)在\(y\)方向上存在梯度，但是在壁面处\(\overline{p}\)在\(y\)方向的梯度为0，因为壁面处脉动速度为0，否则不满足壁面不可穿透的物理条件。

Warning

为什么方程方程(9)中\(\overline{v}=0\)，但是\(\overline{v'v'}\neq 0\)？

可以考虑3个时间步下的序列：\(v_{t1}=1,v_{t1}=2,v_{t1}=3\)。时间平均下的\(\overline{v}=2\)。因此\(v'_{t1}=-1,v_{t1}=0,v_{t1}=1\)。同时有\(\overline{v'v'}=(1+1)/3\neq 2\)。也即\(\overline{v'v'}\neq \overline{v}\)。

方程(9)、(7)、(4)可以进一步进行简化。如果网格仅仅具有\(y\)方向上的维度，且给定流向周期性边界的情况下，自动满足方程(4)。同时，\(\frac{\partial \overline{p}}{ \partial x} \)可以认为是一个标量。

如果进一步调用RANS湍流模型，方程(7)中的雷诺应力\(-\frac{\p \overline{v'u'}}{\p y}\)可以被模化，并写为：

(10)\[ - \frac{\partial}{\partial y}\left(\nu\frac{\partial \overline{u}}{\partial y}\right) = -\frac{\partial \overline{p}}{ \partial x} + \frac{\partial}{\partial y}\left(\nu_t\frac{\partial \overline{u}}{\partial y}\right) \]

即：

(11)\[ - \frac{\partial}{\partial y}\left(\nu_{eff}\frac{\partial \overline{u}}{\partial y}\right) = -\frac{\partial \overline{p}}{ \partial x} \]

很明显，如果给定\(\frac{\partial \overline{p}}{ \partial x} \)的值。方程(11)即为一个泊松方程，可以求解并获得\(\overline{u}\)的值。方程(11)具有如下数学特征：

其类似壁面温度为零的固体导热方程。压力梯度用于提供一个平均的温度，最终温度分布呈现抛物线结构。这和槽道流中的速度型线是一致的。
若压力梯度为零，方程(11)的解为零，这也解释了为什么周期的槽道流为什么需要提供压力梯度驱动力。

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

压力梯度法：若已知压力梯度，则可以直接代入压力梯度值，求解\(\frac{\partial}{\partial y}\left(\nu_t\frac{\partial \overline{u}}{\partial y}\right)\)即可；
速度修正法：若已知平均流速，可以通过附加流速限定性条件，对方程(11)进行封闭；

通常情况下，平均流速比较容易获取，压力梯度的值很难得到。但是对于充分发展的湍流，可以获得压力梯度与壁面剪切力的平衡关系。对方程(11)做\(y\)方向的积分（\(\tau_{xy}=\nu_{eff}\frac{\partial \overline{u}}{\partial y}\)）：

(12)\[ \int_0^h\left(-\frac{\partial \tau_{xy}}{\partial y}\right)\rd y=-\int_0^h\left(\frac{\partial \overline{p}}{ \partial x}\right)\rd y \]

(13)\[ \tau_{xy}^h -\tau_{w}=-\tau_{w}=-\frac{\partial \overline{p}}{ \partial x}h \]

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

Warning

很明显，方程(11)在1维情况下可以直接求解，并不需要方程(9)。确实如此，boundaryFoam并没有对方程(9)进行求解，因此压力\(\overline{p}\)在\(y\)方向上不具有分布。这与DNS时间平均后的数据是不一样的。

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

给定一个猜测的压力梯度，求解传输速度。
计算传输速度和预期速度的差值。
更新压力梯度修正量，速度修正量，并返回第一步重新计算传输速度，使得结果一步步的逼近真实值。

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

(14)\[ -\frac{\partial}{\partial y}\left(\nu_{eff}\frac{\partial \overline{u}}{\partial y}\right)=-\mathrm{gP} \]

首先，对方程(14)做离散有：

(15)\[ A_{\rP}\overline{u}_{\rP}^*+\sum A_{\rN}\overline{u}_{\rN}^*=-\frac{1}{\Delta V_{\rP}}\int \mathrm{gP} ^{n}\rd V=-\mathrm{gP}^{n} \]

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

(16)\[ A_{\rP}\overline{u}_{\rP}^{n+1}+\sum A_{\rN}\overline{u}_{\rN}^{n+1}=-\mathrm{gP}^{n+1} \]

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

(17)\[ A_{\rP}\overline{u}_{\rP}^{'}+\sum A_{\rN}\overline{u}_{\rN}^{'}=-\mathrm{gP}^{'} \]

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

(18)\[ A_{\rP}\overline{u}_{\rP}^{'}=-\mathrm{gP}^{'} \]

看起来像SIMPLE算法？

如果继续回代到方程(15)中有：

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

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

方程(18)有下述特点：

右侧\(-\mathrm{gP}^{'} \)为一个单一值；左侧\(A_\rP,\overline{u}_{\rP}^{'}\)则定义在每个网格点上；
若每个网格的\(A_\rP\)不相同，则\(\overline{u}_{\rP}^{'}\)必然不同；

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

(20)\[ \bar{u}_{mean}=\frac{\sum \overline{u}_\rP V_\rP}{\sum V_\rP} \]

其中\(\bar{u}_{mean}\)为用户给定的值。另一方面，方程(15)获得预测速度\(\overline{u}^*\)的体平均值为

(21)\[ \overline{u}^{*}=\frac{\sum \overline{u}^*_PV_P}{\sum V_P} \]

因此我们近似的将速度修正量理解为：

(22)\[ \overline{u}'\approx \bar{u}_{mean}-\overline{u}^{*}=\bar{u}_{mean}-\frac{\sum \overline{u}^*_PV_P}{\sum V_P} \]

注意这里的\(u^{'}\)为一个单一值而非场。

在这里继续引入假设：每个网格点的速度修正\(\overline{u}_{\rP}'\)均是相同的，即\(\overline{u}'\)。将方程(22)带入到方程(18)左侧有：

(23)\[ A_{\rP}\left(\bar{u}_{mean}-\frac{\sum \overline{u}_\rP^* V_\rP}{\sum V_\rP}\right) \]

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

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

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

(25)\[ -\mathrm{gP}^{*}=-\mathrm{gP}^{n}-\mathrm{gP}^{'}= -\mathrm{gP}^{n}+\frac{\sum V_\rP A_{\rP}}{\sum V_\rP} \left(\bar{u}_{mean}-\frac{\sum \overline{u}_\rP^* V_\rP}{\sum V_\rP}\right) \]

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

(26)\[ \overline{u}_{\rP}^{**}=\overline{u}_{\rP}^{*}+\overline{u}' \]

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

关键代码

boundaryFoam中的速度方程为：

fvVectorMatrix UEqn
(
    divR == gradP
);

其对应方程(15)。

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

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

U += (Ubar - UbarStar);

上述代码即方程(26)。

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

上述代码即方程(25)。

算例在运行后，还会输出相关变量，比如\(\mathbf{R}\):

(27)\[ \mathbf{R}=\frac{2}{3}k\bfI -\nu_t\left(\nabla\bfU+\nabla\bfU^T-\frac{2}{3}( \nabla\cdot\bfU )\bfI\right) \]

在boundaryFoam做简化后，\(\mathbf{R}\)可以写为：

(28)\[\begin{split} \mathbf{R}= \begin{bmatrix} \frac{2}{3}k & -\nu_t\frac{\p \overline{u}}{\p y} & 0\\ -\nu_t\frac{\p \overline{u}}{\p y} & \frac{2}{3}k & 0\\ 0 & 0 & \frac{2}{3}k \end{bmatrix} \end{split}\]

因此，针对boundaryFoam预测的模化后雷诺应力有：

(29)\[ \overline{u'u'}=\overline{v'v'}=\overline{w'w'}=\frac{2}{3}k, \overline{u'v'}=-\nu_t\frac{\p \overline{u}}{\p y} \]