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

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

《MFC数学模型及MATLAB编程程序.docx》由会员分享，可在线阅读，更多相关《MFC数学模型及MATLAB编程程序.docx（14页珍藏版）》请在冰豆网上搜索。

MFC数学模型及MATLAB编程程序

依据建立的（含电子介体和无电子介体（基于导电生物膜））流化床MFC数学模型，编写MATLAB程序，模拟结果显示能够和实验结果对应。

模型包括传质模型，生化反应模型和电化学模型三个部分，能够完整表述微生物燃料电池的内部的传质反应过程。

MATLAB程序如下：

无介体流化床MFC模型程序：

无电子介体模型以mm,kgCOD,s,A,V,L,为单位x_Su为su微生物体积分数c为微生物浓度改以s为单位

clc

clear

k=1;w=1;d=21601;循环次数d-115天内有21601个60秒即21600个时间步长

c=zeros（6,6）;结果保存在c矩阵

e=zeros（6,6）;结果保存在e矩阵

a=zeros（6,6）;作图所需数据结果保存在a矩阵选取特定数据

t=60;时间步长s

Lf=5*10^（-5）;生物膜初始厚度单位m

m=60000;n=7;网格划分n为奇数，生物膜计算区域内划分网格，沿电极方向二维网格垂直于生物膜表面

流化床参数

A=0.00143;阳极电极面积m2，

Vs=3*10^（-3）;流化床硫化速率ms

Qs=12.566*10^（-4）*Vs;流量m3s

电化学参数

k_bio=0.05;生物膜电导率sm

R=91000;总电阻值

R_int=1000;电池内阻

V_anod=-0.27739;阳极初始电势V

E_KA=-0.47;v

V_cat=0.25;阴极电势V

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生物质的密度为80kgCOD*m-3，假设体积为1m3c_Su_f=80*x_Su_f=80*0.1625=13kgCODm-3

c_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-3

c_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-3

L_c_Su=10^（-3）;L_c_Bu=10^（-3）;L_c_Pro=10^（-3）;L_c_Ac_AM=10^（-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-3

L_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列，把每种物质在生物膜内的分布以区域点的形式表示，即点浓度表示小区域浓度

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^（-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，电极表面的电势一致在单个时间步长内不随时间变化特别说明电势计算是从电极表面到生物膜表面即第一列n=1为电极表面

u_v（:

1）=V_anod-E_KA;以上为定义变量，赋初值t=0时，生物膜内电势分布特别说明电势计算是从电极表面到生物膜表面

forj=2:

n-1

u_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

电池整个系统的一个周期内循环，调节d的值来改变计算时间长度

while（k

deltz=Lf*1000（n-1）;deltt=t（m-1）;ltz单位改为mm

tic

各物质的网格内的循环计算

fori=1:

m-1

forj=2:

n-1

u_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））deltz^2-0.0278*x_Su_f*u_Su（i,j）（0.5+u_Su（i,j）））*deltt;x_Su为su微生物体积分数30*8086400反应系数都改为s

u_Bu（i+1,j）=u_Bu（i,j）+（D_Bu*（u_Bu（i,j-1）-2*u_Bu（i,j）+u_Bu（i,j+1））deltz^2-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））deltz^2-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（i,j-1）-2*u_Ac（i,j）+u_Ac（i,j+1））deltz^2-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*（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））deltz^2-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.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））deltz^2+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,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）;生物膜表面无电势差特别说明电势计算是从电极表面到生物膜表面即第一列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为电极表面;

end

toc

计算电流以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;单位V

I=（V_cat-V_anod）R;电路电流计算单位A

I_1=（V_cat-u_v（m,n））R;电路电流计算单位A

V_cell=V_cat-V_anod-I*R_int;电池电压计算单位V

独立于体系之外

V_cell_1=V_cat-V_anod-I_1*R_int;单位V

toc

tic

每次循环保存到矩阵中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（m,1）;

电化学数据保存电流,电压数据到e矩阵中

e（1,k）=R;总电阻

e（2,k）=R-R_int;外接电阻

e（3,k）=I*1000;电流mA

e（4,k）=V_cell*1000;电压mV

e（5,k）=V_cell*I*1000A;功率密度mWm2

e（6,k）=V_anod;单位V

e（10,k）=I_1*1000;电流mA

e（20,k）=V_cell_1*1000;电压mVI

e（30,k）=V_cell_1*I*1000A;功率密度mWm2I

保存数据到a矩阵中，数据是以天为单位用作作图数据

ifk==1

w=k;w=1

保存阳极室内各物质浓度

a（1,w）=u_Su（m,1）;a（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;电流mA

a（13,w）=I*（R-R_int）*1000;电压mV

a（14,w）=I^2*（R-R_int）*10000.00143;率密度mWm2

a（20,w）=double（I_1）*1000;电流mA

a（21,w）=I_1*（R-R_int）*1000;电压mV

a（22,w）=I_1^2*（R-R_int）*10000.00143;率密度mWm2

end

ifmod（k,1440）==0

w=k1440+1;w=2

a（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_int;外接电阻

a（12,w）=double（I）*1000;电流mA

a（13,w）=I*（R-R_int）*1000;电压mV

a（14,w）=I^2*（R-R_int）*10000.00143;率密度mWm2

a（20,w）=double（I_1）*1000;电流mA

a（21,w）=I_1*（R-R_int）*1000;电压mV

a（22,w）=I_1^2*（R-R_int）*10000.00143;率密度mWm2

end

保存数据到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（18,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_Bu（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-1Lf2单位是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*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_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*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-u_Su（m,1））*1000+u_Su（m,1）;不要忘了改乘时间段。

体积为0.001m3t的单位是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为电极表面

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为电极表面

保存数据在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;不要忘了改乘时间段。

（-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.002m3

L_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*（u_Pro（m,1）-L_Pro）*500+L_Pro;不要忘了改乘时间段。

体积为0.002m3

L_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*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