Matlab在求解擴(kuò)散系統(tǒng)之濃度分布中的應(yīng)用

讀書報(bào)告Matlab在求解擴(kuò)散系統(tǒng)之濃度分布

中的應(yīng)用池雨

、問題的提出管中儲(chǔ)放靜止液體B，高度為L=10cm,放置于充滿A氣體的環(huán)境中。假設(shè)與B液體接觸面之濃度為C=0.01mol/m3,且此濃度不A0隨時(shí)間改變而改變，即在操作時(shí)間內(nèi)（h=10天）維持定值。氣體A在液體B中之?dāng)U散系數(shù)為DAB=2X109m2/s。試決定A與B不發(fā)生反應(yīng);AB情況下，氣體A溶于液體B中之流通量（flux）。參考如圖所示的裝TloTlo￡th氣棒血在液休B中的擴(kuò)浪二、知識(shí)背景Fick第一定律：實(shí)驗(yàn)表明，在穩(wěn)態(tài)擴(kuò)散的條件下，單位時(shí)間內(nèi)通過垂直于擴(kuò)散方向的單位截面面積的擴(kuò)散物質(zhì)的通量與濃度梯度成dc正比。數(shù)學(xué)表達(dá)式為：J=-DdxFick第二定律：根據(jù)質(zhì)量平衡關(guān)系即在微小體積中積存的物質(zhì)dcj—j（留入的物質(zhì)量）J1—（留出的物質(zhì)量）J2得出丁=亠」12dtdx因此Fick第二定律的數(shù)學(xué)表達(dá)為：竺=D竺dtdx2三、問題求解根據(jù)題意不同時(shí)間t和距界面厚度不同處x的濃度C=f（z,t）;因氣體A與液體B不發(fā)生反應(yīng)，故其擴(kuò)散現(xiàn)象的質(zhì)量平衡方程根據(jù)Fick

第二定律。-.■-■依題意，其初始及邊界條件為:I.C.Cao（Z,0）=0,Z>0A0QCNAZ（t）二一DABBC.CA（0,t）=CA0,m；占lzJNAZ（t）二一DAB在獲得濃度分布后即可應(yīng)用Fick第一定律求得流通量即四、matlab程序設(shè)計(jì)偏微分方程（PartialDifferentialEquation，簡稱PDE）就是涉及到兩個(gè)自變量以上的微分方程。在化學(xué)工程領(lǐng)域，為了更好的進(jìn)行過程設(shè)計(jì)、優(yōu)化和控制，經(jīng)常需要了解化工設(shè)備（如反應(yīng)器）中的溫度、濃度和速度在不同空間上的分布以及隨時(shí)間的動(dòng)態(tài)變化規(guī)律，因此涉及到許多偏微分方程的問題。Matlab函數(shù)pdepe（）可用于求解偏微分方程，模型為:壬du.du匚仏?。?CJlC7-v（x刃耳,r,虬r心）fit（Jjlor用以解含上述初始值及邊界值條件的偏微分方程MATLAB命令pdepe的用法如下:sol=icfuH.bcfun.tspan.options）若要獲得特定位置及時(shí)間下的解，可配合以pdeval命令。使用格式如下：uoul*t/reMfiZx]=pdeval（m加旳九ui,xaut）第一步：Fick化為標(biāo)準(zhǔn)形式后即f(第一步：Fick化為標(biāo)準(zhǔn)形式后即f(Z，t，CA，二0dc1c(z,t,c,A)二AozDABOc第二步：編寫偏微分方程的系數(shù)向量函數(shù)：function[c,f,s]=ex20_3_2pdefun(z,t,CA,dCAdt)c=1/D;f=dCAdt;s=0;第三步：編寫起始條件：functionCA_i=ex20_3_2ic(z)CA_i=0;第四步：編寫邊界條件：function[pl,ql,pr,qr]=ex20_3_2bc(zl,CAl,zr,CAr,t)globalDABkCA0pl=CAl-CA0;ql=0;pr=0;qr=1/DAB;第五步：取點(diǎn)。例如t=linspace(0,h,100);z=linspace(0,L,10);第六步利用pdepe求解。m=0;sol=pdepe(m,@ex20_3_2pdefuna.@ex20_3_2ic,@ex20_3_2bc,z,t)

sol=pdepe(m,@ex20_3_2pdefuna.@ex20_3_2ic,@ex20_3_2bc,z,t)Sol(j,k,i)第一維代表時(shí)間t,第二維代表空間位置z,第三維代表解向量u的第i個(gè)元素，即ui=sol(:，:,i)第七步顯示結(jié)果。u=sol(:,:,l);surf(x,t,u)title('標(biāo)題')xlabel('位置')ylabel("conc.(mol/mA3)'')zlabel('u')若要獲得特定位置及時(shí)間下的解，可配合以pdeval命令。使用格式如下L本題的程序編寫為：fori=l:length(t)[CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),0);NAz(i)=-dCAdz_i*DAB;End附錄（matlab程序）：functionex20_3_2%擴(kuò)散系統(tǒng)之濃度分布clearclcglobalDABkCA0%給定數(shù)據(jù)a=input('請(qǐng)輸入要求通量處厚度a=')CA0=0.01;L=0.1;DAB=2e-9;h=10*24*3600;%取點(diǎn)t=linspace(0,h,100);z=linspace(0,L,10);m=0;sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t)sol=pdepe(m,@ex20_3_2pdefuna,@ex20_3_2ic,@ex20_3_2bc,z,t)CA=sol(:,:,l);fori=l:length(t)[CA_i,dCAdz_i]=pdeval(m,z,CA(i,:),a);NAz(i)=-dCAdz_i*DAB;endsubplot(211)surf(z,t/(24*3600),CA)title('case(a)')xlabel('length(m)')ylabel('time(day)')zlabel('conc.(mol/mA3)')subplot(212)plot(t/(24*3600),NAz'*24*3600)xlabel('time(day)')ylabel('flux(mol/mA2.day)')%PDE函數(shù)*L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L*function[c,f,s]=ex20_3_2pdefuna(z,t,CA,dCAdz)globalDABkCA0c=1/DAB;f=dCAdz;s=0;%初始條件函數(shù)*L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L*f^.functionCA_i=ex20_3_2ic(z)CA_i=0;Z"V/*L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L*\^f^.%邊界條件函數(shù)Z"V/*L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L**L*