線上求解法最優(yōu)近似格式的探討及應(yīng)用

劉 陽(yáng), 李 凌

(沈陽(yáng)化工大學(xué) 信息工程學(xué)院, 遼寧 沈陽(yáng) 110142)

現(xiàn)代科學(xué)研究過(guò)程中,研究人員常常將研究對(duì)象抽象成數(shù)學(xué)模型,進(jìn)而通過(guò)數(shù)學(xué)建模預(yù)測(cè)研究對(duì)象的變化趨勢(shì).在建立模型和仿真的過(guò)程中常常需要建立數(shù)學(xué)方程并進(jìn)行求解,但在比較復(fù)雜的系統(tǒng)中,所建立模型的數(shù)學(xué)方程經(jīng)常以偏微分方程的形式出現(xiàn).求解偏微分方程時(shí)常常難以得到精確的解析解,所以想對(duì)此類(lèi)系統(tǒng)進(jìn)行全面的數(shù)值分析就需要對(duì)偏微分方程進(jìn)行離散化,化為常見(jiàn)易解的常微分方程形式再求解.線上求解法(method of line,MOL)是一種常用的求解偏微分方程的方法[1-2],該方法是將一個(gè)自變量當(dāng)成連續(xù)變量,再用有限差分法或正交配置法離散剩余的自變量[3],從而將偏微分方程轉(zhuǎn)變?yōu)槌Ｎ⒎址匠探M,將問(wèn)題轉(zhuǎn)化為求解常微分方程組.

在化工過(guò)程動(dòng)態(tài)分析[4]中,基于MATLAB的線上求解法思想有著廣泛應(yīng)用.危鳳等[5]將線上求解法與MATLAB工具相結(jié)合,解決了色譜分離動(dòng)力學(xué)模型的求解與仿真,快速準(zhǔn)確地模擬了色譜分離的過(guò)程;危鳳、沈波等[6]還將基于MATLAB的線上求解法應(yīng)用于奧美拉唑?qū)τ丑w的色譜模擬分離上,對(duì)模擬移動(dòng)床色譜分離技術(shù)制備手性藥物化工過(guò)程的建模仿真做出了較好優(yōu)化;在乙烯環(huán)氧化反應(yīng)工業(yè)領(lǐng)域,石向成等[7]運(yùn)用線上求解法思想,建立了該反應(yīng)的固定床列管反應(yīng)器數(shù)學(xué)模型,對(duì)乙烯環(huán)氧化反應(yīng)的工藝條件做出了分析;鄧家璧等[8]在模擬移動(dòng)床分離過(guò)程的模型預(yù)測(cè)控制方法研究里,結(jié)合線上求解法思想與MATLAB系統(tǒng)辨識(shí)工具箱,設(shè)計(jì)了一個(gè)預(yù)測(cè)控制器對(duì)模擬移動(dòng)床進(jìn)行控制,取得了較好效果.

本文在介紹線上求解法原理的基礎(chǔ)上,通過(guò)使用MATLAB軟件求解工具箱,對(duì)一維Burgers方程進(jìn)行仿真實(shí)驗(yàn),結(jié)合仿真圖像曲線分析線上求解法的最優(yōu)近似格式,并應(yīng)用于苯加氫管式固定床反應(yīng)器的化工反應(yīng)工程仿真建模,驗(yàn)證了該方法的可靠性.

1 線上求解法

線上求解法的求解大致分為兩個(gè)步驟：(1)對(duì)方程中的偏微分項(xiàng)進(jìn)行離散; (2)在離散后的方程上對(duì)時(shí)間變量進(jìn)行積分求解[1].

1.1 偏微分項(xiàng)的離散

偏微分方程

xt=f(t,z,x,xz,xzz),z∈D,t>0.

(1)

邊界條件

a(t,z,x,xz)=0,z∈D,t>0.

(2)

初始條件

x(t=0,z)=x0(z),z∈D∪S,t=0.

(3)

其中：t為時(shí)間且t≥0;z為空間域D中的位置變量;S為空間域D的邊界域;xt為對(duì)時(shí)間變量的一階偏微分項(xiàng);xz為對(duì)空間位置變量的一階偏微分項(xiàng);xzz為對(duì)空間位置變量的二階偏微分項(xiàng).

對(duì)一階偏微分項(xiàng)xz先選取空間網(wǎng)格點(diǎn)zi(i=1,2,…,n),再用泰勒級(jí)數(shù)進(jìn)行展開(kāi),則一階偏微分項(xiàng)的三點(diǎn)中心近似格式為

Δz=zi-zi-1.

(4)

一階偏微分項(xiàng)的五點(diǎn)偏心逆風(fēng)近似格式為

10xzi+3xzi+1)/12Δz+ο(Δz4).

(5)

二階偏微分項(xiàng)的五點(diǎn)中心近似格式為

32xzi+1-2xzi+2)/12Δz2+ο(Δz4).

(6)

在應(yīng)用MATLAB軟件求解工具箱時(shí)可以將式(4)～式(6)用矩陣形式表示.

1.2 時(shí)間變量的積分求解

通過(guò)采用幾種格式近似一階偏微分項(xiàng)和二階偏微分項(xiàng)后,可以通過(guò)MATLAB軟件求解工具箱中的ODE求解器進(jìn)行求解[3],常見(jiàn)調(diào)用格式有ode45、ode23、ode15s等.

2 基于MATLAB的仿真實(shí)驗(yàn)

MATLAB現(xiàn)已在工業(yè)和學(xué)術(shù)界廣泛使用,僅需要較少的編程專(zhuān)業(yè)知識(shí)便可獲得可靠的仿真實(shí)驗(yàn)結(jié)果[9],這為MOL工具箱的開(kāi)發(fā)提供了非常方便的基礎(chǔ).針對(duì)線上求解法在求解偏微分方程時(shí)對(duì)一階偏微分項(xiàng)和二階偏微分項(xiàng)的離散效果,筆者選取了一維Burgers方程作為模型,分別采用幾種不同的離散近似格式,在MATLAB軟件下編寫(xiě)程序調(diào)用ODE/DAE求解器,輸出仿真圖像,結(jié)合圖像進(jìn)行分析探尋并驗(yàn)證離散效果最佳的格式.

2.1 一維Burgers方程

一維Burgers方程,即為簡(jiǎn)化的動(dòng)量平衡方程,如式(7)所示,該方程同時(shí)包含了對(duì)流項(xiàng)與擴(kuò)散項(xiàng).

(7)

式中：x為速度;μ為介質(zhì)黏度.雖然方程是非線性的,但是存在已知的解析解,可用于評(píng)估近似后的數(shù)值解,其解析解為

x(t,z)={0.1exp[-0.05(z-0.5+

4.95t)/μ]+0.5exp[-0.25(z-0.5+

0.75t)/μ]+exp[-0.5(z-

0.375)/μ]}/{exp[-0.05(z-0.5+

4.95t)/μ]+exp[-0.25(z-0.5+

0.75t)/μ]+exp[-0.5(z-0.375)/μ]}.

(8)

2.2 仿真分析

在MATLAB軟件工具下編寫(xiě)程序進(jìn)行仿真實(shí)驗(yàn),再根據(jù)所得圖像進(jìn)行分析.其中μ=0.001;t=0,0.1,0.2,…,1,單位為min;網(wǎng)格數(shù)n=201;解析解輸出圖用虛線顯示,線上求解法近似格式用實(shí)線顯示.先對(duì)式(7)中的二階偏微分項(xiàng)采用五點(diǎn)中心格式進(jìn)行離散近似,再對(duì)一階偏微分項(xiàng)分別進(jìn)行五點(diǎn)偏心逆風(fēng)格式近似、兩點(diǎn)逆風(fēng)格式近似、三點(diǎn)逆風(fēng)格式近似和七點(diǎn)中心格式近似,所得圖像分別如圖1～圖4所示.

圖1中實(shí)線與虛線的貼合度十分完美,即數(shù)值解與解析解非常接近,表明五點(diǎn)偏心逆風(fēng)格式對(duì)一階偏微分項(xiàng)有非常好的離散效果.圖2中顯示一階偏微分項(xiàng)由兩點(diǎn)逆風(fēng)格式近似后陡區(qū)域附近表現(xiàn)出波動(dòng)現(xiàn)象(即解不穩(wěn)定),不適合求解對(duì)流為主的問(wèn)題.采用三點(diǎn)逆風(fēng)格式近似一階偏微分項(xiàng)的圖像如圖3所示.與解析解相比較,三點(diǎn)逆風(fēng)格式近似對(duì)非物理振蕩有比較好的抑制作用,但是陡區(qū)域不能完全貼合解析解,有一定誤差,且在陡區(qū)的中部和后區(qū)的底部有一定振蕩.最后觀察圖4,七點(diǎn)中心格式近似后的圖像前中后部的陡區(qū)和底線與解析解都貼合較好,但是陡區(qū)的頂部都存在振蕩,且在中后部的振蕩幅度和頻率更為明顯.

圖1 五點(diǎn)偏心逆風(fēng)格式近似一階偏微分項(xiàng)與解析解對(duì)比Fig.1 Contrast of five-point biased upwind approximate first-order partial differential terms and analytical solu tions

圖2 兩點(diǎn)逆風(fēng)格式近似一階偏微分項(xiàng)與解析解對(duì)比Fig.2 Contrast of two-point upwind approximate first-order partial differential terms and analytical solutions

圖3 三點(diǎn)逆風(fēng)格式近似一階偏微分項(xiàng)與解析解對(duì)比Fig.3 Contrast of three-point upwind approximate first-order partial differential terms and analytical solutions

圖4 七點(diǎn)中心格式近似一階偏微分項(xiàng)與解析解對(duì)比Fig.4 Contrast of seven-point centered approximate first-order partial differential terms and analytical solutions

對(duì)式(7)中的一階偏微分項(xiàng)給定五點(diǎn)偏心逆風(fēng)格式進(jìn)行離散近似,再對(duì)二階偏微分項(xiàng)分別進(jìn)行五點(diǎn)中心格式近似、兩點(diǎn)逆風(fēng)格式近似、三點(diǎn)逆風(fēng)格式近似,其中五點(diǎn)中心格式近似后圖像仍為圖1,其余兩種格式所得圖像分別如圖5、圖6所示.

圖5 兩點(diǎn)逆風(fēng)格式近似二階偏微分項(xiàng)與解析解對(duì)比Fig.5 Contrast of two-point upwind approximate second-orderpartial differential terms and analytical solutions

圖6 三點(diǎn)逆風(fēng)格式近似二階偏微分項(xiàng)與解析解對(duì)比Fig.6 Contrast of three-point upwind approximate second-order partial differential terms and analytical solutions

結(jié)合圖像可以看出：將兩點(diǎn)逆風(fēng)格式用于近似二階微分項(xiàng)所得的方程曲線與解析解相差很大,僅有陡區(qū)部分貼合,圖像中部產(chǎn)生了巨大的振蕩;三點(diǎn)逆風(fēng)格式方法近似二階偏微分項(xiàng)后的圖像產(chǎn)生了非常大的偏差,從圖6中已看不到解析解的圖像.故可推測(cè)說(shuō)明帶有逆風(fēng)格式的近似方法不適用于近似二階微分項(xiàng).

綜上,在給定條件下通過(guò)對(duì)比幾種近似一維Burgers方程偏微分項(xiàng)的方法,其中采用五點(diǎn)偏心逆風(fēng)格式近似一階偏微分項(xiàng)、五點(diǎn)中心格式近似二階偏微分項(xiàng)后的曲線與解析解的曲線總體貼合度最好,誤差最小.進(jìn)而可以將其推廣到其他領(lǐng)域,在對(duì)方程的偏微分項(xiàng)進(jìn)行近似的時(shí)候,采取此種格式能取得更理想的效果.

3 MOL在化工過(guò)程中的應(yīng)用

化工過(guò)程常常用偏微分方程組形式來(lái)表示,方程除了時(shí)間變量還存在空間自變量,此時(shí)若要對(duì)反應(yīng)過(guò)程進(jìn)行動(dòng)態(tài)分析,方程組求解則要應(yīng)用到MOL思想.其中苯加氫的三區(qū)管式固定床反應(yīng)器的動(dòng)力學(xué)方程組則是一個(gè)典型的偏微分方程組,該模型包括苯和噻吩的物料平衡方程、能量平衡方程以及考慮催化劑失活的方程,形式如式(9)～式(12)所示.

(9)

(10)

(11)

(12)

式中：DT為噻吩擴(kuò)散系數(shù);cT為噻吩的濃度;DB為苯擴(kuò)散系數(shù);cB為苯濃度;λeff為床導(dǎo)熱性系數(shù);R為反應(yīng)器內(nèi)半徑;v為氣體速度;ΔH為焓變;ρc為催化劑系數(shù);ε為床空隙分?jǐn)?shù);ρG為氣體密度;α為壁傳熱系數(shù);rT為噻吩吸附率;rB為加氫率;T為溫度.

根據(jù)以上結(jié)論驗(yàn)證該動(dòng)力學(xué)模型,以五點(diǎn)偏心逆風(fēng)格式近似方程一階偏微分項(xiàng),五點(diǎn)中心格式近似方程二階偏微分項(xiàng),所得反應(yīng)器內(nèi)部的溫度分布曲線如圖7所示.將圖7作為標(biāo)準(zhǔn)圖像.

圖7 標(biāo)準(zhǔn)格式近似圖Fig.7 Standard format approximation

先保持用五點(diǎn)中心格式近似二階微分項(xiàng),對(duì)一階微分項(xiàng)采用七點(diǎn)中心格式近似,所得反應(yīng)器內(nèi)部的溫度分布圖像如圖8所示.與標(biāo)準(zhǔn)圖像對(duì)比可以看出圖8與標(biāo)準(zhǔn)圖像有很大偏差.由于巨大振蕩導(dǎo)致圖像中后部偏離標(biāo)準(zhǔn)圖像特別大,產(chǎn)生了誤差.再將一階偏微分項(xiàng)用五點(diǎn)偏心逆風(fēng)格式近似,對(duì)二階微分項(xiàng)則采用三點(diǎn)逆風(fēng)格式近似,所得反應(yīng)器內(nèi)部的溫度分布圖像如圖9所示.

圖8 七點(diǎn)中心格式近似一階偏微分項(xiàng)Fig.8 Seven-point centered approximate first-order partial differential term

圖9 三點(diǎn)逆風(fēng)格式近似二階偏微分項(xiàng)Fig.9 Three-point upwind approximate second-order partial differential term

與標(biāo)準(zhǔn)圖像對(duì)比分析可知偏心逆風(fēng)格式近似二階偏微分項(xiàng)后圖像曲線單薄,且上部陡區(qū)振蕩明顯,中部下凹曲線與標(biāo)準(zhǔn)圖像相比誤差也較大.

4 結(jié) 論

通過(guò)運(yùn)用線上求解法思想,在給定條件下對(duì)比幾種近似一維Burgers方程偏微分項(xiàng)的方法,得知以五點(diǎn)偏心逆風(fēng)格式近似一階偏微分項(xiàng)、五點(diǎn)中心格式近似二階偏微分項(xiàng)所得曲線具有最理想的效果,并采用苯加氫管式固定床反應(yīng)器這一典型化工過(guò)程進(jìn)行仿真研究,驗(yàn)證了此推論的可靠性.