How to solve a coupled ODE, PDE system with method of lines?

Hi,

I am trying to solve a problem of gas flow through a packed catalyst bed with reaction (gas phase species being adsorbed by the solid catalyst surface). I have the transport equation dCdt = D_ax*d2Cdz2 - u*dCdz + f. I have neglected the dispersion term (D_ax*d2Cdz2), and so there is no second derviative w.r.t. space. The source term 'f' depends on the concentration of the gas C, as well as on the fraction of available sites on the solid catalyst particles (1-q). The gas crosses the control volume boundaries via convection, i.e. has a derivative w.r.t. space however the available adsorption sites are fixed in position, and therefore even though they have a spatial distribution, from the mass balance it follows that theta only has a derivative w.r.t. time. This means I have coupled ODE (balance over the adsorption sites) with a PDE (balance over the gas species):

(1) dC/dt = -u*dC/dz - R

(2) dq/dt = R

Where: R = k*C*(1-q)

The I.C.'s are: (1) C(z,t=0) = 0, and (2) q(z,t=0) = 0. The B.C. used for (1) is a time dependent interpolated concentration step CA0. (2) Does not require a B.C.

I keep getting an error: "Error using odearguments. MYFUNC must return a column vector".... I hope somebody could help me? I have tried reading many other similae questions and answers which have got me to this point, but I really have no idea how to solve this... Thanks a lot!!

This is my script:

Nz = 10; %[-]

dz = L_bed/Nz;

k_opt = 1e6;

CA0 = inlet_ion_current./max(inlet_ion_current) .* max(CO2_alpha); %Defining input B.C. CA0(t) based on data

t_interp_cap = [0:0.5:length(CA0)*2]'; %Time steps belonging to input B.C. CA0(t)

tspan = [0:0.01:t_interp_cap(end)]';

%Initial conditions

C = zeros(Nz,1);

q = zeros(Nz,1);

y_IC = [C;q];

%Solver

[t, SOL] = ode15s(@myfunc,tspan,y_IC,[],Nz,CA0,t_interp_cap,dz,u,k_opt);

%ODE function definition

function dydt = myfunc(t,y,Nz,CA0,CA0_t,dz,u,k1)

global m_cat rho_b A n_CO2_max %(IGNORE) calling parameters

CA0 = interp1(CA0_t,CA0,t,'pchip'); %Interpolating B.C. CA0 to evaluate at tspan times

dydt = zeros(Nz*2,1); %preallocation of space, Nz*2 ODEs (Nz ODEs fror C and NZ for q)

C(1:Nz) = y(1:Nz);

C(1) = CA0; %C(t,z=0) = CA0(t)

q(1:Nz) = y(Nz+1:2*Nz);

for i = 2:Nz

R(i) = k1.*C(i).*(1-q(i)); %Source term: R = f(C,q)

dCdz(i) = 1./(2.*dz).*(C(i)-C(i-1)); %Finite element difference (backward) only calculated for C (dC/dz)

dCdt(i) = -u.*dCdz(i) - R(i).*(n_CO2_max*A*rho_b/m_cat); %Equation (1)

dqdt(i) = R(i); %Equation (2)

end

dydt = [dCdt; dqdt];

end

[Torsten]

Initialize dCdt and dqdt as

dCdt = zeros(Nz,1);

dqdt = zeros(Nz,1);

Further, u seems to be undefined.

dz = L_bed/(Nz-1) instead of dz = L_bed/Nz

Ca0 and Ca0_t have the same length ?

In the definition of dCdz, you have to divide by deltaz, not 2*deltaz. First order upwind for dCdz requires many more points than 10 - you will have to choose a much larger value for Nz to capture the concentration gradient.

Maybe you want to calculate q also in the first grid point ? Because your loop starts at i=2 ...

Maybe there are some more inconsistencies - I can't run the code where I am at the moment.

Answer 1

[Image Analyst]

So what happens if you make it a column vector like it wants by adding this line to the end of your function

dydt = dydt(:); % or dydt = reshape(dydt, [], 1);