劉 陽, 李 凌

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

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

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

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

1 線上求解法

線上求解法的求解大致分為兩個步驟：(1)對方程中的偏微分項進行離散; (2)在離散后的方程上對時間變量進行積分求解[1].

1.1 偏微分項的離散

偏微分方程

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為時間且t≥0;z為空間域D中的位置變量;S為空間域D的邊界域;xt為對時間變量的一階偏微分項;xz為對空間位置變量的一階偏微分項;xzz為對空間位置變量的二階偏微分項.

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

Δz=zi-zi-1.

(4)

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

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

(5)

二階偏微分項的五點中心近似格式為

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

(6)

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

1.2 時間變量的積分求解

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

2 基于MATLAB的仿真實驗

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

2.1 一維Burgers方程

一維Burgers方程,即為簡化的動量平衡方程,如式(7)所示,該方程同時包含了對流項與擴散項.

(7)

式中：x為速度;μ為介質(zhì)黏度.雖然方程是非線性的,但是存在已知的解析解,可用于評估近似后的數(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軟件工具下編寫程序進行仿真實驗,再根據(jù)所得圖像進行分析.其中μ=0.001;t=0,0.1,0.2,…,1,單位為min;網(wǎng)格數(shù)n=201;解析解輸出圖用虛線顯示,線上求解法近似格式用實線顯示.先對式(7)中的二階偏微分項采用五點中心格式進行離散近似,再對一階偏微分項分別進行五點偏心逆風(fēng)格式近似、兩點逆風(fēng)格式近似、三點逆風(fēng)格式近似和七點中心格式近似,所得圖像分別如圖1～圖4所示.

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

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

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

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

圖4 七點中心格式近似一階偏微分項與解析解對比Fig.4 Contrast of seven-point centered approximate first-order partial differential terms and analytical solutions

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

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

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

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

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

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

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

(9)

(10)

(11)

(12)

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

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

圖7 標準格式近似圖Fig.7 Standard format approximation

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

圖8 七點中心格式近似一階偏微分項Fig.8 Seven-point centered approximate first-order partial differential term

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

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

4 結(jié) 論

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