依据建立的（含电子介体和无电子介体(基于导电生物膜)）流化床MFC数学模型，编写MATLAB程序，模拟结果显示能够和实验结果对应。模型包括传质模型，生化反应模型和电化学模型三个部分，能够完整表述微生物燃料电池的内部的传质反应过程。
MATLAB程序如下：无介体MFC模型程序：
%%%%
%无电子介体模型%%%%%%%%%%%%%%%%%以mm,kgCOD,s,A,V,L,为单位%%%%%%%x_Su为su微生物体积分数 c为微生物浓度 改以s为单位
clc
clear
k=1;w=1;d=21601;%%%循环次数d-1%%15天内有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);%%%%%%%流化床硫化速率m/s
Qs=12.566*10^(-4)*Vs;%%%%%%流量 m3/s
%%电化学参数
k_bio=0.05;%%%%生物膜电导率s/m
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%%%%%%生物质的密度为80 kgCOD*m-3，假设体积为1 m3%%%%%%%c_Su_f=80*x_Su_f=80*0.1625=13 kgCODm-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);%%物质传质系数 %mm2/s %%%%%%（根据文献定义）
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时，生物膜内电势分布%特别说明电势计算是从电极表面到生物膜表面!!!!
for j=2:n-1
u_v(1,j)=1/2*(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<d)
deltz=Lf*1000/(n-1);deltt=t/(m-1);%%%%%%%ltz单位改为mm
tic
%%%%%%%%%%各物质的网格内的循环计算————%%%%%%
for i=1:m-1
for j=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*80/86400反应系数都改为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)=1/2*(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)/64/2.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*1000/A;%%%%%%%%%功率密度 %%mW/m2
e(6,k)=V_anod;%%%%%%%单位V
e(10,k)=I_1*1000;%%%%%%%%%电流mA
e(20,k)=V_cell_1*1000;%%%%%%%%%电压mV %%%%I
e(30,k)=V_cell_1*I*1000/A;%%%%%%%%%功率密度 %%mW/m2 %%%I
%%%%%%%%%保存数据到 a 矩阵中，数据是以天为单位 用作作图数据%%%%
if k==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)*1000/0.00143;%%%%%率密度mW/m2
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)*1000/0.00143;%%%%%率密度mW/m2
end
if mod(k,1440)==0
w=k/1440+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)*1000/0.00143;%%%%%率密度mW/m2
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)*1000/0.00143;%%%%%率密度mW/m2
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-1 Lf/2单位是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*1000/86400+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*1000/86400;c(45,k)=r_Bu*t*1000/86400;c(46,k)=r_Pro*t*1000/86400;c(47,k)=r_Ac*t*1000/86400;c(48,k)=r_H*t*1000/86400;c(49,k)=r_CH4*t*1000/86400;
%%%重新计算阳极室内的物质浓度%%%%同时更新矩阵右边界条件%%%%%%%%%
u_Su(:,1) = r_Su*t/0.001/86400+Qs*t*(L_Su-u_Su(m,1))*1000+u_Su(m,1);%%%不要忘了改乘时间段。体积为0.001m3 t的单位是s%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Bu(:,1) = r_Bu*t*1000/86400+Qs*t*(L_Bu-u_Bu(m,1))*1000+u_Bu(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Pro(:,1)= r_Pro*t*1000/86400+Qs*t*(L_Pro-u_Pro(m,1))*1000+u_Pro(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Ac(:,1) = r_Ac*t*1000/86400+Qs*t*(L_Ac-u_Ac(m,1))*1000+u_Ac(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_H(:,1) = r_H*t*1000/86400+Qs*t*(L_H-u_H(m,1))*1000+u_H(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 %n=n为电极表面
u_CH4(:,1) = r_CH4*t*1000/86400+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*500/86400+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的单位是s%%%%%u_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*500/86400+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*500/86400+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*500/86400+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*500/86400+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*500/86400+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_Pro(m,(n-1)/2))+0.4*x_Ac_AM_f*u_Ac(m,(n-1)/2)/(0.15+u_Ac(m,(n-1)/2))+0.615*x_Ac_e_f*u_Ac(m,(n-1)/2)/((0.001+u_Ac(m,(n-1)/2))*(1+exp(-38.28*u_v(m,(n-1)/2)))))-0.08*Lf)*t/86400+Lf;%%%%%%%%%t的单位是s%%%%%lf单位是m%%%%%%%%%%%%%%生物膜内生物的组成分布%%%%%%%
%%%%%生物膜内微生物浓度用作接下来算生物的体积分率%%%%%%假设体积为1m3单位为Kg CODx m-3
c_Su_f=((3*u_Su(m,(n-1)/2)/(0.5+u_Su(m,(n-1)/2))-0.02)*x_Su_f)*80*t/86400+c_Su_f;
c_Bu_f=((1.2*u_Bu(m,(n-1)/2)/(0.2+u_Bu(m,(n-1)/2))-0.02)*x_Bu_f)*80*t/86400+c_Bu_f;
c_Pro_f=((0.52*u_Pro(m,(n-1)/2)/(0.1+u_Pro(m,(n-1)/2))-0.02)*x_Pro_f)*80*t/86400+c_Pro_f;
c_Ac_e_f=((0.615*u_Ac(m,(n-1)/2)/((0.001+u_Ac(m,(n-1)/2))*(1+exp(-38.28*u_v(m,(n-1)/2))))-0.02)*x_Ac_e_f)*80*t/86400+c_Ac_e_f;
c_Ac_AM_f=((0.4*u_Ac(m,(n-1)/2)/(0.15+u_Ac(m,(n-1)/2))-0.02)*x_Ac_AM_f)*80*t/86400+c_Ac_AM_f;
c_H_f=((2.1*u_H(m,(n-1)/2)/(7*10^(-6)+u_H(m,(n-1)/2))-0.02)*x_H_f)*80*t/86400+c_H_f;
%%%%%惰性微生物计算%%%%!!!!！！！计算过程中发现降解氢气产甲烷微生物生长数据异常，故惰性微生物计算省略此部分
c_I_f=0.02*(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f)*t/86400+c_I_f;
%%%%%%%%微生物活性组分分率计算%%%%%%%%
x_Su_f=c_Su_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Bu_f=c_Bu_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Pro_f=c_Pro_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Ac_e_f=c_Ac_e_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Ac_AM_f=c_Ac_AM_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_H_f=c_H_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_I_f=c_I_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
%%阳极室微生物生长模型%%%%%%%%%%%%%%%%%%%%%%%%%%%阳极室微生物生长同时更新微生物的浓度数据用作下一时间步长的循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c_Su=(3*c_Su*u_Su(m,1)/(0.5+u_Su(m,1))-0.02*c_Su)*t/86400+Qs*t*(L_c_Su-c_Su)*1000+c_Su;
c_Bu=(1.2*c_Bu*u_Bu(m,1)/(0.2+u_Bu(m,1))-0.02*c_Bu)*t/86400+Qs*t*(L_c_Bu-c_Bu)*1000+c_Bu;
c_Pro=(0.52*c_Pro*u_Pro(m,1)/(0.1+u_Pro(m,1))-0.02*c_Pro)*t/86400+Qs*t*(L_c_Pro-c_Pro)*1000+c_Pro;
c_Ac_AM=(0.4*c_Ac_AM*u_Ac(m,1)/(0.15+u_Ac(m,1))-0.02*c_Ac_AM)*t/86400+Qs*t*(L_c_Ac_AM-c_Ac_AM)*1000+c_Ac_AM;
c_H=(2.1*c_H*u_H(m,1)/(7*10^(-6)+u_H(m,1))-0.02*c_H)*t/86400+Qs*t*(L_c_H-c_H)*1000+c_H;
c_I=0.02*(c_Su+c_Bu+c_Pro+c_Ac_AM)*t/86400+Qs*t*(L_c_I-c_I)*1000+c_I;
%%%%储槽微生物生长模型%%%%%%同上%%%%%%%%%%%%%
L_c_Su=(3*L_c_Su*L_Su/(0.5+L_Su)-0.02*L_c_Su)*t/86400+Qs*t*(c_Su-L_c_Su)*500+L_c_Su;%%%%%%%积累等于进入减去输出%%%%%%Qs*t*(c_Su-L_c_Su)*500除以体积0.002m3
L_c_Bu=(1.2*L_c_Bu*L_Bu/(0.2+L_Bu)-0.02*L_c_Bu)*t/86400+Qs*t*(c_Bu-L_c_Bu)*500+L_c_Bu;
L_c_Pro=(0.52*L_c_Pro*L_Pro/(0.1+L_Pro)-0.02*L_c_Pro)*t/86400+Qs*t*(c_Pro-L_c_Pro)*500+L_c_Pro;
L_c_Ac_AM=(0.4*L_c_Ac_AM*L_Ac/(0.15+L_Ac)-0.02*L_c_Ac_AM)*t/86400+Qs*t*(c_Ac_AM-L_c_Ac_AM)*500+L_c_Ac_AM;
L_c_H=(2.1*L_c_H*L_H/(7*10^(-6)+L_H)-0.02*L_c_H)*t/86400+Qs*t*(c_H-L_c_H)*500+L_c_H;
L_c_I=0.02*(L_c_I+L_c_Bu+L_c_Pro+L_c_Ac_AM)*t/86400+Qs*t*(c_I-L_c_I)*500+L_c_I;
k=k+1;
toc
end
c;
x=1:n;y=1:m;
figure(1)
subplot(1,7,1)
plot(u_Su(m,:),'r-')
title('Su变化')
subplot(1,7,2)
plot(u_Bu(m,:),'b--')
title('Bu变化')
subplot(1,7,3)
plot(u_Pro(m,:),'g-')
title('Pro变化')
subplot(1,7,4)
plot(u_Ac(m,:),'r-')
title('Ac变化')
subplot(1,7,5)
plot(u_H(m,1:n),'b--')
title('H变化')
subplot(1,7,6)
plot(u_CH4(m,:),'g-')
title('CH4变化')
subplot(1,7,7)
plot(u_v(m,:),'g-')
title('沿生物膜厚度方向电势变化')
接下来无介体模型的Word版和含电子介体的模型Word和MATLAB程序。。
MATLAB程序如下：无介体MFC模型程序：
%%%%
%无电子介体模型%%%%%%%%%%%%%%%%%以mm,kgCOD,s,A,V,L,为单位%%%%%%%x_Su为su微生物体积分数 c为微生物浓度 改以s为单位
clc
clear
k=1;w=1;d=21601;%%%循环次数d-1%%15天内有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);%%%%%%%流化床硫化速率m/s
Qs=12.566*10^(-4)*Vs;%%%%%%流量 m3/s
%%电化学参数
k_bio=0.05;%%%%生物膜电导率s/m
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%%%%%%生物质的密度为80 kgCOD*m-3，假设体积为1 m3%%%%%%%c_Su_f=80*x_Su_f=80*0.1625=13 kgCODm-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);%%物质传质系数 %mm2/s %%%%%%（根据文献定义）
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时，生物膜内电势分布%特别说明电势计算是从电极表面到生物膜表面!!!!
for j=2:n-1
u_v(1,j)=1/2*(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<d)
deltz=Lf*1000/(n-1);deltt=t/(m-1);%%%%%%%ltz单位改为mm
tic
%%%%%%%%%%各物质的网格内的循环计算————%%%%%%
for i=1:m-1
for j=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*80/86400反应系数都改为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)=1/2*(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)/64/2.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*1000/A;%%%%%%%%%功率密度 %%mW/m2
e(6,k)=V_anod;%%%%%%%单位V
e(10,k)=I_1*1000;%%%%%%%%%电流mA
e(20,k)=V_cell_1*1000;%%%%%%%%%电压mV %%%%I
e(30,k)=V_cell_1*I*1000/A;%%%%%%%%%功率密度 %%mW/m2 %%%I
%%%%%%%%%保存数据到 a 矩阵中，数据是以天为单位 用作作图数据%%%%
if k==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)*1000/0.00143;%%%%%率密度mW/m2
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)*1000/0.00143;%%%%%率密度mW/m2
end
if mod(k,1440)==0
w=k/1440+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)*1000/0.00143;%%%%%率密度mW/m2
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)*1000/0.00143;%%%%%率密度mW/m2
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-1 Lf/2单位是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*1000/86400+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*1000/86400;c(45,k)=r_Bu*t*1000/86400;c(46,k)=r_Pro*t*1000/86400;c(47,k)=r_Ac*t*1000/86400;c(48,k)=r_H*t*1000/86400;c(49,k)=r_CH4*t*1000/86400;
%%%重新计算阳极室内的物质浓度%%%%同时更新矩阵右边界条件%%%%%%%%%
u_Su(:,1) = r_Su*t/0.001/86400+Qs*t*(L_Su-u_Su(m,1))*1000+u_Su(m,1);%%%不要忘了改乘时间段。体积为0.001m3 t的单位是s%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Bu(:,1) = r_Bu*t*1000/86400+Qs*t*(L_Bu-u_Bu(m,1))*1000+u_Bu(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Pro(:,1)= r_Pro*t*1000/86400+Qs*t*(L_Pro-u_Pro(m,1))*1000+u_Pro(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_Ac(:,1) = r_Ac*t*1000/86400+Qs*t*(L_Ac-u_Ac(m,1))*1000+u_Ac(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 n=n为电极表面
u_H(:,1) = r_H*t*1000/86400+Qs*t*(L_H-u_H(m,1))*1000+u_H(m,1);%！特别说明物质浓度计算是从生物膜表面到电极表面！！即第一列n=1为生物膜表面 %n=n为电极表面
u_CH4(:,1) = r_CH4*t*1000/86400+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*500/86400+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的单位是s%%%%%u_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*500/86400+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*500/86400+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*500/86400+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*500/86400+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*500/86400+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_Pro(m,(n-1)/2))+0.4*x_Ac_AM_f*u_Ac(m,(n-1)/2)/(0.15+u_Ac(m,(n-1)/2))+0.615*x_Ac_e_f*u_Ac(m,(n-1)/2)/((0.001+u_Ac(m,(n-1)/2))*(1+exp(-38.28*u_v(m,(n-1)/2)))))-0.08*Lf)*t/86400+Lf;%%%%%%%%%t的单位是s%%%%%lf单位是m%%%%%%%%%%%%%%生物膜内生物的组成分布%%%%%%%
%%%%%生物膜内微生物浓度用作接下来算生物的体积分率%%%%%%假设体积为1m3单位为Kg CODx m-3
c_Su_f=((3*u_Su(m,(n-1)/2)/(0.5+u_Su(m,(n-1)/2))-0.02)*x_Su_f)*80*t/86400+c_Su_f;
c_Bu_f=((1.2*u_Bu(m,(n-1)/2)/(0.2+u_Bu(m,(n-1)/2))-0.02)*x_Bu_f)*80*t/86400+c_Bu_f;
c_Pro_f=((0.52*u_Pro(m,(n-1)/2)/(0.1+u_Pro(m,(n-1)/2))-0.02)*x_Pro_f)*80*t/86400+c_Pro_f;
c_Ac_e_f=((0.615*u_Ac(m,(n-1)/2)/((0.001+u_Ac(m,(n-1)/2))*(1+exp(-38.28*u_v(m,(n-1)/2))))-0.02)*x_Ac_e_f)*80*t/86400+c_Ac_e_f;
c_Ac_AM_f=((0.4*u_Ac(m,(n-1)/2)/(0.15+u_Ac(m,(n-1)/2))-0.02)*x_Ac_AM_f)*80*t/86400+c_Ac_AM_f;
c_H_f=((2.1*u_H(m,(n-1)/2)/(7*10^(-6)+u_H(m,(n-1)/2))-0.02)*x_H_f)*80*t/86400+c_H_f;
%%%%%惰性微生物计算%%%%!!!!！！！计算过程中发现降解氢气产甲烷微生物生长数据异常，故惰性微生物计算省略此部分
c_I_f=0.02*(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f)*t/86400+c_I_f;
%%%%%%%%微生物活性组分分率计算%%%%%%%%
x_Su_f=c_Su_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Bu_f=c_Bu_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Pro_f=c_Pro_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Ac_e_f=c_Ac_e_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_Ac_AM_f=c_Ac_AM_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_H_f=c_H_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
x_I_f=c_I_f/(c_Su_f+c_Bu_f+c_Pro_f+c_Ac_e_f+c_Ac_AM_f+c_I_f);
%%阳极室微生物生长模型%%%%%%%%%%%%%%%%%%%%%%%%%%%阳极室微生物生长同时更新微生物的浓度数据用作下一时间步长的循环%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c_Su=(3*c_Su*u_Su(m,1)/(0.5+u_Su(m,1))-0.02*c_Su)*t/86400+Qs*t*(L_c_Su-c_Su)*1000+c_Su;
c_Bu=(1.2*c_Bu*u_Bu(m,1)/(0.2+u_Bu(m,1))-0.02*c_Bu)*t/86400+Qs*t*(L_c_Bu-c_Bu)*1000+c_Bu;
c_Pro=(0.52*c_Pro*u_Pro(m,1)/(0.1+u_Pro(m,1))-0.02*c_Pro)*t/86400+Qs*t*(L_c_Pro-c_Pro)*1000+c_Pro;
c_Ac_AM=(0.4*c_Ac_AM*u_Ac(m,1)/(0.15+u_Ac(m,1))-0.02*c_Ac_AM)*t/86400+Qs*t*(L_c_Ac_AM-c_Ac_AM)*1000+c_Ac_AM;
c_H=(2.1*c_H*u_H(m,1)/(7*10^(-6)+u_H(m,1))-0.02*c_H)*t/86400+Qs*t*(L_c_H-c_H)*1000+c_H;
c_I=0.02*(c_Su+c_Bu+c_Pro+c_Ac_AM)*t/86400+Qs*t*(L_c_I-c_I)*1000+c_I;
%%%%储槽微生物生长模型%%%%%%同上%%%%%%%%%%%%%
L_c_Su=(3*L_c_Su*L_Su/(0.5+L_Su)-0.02*L_c_Su)*t/86400+Qs*t*(c_Su-L_c_Su)*500+L_c_Su;%%%%%%%积累等于进入减去输出%%%%%%Qs*t*(c_Su-L_c_Su)*500除以体积0.002m3
L_c_Bu=(1.2*L_c_Bu*L_Bu/(0.2+L_Bu)-0.02*L_c_Bu)*t/86400+Qs*t*(c_Bu-L_c_Bu)*500+L_c_Bu;
L_c_Pro=(0.52*L_c_Pro*L_Pro/(0.1+L_Pro)-0.02*L_c_Pro)*t/86400+Qs*t*(c_Pro-L_c_Pro)*500+L_c_Pro;
L_c_Ac_AM=(0.4*L_c_Ac_AM*L_Ac/(0.15+L_Ac)-0.02*L_c_Ac_AM)*t/86400+Qs*t*(c_Ac_AM-L_c_Ac_AM)*500+L_c_Ac_AM;
L_c_H=(2.1*L_c_H*L_H/(7*10^(-6)+L_H)-0.02*L_c_H)*t/86400+Qs*t*(c_H-L_c_H)*500+L_c_H;
L_c_I=0.02*(L_c_I+L_c_Bu+L_c_Pro+L_c_Ac_AM)*t/86400+Qs*t*(c_I-L_c_I)*500+L_c_I;
k=k+1;
toc
end
c;
x=1:n;y=1:m;
figure(1)
subplot(1,7,1)
plot(u_Su(m,:),'r-')
title('Su变化')
subplot(1,7,2)
plot(u_Bu(m,:),'b--')
title('Bu变化')
subplot(1,7,3)
plot(u_Pro(m,:),'g-')
title('Pro变化')
subplot(1,7,4)
plot(u_Ac(m,:),'r-')
title('Ac变化')
subplot(1,7,5)
plot(u_H(m,1:n),'b--')
title('H变化')
subplot(1,7,6)
plot(u_CH4(m,:),'g-')
title('CH4变化')
subplot(1,7,7)
plot(u_v(m,:),'g-')
title('沿生物膜厚度方向电势变化')
接下来无介体模型的Word版和含电子介体的模型Word和MATLAB程序。。
