耦合传热描述的是一种热量在流体域和固体域交互传递的物理现象.这一现象常见于许多工程实际应用中，如换热器中流体流动带来的热量传递^[1]，填充床蓄热式热交换器内部的热量传递^[2]，高炉炉缸内部高温铁水和外部冷却循环系统之间的热量传递^[3]等.

求解流固耦合传热问题一般包括两种方法，一种是数值模拟法，另一种是理论分析法.随着计算机技术在计算速度和精度方面的快速发展，数值模拟法以其成本低和设计周期短的特点受到广泛欢迎.目前，国外学者对与耦合传热相关的数值计算方法已进行了较为丰富的研究^[4,5,6,7,8]，而我国在这方面的研究还比较欠缺.

在数值计算中，方法的稳定性、效率和精确性是人们最关心的问题.本文基于Goudonov- Ryabenkii (G-R)稳定性理论，对一维隐式流固耦合稳态传热问题进行了稳定性分析，采用常用的有限体积法(流体域)(FVM)和有限单元法(固体域)(FEM)，应用Dirichlet-Robin组合形式的边界条件，重点研究了交界面移动时的数值方法的稳定性，最终获得一条由耦合系数和移动速度组成的最优曲线，并且证明了当耦合系数和移动速度在这条曲线上取值时，离散的求解域能够达到最快的收敛速度及恒定的稳定性特征.

为简化问题，流体域和固体域的所有物理参数都假设为常数.稳态问题中，交界面移动的情况下，固体的导热微分方程为一维Laplace方程：

流体域的求解需要在时间方向上进行推进，即求解非守恒形式下的N-S方程：

除上述单个子域的控制方程外，在交界面处，需要满足温度和热流的连续性方程组：

Roux等^[9]已经证明Dirichlet-Robin边界条件比Dirichlet-Neumann更容易得到收敛结果并能得到更好的稳定性.因此，本文中直接采用这一组合边界条件，将Dirichlet条件(温度)作为边界施加到流体域，而Robin条件(可视作热流的修正结果)作为边界施加到固体域，此时温度和热流的连续性方程组将变为

在热流计算中暗含了温度梯度与坐标轴正向一致的假设，由于热流方向与温度梯度方向相反，因此可假设热流方向均与坐标轴负向一致.

对偏微分方程(1)和(2)进行数值计算之前，v₀需要将流固系统离散化，原理如图 1所示.

图 1(Figure 1)

Fig. 1

图 1 流固系统离散化原理图 Fig. 1 Schematic of fluid-structure system discretization

如图 1所示，Δx_f和Δx_s分别为流体域和固体域的空间步长；Λ_s为固体域总长；v₀为交界面移动速度；时间步长由流体域的时间步长决定，用Δt表示.

在外部边界的传热计算中，同样采用Robin边界条件：

将FVM应用到流体中，靠近边界处的一点距离边界的长度为空间步长Δx_f的一半^[10]，即点1的坐标为Δx_f /2，其余点之间的距离为Δx_f.

遵照通用的表示方式，用T_jⁿ代表j位置处n时刻的温度值，为求解T_jⁿ⁺¹，对式(1)和式(2)进行差分，交界面处的温度用传统的顺序交错方式^[11]进行迭代求解，迭代过程包括以下4个步骤：

不断执行以上4个步骤直至得到收敛的温度和热流.

由于本文讨论的问题并非周期性变化而且边界条件不可忽略，因此，Von-Neumann稳定性分析方法并不适用.借鉴Giles的研究方法^[12]，本文采用G-R理论，通过寻求系统的正则模态解来判断离散的流固系统(第1节末尾的循环迭代系统)的稳定性.

假设系统有形如公式(7)的正则模态解

G-R理论的完整表述是若离散的系统在 j→±∞ 时，在|R_f|<1，|R_s|>1和|z|>1的条件下有形如式(7)的解，则离散的系统不稳定.将式(7)代入迭代求解步骤中去，经过化简，得到以下两个方程：

依据G-R理论，在|z|≥1复平面内max(|g(z)|)<1是离散系统稳定的一个充要条件，因此，研究max(|g(z)|)在|z|≥1内的变化情况是有必要的.令Z=1/z，则将开放域|z|≥1转化成为有限闭区域|Z|≤1，则式(9)变为

由复变函数基本理论可知，g(Z)在|Z|<1范围内全纯，在边界|Z|=1上连续，则根据最大模原理可知g(Z)只可能在边界上取得最大值，即max(|g(z)|)=max[|g(|z|=1)|].此时，z的模r为1，辐角θ由0到2π变化，即z=e^iθ，代入式(8)和式(9)后可以得到

由式(11)可知，只有当|f(θ)|取最小或最大值时，|g(|z|=1)|才有可能取得最大值.利用一阶导数为0，二阶导数不为0的方法容易发现，|f(θ)|在θ=0或2π处取最小值，在θ=π处取最大值.上述计算过程可在Matlab上实现.因此，max(|g(z)|)只可能发生在z=1(θ=0或2π)或z=-1(θ=π)上，将z=±1代入式(9)，可得

由于|g(z=1)|和|g(z=-1)|的大小关系尚未确定，因此，若令|g(z=1)|=|g(z=-1)|，可得

根据式(14)可得到两条曲线相交时α_flu的值，并定义为α_flu^opt，具体表达式为

结合式(13)和式(15)发现，|g(z=-1)|随α_flu变化而连续变化且在2α_flu^opt处发生单调性质的突变.当α_flu <2α_flu^opt时，|g(z=-1)|单调递减；当α_flu >2α_flu^opt时，|g(z=-1)|单调递增.

由式(12)可知，|g(z=1)|关于α_flu是单调递增函数.令|g(z=1)|=|g(z=-1)|，可得二者交点坐标为

根据函数的单调性质可知，当α_flu <α_flu^opt时，|g(z=-1)|>|g(z=1)|，max (|g(z)|)=|g(z=-1)|；当α_flu >α^opt_flu时，|g(z=-1)|<|g(z=1)|，max(|g(z)|)=|g(z=1)|；当α_flu=α_flu^opt时，|g(z =-1)|=|g(z=1)|，max(|g(z)|)=|g(z=-1)|= |g(z=1)|，且max(|g(z)|)取得最小值，即

用opt标注α_flu及max(|g(z)|)主要有以下原因：

① 由式(9)可知，max(|g(z)|)的物理意义是最大时间放大因子，而当α_flu取α_flu^opt时，max(|g(z)|)取得最小值，意味着由初始温度达到最终温度所需时间最少，说明在这一点系统能获得最快的收敛速度；

② 由式(17)可知，在α_flu取α_flu^opt时，时间放大因子总是小于1，说明在这一点处离散系统是绝对稳定的.

由上述内容可知，max(|g(z)|)准确的表达式为

由式(15)可知，在材料参数已知且求解域完成离散化后，α_flu^opt成为仅仅与交界面的移动速度有关的物理量且为单值对应，这说明每一个速度值都有唯一的最优耦合系数与其对应，以至于有唯一最优的最大时间放大因子与其对应.不难想象，在由v₀，α_flu和max(|g(z)|)三者构成的三维空间中，存在着一条由v₀^opt，α_flu^opt和{max(|g(z)|)}^opt构成的最优曲线，v₀和α_flu在这条最优曲线上取值，必定能得到同等水平(v₀或α_flu固定时，另一个参数变化的情况)中最快的收敛速度和绝对稳定的系统.

第2节中的结论均由理论公式推导而得，因此对于任意数据均成立.任选物理参数(如表 1所示)对以下两个结论进行实验，采用不同的时间步并假设β=1.

表 1

Table 1

表 1(Table 1)

表 1 流体域和固体域物理参数 Table 1 Parameters of fluid and solid domains 流体部分 κ f Δ x f ρ f c f κ f/Δ x f W·K -1·m -1 m kg·m -3 J·kg -1·K -1 W·K -1·m -2 0.02 10 -6 1 1 000 20 000 固体部分 κ s Λ s κ s/ Λ s W·K -1·m -1 m W·K -1·m -2 20 10 -2 2 000

表 1 流体域和固体域物理参数 Table 1 Parameters of fluid and solid domains

①当v₀和α_flu在这条最优曲线上取值时，必定能得到同等水平中最快的收敛速度；

②当v₀和α_flu在这条最优曲线上取值时，必定能得到绝对稳定的系统.

对于第一个结论，根据式(15)，可以把α_flu^opt视作速度v₀的函数(或者反过来，本文中采用公式(15)的表述)，式(18)对任意速度都是适用的.因此任意选取速度值，得出如图 2所示的结果.

图 2(Figure 2)

Fig. 2

图 2 不同时间步和速度下max(|g(z)|)随耦合系数的变化曲线 Fig. 2 Variation of max(|g(z)|) with coupling coefficient under different time steps and velocity (a)—v0=0； (b)—v0=1 m/s; (c)—v0=10 m/s； (d)—v0=100 m/s.

由图 2可知，无论速度和时间步怎样选取，当α_flu取α_flu^opt时，max(|g(z)|)取得最小值，此时，系统能达到最快的收敛速度，由于数据选取的任意性，此结论始终成立.

对于第二个结论.综合式(15)和(17)，可得max(|g(z)|)在这条最优曲线上关于v₀的表达式，即

由式(19)可知，无论v₀取何值，{max(|g(z)|)}^opt恒小于1，依据G-R理论，此时离散流固耦合系统绝对稳定，此结论始终成立.

在Dirichlet-Robin组合边界条件下，考虑交界面移动的情况，对一维隐式流固耦合稳态传热系统进行数值求解，结果表明，存在一条由耦合系数和交界面移动速度组合而成的最优曲线，该曲线具有以下两个重要特征： ①当速度和耦合系数在这条曲线上取值时，离散的流固耦合系统能得到最快的收敛速度； ②当速度和耦合系数在这条曲线上取值时，离散的流固耦合系统能绝对稳定.

本文的结论可以为设计人员解决工程问题时提供合理的参考，不仅可以节约仿真步数，还能避免因参数选取不合理而导致仿真失稳的问题.