钱炜祺,周 宇,*,邵元培

(1.空气动力学国家重点实验室,绵阳 621000;2.中国空气动力研究与发展中心 计算空气动力研究所,绵阳 621000)

0 引 言

表面热流辨识是通过材料边界或内部的温度测量值来反演表面热流,是一类典型的热传导逆问题,在航空航天、冶金制造等领域有广泛的应用。表面热流辨识是数学上的不适定问题,不满足解的稳定性,也就是说即使测点温度存在较小测量白噪声也可能导致解的“爆破”;具体来说,在频域,表面热流的辨识误差在白噪声的影响下随着频率的升高逼近无界。这给辨识精度估计带来了难度,而在工程实际的测热系统设计时,希望能够在对热环境有大致预测的情况下对其辨识结果误差给出较为准确的估计,进而根据辨识结果精度要求来确定传感器的测量位置和测量精度。文献[1]针对这一问题进行了探讨,结果显示,辨识结果误差不仅与防热材料物性参数、测量点位置、测量噪声形式、测量噪声大小有关,而且还与表面热流本身的频域特性有关。将傅立叶数作为相似参数,可以对辨识结果误差进行初步的分析。但是这一结果给出的主要还是定性分析的结论,在定量评估上还存在不足。因此,本文工作是文献[1]的进一步深入探索,首先对给定单一频率的热流辨识误差进行定量分析,然后结合频域分析方法[2-3]对两个和多个给定频率组合情况下的误差规律进行分析,初步建立热流辨识误差的定量分析方法。

1 热流辨识方法及辨识结果误差定义

对典型的一维热传导问题,采用文献[1]中的无量纲方法后,传热方程可写为:

初始条件:t＝0:T＝0

其中,xi(i＝1,P)为温度测点位置;P 为测点数目;为温度测量值;v(t)为测量噪声。通常的热传导正问题计算是给定热流Q(t),计算测点的温度响应,而热流辨识问题则是根据观测方程中的测点温度信息来确定Q(t)。

表面热流辨识方法主要有顺序函数法和共轭梯度法[4-5],本文采用共轭梯度法来对式(1)进行表面热流辨识,共轭梯度法将辨识问题转化为求合适的Q(t)使如下目标函数达极小的优化问题:

对上式做变分后可以得到伴随变量满足的伴随方程,同时可以得到目标函数对Q 的梯度,然后代入具有“停止准则”的共轭梯度法进行优化计算;“停止准则”是一种迭代正则化方法,可以有效抑制辨识结果的高频非物理振荡。具体算法和有效性验证参见文献[5,6]。图1给出了典型算例的辨识结果,此时的表面热流真值取为频率f＝2 Hz的正弦函数,即图1中的“Exact”,测点位于绝热端x＝1,算例中相应的无量纲量k 为1。“Est.(σ＝0)”表示不考虑测点温度测量误差情况下的热流辨识结果,“Est.(σ＝0.05)”表示在测点温度叠加标准差σ＝0.05Tmean白噪声情况下(图2中的“Exp.(σ＝0.05)”)的热流辨识结果,Tmean为测点温度变化平均值。从图中可以看到,在不考虑测量噪声和测量噪声标准差σ＝0.05的情况下,热流辨识结果和真值都符合;在考虑测量噪声情况下,根据热流辨识结果计算出的测点温度(图2中“Fitted(σ＝0.05)”)也与测量值符合,表明辨识算法是有效的。

图1 热流真值与辨识结果对比Fig.1 Comparison of exact and estimated value of heat flux

图2 测点温度测量值与辨识拟合值对比Fig.2 Comparison of measurement temperature

定义辨识结果与真值之间的相对误差E 为:

从图1可以看到,由于传热的延迟特性,测点温度信息对表面热流尾段值不敏感,导致热流尾段值与真值存在一定误差,因而在式(2)计算时只取中间一段来进行分析,即t1＝0.25,t2＝1.25。此外,由文献[1]分析知,在白噪声的影响下辨识结果误差同样具有一定的随机性。为克服这一影响,采用如下方法:首先给定相同均值、相同标准差的N 组白噪声;对叠加这些噪声后的测量结果分别辨识得出相对误差值Ei(i＝1,N);然后对这些相对误差值取平均,将该平均值作为对应这一噪声标准差水平的平均相对误差值,记为Em,即

2 给定频率热流的辨识结果误差建模

采用上节方法,首先针对五组不同频率(f＝1,2,5,8,10 Hz)正弦形式热流在不同测量噪声方差σ＝0.01Tmean、0.02Tmean、0.05Tmean、0.1Tmean情况下的辨识结果平均相对误差进行分析,计算结果如表1所示,图3中给出了频率10 Hz情况下不同噪声水平对应的热流辨识结果误差。

图3 频率10 Hz不同噪声水平对应的热流辨识结果误差Fig.3 Estimation error for 10 Hz heat flux with reapect to noise standard devition

表1 不同频率下对应不同噪声水平的热流辨识结果与给定值平均误差EmTable 1 Mean error of estimated heat flux with different frequencies for different noise standard devition levels

接下来以表1中的f 和σ/Tmean为输入,以Em为输出,采用Kriging响应面模型[7]进行建模,然后利用该模型对频率7 Hz热流在测量噪声水平σ＝0.02Tmean、0.04Tmean、0.06Tmean、0.1Tmean情况下的辨识结果平均相对误差值进行预测,并与实际计算结果进行对比。如图4所示,图中“Model prediction”为Kriging模型预测结果,“Calculated”为实际计算结果。从图中可以看到,模型预测结果与实际计算结果符合较好,数值略偏低。

图4 频率7 Hz热流的辨识结果误差预测结果对比Fig.4 Comparison of calculated and predicted estimation error for 7 Hz heat flux

3 组合频率热流的辨识结果误差估计

对于频率组合热流,选取如下几组热流进行分析:

国家提出的“互联网+”及“大数据”发展战略，使得传统的包装及印刷行业主动或被动地融入其中，为这个古老的加工产业带来了生机与希望。其优势主要体现在遥远的距离被拉近，且囊括了整个加工行业的方方面面。当然，实现“互联网+”及“大数据”战略转型的基础平台之一就是网络搜索引擎的应用。

采用与上一小节相同的方法来计算不同噪声水平σ＝0.01Tmean、0.02Tmean、0.05Tmean情况下的辨识结果平均相对误差值。图5、图6 给出了σ ＝0.05Tmean情况下热流b和热流d的理论值、时域辨识值结果对比和频域的功率谱对比。从频域分析结果可以看到,在低频段,辨识结果与真值较为一致,差异主要体现在高频分量上。由此可知,对于频率组合的热流,低频分量能在辨识结果得到较好地复现,高频分量是导致辨识结果与理论值出现差异的主要原因,因此,辨识结果精度可以通过最高频率热流分量与测量噪声之间的对应关系来进行大致估计。选取两种频率组合的情况进行分析,以热流b为例,将其重写为:

图5 热流b理论值、时域辨识值结果和频域功率谱对比Fig.5 Comparison of exact and estimated value in time and frequency domain for heat flux b

图6 热流d理论值、时域辨识值结果和频域功率谱对比Fig.6 Comparison of exact and estimated value in time and frequency domain for heat flux d

并记Q1和Q2的辨识值为和,则由式(1)的线性可叠加性、函数正交性和式(2)知:

式中E1、E2表示Q1和Q2的辨识结果误差,由前面分析知,E1≈0,所以上式简化为:

对于Q1和Q2,在时域无法分开,但在频域可分开,因此,由Parseval(帕塞瓦尔)定理[8]及图5、图6可知:

式中Δ为离散变换引起的频率变化量。则式(6)可另写为:

其中C 的物理含义为频域中高频分量与低频分量的能量值对比。

对于N 种频率组合(N ≥3)的情况,设f1＜f2＜…＜fN,由图6知,与两种频率组合时类似,N 种频率组合情况下的辨识结果与真值差异主要集中在最高频分量上,此时的误差值为:

图7 热流c理论值、时域辨识值结果和频域功率谱对比Fig.7 Comparison of exact and estimated value in time and frequency domain for heat flux c

表2也给出了对其余几组频率组合情况下采用上述方法的估计结果(表中“估算值”),并与采用式(3)和叠加64组同分布白噪声的热流辨识统计平均结果(表中“计算值”)进行了比较,可以看到,两组结果符合较好。

下面考虑更一般的情况,前面的分析中事先已知了热流中各分量的频率、幅值、相位等信息,且热流均值为0。而在工程实际应用中,这些信息,尤其是热流的幅值和相位信息不是十分明确,此时的估算方法需更多地借助频域分析技术。考虑图8(a)的热流,对其进行频谱分析,得出功率谱如图8(b),由功率谱知该热流的最大频率为10 Hz;由于热流均值不为0,频谱中包含低频分量;同时热流中还包含有频率为4 Hz的分量。由图中频谱曲线积分可知,低频分量对应的C1≈0.125,4 Hz分量对应的C2≈0.445。接下来还需要确定10 Hz热流分量的幅值,采用的方法是先滤掉热流中最大频率f＝10 Hz的成分,再用原始热流值减去滤波后的热流,即得对应频率10 Hz的热流分量幅值,如图9中的“Exact”示,可知,其幅值近似等于0.55。因此,当测量误差σ＝0.005Tmean时,由前面分析知,该误差将全部作用到10 Hz分量上,相当于10 Hz分量对应s＝0.005Tmean/(0.55Tm10Hz)情况,Tm10Hz值在之前已计算出Tm10Hz＝0.01457,测点平均温升Tmean的值由式(1)计算出为Tmean＝0.1414。于是,当σ＝0.005Tmean时,s＝0.0880,在图3中插值得出E3＝0.3396;代入式(9)可估算出此时对组合热流的误差为:

该值与采用式(3)和叠加64组同分布白噪声的热流辨识统计平均值0.2576符合较好。进一步考虑测量误差σ＝0.01Tmean的情况,可采用类似方法计算出s＝0.176,E3＝0.8196,E＝0.6541,与采用式(3)和叠加64组同分布白噪声的热流辨识统计平均值0.6176也符合较好,进一步验证了估计方法的有效性。

综上所述,对于一多频率组合热流,快速分析其辨识结果误差的具体步骤为:

表2 典型频率组合热流的辨识结果平均误差计算结果与分析预测结果对比Table 2 Comparison of calculated and predicted mean error for estimated multiple-frequency combined heat flux

图8 待分析热流及其功率谱Fig.8 Exemplified heat flux and its power spectrum

图9 高频热流分量的幅值计算Fig.9 Amplitude calculation of high frequency component of heat flux

1)对待分析的热流进行频谱分析,获知最大频率以及其余频率的能量占比;

2)通过滤波计算获得对应最大频率热流分量的振幅AmaxHz;

3)针对最大频率的周期热流,计算其测点平均温升TmaxHz及如图3所示的“辨识结果误差～测量误差”对应关系曲线;

4)由式(1)计算待分析热流对应的Tmean;

5)根据测量误差及TmaxHz、Tmean、AmaxHz计算s;

6)由s 值在步骤(3)得出的“辨识结果误差～测量误差”关系曲线中插值得EN值;

7)由式(9)估算出待分析热流的辨识结果误差E。

4 结 论

本文对表面热流的辨识误差进行了深入分析,首先对单一频率的热流辨识误差建立了其与热流频率和测量精度之间的响应面模型。然后重点对多个频率组合情况下的辨识结果误差进行了定量分析。结果显示,频率组合热流中的低频分量在辨识结果中能较好地复现,辨识结果误差可只考虑高频分量的辨识结果误差。于是,由辨识误差定义出发,通过频谱分析方法可计算出组合热流中高频分量的振幅及能量占比,进而建立了组合热流的辨识误差的估算方法,并通过多个算例进行了验证。这一结果一方面有较强的理论意义,揭示了频率组合热流与单一频率热流辨识误差之间的内在关联,深化了对辨识结果时域、频域特性的认识,并且不仅对传热逆问题适用,对其余领域的逆问题误差分析也有一定的参考价值。另一方面,这一结果在工程上也有实用价值,可减少“暴力”的仿真辨识计算,尤其是形成一些数据库后甚至可不进行数值计算,直接查表即可对辨识结果误差进行大致估算。

下一步,方法还需要在以下三个方面进行完善和推广,一是目前分析中热流的组合频率值较为稀疏,在频域中表现为离散谱,便于开展分析,而对于频谱范围较宽、频谱特性为连续谱的情况,需进一步挖掘辨识计算的“停止准则”与热流功率谱之间的内在关联,开展深入研究。二是温度的采样频率对辨识误差的影响。采样频率提高有利于热流辨识精度的提高,但是采样频率与热流辨识精度的定量关系还需进一步分析。三是目前分析中噪声形式只考虑了白噪声形式,对于有色噪声的情况还需要分析拓展。