CN105512433B - 流体-固体的节点化两相流建模方法 - Google Patents

Abstract

一种求解重力驱动无相变两相流流固耦合传热的流体‑固体节点化两相流建模方法。其特征在于忽略对流场的准确计算，用节点热网络思想处理两相流体。具体包括：A)将液相、气相抽象为具有一定质量、体积、温度的节点；B)使用重力方向体积积分计算两相分界面；C)使用一个表面辐射换热模型计算流固耦合壁面辐射热流密度；D)根据流固耦合壁面对流换热热流密度计算气相、液相温度；E)根据流固耦合壁面总热流密度计算固体瞬态导热。针对典型燃油箱内两相流固耦合传热过程，利用本发明的建模方法与Fluent商业软件相场方法计算结果相比，更好的反映了相分界面对传热的影响，计算效率提高一个数量级，满足工程计算需求。

Description

流体-固体的节点化两相流建模方法

技术领域

本发明涉及流体-固体的节点化两相流建模方法。

背景技术

目前，重力驱动简单两相流固耦合传热问题在工程上主要有两种求解方法。

第一种方法用基于扩散分界面的“相场”方法来求解两相流区域的流场，并与固体实时耦合。主要应用于商业CFD(计算流体力学)软件，如FLUENT。“相场”方法用序参数φ_a(x,t)、φ_b(x,t)对两相a、b的空间分布进行建模。在相a占据的区域φ_a(x,t)等于1，在相a不存在的区域φ_a(x,t)等于0，在相a周围的扩散分界面内φ_a(x,t)在1到0之间连续变化。在流场中的任何位置，φ_a(x,t)与φ_b(x,t)之和都等于1。在扩散分界面内，流体的特征参数通过对两相各自的特征参数依据序参数进行插值来获得。所以，流场中流体的特征参数处处连续，在整个求解区域内各相的控制方程是统一的。这种方法适用范围较广，可以处理比较复杂的相变两相流问题，但计算效率较低，而且处理不互溶的简单两相流时有一定误差。第二种方法用节点化的方法对流体、固体进行建模，求解由这些热节点组成的节点热网络。主要应用于商业节点热网络软件，如Flowmaster。在工程应用中，这种算法将固体结构、两相流体抽象成若干个具有一定温度、质量的节点。考虑各节点间的对流、辐射、导热作用，以及环境辐射、油箱内热源，就可以对各节点列出热平衡方程。对整个系统的热平衡微分方程组进行瞬态迭代求解，就可以实现对整个流固耦合传热系统的仿真、求解。虽然这种方法具有高的计算效率，但计算结果的质量较差，无法给出准确的温度场。

发明内容

根据本发明的一个方面，提供了一种流体-固体的节点化两相流建模方法，其特征在于将两相流体中的液相和气相分别抽象为具有温度、质量、体积的节点，并包括：

A)计算当前时刻两相流体中的液相和气相的体积，通过重力方向的体积积分计算出两相分界面；

B)使用对流换热公式，计算流固耦合壁面的对流换热热流密度；

C)使用一个表面辐射换热模型，计算流固耦合壁面的辐射热流密度，并将辐射热流密度与步骤B)中求出的对流换热热流密度相加，获得流固耦合壁面的总热流密度；

D)计算下一时间步内气相和液相获得的总对流换热热量，分别更新气相和液相的加权平均温度；

E)根据流固耦合壁面的总热流密度，计算流固耦合壁面的温度梯度，进行一个时间步的固体瞬态传热计算；

F)下一时间步之后的时刻，重复A)至E)，直到到达瞬态传热终止时刻。

附图说明

图1是根据本发明的一个实施例的流体-固体的节点化两相流建模方法的流程图。

图2是一个具备隔热层和内置电子设备的典型飞机燃油箱模型的纵剖面示意图，用于检验本发明的方法的可靠性。

图3A显示了用根据本发明的流体-固体的节点化两相流建模方法对图2所示的飞机燃油箱模型进行的数值计算的结果；图3B显示了用Fluent软件对图2所示的飞机燃油箱模型进行的计算结果。

具体实施方式

根据本发明的流体-固体的节点化两相流建模方法，将两相流体抽象成具有一定温度、质量、体积的节点，用节点热网络的思想来处理，规避多相流CFD计算。同时，依然用数值传热学的方法来求解固体区域的瞬态传热，保证固体区域温度的计算结果质量。由于计算中将流体节点化处理，所以辐射模型选择不依托流体区域的表面辐射换热模型(S2S)。

如图1所示，根据本发明的一个实施例的流体-固体的节点化两相流建模方法包括：

1)将液相、气相分别抽象为具有温度、质量、体积的“节点”；

2)计算当前时刻两相流体中气相和液相的体积，通过重力方向的体积积分计算出液面高度(两相分界面)；(步骤101、102)

3)使用对流换热公式，计算流固耦合壁面对流换热热流密度；(步骤103)

4)使用表面辐射换热模型(S2S)模型，计算流固耦合壁面的辐射热流密度(步骤104)，并将辐射热流密度与步骤103中求出的对流换热热流密度相加，获得总热流密度；(步骤105)

5)计算下一时间步内气相和液相获得的总对流换热热量，分别更新气相和液相的加权平均温度；(步骤106、107)

6)根据流固耦合壁面的总热流密度计算流固耦合壁面温度梯度，进行一个时间步的固体瞬态传热计算；(步骤108、109)

7)在下一时刻(下一时间步之后的时刻)，重复2)至6)，直到到达瞬态传热终止时刻。

2.1两相区域划分

根据本发明的一个实施例，在如图1中的步骤101、102所表示的处理中，为了确定流固耦合壁面上的对流换热系数，需要对两相流体区域进行划分。

假设液相质量流速很小且液面的波动可以忽略，则两相区域划分可以用液相的液面高度来代替。即首先根据进/出口质量流量获得每一时刻液相体积，然后根据液相体积确定液面高度，即两相分界面。

考虑到固体内几何模型往往比较复杂，液面高度可根据液相体积通过网格单元体积分的方法更简单、更准确地获得。具体方法包括：

1)获得所有体单元的中心点坐标和体积，并根据中心点高度由低到高进行排序。

2)设定初始液面高度，将该高度以下所有体元的体积相加，即可近似得到初始液相体积及质量。

3)根据进/出口质量流量获得各时刻液相质量，并结合液相平均温度确定每一时刻的液相体积。

4)按中心点高度从低到高依次将对应体单元的体积相加，直到第i体元时体积之和首次大于等于该时刻液相体积，则将该体单元网格的中心点高度近似看作当前时刻液面高度，即两相分界面的位置。如果所有体单元网格尺寸均在毫米量级，则这种方法确定的液面高度满足工程设计误差范围。

2.2 S2S辐射模型

在本算法中，两相流区被节点化处理，所以在计算过程中流体网格不会被加载入求解器，需要使用不受流场影响的辐射模型(步骤104)。表面辐射换热模型(S2S)在瞬态计算开始之前，需要使用流体网格来计算辐射角系数，但在瞬态计算中不依赖流场，适合在本方法中使用。

使用S2S模型计算辐射热流时，所需消耗的时间与参与辐射的面单位数目的平方成正比。这意味着，如果以面元为单位计算热辐射，需要消耗大量的时间。所以，可以将若干个面元组合成一个表面束，然后以表面束为单位计算辐射热流。

对于表面束t，其本身的辐射为：

其中，λ_t为表面束t的辐射率，σ为波尔兹曼常数，E_b为黑体辐射力，T_t为表面束t的绝对温度。考虑到投射到表面的辐射G会有(1-λ_t)·G_t的部分被反射，表面束t的有效辐射J_t可以表示为：

J_t＝λ_t·E_b+(1-λ_t)·G_t

进一步的，有有效辐射方程为：

其中F_t,i为表面束t到表面束i得辐射角系数。因为：

所以有效辐射方程的系数矩阵严格对角占优，故可逆。通过求解该方程，可以获得所有表面束的有效辐射热流密度J。最终，表面束t吸收的辐射热流密度的表达式为：

其中Q_rt为表面束t包含的所有面元的辐射热流密度。将辐射热流密度与通过面元进入固体区域的对流换热热流密度相加，就可以得到通过面元进入固体的总热流密度。

由于以表面束为单位计算辐射角系数会造成一定误差，所以在实际计算中，需要在有效辐射方程、辐射热流密度表达式的基础上加以修正。

具体包括，在开始传热计算之前，做如下处理：

1)对参与辐射传热的燃油箱内部壁面进行结块处理，记录每个表面束包含的面元。

2)以表面束为单位，计算各表面束之间的辐射角系数。

3)求解有效辐射方程的系数矩阵的逆矩阵，用于计算各表面束的有效辐射。

在每一个时刻，需要通过如下处理，计算各面元辐射热流，包括：

1)计算表面束的加权平均温度。

2)计算表面束的有效辐射。

3)计算表面束的净辐射热流密度。

4)将表面束的辐射热流密度赋给表面束包含的面元。

2.2.1辐射方程修正

通过表面束操作，辐射面单位的数目大大的降低了。这使得辐射计算的速度大大 提高，但也导致了辐射角系数的计算存在一定误差。对于准封闭的空间(只有出入口与外界 连通)，绝大多数表面束应该满足

为此，将有效辐射方程修正为：

将净辐射热流密度表达式修正为：

对于的表面束，可以近似认为，每一个辐射角系数都存在的相对误差。

由于考虑辐射的流固耦合壁面不一定是严格封闭的，对于的表面束，包含两种情形。在第一种情形下，表面束发射的辐射全部被其它表面束所接收，这样的表面束也应满足可以近似认为，每一个辐射角系数都存在的相对误差。

在第二种情形下，表面束发射的部分辐射通过流体出入口等不封闭处进入了外部环境，这样的表面束应满足而角系数修正相当于将投射辐射修正为：

此时可以近似认为该表达式中的第一项为来自外部环境，通过不封闭处进入的辐射。

2.2.2辐射热流计算

在计算辐射热流时，需要首先获得表面束的加权平均温度。这个温度可以通过表面束包含的面元的温度T_k及面积S_k获得。可以表示成：

获得所有表面束的加权平均温度后，就可以获得各表面束的黑体辐射热流密度：

将该黑体辐射热流密度与有效辐射方程的系数矩阵的逆矩阵相乘，就可以获得各表面束的有效辐射热流密度J_t。最后将有效辐射热流密度代入修正后的净辐射热流密度表达式：

就可以获得通过表面束进入固体区的热流密度。根据表面束的定义，这个热流密度也就是表面束包含的所有面元上的辐射热流密度，将该热流密度赋给各面元即可。

2.3流体加权平均温度更新

根据步骤103中求得的流固耦合壁面对流换热热流密度q_i，分别计算气相、液相获得的总热量Q(步骤106)，如下式所示：

其中Δt为时间步长，s_i为面元面积。

根据进入气相和液相的热量，分别更新这两者的加权平均温度(步骤107)。在瞬态传热过程中，自然对流两相流体可能存在低速的强制流入、流出，而这也会对温度更新产生一定影响。

考虑两相流体的流入、流出，液相、气相分别为流入相和流出相。记当前时刻流出相的体积为V₁、平均温度为T₁，密度ρ₁，比热为C_p1，通过流固耦合壁面进入流出相的总热量为Q_1，

则记流出相在经过Δt时间后的温度为：

考虑流出相流出后，有温度为初始环境温度T₂₀的另一相流体流入，并与原有流体达到等温，记当前时刻流入相的体积为V₂、平均温度为T₂，密度为ρ₂，比热为C_p2，通过流固耦合壁面进入流入相的总热量为Q₂，流出相流出的体积流速为v₀，把流入相经过Δt时间后的温度表示为：

2.4固体瞬态导热计算

在步骤108中，根据步骤105中求得的流固耦合壁面总热流密度q_i0，计算流固耦合壁面面元上的温度梯度

其中，k为导热系数。

在固体区域，可使用通用的固体瞬态导热算法(步骤109)。控制方程形式如下：

其中：

ρ＝密度

k＝导热系数

T＝温度

S_h＝体积热源

控制方程左边第一项表示固体随时间的温度变化，右边两项分别表示传导引起的热流以及固体内部的体积热源。

在有限的时间间隔[t,t+Δt]内，对控制容积作积分，即可将控制方程离散为代数方程。在每个时间步内对代数方程进行迭代求解即可。

3实用算例

为了检验本发明的方法的可靠性，我们选择了具备隔热层和内置电子设备的典型飞机燃油箱模型。在该模型中，燃油箱暴露在机身表面的部位(包括顶面201、水平换热面202、垂直换热面203、底面204)包裹有隔热层(207)，该燃油箱内部有两个设备壳体(205、206)。图2为该燃油箱模型的纵剖面示意图。图中黑色部分为隔热层(207)，未包裹隔热层的区域为绝热壁面(208)，隔热层(207)外有1730K的对流热源(气动热，未显示)。燃油箱与隔热层(207)接触的4个部分，按空间位置由上到下依次命名为：顶面(201)、水平换热面(202)、垂直换热面(203)、底面(204)。燃油箱左侧表面和内部各有一个设备壳体，分别为第一设备(205)和第二设备(206)。位于燃油箱内的第二设备(206)与燃油箱通过肋片(未显示)连接。在初始时刻，燃油箱内充满燃油，固体和流体的温度皆为303K。随时间的推移，燃油匀速流出燃油箱，空出的体积则被空气充满。

本发明人对上述飞机燃油箱模型进行了数值计算，并与Fluent软件的计算结果进行了对比。计算结果如图3A和3B所示。比较两种方法的计算结果，两种方法的求得的1560s温升相对误差在20％以内。水平换热面(202)和垂直换热面(203)两油箱壁面的快速升温过程在根据本发明的方法的计算结果中更明显。

由于气体与液相的比热在同一量级，而密度却相差了3个量级，所以单位体积的气体升高1K所需的热量远小于单位体积的液相。因此，液相升温的速度要小于气体，而获得的热流却高于气体。这也导致了与液相接触的壁面输入气液混合区的热流较大，而自身温度较低。

另一方面，虽然两个设备与主要热源顶面(201)之间同时存在对流换热、辐射传热以及热传导，但经由固体区域(外壁、肋片)的热传导是其中效率最高的换热方式。在温度趋势图中统计的设备温度最高点温度，实际上是热源顶面(201)与设备壁面通过固体热传导换热的结果。

综合以上因素，燃油箱温度变化应符合如下特征：在360s之前，穿过隔热层的热流很小，油箱和其中的流体都没有明显温升。在360s之后，油箱和其中流体的开始明显升温，并且速度逐渐加快。而各条温度曲线也会逐渐分离。顶面(201)的内侧面一直与气体接触，所以顶面(201)的温度是最高的。气体温度以及两个设备上的最高温度紧随其后。水平换热面(202)和垂直换热面(203)的内侧面在初始状态下与液相区接触，而后随着液面的不断降低，逐渐与气体区接触。所以，水平换热面(202)与垂直换热面(203)先后出现温度快速升高的现象。在这两个壁面之后的，温度最低的是内侧面一直与液相区接触的底面(204)和液相。

从fluent软件的方法和根据本发明的节点化两相流建模方法的计算结果可以看出，用两种方法求得的温度曲线的趋势皆符合上述特征。由于根据本发明的方法使用的节点化两相流处理方法基于精确相分界面，Fluent软件使用的“相场”两相流处理方法则基于扩散分界面，所以本算法更好的捕捉到了这一由液面下降引起的快速温升过程。

在计算效率方面，根据本发明的节点化两相流建模方法计算一个工况的用时约为17.3小时，而fluent软件计算一个工况用时约为168小时。简化算法的计算效率相对商业软件提高约10倍。

总的来说，相对于相场CFD方法，本算法的计算效率可以提高至少1个量级，而且可以对不互溶的简单两相流固耦合传热问题获得精度较高的结果。相对于节点热网络法，本算法可以获得固体结构的温度场，更好的满足工程应用需求。

Claims (9)

1.流体-固体的节点化两相流建模方法，其特征在于将两相流体中的液相和气相分别抽象为具有温度、质量、体积的节点，并包括：

A)计算当前时刻两相流体中的液相和气相的体积，通过重力方向的体积积分计算出两相分界面；

B)使用对流换热公式，计算流固耦合壁面的对流换热热流密度；

C)使用一个表面辐射换热模型，计算流固耦合壁面的辐射热流密度，并将辐射热流密度与步骤B)中求出的对流换热热流密度相加，获得流固耦合壁面的总热流密度；

D)计算下一时间步内气相和液相获得的总对流换热热量，分别更新气相和液相的加权平均温度；

E)根据流固耦合壁面的总热流密度，计算流固耦合壁面的温度梯度，进行一个时间步的固体瞬态传热计算；

F)下一时间步之后的时刻，重复A)至E)，直到到达瞬态传热终止时刻。

2.根据权利要求1所述的方法，其特征在于用重力方向体积积分的方法计算简单两相流体的两相分界面包括：

对两相流体的区域进行划分，以确定固体与两相流体接触的壁面上的对流换热系数，

忽略液面的波动，用液相的液面高度代替两相区域划分，即首先根据进/出口质量流量获得每一时刻液相体积，然后根据液相体积确定液面高度，即两相分界面，

根据液相体积，通过网格单元体积分的方法获得液面高度，具体方法包括：

A1)获得所有流体单元的中心点坐标和体积，并根据中心点高度由低到高进行排序，

A2)设定初始液面高度，将该高度以下所有体元的体积相加，从而得到初始液相体积及质量，

A3)根据进/出口质量流量获得各时刻液相质量，并结合液相平均温度确定各时刻的液相体积，

A4)按中心点高度从低到高依次将对应体单元的体积相加，直到第i体元时体积之和首次大于等于该时刻液相体积，则将该体单元网格的中心点高度看作当前时刻液面高度，即两相分界面的位置。

3.根据权利要求1所述的方法，其特征在于使用一个表面辐射换热模型计算流固耦合壁面辐射热流密度包括：

将若干个面元组合成一个表面束，然后以表面束为单位计算辐射热流密度，包括：

把表面束t本身的辐射表示为：

λ_t·E_b＝λ_t·σ·T_t ⁴ (1)

其中，λ_t为表面束t的辐射率，σ为波尔兹曼常数，E_b为黑体辐射力，T_t为表面束t的绝对温度，考虑到投射到表面的辐射G会有(1-λ_t)·G_t的部分被反射，

把表面束t的有效辐射J_t表示为：

J_t＝λ_t·E_b+(1-λ_t)·G_t (2)

采用有效辐射方程：

其中G_t为投射到表面束t的辐射，J_i为投射到表面束i的有效辐射，F_t,i为表面束t到表面束i的辐射角系数，因为：

所以有效辐射方程的系数矩阵严格对角占优，

通过求解方程(3)，获得所有表面束的有效辐射热流密度J，

把表面束t吸收的辐射热流密度的表示为：

其中Q_rt为表面束t包含的所有面元的辐射热流密度，将辐射热流密度与通过面元进入固体区域的对流换热热流密度相加，就可以得到通过面元进入固体的总热流密度。

4.根据权利要求3所述的方法，其特征在于进一步包括在开始瞬态传热计算之前，做如下处理：

对考虑辐射的流固耦合壁面进行结块处理，记录每个表面束包含的面元，

以表面束为单位，计算各表面束之间的辐射角系数，

求解有效辐射方程(3)的系数矩阵的逆矩阵，用于计算各表面束的有效辐射，

从而避免在瞬态传热计算过程中，对流场及流体网格的依赖。

5.根据权利要求3-4之一所述的方法，其特征在于进一步包括：

将有效辐射方程修正为：

将净辐射热流密度表达式修正为：

对于的表面束，认为每一个辐射角系数都存在的相对误差，

由于考虑辐射的流固耦合壁面不一定是严格封闭的，对于的表面束，包含两种情形，在第一种情形下，

在表面束发射的辐射全部被其它表面束所接收的情况下，认为每一个辐射角系数都存在的相对误差，

在表面束发射的部分辐射通过不封闭处进入了外部环境的情况下，将投射辐射修正为：

其中认为表达式(8)中的第一项为来自外部环境，通过不封闭处进入的辐射。

6.根据权利要求3所述的方法，其特征在于进一步包括在瞬态传热计算中的每一个时刻，需要通过如下处理，计算各面元辐射热流密度，包括：

在计算辐射热流密度时，首先获得表面束的加权平均温度，该加权平均温度表示成：

其中，T_k是表面束包含的面元的温度，S_k是表面束包含的面元的面积，

通过获得所有表面束的加权平均温度，获得各表面束的黑体辐射热流密度：

E_b＝σ·T_t ⁴ (10)

将该黑体辐射热流密度与有效辐射方程的系数矩阵的逆矩阵相乘，获得各表面束的有效辐射热流密度J_t，

将有效辐射热流密度代入修正后的净辐射热流密度表达式(7)，从而获得通过表面束进入固体区的辐射热流密度，该热流密度就是表面束包含的所有面元上的热流密度。

7.根据权利要求1所述的方法，其特征在于所述步骤D)包括：

根据步骤B)中求得的流固耦合壁面对流换热热流密度q_i，分别计算气相、液相获得的总热量Q，如：

其中Δt为时间步长，s_i为面元面积，

根据进入气相和液相的热量，分别更新气相和液相的加权平均温度，

考虑两相流体的流入、流出，液相、气相分别为流入相和流出相，记当前时刻流出相的体积为V₁、平均温度为T₁，密度为ρ₁，比热为C_p1，通过流固耦合壁面进入流出相的总热量为Q₁，

则流出相在经过Δt时间后的温度为：

考虑流出相流出后，有温度为初始环境温度T₂₀的另一相流体流入，并与原有流体达到等温，记当前时刻流入相的体积为V₂、平均温度为T₂，密度为ρ₂，比热为C_p2，通过流固耦合壁面进入流入相的总热量为Q₂，流出相流出的体积流速为v₀，把流入相经过Δt时间后的温度表示为：

8.根据权利要求7所述的方法，其特征在于所述步骤E)包括：

根据步骤C)中求得的流固耦合壁面总热流密度q_i0，计算流固耦合壁面面元上的温度梯度：

其中k为导热系数，

采用通用的固体瞬态导热算法，其中控制方程为：

其中：ρ＝密度，显焓k＝导热系数，T＝温度，S_h＝体积热源，T_ref是基准温度，C_p为定压比热容，控制方程(15)左边第一项表示固体随时间的温度变化，右边两项分别表示传导引起的热流以及固体内部的体积热源。

9.根据权利要求8所述的方法，其特征在于包括：

在有限的时间间隔[t,t+Δt]内，对控制容积作积分，从而将控制方程离散为代数方程，

在每个时间步内对代数方程进行迭代求解。