CN105677994B - 流体－固体耦合传热的松耦合建模方法 - Google Patents

Abstract

一种流体－固体耦合传热的松耦合建模方法，其特征在于忽略流场的瞬态变化过程，假设全局瞬态传热过程在“准稳态”流场下进行。具体计算流程包括：A)更新流场：单独以流体作为求解对象，将流固耦合壁面设为流体的固定温度边界，用稳态CFD算法求解流场；B)计算瞬态传热：同时以流体和固体为求解对象，将流固耦合壁面设为传热耦合边界，关闭流体的动量方程和湍流方程，计算瞬态传热直到下次流场更新和/或计算终止；C)重复步骤A)和B)，交替进行流场更新和瞬态传热计算，直到到达瞬态传热终止时刻。针对管内空气强迫对流换热过程，利用本发明的建模方法与Fluent商业软件紧耦合计算结果相比，绝对误差在1K以下，计算效率提高一个数量级，满足工程计算需求。

Description

流体－固体耦合传热的松耦合建模方法

技术领域

本发明涉及一种流体－固体耦合传热的松耦合建模方法。

背景技术

目前，单相流固耦合传热问题在商业软件中多使用瞬态CFD(计算流体力学)方法求解流场并与固体传热实时耦合的计算方法。在通用商业CFD软件FLUENT中，有4中主要的瞬态CFD方法：SIMPLE(压力耦合方程半隐方法)，SIMPLEC(协调一致的压力耦合方程半隐方法)，PISO(压力的隐式算子分割算法)，FSM(有限状态机)。在每一个时间步的迭代中，流体区域的流场与温度场耦合，流体区域的温度场与固体区域的温度场耦合。该算法具有较高的精度和广泛的使用范围，但是计算效率较低。

基本松耦合算法完全规避了瞬态CFD计算，并且完全解除了流体与固体之间的耦合。在该算法中，流体区域与固体区域分别进行求解。其中，流体使用稳态CFD求解，固体使用瞬态热传导求解。通常情况下，在处理流固耦合边界面时，对流体一侧的计算使用温度边界条件(第一类边界条件)，对固体一侧的计算使用热流边界条件(第二类边界条件)。该算法的计算效率比商业软件通用算法高一个量级以上，但精度较低且不符合物理实际。流程如下：

1)更新流场：将流体单独加载入求解器，将流固耦合壁面设为流体的固定温度边界，用稳态CFD算法求解流场；

2)计算瞬态传热：将固体单独加载入求解器，将流固耦合壁面设为固体的固定热流边界，计算瞬态传热直到下次流场更新或计算终止；

3)重复1)、2)两步，交替进行流场更新和瞬态传热计算，直到计算中止。

流固交界面实时耦合的松耦合算法部分规避了瞬态CFD计算，并且完全保持流体与固体之间的耦合。在处理瞬态传热过程时，不是认为流体区域内的所有参数处于准稳态，而是只认为速度、压强等流场参数处于准稳态，进而允许流体的温度在瞬态计算中与固体区域温度耦合变化。该算法算法的计算效率和精度都与商业软件通用算法处于同一量级。流程如下：

1)在一个时间步长t1内同时计算流场和温度场，也就是说在t1的时间步长内对能量方程和动量方程同时进行迭代求解；

2)固定流场，在下一个时间步长t2内不再进行计算流场变化，即不再对流体动量方程进行计算，只对能量方程进行迭代求解一个步长t2，得出温度场；

3)固定温度场不变，在t1时间内继续同时求解流场；

4)重复步骤1)到步骤3)，直到求解结束。

发明内容

根据本发明的一个方面，提供了一种流体－固体耦合传热的松耦合建模方法，其特征在于忽略流场的瞬态变化过程，假设全局瞬态传热过程在“准稳态”流场下进行。具体计算流程包括：

A)更新流场：单独以流体作为求解对象，将流固耦合壁面设为流体的固定温度边界，用稳态CFD算法求解流场；

B)计算瞬态传热：同时以流体和固体为求解对象，将流固耦合壁面设为传热耦合边界，关闭流体的动量方程和湍流方程，计算瞬态传热直到下次流场更新和/或计算终止；

C)重复步骤A)和B)，交替进行流场更新和瞬态传热计算，直到到达瞬态传热终止时刻。

附图说明

图1是根据本发明的一个实施例的流体－固体耦合传热的松耦合建模方法的总体流程图。

图2根据本发明的一个实施例的流体－固体耦合传热的松耦合建模方法的流体与固体间的实时耦合传热计算中使用的分区计算、边界耦合方法的流程图。

图3示意显示了矩形管空气强制对流模型，用于检验根据本发明的流体－固体耦合传热的松耦合建模方法的可靠性。

图4显示了Fluent软件紧耦合的一个实例的计算结果。

图5显示了根据本发明的一个实施例的流体－固体耦合传热的松耦合建模方法的一个应用实例的计算结果。

具体实施方式

根据本发明的松耦合算法一方面完全规避瞬态CFD计算，保证计算效率；另一方面保持瞬态传热计算时流体与固体的耦合，贴近物理实际。根据本发明的一个实施例的算法流程包括：

1)更新流场：单独以流体作为求解对象，将流固耦合壁面设为流体的固定温度边界，用稳态CFD算法求解流场；(步骤101-104)

2)计算瞬态传热：同时以流体和固体为求解对象，将流固耦合壁面设为传热耦合边界，关闭流体的动量方程、湍流方程，计算瞬态传热直到下次流场更新或计算终止；(步骤105-108)

3)重复1)、2)两步，交替进行流场更新和瞬态传热计算，直到到达瞬态传热终止时刻。其总体算法流程图如图1所示。

2.1稳态CFD流场计算

图1中左边一列为流场更新计算(步骤101-104)。在更新流场时(步骤101)，单独以流体作为求解对象，将流体单独加载入求解器。

在步骤102，将流固耦合壁面设为流体区域的温度边界。

在步骤103，开启流体区域的动量方程与湍流方程，进行稳态CFD流场计算。

稳态CFD求解可以使用通用的稳态CFD方法(步骤104)，具体包括：

2.1.1确立控制方程

控制方程的通式为：

▽·(ρVφ)＝▽·(Γ▽φ)+S

其中，ρ表示密度、t表示时间、V表示速度矢量、S表示源项、Γ表示广义扩散系数、φ代表方程变量。具体的，φ在各方程中分别代表x方向速度分量u、y方向速度分量v、z方向速度分量w、湍流动能k、湍流耗散率ε以及温度T。

对于各变量，广义扩散系数Γ为：

u,v,w:Γ＝η+η_t

其中η、η_t分别表示粘性系数和湍流粘性系数，Pr与σ_T分别表示普朗特数与湍流普朗特数，σ_k与σ_ε分别表示k与ε的湍流Schmidt数。

对于各变量，源项S为：

k:S＝ρG_k-ρε

T：S按实际情况而定

其中，c₁、c₂为经验常数，G_k为：

2.1.2用SIMPLE算法求解控制方程

对于经过有限体积法离散的流体控制方程，使用SIMPLE算法求解，实现速度与压力的耦合。具体步骤如下：

1)假定一个速度分布，记为V⁰，以此计算动量离散方程中的系数和常数项；

2)假定一个压力场p*；

3)依次求解各方向上的动量方程，得到V*；

4)求解压力修正方程，得p’；

5)据p’改进速度值；

6)利用改进后的速度场求解那些通过源项无形等与速度场耦合的变量，如果变量并不影响流场，则应在速度场收敛后再求解；

7)利用改进后的速度长重新计算动量离散方程的系数，并用改进后的压力场作为下一层迭代的初值，重复上述步骤，直到获得收敛解。

对于二维结构化网格，速度修正计算式为：

u_e＝u_e*+d_e(p'_P-p'_E)，v_n＝v_n*+d_n(p'_P-p'_N)

压力修正方程为：

a_Pp'_P＝a_Ep'_E+a_Wp'_W+a_Np'_N+a_Sp'_S+b

其中：a_E＝ρ_ed_eΔy，a_W＝ρ_wd_wΔy，a_N＝ρ_nd_nΔx，a_S＝ρ_sd_sΔx

a_P＝a_E+a_W+a_N+a_S

ρ表示密度，下标P表示当前节点，下标E、W、N、S分别表示四周的邻点，下标e、w、n、s分别表示四周的表面，d表示相应方向上控制体表面面积与动量方程系数的商。

在SIMPLE算法中，引入了三方面的简化假定：1.速度场的假定和压力场的假定是各自独立进行的，两者间无任何关系；2.在导出速度修正值计算式时没有计及邻点速度修正值的影响；3.采用线性化了的动量离散方程，即在每一个层次的计算中，动量离散方程中的各个系数及源项b假定均为定值。

在实施SIMPLE算法的过程中，速度与压力的修正值都应作亚松弛处理，但实施的方式有所不同。对压力由于在速度修正值公式中略去了邻点的影响，所解得的p’修正速度是合适的，但对压力修正值本身，则是被夸大了，因而需要亚松弛，直接对其进行亚松弛处理，亦即作为这一迭代层次的解为：

p＝p*+α_pp'

其中α_p为压力松弛因子。

对速度，为限制相邻两层次之间的变化，以利于非线性问题迭代收敛，也要求亚松弛。将亚松弛过程组织到代数方程的求解过程中。将代数方程的对角元系数a_p修正为将源项b修正为其中表示变量上一层次的解，α为压力松弛因子。

在同一层次上代数方程的迭代(内迭代)中，p’方程求解是关键，常常占据内迭代大部分时间，因而终止内迭代常常以该方程为依据，主要有三种方法：简单规定实施交替方向线迭代与块修正运算轮数；规定p’方程余量的范数小于某一数值；规定终止迭代时的范数与初始范数之比小于允许值。

在非线性问题的迭代(外迭代)中，有四种主要的判定终止迭代方法：特征量在连续若干个层次迭代中的相对偏差小于允许值；要求在内节点上连续性方程余量的代数和及节点余量的最大绝对值小于一定的数值；要求连续性方程余量范数的相对值小于允许值；要求在整个求解区域内动量方程余量之和或其范数与参考动量之比小于一定值。

2.2瞬态导热计算

流程图中第二列为瞬态传热计算(步骤105-108)。为保证物理合理性，瞬态传热计算同时以流体和固体为求解对象，将流体与固体同时加载入求解器，将流固耦合壁面设为传热耦合边界，关闭流体动量、湍流方程，使用实时耦合的瞬态数值方法求解传热(步骤108)。

流体与固体间的实时耦合传热计算可以使用分区计算、边界耦合的方法，如图2所示，其操作步骤包括：

1)分别对各个区域中的物理问题建立控制方程。(步骤201)

2)列出对每个区域的边界条件，其中耦合边界条件上的条件取为：耦合边界上温度连续，耦合边界上热流密度连续。(框201)

3)假定耦合边界上的温度分布，作为区域1的边界条件。(步骤202)

4)对其中区域1进行求解，得出耦合边界上的局部热流密度和温度梯度，依据耦合边界上热流密度的连续关系获得区域2的边界条件。(步骤203、204)

5)求解区域2，得出耦合边界上新的局部热流密度和温度梯度，依据耦合边界上热流密度的连续关系获得区域1的边界条件。(步骤205、206)

6)重复第4、5两步计算直到收敛。

当采用无量纲控制方程时，应注意在耦合边界上无量纲温度定义间的一致性，以利于信息的传递。这种计算方法迭代过程收敛的快慢主要取决于耦合边界上的信息传递。

在实时耦合的瞬态传热计算中，处理流体区域及固体区域的具体技术途径包括：

2.2.1流体区域瞬态计算

考虑流场达到准稳态，此时关闭流场的动量、湍流方程，不再对流场进行更新。此时，流体区域控制方程通式为：

其中，φ只代表温度T。其它变量视为定值，其控制方程皆不求解。这样，也就不需要在流场计算中使用SIMPLE算法。

2.2.2固体区域瞬态计算

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

其中：

ρ＝密度

k＝导热率

T＝温度

S_h＝体积热源

固体导热控制方程左边第一项表示固体随时间的温度变化，右边两项分别表示传导引起的热流以及固体内部的体积热源。对于各向异性导热的，其热传导项为▽_·(k_ij▽T)，其中k_ij为导热率张量。

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

2.3流场更新时刻的确定

根据本算法的瞬态传热计算方法，只有当外部边界条件或流体温度场的变化影响了流体区域的稳态流场，并使流固对流换热系数(努赛尔数)发生较大变化时，才需要对流场进行更新。实施中，可以用努赛尔数经验公式来判断努赛尔数是否发生较大变化，需要更新流场。不同的工程问题的努赛尔数经验公式不同，例如管槽内层流强制对流有：

其中，Nu、Re、Pr分别表示努赛尔数、雷诺数和普朗特数，η表示粘性系数，l与d分别表示管长、管径。下角标f、w分别表示流体和管壁。

特别的，当流体为常物性或流体物性在工况温度区间内变化很小时，只需在外边界条件发生较大改变时更新流场。如果外边界条件不发生改变，在初始时刻求得稳态流场后即不需再更新流场。

2.3算法可靠性验证

为了检验本算法的可靠性，我们选择了如图3所示的矩形管空气强制对流模型。矩形管(301)长0.328m，内壁(304)横截面0.0104m*0.0064m，外壁(303)横截面0.013m*0.0083m；流道(302)被内壁(304)所围绕限定。如图3所示。管(301)内通入2m/s，温度400K的空气。固体初始温度为300K，外表面绝热，加热过程持续600s。由于流体与固体都取作常物性，并且入口条件不变，所以松耦合计算只需在初始时刻更新一次流场，然后关闭动量、湍流方程完成瞬态传热计算即可。另外，由于瞬态计算不开启与流场有关的方程，所以使用较长的瞬态时间步长。

Fluent软件紧耦合计算结果与本松耦合算法计算结果如表1、图4和图5所示。计算结果取固体管上距入口10％长度、50％长度、90％长度的三个点上的温度，作为管入口附近、中段、出口附近的温度，单位为K：

表1

从图4、图5的计算结果可以看出，管体的温度沿气流方向递减。随着管体温度逐渐接近400K，其升温速度也在逐渐下降。比较Fluent紧耦合计算结果与本松耦合算法的计算结果，绝对误差在1K以下，满足工程计算需求。

在计算效率方面，本松耦合算法用时2小时，Fluent实时耦合计算用时13.5小时。本算法可以保证计算精度的情况下，将计算效率提高一个量级。

Claims (8)

1.流体－固体耦合传热的松耦合建模方法，其特征在于包括：

A)更新流场：单独以流体作为求解对象，将流固耦合壁面设为流体的固定温度边界，用稳态CFD算法求解流场；

B)计算瞬态传热：同时以流体和固体为求解对象，将流固耦合壁面设为传热耦合边界，关闭流体的动量方程和湍流方程，计算瞬态传热直到下次流场更新和/或计算终止；

C)重复步骤A)和B)，交替进行流场更新和瞬态传热计算，直到到达瞬态传热终止时刻，

其中，忽略了流场的瞬态变化过程，并假设全局瞬态传热过程在“准稳态”流场下进行，

更新流场时单独对流体进行稳态CFD计算，

稳态CFD算法求解流场的处理包括：

确立控制方程，控制方程的通式为：

其中，ρ表示密度、t表示时间、V表示速度矢量、S表示源项、Γ表示广义扩散系数、φ代表方程变量，其中φ在各方程中分别代表x方向速度分量u、y方向速度分量v、z方向速度分量w、湍流动能k、湍流耗散率ε以及温度T，

对于各变量，把广义扩散系数Γ确定为：

u,v,w:Γ＝η+η_t

其中η、η_t分别表示粘性系数和湍流粘性系数，Pr与σ_T分别表示普朗特数与湍流普朗特数，σ_k与σ_ε分别表示k与ε的湍流Schmidt数，

对于各变量，把源项S表示为：

k:S＝ρG_k-ρε

其中，T：S按实际情况而定，c₁、c₂为经验常数，p为压力，G_k为：

2.根据权利要求1所述的建模方法，其特征在于对于经过有限体积法离散的流体控制方程，使用SIMPLE算法求解，实现速度与压力的耦合，包括：

假定一个速度分布，记为V⁰，以此计算动量离散方程中的系数和常数项；

假定一个压力场p*；

依次求解各方向上的动量方程，得到V*；

求解压力修正方程，得p’；

据p’改进速度值；

利用改进后的速度场求解那些通过包括源项和物性的与速度场耦合的变量，如果变量并不影响流场，则应在速度场收敛后再求解；

利用改进后的速度场重新计算动量离散方程的系数，并用改进后的压力场作为下一层迭代的初值，重复上述步骤，直到获得收敛解，

对于二维结构化网格，速度修正计算式为：

u_e＝u_e*+d_e(p'_P-p'_E)，v_n＝v_n*+d_n(p'_P-p'_N)

其中压力修正方程为：

a_Pp'_P＝a_Ep'_E+a_Wp'_W+a_Np'_N+a_Sp'_S+b

其中：a_E＝ρ_ed_eΔy，a_W＝ρ_wd_wΔy，a_N＝ρ_nd_nΔx，a_S＝ρ_sd_sΔx

a_P＝a_E+a_W+a_N+a_S

式中ρ表示密度，下标P表示当前节点，下标E、W、N、S分别表示四周的邻点，下标e、w、n、s分别表示四周的表面，d表示相应方向上控制体表面面积与动量方程系数的商。

3.根据权利要求2所述的建模方法，其特征在于在实施SIMPLE算法的过程中，对速度与压力的修正值作亚松弛处理，

其中

对压力修正值本身直接进行亚松弛处理，即作为这一迭代层次的解为：

p＝p*+α_pp'

其中α_p为对压力修正值进行亚松弛处理的压力松弛因子；

对速度，为限制相邻两层次之间的变化，以利于非线性问题迭代收敛，也进行亚松弛处理，包括将亚松弛过程组织到代数方程的求解过程中，将代数方程的对角元系数a_p修正为将源项b修正为其中表示变量上一层次的解，α为对速度修正值进行亚松弛处理的压力松弛因子。

4.根据权利要求1所述的建模方法，其特征在于计算瞬态传热的处理包括使用分区计算、边界耦合的方法计算流体与固体间的瞬态传热，具体包括：

分别对各个区域中的物理问题建立控制方程，

列出对每个区域的边界条件，其中耦合边界条件上的条件取为：耦合边界上温度连续，耦合边界上热流密度连续，

把耦合边界上的一个假定的温度分布作为区域1的边界条件，

对区域1进行求解，得出耦合边界上的局部热流密度和温度梯度，依据耦合边界上热流密度的连续关系获得区域2的边界条件，

对区域2进行求解，得出耦合边界上新的局部热流密度和温度梯度，依据耦合边界上热流密度的连续关系获得区域1的边界条件，

重复之前最后两个步骤，直到计算结果收敛。

5.根据权利要求4所述的建模方法，其特征在于对流体区域建立只考虑温度的瞬态控制方程：

控制方程通式为：

其中，φ只代表温度T，其它变量视为定值，其控制方程皆不求解。

6.根据权利要求4所述的建模方法，其特征在对固体区域进行只考虑导热的瞬态传热计算，该包括：

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

其中：

ρ＝密度

k＝导热率

T＝温度

S_h＝体积热源

其中，c_p是定压比热，T_ref是基准温度，h为显焓，上述控制方程左边第一项表示固体随时间的温度变化，右边两项分别表示传导引起的热流以及固体内部的体积热源，

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

7.根据权利要求4所述的建模方法，其特征在于流场更新的判别依据包括：

仅当外部边界条件或流体温度场的变化影响了流体区域的稳态流场，并使流固对流换热系数发生较大变化时，才对流场进行更新。

8.根据权利要求7所述的建模方法，其特征在于流场更新的判别包括基于努赛尔数阈值的判别方式：

用努赛尔数经验公式来判断流固对流换热系数的努赛尔数是否发生较大变化。