返回主页

CFD中的能量方程


总能量方程

能量守恒定义为绝热系统的总能量是一个常数。即总的能量不随时间变化,只能从一种形式转换为另一种形式且不能凭空消失。在CFD里面通常只考虑动能(机械能)分子内能(内能)。能量守恒可以表示为$^1$:

流体微团内能量的变化率=流入流体微团的净热流量+体积力和表面力对流体微团做功的功率

单位质量的动能定义为$^1$:
\begin{equation} K=\frac{1}{2}|\mathbf{U}|^2 \label{K} \end{equation}
其中$\mathbf{U}$为速度。单位质量的内能定义为:$e$。那么流体微团动能和内能总和的时间变化率定义为:
\begin{equation} \rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{A} \end{equation}
方程($\ref{A}$)即为"流体微团内能量的变化率"。

流入流体微团的净热流量来自于加热、辐射、热传导等。我们定义$r$为产生的热源。那么流体微团的净热源即为:
\begin{equation} \rho r \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{heatSource} \end{equation}
另外,由于热传导(热流具有方向性)对流体微团的加热为(热通量矢量定义为$\mathbf{q}$):
\begin{equation} - \nabla \cdot \mathbf{q} \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{heatFlux} \end{equation}
进一步的依据傅里叶定律有:
\begin{equation} \mathbf{q}=-k \nabla T \label{fourier} \end{equation}
把方程($\ref{heatFlux}$)和方程($\ref{heatSource}$)加和,因此流入流体微团的净热流量为:
\begin{equation} \rho r \mathrm{d}x\mathrm{d}y\mathrm{d}z - \nabla \cdot \mathbf{q}\mathrm{d}x\mathrm{d}y\mathrm{d}z \label{B} \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}
压力和剪切力对流体微团做功的功率可以表示为$^2$:
\begin{equation} \left( -\nabla\cdot(p\mathbf{U})+ \nabla \cdot(\mathbf{\tau} \cdot \mathbf{U}) \right) \mathrm{d}x\mathrm{d}y\mathrm{d}z \label{C} \end{equation}
结合方程($\ref{A}$),($\ref{B}$),($\ref{g}$),($\ref{C}$)我们有最终的总能量方程:
\begin{equation} \rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} = \rho r - \nabla \cdot \mathbf{q} + \rho \mathbf{g} \cdot \mathbf{U}-\nabla\cdot(p\mathbf{U})+\nabla \cdot(\tau \cdot \mathbf{U}) \label{EEqn} \end{equation}
其中的$\nabla \cdot(\tau \cdot \mathbf{U})=\left(\nabla \cdot \tau \right)\cdot \mathbf{U}+\tau:\nabla \mathbf{U}$。可以看出,方程($\ref{EEqn}$)采用的是非守恒形式。

内能方程

虽然方程($\ref{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} \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}
将方程($\ref{EEqn}$)中减去方程($\ref{K2}$)我们有最终的内能方程:
\begin{equation} \frac{\mathrm{D}e}{\mathrm{D}t} =-p\nabla \cdot \mathbf{U}-\nabla\cdot \mathbf{q} + \tau : \nabla \mathbf{U} + \rho r \label{e} \end{equation}
可以看出,内能方程中不含有体积力项。再一次的,方程($\ref{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}
回到非守恒形式的方程($\ref{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 \mathbf{q} + \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{EEqnconser} \end{equation}
把方程($\ref{hDef}$)代入到方程($\ref{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 \mathbf{q} + \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 \mathbf{q} + \rho \mathbf{g} \cdot \mathbf{U}+\nabla \cdot(\tau \cdot \mathbf{U}) \label{h0Eqn} \end{equation}
需要提及的是,在OpenFOAM的新版本中,传热、可压缩、化学反应求解器求解的主要为守恒形式的总能量方程和比焓方程$^3$。并且,方程($\ref{EEqnconser}$)和方程($\ref{hEqn}$)中的$\nabla \cdot(\tau \cdot \mathbf{U})$以及$\rho \mathbf{g} \cdot \mathbf{U}$被忽略。但是在rhoCentralFoam中,则植入了此项。OpenFOAM官方表示,能量方程的选择在某些情况下是至关紧要的$^3$。比如在进行震波捕获计算的时候,求解总能量方程要比求解内能方程结果精确的多。OpenFOAM中还进一步的假设$\mathbf{q}=-\alpha_\mathrm{eff}\nabla e$或$\mathbf{q}=-\alpha_\mathrm{eff}\nabla h$,这样方程($\ref{EEqnconser}$)以及方程($\ref{hEqn}$)则简化为(若无热源$\rho 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}
能量守恒相对于动量守恒要复杂一些。在不可压缩流动中,唯一重要的是内能方程。在考虑传热的情况下,内能要比热能小得多$^4$。只要温度对流体的影响不是很重要,那么能量方程可以放在压力速度耦合之后进行求解。在这种情况下是一种单向耦合。同时,需要注意的是内能方程实际上来源于动量方程。但是总能量方程以及焓方程却可以单独推导。在某些情况下求解内能方程可能会带来一些问题$^4$。同时,对于不可压缩流体,能量的增加很少见,能量的耗散却比较明显。对于可压缩流,双方都比较重要。

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$表示流体的体膨胀系数。考虑可压缩并附加重力的NS方程有:
\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假定的动量方程为不可压缩的。

参考文献

[1]. 计算流体力学基础及其应用
[2]. An introduction to computational fluid dynamics: the finite volume method
[3]. http://cfd.direct/openfoam/energy-equation
[4]. Computational Methods for Fluid Dynamics

2018.01.09:添加内能部分的描述
2017.07.17:添加boussinesq假定
2016.12.13:依据反馈添加方程\eqref{h0Eqn}
2016.11.16:统一公式符号类型,添加方程\eqref{Efinal}以及方程\eqref{hfinal}


东岳流体dyfluid.com
勘误、讨论、补充内容请前往CFD中国