chtMultiRegionFoam解析


李子丰


1. 引言

chtMultiRegionFoam是一个多计算域求解器。可用于同时求解流体换热以及固体导热。类似的应用非常广泛,比如电子元件的散热,管道系统的换热等。其他传热类求解器仅仅可以用于求解流体域的换热。chtMultiRegionFoam可以同时将固体域与流体域的热交换耦合起来,同时考虑固体与流体间的换热效果对散热的影响,使温度场满足热力学定律;此求解器解决的这种问题被称为共轭传热问题。chtMultiRegionFoam类似于其他的瞬态求解器,我们需要对系统中每个变量的方程依次求解,求解思路大致为:首先使用前一次迭代的固体温度求解流体方程以定义流体温度的边界条件,而后使用前一次迭代的流体温度求解固体方程,以定义固体温度的边界条件;运行这个迭代过程一直执行到收敛。

2. 控制方程与算法

我们分为三个部分来探讨chtMultiRegionFoam这一求解器

1.流体部分控制方程(solveFluid.H)

2.固体部分控制方程(solveSolid.H)

3.温度传输的流固耦合(turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C)

2.1 流体控制方程

我们此处不考虑反应故省略物料输运方程(YEqn.H)。流体中的质量方程,动量方程,以及能量守恒方程与rhoPimpleFoam类似,如下所示:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \right)=0 \label{NS1} \end{equation}

\begin{equation} \frac{\partial \rho \mathbf{U} }{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \mathbf{U} \right)-\nabla \cdot \tau=-\nabla p \label{NS2} \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}

附加状态方程

\begin{equation} \rho=\psi p \label{bar} \end{equation}

其中的$\mathbf{U}$为速矢量,$p$为压力,$\rho$为密度,$\tau$为剪切应力,$h$为比焓,$K$表示机械能,$\alpha_{\mathrm{eff}}$表示有效导热系数,其等于$\frac{\nu_t}{\mathrm{Pr}}$,$\psi=1/RT$表示可压缩性。对于能量方程的的讨论详见CFD中的能量方程。首先对方程\eqref{NS2}中的时间项进行对速度$\bfU$关于时间t的欧拉全隐离散有:

\begin{equation} \int \int {\frac{{\partial\rho \mathbf{U}}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {\rho_\rP^{t+\Delta t}\mathbf{U}_\mathrm{P}^* -\rho_\rP^{t} \mathbf{U}_\mathrm{P}^t} \right)\;} {V_\mathrm{P}}, \label{mom1} \end{equation}

其中$\bfU_\rP^*$为待求速度,但方程\eqref{mom1}中继续存在未知量$\rho_\rP^{t+\dt}$,因此,在离散动量方程之前,需要对密度方程\eqref{NS1}进行更新获得$\rho_\rP^{t+\dt} $。但由于方程\eqref{NS1}中的速度为当前的速度,因此求解的密度为显性解$\rho_\rP^{t+\dt}$。随后,对流项隐性离散:

\begin{equation} \int \int {\nabla \cdot \left( {\rho\mathbf{U}\mathbf{U}} \right)\mathrm{d}V\mathrm{d}t = \int \int {\rho\mathbf{U}\mathbf{U}\cdot\mathrm{d}\bfS\mathrm{d}t} = } \sum {{{\left( {{\rho^t\mathbf{U}^*}{\mathbf{U}^t}} \right)}_f}} \cdot\bfS_f\Delta t = \Delta t\sum {{F_f^t \mathbf{U}_f^*} } , \label{mom2} \end{equation} 其中的$F_f^t$表示质量通量,其等于$\rho_f\bfU_f\cdot\bfS_f$。方程\eqref{NS2}中的$\nabla\cdot\tau$ 通过湍流模型进行模化,在此不做离散处理。压力项显性离散: \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^t\bfS_f\right), \label{gradp} \end{equation} 将方程\eqref{mom1}、\eqref{mom2}以及\eqref{gradp}整理有(此处未考虑粘性拉普拉斯项): \begin{equation} \left( {\rho_\rP^{t+\dt}\mathbf{U}_\mathrm{P}^* -\rho_\rP^t \mathbf{U}_\mathrm{P}^t} \right)\; {V_\mathrm{P}} + \Delta t\sum F_f^t \mathbf{U}_f^* =-\Delta t\sum \left(p_f^t\bfS_f\right) \label{predic} \end{equation} 方程\eqref{predic}中$\mathbf{U}_f^*$被定义为面上的预测速度。面速度需要从体心速度进行插值来获得,在此步可以引入各种插值格式。假设在均一网格上使用中心线性格式: \begin{equation} \mathbf{U}_f^* = \frac{\bfU_\rP^* + \bfU_\rN^*}{2}. \label{linear} \end{equation} \begin{equation} p_f^t = \frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}. \label{linear2} \end{equation} 将方程\eqref{linear},\eqref{linear2}代入到方程\eqref{predic}有 \begin{equation} \left( {\rho_\rP^{t+\dt}\mathbf{U}_\mathrm{P}^* -\rho_\rP^t \mathbf{U}_\mathrm{P}^t} \right)\; {V_\mathrm{P}} + \Delta t\sum F_f^t \frac{\bfU_\rP^* + \bfU_\rN^*}{2} =-\Delta t\sum \left(\frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}\bfS_f\right) \end{equation} 整理有: \begin{equation} \left(\frac{\rho_\rP^{t+\dt}}{\Delta t}+\frac{\sum F_f^t}{2V_\rP} \right)\bfU_\rP^* + \frac{\sum F_f^t}{2V_\rP}\bfU_\rN^* = \frac{\rho_\rP^{t}}{\Delta t}\bfU_\rP^t-\frac{1}{V_\rP}\sum \left(\frac{p_\mathrm{P}^t + p_\mathrm{N}^t}{2}\bfS_f\right) \label{predic2} \end{equation} 将上述方程简写为 \begin{equation} {A_\mathrm{P}}\mathbf{U}_\mathrm{P}^*{\rm{ + }}\sum {A_\mathrm{N}\mathbf{U}_\mathrm{N}^*} = S_\mathrm{P}^t-\frac{1}{V_\rP}\sum \frac{{p_\mathrm{P}^t + p_\mathrm{N}^t}}{2}\bfS_f, \label{apanmom} \end{equation} 其中 \begin{equation}\label{Ap} A_\mathrm{P}=\frac{\rho_\rP^{t+\dt}}{\Delta t}+\frac{\sum F_f^t}{2V_\rP}, \end{equation} \begin{equation} A_\mathrm{N}= \frac{{F_f^t}}{2V_\rP} , \end{equation} \begin{equation} S_\mathrm{P}^t= \frac{\rho_\rP^{t}}{{\Delta t}}\mathbf{U}_\mathrm{P}^t. \end{equation} 上述方程中$\rho_\rP^{t+\dt}$为已知量,因为其已经通过密度方程获得显性解。求解方程\eqref{apanmom}即可获得速度$\bfU^*$。此处的速度并不满足质量守恒。在求解预测速度$\bfU^*之$后,我们将其用于求解能量守恒方程(焓方程),由于能量守恒方程不参与下文质量守恒方程所约束的压力速度耦合过程,求解焓与温度等过程此处进行省略。

2.2 压力泊松方程
压力基求解器需要利用连续性方程构建压力泊松方程,进而修正速度。在获得速度$\bfU^*$后,利用$A_\rP$与$A_\rN$可以组建$\bfHbyA^*$速度收敛量(这里并未考虑压力的作用): \begin{equation} \bfHbyA_\rP^{*} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}^{*}} + S_\rP^t} \right) \label{hbya2} \end{equation} 同时有: \begin{equation} \bfU_\rP^{*}=\bfHbyA_\rP^{*} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t}\bfS_f \label{Up2} \end{equation} \begin{equation} \bfU_f^{*}=\bfHbyA_f^{*} - \frac{1}{{{A_{\mathrm{P},f}}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t}\bfS_f\right)_f \label{Up3} \end{equation} 现考虑连续性方程,其最终的离散形式为: \begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f (\rho_f^t+\rho_f')(\mathbf{U}_f^{*}+\bfU_f') \cdot \bfS_f=0. \label{cont0} \end{equation} 参考rhoSimpleFoam解析,其最终可以写为 \begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \left(\rho_f^t\bfU_f^{*} + \rho_f^n\bfU_f' + \frac{1}{RT}p_f'\bfU_f^{*} \right)\cdot\bfS_f=0. \label{cont} \end{equation} 略去邻点影响有: \begin{equation}\label{uFf} \bfU_f'=-\frac{1}{V_\rP}\frac{1}{{{A_f^{n}}}} \left(\sum_f p_f’\bfS_f\right)_f \end{equation} \begin{equation} \bfU_\rP^{t+\dt}=\bfHbyA_\rP^{*} - \frac{1}{{{A_\mathrm{P}}}} \frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f \label{Up22} \end{equation} 将方程\eqref{uFf}、\eqref{Up3}同时代入至\eqref{cont}有: \begin{multline} \frac{\rho^{t+\dt}-\rho^t}{\dt} \\ +\sum_f \left(\rho_f^t\left(\bfHbyA_f^{*} - \frac{1}{{{A_{\mathrm{P},f}}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t}\bfS_f\right)_f\right) - \rho_f^n\frac{1}{V_\rP}\frac{1}{{{A_f^{n}}}} \left(\sum_f p_f’\bfS_f\right)_f + \frac{1}{RT}p_f'\bfU_f^{*} \right)\cdot\bfS_f \\ =0. \label{cont2} \end{multline} 引入修正压力: \begin{equation} p_f^{t+\dt}=p_f^t+p_f' \end{equation} 整理方程\eqref{cont2}有 \begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^n \mathbf{HbyA}^{*}_f\cdot\bfS_f + \sum_f \frac{1}{RT}p_f'\bfU_f^{*}\cdot\bfS_f = \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \right)\cdot\bfS_f \label{press} \end{equation} 由于 \begin{equation} \frac{1}{RT}=\frac{\gamma}{c^2}=\frac{\p\rho}{\p p}. \label{RT} \end{equation} 在马赫数较小的情况下,c趋向于无穷大,RT项可以省略。有: \begin{equation} \frac{\rho^{t+\dt}-\rho^t}{\dt}+\sum_f \rho_f^n \mathbf{HbyA}^{*}_f\cdot\bfS_f = \sum_f \left(\frac{\rho_f^n}{{{A_f^n}}} \left(\frac{1}{V_\rP}\sum_f p_f^{t+\dt}\bfS_f\right)_f \right)\cdot\bfS_f \label{pres} \end{equation} 其即为亚音速情况下连续形式的压力泊松方程: \begin{equation} \frac{\p\rho}{\p t}+\nabla \cdot \left(\rho\mathbf{HbyA} \right)=\nabla \cdot \left( \frac{\rho}{{{A}}}\nabla p \right) \label{p2pre} \end{equation} 求解器中的迭代过程可以表示为下面几个步骤:

1.依据初始速度组建通量,显性离散方程\eqref{NS1},获得下一个时间步的密度$\rho^{t+\dt}$;

2.通过方程\eqref{apanmom},获得$\bfU^*$以及$\mathbf{HbyA}^*$;

3.求解方程\eqref{pres},获得压力$p^{t+\dt}$,注意,$p^{t+\dt}$并不是精准的(忽略了临点影响);

4.通过方程\eqref{Up22}获得速度$\bfU^{t+\dt}$以及通量,再次通过方程\eqref{NS1}更新密度($\rho^{t+\dt}$并不是精准的);

5.依据状态方程,以及压力$p^{t+\dt}$,更新$\rho^{t+\dt}$,通过$\rho^*$与$\rho^{t+\dt}$,求得连续性误差;

6.回到第一步迭代求解几次,在时间步内收敛的情况下,$\rho^* \rightarrow \rho^{t+\dt}$

对比Piso与Pimple的迭代逻辑的区别在于:Pimple逻辑中:while (pimple.1oop())在时间步中增加了压力速度耦合的求解,即在Piso逻辑外增加了循环,有利于大时间步长瞬态计算的收敛,(例如在大计算域范围中的共轭传热计算),注:在fvSolution中指定nCorrectors为1即为Piso逻辑。

2.3 固体部分控制方程

在固体区域由于不涉及任何流动,故只有热传导过程,在本文中我们忽略固体内空隙的热传递及热辐射过程,有固体的热传导方程为:

\begin{equation} \frac{\partial}{\partial t}(\rho h)=\nabla \alpha_{e f f} \nabla h \end{equation}

2.4 温度传输的流固耦合

温度传输的流固耦合(共轭传热)本质上来讲是传热学(杨世铭, 陶文铨. 《传热学(第4版)》;P45(2-17))之中的第三类边界条件;这种边界条件在OpenFoam-9中由src中turbulentTemperatureCoupledBaffleMixed文件夹下的turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.C定义


固体和气体区域界面上的两个单元格示意图


如图所示,在有限体积法离散的网格单元中,我们假设左边网格1为流体网格,右边网格2为固体网格, $T_{c}$为储存在网格中心的温度, $ T_{p} $为两相邻网格共享面上的温度, $ q_{1}^{\prime \prime} $为流出网格1的热通量 ,$ q_{2}^{\prime \prime} $为流入网格2的热通量。由第三类边界条件可得到如下关系:

\begin{equation} T_{p, 1}=T_{p, 2}=T_{p}, \end{equation} \begin{equation} q_{1}^{\prime \prime}=q_{2}^{\prime \prime}=q^{\prime \prime}. \end{equation}

由傅里叶定律可得:

\begin{equation} k_{1} \Delta_{1}\left(T_{c, 1}-T_{p}\right)=k_{2} \Delta_{2}\left(T_{p}-T_{c, 2}\right) . \end{equation}

其中 $ k $为传热系数,进而我们可以得到网格边界上的温度与热通量如下所示:

\begin{equation} T_{p}=\frac{k_{1} \Delta_{1} T_{c, 1}+k_{2} \Delta_{2} T_{c, 2}}{k_{1} \Delta_{1}+k_{2} \Delta_{2}} \end{equation}

\begin{equation} q^{\prime \prime}=k_{1} \Delta_{1}\left(T_{c, 1}-T_{p}\right)=k_{2} \Delta_{2}\left(T_{p}-T_{c, 2}\right) \end{equation}

3. 验证算例

本算例为OKS课程基本算例之一,算例来源于某电路板散热项目,基于此项目的原型进行简化。算例如下图所示:


在生成网格的过程中,可以首先生成整个计算域的网格,然后通过topoSet定义固体区域的网格,剩余的部分即为流体区域的网格。在通过blockMesh生成整个计算域网格之后,需要修改system/topoSetDict中的参数定义固体区域,如下所示:

actions
(
    {
        name    cs;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (1 -0.5 -0.1) (8 0.5 0.56);
        }
    }
    {
        name    solid;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo { set cs; }
    }
);

同时,需要在constant文件夹中,通过regionProperties文件来确定我们所定义的区域

regions
(
   fluid (fluid)//括号中为0文件夹中我们定义的流体名称,同时fliud子区域的边界将获得与其所属整体边界相同的名称
   solid (solid)//括号中为0文件夹中我们定义的固体名称,同时solid子区域的边界将获得与其所属整体边界相同的名称
);

初始条件设置:首先设置液体区域的初始场,其位置为/0/fluid

  • 初始场p:内部即均一分布uniform 1e5
  • 初始场速度U:入口速度为0.2,即均一分布uniform (0.2 0 0)
  • 流固耦合区域p:wall为calculated条件;
  • 流固耦合区域U:wall为noSlip无滑移边界;
  • 流固耦合区域T:整体温度场设置为uniform 297由上文所述,wall的类型设置为compressible类型下的turbulentTemperatureCoupledBaffleMixed
  • wall
    {
        type            compressible::turbulentTemperatureCoupledBaffleMixed;
        kappa           kappa;
        Tnbr            T;
        value           $internalField;
    }
    

    设置固体区域的初始场,其位置为:/0/solid

  • 由于固体部分只涉及导热计算,设置温度场T整体为:uniform 297,wall为compressible::turbulentTemperatureCoupledBaffleMixed;类型;并设置其底部区域topAndBottom拥有初始温度uniform 313进而在固体内产生温度梯度。
  • 流体与固体的物性参数在constant文件夹下指定:

  • 在constant/fluid/momentumTransport中粘度模型选择层流,即laminar
  • 在constant/fluid/thermophysicalProperties中我们指定流体物性参数如下;
  • mixture
    {
        specie
        {
            nMoles      1;
            molWeight   18;         // [g/mol]
        }
        equationOfState
        {
            rho         998.19;     // [kg/m^3]
        }
        thermodynamics
        {
            Cv          4150;       // [J/kg/K]
            Hf          0;
        }
        transport
        {
            mu          1.0005e-3;  // [kg/m/s]
            Pr          6.31;
        }
    }
  • 在constant/solid/thermophysicalProperties中我们指定固体物性参数如下:
  • mixture
    {
        specie
        {
            nMoles      1;
            molWeight   63.5;   // [g/mol]
        }
    
        transport
        {
            kappa   380;        // [W/m/K]
        }
    
        thermodynamics
        {
            Hf      0;
            Cv      385;        // [J/kg/K]
        }
    
        equationOfState
        {
            rho     8940;       // [kg/m^3]
        }
    }
    

    物性参数的各个符号及计算公式详见OpenFOAM用户指南-热物理模型。算例可以通过输入./Allrun直接得到计算结果。可以看到固体内部的温度根据导热定律传递,而到达固液界面时,根据第三类边界条件进行固体与液体之间的换热。同时液体区域的问题也会影响固体区域的问题。固体的高温也会导致液体的问题上升。下面这个动图即为本算例模拟的动态过程:


    点击下载电子元件传热算例

    东岳流体 2014 - 2021
    勘误、讨论、补充内容请前往CFD中文网