icoFoam解析


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

1. 引言

计算流体力学 (CFD) 采用计算机模拟流体的流动行为。流体动力学包含各种类型,如可压缩流动、不可压缩流动、多相流动以及化学反应、组分输运等。所有的这些流动性为都可以用偏微分方程来描述。例如,不可压缩牛顿流体的控制方程可以表示为 (Ferziger and Peric, 2012):

\begin{equation} \nabla\cdot\bfU=0, \label{C} \end{equation}
\begin{equation} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nu \nabla \mathbf{U}), \label{mom} \end{equation}

其中$\rho$为密度,$\bfU$为速度,$p$为压力,$\nu$为粘度。该方程具有以下特点:

  • 方程\eqref{mom}中左边第二项是关于$\bfU$乘积的偏导数,这种未知量和未知量乘积的问题构成非线性问题,CFD对非线性问题需要特殊处理。另一方面,非线性的双曲问题的解可能会存在间断(如激波)。激波通常存在于高超声速的欧拉问题求解中。同时,非线性项也是湍流在数学方程中的体现;
  • 方程\eqref{mom}的数学特征为抛物线。不同数学特征的问题需要调用不同的时间离散格式,隐性时间格式更有利于求解抛物线问题。若方程\eqref{mom}中省略若干项则会改变方程的数学特征,例如若将方程右侧置为$0$,则变为双曲特征的欧拉方程。欧拉方程得益于其双曲特性,可采用迎风类显性算法推进,各种基于有限体积法的高分辨率格式因此而生(交错网格中心格式、中心-迎风格式等);
  • 方程\eqref{C}和\eqref{mom}中存在未知量压力$p$,同时压力和速度是耦合在一起的(二者相互影响),但并没有单独的压力方程。这导致压力的求解需要特殊的策略。这也是CFD中压力基SIMPLE/PISO算法、密度基、耦合/解耦算法要处理的问题。另一方面,方程\eqref{mom}右侧存在压力梯度项,所有类似的梯度项都需要特殊的数值处理否则会引起数值震荡;
  • 方程\eqref{C}和\eqref{mom}为宏观方程,调用了宏观假定,其起源为玻尔兹曼方程(又名动理学方程)。在更底层的介尺度研究领域,方程\eqref{C}和\eqref{mom}也即从介尺度模型演化的宏观二阶矩模型。在无压力无粘性的条件下具备弱双曲特征。由于失去了高阶矩的统计学特征,因此方程\eqref{C}和\eqref{mom}在某些情况下是不适用的;

在考虑适当的简化情况下,方程\eqref{C}和\eqref{mom}可被称之为Navier-Stokes (N-S) 方程。除了上述数学特征,求解N-S方程可分为不同的方法如有限体积法、有限差分法以及谱方法等。有限体积法得益于其天然守恒的数学特征,被大量的植入于CFD软件(如ANSYS Fluent、OpenFOAM等)中用于流体动力学计算模拟。

icoFoam为OpenFOAM中一个基于有限体积法的,用于求解不可压缩牛顿流体N-S方程的求解器。它可用于模拟层流流动或流场剪切力引致的湍流准直接模拟 (quasi-DNS),也被称之为准DNS。在这里需要注意的是,得益于谱方法以及有限差分法高阶精度的易用性,DNS计算通常采用谱方法或有限差分法。有限体积法由于本身在计算通量时引入的中点插值法则,大大限制了有限体积法的精度。然而Komen等通过研究圆管内的湍流流动,认为在网格分辨率足够的情况下,使用有限体积法同样可以获得高质量的平均速度以及脉动速度流场 (Komen et al., 2014)。Komen等进一步的使用有限体积法进行准DNS计算,并与有限体积法大涡模拟进行标定,认为大涡模拟方法会引致一定的数值耗散 (Komen et al., 2017)。Axtmann与Rist也使用OpenFOAM进行了准DNS直接模拟,并研究了编译器与并行计算的特性 (Axtmann and Rist, 2016) 。liu等通过icoFoam对开口顶盖驱动流进行了准直接模拟 (Liu et al., 2016) ,进一步验证了使用OpenFOAM进行准直接模拟的可行性。

在icoFoam中,压力和速度的耦合通过瞬态PISO算法进行计算 (Issa, 1986) ,且没有考虑其他体积力(如重力等)。除去三个基本求解器laplacianFoam, potentialFoam, scalarTransportFoam,icoFoam为一个最基本的求解器范例,可用于理解 OpenFOAM中压力速度耦合的框架策略。需要注意的是对于流场内温度变化较大,但可忽略浮力以及重力的可压缩求解器,可选用rhoPimpleFoam。附加重力的可压缩浮力驱动流求解器可选用buoyantPimpleFoam。本文从最基本的N-S方程在笛卡尔网格下的有限体积离散开始推导,介绍icoFoam求解器中的算法并和OpenFOAM植入的代码相对应,可用于理解N-S方程在有限体积法中的离散以及瞬态PISO算法在OpenFOAM中的实现。需要注意的是,由于为不可压缩流体求解器,下文中$p$的单位为压力除以密度的单位,即$\mathrm{m}^2\cdot \mathrm{s}^{-2}$。在下文中,$p$也对应方程\eqref{mom}中的$p/\rho$。

2. 控制方程与算法
2.1 有限体积法离散

首先对方程\eqref{mom}进行对速度$\bfU$关于时间$t$的欧拉全隐离散有 (Jasak, 1996):

\begin{equation} \int \int {\frac{{\partial \mathbf{U}}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {\mathbf{U}_\mathrm{P}^r - \mathbf{U}_\mathrm{P}^n} \right)\;} {V_\mathrm{P}}, \label{mom1} \end{equation}
\begin{equation} \int \int {\nabla \cdot \left( {\mathbf{U}\otimes\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\mathbf{U}\otimes\mathbf{U}\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\mathbf{U}^n}\otimes{\mathbf{U}^r}} \right)}_f}} \bfS_f\Delta t = \Delta t\sum {\left( {F_f^n \mathbf{U}_f^r} \right)} , \label{mom2} \end{equation}
\begin{equation} {\int \int \nabla \cdot \left( {\nu \nabla \mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\nu \nabla \mathbf{U}\mathrm{d}\bfS\mathrm{d}t = \sum {\left( {\nu \nabla {\mathbf{U}^r}} \right)} } _f}\bfS_f\Delta t = \Delta t\sum {\left( \nu \nabla\mathbf{U}_f^r\cdot\bfS_f \right)} , \label{mom3} \end{equation}
\begin{equation} \int \int \nabla p \mathrm{d}V\mathrm{d}t=\int\int p \mathrm{d}\bfS\mathrm{d}t=\Delta t\sum \left(p_f^n\bfS_f\right), \label{gradp} \end{equation}

其中上标$^n$表示为当前的时间步(已知),上标$^r$表示预测时间步(待求),下标$_f$表示网格单元面上的值,$V_\rP$表示网格单元体积,$\bfS$表示网格单元的各个面的面矢量,$\Delta t$表示时间步长,$F_f$为通量,$\nu$为动力粘度。方程\eqref{mom3}中$\nabla\mathbf{U}_f^r\cdot\bfS_f$表示面法向梯度,其在正交网格中进一步可以表示为:

\begin{equation} \nabla\mathbf{U}_f^r\cdot\bfS_f= \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^r - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}, \label{sngrad} \end{equation}

其中$d$表示网格单元体心之间的距离,下标$_\mathrm{N}$表示相邻网格单元的速度,下标$_\mathrm{P}$表示当前网格单元的速度。面法向梯度在三维情况下为一个矢量。在这里需要注意的是,方程\eqref{sngrad}只对矩形网格精准,网格非正交性增加,需要特定的数值处理提高其计算精准性 (Demirdžić,2015) 。将方程\eqref{mom1}-\eqref{sngrad}代入到方程\eqref{mom}有:

\begin{equation} \frac{{\mathbf{U}_\mathrm{P}^r - \mathbf{U}_\mathrm{P}^n}}{{\Delta t}}{V_\mathrm{P}} + \sum {\left( {F_f^n \mathbf{U}_f^r} \right)} = \sum {\left( {\nu \nabla\mathbf{U}_f^r\cdot\bfS_f} \right)} -\sum \left(p_f^n\bfS_f\right). \label{momF} \end{equation}

需要注意的是,方程\eqref{mom}中的左边第二项(对流项)是非线性的。在求解的过程中,要么选择非线性求解器,要么将对流项线性化。由于非线性求解器非常复杂,因此在OpenFOAM均采用线性化处理。具体的,在对方程\eqref{mom2}中的左边第二项对流项离散中,其中的通量$F_f$采用当前已知时间步的速度$\mathbf{U}^n$来计算,同时保留预测速度$\mathbf{U}^r_\mathrm{P}$为未知量。这种将高阶非线性项降为一阶线性项的过程即为线性化操作 (Versteeg and Malalasekera, 2007; 陶文铨, 2001) 。线性化的一个问题即为通量速度的信息有所滞后 (Jasak, 1996) 。

2.2. 动量预测

方程\eqref{momF}中$\mathbf{U}_f^r$被定义为面上的预测速度,然而在计算中求得的均为网格单元体心上的速度。因此需要从体心速度进行插值来获得面速度,在此步可以引入各种插值格式。假设使用中心线性格式:

\begin{equation} \mathbf{U}_f^r = \frac{{\mathbf{U}_\mathrm{P}^r + \mathbf{U}_\mathrm{N}^r}}{2}. \label{linear} \end{equation}
\begin{equation} p_f^n = \frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}. \label{linear2} \end{equation}

将方程\eqref{sngrad},\eqref{linear}代入到方程\eqref{momF}有

\begin{equation} \frac{{\mathbf{U}_\mathrm{P}^r - \mathbf{U}_\mathrm{P}^n}}{{\Delta t}}{V_\mathrm{P}} + \sum {\left( {F_f^n \frac{{\mathbf{U}_\mathrm{P}^r + \mathbf{U}_\mathrm{N}^r}}{2}} \right)} = \sum {\left( {\nu \left| \bfS_f \right|\frac{{\mathbf{U}_\mathrm{N}^r - \mathbf{U}_\mathrm{P}^r}}{{\left| d \right|}}} \right)}-\sum \left(\frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f\right), \end{equation}
\begin{equation} \left( {\frac{{{V_\mathrm{P}}}}{{\Delta t}} + \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right)\mathbf{U}_\mathrm{P}^r = - \sum {\left( {\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)\mathbf{U}_\mathrm{N}^r} \right)} + \frac{{{V_p}}}{{\Delta t}}\mathbf{U}_\mathrm{P}^n-\sum \left(\frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f\right). \end{equation}

将上式简写为

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^r{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^r} = S_\mathrm{P}^n-\sum \left(\frac{{p_\mathrm{P}^n + p_\mathrm{N}^n}}{2}\bfS_f\right), \label{apanmom} \end{equation}

其中

\begin{equation} A_\mathrm{P}=\left( {\frac{{{V_p}}}{{\Delta t}} + \sum {\left( {\frac{{F_f^n}}{2}} \right) + \sum {\left( {\nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right)} } } \right), \end{equation}
\begin{equation} A_\mathrm{N}=\left( {\frac{{F_f^n}}{2} - \nu \frac{{\left| \bfS_f \right|}}{{\left| d \right|}}} \right), \end{equation}
\begin{equation} S_\mathrm{P}^n=\frac{{{V_p}}}{{\Delta t}}\mathbf{U}_\mathrm{P}^n. \end{equation}

可以看出在某个时间步内$A_\mathrm{P}$、$A_\mathrm{N}$和$S_\mathrm{P}^n$保持不变。求解方程\eqref{apanmom}即可获得预测速度$\mathbf{U}^r_\mathrm{P}$。方程\eqref{apanmom}即为OpenFOAM中的动量预测方程。需要注意的是,动量预测步骤并不是必须的,这主要和速度压力耦合求解策略有关。如果不调用动量预测步骤,则$\mathbf{U}^r_\mathrm{P}=\mathbf{U}^n_\mathrm{P}$。

2.3. 压力泊松方程

N-S方程求解的关键问题之一在于并没有压力的方程出现。仔细分析得知上文中求解动量方程获得了速度,但连续性方程并不是关于压力的方程。压力泊松方程的构建正式基于连续性方程而来,其通过对动量方程继续做散度并相加获得。考虑最终收敛的情况下,对连续性方程$\nabla\cdot\bfU=0$进行离散后的形式为

\begin{equation} \sum \left(\mathbf{U}_{\mathrm{P},f}^{n+1} \cdot \bfS_f\right)=0. \label{cont} \end{equation}

如果能用压力表示方程\eqref{cont}中的$\mathbf{U}_{\mathrm{P},f}^{n+1}$,是不是就可以得出压力泊松方程了呢?答案是肯定的。首先,方程\eqref{apanmom}可以写为:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{r}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}} = S_\mathrm{P}^n-\nabla p_\mathrm{P}^{n}, \label{apanmomr} \end{equation}

收敛情况下为:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{n+1}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{n+1}} = S_\mathrm{P}^n-\nabla p_\mathrm{P}^{n+1}, \label{apanmomNew} \end{equation}

若将方程\eqref{apanmomNew}减去\eqref{apanmomr}有:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}'+\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}'} = -\nabla p_\mathrm{P}', \label{simple} \end{equation}

其中$\mathbf{U}_\mathrm{P}'=\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U}_\mathrm{P}^{r}$(其他同理)。在这里未调用任何略去临点的假定,方程\eqref{apanmomNew}是严格成立的(精准的)。对方程\eqref{apanmomNew}移项有

\begin{equation} \mathbf{U}_\mathrm{P}^{n+1} = \mathbf{HbyA}^{n+1}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\nabla {p_\rP^{n+1}}. \end{equation}

其中

\begin{equation} \mathbf{HbyA}^{n+1}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{n+1}} + S_\mathrm{P}^n} \right). \label{hbya} \end{equation}

类似的,通过Rhie-Chow插值,面上的速度可以通过下面的公式获得 (Rhie and Chow, 1983):

\begin{equation} \mathbf{U}_{\mathrm{P},f}^{n+1} = \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}}. \label{uf} \end{equation}

这样,当前网格面上的速度$\bfU_\rP$即和面法向压力梯度联立起来。面法向压力梯度进一步和当前网格与相邻网格的压力有关,即1-$\delta$压差,可用来消除棋盘分布。将方程\eqref{uf}代入到方程\eqref{cont}有:

\begin{equation}\label{eg} \sum \left( \left( \mathbf{HbyA}_f^{n+1} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{n+1}} \right) \cdot \bfS_f \right)=0, \end{equation}
\begin{equation} \sum {\left( {\mathbf{HbyA}_f^{n+1}} \cdot \bfS_f\right) } = \sum {\left( {\frac{1}{{A_{\mathrm{P},f}}}{\nabla _f}{p^{n+1}}} \cdot \bfS_f\right) }, \label{jim} \end{equation}

方程\eqref{jim}即下述方程的离散形式:

\begin{equation} \nabla \cdot (\mathbf{HbyA}^{n+1}) = \nabla \cdot\left(\frac{1}{A_{\mathrm{P}}} \nabla p^{n+1}\right). \label{poss} \end{equation}

方程\eqref{poss}即为压力泊松方程。若有$\mathbf{HbyA}^{n+1}$,则可求得收敛的压力。

2.4. 迭代算法

可以看出,最后求解的速度$\mathbf{U}_\mathrm{P}^{n+1}$和压力$p^{n+1}$应该符合方程\eqref{apanmom}和\eqref{poss},分别对应动量方程和连续性方程。然而目前,我们只有通过2.1节动量预测步骤求出来的$\mathbf{U}_\mathrm{P}^r$。$\mathbf{HbyA}_\mathrm{P}^{n+1}$和$p^{n+1}$均为未知的。压力泊松方程不能求解。1986年Issa提出了PISO (Pressure Implicit with Splitting of Operators) 算法可以解决这个问题 (Issa, 1986) 。PISO算法在提出时主要针对不可压缩非稳态计算,其为一种在时间补内非迭代的算法,在随后也被拓展用于稳态问题中。PISO算法通过对当前时间步内的多次修正来获得最终的$\mathbf{U}_\mathrm{P}^{n+1}$和$p^{n+1}$。

依据PISO算法,在方程\eqref{simple}的基础上引入略去临点影响的假定有:

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}' = -\nabla p_\mathrm{P}', \label{uprime} \end{equation}

将方程\eqref{uprime}和\eqref{apanmomr}相加有

\begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^{*}{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^{r}} = S_\mathrm{P}^n-\nabla p_\mathrm{P}^{*}, \label{pisoapan} \end{equation}

其中$\mathbf{U}_\mathrm{P}^{*}=\mathbf{U}_\mathrm{P}^{r}+\mathbf{U}_\mathrm{P}'$,$p_\mathrm{P}^{*}=p_\mathrm{P}^{n}+p_\mathrm{P}'$。相比方程\eqref{apanmomNew},方程\eqref{pisoapan}由于忽略了临点的影响并不是精准的。对方程\eqref{pisoapan}进行移项有

\begin{equation} \mathbf{U}_\mathrm{P}^* = \mathbf{HbyA}^{r}_\mathrm{P} - \frac{1}{{{A_\mathrm{P}}}}\nabla {p^{*}}. \end{equation}
\begin{equation} \mathbf{HbyA}^{r}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{r}} + S_\mathrm{P}^n} \right). \label{hbya2} \end{equation}
\begin{equation} \mathbf{U}_{\mathrm{P},f}^{*} = \mathbf{HbyA}_f^{r} - \frac{1}{{{A_{\mathrm{P},f}}}}\nabla_f {p^{*}}. \label{pisouf} \end{equation}

方程\eqref{pisouf}中可参考方程\eqref{eg}-\eqref{poss}的步骤组建压力泊松方程,即

\begin{equation} \nabla \cdot (\mathbf{HbyA}^{r}) = \nabla \cdot\left(\frac{1}{A_{\mathrm{P}}} \nabla p^{*}\right). \label{pisoposs} \end{equation}

对方程\eqref{pisoposs}求解后有$p^{*}$。将$p^*$回代到方程\eqref{pisouf}的时候,有$\mathbf{U}_{\mathrm{P}}^*$。但需要注意的是,方程\eqref{pisouf}并不是精准的。因此$p^{*}$和$\mathbf{U}_{\mathrm{P}}^*$并不是最终时间步$n+1$的结果。通过求解多次的压力泊松方程,这种滞后性逐渐消除并可最终忽略 (Issa, 1986)。总而言之,PISO算法中的的迭代过程可以表示为下面几个步骤:

  1. 依据初始条件求解预测速度$\mathbf{U}^r$,组建$\mathbf{HbyA}^r$;
  2. 求解方程\eqref{pisoposs}求解压力$p^*$;
  3. 通过方程\eqref{pisouf}依据压力$p^*$以及$\mathbf{HbyA}^r$更新速度有$\mathbf{U}^{*}$;
  4. 回到第二步迭代求解几次即可;
3. 代码分析
#include "fvCFD.H"
//有关头文件等参考laplacianFoam, potentialFoam, scalarTransportFoam解析

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        solve(UEqn == -fvc::grad(p));
        //对应方程(12)

        for (int corr=0; corr< nCorr; corr++)
        {
            volScalarField rAU(1.0/UEqn.A());
            //方程(13)中A_p的倒数

            volVectorField HbyA("HbyA", U);
            HbyA = rAU*UEqn.H();
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::ddtPhiCorr(rAU, U, phi)
            );
            //方程(25)中左边的第一项

            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );
                //方程(26)

                pEqn.solve();

                if (nonOrth == nNonOrthCorr)
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }
            U = HbyA - rAU*fvc::grad(p);
            //方程(23)
//方程\eqref{mom1}的离散格式
ddtSchemes
{
    default         Euler;
}
//方程\eqref{gradp}的离散格式
gradSchemes
{
    default         Gauss linear;
}
//方程\eqref{mom2}中$\bfU_f^n$的离散格式
divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}
//方程\eqref{mom3}的离散格式
laplacianSchemes
{
    default         Gauss linear orthogonal;
}
//例如方程\eqref{mom2}中$F_f^n$的离散格式
interpolationSchemes
{
    default         linear;
}
//例如方程\eqref{sngrad}的离散格式
snGradSchemes
{
    default         orthogonal;
}
4. 验证算例

icoFoam自带若干算例。在OpenFOAM-6中,自带的算例主要为以下几个:

上述算例在《OpenFOAM用户指南》已经进行了非常详细的阐述。在此不做过多介绍。除上述算例外,icoFoam也可以进行直接模拟。icoFoam不附加任何湍流模型直接求解N-S方程,并采用二阶有限体积法进行离散求解。因此可用作为一个直接模拟求解器进行直接模拟。一些文章中更严谨的将其称为准直接模拟(quasi-DNS)。另外,在这里强调OpenFOAM中的icoFoam、pisoFoam以及pimpleFoam在不使用湍流模型的情况下是相同的。

本算例对标相关文献的研究工作 (Komen et al., 2017)。


采用icoFoam模拟的壁面湍流

4.1. 模型设置

算例网格为采用blockMesh生成的3D纯六面体结构网格。左侧为周期速度进口,右侧为周期速度出口,上下为用于附加剪切的壁面。需要注意的是,对于DNS以及LES,流场本身的湍流结构以及进口边界条件非常敏感。在不合理的初始条件下,二阶精度的有限体积法可能并不会捕获精细的涡结构。因此初始场应该尽可能的精准。本算例初始场通过boxTurb给定湍流涡结构之后,运行了200个周期,作为稳定的初始场。$\mathrm{Re}_\tau\approx 180$。

constant文件夹的物性主要设置如下:

4.2. 算例运行

相关算例运行首先需要执行网格转换命令:

blockMesh

网格生成完毕,后直接运行直接模拟求解器运行即可:

pisoFoam

在这里提醒,pisoFoam是icoFoam的进阶版本。由于icoFoam中不支持fvOptions,因此本算例采用pisoFoam进行计算。本算例比较耗费计算机资源,采用12核并行计算大约需要10个小时。在运行之后,针对平均速度UMean即可以求出u+以及y+(如下图)。


点击下载

因为初始场信息已经经过计算200个周期,因此文件较大,约65m。

参考文献

Komen, E., Shams, A., Camilo, L., Koren, B., 2014. Quasi-DNS capabilities of OpenFOAM for different mesh types. Computers & Fluids, 96, 87-104.

Komen, E., Camilo, L. H., Shams, A., Geurts, B. J., Koren, B., 2017. A quantification method for numerical dissipation in quasi-DNS and under-resolved DNS, and effects of numerical dissipation in quasi-DNS and under-resolved DNS of turbulent channel flows. Journal of Computational Physics, 345, 565-595.

Axtmann, G., Rist, U., 2016. Scalability of OpenFOAM with Large Eddy Simulations and DNS on High-Performance Systems. In High Performance Computing in Science and Engineering, 16, 413-424. Springer.

Liu, Q., Gómez, F., Pérez, J. M., Theofilis, V., 2016. Instability and sensitivity analysis of flows using OpenFOAM. Chinese Journal of Aeronautics, 29, 316-325.

Issa, R. I., 1986. Solution of the implicitly discretised fluid flow equations by operator-splitting. Journal of Computational Physics, 62, 40-65.

Ferziger, J. H., Peric, M., 2012. Computational methods for fluid dynamics. Springer Science Business Media.

Jasak, H., 1996. Error analysis and estimation for finite volume method with applications to fluid flow. PhD Thesis. Imperial College London.

Demirdžić, I., 2015. On the discretization of the diffusion term in finite-volume continuum mechanics. Numerical Heat Transfer, Part B: Fundamentals, 68, 1-10.

Versteeg, H. K., Malalasekera, W., 2007. An introduction to computational fluid dynamics: the finite volume method. Pearson Education.

陶文铨, 2001. 数值传热学. 西安交通大学出版社.

Rhie, C. M., Chow, W. L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal, 21, 1525-1532.

更新历史
2019.10.12更正一些符号 | 2019.6.12去掉一些方程中的$1/V_\rP$ | 2018.12.08丁长友:修正方程\eqref{cont}中的$\mathbf{U}_{f}$为$\mathbf{U}_{\mathrm{P},f}$,方程\eqref{sngrad}中$\mathbf{U}_{N}^n$为$\mathbf{U}_{N}^r$ | 2018.11.07更新方程索引 | 2018.08.22于鸿翔:修正方程7-10、14中的$\mathbf{U}_\mathrm{N}^n$为$\mathbf{U}_\mathrm{N}^r$ | 2018.07.05尹胜强:修正方程26、27的对应关系 | 2018.05.21JimShi:将方程22中的$\bfS$移至括号内 | 2018.05.08去掉方程\eqref{poss}中的下标$_f$ | 2018.04.15冯强:增加参考文献更正错别字 | 2017.08.13增加方程\eqref{uprime}

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