2. 中国空气动力研究与发展中心计算空气动力学研究所, 四川绵阳 621000
2. Computational Aerodynamics Research Institute, China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China
在航空航天工程中,由于不同的功能需求,涡轮 发动机、火箭发动机的表面涂层材料,航天器的电磁窗口、观察窗,以及宇航飞行器的表面需要采用半透明材料,甚至某些飞行器大面积防热材料(如陶瓷材料、玻璃钢)本身就是半透明材料。这类材料的传热问题不仅包含热传导过程,还包含辐射效应,其控制方程为耦合热传导/辐射的复合换热微分方程。目前,国内外对这类传热正问题的研究较多[1, 2, 3, 4],对相关的传热逆问题也开展了一些研究工作,如Daouas N利用热线温度数据辨识了半透明材料的热传导系数、热容和辐射特性参数[5];Li H Y推导了辨识二维传导/辐射方程中辐射源项的灵敏度算法[6];Park H M和Yoon T Y建立了二维传导/辐射方程中辐射参数辨识的共轭梯度法[7];刘林华、谈和平提出了一种由一侧边界出射辐射强度辨识另一侧边界稳态入射辐射热流密度的方法;还提出了一种在介质辐射特性已知条件下,由壁面入射辐射热流的测量值辨识燃烧室内三维温度场的方法[8, 9]。总体来看,这些方法大多针对的是稳态问题和辐射参数辨识问题,对测点位置选取的分析较少。本文针对一维传导/辐射耦合传热问题,首先实现了同时辨识材料热传导系数和辐射特性参数的方法,并分析了在单一测点情况下测点位置对辨识结果的影响;同时,针对实际工程应用中常遇到的表面热流时间变化历程辨识问题建立了辨识算法,分析了测点位置和测量精度影响。
1 传导/辐射耦合传热问题求解方法一维传导/辐射耦合传热系统的控制方程为:
式中:I(x,μ)为辐射强度;Ib表示黑体辐射强度;μ=cosθ为方向余弦;n为介质的折射率;σ为斯蒂芬-波尔兹曼(Stephen-Boltzmann)常数,值为5.67×10-8W/(m2K4);ka为介质的吸收系数;ks为介质的散射系数;p为散射相函数。
对式(1)、(2)进行联立求解,式(1)中的微分方程可沿用求解热传导问题的有限控制体积法[10]求解,每一时刻的辐射强度积分项作为源项通过求解式(2)获取;式(2)称为RTE(Radiative Transfer Equation)方程,当温度场已知时可采用参考文献[8]中给出的基于有限体积法的计算方法来进行求解。用上述计算方法对文献[2]中给出的如下算例进行考核计算:设有2个平行黑平板之间充满着传导辐射性材料的介质层(见图 1),平板1的温度为T1,平板2的温度为T2,平板间距为D,介质的n=1,ks=0;定义无量纲量:N1=kka/(4σT13),称为Stark(斯塔克)数,又称传导-辐射参数。当N1→∞时,此时热传导占主导地位,方程(1)、(2)解的温度为线性分布,而当N1→0时,此时热辐射占主导地位。图 2给出了此时的计算结果(图中“Calc.”)与文献[2]中结果(图中“Ref.”)的比较,可以看到2者符合较好,校核了计算方法的有效性。
2 热传导系数和辐射吸收系数的辨识热传导系数和辐射吸收系数的辨识是指传热问题(1)、(2)式中的热传导系数和辐射吸收系数值未知,需要由某些测点上的温度测量值来确定。此时采用的求解方法是将该问题转化为优化问题,即选取合适的k和ka使如下目标函数达极小:
式中:xmi表示第i个测点位置,T为该测点温度的计算值,为测量值,tf为传热时间。采用最速下降法对其进行优化计算,即
式中:符号“^”表示估计值,下标“l”和“l+1”表示迭代步数;β为步长;目标函数对参数的梯度通过差分来计算,即:
Δ为摄动量。
下面给出热传导系数和辐射参数的辨识算例,设防热层为半透明材料,其中热物性参数ρ=1000kg/m3,Cp=1000J/(kg·K),k=0.5W/(m·K),厚度L=0.01m,初始温度300K,半透明材料吸收系数ka=0.5m-1,散射系数为0.5m-1,折射率1.0,散射相函数各向同性。表面热流随时间变化的规律Q(t)已知,如图 3所示,在绝热端(x=0.01m,图 4中的“Right measurement position”) 布置温度测点,利用该测点的温度测量结果来对热传导系数和辐射参数进行辨识。利用图 3中热流可仿真计算出测点温度变化如图 5中“Tcal.”所示,首先以此计算值作为测量值,取参数辨识初值k=0.1,ka=0.1来进行辨识计算,结果如表 1所示。可以看到,辨识结果与真值符合。
k | ka | |||
均值 | 标准差 | 均值 | 标准差 | |
0 | 0.49987 | - | 0.49808 | - |
σ= 0.01 | 0.49987 | 0.00144 | 0.49973 | 0.02318 |
σ= 0.05 | 0.49955 | 0.00710 | 0.50328 | 0.11561 |
σ= 0.1 | 0.49919 | 0.01416 | 0.51288 | 0.23443 |
σ= 0.05, ka=5.0 | 0.49973 | 0.00679 | 5.01838 | 0.35744 |
接下来通过以下方式在温度计算结果上叠加测量噪声作为测量值来辨识k,ka。
v(t)是均值为零标准差σ=0.01、0.05、0.1的白噪声。
测量白噪声会导致辨识结果随机散布[11],辨识结果散布可以通过重复多次试验获得。这里选取了64组相同标准差的噪声按上面的方法生成64组测量数据来做辨识计算。图 5给出了叠加不同测量噪声后的测量值。图 6给出了k和ka随迭代步的收敛过程。表 1给出了64组辨识结果的均值及标准差,可以看到:辨识结果均值与真值符合,但测量噪声越大,辨识结果的标准差越大,两者线性相关。
表 1还给出了半透明板光学厚度增加情况下的热物性参数辨识情况。算例中半透明材料吸收系数ka=5.0(1/m);在绝热端(x=0.01m)布置温度测点,在测点温度计算值上叠加64组相对误差5%的白噪声生成温度测量值,图 5给出了其中一组测量值。与吸收系数较小的辨识结果相比较可以看到:导热系数辨识结果没有显著区别;吸收系数辨识结果的标准差变大,但是相对误差较小。
下面考虑不同测点位置对辨识结果影响。考虑ka=0.5(1/m)的情况,针对上面的算例再分别选取测点位置处于左端加热面(x=0m,图 4中的Left measurement position”)和材料内部中点(x=0.005m,图 4中的“Middle measurement position”)。并在测点温度计算结果上叠加相对误差5%的白噪声当作温度测量值。同样选取了64组不同形式相同标准差的噪声来进行分析。表 2给出了64组辨识结果的均值与标准差。一般来说温度相对辨识参数的灵敏度越大,辨识结果越好。图 7给出了用测量噪声标准差正规化后的灵敏度值,可以看到同一测点导热系数灵敏度远大于吸收系数的灵敏度,导热系数辨识结果也好于辐射吸收系数。
k | ka | |||||
均值 | 标准差 | CR界 | 均值 | 标准差 | CR界 | |
左端测点 | 0.50189 | 0.01352 | 0.01041 | 0.49473 | 0.08071 | 0.08042 |
中间测点 | 0.49918 | 0.01850 | 0.01316 | 0.50381 | 0.06724 | 0.07354 |
右端测点 | 0.49955 | 0.00710 | 0.00580 | 0.50328 | 0.11561 | 0.01070 |
同时辨识2个参数时,2个参数灵敏度之间的相关性也会影响辨识结果。Cramer-Rao界[12]通过对信息矩阵求逆来估计参数辨识结果的精度,信息矩阵对角线元素反映参数的灵敏度,非对角线元素则反映了参数之间的相关性。通过图 7可以看到左端测点辐射吸收系数灵敏度比其他测点辐射吸收系数灵敏度都要大;但是实际计算结果(见表 2)左端测点辐射吸收系数辨识结果却比中间测点差,这与表 2中Cramer-Rao界的预测是相一致的;这说明了相关性对辨识精度的影响。表 2中给出了本算例中不同测点位置对应辨识结果的Cramer-Rao界。由表中可以看到:通过重复多次数值试验获得辨识标准差与Cramer-Rao界比较一致;基于右端测点温度辨识出的热传导系数标准差最小,辐射吸收系数标准差较大;基于中部温度辨识出的辐射吸收系数标准差较小,热传导系数标准差较大。
3 表面热流辨识表面热流辨识问题是在表面热流Q(t)未知的情况下,由测量方程中的信息来辨识表面热流Q(t)。采用的求解方法同样是将其转化为优化问题,即选取合适的Q(t)使如下目标函数达极小,而传热模型式(1)、(2)则视作为对未知参数的约束。
首先利用拉格朗日乘数法将其转化为使如下目标函数达极小的无约束优化问题:
对此式分部积分再做变分可得伴随变量满足的控制方程:
而目标函数对Q的梯度为:
由此梯度出发即可进行后续的优化计算,本文中采用的是共轭梯度法,即:
该优化算法的计算停止准则为:
下面给出算例,算例的基本参数与上节相同。考虑ka=0.5m-1的情况,边界热流随时间变化的规律仍为图 3,亦即图 8中的“Exact”,测点同样先取在绝热端,利用此热流可仿真计算出测点温度变化,如图 5中的“Tcal.”所示。以此温度变化历程为实测值,用上述方法对边界热流进行辨识,结果如图 8中的虚线“Estimated(σ=0)”所示,从图中可以看到:辨识结果与热流的给定值符合一致,表明了辨识方法的正确性。
下面考虑噪声影响,测点温度计算值上分别叠加64组相对误差2%,5%,10%的白噪声生成温度测量值,分别进行辨识计算,图 8中给出了其中一组的辨识结果。若定义单组辨识结果与真值的偏差及多组结果的平均值为:
平均值。
则该平均值可作为对应这一噪声标准差水平的辨识结果偏差值。对于本算例,相对误差2%、5%、10%的辨识结果平均偏差值分别为0.1046、0.2068和0.2933。图 8中还给出了提高测量数据采样频率(即测量数据时间间隔由1s减少至0.2s)的辨识结果。其相对误差5%辨识结果平均偏差值为0.1882。通过对比辨识结果可以发现:测量噪声越大辨识误差越大;提高测量数据采样频率可以减小辨识误差。
接下来分析测点位置的影响,进一步将测点选择在距左端加热面0.0025m处、中点即x=0.005m处。在测点温度计算值上均叠加64组相对误差5%的白噪声生成温度测量值,图 9给出了2个测点其中1组测量值。用前面介绍的方法对边界热流进行辨识,可计算出x=0.005、0.0025m的平均相对偏差值Em分别为0.1731和0.1156,其中1组辨识结果如图 10所示,并与测点取在绝热端(x=0.01m)的辨识结果进行了对比,从结果中可看到:与只考虑热传导的情况类似,测点越靠近加热面,辨识结果的精度越高。
下文分析半透明板光学厚度增加的情况下的表面热 流辨识,考虑半透明材料吸收系数ka=5.0m-1的情况。 在测点温度计算值上均叠加64组相对误差5%的白噪声生成温度测量值,图 9给出了其中1组测量值。其他计算条件与上面保持一致,可计算出x=0.01、0.005、0.0025m的平均相对偏差值Em分别为0.2090557、0.1582443和0.1285126。与吸收系数较小的辨识结果相比较可以看到:x=0.01m的辨识结果没有显著变化;中间测点辨识结果稍好;0.0025m的测点结果稍差。
4 小 结本文针对耦合传导/辐射的传热问题开展了相关的辨识问题研究,首先建立了耦合传导/辐射的传热问题的数值方法,进行了算例验证。然后建立了同时辨识材料热物性参数和辐射特性参数的算法和耦合传导/辐射后的非稳态表面热流密度辨识算法;分析了测点位置、测量误差、吸收系数对辨识结果的影响,结果表明:同时辨识材料热物性参数和辐射吸收系数时,利用右端测点的温度辨识出热传导系数值的精度较高,利用中间测点温度辨识出的辐射吸收特性系数精度较高。而在较大的吸收系数情况下,辨识的所得吸收系数相对误差较小。对于表面热流的辨识,测点越靠近加热面,辨识结果的精度越高;吸收系数对中间测点的辨识结果影响大于其他位置测点的辨识结果。
[1] | 范绪箕.气动加热与热防护系统[M].北京:科学出版社, 2004. |
[2] | 姜贵庆,刘连元.高速气流传热与烧蚀热防护[M].北京:国防工业出版社, 2003. |
[3] | Modest M F. Radiative heat transfer[M]. San Diego:Academic Press, 2003. |
[4] | Asllanaj F, Milandri A, Jeandel G, et al. A finite difference solution of nonlinear system s of radiative-conductive heat transfer equations[J]. International Journal for Numerical Methods in Engineering, 2002, 54:1649-1668. |
[5] | Daouas N, Fguiri A, Radhouani M S. Solution of a coupled inverse heat conduction-radiation problem for the study of radiation effects on the transient hot wire measurements[J]. Experimental Thermal and Fluid Science, 2008, 32:1766-1778. |
[6] | Li H Y. Inverse radiation problem in two-dimensional rectangular media[J]. Journal of Thermophysics and Heat Transfer, 1997, 11(4):556-561. |
[7] | Park H M, Yoon T Y. Solution of the inverse radiation problem using a conjugate gradient method[J]. International Journal of Heat and Mass Transfer, 2000, 43:1767-1776. |
[8] | 刘林华,谈和平,余其铮.半透明平板边界入射辐射热流密度的反问题[J].计算物理, 1999, 16(3):235-242. Liu Linhua, Tan Heping, Yu Qizheng. Inverse problem of boundary incident radiation heat flux in semitransparent planar slab[J]. Chinese Journal of Computation Physics, 1999, 16(3):235-242. |
[9] | 刘林华,谈和平,余其铮.燃烧室内三维温度场的辐射反问题[J].燃烧科学与技术, 1999, 5(1):62-69. Liu Linhua, Tan Heping, Yu Qizheng. Inverse radiation problem of three dimensional temperature field in furnace[J]. Journal of Combustion Science and Technology, 1999, 5(1):62-69. |
[10] | Beck J V, Blackwell B, Jr Clair C R S. Inverse heat conduction ill-posed problems[M]. John Wiley&Sons. New York, 1985. |
[11] | 钱炜祺,周宇,邵元培,等.表面热流可辨识性初步分析[J].实验流体力学, 2013, 27(4):17-22. Qian Weiqi, Zhou Yu, Shao Yuanpei, et al. A preliminary analysis of identifiability of surface heat flux[J]. Journal of Experiments in Fluid Mechanics, 2013, 27(4):17-22. |
[12] | 蔡金狮,汪清,王文正.飞行器系统辨识学[M].北京:国防工业出版社,北京, 2003. |