• 返回主页
  • CFD中的能量方程

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

    总能量方程

    能量守恒定义为绝热系统的总能量是一个常数。即总的能量不随时间变化,只能从一种形式转换为另一种形式且不能凭空消失。在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}

    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]. 安德森, 吴颂平, 刘赵淼. 计算流体力学基础及其应用[M]. 机械工业出版社, 2007.
    [2]. Versteeg H.K, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method[M]. Pearson Education, 2007.
    [3]. Energy Equation in OpenFOAM. http://cfd.direct/openfoam/energy-equation


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