return 0;
wmake

CFD中的能量方程

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


能量守恒定义为绝热系统的总能量是一个常数。能量的单位为$\mathrm{J}=\mathrm{N}\cdot \mathrm{m}=\frac{\mathrm{kg}\cdot \mathrm{m}}{\mathrm{s}^2}\mathrm{m}$,在CFD中,通常用$E\left[\frac{\mathrm{m}^2}{\mathrm{s}^2}\right]$来表示比能(specific energy),其表示每单位质量的总能量,其值不随时间变化。比能通常由三方面构成,即单位质量的动能(kinetic energy)以及单位质量的内能(internal energy)。流体微元能量的变化可来源于自己产生/消失的能量、流入/流出的能量、以及做功增加/减少的能量。

总能量方程

首先,定义单位质量的动能为$K$:

\begin{equation} K=\frac{1}{2}|\mathbf{U}|^2 \label{K} \end{equation}

其中$\mathbf{U}$为速度。单位质量的内能定义为$e$。这样,流体微元的比能$E$可表示为

\begin{equation} E=K+e \label{E} \end{equation}

比能每单位体积的时间变化率即为每单位质量的比能变化率乘以密度:

\begin{equation} \rho\frac{\rD E}{\rD t}\left[\frac{\mathrm{kg}}{\mathrm{s}^3\cdot\mathrm{m}}\right]=\rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} \label{A} \end{equation}

对于体积为$\mathrm{d}x\mathrm{d}y\mathrm{d}z$的流体微元,其比能的时间变化率为

\begin{equation} \rho\frac{\rD E}{\rD t}\mathrm{d}x\mathrm{d}y\mathrm{d}z\left[\frac{\mathrm{kg}\cdot\mathrm{m}^2}{\mathrm{s}^3}\right]=\rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{AA} \end{equation}

方程\eqref{AA}构成能量方程左边的部分。接下来考虑能量方程右边的部分。正如之前所说的,流体微元能量的变化可能来自于自发产生的能量、流入流出的能量以及做功产生消失的能量。首先考虑内部产生的能量。定义$r\left[\frac{\mathrm{m}^2}{\mathrm{s}^3}\right]$为比热源,其表示每单位时间、每单位质量的能量产生。那么流体微元的净比热源可表示为

\begin{equation} \rho r \mathrm{d}x\mathrm{d}y\mathrm{d}z\left[\frac{\mathrm{kg}\cdot\mathrm{m}^2}{\mathrm{s}^3}\right] \label{heatSource} \end{equation}

下一步,考虑流入流出的能量。热传导(热流具有方向性)对流体微元的加热为(热通量矢量定义为$\mathbf{q}\left[\frac{\mathrm{w}}{\mathrm{m}^2}=\frac{\mathrm{kg}}{\mathrm{s}^3}\right]$)

\begin{equation} - \nabla \cdot \mathbf{q} \mathrm{d}x\mathrm{d}y\mathrm{d}z\left[\frac{\mathrm{kg}\cdot\mathrm{m}^2}{\mathrm{s}^3}\right] \label{heatFlux} \end{equation}

依据傅里叶定律,$\mathbf{q}$和温度$T$的关系为

\begin{equation} \mathbf{q}=-k \nabla T \label{fourier} \end{equation}

接下来考虑做功引起的能量变化。单位时间内做功的变化为功率。作用在流体微元上的力,对其做功的功率等于这个力乘以速度在此力作用方向上的分量。因此压力的功率为:

\begin{equation}\label{pressure} -\nabla\cdot(p\mathbf{U})\mathrm{d}x\mathrm{d}y\mathrm{d}z= -\left(\frac{\p up}{\p x}+\frac{\p vp}{\p y}+\frac{\p wp}{\p z}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z \end{equation}

剪切力对流体微团做功的功率为:

\begin{equation}\label{stress} \nabla \cdot(\mathbf{\tau} \cdot \mathbf{U})\mathrm{d}x\mathrm{d}y\mathrm{d}z \end{equation}

重力矢量$\mathbf{g}$以体积力的方式作用于流体微元,引起重力势能的变化。单位时间内,重力做功的功率为

\begin{equation} \rho \mathbf{g} \cdot \mathbf{U} \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{g} \end{equation}

结合方程\eqref{AA}、\eqref{heatSource}、\eqref{heatFlux}、\eqref{pressure}、\eqref{stress}、\eqref{g}有最终的总能量方程:

\begin{equation} \rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} = \rho r -\nabla\cdot\bfq -\nabla\cdot(p\mathbf{U})+\nabla \cdot(\tau \cdot \mathbf{U})+ \rho \mathbf{g} \cdot \mathbf{U} \label{EEqn} \end{equation}

其展开形式为

\begin{equation} \rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} = \rho r-\frac{\p \bfq_x}{\p x}-\frac{\p \bfq_y}{\p y}-\frac{\p \bfq_z}{\p z}-\left(\frac{\p up}{\p x}+\frac{\p vp}{\p y}+\frac{\p wp}{\p z}\right) \\ +\frac{\p u \tau_{xx}}{\p x}+\frac{\p u \tau_{yx}}{\p y}+\frac{\p u \tau_{zx}}{\p z} +\frac{\p v \tau_{xy}}{\p x}+\frac{\p v \tau_{yy}}{\p y}+\frac{\p v \tau_{zy}}{\p z} +\frac{\p w \tau_{xz}}{\p x}+\frac{\p w \tau_{yz}}{\p y}+\frac{\p w \tau_{zz}}{\p z}+ \rho \mathbf{g} \cdot \mathbf{U} \label{EEqnP} \end{equation}

其中的$\nabla \cdot(\tau \cdot \mathbf{U})=\left(\nabla \cdot \tau \right)\cdot \mathbf{U}+\tau:\nabla \mathbf{U}$。从方程\eqref{EEqn}左边可以看出其采用的是非守恒形式,且目前尚未封闭。

内能方程

虽然方程\eqref{EEqn}为最终的能量方程。但是在CFD中通常的做法是抽离动能项来获得一个内能方程。内能方程可以从能量方程中减去动能(机械能)方程来获得。首先有动量方程:

\begin{equation} \rho \frac{\mathrm{D}\mathbf{U}}{\mathrm{D}t}=\nabla \cdot \tau + \rho \mathbf{g} - \nabla p \label{mom} \end{equation}

将动量方程中的每个速度分量方程乘以速度分量并加和有:

\begin{equation} \rho \frac{\mathrm{D} \left[\frac{1}{2}\left(u_x^2+u_y^2+u_z^2 \right)\right]}{\mathrm{D}t} =\left(\nabla \cdot \tau \right) \cdot \mathbf{U} + \rho \mathbf{g} \cdot \mathbf{U} -\mathbf{U} \cdot \nabla p \label{mom2} \end{equation}

即:

\begin{equation} \rho\frac{\mathrm{D}K}{\mathrm{D}t} =\left(\nabla \cdot \tau \right) \cdot \mathbf{U} + \rho \mathbf{g} \cdot \mathbf{U} -\mathbf{U} \cdot \nabla p \label{K2} \end{equation}

将方程\eqref{EEqn}中减去方程\eqref{K2}有内能方程:

\begin{equation} \rho\frac{\mathrm{D}e}{\mathrm{D}t} =-p\nabla \cdot \mathbf{U}-\nabla\cdot\bfq + \tau : \nabla \mathbf{U} + \rho r \label{e} \end{equation}

可以看出,内能方程中不含有体积力项。同样,方程\eqref{e}采用的是非守恒形式,同样尚未封闭。

焓方程

对于可压缩流体,通常把总能量方程简化为焓方程。首先有比焓$h$以及总比焓$h_0$的定义:

\begin{equation} h=e+\frac{p}{\rho} \label{hDef} \end{equation}
\begin{equation} h_0=h+K=e+\frac{p}{\rho}+K \label{h0Def} \end{equation}

回到非守恒形式的方程\eqref{EEqn},有守恒形式的总能量方程(调用连续性方程):

\begin{equation} \frac{\partial \rho e}{\partial t}+\nabla \cdot (\rho \mathbf{U} e) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K)= -\nabla\cdot(p\mathbf{U})+\rho r -\nabla\cdot\bfq + \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{EEqnconser} \end{equation}

把方程\eqref{hDef}代入到方程\eqref{EEqnconser}中有最终的守恒形式的比焓方程:

\begin{equation} \frac{\partial \rho h}{\partial t}+\nabla \cdot (\rho \mathbf{U} h) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K) =\frac{\partial p}{\partial t}+ \rho r -\nabla\cdot\bfq + \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{hEqn} \end{equation}

以及守恒形式的总比焓方程:

\begin{equation} \frac{\partial \rho h_0}{\partial t}+\nabla \cdot (\rho \mathbf{U} h_0) =\frac{\partial p}{\partial t}+ \rho r -\nabla\cdot\bfq+ \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{h0Eqn} \end{equation}

再一次的,总比焓方程尚未封闭。

封闭

在开源CFD软件中,假设$\mathbf{q}=-\alpha_\mathrm{eff}\nabla e$或$\mathbf{q}=-\alpha_\mathrm{eff}\nabla h$,这样方程\eqref{EEqnconser}以及方程\eqref{hEqn}则简化为(无热源$r=0$):

\begin{equation} \frac{\partial \rho e}{\partial t}+\nabla \cdot (\rho \mathbf{U} e) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K)- \nabla \cdot (\alpha_\mathrm{eff}\nabla e)= -\nabla\cdot(p\mathbf{U}) \label{Efinal} \end{equation}
\begin{equation} \frac{\partial \rho h}{\partial t}+\nabla \cdot (\rho \mathbf{U} h) + \frac{\partial \rho K}{\partial t}+\nabla \cdot (\rho \mathbf{U} K) - \nabla \cdot (\alpha_\mathrm{eff}\nabla h) =\frac{\partial p}{\partial t} \label{hfinal} \end{equation}

其中的$K$通过方程\eqref{K}更新,且其中的速度项通过分离/耦合式算法求解,因此得以封闭。

可见,总能量方程、内能方程和焓方程在求解的时候哪一个最优呢?目前从经验来看,传热、可压缩、化学反应求解器使用的主要为守恒形式的总能量方程和比焓方程。并且,方程\eqref{EEqnconser}和方程\eqref{hEqn}中的$\nabla \cdot(\tau \cdot \mathbf{U})$以及$\rho \mathbf{g} \cdot \mathbf{U}$通常被忽略。但是在rhoCentralFoam中,则植入了此项。一些文献表示,能量方程的选择在某些情况下是至关紧要的。比如在进行激波捕获计算的时候,求解总能量方程要比求解内能方程结果精确的多。在不可压缩流动中,唯一重要的是内能方程。在考虑传热的情况下,内能要比热能小得多。只要温度对流体的影响不是很重要,那么能量方程可以放在压力速度耦合之后进行求解。在这种情况下是一种单向耦合。同时,需要注意的是内能方程实际上来源于动量方程。但是总能量方程以及焓方程却可以单独推导。在某些情况下求解内能方程可能会带来一些问题。同时,对于不可压缩流体,能量的增加很少见,能量的耗散却比较明显。对于可压缩流,双方都比较重要。

Boussinesq假定

下面考虑传热领域的Boussinesq假定,而不是湍流以及本构方程领域的Boussinesq近似/假定。

在传热领域内,Boussinesq假定认为在流动中温度的变化是非常小的,因此对应的密度变化也非常小。所以在控制方程中,除了浮力项$\rho\bfg$,其他项的密度被认为是常数。Boussinesq假定使得方程的非线性特性降低。但由于Boussinesq假定的限定条件,也使得调用Boussinesq假定的求解器存在一定的限制。在工程中,Boussinesq假定主要用于室温下的液体对流、建筑物对流、气相分散等。在温度变化比较大的情况下,不建议使用Boussinesq假定。Boussinesq假定认为流体的密度可以这样计算:

\begin{equation}\label{boussi} \rho=\rho_\rref\left(1-\beta\left(T-T_\rref\right)\right) \end{equation}

其中$\rho$表示流体的密度,$\rho_\rref$表示流体的参考密度,$T$表示流体的温度,$T_\rref$表示流体的参考温度,$\beta$表示流体的体膨胀系数。现考虑可压缩流体并附加重力的N-S方程有:

\begin{equation} \frac{\partial \rho\mathbf{U}}{\partial t}+\nabla \cdot (\rho\mathbf{U} \otimes\mathbf{U})=-\nabla p+\nabla \cdot(\mu (\nabla \mathbf{U}+\nabla\bfU^\bfT))-\frac{2}{3}\mu(\nabla\cdot\bfU)\bfI+\rho\bfg \end{equation}

依据Boussinesq假定,认为浮力项外的密度为常数,因此提出$\rho$有:

\begin{equation}\label{MM} \rho\frac{\partial \mathbf{U}}{\partial t}+\rho\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla p+\nabla \cdot(\mu (\nabla \mathbf{U}+\nabla\bfU^\bfT))-\frac{2}{3}\mu(\nabla\cdot\bfU)\bfI+\rho\bfg \end{equation}

对于连续性方程,同样的,如果认为密度为常数,有连续性方程:

\begin{equation}\label{cont} \nabla\cdot\bfU=0 \end{equation}

把连续性方程\eqref{cont}带入到动量方程\eqref{MM}中有:

\begin{equation}\label{MM2} \rho\frac{\partial \mathbf{U}}{\partial t}+\rho\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla p+\nabla \cdot(\mu \nabla \mathbf{U})+\rho\bfg \end{equation}

将方程\eqref{boussi}代入到\eqref{MM2}中:

\begin{equation}\label{MM3} \rho\frac{\partial \mathbf{U}}{\partial t}+\rho\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla p+\nabla \cdot(\mu \nabla \mathbf{U})+\left(\rho_\rref\left(1-\beta\left(T-T_\rref\right)\right)\right)\bfg \end{equation}

令:

\begin{equation} \rho_k=\frac{\left(\rho_\rref\left(1-\beta\left(T-T_\rref\right)\right)\right)}{\rho} \end{equation}

有:

\begin{equation}\label{MM4} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U} \otimes\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nu \nabla \mathbf{U})+\rho_k \bfg \end{equation}

方程\eqref{MM4}即为考虑Boussinesq假定的动量方程。可以看出调用Boussinesq假定的动量方程为不可压缩的。在调用Boussinesq假定的情况下,如果温差较低的情况下(如空气的15°温差内),误差将在1%以下。如果温差过大,调用Boussinesq假定可能会带来较大的误差。

更新历史
2018.12.12大修

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