MFC数学模型及MATLAB编程程序.docx

1、MFC数学模型及MATLAB编程程序依据建立的（含电子介体和无电子介体(基于导电生物膜)）流化床MFC数学模型，编写MATLAB程序，模拟结果显示能够和实验结果对应。模型包括传质模型，生化反应模型和电化学模型三个部分，能够完整表述微生物燃料电池的内部的传质反应过程。MATLAB程序如下：无介体流化床MFC模型程序：无电子介体模型以mm,kgCOD,s,A,V,L,为单位x_Su为su微生物体积分数c为微生物浓度改以s为单位clccleark=1;w=1;d=21601;循环次数d-115天内有21601个60秒即21600个时间步长c=zeros(6,6);结果保存在c矩阵e=zeros(6,

2、6);结果保存在e矩阵a=zeros(6,6);作图所需数据结果保存在a矩阵选取特定数据t=60;时间步长sLf=5*10(-5);生物膜初始厚度单位mm=60000;n=7;网格划分n为奇数，生物膜计算区域内划分网格，沿电极方向二维网格垂直于生物膜表面流化床参数A=0.00143;阳极电极面积m2，Vs=3*10(-3);流化床硫化速率msQs=12.566*10(-4)*Vs;流量 m3s电化学参数k_bio=0.05;生物膜电导率smR=91000;总电阻值R_int=1000;电池内阻V_anod=-0.27739;阳极初始电势VE_KA=-0.47;vV_cat=0.25;阴极电势V

3、V_cell=0.5;初始电压值V生物膜内生物组成(可自己定义)x_Su_f降解葡萄糖微生物分率;x_Bu_f降解丁酸微生物分率;x_Pro_f降解丙酸微生物分率;x_Ac_e_f降解乙酸产电微生物分率;x_Ac_AM_f降解乙酸产甲烷微生物分率;x_H_f降解氢气微生物分率;x_I_f惰性微生物分率x_Su_f=0.1625;x_Bu_f=0.1625;x_Pro_f=0.1625;x_Ac_e_f=0.1625;x_Ac_AM_f=0.1625;x_H_f=0.1625;x_I_f=0.025; 生物膜内微生物浓度kgCODm-3生物质的密度为80 kgCOD*m-3，假设体积为1 m3c

4、_Su_f=80*x_Su_f=80*0.1625=13 kgCODm-3c_Su_f=13;c_Bu_f=13;c_Pro_f=13;c_Ac_e_f=13;c_Ac_AM_f=13;c_H_f=13;c_I_f=2; 阳极室微生物初始浓度（根据文献定义）单位kgCOD*m-3c_Su=10(-3);c_Bu=10(-3);c_Pro=10(-3);c_Ac_AM=10(-3);c_H=10(-3);c_I=10(-5); 储槽微生物初始浓度（与阳极室内浓度相同）单位kgCOD*m-3L_c_Su=10(-3);L_c_Bu=10(-3);L_c_Pro=10(-3);L_c_Ac_AM=1

5、0(-3);L_c_H=10(-3);L_c_I=10(-5);物质传质系数 mm2s （根据文献定义）D_Su=5.6944*10(-4);D_Bu=8.6806*10(-4);D_Pro=9.5255*10(-4);D_Ac=1.0891*10(-3);D_CH4=1.8403*10(-4);D_H=4.08*10(-3);储槽中物质初始浓度各物质浓度为kgCODm-3L_Su=3;L_Bu=10(-4);L_Pro=10(-4);L_Ac=10(-4);L_H=4*10(-8);L_CH4=10(-8);对生物膜内划分区域zeros(m,n)即m行，n列 ，把每种物质在生物膜内的分布以区

6、域点的形式表示，即点浓度表示小区域浓度u_Su=zeros(m,n);u_Bu=zeros(m,n);u_Pro=zeros(m,n);u_Ac=zeros(m,n);u_H=zeros(m,n);u_CH4=zeros(m,n);u_v=zeros(m,n);u_v=zeros(m,n)表示电势分布生物膜内物质初始时刻浓度kgCODm-3即区域矩阵第一行的赋值，即 m=1（t=0）,沿生物膜厚度方向物质的浓度相同，反应未开始即无浓度梯度u_Su(1,:)=3;u_Bu(1,:)=10(-4);u_Pro(1,:)=10(-4);u_Ac(1,:)=10(-4);u_H(1,:)=4*10(-

7、8);u_CH4(1,:)=10(-8); 阳极室物质初始浓度kgCODm-3即生物膜靠近阳极室宏观液相区一侧表面浓度等同于液相主题浓度，区域矩阵第一列的赋值，即 n=1,沿生物膜平面方向物质的浓度相同，注意特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面u_Su(:,1)=3;u_Bu(:,1)=10(-4);u_Pro(:,1)=10(-4);u_Ac(:,1)=10(-4);u_H(:,1)=4*10(-8);u_CH4(:,1)=10(-8); 矩阵右边界条件阳极电极表面电势 V 即n=1，电极表面的电势一致 在单个时间步长内不随时间变化特别说明电势计算是从电极表

8、面到生物膜表面即第一列n=1为电极表面u_v(:,1)=V_anod-E_KA;以上为定义变量，赋初值t=0时，生物膜内电势分布特别说明电势计算是从电极表面到生物膜表面for j=2:n-1u_v(1,j)=12*(u_v(1,j-1)+u_v(1,j+1)-(90856*x_Ac_e_f*u_Ac(1,j)(u_Ac(1,j)+0.001)*(1+exp(-38.28*u_v(1,j)+781.76*x_Ac_e_f(1+exp(-38.28*u_v(1,j)*(Lf(n-1)2(2*k_bio);已验证正确性。u_v(1,n)=u_v(1,n-1);end电池整个系统的一个周期内循环，调节

9、d的值来改变计算时间长度while(kd)deltz=Lf*1000(n-1);deltt=t(m-1);ltz单位改为mmtic各物质的网格内的循环计算for i=1:m-1for j=2:n-1u_Su(i+1,j)=u_Su(i,j)+(D_Su*(u_Su(i,j-1)-2*u_Su(i,j)+u_Su(i,j+1)deltz2-0.0278*x_Su_f*u_Su(i,j)(0.5+u_Su(i,j)*deltt; x_Su为su微生物体积分数30*8086400反应系数都改为su_Bu(i+1,j)=u_Bu(i,j)+(D_Bu*(u_Bu(i,j-1)-2*u_Bu(i,j)+

10、u_Bu(i,j+1)deltz2-0.0185*x_Bu_f*u_Bu(i,j)(0.2+u_Bu(i,j)+0.00325*x_Su_f*u_Su(i,j)(0.5+u_Su(i,j)*deltt;u_Pro(i+1,j)=u_Pro(i,j)+(D_Pro*(u_Pro(i,j-1)-2*u_Pro(i,j)+u_Pro(i,j+1)deltz2-0.01204*x_Pro_f*u_Pro(i,j)(0.1+u_Pro(i,j)+0.00675*x_Su_f*u_Su(i,j)(0.5+u_Su(i,j)*deltt; u_Ac(i+1,j)=u_Ac(i,j)+(D_Ac*(u_Ac(

11、i,j-1)-2*u_Ac(i,j)+u_Ac(i,j+1)deltz2-0.00837*x_Ac_e_f*u_Ac(i,j)(0.001+u_Ac(i,j)*(1+exp(-38.28*u_v(i,j)-. . 0.007407*x_Ac_AM_f*u_Ac(i,j)(0.15+u_Ac(i,j)+0.01025*x_Su_f*u_Su(i,j)(0.5+u_Su(i,j)+0.0139*x_Bu_f*u_Bu(i,j)(0.2+u_Bu(i,j)+0.00658667*x_Pro_f*u_Pro(i,j)(0.1+u_Pro(i,j)*deltt; 生物膜电势u_v(i+1,j)=12*(

12、u_v(i,j-1)+u_v(i,j+1)-(90856*x_Ac_e_f*u_Ac(i,j)(u_Ac(i,j)+0.001)*(1+exp(-38.28*u_v(i,j)+781.76*x_Ac_e_f(1+exp(-38.28*u_v(i,j)*(deltz*0.001)2(2*k_bio);Udeltz单位是mm，需要转换成m u_H(i+1,j)=u_H(i,j)+(D_H*(u_H(i,j-1)-2*u_H(i,j)+u_H(i,j+1)deltz2-0.0324*x_H_f*u_H(i,j)(7*10(-6)+u_H(i,j)+0.00477*x_Su_f*u_Su(i,j)(0

13、.5+u_Su(i,j)+0.00348*x_Bu_f*u_Bu(i,j)(0.2+u_Bu(i,j)+0.00497*x_Pro_f*u_Pro(i,j)(0.1+u_Pro(i,j)*deltt;u_CH4(i+1,j)=u_CH4(i,j)+(D_CH4*(u_CH4(i,j-1)-2*u_CH4(i,j)+u_CH4(i,j+1)deltz2+0.00704*x_Ac_AM_f*u_Ac(i,j)(0.15+u_Ac(i,j)+0.0305*x_H_f*u_H(i,j)(7*10(-6)+u_H(i,j)*deltt;end左边界（电极表面）条件u_Su(i+1,n)=u_Su(i+1

14、,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_Bu(i+1,n)=u_Bu(i+1,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_Pro(i+1,n)=u_Pro(i+1,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_Ac(i+1,n)=u_Ac(i+1,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_v(i+1,n)=u_v(i+1,n-1);生物膜表面无电势差特

15、别说明电势计算是从电极表面到生物膜表面即第一列n=1为电极表面 n=n为生物膜表面u_H(i+1,n)=u_H(i+1,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_CH4(i+1,n)=u_CH4(i+1,n-1); 特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面;endtoc计算电流以A,V,mol为单tic计算阳极电极表面电势V_anod=0.187-0.0032653*log(u_Ac(m,1)642.5e-68);单位V更新阳极表面电势u_v(:,1)=V_anod-E_KA;单位VI

16、=(V_cat-V_anod)R;电路电流计算单位AI_1=(V_cat-u_v(m,n)R;电路电流计算单位AV_cell=V_cat-V_anod-I*R_int;电池电压计算单位V独立于体系之外V_cell_1= V_cat-V_anod-I_1*R_int;单位Vtoctic每次循环保存到矩阵中m,kgCOD,s,A,V,L,为单位监测数据c(1,k)=u_Su(m,1);c(2,k)=u_Bu(m,1);c(3,k)=u_Pro(m,1);c(4,k)=u_Ac(m,1);c(5,k)=u_H(m,1);特别说明第一列n=1为生物膜表面即阳极室内溶液主体浓度c(6,k)=u_CH4(

17、m,1);电化学数据保存电流,电压数据到 e 矩阵中e(1,k)=R;总电阻e(2,k)=R-R_int;外接电阻e(3,k)=I*1000;电流mAe(4,k)=V_cell*1000;电压mVe(5,k)=V_cell*I*1000A;功率密度 mWm2e(6,k)=V_anod;单位Ve(10,k)=I_1*1000;电流mAe(20,k)=V_cell_1*1000;电压mV Ie(30,k)=V_cell_1*I*1000A;功率密度 mWm2 I保存数据到 a 矩阵中，数据是以天为单位 用作作图数据if k=1w=k;w=1保存阳极室内各物质浓度a(1,w)=u_Su(m,1);a

18、(2,w)=u_Bu(m,1); a(3,w)=u_Pro(m,1); 别说明第一列n=1为生物膜表面即阳极室内溶液主体浓度a(4,w)=u_Ac(m,1);a(5,w)=u_Su(m,1)+u_Bu(m,1)+u_Pro(m,1)+u_Ac(m,1);a(6,w)=u_H(m,1);a(7,w)=u_CH4(m,1);a(10,w)=R;电阻a(11,w)=R-R_int;外接电阻a(12,w)=double(I)*1000;电流mAa(13,w)=I*(R-R_int)*1000;电压mVa(14,w)=I2*(R-R_int)*10000.00143;率密度mWm2a(20,w)=dou

19、ble(I_1)*1000;电流mAa(21,w)=I_1*(R-R_int)*1000;电压mVa(22,w)=I_12*(R-R_int)*10000.00143;率密度mWm2endif mod(k,1440)=0w=k1440+1;w=2a(1,w)=u_Su(m,1);a(2,w)=u_Bu(m,1);a(3,w)=u_Pro(m,1);a(4,w)=u_Ac(m,1);a(5,w)=u_Su(m,1)+u_Bu(m,1)+u_Pro(m,1)+u_Ac(m,1);a(6,w)=u_H(m,1);a(7,w)=u_CH4(m,1);a(10,w)=R;电阻a(11,w)=R-R_in

20、t;外接电阻a(12,w)=double(I)*1000;电流mAa(13,w)=I*(R-R_int)*1000;电压mVa(14,w)=I2*(R-R_int)*10000.00143;率密度mWm2a(20,w)=double(I_1)*1000;电流mAa(21,w)=I_1*(R-R_int)*1000;电压mVa(22,w)=I_12*(R-R_int)*10000.00143;率密度mWm2end保存数据到C矩阵中：阳极室微生物浓度c(12,k)=c_Su; c(13,k)=c_Bu; c(14,k)=c_Pro;c(15,k)=c_Ac_AM; c(17,k)=c_H; c(1

21、8,k)=c_I;监测微生物膜内各微生物组分c(19,k)=c_Su_f;c(20,k)=c_Bu_f;c(21,k)=c_Pro_f;c(22,k)=c_Ac_AM_f;c(23,k)=c_Ac_e_f;c(24,k)=c_H_f;c(25,k)=c_I_f;微生物体积分率c(26,k)=x_Su_f;c(27,k)=x_Bu_f;c(28,k)=x_Pro_f;c(29,k)=x_Ac_AM_f;c(30,k)=x_Ac_e_f;c(31,k)=x_H_f;c(32,k)=x_I_f;生物膜厚度变化c(33,k)=Lf;电极表面检测c(34,k)=u_Su(m,n);c(35,k)=u_B

22、u(m,n);c(36,k)=u_Pro(m,n);c(37,k)=u_Ac(m,n);c(38,k)=u_H(m,n);c(39,k)=u_CH4(m,n);阳极室组分变化更新右边界条件m,kgCOD,day,A,V,L为单位阳极室内反应速率表达式r_Su=-0.03*c_Su*u_Su(m,1)(0.5+u_Su(m,1);单位kgCOD,day-1 Lf2单位是m 体积是1L 边界条件第一列和第二列数据是一样的r_Bu=-0.02*c_Bu*u_Bu(m,1)(0.2+u_Bu(m,1)+0.00351*c_Su*u_Su(m,1)(0.5+u_Su(m,1);r_Pro=-0.013*

23、c_Pro*u_Pro(m,1)(0.1+u_Pro(m,1)+0.00729*c_Su*u_Su(m,1)(0.5+u_Su(m,1);r_Ac=-0.008*c_Ac_AM*u_Ac(m,1)(0.15+u_Ac(m,1)+0.01107*c_Su*u_Su(m,1)(0.5+u_Su(m,1)+0.01504*c_Bu*u_Bu(m,1)(0.2+u_Bu(m,1)+0.0071136*c_Pro*u_Pro(m,1)(0.1+u_Pro(m,1);r_H=0.00513*c_Su*u_Su(m,1)(0.5+u_Su(m,1)+0.00376*c_Bu*u_Bu(m,1)(0.2+u_

24、Bu(m,1)+0.0052546*c_Pro*u_Pro(m,1)(0.1+u_Pro(m,1);-0.035*c_H*u_H(m,n)(7*10(-6)+u_H(m,n); c(77,k)=-0.035*c_H*u_H(m,1)(7*10(-6)+u_H(m,1)+0.00513*c_Su*u_Su(m,1)(0.5+u_Su(m,1)+0.00376*c_Bu*u_Bu(m,1)(0.2+u_Bu(m,1)+0.0052546*c_Pro*u_Pro(m,1)(0.1+u_Pro(m,1)-4589*10(-8)*(u_H(m,1)-u_H(m,n)Lf*t*0+u_H(m,1); Qs

25、*t*(L_H-u_H(m,n)*1000影响较小r_CH4=0.0076*c_Ac_AM*u_Ac(m,1)(0.15+u_Ac(m,1)+0.0329*c_H*u_H(m,1)(7*10(-6)+u_H(m,1);ltz单位改为mm保存在c矩阵中监测反应速率数据c(44,k)=r_Su*t*0;c(45,k)=r_Bu*t*0;c(46,k)=r_Pro*t*0;c(47,k)=r_Ac*t*0;c(48,k)=r_H*t*0;c(49,k)=r_CH4*t*0;重新计算阳极室内的物质浓度同时更新矩阵右边界条件u_Su(:,1) = r_Su*t0.00186400+Qs*t*(L_Su-

26、u_Su(m,1)*1000+u_Su(m,1);不要忘了改乘时间段。体积为0.001m3 t的单位是s特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_Bu(:,1) = r_Bu*t*0+Qs*t*(L_Bu-u_Bu(m,1)*1000+u_Bu(m,1);特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_Pro(:,1)= r_Pro*t*0+Qs*t*(L_Pro-u_Pro(m,1)*1000+u_Pro(m,1);特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极

27、表面u_Ac(:,1) = r_Ac*t*0+Qs*t*(L_Ac-u_Ac(m,1)*1000+u_Ac(m,1);特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_H(:,1) = r_H*t*0+Qs*t*(L_H-u_H(m,1)*1000+u_H(m,1);特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n为电极表面u_CH4(:,1) = r_CH4*t*0+Qs*t*(L_CH4-u_CH4(m,1)*1000+u_CH4(m,1);特别说明物质浓度计算是从生物膜表面到电极表面即第一列n=1为生物膜表面 n=n

28、为电极表面保存数据在c矩阵中同时监测数据储槽各物质浓度c(52,k)=L_Su;c(53,k)=L_Bu; c(54,k)=L_Pro; c(55,k)=L_Ac; c(56,k)=L_H;c(57,k)=L_CH4;储槽微生物浓度c(62,k)=L_c_Su;c(63,k)=L_c_Bu;c(64,k)=L_c_Pro;c(65,k)=L_c_Ac_AM;c(67,k)=L_c_H;c(68,k)=L_c_I;储槽中的各物质浓度浓度计算同时更新物质浓度L_Su=(-0.03*L_c_Su*L_Su(0.5+L_Su)*2*t*+Qs*t*(u_Su(m,1)-L_Su)*500+L_Su;不

29、要忘了改乘时间段。(-0.03*L_c_Su*L_Su(0.5+L_Su)*2体积为0.002m3 -0.03是体积为0.001m3 前面时间单位是天 t的单位是su_Bu(m,n)是反应之后的阳极室浓度L_Bu=(-0.02*L_c_Bu*L_Bu(0.2+L_Bu)+0.00351*L_c_Su*L_Su(0.5+L_Su)*2*t*+Qs*t*(u_Bu(m,1)-L_Bu)*500+L_Bu;不要忘了改乘时间段。体积为0.002m3L_Pro=(-0.013*L_c_Pro*L_Pro(0.1+L_Pro)+0.00729*L_c_Su*L_Su(0.5+L_Su)*2*t*+Qs*t

30、*(u_Pro(m,1)-L_Pro)*500+L_Pro;不要忘了改乘时间段。体积为0.002m3L_Ac=(-0.008*L_c_Ac_AM*L_Ac(0.15+L_Ac)+0.01107*L_c_Su*L_Su(0.5+L_Su)+0.01504*L_c_Bu*L_Bu(0.2+L_Bu)+0.0071136*L_c_Pro*L_Pro(0.1+L_Pro)*2*t*+Qs*t*(u_Ac(m,1)-L_Ac)*500+L_Ac;L_H=(-0.035*L_c_H*L_H(7*10(-6)+L_H)+0.00513*L_c_Su*L_Su(0.5+L_Su)+0.00376*L_c_Bu

31、*L_Bu(0.2+L_Bu)+0.0052546*L_c_Pro*L_Pro(0.1+L_Pro)*2*t*+Qs*t*(u_H(m,1)-L_H)*500+L_H;L_CH4=(0.0076*L_c_Ac_AM*L_Ac(0.15+L_Ac)+0.0329*L_c_H*L_H(7*10(-6)+L_H)*2*t*+Qs*t*(u_CH4(m,1)-L_CH4)*500+L_CH4;更新生物膜区域矩阵上边界条件 即作下一时间步长的初始值u_Su(1,:)=u_Su(m,:);u_Bu(1,:)=u_Bu(m,:);u_Pro(1,:)=u_Pro(m,:);u_Ac(1,:)=u_Ac(m,:);u_v(1,:) = u_v(m,:);电势初值更新u_H(1,:) = u_H(m,:);u_CH4(1,:) = u_CH4(m,:);生物膜微生物生长模型生物膜厚度生长Lf=(Lf*(3*x_Su_f*u_Su(m,(n-1)2)(0.5+u_Su(m,(n-1)2)+1.2*x_Bu_f*u_Bu(m,(n-1)2)(0.2+u_Bu(m,(n-1)2)+0.52*x_Pro_f*u_Pro(m,(n-1)2)(0.1+u_Pr