非稳态共轭传热中宽频热波动的多尺度时域方法

《International Journal of Heat and Mass Transfer》：A multiscale time-domain method for broadband thermal fluctuations in unsteady conjugate heat transfer

　　

•基于网格的截止准则根据扩散尺度对宽频热波动进行划分。•低频波动通过多维方法求解，高频波动则采用递归卷积方法处理。•复频域传递函数拟合适用最大模原理实现误差控制。•该方法在时域中处理周期性、宽频和非周期性热扰动。•热收敛时间降低至同步参考算例的2.0–8.04%。1. 引言共轭传热是一种耦合方法，能够求解涉及固体结构与流体流动相互作用的复杂热问题[1,2]。CHT方法在对流换热系数与壁面条件强耦合的场景中更为适用[3]，例如非均匀壁温分布、固体内部热传导以及近壁区域的湍流输运过程。随着计算机硬件的不断进步，对能够解析湍流的计算流体动力学（CFD）模拟（如大涡模拟（LES））以及非定常共轭传热计算的需求日益增长[4,5,6,7,8]。流体对流时间尺度与固体导热时间尺度之间的差异可能非常大，可跨越四个数量级以上[9]。因此，为准确解析流体扰动并确保固体对流体诱导扰动产生足够的热响应，瞬态CHT计算需要采用极小的时间步长，从而需要极其大量的时间推进迭代。非定常CHT本质上是一个多尺度问题，固体的热响应受到物理和数值因素的综合影响。这些因素包括固体的热物性、热载荷的特征时间尺度和频谱内容、由此产生的热穿透深度、固体域中的网格[10,11]和时间步长分辨率[12,13,14]，以及流体求解器和固体求解器之间使用的界面耦合过程[15,16,17,18]。一般来说，短时间尺度波动与高度局域化的热响应相关，而长时间尺度波动则影响固体中相当大的区域[19]。由于涉及的时间尺度范围很广，不同特征时间尺度的热波动可能需要不同的数值处理。本研究聚焦于该多尺度问题的框架层面，旨在开发一种高效的非定常CHT计算方法。对于短时间尺度的波动，其相关高频导致非常小的热穿透深度。根据Hickling和He[9]的研究结果，需要足够精细的固体网格才能准确捕捉这些波动，这增加了额外的计算负担，并促使开发了无需网格细化即可处理高频波动的计算方法。Schultz[20]基于一维半无限假设推导了固体表面温度的传递函数。该传递函数建立了热流波动与固体表面温度波动之间的映射关系，能够捕捉单位热流波动对应的温度波动的振幅衰减和相位滞后。He等人[9,21]首次将该方法应用于非定常CHT模拟，证明了其在RANS（雷诺平均Navier-Stokes）和LES方法所获得的流场中具有优异的适用性。该CHT方法在频域而非传统时域中运行，非常适合计算影响区域较小的高频波动，尤其适用于湍流诱导的温度波动。对于长时间尺度的波动，其相关低频导致较高的热穿透深度，使得温度波动的影响区域较大。这类波动主要由相干流体结构诱导，常见于旋转部件和燃烧不稳定性现象中[22]。Hickling和He[23,24]提出了一种多尺度频域方法来解析低频波动。该方法通过求解给定频率下的频域热传导方程，能够准确计算指定频率处的温度波动，使其适合与高保真流体模拟结合使用。Koren等人[12,13]将固体温度分解为时间平均分量和波动分量。前者在稳态下求解，旨在人为加速固体时间平均温度向准稳态的收敛；后者在瞬态下以小自适应时间步长求解，以准确捕捉固体的非定常温度波动特征。此后，该方法已被扩展至涉及多介质的共轭传热问题[25]。对于发生在极长时间尺度上的波动，例如发动机启动期间的温升，固体与流体之间显著的时间尺度差异使得流体可被视为处于准稳态[26,27,28,29]。基于这一假设，已经开发了大量相关算法[30,31]来求解固体内时间平均温度的演化。在此基础上，He和Fadl[32]提出了一种以三时间形式表述的流体域多尺度框架。该方法在加速固体温度演化的同时保持真实的瞬态演化行为，并能同时捕捉由固体温度变化诱导的瞬态流体特征。数据驱动和物理信息方法最近的进展也为加速热模拟和CHT模拟提供了替代方案。基于PINN（物理信息神经网络）的方法已被应用于共轭传热问题，尽管其在多区域和多物理场系统中的应用可能涉及训练收敛和稳定性方面的挑战[33,34]。算子学习方法（如DeepONet）在离线训练阶段之后可以提供参数化热传导解的快速预测[35]。降阶模型也通过利用高保真解数据构建的低维表示，在非定常CHT问题上展示了显著的计算节省[36,37]。这些方法对于重复或参数化评估前景广阔。本研究聚焦于两个主要目标：第一，加速时间平均固体温度场向稳态的演化；第二，准确解析由流体扰动诱导的固体温度波动，包括宽频湍流、相干流结构和旋转部件。这两个方面在工程应用中具有重要的实际意义，特别是在LES与CHT的耦合方面。已开发了多种方法来应对这些挑战。Duchaine等人[4,5]在固体域中采用去同步化方法来加速温度演化，实现固体和流体温度的同步收敛，该方法也被商业软件Fluent用于CHT加速[38]。然而，该方法会导致固体中波动频率的低估，从而导致固体温度波动的错误放大[13]。另一些研究者采用了一种策略，即在初始阶段使用较大的固体时间步长获得近似温度场，随后以一致的时间步长进行时间推进。有研究者指出该方法受到初始瞬态问题的限制[23]。Koren等人[12,13]提出的算法采用稳态方法求解时间平均温度，从而在一定程度上缓解了初始瞬态问题。通过自适应调节流体-固体耦合周期，可以实现低频波动的有效计算。然而，当应用于特征时间尺度较小的高频波动时，固体时间步长受到极大的限制，时域方法不再提供显著的计算优势。Hickling和He[23,24]提出的频域算法直接在主导频率下求解热传导方程，能够直接获得准稳态结果，从而解决初始瞬态问题。此外，对于热穿透深度较小的高频热流波动，该方法利用温度传递函数直接确定壁面温度波动，无需采用极小时间步长进行时域推进。然而，对于一般流场配置，主导频率事先并不已知，通常需要通过先前的数值模拟来确定。在没有明显主导频率的情况下，或者主导频率在不同区域有所变化时，需要在多个频率下求解热传导方程，从而增加了计算负担。此外，应用于主导频率范围之外低频区域的温度传递函数基于一维假设推导。将其适用性扩展到热穿透深度显著且横向热传导不可忽略的低频区域，还需要进一步的验证。为克服这些局限，本研究开发了一种面向宽频热波动的共轭传热多尺度时域框架。温度场被分解为时间平均分量、低频波动分量和高频波动分量。时间平均温度采用稳态公式求解，以加速向热平衡的演化。然后，根据热穿透深度和固体网格分辨率确定的截止频率对完整的波动频谱进行划分。所有低于截止频率的波动（包括非主导低频分量）通过多维瞬态固体求解器求解，而高频分量则采用局部递归卷积壁面模型处理。因此，截止频率作为将波动频谱在多维数值解和局部一维传递函数模型之间进行分配的基于网格的准则，而不仅仅是一个信号处理参数。本研究进一步将基于递归卷积的传递函数方法[23]扩展到复频域，其中利用最大模原理结合边界误差估计，在指定区域内控制近似误差。该复域公式旨在为有限持续时间和非周期性瞬态响应提供更优的误差控制，同时保持递归卷积的计算效率。本文的结构组织如下。第2节详细介绍了所提出的方法，包括温度在不同时间尺度上的分解、区分不同尺度波动的准则及相应的数值方法，以及整体计算流程。第3节展示了该方法的数值评估，涵盖单介质固体、多介质固体和三维涡轮叶片中的各种典型工况。第4节总结全文结论。2. 方法2.1. 时间尺度分解方法根据第1节的描述，本研究建立的多尺度非定常CHT模拟方法需要实现以下目标：（1）鉴于固体传热时间尺度大的特点，应加速固体温度演化，以缓解初始瞬态的影响。（2）该方法应能准确捕捉跨多个时间尺度的热波动。对于低频波动，应充分考虑由较大热穿透深度引起的横向热扩散。对于高频波动，应采用温度传递函数的基本概念来建立热流波动与温度波动之间的映射关系，同时规避频域方法固有的频率选择要求。（3）为捕捉波动壁温对流体的影响，应建立在线时域方法。高频与低频区域之间的边界应明确定义，以便在不同频率范围内应用不同的数值方法。如图1所示，遵循先前研究者[12,23,25]采用的方法，非定常温度T被分解为三个不同的时间尺度分量：（1）T = T? + T′ + T″其中T?为时间平均温度，T′和T″分别为低频波动分量和高频波动分量。忽略热物性变化，考虑瞬态固体热传导方程：（2）?T/?t = α?2T其中α = k/(ρc)为热扩散率，k、ρ和c分别表示固体的导热系数、密度和比热容。将式（1）的温度分解代入式（2），得到时间平均温度T?和波动温度T′、T″的控制方程：（3）?T?/?τ = α?2T?（4）?T′/?t = α?2T′（5）?T″/?t = α?2T″其中τ为伪时间，用于保持收敛率并提高求解的数值稳定性。流体域由非定常Navier-Stokes（N-S）方程控制。

通量前推温度回推（FFTB）方法[17]被用于施加固-液界面的边界条件，其中相关瞬时量被分解为时间平均分量和波动分量。所有时间平均量均可通过移动平均法计算：

(6) φ???? = (1/(n+1))φ??? + (n/(n+1))φ??

其中上标n表示时间步。平均运算以递归方式实现，无固定平均窗口。当有效样本数变得过大时，将其重置为较小的值，同时保留当前的平均值，以对新采样数据保持足够的权重。此外，界面热流密度也被分解为三个分量：时间平均部分q?和波动部分q′、q″：

(7) q = q? + q′ + q″

2.2 温度波动求解

2.2.1 低频波动求解

当波动传播进入固体内部时，其振幅逐渐衰减。本研究将热穿透深度δ?定义为从表面量起振幅衰减到其表面值e?1处的距离[19]，其表达式为：

(8) δ? = √(α/(πf))

可以看出，当频率f相对较低时，热穿透深度δ?变得较大，横向热传导不可忽略。在本研究中，低频温度波动T′通过直接求解式(4)所示的能量方程获得。一方面，时域计算天然涵盖了该频段内所有频率的波动信号，无需像频域方法[23]那样选择主导频率。另一方面，由于低频波动对应的时间尺度相对较大，固体域的时间步可以相应地设置得较大。因此，固-液界面处的耦合周期可以适当地增加，从而降低固体中热传导相关的计算成本。

波动热流密度可通过从瞬时热流密度q中减去时间平均热流密度q?获得，它同时包含低频和高频分量，两者之间的截止频率记为f??。在每个时间步，采用二阶Butterworth滤波器在线提取频率低于f??的波动热流密度q′。先前的研究[39]已证明该滤波算法具有很高的精度。随后，q″可通过从总波动热流密度中减去q′获得。滤波主要重新分配截止频率附近两个分支之间的频谱内容，对最终温度响应的影响较小。

2.2.2 高频波动求解

对于高频温度波动，直接求解式(5)所示的热传导方程并非最优选择。在许多实际应用中，波动频率非常高，例如在爆震发动机中可能存在量级达数百万赫兹的主导频率[40]，这需要极小的时间步来求解式(5)，并且需要频繁地在固-液界面处交换信息，从而产生显著的额外计算开销。

现有研究表明，非同步耦合方案可以显著加速固体温度的收敛率，但代价是高估了表面温度波动的振幅[4]。如图2所示，对于一维问题，放大因子为√γ，其中γ = Δt?/Δt_f表示固体时间步增加的因子。对于三维问题，放大因子与横向扩散有关，缺乏解析表达式。此外，固体内部的相位分布也会发生畸变。

图2. 非同步耦合方案引起的表面温度振幅高估。

鉴于高频热流密度波动引起的热穿透深度非常小，采用一维假设来模拟高频波动。将一维高频波动能量方程即式(5)变换到复频域，并假设初始波动温度为零，复频域中的热传导方程可以表示为[41]：

(9) d2T?/dx2 ? sT?/α = 0

其中s = σ + iω表示复频域变量，T?(s) = ∫?^∞ T(t)e???dt。基于一维半无限假设，边界条件设定如下：

(10) x = 0: T? = T?_w
x = ∞: dT?/dx = 0，且T? = 0

上述二阶常微分方程（ODE）的解析解为：

(11) T?(x, s) = T?_w e?√(s/α) x

当在x = 0处施加Neumann边界条件q?_w = ?k dT?/dx|_{x=0}时，可以获得T?_w与q?_w之间的关系：

(12) q?_w = C?_Tq · T?_w，C?_Tq = k√(s/α)

其中C?_Tq表示温度传递函数，建立了固体表面处热流密度波动q?_w与温度波动T?_w之间的映射关系。

固体表面温度的时间变化通过逆Laplace变换获得：

(13) T_w(t) = L?1{T?_w(s)} = √(α/k) · √(π) · ∫?? q(τ)/√(t?τ) dτ

可以看出，通过这种方式获取每个时间步的壁面温度T_w(t)需要知道所有先前时刻的热流密度q(τ)，这意味着必须存储整个热流密度的历史。因此，存储成本随着计算的推进而增加，从而对这种方法构成了限制。该方法的另一个问题是，当考虑存在热障涂层（TBC）等特殊情形时，壁面温度随时间变化的表达式可能难以明确给出。

2.2.3 递归卷积法

采用适合递归卷积的传递函数来拟合温度传递函数：

(14) 1/C?_Tq = 1/(k√(αs)) ≈ Σ(m=1到M) b?/(s+a?)

其中M为拟合阶数，a?和b?均为实常数，且a? > 0以保证传递函数的稳定性。时域传递函数K_qT(t)可以表示为：

(15) K_qT(t) = L?1{1/C?_Tq} ≈ L?1{Σ(m=1到M) b?/(s+a?)} = Σ(m=1到M) b? e????

T_w(t)与传递函数K_qT(t)之间的关系如下：

(16) T_w(t) = ∫?? K_qT(t?τ)q(τ)dτ = Σ(m=1到M) b? ∫?? e???????? q(τ)dτ

定义z?(t) = b? ∫?? e???????? q(τ)dτ，则上式可以表示为：

(17) T_w(t) = Σ(m=1到M) z?(t)

对z?(t)应用Leibniz法则求导得到?? = b?(q(t) ? ∫?? e???????? q(τ)dτ)：

(18) ?? = b?(q(t) ? a? ∫?? e???????? q(τ)dτ)

重新排列上式得到：

(19) ?? + a? z? = b? q(t)

该一阶ODE可以通过递归计算求解，而无需存储壁面热流密度的完整历史。上述ODE的解通过积分因子法获得。将方程两边乘以e???得到：

(20) e??? ?? + e??? a? z? = e??? b? q(t) → d(e??? z?)/dt = e??? b? q(t)

将上式关于时间从0积分到t得到：

(21) e??? z?(t) ? z?(0) = b? ∫?? e??? q(τ)dτ

由定义可知z?(0) = 0。因此，z?(t)的表达式为：

(22) z?(t) = b? ∫?? e???????? q(τ)dτ

当时间步推进到t + Δt时：

(23) z?(t + Δt) = b? e??????Δ?? (∫?? e??? q(τ)dτ + ∫_t^{t+Δt} e??? q(τ)dτ) = e???Δ? z?(t) + b? ∫?^{Δt} e????Δ???? q(s+t)ds，其中s = τ ? t

假设在每个时间步内热流密度保持为q(t)不变，则式(23)可以简化为：

(24) z?(t + Δt) = e???Δ? z?(t) + q(t) · b?/a? · (1 ? e???Δ?)

该表达式使得壁面温度T_w(t)能够递归更新，仅需存储2M个拟合系数a?和b?，而无需保留壁面热流密度的完整历史。由相同材料组成的固体壁面可以共享同一组拟合系数。

2.2.4 递归卷积误差分析

在复频域中，σ代表指数衰减率，而ω代表波动的角频率。对于瞬态共轭热传递（CHT）计算，通常只关注特定区域D。区域D的范围由问题的时间尺度ω_min、ω_max、σ_min和σ_max确定。因此，只需确保在区域D内具有足够的拟合精度。

拟合误差函数E(s)定义如下：

(25) E(s) = 1/(k√(αs)) ? Σ(m=1到M) b?/(s+a?)

式(14)中所示的传递函数定义在复频域中，只能在有限个点上进行拟合。幸运的是，根据复分析中的最大模原理[42]，如果一个复值函数f在复平面的一个连通开子集上是全纯函数且不是常数函数，那么|f|的最大值必然在边界上取得。本研究考虑的定义在σ > 0且ω > 0区域的复值函数E(s)是全纯函数，选取√s的主支，从而满足最大模原理的适用条件。因此，控制区域D边界（记为d）上的拟合误差就能控制整个区域的拟合精度，如图3所示。

图3. 感兴趣区域和拟合点。

根据最大模原理，区域D内的拟合误差与其边界d上的最大拟合误差之间的关系为：

(26) |E(s)|_{s∈D} ≤ max_{s∈d} |E(s)|

仅在边界上选择了有限个点进行拟合，这不足以获得整个边界上的最大拟合误差。因此，需要利用误差函数的导数来确定边界上的最大误差。对于边界上相邻拟合点p?和p?之间形成的长度为h_i的线段d_i，该边界线段上拟合误差上界的表达式为：

(27) |E(s)|_{s∈d_i} ≤ max(|E(p?)|, |E(p?)|) + (h_i/2) · max_{s∈d_i} |E′(s)|

对于式(25)中给出的E(s)，其导数表示为：

(28) |E′(s)|_{s∈d_i} = |?√α/(2k) · s?3?2 + Σ(m=1到M) b?/(s+a?)2| ≤ |√α/(2k) · d??3?2| + Σ(m=1到M) |b?|/d?2

其中d?表示该线段到原点的最小距离，d?表示该线段到极点?a?的最小距离。边界拟合误差和区域拟合误差的上界由所有边界线段拟合误差上界的最大值给出：

(29) |E(s)|_{s∈D} ≤ max_{s∈d} |E(s)| ≤ max_{d_i∈d} (sup_{s∈d_i} |E(s)|)

由式(15)可知，a?的量纲为1/t。因此，在进行拟合之前，根据问题的时间尺度确定a?。参数a?在对数均匀空间中分布于[a_min, a_max]范围内，其中a_min和a_max的取值如下：

(30) a_min = min(ω_min/(2π), 1/τ_max)
a_max = max(ω_max/(2π), 1/τ_min)

其中ω_min和ω_max分别为感兴趣范围内的最小和最大角频率，τ_max设置为预期的计算总时长，τ_min设置为流体时间步。随后，求解最小二乘问题以最小化绝对误差来确定b?。本研究与Hickling和He[23]的工作在计算高频温度波动方面最显著的区别在于：他们仅沿虚轴拟合传递函数，而本研究沿感兴趣时间尺度的边界进行拟合。这种方法的优势在于，借助最大模原理，可以保证在整个感兴趣域内的拟合精度，而不仅仅是在其边界上。

随后测试该方法在固体表面近似传递函数的能力。选择钛作为固体材料，其热物性参数列于表1中。复频域中的拟合范围定义为τ_min = 1×10?? s，τ_max = 10 s，ω_max = 1×10? rad/s，ω_min = 1 rad/s。表2列出了拟合区域内最大相对误差随区域每侧拟合点数N和拟合函数阶数M的变化情况。可以观察到，随着阶数M的增加，拟合精度不断提高。当M = 16且N = 1000时，整个域内的相对误差小于2%，表明精度很高。可以调整M使得相对拟合误差E_r(s)保持在预设阈值E_r,tol以下。E_r(s)定义为：

(31) E_r(s) = sup_{s∈D} |?(s) ? ?_M(s)|/|?(s)| ≤ E_r,tol

其中?(s) = 1/C?_Tq(s)，?_M(s)为拟合获得的传递函数。M的选择准则如下：

(32) M_fit = min{M : E_r(s) ≤ E_r,tol}

表1. 固体材料的热物性参数。

材料 | ρ (kg/m3) | c (J/(kg·K)) | k (W/(m·K)) | α (m2/s)
---|---|---|---|---
钛 | 4500 | 570 | 7.3 | 2.85×10??
Kapton | 1420 | 1090 | 0.155 | 1×10??

表2. 最大相对拟合误差随拟合点数和拟合阶数的变化。

拟合阶数 | N = 1000 | N = 5000 | N = 10000 | N = 50000
---|---|---|---|---
M = 4 | 1.7223 | 1.7237 | 1.7238 | 1.7240
M = 8 | 0.1295 | 0.1306 | 0.1308 | 0.1308
M = 16 | 0.0152 | 0.0047 | 0.0046 | 0.0046
M = 24 | 0.0147 | 0.0098 | 0.0073 | 0.0039

关于拟合阶段的参数选择，M和N主要决定拟合精度，当拟合精度不足时应适当增加。同时，τ_min和τ_max定义了卷积系统能够保存在其"记忆"中的时间范围。τ_min可设为流体时间步，而τ_max可设为预期的总计算时间。

2.2.5 网格分辨率要求

本节考察固体网格分辨率对所提方法结果的影响。考虑到该方法适用于三维问题，需要考虑壁面法向和横向的网格尺寸要求。定义参数β_normal = Δx_normal/δ和β_span = Δx_span/δ，分别表示壁面法向和横向的网格间距与式(8)所给热穿透深度的比值。

高频波动与低频波动之间的截止频率定义为f_c。频率高于f_c的热流密度波动由上述递归卷积法处理，而频率低于f_c的波动则通过求解波动温度的瞬态控制方程来解析。需要满足两个要求。

首先，数值解的准确性取决于网格布置，因此，必须指定足以分辨热穿透深度的最大网格尺寸。其次，温度传递函数的推导基于一维假设，当展向热扩散不可忽略时，该假设将失效。如果目标是使递归卷积结果与数值解尽可能接近，第二个要求施加了一个约束，即壁面上的展向网格尺寸不应过小。当展向网格尺寸趋近于零时，展向热扩散不可避免地变得不可忽略。因此，需要明确定义固体表面允许的最小展向网格尺寸。

假设壁面与第一层网格单元中心之间呈线性温度变化，输入热流呈角频率为 $\omega$ 的单频正弦波动，并使用有限体积法 (FVM) 对热传导方程进行离散化。随着网格尺寸的变化，数值计算得到的壁温波动与由方程 (12) 中的传递函数 $\hat{C}_{Tq}$ 推导出的波动之间的比较如图 4 所示。$\hat{A}(\hat{T}_{ref}')$ 表示通过传递函数计算得到的波动幅值。可以看出，当网格间距小于热穿透深度时，数值结果与由传递函数推导出的温度波动吻合良好；然而，当 $\beta_{normal}$ 大于 2 时，数值结果显著高估了壁温波动的幅值。

图 4. 壁面法向网格分辨率不足导致的壁温波动高估。

考虑到方程 (9) 中给出的一维频域热传导方程解析解以及温度传递函数 $\hat{C}_{Tq}$，通过一阶外推得到的壁温表达式如下：
(33) $\hat{T}_{w,num} = \hat{T}_1 + \frac{\Delta x_{normal}}{2} k \hat{q}_w = \hat{T}_{w,exact} - \frac{\sqrt{i\omega\alpha} \Delta x_{normal}}{2} + \frac{\Delta x_{normal}}{2} \sqrt{i\omega\alpha} \hat{T}_{w,exact}$
其中 $\hat{T}_{w,num}$ 是数值解，$\hat{T}_{w,exact}$ 是精确解，$\hat{T}_1$ 是第一层网格单元中心的温度，$\Delta x_{normal}$ 是壁面法向上第一层网格的尺寸。定义 $W = |\hat{T}_{w,num} / \hat{T}_{w,exact}|$ 为由壁面法向网格分辨率引起的壁温波动幅值的高估因子。则 $W$ 的表达式如下：
(34) $W = |e^{-\sqrt{i\omega\alpha} \frac{\Delta x_{normal}}{2}} + \frac{\Delta x_{normal}}{2} \sqrt{i\omega\alpha}| = |e^{-\sqrt{2i} \frac{\beta_{normal}}{2}} + \sqrt{2i} \frac{\beta_{normal}}{2}|$
图 5 说明了由数值计算和解析表达式得到的壁温幅值高估因子与参数 $\beta_{normal}$ 之间的关系。可以观察到两者吻合良好。可以看出，当采用一阶外推法重构壁温时，采用网格间距准则 $\beta_{normal} < 1.1$ 可确保壁温幅值的误差保持在 10% 以下，而 $\beta_{normal} < 0.9$ 可保证误差小于 5%。

图 5. 随 $\beta_{normal}$ 变化的壁温波动幅值高估因子：数值结果与解析解的比较。

上述温度传递函数基于一维假设推导得出，其在三维共轭传热 (CHT) 问题中的适用范围需要澄清。鉴于难以获得其三维形式的解析解，这部分工作通过数值计算进行。计算域为一个三维立方体，在上表面的单个单元面上施加波动热流，而其他所有边界均为绝热。立方体足够大，以消除边界效应对传热的影响。假设施加在网格上的波动热流表示为 $\hat{q}_{ref}$。通过壁面传递函数计算得到的温度波动为 $\hat{T}_{ref} = \hat{q}_{ref}\hat{C}_{Tq} = \hat{q}_{ref}/(k\sqrt{i\omega/\alpha})$，热穿透深度为方程 (8) 中给出的 $\delta_P = \sqrt{2\alpha/\omega}$。以上述量为参考值，频域热传导方程及其边界条件被无量纲化为 $\hat{T}^* = \hat{T}/\hat{T}_{ref}, x^* = x/\delta_P$：
(35) $\nabla^2 \hat{T}^* - 2i\hat{T}^* = 0$
(36) $\frac{\partial \hat{T}^*}{\partial n^*} = -\sqrt{2i}$
其中 $n^*$ 是热斑网格单元壁面的法向无量纲坐标。对于上述控制方程和边界条件，在没有几何和边界效应的情况下，解仅取决于两个无量纲参数：距热斑网格的无量纲距离 ($r^*$) 和壁面处的无量纲展向热斑网格单元间距 ($\beta_{span}$)。$\beta_{span}$ 表示展向长度尺度与热穿透深度之比。随着 $\beta_{span}$ 减小，热穿透深度变得与展向长度尺度相当，允许热扰动更显著地扩散到邻近区域，使展向扩散变得越来越重要。相反，较大的 $\beta_{span}$ 对应于更局部化、近似一维的热响应。

求解了上述问题，$\hat{T}^*$ 随 $r^*$ 和 $\beta_{span}$ 变化的情况如图 6 所示。可以观察到，随着 $\beta_{span}$ 增加，数值结果趋向于传递函数预测的值。确保展向网格参数满足 $\beta_{span} > 2$ 可保持两者之间的相对误差在 4% 以下。这并不意味着当 $\beta_{span} < 2$ 时温度传递函数变得不适用。在这种条件下，传递函数的准确性与热流分布有关。如果热流分布近似局部一维，传递函数仍然有效。相反，如果热流分布表现出显著的不均匀性，传递函数可能会引入大得多的误差。

图 6. 随固体壁面展向网格尺寸变化的展向热扩散强度。

因此，对于给定的固体表面网格尺寸，可以使用准则 $\beta_{span} > 2$ 确定满足一维假设的最小频率 $f_c$。由于该准则仅取决于固体网格尺寸、材料热扩散率和频率 $f_c$，因此可以在耦合模拟之前确定，而不需要预先知道表面热流分布。随后，根据准则 $\beta_{normal} < 0.9$ 设计壁面法向的网格尺寸，使得所得到的网格分辨率足以准确捕捉频率低于 $f_c$ 的波动的热穿透深度。由于 $f_c$ 的确定与网格尺寸有关，每当固体网格尺寸发生变化时（例如涉及移动边界或旋转部件的情况），都应重新评估 $f_c$ 的值。

### 2.3. 瞬态共轭传热方法的计算流程

在本节中，前述各节描述的方法被整合为一个统一的瞬态共轭传热计算方法，并给出了详细的求解流程。提出的方法以松耦合方式实现，求解传递策略如图 7 所示。流体求解器和固体求解器中使用的时间步长不同。高频温度波动使用递归卷积法解析，其计算成本非常低，每个网格面仅需 $M$ 次代数运算。因此，它可以在每个流体时间步长进行求解，从而实现高频扰动的高分辨率捕捉。相比之下，时均温度和低频波动以数个流体时间步长为间隔进行求解，两者的求解间隔可以不同。需要注意的是，用于解析低频温度波动的时间步长必须足够小，以捕捉低频带的上限频率 $f_c$。流体求解器每个时间步长推进所获得的壁面热流经过移动平均和滤波处理，生成更新的热流 $\bar{q}, q'$ 和 $q''$ 作为固体求解器的边界条件。反过来，固体求解器的每次响应产生更新的壁温 $T = \bar{\bar{T}} + T' + T''$，作为流体求解器的边界条件。提出方法的计算流程如图 8 所示。

图 7. 瞬态共轭传热方法中的求解传递策略。
图 8. 瞬态共轭传热方法的计算流程。

### 3. 结果与讨论

在本节中，使用四个选定的算例对所提出的多尺度时域共轭传热计算方法进行评估。第一个算例是单介质固体中的一维非稳态传热。第二个算例是具有热障涂层 (TBC) 的固体中的一维传热。第三个算例是涡轮叶片中的三维传热。第四个算例是圆柱绕流的非稳态共轭传热问题。前两个算例的示意图如图 9 所示。由于与流体求解器耦合会引入大量额外的计算成本，且在受控数值算例中流体求解器的主要作用是提供波动热流载荷，因此前三个算例采用预置的波动热流边界条件。对于这三个算例，采用预置的波动热流边界条件来表示具有不同频谱特性的热载荷，包括离散频率和宽带扰动，同时避免了显式求解流体域的额外计算成本。

图 9. 一维瞬态传热问题示意图：(a) 单介质固体；(b) 带热障涂层的固体。

#### 3.1. 单介质固体中的一维非稳态传热

一维固体长度为 $L=5$ mm，由钛制成，热物性如表 1 所列。在热侧施加波动热流 $q_w$，而冷侧保持固定温度 $T_c=0$。热流 $q_w$ 由时均分量 $\bar{q}$、低频波动分量 $q'$ 和高频波动分量 $q''$ 组成。时均热流为 $1.2\times10^5$ W/m$^2$。低频波动由 20 Hz 的单频波动组成，幅值为 $1\times10^5$ W/m$^2$。高频波动由叠加频率为 500、1000、5000 和 10000 Hz 的波动构成，每个分量的幅值均为 $3\times10^5$ W/m$^2$，且每个分量的相位随机。低频和高频之间的截止频率 $f_c$ 设置为 50 Hz，拟合阶数 $M$ 为 16。流体时间步长 $\Delta t_f$ 假设为 $1\times10^{-5}$ s。时均温度 $\bar{\bar{T}}$ 每 100 个流体时间步长更新一次，低频波动温度 $T'$ 每 200 个流体时间步长更新一次，确保在低频范围内，每个波动周期至少由 10 个时间步长解析。

滤波过程获得的时均和波动热流均方根 (RMS) 值的时域演化如图 10 所示。可以看出，仅在 0.1 s 内，移动平均法和巴特沃斯滤波器就能解析出相当准确的时均热流和低频波动热流。通过与非同步方法和同步方法的比较，展示了所提出方法在收敛速度方面的改进。固体表面温度的收敛历史如图 11 所示。收敛时间定义为时均温度与收敛解之间的相对误差降至 1% 以下所需的时间。基于此定义，表 3 列出了四种方法收敛解的收敛时间和温度波动的均方根值。可以观察到，本研究构建的多尺度方法的收敛时间仅为同步方法的 2%。这里，收敛时间指的是达到热收敛所需的耦合时间步长对应的模拟时间。此外，其收敛时间比 $\gamma=10$ 的非同步方法更短，且仅比 $\gamma=100$ 的方法稍长。同时，多尺度方法获得的热波动均方根值与同步结果的偏差仅为 4.2%，这远小于 $\gamma=10$ 的非同步情况的 217% 偏差和 $\gamma=100$ 情况的 903% 偏差。

图 10. 不同时间尺度下分解热流分量的均方根值。
图 11. 固体表面温度的收敛历史。

表 3. 四种方法的收敛时间和温度波动均方根值。

| 方法 | 同步 | 非同步 $\gamma=10$ | 非同步 $\gamma=100$ | 多尺度方法 |
| :--- | :--- | :--- | :--- | :--- |
| 收敛时间 (s) | 16.44 | 41.32 | 0.12 | 0.33 |
| RMS 温度 (K) | 1.834 | 5.818 | 18.404 | 1.756 |

不同位置的收敛温度如图 12 所示。可以看出，固体表面和固体内部的温度波动都能以合理的精度计算出来。

关于表面温度，如式（29）所示，拟合区域内误差的上界可以直接利用导数分布和最大模定理获得。对于固体内部区域，传递函数具有更复杂的表达形式，难以直接获得拟合区域内误差的上界。然而，最大模定理仍然适用，通过增加拟合点的数量同样可以控制拟合区域内误差。下载：下载高分辨率图像（250KB）下载：下载全尺寸图像图12.（a）x = 0处和（b）x = 0.1 mm处的收敛温度波动。在实际应用中，可能会发生突发故障，从而向流场引入扰动。这些扰动是随机的且非周期性的，因此有必要评估所提出的多尺度方法处理此类扰动的能力。保持时间平均热通量和低频波动热通量与本案例一致，高频波动设定如下：(37) ??″(??) = {|??| sin(2π???? + ??) · 1/2 [1 + cos(π/2 · (?? ? ????)/????)], |?? ? ????| ≤ ????/2; 0, |?? ? ????| > ????/2 其中????为施加扰动时间段的中点，????为扰动的持续时间。热通量扰动的分布和计算得到的温度如图13所示。由于拟合区域是沿??轴和??轴定义的，拟合函数能够处理非周期性扰动，并能准确计算非周期性温度波动中的尖峰。下载：下载高分辨率图像（212KB）下载：下载全尺寸图像图13.（a）热通量扰动的分布。（b）计算得到的温度扰动。在高保真计算中，湍流对壁面热通量的影响通常是高频且宽带的。因此，应检验本方法计算此类宽带扰动的能力。低频波动热通量和时间平均热通量保持与上述计算一致。高频热通量通过叠加频率范围150–1000 Hz内的波动来构建，频率间隔d??为0.1 Hz。各频率下热通量波动的幅值如图14(a)所示，相位设为随机值。计算得到的温度分布如图14(b)所示。可以观察到，多尺度方法获得的结果与同步方法的结果吻合良好，证明了所提出方法在高频宽带问题中的适用性。下载：下载高分辨率图像（230KB）下载：下载全尺寸图像图14.（a）热通量幅值的频域分布。（b）宽带扰动下的温度分布对比。3.2 带热障涂层的固体一维非稳态传热 在燃气轮机热部件（如燃烧室和涡轮叶片）的制造过程中，通常会在金属基体上施加热障涂层以防止壁面烧蚀。本算例考察所提出方法在带有热障涂层的多种材料热传导问题中的适用性。Kapton被用作热障涂层材料，其热物性参数列于表1。热障涂层厚度????????为0.5 mm。所有边界条件与前序算例保持一致。针对多种材料热传导问题，He [21]基于一维半无限假设建立了传递函数：(38) ???????,?????? = ??1??1[??^(????1(1+??)) ? ??^(?????1(1???))] / [??^(????1(1+??)) + ??^(?????1(1???))] ??1 = √(??/??1), ??2 = √(??/??2), ?? = ??2??2 / ??1??1 其中下标1表示热障涂层层的参数，下标2表示金属基体的参数。采用同步方法、不同步方法以及本研究提出的多尺度方法计算带有热障涂层的固体中的瞬态热传导。固体表面温度的收敛历史如图15所示，收敛时间及最终波动温度的均方根值列于表4。可以观察到，所提出方法在0.61 s内达到收敛状态，仅为同步方法所需时间的6%，为?? = 10时不同步方法所需时间的53.5%，但长于?? = 100时的收敛时间。然而，所提出方法的均方根误差仅为0.7%，远小于?? = 10时不同步方法的206%误差和?? = 100时的938%误差。下载：下载高分辨率图像（203KB）下载：下载全尺寸图像图15. 带热障涂层固体表面温度的收敛历史。表4. 带热障涂层固体中四种方法的收敛时间与温度波动均方根值。 同步方法 不同步??=10 不同步??=100 多尺度方法 收敛时间(s) 10.13 1.14 0.15 0.61 均方根温度(K) 16.73 51.19 173.71 16.84 收敛后获得的固体表面温度和深度0.01 mm处的温度如图16所示。选择如此浅的位置是因为热障涂层材料的热扩散率非常小，导致高频波动幅值迅速衰减，在极短距离内衰减至零。可以看到，同步方法和多尺度方法在两处的温度分布计算结果吻合良好，证明了所提出方法在多种材料传热问题中的适用性。下载：下载高分辨率图像（262KB）下载：下载全尺寸图像图16. 带热障涂层固体中（a）x = 0处和（b）x = 0.01 mm处的收敛温度波动。3.3 涡轮叶片非稳态传热 前述算例均为一维的，未考虑三维情况下的几何复杂性或可能存在的展向扩散现象。本节以涡轮叶片为计算对象，考察所提出方法在三维问题中的适用性。叶片材料为钛，几何与边界条件设置如图17(a)所示。等温边界固定为600 K。下载：下载高分辨率图像（132KB）下载：下载全尺寸图像图17.（a）涡轮叶片网格与边界条件设置，以及（b）热通量放大系数的分布。在热通量边界上施加空间非均匀的宽带热通量。时间平均热通量为4×10? W/m2，波动分量覆盖30至300 Hz的频率范围，频率间隔为1 Hz。为表征类湍流的宽带分布，波动热通量的功率谱密度遵循??^(?5/3)依赖关系，并为各个频率分量赋予确定性随机相位。波动热通量的总均方根值为4×10? W/m2。上述热通量乘以依赖于几何位置的热通量放大系数，该放大系数的表达式如图17(b)给出。取???? = 80 Hz作为高低频之间的截止频率，该频率下的热渗透深度为0.106 mm。壁面上的展向网格间距设为0.5 mm，壁面法向网格间距设为0.06 mm，网格尺度满足??_normal < 0.9和??_span > 2。选取叶片前缘处的一点作为监测点，以评估所提出方法的收敛加速效果和精度，报告的收敛时间基于该监测点的时间平均温度定义。监测点温度的收敛历史如图18(a)所示。与一维算例一致，所提出方法在三维构型中展现出更快的收敛速度，仅在0.23 s内实现时间平均温度的收敛，而同步方法需要2.86 s。多尺度方法的收敛时间仅为同步方法的8.04%。收敛后的温度分布如图18(b)所示。两种方法获得的涡轮叶片表面波动温度均方根值如图19所示。两种结果的分布趋势总体吻合良好。进一步对表面温度波动均方根值进行逐点对比，结果如图20所示。多尺度与同步解之间的相对误差为2.86%，最大绝对差为0.49 K，表明在整个加热表面上吻合密切。这些结果证明了本研究建立的多尺度时域计算方法在三维问题中的适用性。下载：下载高分辨率图像（253KB）下载：下载全尺寸图像图18.（a）监测点温度波动的收敛历史，以及（b）监测点收敛温度分布。下载：下载高分辨率图像（144KB）下载：下载全尺寸图像图19. 同步方法和多尺度方法获得的涡轮叶片表面波动温度均方根值。下载：下载高分辨率图像（121KB）下载：下载全尺寸图像图20. 同步方法和多尺度方法获得的表面温度波动均方根值对比。3.4 圆柱绕流的非稳态共轭传热 所提出多尺度方法的加速性能已在前述涉及一维和三维固体热传导问题的数值算例中得到验证。本小节针对圆柱绕流传热的三维问题进行共轭传热计算。流体域采用LES方法求解，主要目的在于考察所提出的多尺度CHT方法获得的结果是否与传统同步耦合方案的结果一致。计算构型如图21所示。模拟在雷诺数Re = 3920下进行。圆柱内表面作为等温加热壁面。入口温度为300 K，壁面温度为310 K。圆柱尾迹区域的周期性涡脱落引发波动热通量。采用WALE（壁面适应局部涡黏度）模型作为亚格子模型。将计算得到的圆柱表面时间平均压差系数与实验测量值[43]进行对比，结果如图22所示，两者吻合良好。该一致性验证了LES模拟在解析近壁流行为方面的准确性。下载：下载高分辨率图像（70KB）下载：下载全尺寸图像图21. 计算域构型与边界条件。下载：下载高分辨率图像（133KB）下载：下载全尺寸图像图22. 圆柱表面时间平均压差系数与实验测量值[43]的对比。在评估波动CHT响应之前，首先利用多尺度过程建立时间平均固体温度场。时间平均壁面热通量使用式(6)计算，相应的时间平均固体温度通过伪时间步长推进。在此阶段，波动温度分量尚未求解。时间平均固体温度在LES计算大约经历15个物理涡脱落后达到充分收敛状态。温度随后通过??_nd = (T ? T_in)/(T_heat ? T_in)进行无量纲化，得到的无量纲温度分布如图23所示。下载：下载高分辨率图像（64KB）下载：下载全尺寸图像图23. 瞬时流体与时间平均固体无量纲温度的分布。在时间平均固体温度场收敛之后，激活波动温度计算。传统同步CHT计算与所提出的多尺度CHT计算均从同一收敛的时间平均温度场初始化。因此，该阶段的对比旨在评估两种耦合方案所预测的波动热响应的一致性。为捕捉更精细的流场特征，本算例中壁面上采用相对精细的展向网格。根据??_span > 2的准则，确定分离高频与低频分量的截止频率为850 Hz。波动热通量通过每个时间步长上瞬时热通量减去时间平均热通量获得。然后施加滤波程序将热通量分解为低频和高频分量。低频热通量输入热传导求解器，而高频热通量采用递归卷积方法处理。在多尺度耦合算例中，时间平均温度和低频波动温度的配置使得在截止频率处的波动，每个周期进行十次耦合。对于采用CFL数0.9确定的全局流体时间步，该配置对应每18个流体时间步进行一次流体-固体耦合。在每个耦合时间步，固体求解器求解低频温度波动方程以获得波动温度的增量。该增量随后被等分为若干部分，在随后的非耦合时间步中逐步添加。因此，固体内的低频波动温度在非耦合区间内得以持续更新。同步方法和多尺度方法获得的固体表面温度波动的无量纲均方根值分布如图24所示。

可以观察到，多尺度方法的结果与同步方法的结果高度一致，两种方法之间表面温度波动均方根的最大相对误差仅为1.45%，表明多尺度方法在本文问题中能产生时间一致的结果。下载：下载高分辨率图片（92KB）下载：下载完整尺寸图片 图24. (a)同步方法和(b)多尺度方法获得的无量纲温度波动的均方根值。4. 结论 针对宽带热波动条件下的非稳态共轭传热，本文开发了一种多尺度时域方法。固体温度被分解为时间平均、低频和高频分量。平均温度使用伪时间过程进行加速，低频波动通过瞬态热传导方程求解，高频波动则利用递归卷积壁模型在线重构。主要发现总结如下：• 利用热穿透深度建立了一种基于网格的截止频率准则。侧向热扩散可能具有重要作用的低频分量被保留在多维固体求解中，而当一维假设成立时，高频分量则使用局部壁传热函数进行处理。• 壁传热函数在预设的复频率域上进行了近似。利用最大模原理和边界导数分布来评估和控制该域内的拟合误差。所得到的递归卷积公式可在时域中直接处理周期性、宽带和非周期性热流扰动。• 该方法在保持满意的波动精度的同时，大幅减少了热收敛所需的模拟物理时间。对于单一介质情形，收敛时间从16.44s降至0.33s，温度波动均方根偏差为4.2%。对于TBC（热障涂层）情形，收敛时间从10.13s降至0.61s，均方根误差约为0.7%。• 对于三维涡轮叶片情形，收敛时间从2.86s降至0.23s，仅为同步方法所需时间的8.04%。表面温度波动均方根的点间比较结果为相对误差2.86%，最大绝对差值为0.49K。完全耦合的圆柱体情形进一步表明，表面温度波动均方根的最大相对误差仅为1.45%，证实了所提出框架的时间一致性。总体而言，该方法将多维低频热传导与高效高频递归卷积处理相结合，从而在不损失宽带壁面温度响应的情况下加速缓慢的热演化过程。其适用性取决于局部一维假设的有效性、固体网格分辨率以及截止频率、拟合范围和滤波参数的选取。文中特别包含了一个初步的完全耦合LES-CHT圆柱体算例，以验证所提出的多尺度处理满足非稳态CHT耦合的时间一致性要求。对于具有显著热反馈的强耦合构型，壁面温度响应中的误差或时间延迟可能会反馈到流体侧热流中，因此需要进一步评估耦合一致性和稳定性。未来工作将聚焦于此类强耦合LES-CHT问题、温度相关物性以及复杂多层构型。CRediT作者贡献声明 徐波：调查、形式分析、概念设计、方法学、验证、初稿撰写。蔡飞雪：调查、概念设计。郭欣欣：概念设计、调查。李小鹏：方法学、调查。姚敏：概念设计、资金获取、方法学、验证、审阅与编辑。任竹音：资金获取、概念设计、方法学、监督、审阅与编辑。徐波|蔡飞雪|郭欣欣|李小鹏|姚敏|任竹音