return 0;
wmake

boundaryFoam解析

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

boundaryFoam为一个充分发展一维槽道湍流求解器。主要用于模拟充分发展的槽道流流中的径向湍流变量分布。如果用户要模拟一个充分长的管道,要获得充分发展的湍流,则需要此管道足够的长,这样会大大增加网格数量并大大增加计算时间。但通过boundaryFoam,仅仅需要几十个网格就可以完成求解。本文分析boundaryFoam中的方程和代码,很遗憾目前没有发现使用boundaryFoam进行工程应用案例和学术论文。

如果打不开图像,请右键在新标签页打开图像后刷新几次

图1. 充分发展一维槽道湍流示意图。

方程离散以及求解

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

\begin{equation}\label{mom} \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} = -\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{equation}

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

因此,方程\eqref{mom}可以化简为:

\begin{equation}\label{mom2} -\frac{\partial}{\partial y}\left(\nu_t\frac{\partial u_1}{\partial y}\right)=-\frac{\partial p}{ \partial x} \end{equation}

仔细观察方程\eqref{mom2},且其中有两个未知数$u_1$和$p$并不封闭。通常情况下,可以结合连续性方程进行处理。但在一维问题下,连续性方程自然满足。在这里需要添加另一个条件,即用户需要给定平均的流动速度(通量)。boundaryFoam即通过用户给定的平均流动速度和方程\eqref{mom2}来进行封闭求解。为方便讨论,下文中使用标量$T$表示速度分量$u_1$并略去湍流粘度。若不考虑压力梯度,方程\eqref{mom2}即为:

\begin{equation}\label{lap} -\frac{\partial}{\partial y}\left(\frac{\partial T}{\partial y}\right)=0 \end{equation}

其为一个沿$y$方向的抛物线一维热传导方程。依据不同的边界条件具有不同的分布。接下来,将压力梯度项重写为$\gradP$并考虑进来有:

\begin{equation}\label{lap2} -\frac{\partial}{\partial y}\left(\frac{\partial T}{\partial y}\right)=\gradP \end{equation}

若将$\gradP$看做热源,方程\eqref{lap2}则变为附加热源向的一维热传导方程。在boundaryFoam中,即为求解这样的一个方程,将其看做热传导方程,理解起来会变得更加明晰。在热传导方程中,给定$\gradP$(热源)之后就可以求解$T$。但往往这个压力梯度是未知的,因此boundaryFoam采用了另一种求解策略,即调用上文中用户给定的平均流动速度。

若用户定义在此一维问题上的$T$的平均值为给定值$T_\mathrm{known}$,则在求解后应有:

\begin{equation}\label{phi} \frac{\sum T^{n+1}\Delta V}{\sum\Delta V}=T_\mathrm{known} \end{equation}

其中$T^{n+1}$表示收敛后的值。方程\eqref{phi}中即表示$T$的体平均值满足用户给定的$T_\mathrm{known}$,且存在径向分布。在求解方程\eqref{lap2}的时候,需要结合方程\eqref{phi}迭代进行,步骤如下:

  1. 首先,用户给定$\gradP^n$,并求解下面的方程来获得$T^n$:

  2. \begin{equation}\label{momF1} -\frac{\partial}{\partial y}\left(\frac{\partial T^n}{\partial y}\right)=\gradP^n \end{equation}
  3. 由于$T^n$和$T_\mathrm{known}$的关系并不满足方程\eqref{phi},因此通过方程\eqref{phi}来获得修正变量$T'$:

  4. \begin{equation}\label{momF4} T'=T_\mathrm{known}-\frac{\sum T^{n}\Delta V}{\sum\Delta V} \end{equation}
  5. 更新变量$T^{n+1}$:

  6. \begin{equation}\label{momF5} T^{n+1}=T+T'=T+\left(T_\mathrm{known}-\frac{\sum T^{n}\Delta V}{\sum\Delta V}\right) \end{equation}
  7. 下一步更新变量$\gradP^{n+1}$。在更新的时候,需要获得修正变量$\gradP'$。考虑在最终稳态的情况下有:

    \begin{equation}\label{momF2} -\frac{\partial}{\partial y}\left(\frac{\partial T^{n+1}}{\partial y}\right)=\gradP^{n+1} \end{equation}

    方程\eqref{momF2}减去方程\eqref{momF1}有:

    \begin{equation}\label{momF3} -\frac{\partial}{\partial y}\left(\frac{\partial T'}{\partial y}\right)={\gradP'} \end{equation}

    将$T'$和方程\eqref{momF3}结合可用于获得$\gradP'$以及$\gradP^{n+1}$:

  8. \begin{equation}\label{momF6} \gradP^{n+1}=\gradP^{n}+\gradP'=\gradP^{n}-\frac{\partial}{\partial y}\left(\frac{\partial T'}{\partial y}\right) \end{equation}
  9. 随后依据$\gradP^{n+1}$进入到步骤1进行下一个时间步的迭代,直至满足方程\eqref{phi}即收敛。

代码分析

主要分析boundaryFoam.C中的主要代码部分:

//方程\eqref{lap2}
fvVectorMatrix UEqn
(
    divR == gradP + fvOptions(U)
);

UEqn.solve();

...
 
//方程\eqref{phi}左边的第一项     
dimensionedVector UbarStar = flowMask & U.weightedAverage(mesh.V());

//方程\eqref{momF5}
U += (Ubar - UbarStar);

//方程\eqref{momF6} 
gradP += (Ubar - UbarStar)/(1.0/UEqn.A()().weightedAverage(mesh.V());

更新历史
2018.06.28 大修

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网

');