工程數(shù)學(xué)作業(yè)第十三次方健(1217)

1、姓名：方健 學(xué)號：652081701073 專業(yè)：化學(xué)工程 班級：化工1503問題：利用MATLAB對與化學(xué)工程學(xué)科的實際問題進行解決，比如使用正交配置法求解有軸向擴散的固定床反應(yīng)器中的催化劑的溫度和濃度分布，柱形固定床反應(yīng)器中的催化的溫度和濃度方程為：-Tr=1=BiwT1,z-Twz, -cr=1=0Tr,0=T0, cr,0=c0其中R(c,T)為催化反應(yīng)的速率方程，其形式為：分析：從上面的問題可以知道，以上是一個偏微分方程的數(shù)值求解問題，由于該方程含有初值條件，還有左右邊界條件，所以應(yīng)該使用pdepe函數(shù)進行求解，于是先學(xué)習(xí)了pdepe函數(shù)的具體用法如下：pdefun(函數(shù)格式描述)：

2、c,f,s=pdefun(x,t,u,du)pdebc（邊界條件描述）：pa,qa,pb,qb=pdebc(x,t,u,du) l:左邊界，r:右邊界調(diào)用格式：sol=pdepe(m, pdepe, icfun, bcfun, xmesh, tspan)sol=pdepe(m, pdepe, icfun, bcfun, xmesh, tspan, options)輸入?yún)?shù)：1. m表示定解問題的對稱性，m=0, 1, 2分別代表平板、圓柱和球形。2. pdefun是描述PDE問題的函數(shù)，輸出變量為標(biāo)準(zhǔn)形式中的c, f和s，其格式為c, f,s=pdefunc(x,t,u,dudx)。輸入變量x

3、和t為標(biāo)量，u和dudx是向量。u和dudx分別是問題的解u(x,t)和它對x的偏導(dǎo)數(shù)的近似。c, f, s是列向量。3. icfun描述定解問題初始條件的函數(shù)，其格式為u=icfun(x)。icfun計算和返回解的初始值。4. bcfun是描述定解問題邊界條件的函數(shù)，格式為pl, ql, pr, qr=bcfun(xl, ul, xr, ur, t)。其中ul是在左邊界xl=a處的近似解，ur是在右邊界xr=b處的近似解。pl和ql是p和q在xl處的列向量值，同樣pr和qr是p和q在xr處的列向量值5. xmesh是空間變量x的網(wǎng)點向量，pdepe不會自動選擇。xmesh=x0, x1, ,

4、 xn，滿足x0<x1<<xn，且xmesh的長度必須大于3，xmesh(1)=a和xmesh(end)=b。pdepe的求解效率與xmesh的選擇好壞關(guān)系很大。通常在梯度較大處應(yīng)加密網(wǎng)格。通過上面的學(xué)習(xí)，應(yīng)該先將原方程組進行變形，將溫度變量T使用v1代替,用v2代替濃度c，根據(jù)以上的函數(shù)的調(diào)用特點，以及原方程組的組成，可以將原方程組改寫成： 其中：于是編寫偏微分方程的函數(shù)如下，fjpdepe.m文件如下：function c,f,s=fjpdepe(r,z,v,dv)c=1;1;lamb=20;F=(1-v(2)*exp(lamb-lamb/v(1);s=F*0.2;0.3

5、;f=r*dv(1);r*dv(2);end將邊界條件改寫為：下邊界： 上邊界： 下邊界在r=0處和下邊界在r=1處，編寫fjbc.m文件如下：function pl,ql,pr,qr=fjbc(rava,rb,vb,r,z)% l表示下邊界，r表示上邊界pl=0;0;ql=1;1;pr=r*(vb(2)-0.92);0;qr=1;1;初值條件改寫為： 下面編寫初試條件fjic.m文件如下：function v0=fjic(r)v0=1,0;最后進行計算函數(shù)主函數(shù)的編寫，編寫fj_pdepe.m文件如下：首先給出繪制圖形的區(qū)域，將微分方程組的結(jié)果繪制稱為圖象，其中定義r的范圍在（0,2），z的

6、范圍在（0.3）于是編寫程序如下：clear allclcr=0:0.05:2; z=0:0.05:3;%確定繪圖區(qū)域m=1;%其中m為問題之對稱參數(shù)，求解方程組類型sol = pdepe(m,fjpdepe,fjic,fjbc,r,z);v1 = sol(:,:,1);%第一個狀態(tài)之?dāng)?shù)值解輸出v2 = sol(:,:,2);%第二個狀態(tài)之?dāng)?shù)值解輸出%*%繪圖程序figure(1) %繪制T(r,z)圖形surf(r,z,v1)title('T(r,z)')xlabel(' r')ylabel('z')figure(2) %繪制T(r,z)圖形s

7、urf(r,z,v2)title('c(r,z)')xlabel(' r')ylabel(' z')figure(3) %繪制溫度T隨高度z變化圖形plot(z,v1)title('z,T')xlabel('z')ylabel('T')figure(4) %繪制濃度c隨高度z變化分布圖形plot(z,v2)title('z,c')xlabel('z')ylabel('c')于是在MATLAB的command window中輸入>>fj_pd

8、epe運行結(jié)果如下：濃度隨床層高度變化曲線溫度隨床層高度變化曲線圖形3圖形4結(jié)果分析與討論：從上面的結(jié)果中可以看出，利用pdepe函數(shù)以及對于微分方程的邊界條件進行分析，可以得到了反應(yīng)過程中的濃度以及溫度隨著高度的變化過程，從第一張圖中可以看出，濃度隨著高度變化的曲線，從圖中可以看出來，在高度等于0.3處出現(xiàn)突變，之后濃度很快上升達到了最大，其中隨著床層的高度增加溫度也是隨著下降，從圖像中溫度的變化可以知道該反應(yīng)是放熱反應(yīng)，所以，但是隨著后面的濃度減少導(dǎo)致溫度下降。MATLAB 的PDE求解器pdepe的算法，主要是將原來的橢圓型和拋物線型偏微分方程轉(zhuǎn)化為一組常微分方程。此轉(zhuǎn)換的過程是基于使用者所指定的mesh點，以二階空間離散化(spatial discretization)技術(shù)，然后以ode15s 的指令求解。采用ode15s 的ode 解法，主要是因為在離散化的過程中，橢圓型偏微分方程被轉(zhuǎn)化為一組代數(shù)方程，而