recdiffq9.m

反应扩散问题长期在研究领域中c++去实现, 现在我们用MATLAB也可以

function rec_diff_d2q9
omega=0.5063; 
%定义网格数目
nx=256; 
ny=256;
%反应的物性常数
n=2;
FLD=96487;
A=270;                 %单位cm2
R=8.3144;
N=60000;
k01=3.48e+15;
k0_1=1.41e+15;
Ea1=4e+5;
Ea2=4.066e+5;
S298=-44.397;
syms T1;
T=input('input the temperature : ')

CPH2O=(30.3794+9.6212*10^(-3)*T1+1.1848*10^(-6)*T1^2)/T1;
SRH2O=int(CPH2O,298,T);
SRH2O=eval(SRH2O);

CPH2=(29.0856-0.8373*10^(-3)*T1+2.0138*10^(-6)*T1^2)/T1;
SRH2=int(CPH2,298,T);
SRH2=eval(SRH2);

CPO2=(25.8911+12.9874*10^(-3)*T1-3.8644*10^(-6)*T1^2)/T1;
SRO2=int(CPO2,298,T);
SRO2=eval(SRO2);

S=S298+SRH2O-SRH2-0.5*SRO2
n_T=1;
ka=k01*exp(-Ea1/R/T);
kd=k0_1*exp(-Ea2/R/T);
po2=9.9167;                 %单位mol/s
ph2o=0;
c=1/9;      

den=15.2250;                %单位mol/s
%初始条件随机正态分布
iniden=normrnd(den,0.1,nx,ny);
for k=1:9
    for x=1:nx
        for y=1:ny
            F(x,y,k)=iniden(x,y)/9;
            DENSITY(x,y)=iniden(x,y);
        end
    end
end
FEQ=F;
ts=0;
%step=10000;
step=input('input the step : ')
DENSITYave=repmat(0,[1 step]);
denj=repmat(den,[1 step]);
i=1;
ii=1;
while(i<=step) 
    while (ts<step)
    % Propagate
    F(:,:,4)=F([2:nx 1],[ny 1:ny-1],4);
    F(:,:,3)=F(:,[ny 1:ny-1],3);
    F(:,:,2)=F([nx 1:nx-1],[ny 1:ny-1],2);
    F(:,:,5)=F([2:nx 1],:,5);
    F(:,:,1)=F([nx 1:nx-1],:,1);
    F(:,:,6)=F([2:nx 1],[2:ny 1],6);
    F(:,:,7)=F(:,[2:ny 1],7); 
    F(:,:,8)=F([nx 1:nx-1],[2:ny 1],8);
  %  DENSITY=sum(F,3);
    for ic=1:9;
            FEQ(:,:,ic)=c*DENSITY;
    end
    F=omega*FEQ+(1-omega)*F;
    force=ka*DENSITY*sqrt(po2)-kd*ph2o;
    po2=po2-0.0010;                           %1073k
    ph2o=ph2o+0.0020;                         %1073k
 %   po2=po2-0.0005;                           %1173k
 %   ph2o=ph2o+0.0010;                         %1173k
    forceave=-force/9;
    if forceave>=0
  %      break
       ii=ii+1;
       forceave=0;
    end
    for ic=1:9;
            F(:,:,ic)= F(:,:,ic) + forceave;
    end
    DENSITY=sum(F,3);
    B=reshape(DENSITY,nx*ny,1);
    Bave=mean(B);
    DENSITYave(1,i)=Bave;
    j(1,i)=(denj(1,i)-DENSITYave(1,i))*n*FLD/A/N;
    i=i+1;
    ts=ts+1;
    end 
end
rec_t(1,n_T)=1/kd;
rec_equ_t(1,n_T)=rec_t(1,n_T)*(step-ii)/step;
j_T(1,n_T)=j(1,step);
DENSITYave_T(1,n_T)=DENSITYave(1,step);
n_T=n_T+1;
%qq=j(1,step)
q=T*S*j(1,step)/n/FLD

figure;                       
[X,Y]=meshgrid(50/256:50/256:50,0.01/256:0.01/256:0.01);mesh(X,Y,DENSITY(1:nx,1:ny))
title(['Concentration field after ',num2str(rec_t),'\deltat']);
hold on;
i=[0:rec_t/step:rec_t-rec_t/step];
figure;
plot(i,DENSITYave);
grid on;
title('concentration distribution of hydrogen');
xlabel('time(s)');
ylabel('concentration(mol)');
figure;
plot(i,j);
grid on;
title('current density distribution');
xlabel('time(s)');
ylabel('current density(A/cm2)');