描述金屬氫化物反應(yīng)器中反應(yīng)-傳熱耦合現(xiàn)象的理論模型構(gòu)建
2025-04-30 00:00:00王京李瑞卿代敏楊福勝JasminaGrbovi?NOVAKOVI?吳震張早校
西安交通大學(xué)學(xué)報 2025年5期
摘要：為了完善金屬氫化物反應(yīng)器理論模型，實現(xiàn)精確描述反應(yīng)器內(nèi)部反應(yīng)過程，利用經(jīng)典的Stefan問題，描述金屬氫化物反應(yīng)器內(nèi)部的氫熱耦合現(xiàn)象，將之簡化為反應(yīng)鋒面推進(jìn)問題。通過結(jié)合熱勢理論和過増元等的熱質(zhì)理論，提出了一套新的基于熱勢理論的修正求解方法。該方法在求解耗散函數(shù)過程中，不再依照熱位移場與內(nèi)熱源的關(guān)系將熱位移場區(qū)分為兩個部分分別求導(dǎo)，而是對熱位移場整體進(jìn)行求導(dǎo)。結(jié)果發(fā)現(xiàn)，與文獻(xiàn)結(jié)果對比，熱勢理論修正方法求得的近似解在結(jié)構(gòu)上涵蓋了現(xiàn)有的理論模型，且具有較高的精度，與模擬結(jié)果的平均誤差為0.6%。研究表明，利用反應(yīng)鋒面推進(jìn)模型描述反應(yīng)器中氫熱耦合現(xiàn)象是合理的，所構(gòu)建的理論模型可以對不同工況下的反應(yīng)過程進(jìn)行合理及精確的預(yù)測。
關(guān)鍵詞：金屬氫化物反應(yīng)器；氫熱耦合；熱勢理論
中圖分類號：TK02 文獻(xiàn)標(biāo)志碼：A
DOI：10.7652/xjtuxb202505007 文章編號：0253-987X（2025）05-0064-11
Theoretical Model Building for Describing the Reaction-Heat Coupling
Phenomenon in Metal Hydride Reactors
WANG Jing1， LI Ruiqing1， DAI Min2， YANG Fusheng1，
Jasmina Grbovi? NOVAKOVIC′3， WU Zhen1， ZHANG Zaoxiao1
（1. School of Chemical Engineering and Technology， Xi’an Jiaotong University，" Xi’an 710049， China;2. Chongqing Jialing Special Equipment Company Limited， Chongqing 400032， China;
3. Vincˇa Institute of Nuclear Science， University of Belgrade， Belgrade 11000， Serbia）
Abstract：In order to improve the theoretical model of metal hydride reactors and achieve an accurate description of the reaction process of the reactor， the classical Stefan problem is utilized to describe the hydrogen-thermal coupling phenomenon inside the metal hydride reactor， simplifying it into a reaction front advance problem. By combining the thermal potential theory and the caloric theory proposed by Guo Zengyuan et al， a new modified method based on the “thermal potential” theory is introduced. In this method， when solving for the dissipation function， instead of differentiating the thermal displacement field into two parts based on its relationship with the internal heat source， the entire thermal displacement field is differentiated as a whole. As a result， it was found that compared with the literature results， the approximate solution obtained by the modified method of thermal potential theory encompassed the existing theoretical models structurally and exhibited high accuracy， with an average error of 0.6% compared to the simulation results. Research has shown that utilizing the reaction front advance model to describe hydrogen-heat coupling phenomena in the reactor is reasonable， and the proposed theoretical model can provide reasonable and accurate predictions for the reaction process under different operating conditions.
Keywords：metal hydride reactor; hydrogen thermal coupling; thermal potential
全球氣候問題日益嚴(yán)重，人們開始尋求應(yīng)用新型的可再生能源技術(shù)，高效的儲能是發(fā)展新能源過程中重要的一環(huán)。氫被視作一種具有前景的儲能介質(zhì)，因為它的能量密度高、綠色環(huán)保、獲得方式多樣。目前，儲氫技術(shù)存在多種技術(shù)方案，多種金屬或合金（即儲氫合金）具有通過可逆反應(yīng)大量吸收/釋放氫氣的能力，吸氫產(chǎn)物稱為金屬氫化物，相應(yīng)的固態(tài)儲氫技術(shù)也是備選方案之一。金屬氫化物具有較高的體積儲能密度、可長期安全儲能，因此受到眾多研究者關(guān)注。但是金屬氫化物材料顆粒導(dǎo)熱系數(shù)低（0.1～1 W/（m·K）），而反應(yīng)焓值較大（一般大于等于30 kJ/mol）［1］，在裝填材料的反應(yīng)器中常發(fā)生明顯的溫度突然增加/降低的現(xiàn)象，阻礙了反應(yīng)的進(jìn)一步進(jìn)行。金屬氫化物儲氫技術(shù)的推廣應(yīng)用受制于反應(yīng)器換熱能力，因此研究人員需要在反應(yīng)器設(shè)計中著重考慮更為高效的熱管理方案［2］，提出了許多傳熱強化的方法并應(yīng)用于金屬氫化物反應(yīng)器中。
目前，增強反應(yīng)器傳熱的手段可分為兩類：一類是增加金屬氫化物床層的導(dǎo)熱系數(shù)，例如Klein和Groll對膨脹石墨/金屬氫化物壓塊的使用［3］；另一類是增加反應(yīng)器中的傳熱面積，例如增加翅片、換熱管和夾套［4-5］。強化傳熱的技術(shù)手段已十分豐富，但在熱設(shè)計過程中尚缺乏簡單準(zhǔn)確的理論工具作為設(shè)計指導(dǎo)。Guo等［6］類比電學(xué)中的“電勢能”概念提出了新的物理量“火積”，用來衡量物體對外傳熱的能力。傳熱過程中的不可逆損失可以用火積耗散來進(jìn)行描述，對于特定傳熱問題更小的火積耗散代表著更小的不可逆損失和更高的傳熱效率，可以用來指導(dǎo)傳熱的優(yōu)化設(shè)計。Zhang等［7］首次將火積概念引入金屬氫化物反應(yīng)器熱設(shè)計。目前，描述氫化物反應(yīng)-傳熱過程的理論方法普遍是基于穩(wěn)態(tài)假設(shè)，例如可接受包絡(luò)法［8］、無量綱熱導(dǎo)法［9］和傳熱控制反應(yīng)速率法［10］。這些方法可以方便地進(jìn)行傳熱分析，并以平均溫度作為優(yōu)化目標(biāo)進(jìn)行傳熱優(yōu)化，但也存在如下不足：
（1）利用平均溫度作為目標(biāo)只能進(jìn)行傳熱優(yōu)化，并不能簡單等價于反應(yīng)優(yōu)化；
（2）基于穩(wěn)態(tài)假設(shè)的方法不能描述反應(yīng)器內(nèi)實際發(fā)生的瞬態(tài)過程，無法判斷反應(yīng)的起止時間和反應(yīng)速率的實時變化情況。
為此有研究者開發(fā)了瞬態(tài)模型，例如Marty等［11］提出的反應(yīng)鋒面模型和Wang等［12］提出的判斷床層特征尺寸的熱擴散模型。受到前人求解Stefan問題的啟發(fā)，Marty等選擇將反應(yīng)過程簡化為反應(yīng)鋒面由換熱壁面向絕熱壁面推進(jìn)的過程。雖然，Marty等并未說明自己采用何種方法得到模型結(jié)果，但是從過往研究來看其結(jié)果可以通過準(zhǔn)靜態(tài)方法方便地求解得到。準(zhǔn)靜態(tài)方法是對于偏微分方程的一種簡化方法，但是其應(yīng)用范圍是受限的。熱擴散模型假設(shè)第一階段反應(yīng)快速發(fā)生、床層溫度快速上升，第二階段反應(yīng)受到傳熱控制，具體持續(xù)時間由熱擴散時間決定。但是需要指出，對于金屬氫化物反應(yīng)過程來說，影響傳熱速率的重要因素應(yīng)該包括床層熱容和金屬氫化物材料的反應(yīng)焓，兩者共同決定了反應(yīng)區(qū)域由換熱壁面向絕熱壁面推進(jìn)單位距離時反應(yīng)器需要的熱量。反應(yīng)鋒面模型忽略了床層熱容，假定矩形反應(yīng)器反應(yīng)鋒面到換熱面的溫度分布為線性，極大簡化了推導(dǎo)過程，但是溫度分布線性假設(shè)是否合理是缺乏討論的。熱擴散模型則完全忽略了金屬氫化物反應(yīng)焓的影響，造成了極大的誤差。綜上所述，影響傳熱速率的兩個重要因素（熱容和反應(yīng)焓）沒有在現(xiàn)有模型中同時得到充分討論，忽略其中一個來簡化推導(dǎo)過程將影響預(yù)測效果，另外基于偏微分方程推導(dǎo)的結(jié)果也難以直接利用實施優(yōu)化。
眾所周知，物理問題的描述可以通過偏微分方程，也可以通過變分原理（最小作用量原理）。變分原理相較于偏微分方程能夠?qū)ξ锢韱栴}進(jìn)行整體描述，其變分積分可以是有意義的物理量，對物理過程的描述更簡潔且切中要害［13］。Biot在1954年及之后陸續(xù)發(fā)表一系列文章［14-16］，通過在傳熱學(xué)領(lǐng)域引入力學(xué)中哈密頓原理構(gòu)造了變分形式的熱傳導(dǎo)方程，得到了以廣義坐標(biāo)表示的導(dǎo)熱問題的拉格朗日方程。為此Biot引入的概念包括：熱位移場H，廣義熱力Q，熱勢E，耗散函數(shù)D等。熱勢耗散這一概念描述了物體中“傳熱能力的耗散”。以熱勢為基礎(chǔ)的變分原理研究啟發(fā)了后續(xù)一系列對于熱傳導(dǎo)過程的研究，例如熱質(zhì)理論［17］、火積理論［18］和傳熱優(yōu)化研究［19-20］。
事實上，金屬氫化物反應(yīng)器中狹窄反應(yīng)區(qū)域隨時間前進(jìn)的反應(yīng)-傳熱耦合現(xiàn)象可以看作一類受到傳熱控制的自由邊界問題，而相變過程中融化或者凝固的自由邊界問題（一般稱之為Stefan問題），已經(jīng)被廣泛討論。這類問題中存在一個或多個界面，把不同熱物性的區(qū)域分開，界面隨時間變化而運動，在界面處能量被釋放或者吸收。Stefan問題是數(shù)學(xué)和物理領(lǐng)域經(jīng)久不衰的研究對象，對于這一問題的解法也是相當(dāng)豐富的［21-23］。從金屬氫化物反應(yīng)-傳熱耦合機理來看，其反應(yīng)鋒面前進(jìn)的過程可以被視作一類廣義Stefan問題。基于熱勢概念的拉格朗日方程在求解多種實際傳熱問題（包括Stefan問題）的近似解中取得了令人滿意的效果［24］，但是用以求解金屬氫化物反應(yīng)鋒面前進(jìn)問題存在一定爭議。
金屬氫化物中存在的是一種典型的含內(nèi)熱源的自由邊界問題，該內(nèi)熱源是特殊的，在一維過程中可以近似認(rèn)為是隨反應(yīng)鋒面前進(jìn)的移動熱源。過增元團(tuán)隊的宋柏等［25］認(rèn)為Biot的熱勢理論在處理含源問題時存在缺陷。Biot將含源問題的熱位移場分為兩個部分，并認(rèn)為其中一部分與內(nèi)熱源無關(guān)，而求解耗散函數(shù)時必須從與內(nèi)熱源無關(guān)的熱位移場中求得。宋柏等對此提出了疑問，并且構(gòu)建了一個隨時間變化的內(nèi)熱源傳熱問題，并利用熱質(zhì)理論進(jìn)行了變分計算，求解結(jié)果較為符合解析解，證明Biot在熱勢理論中處理內(nèi)熱源問題的思路存在不足。值得注意的是，Biot對內(nèi)熱源的處理思路并未被其他研究者應(yīng)用到處理相變潛熱中。其他研究者在利用熱勢理論處理Stefan問題中相變潛熱的方法與宋柏等利用熱質(zhì)理論處理內(nèi)熱源問題的方法是相近的。
本文在構(gòu)建反應(yīng)模型時同時考慮反應(yīng)焓值和床層熱容，并通過討論兩者間相對大小來選擇合適的溫度分布曲線。在求解理論模型的過程中，將處理內(nèi)熱源的思路應(yīng)用到熱勢理論中，去求解含移動點熱源的金屬氫化物反應(yīng)鋒面前進(jìn)問題，并與熱勢的原始方法得到的結(jié)果進(jìn)行比較。除此之外，本文還在模型求解中應(yīng)用基于偏微分方程的準(zhǔn)靜態(tài)方法，將其結(jié)果與熱勢理論修正方法進(jìn)行比較，同時討論了床層熱容、反應(yīng)焓值和溫度分布曲線對兩種方法的影響。由于利用了熱勢理論并進(jìn)行了求解方法的修正，本模型能夠綜合考慮床層熱容、反應(yīng)焓值和溫度分布曲線等對鋒面推進(jìn)過程的影響，得到的理論結(jié)果符合模擬結(jié)果。
1 理論推導(dǎo)
金屬氫化物在反應(yīng)過程中可以分為兩個階段［12］。第一階段由初始溫度和初始壓力控制，近似呈現(xiàn)反應(yīng)器整體同步反應(yīng)的狀態(tài)，床層溫度可視為快速地均勻升溫（或降溫）至平衡溫度；第二階段，受到傳熱限制，反應(yīng)集中在部分溫度略高于（或略低于）平衡溫度的狹窄區(qū)域，隨著時間不斷向內(nèi)部延伸，如圖1所示。
出于簡化考慮，將第二階段反應(yīng)區(qū)域視作厚度無限薄的反應(yīng)鋒面，且溫度無限接近平衡溫度。此時，可以利用Biot引入的廣義坐標(biāo)滲透距離來描述反應(yīng)鋒面位置，進(jìn)而利用變分原理對反應(yīng)鋒面推進(jìn)這一導(dǎo)熱和反應(yīng)耦合過程進(jìn)行理論分析。
1.1 基本假設(shè)
反應(yīng)過程如圖2所示。本文基本假設(shè)包括：
（1）熱量僅沿垂直于換熱面的方向由換熱面?zhèn)鬟f到絕熱面；
（2）反應(yīng)速率僅受熱量傳遞速率控制，任一時刻傳遞到反應(yīng)鋒面處的熱量僅造成氫氣的釋放，不造成反應(yīng)鋒面處的溫度上升；
（3）反應(yīng)鋒面無限窄，在反應(yīng)鋒面前床層完全無反應(yīng)，反應(yīng)鋒面后床層完全反應(yīng)；
（4）反應(yīng)鋒面前床層溫度保持為第一階段結(jié)束時吸氫/放氫平衡溫度，反應(yīng)鋒面后床層溫度由平衡溫度逐漸變化為壁面溫度；
（5）忽略床層中氫氣擴散的時間，認(rèn)為床層中氫氣壓力梯度為0。
出于簡化考慮，在一開始的推導(dǎo)過程中，反應(yīng)鋒面后的床層溫度分布假設(shè)為線性，其合理性將在1.5 節(jié)中詳細(xì)討論，并在3.1節(jié)考慮更為復(fù)雜的溫度分布。在本文的推導(dǎo)過程中，認(rèn)為反應(yīng)器中只發(fā)生了吸氫反應(yīng)。
1.2 邊界條件
模型左側(cè)換熱壁面為第一類邊界條件（溫度恒定）
tw=t0+teq （1）
式中：tw為壁面溫度；t0為過余溫度，即換熱邊界處壁面溫度減去平衡溫度；teq為材料平衡溫度，該參數(shù)可由金屬氫化物材料的范特霍夫方程求得。模型右側(cè)絕熱壁
?t?xxmax=0 （2）
式中：t為溫度；x為在熱量傳遞方向上的位移量；xmax為最大床層厚度。
1.3 線性溫度分布假設(shè)下的Biot原始熱勢理論
Stefan問題描述的是一種非線性現(xiàn)象。通常研究者通過假設(shè)各種溫度曲線來解決這個問題，并使用熱平衡方程來確定界面的運動［21］。如1.1節(jié)提到的，反應(yīng)鋒面后的溫度分布假設(shè)為線性分布。在構(gòu)建金屬氫化物反應(yīng)-傳熱耦合模型的研究中，線性溫度分布假設(shè)是被采用過的一種溫度分布情況［16］。線性溫度分布假設(shè)可以被描述為
t=t01－xδ （3）
式中：δ為滲透距離，是時間的函數(shù)，代表了反應(yīng)鋒面的瞬時位置。不同時間下的溫度分布見圖3，圖中溫度曲線的轉(zhuǎn)折處的橫坐標(biāo)即為滲透距離δ1－δ2，xwall為壁面位置。吸氫過程中壁面溫度低于平衡溫度，過余溫度t0為負(fù)值，放氫過程反之。
已知Biot得到的熱傳導(dǎo)過程的拉格朗日方程如下
?E?δ+?D?=Q （4）
式中：E為熱勢，描述物體對外傳熱的能力，其作用與位勢（位能）相類似；D為在純粹的熱的課題（不考慮熱功轉(zhuǎn)換）中的耗散函數(shù)；Q為廣義熱力，是溫度t在熱流場中虛位移方向上所作的功；為滲透距離對時間的導(dǎo)數(shù)。式（4）為Biot首先提出的在傳熱過程中的變分方程，公式形式類比分析力學(xué)中的拉格朗日方程，描述了熱傳遞過程中的“熱的勢能”的積累（?E/?δ）、耗散（?D/?σ）與傳遞（Q）。為求解這一方程，首先引入熱位移場H，其定義為在任意點關(guān)于時間的導(dǎo)數(shù)為熱流密度
?H?τ=q （5）
式中：τ為時間；q為熱流密度。
帶有H的能量守恒方程可以描述為
cρt+′=－divH （6）
式中：c為床層比熱容；ρ為床層密度；′為內(nèi)熱源強度，其意義為滲透距離前進(jìn)單位長度時，床層中反應(yīng)需要吸收的熱量。能量守恒方程右手邊為熱位移場的散度，由于熱位移場為熱流在時間上的積分，該項代表此單位空間熱流的散度在時間上的積分。根據(jù)能量守恒方程可得熱位移場H為
H=12ρct0δ1－xδ2+′δ1－xδ （7）
觀察式（7）可以發(fā)現(xiàn)，熱位移場依據(jù)與內(nèi)熱源的關(guān)系分為兩個部分
H=H++H* （8）
式中：H+為與內(nèi)熱源無關(guān)的熱位移場，代表受床層熱容影響而發(fā)生改變的熱位移場；H*為內(nèi)熱源造成的熱位移場，代表受內(nèi)熱源影響而發(fā)生改變的熱位移場。據(jù)此，可知
H+=12ρct0δ1－xδ2 （9）
H*=′δ1－xδ （10）
式（10）表示在滲透距離處，熱通量等于反應(yīng)吸收的熱量。Biot認(rèn)為內(nèi)熱源造成的熱位移場H*不能對廣義坐標(biāo)即滲透距離δ求偏導(dǎo)。因此，耗散函數(shù)與無內(nèi)熱源的耗散函數(shù)產(chǎn)生區(qū)別，含內(nèi)熱源耗散函數(shù)項可以描述如下
?D?=?D+?－Q+ （11）
D+=12∫1k（H·+）2dV （12）
Q+=－∫1k?H+?δH·*dV （13）
式中：D+為由熱位移場H+引起的耗散函數(shù)；Q+為由內(nèi)熱源造成的廣義熱力。為了計算的簡單，這里需要進(jìn)行一個代換
?H+?δ=?H·+? （14）
綜合式（11）～（14），可以求得
?D?=215k ρ2c2t20+13k ρct0′δ （15）
除此之外，熱勢的定義和它的求解過程如下
E=12∫ρct2dV=12∫δ0ρct201－xδ2dx=16ρct20δ（16）
熱勢對滲透距離的偏導(dǎo)，可得
?E?δ=16ρct20 （17）
廣義熱力亦可求得
Q=∫t?H?δdS=∫t0?H?δdS=12ρct20+′t0（18）
將式（16）～（18）代入式（4），可得
δ（τ）=2k（（1/3）ρct20+q·′t0）（2/15）ρ2c2t20+（1/3）ρct0q·′τ （19）
1.4 線性溫度分布假設(shè)下的熱勢理論修正方法
宋柏等［25］基于熱質(zhì)理論［26］提出一套新的拉格朗日方程，對導(dǎo)熱問題進(jìn)行了一系列分析，其具體形式與Biot的方程存在較大差異。但是，兩者同樣利用拉格朗日方程求解內(nèi)熱源導(dǎo)熱過程，而且結(jié)構(gòu)相類似，宋柏等對耗散函數(shù)的討論仍有借鑒意義。宋柏等認(rèn)為熱位移場分為與內(nèi)熱源有關(guān)和無關(guān)的兩部分沒有必要，這一處理的具體意義不明，故對熱位移場進(jìn)行整體求導(dǎo)，求解了隨時間變化的含內(nèi)熱源導(dǎo)熱過程，取得了較好的效果。其實從本文1.3節(jié)的推導(dǎo)過程中可以發(fā)現(xiàn)，從H*的公式看，它與H+一樣包含了滲透距離δ，因此從數(shù)學(xué)上講，H*對δ進(jìn)行偏導(dǎo)計算不存在障礙。本文使用宋柏的思路對熱勢理論處理內(nèi)熱源的方式進(jìn)行修正。
熱勢理論修正方法下的熱位移場與原始熱勢理論不存在區(qū)別，但是前者需要對熱位移場進(jìn)行整體求導(dǎo)，導(dǎo)致耗散函數(shù)不同
H·=?H?τ=12ρct01－xδ2+
ρct01－xδxδ+′（20）
D=12∫1kH·2dV=
115kρ2c2t20+13kρct0′+12k′22δ（21）
進(jìn)一步，還會導(dǎo)致拉格朗日方程中耗散項出現(xiàn)差別
?D?=215kρ2c2t20+23kρct0′+1k′2δ（22）
除了耗散項之外，熱勢理論修正方法熱勢E和廣義熱力Q的方程與原始熱勢理論是一樣的。將此3項代入拉格朗日方程，最終結(jié)果如下
δ（τ）=2k（（1/3）ρct20+q·′t0）（2/15）ρ2c2t20+（2/3）ρct0q·′+q·′2τ（23）
從結(jié)果上看，主要區(qū)別在于根號下分式分母不同，熱勢理論修正方法相較于原始熱勢理論多出了13ρct0′+′2，即內(nèi)熱源產(chǎn)生的熱位移場對滲透距離δ進(jìn)行求導(dǎo)的結(jié)果。Lardner等在處理融化相變中的Stefan問題時［24］，構(gòu)建的熱位移場結(jié)構(gòu)與本文的相類似：一部分是相變造成的熱位移場，熱量用于相變潛熱；另一部分是與相變無關(guān)的熱位移場，兩部分均采取對δ求導(dǎo)的同一處理方法。Biot引用了Lardner的工作［16］并未對其處理方法提出異議。從傅里葉導(dǎo)熱定律求解問題的傳統(tǒng)方法來看，由反應(yīng)造成的熱源和潛熱造成的熱源在性質(zhì)上是相同的，在處理方法上不應(yīng)該存在區(qū)別。但是，在利用變分原理求解含內(nèi)熱源的問題時，原始熱勢理論在處理反應(yīng)造成的熱源和潛熱造成的熱源上產(chǎn)生了差異。這種處理方法的差異在內(nèi)熱源數(shù)值較大時將引起大的計算誤差，必須進(jìn)行合理性的討論。
1.5 線性溫度分布假設(shè)下的準(zhǔn)靜態(tài)方法
除了前面所述的利用變分原理近似求解金屬氫化物反應(yīng)鋒面前進(jìn)這類廣義Stefan問題，利用偏微分方程同樣可以進(jìn)行Stefan問題的近似求解。在求解問題之前，估計并假設(shè)其溫度分布狀態(tài)是同樣重要的。需要注意的是，在利用微分方程求解經(jīng)典相變Stefan問題時，通常認(rèn)為：在截面恒定的情況下，當(dāng)固體層中的顯熱與界面處吸收的潛熱相比較小時，平面固體中的溫度分布幾乎是線性的［21］。經(jīng)典Stefan問題中引入了Stefan數(shù)（Sste）來定義固體層中顯熱與潛熱的比值［21］，其公式如下
Sste=ct0h （24）
式中：h為融化潛熱。當(dāng)顯熱與潛熱相比較小（Sstelt;1）時，導(dǎo)熱微分方程中的瞬態(tài)項可忽略，傳導(dǎo)過程幾乎與穩(wěn)態(tài)相同。忽略瞬態(tài)項意味著固液混合物中溫度近似處于穩(wěn)態(tài)分布（在一維問題中代表線性分布）。這一條件下，時間的影響只能通過相變界面的緩慢移動來體現(xiàn)，Stefan問題呈現(xiàn)準(zhǔn)穩(wěn)態(tài)形式。已有研究將準(zhǔn)靜態(tài)方法應(yīng)用到含有內(nèi)熱源的Stefan問題求解中［27-28］，這一思路同樣可以應(yīng)用到金屬氫化物反應(yīng)過程中的廣義Stefan問題。類比Stefan數(shù)，本文構(gòu)建無量綱參數(shù)“廣義Stefan數(shù)”（Sgste）描述反應(yīng)鋒面前進(jìn)過程中顯熱與反應(yīng)潛熱的比值，定義如下
Sgste=ct0MMHΔH （25）
式中：ΔH為金屬氫化物反應(yīng)焓；MMH為金屬氫化物的摩爾質(zhì)量。 通過計算可以發(fā)現(xiàn)，在金屬氫化物領(lǐng)域，普遍小于1。 以LaNi5為例，取溫差t0為10℃，此時廣義Stefen數(shù)為0.06，遠(yuǎn)小于1。說明這種情況下利用線性溫度分布作為模型基本假設(shè)是合理的。
反應(yīng)發(fā)生后，滲透距離內(nèi)床層中積累的顯熱Q+為
Q+=∫ρctdV=∫δ0ρct01－xδdx=12ρct0δ（26）
滲透距離前進(jìn)dδ時，床層溫度分布發(fā)生變化，床層中增加的顯熱dQ+可以求得
dQ+=∫δ+dδ0ρct01－xδdx－∫δ0ρct01－xδdx=
12ρct0dδ（27）
同時，厚度為dδ的床層發(fā)生了反應(yīng)，這部分反應(yīng)需要的熱量為
dQ*=′dδ （28）
由傅里葉導(dǎo)熱定律可知，外界提供的熱流近似為
Q·=kdtdx=kt0δ=12ρct0+′dδdτ（29）
隨著反應(yīng)的進(jìn)行，熱流在時間上的積累等于滲透距離前進(jìn)dδ所需要的熱量總和，即
Q·dτ=dQ*+dQ+ （30）
最終可以求得跟變分方法類似的結(jié)果為
δ=2kt0（1/2）ρct0+q·′τ （31）
2 參數(shù)與模擬驗證結(jié)果
2.1 金屬氫化物反應(yīng)器模型構(gòu)建與LaNi5物性參數(shù)
LaNi5作為應(yīng)用最為廣泛的儲氫合金之一，本文選取其做為數(shù)值模型驗證的參考材料，模型由COMSOL Multiphysics 5.6構(gòu)建。本文模型的離散化方法采用向后微分公式（BDF），方程的主要參數(shù)采用二階離散化精度，離散計算域采用非結(jié)構(gòu)化三角網(wǎng)格，單元格數(shù)為1 282，模型的相對容差為10－4，方程的收斂準(zhǔn)則由相對容差控制。模型計算結(jié)果與實驗結(jié)果驗證后，說明本文使用的求解器設(shè)置是可接受的。考慮到實際的瞬態(tài)反應(yīng)過程，在制定模型時做了以下假設(shè)：①模型中忽略對流換熱與輻射換熱，床層和壁面之間的接觸熱阻為0；②金屬氫化物的熱物性保持不變且各向同性；③氫氣與金屬氫化物材料保持局部熱平衡；④忽略金屬氫化物反應(yīng)器中氫氣流動引起的壓力波動。
LaNi5材料的反應(yīng)動力學(xué)方程為
=kexp－EaRtlnp0peq（ρsat－ρs）（32）
式中：為氫氣釋放速率；k為指前因子；Ea為活化能；R為理想氣體常數(shù)；p0為吸氫壓力；peq為平衡壓力；ρsat為飽和狀態(tài)的金屬氫化物密度；ρs為當(dāng)下金屬氫化物密度。該材料的平衡壓力可由范特霍夫方程求解
lnpeqpref=a－bT （33）
式中：pref為參考壓力，取值為1 MPa；a為范特霍夫公式中常數(shù)，取值為10.7；b為范特霍夫公式中常數(shù)，取值為3 704.6 K。數(shù)值模型的其他主要控制方程等具體細(xì)節(jié)可參考團(tuán)隊已發(fā)表文獻(xiàn)［29］，LaNi5物性參數(shù)［30-31］見表1。
表1 LaNi5物性參數(shù)
Table 1 LaNi5 physical property parameter參數(shù)數(shù)值吸收速率常數(shù)A/s－159.187吸氫反應(yīng)活化能Ea/（J·mol－1）21 179.600反應(yīng)焓ΔH/（J·mol－1）30 100飽和金屬氫化物密度ρsat/（kg·m－3）8 526儲氫合金密度ρemp/（kg·m－3）8 400材料有效導(dǎo)熱系數(shù)k/（W·m－1·K－1）1.087金屬氫化物材料比定壓熱容cps/（J·kg－1·K－1）419床層孔隙率ε0.500摩爾質(zhì)量MMH/（g·mol－1）432.373
2.2 參數(shù)驗證
將表1參數(shù)代入式（19）、（23）、（31），推導(dǎo)得到的結(jié)果均可寫為如下形式
δ（τ）=βτ （34）
式中：β為Stefan問題解的系統(tǒng)。
可見，式（34）符合傳統(tǒng)Stefan問題中解的一般形式。當(dāng)t0分別為10、20、40 K時，討論3種方法得到的β，具體結(jié)果可以參考表2。
隨著t0增加，熱勢理論修正方法和準(zhǔn)靜態(tài)方法得到的β均逐漸增大，這一趨勢符合物理規(guī)律：溫差t0越大，推動鋒面前進(jìn)的勢差越大，鋒面前進(jìn)速率越快，表現(xiàn)為β越大。原始熱勢理論得到的結(jié)果與此趨勢相反，證明其不符合物理機理，是錯誤的。從數(shù)值上看，熱勢理論修正方法與準(zhǔn)靜態(tài)方法得到的結(jié)果接近，且均與原始熱勢理論結(jié)果相差巨大，幾乎達(dá)到一個數(shù)量級的差別。通過數(shù)值模擬金屬氫化物吸氫過程，選取反應(yīng)速率最高的位置作為假想的反應(yīng)鋒面，將模擬結(jié)果與推導(dǎo)結(jié)果比較，最終結(jié)果如圖4所示。
圖4表明，熱勢理論修正方法和準(zhǔn)靜態(tài)方法可以較為精確地衡量鋒面前進(jìn)的情況，而原始熱勢理論的結(jié)果大大偏離了實際情況。從參數(shù)分析中發(fā)現(xiàn)，單位長度反應(yīng)內(nèi)熱源的數(shù)值較大，約為單位長度床層顯熱（床層熱容和溫差決定）的20倍，這一點即反映在廣義Stefan數(shù)中。修正方法中如實體現(xiàn)了反應(yīng)內(nèi)熱源的巨大影響，基本可忽略床層顯熱變化對滲透距離前進(jìn)速率的影響；Biot原始熱勢理論缺少了13ρct0′+′2項，內(nèi)熱源的影響沒有充分考慮，預(yù)測耦合過程由床層顯熱主導(dǎo)，其結(jié)果偏離了實際的情況。
3 結(jié)果與討論
3.1 溫度分布假設(shè)對結(jié)果的影響
與1.3.1節(jié)不同，這一節(jié)中反應(yīng)鋒面到換熱邊界之間的床層的溫度分布假設(shè)為拋物線型分布。從文獻(xiàn)中可知，研究人員用熱勢理論求解無內(nèi)熱源的問題時，采用的即為拋物線型溫度分布［16］。在本研究中，可以認(rèn)為內(nèi)熱源強度趨近于0時，溫度分布為拋物線型，描述為
t=t01－xδ2 （35）
不同時間下的溫度分布可以表示為圖5。可以確定，當(dāng)溫差不為0而內(nèi)熱源強度趨近于0時，廣義Stefan數(shù)將趨近于正無窮。嘗試將熱勢理論修正方法和準(zhǔn)靜態(tài)方法的溫度假設(shè)更換為拋物線假設(shè)并重新推導(dǎo)。
3.1.1 拋物線型溫度分布假設(shè)下的熱勢理論修正方法
根據(jù)能量守恒方程可得熱位移場H為
H=13ρct0δ1－xδ3+′δ1－xδ （36）
熱位移場對滲透距離求導(dǎo)為
H·=?H?τ=
13ρct01－xδ3+ρct0xδ1－xδ2+′ （37）
耗散函數(shù)和拉格朗日方程耗散項可求得
D=12∫1kH·2dV=
13630kρ2c2t20+16kρct0′+12k′22δ（38）
?D?=13315kρ2c2t20+13kρct0′+1k′2δ （39）
除此之外，熱勢的定義和它的求解過程如下
E=12∫ρct2dV=12∫δ0ρct201－xδ4dx=110ρct20δ（40）
熱勢對滲透距離的偏導(dǎo)，可得
?E?δ=110ρct20 （41）
廣義熱力亦可求得
Q=∫t?H?δdS=∫t0?H?δdS=13ρct20+′t0（42）
將式（40）～（42）代入拉格朗日方程式（4），可得
δ（τ）=2k（7/30）ρct20+q·′t0（13/315）ρ2c2t20+（1/3）ρct0q·′+q·′2τ（43）
3.1.2 拋物線型溫度分布假設(shè)下的準(zhǔn)靜態(tài)方法
需要強調(diào)的是，假設(shè)內(nèi)熱源強度趨近于0時廣義Stefen數(shù)趨近于正無窮，嚴(yán)格意義上準(zhǔn)靜態(tài)方法并不適用這種情況。為了討論拋物線型溫度分布假設(shè)在問題求解中的重要性，本節(jié)進(jìn)行了這一并不精確的推導(dǎo)，借此分析溫度分布假設(shè)對結(jié)果的影響。事實上，在修正方法中應(yīng)用拋物線型溫度分布同樣存在一定不嚴(yán)謹(jǐn)?shù)膯栴}。因為對于拋物線型溫度分布來說，在反應(yīng)鋒面（滲透距離）處溫度梯度（曲線斜率）為0。基于傅里葉導(dǎo)熱定律進(jìn)行分析會發(fā)現(xiàn)，反應(yīng)鋒面該處的熱流密度亦為0。從這一角度看，反應(yīng)區(qū)域并沒有進(jìn)行任何熱量的傳遞，反應(yīng)鋒面也不應(yīng)該推進(jìn)。研究者討論過這一問題［24］。這是由于，在變分公式中，熱位移場 H 的時間導(dǎo)數(shù)表示熱通量而不是溫度梯度。傅里葉傳導(dǎo)定律在變分原理中并不直接使用，而是由變分原理近似。拋物線型溫度分布假設(shè)是出于近似計算的需要，同時研究者們也研究過指數(shù)型、三次方型曲線等［32］。總的來說，拋物線型溫度分布應(yīng)用最為廣泛、計算方便、精度較高。
與1.5節(jié)推導(dǎo)過程類似，反應(yīng)發(fā)生后滲透距離內(nèi)床層積累顯熱Q+為
Q+=∫ρctdV=∫δ0ρct01－xδ2dx=13ρct0δ（44）
當(dāng)滲透距離前進(jìn)dδ時，床層溫度分布發(fā)生變化，床層中增加的顯熱dQ+為
dQ+=∫δ+dδ0ρct01－xδ2dx－∫δ0ρct01－xδ2dx=
13ρct0dδ（45）
同時，厚度為dδ的床層發(fā)生了反應(yīng)，對應(yīng)需要的熱量為
dQ*=′dδ （46）
由傅里葉導(dǎo)熱定律可知，外界提供的熱流近似為
Q·=kdtdx=k2t0δ=13ρct0+′dδdτ（47）
最終可以求得結(jié)果
δ=2kt0（1/3）ρct0+q·′τ（48）
單純比較不同方法、不同溫度分布假設(shè)下求得的4個滲透距離的表達(dá)式（式（23）、式（31）、式（43）和式（48））是困難的，需要結(jié)合影響傳熱過程的重要參數(shù)來討論其合理性。對傳熱過程影響最明顯的參數(shù)包括材料的反應(yīng)熱和床層顯熱，同時這也是直接決定廣義Stefan數(shù)大小的兩個關(guān)鍵參數(shù)。
3.2 內(nèi)熱源強度對結(jié)果的影響
從前人的研究可知，在內(nèi)熱源強度趨近于0時或者說廣義Stefen數(shù)趨近于正無窮時，求近似解采用拋物線型溫度分布進(jìn)行計算是合理的。將4個結(jié)果的內(nèi)熱源強度取值為0，即將式（23）、式（31）、式（43）和式（48）中的′取值為0，可得線性溫度分布下熱勢理論修正方法及準(zhǔn)靜態(tài)方法的δ表達(dá)式
δ=2.24ατ （49）
δ=2.00ατ （50）
拋物線型溫度分布下，熱勢理論修正方法及準(zhǔn)靜態(tài)方法的δ表達(dá)式
δ=3.36ατ （51）
δ=3.46ατ （52）
式中：α為熱擴散系數(shù)，等價于k/ρc。
從式（49）～式（52）可以看出，在無內(nèi)熱源的假設(shè)下，采用線性溫度分布計算的兩種方法的結(jié)果差距較大，而采用拋物線型溫度分布得到的兩種方法的結(jié)果差距相對較小。從前人研究可知，原始熱勢理論的拋物線型溫度分布的精度較高，而去掉內(nèi)熱源影響后，原始熱勢理論和熱勢理論修正方法在機理上是等價的，因此可以認(rèn)為式（51）的結(jié)果是合理的近似解。無內(nèi)熱源假設(shè)下原本的Stefan問題可以被視為另一簡單問題A，即原有均勻溫度teq的半無限大平板，其邊界溫度突然變?yōu)閠w且保持恒定，求解溫度變化情況。問題A存在精確解［33］，任意時刻、任意位置的溫度分布的表達(dá)式如下
t=teq+t02π∫x2ατ0e-η2dη（53）
將各近似解與精確解比較，可求得近似解存在的誤差。式（49）～式（52）誤差分別為11.3%、15.7%、1.8%和1.5%，可見采用線性溫度分布假設(shè)的兩個解誤差較大，都達(dá)到了10%以上，而且兩者之間的差距也很大。采用正確溫度分布假設(shè)的兩個解誤差均小于2%，是合理的近似解。綜上所述，在內(nèi)熱源趨近于0這一特殊情況下，選取了合適溫度分布假設(shè)的熱勢理論修正方法和準(zhǔn)靜態(tài)方法是合理的。需要注意的是，Wang等［12］得出的判斷床層特征尺寸的熱擴散模型在反應(yīng)第二階段忽略了反應(yīng)焓值對反應(yīng)過程的影響，其結(jié)果（δ=ατ）與式（49）～式（52）在形式上是一致的。求解過程極為簡化，并未對溫度分布假設(shè)等做出討論，與精確解的誤差達(dá)到48.0%，遠(yuǎn)超過本文討論的4種近似解結(jié)果。但是，拋開其精度問題，該模型與問題A的描述是一致的，在金屬氫化物理論模型構(gòu)建上啟發(fā)了新的思路。
3.3 床層顯熱對結(jié)果的影響
從前人研究可知，在床層顯熱為0時或者說廣義Stefen數(shù)趨近于0時，采用線性溫度分布假設(shè)是合理且精確的，應(yīng)用準(zhǔn)靜態(tài)方法得到的解可以視作精確解［21］。在這一節(jié)，將4個結(jié)果中的床層顯熱取值為0，即將式（23）、式（31）、式（43）和式（48）中的床層熱容c取值為0，可得線性溫度分布下熱勢理論修正方法及準(zhǔn)靜態(tài)方法的δ（τ）為
δ（τ）=1.41kt0q·′τ （54）
δ（τ）=1.41kt0q·′τ （55）
拋物線型溫度分布下熱勢理論修正方法及準(zhǔn)靜態(tài)方法的δ（τ）為
δ（τ）=1.41kt0q·′τ （56）
δ（τ）=2kt0q·′τ （57）
前面已經(jīng)提到，應(yīng)用線性溫度分布假設(shè)的準(zhǔn)靜態(tài)方法得到的式（55）是該情況下的精確解，因此4個結(jié)果誤差分別為0、0、0和41.8%。比較發(fā)現(xiàn)，應(yīng)用線性溫度分布假設(shè)得到的兩個解相等，說明溫度分布假設(shè)符合真實情況時得到的結(jié)果合理；應(yīng)用拋物線型溫度分布假設(shè)的熱勢理論修正方法得到的結(jié)果也與精確解相同。這說明熱勢理論修正方法在顯熱等于0或接近于0的情況下，即使選擇了不準(zhǔn)確的溫度分布假設(shè)仍可以得到準(zhǔn)確的解。結(jié)合3.2節(jié)分析結(jié)果可以發(fā)現(xiàn)，熱勢理論修正方法選用兩種溫度分布假設(shè)得到的近似解相差不大，而準(zhǔn)靜態(tài)方法只有選擇準(zhǔn)確合理的溫度分布假設(shè)時才能得到合理的近似解。Marty等［11］提出的反應(yīng)鋒面模型同樣忽略了床層熱容對反應(yīng)過程的影響，其結(jié)果與線性溫度分布的準(zhǔn)靜態(tài)方法得到的解（式（55））是等價的。
金屬氫化物反應(yīng)器中存在的反應(yīng)鋒面推進(jìn)現(xiàn)象更為接近廣義Stefen數(shù)趨近于0的情況，但是床層熱容仍可能對反應(yīng)情況存在影響，討論床層熱容為0假設(shè)下得到的近似結(jié)果與圖4中模擬結(jié)果之間的差距是有必要的。將表1中LaNi5材料的主要參數(shù)代入式（55）求得β為3.011×10－4 m/s1/2。該β求得的反應(yīng)鋒面位置與模擬結(jié)果平均誤差為3.1%，未做該簡化假設(shè)的熱勢理論修正方法近似解平均誤差為0.6%，準(zhǔn)靜態(tài)方法近似解平均誤差為0.6%，這兩種求解方法的精度要高于采用了忽略熱容假設(shè)的近似解結(jié)果。
有必要從物理角度分析內(nèi)熱源強度、床層顯熱和兩者之間的相對大小（廣義Stefen數(shù)）在反應(yīng)器氫熱耦合現(xiàn)象中的作用。從熱量分配的角度看，外界熱源提供的熱量（或冷量）一部分供給床層內(nèi)熱源用于脫氫反應(yīng)（或吸氫反應(yīng)），一部分被床層吸收導(dǎo)致其自身溫度上升（或下降）。當(dāng)外界熱源提供的熱流一定的情況下，單位時間內(nèi)反應(yīng)熱和床層顯熱需要的總熱量恒定。廣義Stefen數(shù)的大小體現(xiàn)了兩者之間熱量分配的情況。當(dāng)內(nèi)熱源強度不變，如果床層熱容越小，那么床層顯熱對傳熱過程的影響越小，外界熱源提供的熱量幾乎全部流至反應(yīng)區(qū)域供給吸放氫反應(yīng)。對于矩形反應(yīng)器來說，此時床層的溫度分布曲線接近平直的線性分布，第一類邊界條件下反應(yīng)區(qū)域的移動將更為快速。這反映了金屬氫化物反應(yīng)特性和熱物性對傳熱過程的影響。
在廣義Stefen數(shù)趨近于0和趨近于正無窮兩種特殊情況下，熱勢理論修正方法在選擇了合理溫度分布假設(shè)時，得到的結(jié)果都是合理的；當(dāng)選擇了不合理的溫度分布假設(shè)時，得到的結(jié)果精度仍比準(zhǔn)靜態(tài)方法要高。綜上所述，熱勢理論修正方法對不同工況條件都具有較好的適應(yīng)性和精確性，是求解金屬氫化物反應(yīng)器反應(yīng)過程理論模型的合適方法。
4 結(jié) 論
（1）熱勢理論修正方法求解的理論模型的結(jié)果與數(shù)值模型結(jié)果相比平均誤差為0.6%，證明該理論模型具有準(zhǔn)確性和合理性。該理論模型可以應(yīng)用到吸放氫速率估計、儲氫系統(tǒng)設(shè)計和未來的導(dǎo)熱優(yōu)化中。原始熱勢理論不能完整地描述內(nèi)熱源作用對滲透距離的影響，偏離了實際物理過程，求解得到的結(jié)果證明Biot將與內(nèi)熱源有關(guān)的熱位移場單獨處理的方法并不合理。
（2）理論模型的推導(dǎo)結(jié)果顯示鋒面位置與時間的平方根、導(dǎo)熱系數(shù)的平方根成正比，表明熱勢理論修正方法所求得的解符合傳統(tǒng)的Stefan問題中解的一般形式。
（3）分別利用熱勢理論修正方法和準(zhǔn)靜態(tài)方法在不同溫度分布假設(shè)下求解理論模型，4個求解結(jié)果記作：熱勢理論修正方法-線性溫度分布，準(zhǔn)靜態(tài)方法-線性溫度分布，熱勢理論修正方法-拋物線型溫度分布，準(zhǔn)靜態(tài)方法-拋物線型溫度分布。廣義Stefen數(shù)趨近于正無窮時，4個結(jié)果誤差分別為11.3%、15.7%、1.8%和1.5%。廣義Stefen數(shù)趨近于0時，4個結(jié)果誤差分別為0、0、0和41.8%。
（4）結(jié)論3說明，廣義Stefen數(shù)趨近于正無窮時，拋物線型溫度分布假設(shè)是合理的；廣義Stefen數(shù)趨近于0時，線性溫度分布假設(shè)是合理的。對于金屬氫化物材料的反應(yīng)器來說，其廣義Stefen數(shù)的值普遍小于1，線性溫度分布假設(shè)是合理的，此時床層顯熱對反應(yīng)鋒面推進(jìn)過程影響有限。以LaNi5材料為例，如果忽略床層顯熱利用熱勢理論修正方法求解理論模型，其結(jié)果與模擬結(jié)果相比誤差為3.1%。
（5）熱勢理論修正方法和準(zhǔn)靜態(tài)方法在選擇合理的溫度分布假設(shè)時才能得到精確的近似解。當(dāng)選擇的溫度分布假設(shè)與實際溫度分布情況偏差較大的情況下，兩種方法得到的結(jié)果誤差較大。此時，熱勢理論修正方法的近似解誤差總體小于準(zhǔn)靜態(tài)方法，因此可以在不確定實際溫度分布情況下進(jìn)行近似求解，但是其計算過程是復(fù)雜的。
參考文獻(xiàn)：
［1］SREEDHAR I， KAMANI K M， KAMANI B M， et al. A bird’s eye view on process and engineering aspects of hydrogen storage ［J］. Renewable and Sustainable Energy Reviews， 2018， 91： 838-860.
［2］AFZAL M， MANE R， SHARMA P. Heat transfer techniques in metal hydride hydrogen storage： a review ［J］. International Journal of Hydrogen Energy， 2017， 42（52）： 30661-30682.
［3］KLEIN H P， GROLL M. Heat transfer characteristics of expanded graphite matrices in metal hydride beds ［J］. International Journal of Hydrogen Energy， 2004， 29（14）： 1503-1511.
［4］BAI Xiaoshuai， YANG Weiwei， ZHANG Wanyu， et al. Hydrogen absorption performance of a novel cylindrical MH reactor with combined loop-type finned tube and cooling jacket heat exchanger ［J］. International Journal of Hydrogen Energy， 2020， 45（52）： 28100-28115.
［5］ASKRI F， BEN SALAH M， JEMNI A， et al. Optimization of hydrogen storage in metal-hydride tanks ［J］. International Journal of Hydrogen Energy， 2009， 34（2）： 897-905.
［6］GUO Zengyuan， ZHU Hongye， LIANG Xingang. Entransy： a physical quantity describing heat transfer ability ［J］. International Journal of Heat and Mass Transfer， 2007， 50（13/14）： 2545-2556.
［7］ZHANG Shiwei， YANG Fusheng， ZHOU Lu， et al. A novel multilayer fin structure for heat transfer enhancement in hydride-based hydrogen storage reactor ［J］. International Journal of Energy Research， 2018， 42（12）： 3837-3850.
［8］CORGNALE C， HARDY B J， TAMBURELLO D A， et al. Acceptability envelope for metal hydride-based hydrogen storage systems ［J］. International Journal of Hydrogen Energy， 2012， 37（3）： 2812-2824.
［9］VISARIA M， MUDAWAR I， POURPOINT T， et al. Study of heat transfer and kinetics parameters influencing the design of heat exchangers for hydrogen storage in high-pressure metal hydrides ［J］. International Journal of Heat and Mass Transfer， 2010， 53（9/10）： 2229-2239.
［10］YANG Fusheng， MENG Xiangyu， DENG Jianqiang， et al. Identifying heat and mass transfer characteristics of metal hydride reactor during adsorption： parameter analysis and numerical study ［J］. International Journal of Hydrogen Energy， 2008， 33（3）： 1014-1022.
［11］MARTY P， DE RANGO P， DELHOMME B， et al. Various tools for optimizing large scale magnesium hydride storage ［J］. Journal of Alloys and Compounds， 2013， 580， Supplement 1： S324-S328.
［12］WANG Chunsheng， BRINKERHOFF J. Is there a general time scale for hydrogen storage with metal hydrides or activated carbon？ ［J］. International Journal of Hydrogen Energy， 2021， 46（21）： 12031-12034.
［13］程新廣， 李志信， 過增元. 熱傳導(dǎo)中的變分原理 ［J］. 工程熱物理學(xué)報， 2004， 25（3）： 457-459.
CHENG Xinguang， LI Zhixin， GUO Zengyuan. Variational principles in heat conduction ［J］. Journal of Engineering Thermophysics， 2004， 25（3）： 457-459.
［14］BIOT M A. Theory of stress-strain relations in anisotropic viscoelasticity and relaxation phenomena ［J］. Journal of Applied Physics， 1954， 25（11）： 1385-1391.
［15］BIOT M A. New methods in heat flow analysis with application to flight structures ［J］. Journal of the Aeronautical Sciences， 1957， 24（12）： 857-873.
［16］BIOT M A. Variational principles in heat transfer ［M］. Oxford， UK： Oxford University Press， 1970.
［17］宋柏， 吳晶， 過增元. 基于熱質(zhì)理論的Hamilton原理 ［J］. 物理學(xué)報， 2010， 59（10）： 7129-7134.
SONG Bai， WU Jing， GUO Zengyuan. Hamilton’s principle based on thermomass theory ［J］. Acta Physica Sinica， 2010， 59（10）： 7129-7134.
［18］程新廣. （火積）及其在傳熱優(yōu)化中的應(yīng)用 ［D］. 北京： 清華大學(xué)， 2004.
［19］夏再忠. 導(dǎo)熱和對流換熱過程的強化與優(yōu)化 ［D］. 北京： 清華大學(xué)， 2001.
［20］程新廣， 夏再忠， 李志信， 等. 導(dǎo)熱優(yōu)化： 熱耗散與最優(yōu)導(dǎo)熱系數(shù)場 ［J］. 工程熱物理學(xué)報， 2002， 23（6）： 715-717.
CHENG Xinguang， XIA Zaizhong， LI Zhixin， et al. Optimization of heat conduction： thermal dissipation and optimal thermal conductivity distribution ［J］. Journal of Engineering Thermophysics， 2002， 23（6）：
715-717.
［21］LOCK G S H. Latent heat transfer： an introduction to fundamentals ［M］. Oxford： Oxford University Press， 1994.
［22］姜禮尚. 二相Stefan問題：Ⅰ ［J］. 數(shù)學(xué)學(xué)報， 1963， 13（4）： 631-646.
［23］姜禮尚. 二相Stefan問題：Ⅱ ［J］. 數(shù)學(xué)學(xué)報， 1964， 14（1）： 33-49.
［24］LARDNER T J. Approximate solutions to phase-change problems ［J］. AIAA Journal， 1967， 5（11）： 2079-2080.
［25］宋柏， 吳晶， 過增元. 基于拉格朗日方程的含源導(dǎo)熱問題分析 ［J］. 工程熱物理學(xué)報， 2011， 32（1）： 141-144.
SONG Bai， WU Jing， GUO Zengyuan. Analysis of heat conduction problem involving heat source using Lagrangian equation ［J］. Journal of Engineering Thermophysics， 2011， 32（1）： 141-144.
［26］WU Jing， GUO Zengyuan， SONG Bai. Application of Lagrange equations in heat conduction ［J］. Tsinghua Science amp; Technology， 2009， 14（S2）： 12-16.
［27］CREPEAU J， SIAHPUSH A. Approximate solutions to the Stefan problem with internal heat generation ［J］. Heat and Mass Transfer， 2008， 44（7）： 787-794.
［28］MCCORD D， CREPEAU J， SIAHPUSH A， et al. Analytical solutions to the Stefan problem with internal heat generation ［J］. Applied Thermal Engineering， 2016， 103： 443-451.
［29］WANG Jing， LIU Jiaxuan， WU Zhen， et al. A novel design for fin profile in metal hydride reactor towards heat transfer enhancement： considering the limitation of the entransy theory in practical application ［J］. International Journal of Heat and Mass Transfer， 2024， 235： 126221.
［30］BAO Zewei， YANG Fusheng， WU Zhen， et al. Optimal design of metal hydride reactors based on CFD-Taguchi combined method ［J］. Energy Conversion and Management， 2013， 65： 322-330.
［31］MACDONALD B D， ROWE A M. Impacts of external heat transfer enhancements on metal hydride storage tanks ［J］. International Journal of Hydrogen Energy， 2006， 31（12）： 1721-1731.
［32］LARDNER T J. Biot’s variational principle in heat conduction ［J］. AIAA Journal， 1963， 1（1）： 196-206.
［33］ECKERT E R G. Analysis of heat and mass transfer ［M］. ［S.l.］： McGraw-Hill， 1972.
（編輯 杜秀杰）
