基于Levenberg-Marquardt算法的材料溫度相關(guān)熱物性參數(shù)辨識(shí)

竇怡彬，陳 昂，陸云超，劉陸廣，李宗陽(yáng)

（上海機(jī)電工程研究所，上海 201109）

0 引 言

隨著導(dǎo)彈向著高速、遠(yuǎn)程的方向發(fā)展，氣動(dòng)加熱及熱防護(hù)設(shè)計(jì)已經(jīng)成為高超聲速飛行器研制過(guò)程中一項(xiàng)重要且必不可少的工作。在方案設(shè)計(jì)階段，通過(guò)對(duì)彈體結(jié)構(gòu)大面積部位開(kāi)展傳熱分析，并根據(jù)得到的溫度響應(yīng)給出初步的防熱方案。結(jié)構(gòu)的溫度響應(yīng)受材料熱物性參數(shù)影響較大，而這些參數(shù)往往難以獲得或者信息不全面、不準(zhǔn)確，只能利用類(lèi)似材料的參數(shù)近似求解。因此需要開(kāi)展物性參數(shù)辨識(shí)研究，為后續(xù)防熱設(shè)計(jì)和熱傳導(dǎo)模型修正提供依據(jù)。

熱傳導(dǎo)逆問(wèn)題（inverse heat conduction problem，IHCP），是指利用實(shí)驗(yàn)手段測(cè)得物體內(nèi)部或邊界上某點(diǎn)或某些點(diǎn)上的溫度及其隨時(shí)間的變化歷程，通過(guò)求解傳熱微分方程來(lái)反演物體邊界熱流、材料熱傳導(dǎo)系數(shù)或物體內(nèi)部的熱源分布等參數(shù)。國(guó)內(nèi)外許多學(xué)者提出了大量的逆問(wèn)題辨識(shí)方法并應(yīng)用于不同的工程領(lǐng)域。錢(qián)煒祺等提出了基于靈敏度方法和順序函數(shù)法的熱導(dǎo)率辨識(shí)和表面熱流辨識(shí)方法。何開(kāi)鋒等總結(jié)了高超聲速飛行器氣動(dòng)力/熱參數(shù)辨識(shí)國(guó)內(nèi)外研究現(xiàn)狀和辨識(shí)方法。周宇等基于最速下降法研究了耦合傳導(dǎo)/輻射情況下的一維半透明材料的傳熱系數(shù)和表面熱流辨識(shí)問(wèn)題。朱燕偉詳細(xì)研究了蜂窩增強(qiáng)低密度燒蝕材料的熱物性參數(shù)和表面熱流辨識(shí)問(wèn)題。Duda給出了瞬態(tài)多維熱傳導(dǎo)逆問(wèn)題的一般求解方法。Yang等研究了剎車(chē)盤(pán)表面熱流辨識(shí)問(wèn)題。Xie等研究了氣凝膠材料吸熱反應(yīng)參數(shù)及物理參數(shù)的反演問(wèn)題。Kamalpreet 等進(jìn)行了熱成型模具淬火加熱階段超高強(qiáng)度鋼比熱容的反演分析。

熱傳導(dǎo)逆問(wèn)題主要有兩類(lèi)求解方法：基于梯度的方法和基于隨機(jī)的方法。隨機(jī)方法主要指隨機(jī)智能優(yōu)化算法，這類(lèi)方法全局最優(yōu)求解能力強(qiáng)，但是計(jì)算量巨大且辨識(shí)參數(shù)數(shù)量不能過(guò)多。基于梯度的方法具有求解精度高、速度快等優(yōu)點(diǎn)，但是容易陷入局部最優(yōu)解。L-M 算法是一種基于梯度的算法，在熱傳導(dǎo)逆問(wèn)題中被廣泛應(yīng)用。Cui對(duì)L-M 算法中不同形式阻尼因子對(duì)計(jì)算收斂性和收斂時(shí)間進(jìn)行了研究，提出了基于殘差的阻尼因子確定方式。

綜合國(guó)內(nèi)外熱傳導(dǎo)逆問(wèn)題研究可知，在熱物性參數(shù)研究方面，大部分研究集中在熱傳導(dǎo)的反演上，對(duì)于同時(shí)反演熱導(dǎo)率和比熱容的研究相對(duì)較少。周宇基于伴隨方程方法對(duì)熱物性參數(shù)按溫度區(qū)間離散為不同的常數(shù)進(jìn)行優(yōu)化辨識(shí)，其他相關(guān)文獻(xiàn)對(duì)熱物性參數(shù)的研究往往以第二類(lèi)邊界條件為例，和工程研制中常用的石英燈加熱試驗(yàn)或者電弧風(fēng)洞試驗(yàn)中的耦合邊界條件有出入，無(wú)法直接使用。本文以L-M 算法為基礎(chǔ)，推導(dǎo)了任意截面積一維固體熱傳導(dǎo)數(shù)學(xué)模型，耦合了第二類(lèi)/第三類(lèi)/黑體輻射邊界條件，并給出了隱式求解過(guò)程。其中熱物性參數(shù)定義為隨溫度變化的多項(xiàng)式函數(shù)，給出了熱導(dǎo)率和比熱容辨識(shí)的整個(gè)過(guò)程，并通過(guò)2 個(gè)算例驗(yàn)證了本文方法的有效性。本文的方法可以適用于石英燈加熱試驗(yàn)或者電弧風(fēng)洞試驗(yàn)，拓展了溫度相關(guān)熱物性參數(shù)辨識(shí)方法的使用范圍。

1 計(jì)算模型

1.1 一維熱傳導(dǎo)數(shù)學(xué)模型

任意面積下一維固體熱傳導(dǎo)方程的數(shù)學(xué)表達(dá)式為

式中：、c、、、?、和分別為材料密度、比熱容、定壓條件、熱導(dǎo)率、內(nèi)熱源、溫度和一維方向坐標(biāo)；取0、1和2分別對(duì)應(yīng)笛卡爾坐標(biāo)系、圓柱坐標(biāo)系和球坐標(biāo)系。采用有限體積法對(duì)一維結(jié)構(gòu)進(jìn)行離散，具體如圖1所示。

圖1 控制體示意圖Fig.1 The diagram of control volume

式中：卡方函數(shù)()為實(shí)值有界函數(shù)，為待辨識(shí)參數(shù)向量；σ是溫度測(cè)量誤差；是測(cè)量誤差倒數(shù)構(gòu)成的協(xié)方差矩陣。式（17）是非線性?xún)?yōu)化問(wèn)題，需要采用迭代方法求解得到。

1.2 Levenberg-Marquardt算法

Levenberg-Marqurdt 算法被廣泛用于求解線性病態(tài)問(wèn)題，該算法在最優(yōu)點(diǎn)附近退化為Gauss-Newton法，在遠(yuǎn)離最優(yōu)點(diǎn)的時(shí)候具有最速下降法的特點(diǎn)，其具體表達(dá)式為

式中：上標(biāo)表示迭代步數(shù)；=??為雅各比矩陣；為對(duì)角矩陣；為阻尼因子；h為增量步長(zhǎng)。對(duì)角矩陣的目的是通過(guò)使得矩陣對(duì)角占優(yōu)來(lái)抑制由于逆問(wèn)題的病態(tài)特性導(dǎo)致矩陣近似奇異而引起的數(shù)值振蕩和不穩(wěn)定。不同版本L-M 方法的差別主要體現(xiàn)在對(duì)角矩陣和阻尼因子的選擇上，本文采用如下形式矩陣

圖2給出了基于L-M算法的熱物性參數(shù)逆問(wèn)題計(jì)算流程圖。

圖2 逆問(wèn)題求解流程圖Fig.2 Flowchart of solving IHCP

辨識(shí)參數(shù)的漸進(jìn)標(biāo)準(zhǔn)差可以通過(guò)協(xié)方差矩陣來(lái)估計(jì)，協(xié)方差矩陣可以用近似Hession 矩陣來(lái)表示

1.3 置信區(qū)間分析

采用F-檢驗(yàn)方法確定辨識(shí)參數(shù)的置信區(qū)間，在誤差曲面上最優(yōu)辨識(shí)參數(shù)對(duì)應(yīng)誤差面上的最小值，每次固定一個(gè)參數(shù)p，其他參數(shù)可調(diào)節(jié)。對(duì)模型采用L-M算法進(jìn)行優(yōu)化得到新的卡方函數(shù)，建立如下分布函數(shù)

2 數(shù)值仿真及結(jié)果分析

2.1 有噪聲狀態(tài)熱傳導(dǎo)率逆問(wèn)題

以TC4材料為例，通過(guò)查閱手冊(cè)得到該材料的熱傳導(dǎo)系數(shù)如表1 所示，密度=4 440 kg/m。對(duì)一塊4 mm 厚的TC4 鈦合金板材，對(duì)外表面施加150 kW/m熱流，計(jì)算試片背溫。對(duì)背溫?cái)?shù)據(jù)增加0均值白噪聲，其具體表達(dá)式為

表1 TC4熱物性參數(shù)Tab.1 Thermophysical parameters of TC4

式中：為測(cè)量誤差；為［-1，1］內(nèi)正態(tài)分布的隨機(jī)數(shù)；?為含測(cè)量噪聲的溫度。計(jì)算中取=0.01/0.03/0.05。不同的測(cè)量噪聲會(huì)導(dǎo)致辨識(shí)結(jié)果隨機(jī)散布，辨識(shí)結(jié)果散布程度可以通過(guò)多次重復(fù)試驗(yàn)獲得。對(duì)每一個(gè)取值構(gòu)造50組含有不同測(cè)量噪聲的溫度數(shù)據(jù)進(jìn)行辨識(shí)，圖3給出了其中一組溫度數(shù)據(jù)。

圖3 不同隨機(jī)測(cè)量誤差下的溫度Fig.3 Temperatures with different random measurement errors

表2 給出了不同下辨識(shí)參數(shù)的平均值，其中最大相對(duì)誤差為1.6%，該誤差較小，可以認(rèn)為參數(shù)辨識(shí)結(jié)果是無(wú)偏的。表3 給出了辨識(shí)結(jié)果的偏離程度估計(jì)，通過(guò)對(duì)辨識(shí)結(jié)果計(jì)算樣本標(biāo)準(zhǔn)差和平均漸進(jìn)標(biāo)準(zhǔn)差，反映出參數(shù)估計(jì)的不確定度。整體上平均漸進(jìn)標(biāo)準(zhǔn)差和樣本標(biāo)準(zhǔn)差趨勢(shì)相同，和具有較明顯的正相關(guān)性，隨著增大而增大。而平均值相對(duì)誤差和之間沒(méi)有非常明顯的線性相關(guān)性。圖4 和圖5 分別給出了=0.01 狀態(tài)下熱導(dǎo)率辨識(shí)結(jié)果和漸進(jìn)標(biāo)準(zhǔn)誤差分布。

圖4 熱導(dǎo)率辨識(shí)結(jié)果（ξ=0.01）Fig.4 Identified thermal conductivities（ξ=0.01）

圖5 熱導(dǎo)率辨識(shí)結(jié)果漸進(jìn)標(biāo)準(zhǔn)誤差（ξ=0.01）Fig.5 Asymptotic standard errors of identified thermal conductivities（ξ=0.01）

表2 熱導(dǎo)率辨識(shí)結(jié)果和相對(duì)誤差Tab.2 Identified thermal conductivities and errors

表3 熱導(dǎo)率辨識(shí)結(jié)果的標(biāo)準(zhǔn)差和漸進(jìn)標(biāo)準(zhǔn)差Tab.3 Standard deviations and asymptotic standard errors of identified thermal conductivities

實(shí)際測(cè)試中往往只有一組試驗(yàn)的測(cè)量數(shù)據(jù)，無(wú)法通過(guò)多次平均得到較為準(zhǔn)確的辨識(shí)結(jié)果，這時(shí)需要對(duì)辨識(shí)結(jié)果進(jìn)行置信度分析或者給出指定置信度下的參數(shù)區(qū)間。從=0.05 構(gòu)造的響應(yīng)數(shù)據(jù)中取其中一個(gè)響應(yīng)結(jié)果，采用F-檢驗(yàn)方法分析其68%和95%置信區(qū)間，并和1 倍漸進(jìn)標(biāo)準(zhǔn)差下的置信區(qū)間結(jié)果進(jìn)行比較，結(jié)果如表4 所示?？梢园l(fā)現(xiàn)：1 倍漸進(jìn)標(biāo)準(zhǔn)差下的置信區(qū)間和F-檢驗(yàn)68%置信區(qū)間較為接近。

表4 熱導(dǎo)率置信區(qū)間Tab.4 Confidence interval of identified thermal conductivities

圖6 給出了熱導(dǎo)率辨識(shí)結(jié)果和置信區(qū)間，采用F-檢驗(yàn)得到的68%、95%置信區(qū)間以及1 倍漸進(jìn)標(biāo)準(zhǔn)差得到的置信區(qū)間都將真值包括在區(qū)間范圍內(nèi)。從表4和圖6可以看出：6個(gè)熱導(dǎo)率參數(shù)中的不確定度最大，且不確定度的非線性最強(qiáng)。這主要是因?yàn)楸疚挠玫谋孀R(shí)數(shù)據(jù)為階躍響應(yīng)函數(shù)，其溫度響應(yīng)值落在和對(duì)應(yīng)的溫度區(qū)間范圍內(nèi)的情況較少導(dǎo)致辨識(shí)結(jié)果的不確定度較大。F-檢驗(yàn)和漸進(jìn)標(biāo)準(zhǔn)差都可以反映出辨識(shí)參數(shù)的不確定度，但F-檢驗(yàn)還可以反映出不確定度的非線性程度，即置信區(qū)間關(guān)于辨識(shí)結(jié)果的不對(duì)稱(chēng)程度。

圖6 辨識(shí)結(jié)果和置信區(qū)間Fig.6 Identified results and confidence interval

2.2 有噪聲狀態(tài)熱傳導(dǎo)率和比熱容逆問(wèn)題

2.1節(jié)對(duì)熱導(dǎo)率辨識(shí)進(jìn)行了研究，本節(jié)主要考核程序同時(shí)對(duì)熱導(dǎo)率和比熱容進(jìn)行辨識(shí)的能力。仍然以TC4材料為例，對(duì)一塊4.0 mm厚的TC4鈦合金板材，在其外表面施加150 kW/m的熱流，持續(xù)時(shí)間100 s，以試片背溫作為辨識(shí)輸入?yún)?shù)，同時(shí)辨識(shí)熱導(dǎo)率和比熱容隨溫度變化關(guān)系。以常溫狀態(tài)辨識(shí)得到的材料參數(shù)作為計(jì)算初始條件，常溫狀態(tài)下辨識(shí)得到熱導(dǎo)率為8 W/（m·K），比熱容為660 J/（kg·K）。以=0.01時(shí)的50組數(shù)據(jù)進(jìn)行辨識(shí)，得到的辨識(shí)結(jié)果如表5和表6所示?？梢钥闯?，～和c這5個(gè)參數(shù)的漸進(jìn)標(biāo)準(zhǔn)差與辨識(shí)參數(shù)的比值較大（即不確定度較大），且相對(duì)誤差也較大，說(shuō)明這些參數(shù)之間的相關(guān)性較強(qiáng)，可辨識(shí)性較低，屬于不可辨識(shí)參數(shù)。

表5 熱導(dǎo)率辨識(shí)結(jié)果和誤差Tab.5 Identified thermal conductivities and errors

表6 比熱容辨識(shí)結(jié)果和誤差Tab.6 Identified specific heat capacities and errors

去掉上述5個(gè)不可辨識(shí)參數(shù)后，重新進(jìn)行計(jì)算，得到辨識(shí)結(jié)果和相對(duì)誤差如表7和表8所示。熱導(dǎo)率平均值相對(duì)誤差最大為13.7%，比熱容平均值相對(duì)誤差最大為2.22%。整體來(lái)看比熱容的辨識(shí)結(jié)果精度更高，這主要是因?yàn)榈目杀孀R(shí)性相比于比熱容參數(shù)來(lái)說(shuō)更差，對(duì)的辨識(shí)引起熱傳導(dǎo)率參數(shù)辨識(shí)結(jié)果整體誤差增大。

表7 熱導(dǎo)率辨識(shí)結(jié)果和誤差Tab.7 Identified thermal conductivities and errors

表8 比熱容辨識(shí)結(jié)果和誤差Tab.8 Identified specific heat capacities and errors

表9和表10 則給出了辨識(shí)結(jié)果的標(biāo)準(zhǔn)差和漸進(jìn)標(biāo)準(zhǔn)差，隨著數(shù)據(jù)噪聲水平的增大而增大，反映了參數(shù)的不確定度和具有較明顯的正相關(guān)性。而平均值相對(duì)誤差和之間沒(méi)有非常明顯的線性相關(guān)性，這點(diǎn)和只辨識(shí)熱導(dǎo)率得到的結(jié)論一樣。

表9 熱導(dǎo)率辨識(shí)結(jié)果的標(biāo)準(zhǔn)差和漸進(jìn)標(biāo)準(zhǔn)差Tab.9 Standard deviations and asymptotic standard errors of identified thermal conductivities

表10 比熱容辨識(shí)結(jié)果的標(biāo)準(zhǔn)差和漸進(jìn)標(biāo)準(zhǔn)差Tab.10 Standard deviations and asymptotic standard errors of identified specific heat capacities

圖7 和圖8 是噪聲水平=0.01 時(shí)的參數(shù)辨識(shí)結(jié)果和漸進(jìn)標(biāo)準(zhǔn)差分布結(jié)果（不可辨識(shí)參數(shù)賦值情況為熱導(dǎo)率賦值8 W/（m·K），比熱容賦值660 J/（kg·K））。對(duì)比真值可以看出辨識(shí)結(jié)果是有偏的。圖9 是噪聲水平=0.01 時(shí)的參數(shù)辨識(shí)結(jié)果（不可辨識(shí)參數(shù)賦真值），平均值相對(duì)誤差分別為0.25%、0.71%、0.000 4%、0.009 6%、0.006 1%、0.027 7%和0.036 3%，為無(wú)偏估計(jì)。由此可見(jiàn)，是不可辨識(shí)參數(shù)的賦值誤差導(dǎo)致了估計(jì)結(jié)果的有偏。2.1 節(jié)中也是因?yàn)榻o比熱容賦值為真值才保證了熱導(dǎo)率辨識(shí)結(jié)果是無(wú)偏的。

圖7 熱導(dǎo)率和比熱容辨識(shí)結(jié)果（ξ=0.01）Fig.7 Identified thermal conductivities and specific heat capacities（ξ=0.01）

圖8 熱導(dǎo)率和比熱容辨識(shí)結(jié)果漸進(jìn)標(biāo)準(zhǔn)誤差（ξ=0.01）Fig.8 Asymptotic standard error of identified thermal conductivities and specific heat capacities（ξ=0.01）

圖9 熱導(dǎo)率和比熱容辨識(shí)結(jié)果（ξ=0.01，不可辨識(shí)參數(shù)賦真值）Fig.9 Identified thermal conductivities and specific heat capacities（ξ=0.01，unidentified parameters assigned true value）

取=0.01時(shí)的一組數(shù)據(jù)進(jìn)行辨識(shí)，得到辨識(shí)結(jié)果并用F-檢驗(yàn)方法計(jì)算置信區(qū)間，結(jié)果如表11～12、圖10～11所示。辨識(shí)結(jié)果為有偏估計(jì)，置信區(qū)間無(wú)法將真值包括在區(qū)間范圍內(nèi)。68%置信區(qū)間和1倍漸進(jìn)標(biāo)準(zhǔn)差置信區(qū)間基本相同。同時(shí)熱導(dǎo)率參數(shù)的不確定度比比熱容的不確定度更大。但是所有參數(shù)的不確定度的非線性程度都較低，特別是比熱容參數(shù)基本沒(méi)有體現(xiàn)出非線性。

表11 熱導(dǎo)率置信區(qū)間Tab.11 Confidence interval of identified thermal conductivities

表12 比熱容置信區(qū)間Tab.12 Confidence interval of identified specific heat capacities

圖10 熱導(dǎo)率辨識(shí)結(jié)果和置信區(qū)間Fig.10 Identified results and confidence interval of thermal conductivities

圖11 比熱容辨識(shí)結(jié)果和置信區(qū)間Fig.11 Identified results and confidence interval of specific heat capacities

3 結(jié)束語(yǔ)

本文推導(dǎo)了基于耦合邊界條件的一維熱傳導(dǎo)方程，基于Levenberg-Marquardt（L-M）算法給出了一維熱傳導(dǎo)方程中熱物性參數(shù)辨識(shí)方法，并通過(guò)兩個(gè)數(shù)值算例驗(yàn)證了本文方法的可行性和準(zhǔn)確性，形成以下結(jié)論：

1）對(duì)于L-M算法辨識(shí)熱物性參數(shù)，固定不變的參數(shù)或者不可辨識(shí)參數(shù)的賦值會(huì)影響辨識(shí)結(jié)果的精度，即當(dāng)固定參數(shù)賦值為真值時(shí)，其余參數(shù)辨識(shí)結(jié)果是無(wú)偏的；當(dāng)固定參數(shù)賦值不為真值時(shí)，其余參數(shù)辨識(shí)結(jié)果是有偏的。

2）辨識(shí)結(jié)果的平均值和真值的相對(duì)誤差與噪聲水平?jīng)]有明顯的正相關(guān)性；辨識(shí)參數(shù)的樣本標(biāo)準(zhǔn)差和平均漸進(jìn)標(biāo)準(zhǔn)差隨著噪聲水平的增大而增大，即辨識(shí)結(jié)果的不確定度與噪聲水平具有明顯的正相關(guān)性。

3）只辨識(shí)熱導(dǎo)率情況下，6 個(gè)熱導(dǎo)率參數(shù)中的不確定度最大，不確定度的非線性程度最高。這主要是因?yàn)楸疚挠玫谋孀R(shí)數(shù)據(jù)為階躍響應(yīng)函數(shù)，其溫度響應(yīng)值落在和對(duì)應(yīng)的溫度區(qū)間范圍內(nèi)的情況較少導(dǎo)致辨識(shí)結(jié)果的不確定度較大。

4）對(duì)于同時(shí)辨識(shí)熱導(dǎo)率和比熱容的情況，熱導(dǎo)率和比熱容參數(shù)之間存在一定的相關(guān)性，導(dǎo)致部分參數(shù)不可辨識(shí)。

5）對(duì)于同時(shí)辨識(shí)熱導(dǎo)率和比熱容的情況，比熱容辨識(shí)結(jié)果精度高于熱導(dǎo)率辨識(shí)結(jié)果精度，這主要是因?yàn)闊釋?dǎo)率的可辨識(shí)性更差，更大的不確定性增大了參數(shù)辨識(shí)結(jié)果的誤差。

6）使用F-檢驗(yàn)方法計(jì)算辨識(shí)結(jié)果的置信區(qū)間可以有效反映出辨識(shí)結(jié)果不確定度的非線性程度，即置信區(qū)間關(guān)于辨識(shí)結(jié)果是不對(duì)稱(chēng)的，而1 倍漸進(jìn)標(biāo)準(zhǔn)差方法得到的置信區(qū)間則無(wú)法反映這種特征。

后續(xù)可繼續(xù)研究參數(shù)的可辨識(shí)性，確認(rèn)最優(yōu)的辨識(shí)參數(shù)個(gè)數(shù)；研究辨識(shí)輸入信號(hào)和參數(shù)可辨識(shí)性的關(guān)系，設(shè)計(jì)最優(yōu)辨識(shí)信號(hào)。