基于空間基函數(shù)優(yōu)化組合的微波加熱溫度分布數(shù)值求解

楊 彪 馬紅濤杜 婉劉 承高 皓

(1.昆明理工大學(xué)信息工程與自動化學(xué)院,云南昆明 650500;2.昆明理工大學(xué)云南省人工智能重點實驗室,云南昆明 650500;3.昆明理工大學(xué)非常規(guī)冶金教育部重點實驗室,云南昆明 650093)

1 引言

近些年來隨著微波加熱技術(shù)的成熟,微波加熱在工業(yè)方面取得了廣泛的應(yīng)用.與傳統(tǒng)的加熱方式相比微波加熱具有高效率,清潔無污染的優(yōu)點.用于食品加工、木材干燥以及化工產(chǎn)業(yè)等[1-3].微波加熱與傳統(tǒng)熱工過程有很大不同,傳統(tǒng)加熱中熱量從媒質(zhì)外部通過對流、傳導(dǎo)、熱輻射進(jìn)入媒質(zhì)內(nèi)部,而微波加熱則是通過電磁場使得被加熱媒質(zhì)內(nèi)部極性分子在短時間內(nèi)快速相互作用從而產(chǎn)生大量的熱.微波加熱在加熱效率快速高效的同時也伴隨著加熱溫度分布不均勻的缺點,熱點和冷點的出現(xiàn)會導(dǎo)致局部過熱的現(xiàn)象,當(dāng)局部溫度超過材料臨界溫度時就會產(chǎn)生熱失控[4].熱失控現(xiàn)象是阻礙微波加熱工業(yè)應(yīng)用發(fā)展的重大阻礙.研究求解微波加熱過程中媒質(zhì)的全局溫度分布是預(yù)防和解決熱失控現(xiàn)象的首要問題.

微波加熱溫度模型由熱傳遞方程、邊界條件、初始溫度條件組成,熱傳遞方程可由一類典型的拋物線型偏微分方程(partial differential equations,PDE)表示,其具有強(qiáng)烈的時空耦合特性及時間與空間無限維的特點,并且熱傳遞模型的非齊次耗散功率項的存在導(dǎo)致微波加熱溫度模型具有強(qiáng)烈的非線性特征,再加上熱動力學(xué)場非齊次邊界條件和初始溫度的影響,利用解析方法獲得全局溫度分布的表達(dá)式幾乎是不可能的,因此不少學(xué)者對微波加熱過程中的溫度分布求解提出了多種數(shù)值計算方法.Yang等[5]通過COMSOL建立了多饋微波裝置流體加熱的三維有限元數(shù)值模型,模擬了不同功率、初始溫度等參數(shù)對加熱溫度的影響,最終選擇出最佳加熱管半徑.Oliveira等[6]利用有限元法與分析技術(shù)相結(jié)合解決了加載陶瓷樣品的單模填充器內(nèi)部熱電磁問題.Ma[7]等利用時域有限差分法(finite difference time domain method,FDTD)描述微波加熱德拜物質(zhì)的溫度分布.此外還有許多數(shù)值方法,比如傳輸線矩陣法(transmission line matrix,TLM)、矩量法(method of moments,MOM)、有限體積法(finite volume method,FVM)等[8-10].雖然上述數(shù)值方法能有效的解決微波加熱過程中的溫度分布問題,但是在計算過程中計算速度、精度和實時性很難取得良好的平衡,因此需尋求一種新的求解微波加熱溫度分布的方法.

Zhong等[11]利用輔助函數(shù)提出一個中間模型,將PDE溫度模型的非齊次Dirichlet邊界條件轉(zhuǎn)化為齊次Dirichlet邊界條件,便于求得PDE模型空間微分算子的特征值與特征函數(shù),并將特征函數(shù)作為全局空間基函數(shù)[12],以空間基函數(shù)和其對應(yīng)時間系數(shù)的級數(shù)和的形式將PDE模型進(jìn)行時空分離,結(jié)合伽遼金權(quán)重殘差法(Galerkin weighted residual method,GWRM)[13]和譜方法成功將無限維PDE模型轉(zhuǎn)化為近似有限維ODE模型.并將模型的數(shù)值仿真結(jié)果與傳統(tǒng)數(shù)值方法求解結(jié)果進(jìn)行對比驗證了該近似模型的有效性.雖然該方法將無限維模型近似為有限維模型,大大減小了模型計算的復(fù)雜度,但是在該方法中有限維ODE模型階數(shù)的選取有較大的主觀性,并且根據(jù)理論分析判斷模型階數(shù)的選取對全局溫度分布的求解有較大的影響,當(dāng)階數(shù)取較大數(shù)值時,有限維ODE溫度模型維數(shù)較高能保證全局溫度分布求解的精度但會使計算時間過長實時性降低.當(dāng)階數(shù)取較小數(shù)值時,實時性增強(qiáng)但精度會降低.尋找一個合適的方法使得ODE溫度模型溫度分布求解的計算速度和計算精度達(dá)到一個更好的平衡是本文的主要目的.

受到Zhong有限維微波加熱ODE模型的啟發(fā),本文在第3節(jié)中通過尋找一個最優(yōu)的空間基函數(shù)轉(zhuǎn)換矩陣將其選取的全局空間基函數(shù)進(jìn)行線性組合[14],轉(zhuǎn)換矩陣的取值可通過最小化誤差函數(shù)近似得到,誤差函數(shù)表示為優(yōu)化組合后的空間基函數(shù)對應(yīng)的ODE模型與初始空間基函數(shù)對應(yīng)的ODE模型之間預(yù)測輸出的差值.通過優(yōu)化組合后的空間基函數(shù)將有限維ODE溫度模型進(jìn)一步降維.第4節(jié)中,本文首先求解了初始高維ODE模型不同階數(shù)取值下的溫度分布與求解時間,分析了階數(shù)取值的不同對溫度分布求解的影響,并通過不同模型階數(shù)下某些點處溫度的趨勢走向選取了合適的高維ODE模型階數(shù)求得近似標(biāo)準(zhǔn)溫度分布.其次,同樣的將進(jìn)一步降維后的低維ODE模型在不同階數(shù)取值的情況下進(jìn)行數(shù)值仿真求解,得到其溫度分布和計算時間,分別將初始高維ODE模型不同階數(shù)取值下的溫度分布和低維ODE模型在不同階數(shù)取值下溫度分布與近似標(biāo)準(zhǔn)溫度分布進(jìn)行對比.最后,利用有多物理場仿真軟件(COMSOL),即有限元方法(finite element method,FEM)再次驗證了模型的準(zhǔn)確性,充分證明了本文方法的有效性.

2 微波加熱溫度模型

2.1 無限維微波加熱PDE溫度模型

在微波加熱過程中,電磁場隨時間空間變化的規(guī)律可由Maxwell方程表示,麥克斯韋方程組根據(jù)電磁理論,變化的電場和變化的磁場相互激發(fā),相互聯(lián)系形成統(tǒng)一的電磁場.Maxwell 方程組可寫成如下形式[15]:

其中:Je,Jm,ρe,ρm分別表示電流密度、磁流密度、電荷密度和磁荷密度;σm是磁電阻率(Ω/m)對應(yīng)磁場的損耗.在無源媒質(zhì)中ρe=0和ρm=0,在無電或無磁損耗的媒質(zhì)中可令σ或σm等于0.H表示磁場強(qiáng)度,E表示電場強(qiáng)度,B表示磁感應(yīng)強(qiáng)度,D表示電位移矢量.

微波加熱過程中,被加熱媒質(zhì)中的電磁場計算可由Helmholtz方程[16]表示:

其中:k0=表示媒質(zhì)的相對磁導(dǎo)率,ε(T)表示媒質(zhì)的復(fù)介電常數(shù),σ表示媒質(zhì)的電導(dǎo)率,E和ω分別入射電場強(qiáng)度和入射電磁波角頻率,k0表示自由空間的波數(shù),ε0=8.854×10?14F/cm,μ0=4π×10?7H/m表示真空磁導(dǎo)率,j表示虛數(shù)單位.

在笛卡爾坐標(biāo)系下,根據(jù)熱動力學(xué)第一定律及傅里葉熱傳導(dǎo)定律,在一維方向上微波加熱媒質(zhì)過程可由下述含非齊次項的偏微分方程進(jìn)行描述[17]:

非齊次Dirichlet邊界條件:

初始溫度條件:

其中:ρ,Cp,k分別表示加熱媒質(zhì)的密度,比熱容和導(dǎo)熱系數(shù).hc表示熱傳系數(shù),Ta,Tb表示環(huán)境溫度.Pabs表示耗散功率,微波功率的計算可采用朗伯(Lambert)定律,朗伯定律[18]以指數(shù)衰減的形式表示被加熱媒質(zhì)內(nèi)部微波能量吸收的情形,在對被加熱媒質(zhì)的全局溫度分布趨勢的研究上在很大程度上降低了電磁場計算的復(fù)雜性.Lambert定理所適用的條件為:電磁場垂直入射到媒質(zhì)表面,介電常數(shù)不隨溫度變化,并且媒質(zhì)長度是無限的.然而上述使用條件僅僅是在理想狀態(tài)下,而在實際應(yīng)用中,若媒質(zhì)的長度大于或者等于趨膚深度的2倍,就可以忽略透射波和反射波的影響,近似地將耗散功率寫成Lambert定理的形式可表示為

其中:f表示微波頻率,ε′′表示相對介電損耗,ε′表示被加熱媒質(zhì)的相對介電系數(shù).E+表示入射電場強(qiáng)度,Pint表示入射功率密度,Dp為趨膚深度,表示入射功率從媒質(zhì)表面衰減至表面功率值的1/e時的距離,c=3×108m/s為微波在真空環(huán)境中的傳播速度.

2.2 有限維微波加熱ODE溫度模型

由于上述PDE模型具有非齊次邊界條件的制約,難以獲取空間微分算子的特征譜和特征函數(shù),Zhong利用構(gòu)造輔助函數(shù)將模型(3)-(5)等價為如下中間模型[11]:

其中:λi表示特征值,φi(z)表示相對應(yīng)的特征函數(shù).選擇特征函數(shù)φi(z)作為空間基函數(shù),式(9)中的時空耦合變量Θ(z,t)可以能展開成一個無限維空間基函數(shù)的集合和其對應(yīng)的時間系數(shù)的級數(shù)和的形式:

可通過下式進(jìn)行最小化:

結(jié)合式(14),將式(15)帶入式(16)可得到時空變量分離的無限維微波加熱ODE溫度模型:

3 空間基函數(shù)變換

空間基函數(shù)的選擇對于偏微分方程系統(tǒng)的降維十分關(guān)鍵,能極大程度的影響模型的精度.第2節(jié)中通過選取空間微分算子的特征函數(shù)作為降維的全局基函數(shù),其特征函數(shù)是一組無限的光滑的全局正交基函數(shù),非線性偏微分方程動態(tài)系統(tǒng)的主要信息存在于有限維特征譜對應(yīng)的有限維空間基函數(shù)中.本小節(jié)通過將對上一節(jié)中有限維模型的n+1階空間基函數(shù)進(jìn)行線性組合,得到新的m+1(m ＜n)階空間基函數(shù),對有限維微波加熱ODE溫度模型進(jìn)一步降維,得到階數(shù)更低的有限維模型.

3.1 基函數(shù)轉(zhuǎn)換矩陣

定義一個基函數(shù)轉(zhuǎn)換矩陣R(n+1)×m+1),矩陣中每個位置的數(shù)值都是n+1階空間基函數(shù)進(jìn)行線性組合,得到新的m+1(m ＜n)階空間基函數(shù)的線性組合系數(shù),可由下式表達(dá):

基于m+1階新空間基函數(shù)可得到新的有限維微波加熱ODE溫度模型,如下:

3.2 基函數(shù)轉(zhuǎn)換矩陣的求解

為求解空間基函數(shù)的轉(zhuǎn)換矩陣,現(xiàn)給出如下誤差函數(shù):

此誤差函數(shù)表示基于n+1階空間基函數(shù)的ODE溫度模型與基于m+1階空間基函數(shù)的ODE溫度模型的溫度分布的時空積分的平方差,其中tmax表示最大加熱時間,z0表示被加熱媒質(zhì)深度,In+1表示n+1階單位矩陣,從上式可以看出,誤差函數(shù)只與初始n+1階空間基函數(shù)對應(yīng)的時間系數(shù)相關(guān).

空間基函數(shù)轉(zhuǎn)換矩陣可通過最小化誤差函數(shù)e(t)近似得到,本文對e(t)的最小化采用隨機(jī)優(yōu)化的方法來計算求解.粒子群算法[20](particle swarm optimization,PSO)是受自然界中鳥群覓食的啟發(fā)的一種群體智能優(yōu)化算法,并且是高精度的隨機(jī)優(yōu)化算法.但是傳統(tǒng)PSO算法在尋優(yōu)過程中存在陷入局部最優(yōu)和早熟收斂的問題.為更準(zhǔn)確快速的求得基函數(shù)轉(zhuǎn)換矩陣R,本文采用自適應(yīng)變異粒子群優(yōu)化算法(adaptive mutation particle swarm optimization,AMPSO)[21],即為粒子群優(yōu)化算法增加了一個自適應(yīng)變異機(jī)制,即每次算法計算周期,都有一定的概率隨機(jī)改變某一個個體的數(shù)值,從而提高整個粒子群的尋優(yōu)能力.

在粒子群優(yōu)化算法的思想中,每個粒子有兩個屬性,一個是粒子所在的位置,一個是粒子的速度,粒子的運動就是尋找最優(yōu)解的運動.采用粒子群算法進(jìn)行本次優(yōu)化時,每個粒子的位置有(n+1)(m+1)個,其中每個粒子的位置表示為

t表示當(dāng)前時刻,t ?1表示前一時刻,t+1表示后一時刻,ω0為權(quán)重系數(shù),c1,c2為學(xué)習(xí)因子,r1,r2為[0,1]的隨機(jī)數(shù),T表示粒子每次運動時間.粒子根據(jù)自身先前速度以及自身位置與個體最佳位置和群體最優(yōu)位置的距離更新自身速度,之后運動到新的位置.

基于上述PSO算法,優(yōu)化求解空間基函數(shù)轉(zhuǎn)換矩陣R(n+1)×(m+1)的算法流程如下:

4 仿真分析及驗證

4.1 數(shù)值仿真及分析

本小節(jié)首先通過MATLAB數(shù)值仿真驗證分析有限維ODE溫度模型階數(shù)的取值對微波加熱全局溫度分布求解的影響,并選取合適的階數(shù)的模型將其溫度分布作為近似標(biāo)準(zhǔn)溫度分布.接著將ODE溫度模型(23)-(25)與模型(19)-(21)在計算速度和溫度分布誤差方面進(jìn)行比較,證明了基于m+1階空間基函數(shù)對應(yīng)的模型的有效性.最后利用FEM驗證了m+1階空間基函數(shù)對應(yīng)的模型的準(zhǔn)確性.因為去離子水具有熱力學(xué)參數(shù)隨溫度變化很小的特性,為更方便的表示本文模型的有效性,假設(shè)離子水的熱力學(xué)參數(shù)均不隨溫度發(fā)生變化,介電系數(shù)視作恒定,其相應(yīng)熱力學(xué)參數(shù)如表1所示.假設(shè)被加熱去離子水的長度大于其2倍的趨膚深度并充滿長波導(dǎo),其加熱的示意圖如圖1所示.

表1 去離子水恒定熱力學(xué)參數(shù)Table 1 Constant thermodynamic parameters of deionized water

圖1 長波導(dǎo)中加熱去離子水示意圖Fig.1 Schematic diagram of heated deionized water column in long waveguide

當(dāng)微波頻率為2.45 GHz,非齊次邊界條件f1=1°C/cm,f2=?1°C/cm,初始溫度為20°C,入射功率密度為20 W/cm2,加熱時間tmax=10 s,加熱物料深度z0=10 cm時,分別求取n=50,n=47,n=45,n=43,n=23,n=3時有限維ODE溫度模型的溫度分布,其結(jié)果如圖2(a)-(f)所示.通過對比可以看出模型階數(shù)選取對微波加熱全局溫度分布求解有很大的影響,選取階數(shù)越大,計算時間增加溫度分布越精確,相反,計算時間減少溫度分布的誤差越明顯.

圖2(a) n=50溫度分布Fig.2(a) n=50 temperature distribution

圖2(b) n=47溫度分布Fig.2(b) n=47 temperature distribution

圖2(c) n=45溫度分布Fig.2(c) n=45 temperature distribution

圖2(d) n=43溫度分布Fig.2(d) n=43 temperature distribution

圖2(e) n=23溫度分布Fig.2(e) n=23 temperature distribution

圖2(f) n=3溫度分布Fig.2(f) n=3 temperature distribution

從圖2(a)中可以明顯看到在微波加熱過程中饋口位置吸收大量微波能量,迅速升溫,隨著介質(zhì)深度的增加微波功率密度呈指數(shù)型衰減,升溫速度也隨之改變,其溫度變化符合Lambert定律的溫度變化走向,并與Zhong進(jìn)行對比保證了基于初始n+1階空間基函數(shù)對應(yīng)的ODE溫度模型的有效性.取不同模型階數(shù)下t=4 s,z=1.5 cm;t=5 s,z=2 cm;t=6 s,z=1 cm 三點處的溫度趨勢如圖3(a)-(c)所示.

圖3(a) 不同階數(shù)下t=4 s,z=1.5 cm處溫度趨勢Fig.3(a) Temperature trend at t=4 s,z=1.5 cm under different orders

圖3(b) 不同階數(shù)下t=5 s,z=2 cm處溫度趨勢Fig.3(b) Temperature trend at t=5 s, z=2 cm under different orders

圖3(c) 不同階數(shù)下t=6 s,z=1 cm處溫度趨勢Fig.3(c) Temperature trend at t=6 s, z=1 cm under different orders

結(jié)合圖2(a)-(c)可以看出當(dāng)溫度模型n=47時,溫度逐漸趨于一致,只發(fā)生微小變化,其溫度分布滿足初始溫度條件,因此將n=47時模型的溫度分布作為近似標(biāo)準(zhǔn)溫度分布并得到n=47時模型溫度分布的計算時間為0.86 s.

在上述被加熱媒質(zhì)熱力學(xué)參數(shù)、初始溫度條件、非齊次邊界條件、入射功率密度不變的情況下,分別求解n=3,n=8,n=13,n=18,n=23時有限維模型的溫度分布和計算時間.在初始高維有限維ODE模型n=47的基礎(chǔ)上通過空間基函數(shù)轉(zhuǎn)換矩陣將有限維ODE溫度模型進(jìn)一步降維,并分別求解m=3,m=8,m=13,m=18,m=23時模型的溫度分布與求解時間.m,n不同取值情況下模型溫度分布與標(biāo)準(zhǔn)近似溫度分布最大相對誤差如圖4(a)所示,m,n不同取值情況下模型計算時間如4(b)所示.

圖4(a) m,n不同取值下與近似標(biāo)準(zhǔn)溫度分布誤差Fig.4(a) Error of approximate standard temperature distribution with different values of m and n

圖4(b) m,n不同取值下的計算時間Fig.4(b) Calculation time under different values of m and n

m=3時模型的溫度分布如圖5(a)所示,并分別繪制n=3,m=3時模型的溫度分布與近似標(biāo)準(zhǔn)溫度分布的誤差曲面圖,如圖5(b)-(c)所示.

圖5(a) m=3溫度分布Fig.5(a) m=3 temperature distribution

圖5(b) n=3與近似標(biāo)準(zhǔn)溫度分布誤差Fig.5(b) Error of approximate standard temperature distribution with n=3

圖5(c) m=3與近似標(biāo)準(zhǔn)溫度分布誤差Fig.5(c) Error of approximate standard temperature distribution with m=3

結(jié)合圖4(a)-(b),5(b),5(c)可以看出m=3時最大誤差為1.64%.而不經(jīng)空間基函數(shù)優(yōu)化組合直接取n=3時模型的溫度分布與標(biāo)準(zhǔn)近似溫度分布之間的誤差較大,最大誤差為?13.73%.并且m=3時模型計算時間和n=47時相比提升了72.2%.

4.2 模型驗證

為了進(jìn)一步驗證所提出ODE溫度模型(23)-(25)的準(zhǔn)確性,在加熱條件與有限維ODE溫度模型完全相同的情況下,使用FEM對微波加熱PDE溫度模型(3)-(5)進(jìn)行溫度分布求解.

4.2.1 模型參數(shù)與邊界條件

模型由微波反應(yīng)腔體、去離子水及微波入射端三部分組成.反應(yīng)腔體與微波入射端為長方體.為便于計算過程中的網(wǎng)格劃分,去離子水并未充滿整個反應(yīng)腔體,與腔體壁之間留有極小的距離.去離子水與腔壁之間的介質(zhì)為空氣.包裹被加熱去離子水的材料為聚四氟乙烯(PTFE),該材料介電損耗角很小,為透波材料,進(jìn)行數(shù)值計算時可忽略其對微波加熱過程的影響.腔體與WR340型標(biāo)準(zhǔn)波導(dǎo)使用銅材質(zhì),腔體規(guī)格為100 mm×86.4 mm×43.2 mm,饋送微波頻率為2.45 GHz.系統(tǒng)從初始溫度20°C開始加熱,加熱時間為10 s.銅與空氣的熱力學(xué)參數(shù)[22]為COMSOL內(nèi)置參數(shù),被加熱去離子水熱力學(xué)參數(shù)如上述表1所示,微波加熱去離子水幾何模型如圖6(a)所示.

圖6(a) 微波加熱去離子水幾何模型Fig.6(a) Geometric model of deionized water heated by microwave

邊界條件分為電磁邊界條件和熱傳導(dǎo)邊界條件[22].針對電磁邊界條件,反應(yīng)腔壁和波導(dǎo)壁材質(zhì)為銅,可視為完美導(dǎo)體,可由下列公式表達(dá):

其中腔壁表面的單位法向量.

熱傳導(dǎo)邊界條件可由下式表達(dá):

其中Tair表示反應(yīng)腔體內(nèi)空氣介質(zhì)的溫度,h表示空氣的熱傳系數(shù)(10 W/(m2·K)).同時,假設(shè)去離子水與腔體底部的接觸邊界熱絕緣,即

4.2.2 數(shù)值計算結(jié)果與分析

利用第4.2.1小節(jié)中所提數(shù)值計算模型加熱去離子水10 s,其整體表面溫度分布與多水平切面溫度分布如圖6(b)所示.

圖6(b) 整體表面溫度分布(上)與多水平切面溫度分布(下)Fig.6(b) Overall surface temperature distribution(left)and multi-level section temperature distribution(right)

選取x=?15 mm時的水平切面,繪制z=0.5 cm,z=2 cm,z=4.5 cm,三點處10 s內(nèi)的溫升曲線并與m=3時ODE溫度模型的溫升曲線進(jìn)行對比,如圖7(a)-7(c)所示.

通過溫升曲線對比圖7(a)-7(c)可知,經(jīng)空間基函數(shù)轉(zhuǎn)換后的低維ODE溫度模型與有限元模型相比也能精確的描述微波加熱媒質(zhì)的溫升過程.因此,從數(shù)值分析的角度可驗證本文所提出的微波加熱ODE溫度模型是正確且有效的.

圖7(a) z=0.5 cm時溫升曲線對比Fig.7(a) Comparison of temperature rise curves at z=0.5 cm

圖7(c) z=4.5 cm時溫升曲線對比Fig.7(c) Comparison of temperature rise curves at z=4.5 cm

5 結(jié)論

本文在有限維ODE模型的基礎(chǔ)上,分析了有限維ODE溫度模型階數(shù)的選取對被加熱媒質(zhì)溫度分布求解的影響,選取n=47時的模型的溫度分布作為近似標(biāo)準(zhǔn)溫度分布,并且利用空間基函數(shù)優(yōu)化組合的方法對有限維ODE溫度模型進(jìn)一步降維,得到新的低維微波加熱ODE溫度模型.經(jīng)仿真實驗證明,優(yōu)化組合后的m+1階空間基函數(shù)對應(yīng)的溫度模型求解得到的溫度分布與標(biāo)準(zhǔn)近似溫度分布之間的誤差維持在很小的范圍內(nèi),m=3時最大誤差為1.64%.而不經(jīng)空間基函數(shù)優(yōu)化組合直接取n=3時模型的溫度分布與標(biāo)準(zhǔn)近似溫度分布之間的誤差較大,最大誤差為?13.73%.并且m=3時模型計算時間和n=47時相比提升了72.2%.并且通過多物理場耦合FEM 求解微波加熱PDE溫度模型,進(jìn)一步驗證了優(yōu)化后的低維ODE溫度模型的準(zhǔn)確性.本方法能在降低少許精度的情況下,大幅度提升溫度分布求解速度,極大提高微波加熱溫度采集的實時性,具有較高的工業(yè)應(yīng)用價值.