線上求解法最優近似格式的探討及應用

2022-04-12 01:42:28劉陽,李凌

沈陽化工大學學報 2022年6期

劉陽, 李凌

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

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

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

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

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先選取空間網格點zi(i=1,2,…,n),再用泰勒級數進行展開,則一階偏微分項的三點中心近似格式為

Δz=zi-zi-1.

(4)

一階偏微分項的五點偏心逆風近似格式為

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

(5)

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

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

(6)

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

1.2 時間變量的積分求解

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

2 基于MATLAB的仿真實驗

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

2.1 一維Burgers方程

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

(7)

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

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軟件工具下編寫程序進行仿真實驗,再根據所得圖像進行分析.其中μ=0.001;t=0,0.1,0.2,…,1,單位為min;網格數n=201;解析解輸出圖用虛線顯示,線上求解法近似格式用實線顯示.先對式(7)中的二階偏微分項采用五點中心格式進行離散近似,再對一階偏微分項分別進行五點偏心逆風格式近似、兩點逆風格式近似、三點逆風格式近似和七點中心格式近似,所得圖像分別如圖1～圖4所示.

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

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

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

圖3 三點逆風格式近似一階偏微分項與解析解對比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)中的一階偏微分項給定五點偏心逆風格式進行離散近似,再對二階偏微分項分別進行五點中心格式近似、兩點逆風格式近似、三點逆風格式近似,其中五點中心格式近似后圖像仍為圖1,其余兩種格式所得圖像分別如圖5、圖6所示.

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

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

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

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

3 MOL在化工過程中的應用

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

(9)

(10)

(11)

(12)

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

根據以上結論驗證該動力學模型,以五點偏心逆風格式近似方程一階偏微分項,五點中心格式近似方程二階偏微分項,所得反應器內部的溫度分布曲線如圖7所示.將圖7作為標準圖像.

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

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