Method of Lines multiple PDE system

I have solved a coupled PDE system of two simultaneous PDEs, I am trying to solve a further system of four simultaneous PDEs. The code below shows the solution for a two PDE system but I can't figure out a way of solving it with two further equations.

% Clear previous files clear all clc

% Parameters shared with other routines global ncall ncase n z ndss

% Two cases : ncase = 1, second order finite differences % ncase = 2, fourth order finite differences

for ncase = 1:2

    % Initial condition 
    n = 201; % Number of lines in the space grid 
    t0 = 0.0;
    y0 = initial_1(t0);

    % Independent variable for ODE integration
    tf = 1000.0; % Simulation time
    nout = 101; % Number of lines in the time grid
    tout = linspace(t0,tf,nout); % Vector of time grid
    ncall = 0;

    % ODE integration ;
    reltol=1.0e-8; abstol=1.0e-8; % Tolerance level
    options=odeset('RelTol',reltol,'AbsTol',abstol);
    if(ncase==1), ndss=2; % ndss = 2, 4, 6, 8 or 10 required
    [t,y]=ode15s(@pde_1,tout,y0,options); end
    if(ncase==2), ndss=4; % ndss = 2, 4, 6, 8 or 10 required
    [t,y]=ode15s(@pde_1,tout,y0,options); end

% One vector to two vectors

Ca = zeros(1,length(n)); % Preallocating Ca for faster processing Xb = zeros(1,length(n)); % Preallocating Xb for faster processing

for it=1:nout for i=1:n

Ca(it,i) = real(y(it,i)); Xb(it,i) = real(y(it,i+n));

end end

% Plot solutions if(ncase==1)

% Parametric plots figure(1); plot(z,Ca(nout,:),'-r') hold on plot(z,Xb(nout,:),'-b') hold on axis([0 10 0 1.01]) xlabel('z [m]') ylabel('CA/CAO or XB')

% Surface plots figure(2); surf(real(Ca)); axis tight xlabel('z grid number'); ylabel('t grid number'); zlabel('Ca'); title('Concetration of Fluid A'); colormap summer rotate3d on;

figure(3); surf(real(Xb)); axis tight xlabel('z grid number'); ylabel('t grid number'); zlabel('Xb'); title('Conversion of Solid B'); colormap winter rotate3d on;

end end

function yt=pde_1(~,y) %(t,y)

%Parameters of fluid A and solid B in the reactor a = 2; %Stoichiometric coefficient of fluid A b = 1; %Stoichiometric coefficient of solid B rho_b = 2960; %Density of solid B in kg/m^3 Rs = 1*10^-4; %Initial radius of solid B in m Mb = 84.3; %Molecular mass of solid B e = 0.4; %Void fraction u = 0.25; %Superficial velocity in m/s kt = 1*10^-4; %Fluid to solid mass transfer coefficient in m/s De = 1*10^-9; %Effective diffusion coefficient of fluid A in m/s k = 8*10^-6; %First order surface rate constant in m/s Np = (3*(1 - e))/(4*pi*(Rs^3)); %Number of particles suspended in a unit volume of reactor

% Function pde_1 defines the ODEs in the MOL solution of two nonlinear PDEs

global ncall ncase Ca Xb Cat Xbt Caz Cazz Xbz Xbzz n z

% One vector to two vectors for i=1:n

Ca(i)=y(i); Xb(i)=y(i+n);

end

% Boundary conditions at x = 0 Ca(1)=1.0; Cat(1)=0.0; Xb(1)=0.0; Xbt(1)=0.0;

% First order spatial derivatives zl=z(1); zu=z(n);

% Three point centered differences if(ncase==1), Caz=dss002(zl,zu,n,Ca); end if(ncase==1), Xbz=dss002(zl,zu,n,Xb); end

% Five point centered differences if(ncase==2), Caz=dss004(zl,zu,n,Ca); end if(ncase==2), Xbz=dss004(zl,zu,n,Xb); end

% Boundary conditions at z = 10 Caz(n)=0; Xbz(n)=0;

% Second order spatial derivatives % Three point centered differences if(ncase==1), Xbzz=dss002(zl,zu,n,Xbz); end

% Five point centered differences if(ncase==2), Xbzz=dss004(zl,zu,n,Xbz); end

% Array Caz is used as temporary storage in the calculation % of the term ((Xb - 1)*Ca ) which is finally stored in array Cazz z z for i=1:n Xbz(i)=(real(Xb(i))-1.0)*real(Caz(i)); end if(ncase==1), Cazz=dss002(zl,zu,n,Xbz); end if(ncase==2), Cazz=dss004(zl,zu,n,Xbz); end

% Two PDEs for i=2:n

kdash = (1/((1/kt) + ((Rs/De)*(((1 - real(Xb(i)))^(-1/3)) - 1)) + (((1 - real(Xb(i)))^(-2/3))/k))); Xbt(i) = (((3*Mb*b)/(Rs*rho_b*a))*kdash*real(real(Ca(i)))); Cat(i) = ((-4*pi*(Rs^2)*kdash*real(Ca(i))*Np) - u*real(Caz(i)))/e;

end

% Two vectors to one vector

yt = zeros(1,length(n));

for i=1:n

yt(i) = real(Cat(i)); yt(i+n) = real(Xbt(i));

end

yt=yt';

% Increment calls to pde_1 ncall = ncall+1;

function y=initial_1(~) % (t0)

% Function inital_1 is called by the main program to define % the initial conditions in the MOL solution of two nonlinear % PDEs

global ncall Ca Xb n z

% Spatial increment dz = 10.0/(n-1);

% Values of x along the spatial grid z = 0.0:dz:10.0;

% Initial conditions y = zeros(1,length(n));

for i=1:n

Ca(i) = 1.0; Xb(i) = 0.0; y(i) = Ca(i); y(i+n) = Xb(i);

end

% Initialize calls to pde_1 ncall = 0;

function [Caz]=dss002(zl,zu,n,Ca)

% Compute the spatial increment dz=(zu-zl)/(n-1); r2fdz=1./(2.*dz); nm1=n-1;

% Equation (1) (note - the rhs of the finite difference approxi- % tions, equations (1), (2) and (3) have been formatted so that % the numerical weighting coefficients can be more easily associ- % ated with the Bickley matrix listed above) Caz(1)=r2fdz*(-3.*real(Ca(1))+4.*real(Ca(2))-1.*real(Ca(3)));

% Equation (2) for i=2:nm1

     Caz(i)=r2fdz*(-1.*real(Ca(i-1))+0.*real(Ca(i))+1.*real(Ca(i+1)));

   end

% Equation (3) Caz(n)=r2fdz*(1.*real(Ca(n-2))-4.*real(Ca(n-1))+3.*real(Ca(n)));

 function [Caz]=dss004(zl,zu,n,Ca)

% Compute the spatial increment dz=(zu-zl)/(n-1); r4fdz=1./(12.*dz); nm2=n-2;

% Equation (1) (note - the rhs of equations (1), (2), (3), (4) % and (5) have been formatted so that the numerical weighting % coefficients can be more easily associated with the Bickley % matrix above) Caz(1)=r4fdz*(-25.*real(Ca(1)) +48.*real(Ca(2))-36.*real(Ca(3))+16.*real(Ca(4))-3.*real(Ca(5)));

% Equation (2) Caz(2)=r4fdz*(-3.*real(Ca(1)) -10.*real(Ca(2)) +18.*real(Ca(3))-6.*real(Ca(4))+1.*real(Ca(5)));

% Equation (3) for i=3:nm2

     Caz(i)=r4fdz*(+1.*real(Ca(i-2)) -8.*real(Ca(i-1)) +0.*real(Ca(i))+8.*real(Ca(i+1))-1.*real(Ca(i+2)));

   end

% Equation (4) Caz(n-1)=r4fdz*(-1.*real(Ca(n-4)) +6.*real(Ca(n-3)) -18.*real(Ca(n-2)) +10.*real(Ca(n-1)) +3.*real(Ca(n)));

% Equation (5) Caz(n)=r4fdz*(3.*real(Ca(n-4)) -16.*real(Ca(n-3)) +36.*real(Ca(n-2)) -48.*real(Ca(n-1)) +25.*real(Ca(n)));

Thanks.

Answer 1

[Albert Shikongo]

email me at [e-mail]