1.本发明涉及传热学反演技术领域，尤其涉及一种瞬态分布热负荷交替方向加权模糊反演方法。
背景技术：2.传热学反问题(inverse heat transfer problem,ihtp)是指根据传热系统内部或表面的部分温度信息反求系统的某些未知特征参量，如边界条件、热物性参数、几何形状、初始条件以及源项等。传热学反问题在航天、动力工程、材料加工处理、生物工程和无损检测等工程领域具有广泛的应用。
3.传热学反问题通常是在hadamard意义下的不适定(ill-posed)问题，即解的存在性、唯一性和稳定性三个条件不能同时满足。由于反问题的不适定性，对于经典的反演算法如共轭梯度法(conjugate gradient method,cgm)、levenberg-marquard方法(l-mm)，当测点数目减少或测量误差增加时，反演结果常会严重恶化。
4.模糊推理基于模糊理论，对输入信息具有明显的抗干扰能力，推理过程具有良好的鲁棒性和容错能力；能够有效利用不精确、不确定和不完备信息进行推理和决策；能够综合运用定性知识(包括经验知识)和定量知识完成推理过程，并具有较低的计算成本。模糊推理方法的这些特点无疑能够为解决包括传热学反问题在内的不适定问题提供实质性的帮助。尽管分散模糊推理方法在传热反问题中已有研究，然而这种方法现在还没有有效地应用于求解非稳态传热反问题。应用分散模糊推理方法求解非稳态传热反问题的难点在于，针对非稳态反问题构造恰当的综合加权方式，不合适的加权方式，不仅会严重影响收敛速度，甚至导致反演结果发散。
技术实现要素：5.为了解决背景技术中的现有技术问题，本发明的目的是提供一种能够明显降低对测点数目的要求，对边界热负荷的初始猜测值不敏感，对测量误差具有良好的抗干扰能力的瞬态分布热负荷反演方法。
6.为了实现上述目的，本发明采用了如下技术方案：
7.一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，包括以下步骤：
8.步骤1:建立被测对象的传热模型，在被测表面布置nm个温度测量点，对瞬态分布热负荷检测问题进行初始化，给定边界热负荷的初始猜测值矩阵设定反演过程迭代次数标识的初值为n＝0；
9.步骤2:根据的猜测值求解传热正问题获得测点处的温度计算值矩阵计算温度误差矩阵判断en是否满足迭代停止条件，如果满足则迭代停止，则输出分布热负荷；否则，进入步骤3；
10.步骤3:将en作为模糊推理模块fim的输入，进行分散模糊推理，得到模糊推理结果矩阵δu
t
；
11.步骤4:对δu
t
进行时间方向加权综合，获得热负荷的补偿矩阵δq
t
′
，通过插值函数进行空间方向插值，热负荷补偿量变为δq
t
；
12.步骤5:根据式对待反演点热负荷的当前猜测值进行更新，获得仅考虑了时间加权的热流估计值
13.步骤6:基于求解传热正问题获得测点处的温度计算值矩阵并计算温度误差矩阵
14.步骤7:将e
n+1/2
作为模糊推理模块fim的输入，进行模糊推理，对模糊推理结果矩阵δus进行空间方向加权综合，获得热流的补偿矩阵δqs；
15.步骤8:根据式对热流的当前猜测值矩阵进行刷新，获得既考虑了时间加权又考虑了空间加权的热流估计值
16.步骤9:置n＝n+1；返回步骤2继续循环直到满足其中停机准则。
17.进一步地，步骤1中，所述传热模型的数学描述如下：
[0018][0019]
t(x,y,t)＝t
0 0≤x≤l
x
,0≤y≤ly,t＝0
[0020][0021][0022][0023][0024]
其中，x＝(x,y)为空间坐标向量；t为温度；ρ为密度；c
p
为比热容；λ为导热系数；h为对流换热系数；t
amb
为环境温度；t0为初始温度分布；tf为仿真结束时间；
[0025]
在传热模型的(xm,ym)处布置nm个温度传感器，测量的温度为y(xm,ym,tk)(m＝1,2,
…
,nm)，反演的未知边界热负荷分布q(x,t)，热负荷离散表示为：
[0026][0027]
进一步地，步骤2中，
[0028]
[0029][0030]
其中表示边界热流为时，根据传热正问题计算获得的计算温度矩阵；
[0031]
所述迭代停止的标准如下：
[0032][0033]
式中的ε为给定的小正数，根据偏差原理，取ε＝nmn
t
σ2；当不考虑测量误差时(即当σ＝0时)，取ε＝0.001℃2，σ为测量误差的标准差的估计值。
[0034]
进一步地，步骤3中，所述模糊推理模块fim使用nm×nt
个分散的模糊推理单元fiui(i＝1,2,...,nm×nt
)，用于实现由温度误差矩阵中任意元素ei到模糊推理结果矩阵任意元素δui的推理过程；
[0035]
取ei和δui的模糊集论域分别为[-pe,pe]和[-pu,pu]，将ei和δui论域分别划分为7个模糊集合{a1,a2,
…
,a7}和{b1,b2,
…
,b7}，这些模糊集合对应的语言值分别为：nb(负大)，nm(负中)，ns(负小)，ze(零)，ps(正小)，pm(正中)，pb(正大)；采用三角形隶属度函数确定各模糊集合的隶属度；
[0036]
fiui的模糊推理规则如下：
[0037][0038]
nm×nt
个推理结果组成推理结果矩阵，即
[0039][0040]
进一步地，于步骤4中，时间方向加权综合采用卡方分布模型构造；
[0041]
卡方分布概率密度函数为：
[0042][0043]
其中nf为卡方分布的自由度，它是一个正整数；调整卡方分布自由度可以改变温度对热流的响应模式；tk＝kdt(k＝1,2...,n
t
)，为离散时间；γ(x)为伽马函数，即
[0044]
基于卡方分布的时间加权系数模型为：
[0045][0046]
归一化后可得加权系数：当i＜j时，也就是时间加权矩阵λ
t
为一个下三角矩阵；
[0047]
通过时间域的加权综合公式[δq
t
′
]
nm
×
nt
＝[δu
t
]
nm
×
nt
[λ
t
]
nt
×
nt
，产生时间域的推理补偿矩阵，由于空间上反演点个数与测量点个数不相等，所以在空间上采用插值的方法将矩阵[δq
t
′
]
nm
×
nt
转化为[δq
t
]
ne
×
nt
。
[0048]
进一步地，步骤6，主要包括：
[0049]
基于求解步骤1中建立的传热正问题模型获得测点处的温度计算值矩阵并计算温度误差矩阵如下：
[0050][0051]
进一步地，步骤7中，其中，将e
n+1/2
作为模糊推理模块fim的输入，进行模糊推理，推理方法同步骤4；
[0052]
其中，行空间方向加权综合，主要包括：
[0053]
空间域加权综合矩阵λs采用基于正态分布的综合加权模型，空间加权综合矩阵λs元素为
[0054][0055]
式中的θ＞0，为正态分布方差系数，xi和xj分别为反演点和测点在x方向的位置。
[0056]
本发明至少具备以下有益效果：
[0057]
本发明能够利用有限的局部温度测量信息，根据所建立的瞬态分布热负荷时空交替方向反演方法，对广泛存在的不易或不能直接测量的边界热负荷进行估算、检测。
[0058]
本发明在实际使用中，能够明显降低反演过程对测点数目的要求，对热负荷分布的初始猜测值不敏感，对测量温度误差具有良好的抗干扰能力。
附图说明
[0059]
图1为本发明方法的系统示意图。
[0060]
图2为本发明实施例中二维非稳态导热系统示意图。
[0061]
图3为ei对于模糊集合a
l
的隶属度。
[0062]
图4为δui对于模糊集合b
l
的隶属度。
[0063]
图5为xe＝0.065m处不同初始猜测值的本发明的热流反演结果。
[0064]
图6为xe＝0.065m处不同初始猜测值的cgm热流反演结果。
[0065]
图7为t＝105s，nm＝15时本发明方法与cgm的热流反演结果对比。
[0066]
图8为t＝105s，nm＝10时本发明方法与cgm的热流反演结果对比。
[0067]
图9为t＝105s，nm＝5时本发明方法与cgm的热流反演结果对比。
[0068]
图10为t＝105s，nm＝15，σ＝0.05℃时本发明方法与cgm的热流反演结果对比。
[0069]
图11为t＝105s，nm＝15，σ＝0.1℃时本发明方法与cgm的热流反演结果对比。
[0070]
图12为t＝105s，nm＝5，σ＝0.05℃时本发明方法与cgm的热流反演结果对比。
[0071]
图13为t＝105s，nm＝5，σ＝0.1℃时时本发明方法与cgm的热流反演结果对比。
具体实施方式
[0072]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0073]
以二维非稳态钢坯导热过程边界热流反演为例，参照图1和2，具体实施例的一种瞬态分布热负荷交替方向加权模糊反演方法包括如下步骤：
[0074]
1)建立二维边界热流反演的传热模型(图2)，在被测表面布置nm个温度测量点，对瞬态分布热负荷检测问题进行初始化，给定边界热负荷的初始猜测值矩阵设定反演过程迭代次数标识的初值为n＝0；
[0075]
其中的传热模型数学描述由控制方程、边界条件以及初始条件组成，
[0076][0077]
t(x,y,t)＝t
0 0≤x≤l
x
,0≤y≤ly,t＝0
[0078][0079][0080][0081][0082]
其中，x＝(x,y)为空间坐标向量；ρ为密度；c
p
为比热容；λ为导热系数；h为对流换热系数；t
amb
为环境温度；t0为初始温度分布；tf为仿真结束时间。
[0083]
在传热模型的(xm,ym)处布置nm个温度传感器，测量的温度为y(xm,ym,tk)(m＝1,2,
…
,nm)，反演的未知边界热负荷分布q(x,t)，热负荷离散表示为
[0084]
[0085]
其中，ne为空间反演点个数。
[0086]
2)根据的猜测值求解传热正问题获得测点处的温度计算值矩阵计算温度误差矩阵判断en是否满足迭代停止条件，如果满足则迭代停止，则输出分布热负荷；否则，进入步骤3)
[0087][0088][0089]
其中表示边界热流为时，根据传热正问题计算获得的计算温度矩阵。
[0090]
所述迭代停止的标准如下：
[0091][0092]
式中的ε为给定的小正数，根据偏差原理，取ε＝nmn
t
σ2；当不考虑测量误差时(即当σ＝0时)，取ε＝0.001℃2，σ为测量误差的标准差的估计值。
[0093]
3)将en作为模糊推理模块fim的输入，进行分散模糊推理。模糊推理模块fim使用了nm×nt
个分散的模糊推理单元fiui(i＝1,2,
…
,nm×nt
)，实现由温度误差矩阵中任意元素ei到模糊推理结果矩阵任意元素δui的推理过程。
[0094]
取ei和δui的模糊集论域分别为[-pe,pe]和[-pu,pu]。将ei和δui论域分别划分为7个模糊集合{a1,a2,
…
,a7}和{b1,b2,
…
,b7}，这些模糊集合对应的语言值分别为：nb(负大)，nm(负中)，ns(负小)，ze(零)，ps(正小)，pm(正中)，pb(正大)。采用三角形隶属度函数确定各模糊集合的隶属度γ
al
(ei)和γ
bl
(δui)，如图3和图4。fiui的模糊推理规则见表1。
[0095]
表1 fiui的模糊推理规则
[0096][0097]
nm×nt
个推理结果组成推理结果矩阵，即
[0098][0099]
4)对δu
t
进行时间方向加权综合，获得热负荷的补偿矩阵δq
′
t
，通过插值函数进行空间方向插值，热负荷补偿量变为δq
t
，时间域加权综合矩阵λ
t
采用卡方分布模型构造。
卡方分布概率密度函数为：
[0100][0101]
其中nf为卡方分布的自由度，它是一个正整数。调整卡方分布自由度可以改变温度对热流的响应模式。tk＝kdt(k＝1,2...,n
t
)，为离散时间。γ(x)为伽马函数，即
[0102]
基于卡方分布的时间加权系数模型为：
[0103][0104]
归一化后可得加权系数：当i＜j时，也就是时间加权矩阵λ
t
为一个下三角矩阵。这是因为热流只会对热流作用以后的时刻有影响，对热流作用以前的时刻没有影响。
[0105]
通过时间域的加权综合公式[δq
′
t
]
nm
×
nt
＝[δu
t
]
nm
×
nt
[λ
t
]
nt
×
nt
，产生时间域的推理补偿矩阵，由于空间上反演点个数与测量点个数不相等，所以在空间上采用插值的方法将矩阵[δq
′
t
]
nm
×
nt
转化为[δq
t
]
ne
×
nt
。
[0106]
5)根据式对待反演点热负荷的当前猜测值进行更新，获得仅考虑了时间加权的热流估计值
[0107]
6)基于求解步骤1)中建立的传热正问题模型获得测点处的温度计算值矩阵并计算温度误差矩阵
[0108][0109]
7)与步骤4)中的推理方法一样，将e
n+1/2
作为模糊推理模块fim的输入，进行模糊推理，获得模糊推理结果矩阵δus。nm×nt
个推理结果组成推理结果矩阵，即
[0110][0111]
空间域加权综合矩阵λs采用基于正态分布的综合加权模型，空间加权综合矩阵λs元素为
[0112][0113]
式中的θ＞0，为正态分布方差系数，xi和xj分别为反演点和测点在x方向的位置。
[0114]
根据式[δqs]
ne
×
nt
＝[λs]
ne
×
nm
[δus]
nm
×
nt
进行空间方向加权综合，获得热流的补偿矩阵δqs；
[0115]
8)根据式对热流的当前猜测值矩阵进行刷新，获得既考虑了时间加权又考虑了空间加权的热流估计值
[0116]
9)置n＝n+1；返回步骤2)继续循环直到满足其中停机准则。
[0117]
下面给出一个采用本发明方法实现二维非稳态钢坯边界热流的反演验证实例。并将本发明方法的反演结果与共轭梯度法(cgm)的反演结果进行了比较以验证本发明方法的有效性。
[0118]
钢坯的几何尺寸参数l
x
＝0.15m，ly＝0.075m。热物性参数ρ＝7854kg/m3，c
p
＝445j/kg
·
℃，λ＝43.5w/m
·
℃。环境温度t
amb
＝20℃，对流换热系数h＝60w/m2·
℃。初始温度t0＝20℃。仿真结束时间tf＝120s。假定钢坯上表面热流分布如下所示：
[0119]
q(x,t)＝100+50sin(πx/l
x
)sin(πt/tf)kw/m2[0120]
在实际工程问题中，由于人为或非人为因素，热电偶的温度测量结果存在不可避免的测量误差。数值模拟计算中，通常按下式计算热电偶的温度“测量值”矩阵y：
[0121]
y＝t(q)+σω
[0122]
式中，t(q)表示由准确热流通过计算正问题得到的温度精确值，σ为测量误差的标准差，ω为99％置信度对应置信区间为[-2.576,2.576]的服从标准正态分布的随机数矩阵。
[0123]
本发明方法和cgm反演热流时，采用相同的迭代停止标准。
[0124]
模糊推理单元的输入误差模糊论域为pe＝1℃，输出热流模糊论域为pu＝500w/m2，空间加权正态分布方差θ＝0.1，时间加权卡方分布的自由度nf＝5。
[0125]
取nm＝15，σ＝0℃。分别取100kw/m2和125kw/m2(i＝1,2,
…
,ne；k＝1,2,
…
,n
t
)，在不考虑测量误差的情况下，考察热流初值选取对反演结果的影响。本发明方法和cgm对应的反演结果如图5和图6。
[0126]
图5表明，对于所给的热流初值，本发明方法均能得到良好的反演效果，而对cgm，在时反演结果良好，当和时，反演结果图出现微小的偏移，(见图6)。可见，本发明的方法对初值的选取不敏感。
[0127]
取σ＝0℃。分别取温度传感器数目nm＝15、nm＝10和nm＝5，在不考虑测量误差的情况下讨论测点数目对反演的影响，并与cgm的反演结果进行对比。反演结果分别如图7-9所示。从图7可以看出，当nm＝15时，本发明方法和cgm均能得到良好的反演效果。当测点数减少至nm＝10和nm＝5时，本发明方法在nm＝10及nm＝5时仍能得到较为准确的反演结果；但cgm的反演结果逐渐恶化，测点数越少，反演结果越差。
[0128]
可见，与cgm相比，本发明方法在温度测点数目明显减少时，仍然能够获得稳定有
效的反演结果，明显降低了对温度测点数目的依赖性。
[0129]
取测点数目nm＝15，热流初值测量误差的标准差分别取为σ＝0.05℃和σ＝0.1℃，考察测量误差对反演结果的影响。图10和图11分别绘出了t＝105s时的本发明和cgm反演热流。观察可知，cgm的反演结果震荡明显强于本发明方法。
[0130]
进一步，减少测点数目至5，分别取测量误差的标准差σ＝0.05℃和σ＝0.1℃，反演结果如图12-13所示。由反演结果可知，对于本发明方法，当测点数目减少时，反演结果的平均相对误差几乎不变，但是cgm的反演误差显著增加。这也证明了本发明方法对测量数目不敏感。
[0131]
以上的反演结果说明，相比cgm，本发明方法增强了对测量误差的抗干扰能力。
[0132]
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是本发明的原理，在不脱离本发明精神和范围的前提下本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明的范围内。本发明要求的保护范围由所附的权利要求书及其等同物界定。
技术特征：1.一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，包括以下步骤：步骤1:建立被测对象的传热模型，在被测表面布置n
m
个温度测量点，对瞬态分布热负荷检测问题进行初始化，给定边界热负荷的初始猜测值矩阵设定反演过程迭代次数标识的初值为n＝0；步骤2:根据的猜测值求解传热正问题获得测点处的温度计算值矩阵计算温度误差矩阵判断e
n
是否满足迭代停止条件，如果满足则迭代停止，则输出分布热负荷；否则，进入步骤3；步骤3:将e
n
作为模糊推理模块fim的输入，进行分散模糊推理，得到模糊推理结果矩阵δu
t
；步骤4:对δu
t
进行时间方向加权综合，获得热负荷的补偿矩阵δq
′
t
，通过插值函数进行空间方向插值，热负荷补偿量变为δq
t
；步骤5:根据式对待反演点热负荷的当前猜测值进行更新，获得仅考虑了时间加权的热流估计值步骤6:基于求解传热正问题获得测点处的温度计算值矩阵并计算温度误差矩阵步骤7:将e
n+1/2
作为模糊推理模块fim的输入，进行模糊推理，对模糊推理结果矩阵δu
s
进行空间方向加权综合，获得热流的补偿矩阵δq
s
；步骤8:根据式对热流的当前猜测值矩阵进行刷新，获得既考虑了时间加权又考虑了空间加权的热流估计值步骤9:置n＝n+1；返回步骤2继续循环直到满足其中停机准则。2.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤1中，所述传热模型的数学描述如下：t(x,y,t)＝t
0 0≤x≤l
x
,0≤y≤l
y
,t＝0,t＝0,t＝0,t＝0其中，x＝(x,y)为空间坐标向量；t为温度；ρ为密度；c
p
为比热容；λ为导热系数；h为对流换热系数；t
amb
为环境温度；t0为初始温度分布；t
f
为仿真结束时间；在传热模型的(x
m
,y
m
)处布置n
m
个温度传感器，测量的温度为y(x
m
,y
m
,t
k
)(m＝1,2,
…
,
n
m
)，反演的未知边界热负荷分布q(x,t)，热负荷离散表示为：3.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤2中，于，步骤2中，其中表示边界热流为时，根据传热正问题计算获得的计算温度矩阵；所述迭代停止的标准如下：式中的ε为给定的小正数，根据偏差原理，取ε＝n
m
n
t
σ2；当不考虑测量误差时(即当σ＝0时)，取ε＝0.001℃2，σ为测量误差的标准差的估计值。4.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤3中，所述模糊推理模块fim使用n
m
×
n
t
个分散的模糊推理单元fiu
i
(i＝1,2,
…
,n
m
×
n
t
)，用于实现由温度误差矩阵中任意元素e
i
到模糊推理结果矩阵任意元素δu
i
的推理过程；取e
i
和δu
i
的模糊集论域分别为[-p
e
,p
e
]和[-p
u
,p
u
]，将e
i
和δu
i
论域分别划分为7个模糊集合{a1,a2,
…
,a7}和{b1,b2,
…
,b7}，这些模糊集合对应的语言值分别为：nb(负大)，nm(负中)，ns(负小)，ze(零)，ps(正小)，pm(正中)，pb(正大)；采用三角形隶属度函数确定各模糊集合的隶属度；fiu
i
的模糊推理规则如下：n
m
×
n
t
个推理结果组成推理结果矩阵，即
5.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤4中，时间方向加权综合采用卡方分布模型构造；卡方分布概率密度函数为：其中n
f
为卡方分布的自由度，它是一个正整数；调整卡方分布自由度可以改变温度对热流的响应模式；t
k
＝kdt(k＝1,2...,n
t
)，为离散时间；γ(x)为伽马函数，即基于卡方分布的时间加权系数模型为：归一化后可得加权系数：当i＜j时，也就是时间加权矩阵λ
t
为一个下三角矩阵；通过时间域的加权综合公式[δq
′
t
]
nm
×
nt
＝[δu
t
]
nm
×
nt
[λ
t
]
nt
×
nt
，产生时间域的推理补偿矩阵，由于空间上反演点个数与测量点个数不相等，所以在空间上采用插值的方法将矩阵[δq
′
t
]
nm
×
nt
转化为[δq
t
]
ne
×
nt
。6.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤6，主要包括：基于求解步骤1中建立的传热正问题模型获得测点处的温度计算值矩阵并计算温度误差矩阵如下：7.根据权利要求1所述的一种瞬态分布热负荷交替方向加权模糊反演方法，其特征在于，步骤7中，其中，将e
n+1/2
作为模糊推理模块fim的输入，进行模糊推理，推理方法同步骤4；其中，行空间方向加权综合，主要包括：空间域加权综合矩阵λ
s
采用基于正态分布的综合加权模型，空间加权综合矩阵λ
s
元素为
式中的θ＞0，为正态分布方差系数，x
i
和x
j
分别为反演点和测点在x方向的位置。
技术总结本发明涉及传热学反演技术领域，尤其涉及一种瞬态分布热负荷交替方向加权模糊反演方法。包括：步骤1：建立被测对象的传热模型，在被测表面布置N
技术研发人员：万世斌 熊刚 徐鹏 杨佳 万友朋
受保护的技术使用者：重庆理工大学
技术研发日：2022.03.29
技术公布日：2022/7/5
