Solving a system a coupled ODE and PDE

In order to study the evolution of the temperature and the moisture content inside an aerated bioreactor I want to solve a system of 4 equations (2 ODE and 2 PDE) that model the mass and energy transfers (see attached file).

For the PDEs I used a finite difference method in order to discretize the spatial derivative on N points and thus obtaining a system of N ODEs.

I then multiply the ODE equation with a vector of size N to have N ODEs (same size as the PDEs).

The function I generated for the system is the following:

function [dydt] = ode_syst2(t, Tg, Ts, Fig, Fis)

%Parameters:
Cpg = 1005; %(J/kg of dry air*°C) Heat capacity of dry air
Cps = 1455; %(J/kg of dry solid*°C) Heat capacity of the solid
Cpv = 1791; %(J/kg of vapor*°C) Heat capacity of the water vapor
Cpw = 4184; %(J/kg of water*°C) Heat capacity of the liquid water
he = 16000; %(J/(s*m^3 of bed*°C) Solid heat transfer coefficient
G = 0.114; %(kg of dry air/(m^2*s)) Air flow rate through the bed
L = 0.4; %(m) Overall bed height
ke = 0.06; %(kg water/s*m^3 of bed)
P = 91325; %(Pa) Total pressure within the system
epsilon = 0.561; %(m^3 of void/m^3 of bed) Porosity of the bed
lambda = 2414300; %(J/kg water) Enthalpy of vaporization of water
rog = 1.14; %(kg dry air/m^3) Density of the air phase
S = 190; %(kg solid DM/m^3) Bulk density of the bed (Dry basis) (For all z, S0 = 190 kg/m^3 solid (DM))
Pwsatin = (133.322*(exp(18.3036-(3816.44/(Tgin+273.15-46.13))))); %Saturation pressure of water of the air at the ir inlet
Figin = (0.62413*(Pwsatin/(P-Pwsatin))); % (kg of water/kg of dry solid) Humidity of the air phase at the air inlet
Tgin = 30; %Temperature at the air inlet

N = 21;
dx = (L-1)/(N-1);
vec = ones(1,N); %spatial vector representing the lenght of the bioreactor vec(1) = entrance, vec (N) = exit

%Energy balance over the gas phase
dTgdx = ones(1,N);
dTgdx (1) = Tgin; %Boundary condition
for x = 2:N
    dTgdx (x) = ((-(he*(G/0.095))*((Tg*t)-(Ts*t))-(G*(Cpg+Cpv*Fig(t))))*((Tg*(x) - Tg*(x-1))/(dx)))/(epsilon*rog(Cpg+(Cpv*(Fig*t))));
end
dTgdt = 1*dTgdx';

%Energy balance over the solid phase
%All the values of Ts are equal along the lenght of the bioractor at a given time t
dTsdx = vec*(((he*(G/0.095))*((Tg*t)-(Ts*t)))-((ke*(G/0.095)*(((Fis*t)/1.5)))*lambda*((Fis*t)-((0.095*(((((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))/(1-(((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))))^0.488))))))/(S(Cps+(Cpw*(Fis*t))));
dTsdt = 1*dTsdx';

%Mass balance over the gas phase
dFigdx = ones(1,N);
dFigdx (1) = Figin; %Boundary condition
for x = 2:N
    dFigdx (x) = ((-G*(((Fig*x) - Fig*(x-1))/(dx)))+(((ke*(G/0.095)*(((Fis*t)/1.5)))*((Fis*t)-((0.095*(((((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))/(1-(((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))))^0.488))))))/(epsilon*rog));
end
dFigdt = 1*dFigdx';

%Mass balance over the solid phase
%All the values of Fis are equal along the lenght of the bioractor at a given time t
dFisdx = vec*((-(ke*(G/0.095)*(((Fis*t)/1.5))))*((Fis*t)-((0.095*(((((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))/(1-(((Fig*t)*P)/((Fig*t)*(((133.322*(exp(18.3036-(3816.44/((Tg*t)+273.15-46.13))))))+0.62413)))))^0.488))))/S);
dFisdt = 1*dFisdx';


dydt = [dTgdt; dTsdt; dFigdt; dFisdt];
end

Afer that, I tried to solve the system using the following code. We create vector of size N for the 4 ODEs that contain each N equations.

%Initial condition
Pwsatin = (133.322*(exp(18.3036-(3816.44/(Tgin+273.15-46.13))))); %Saturation pressure of water of the air at the ir inlet
Figin = (0.62413*(Pwsatin/(P-Pwsatin))); % (kg of water/kg of dry solid) Humidity of the air phase at the air inlet
Tgin = 30; %Temperature at the air inlet
N=21;

Tg0 = Tgin * ones(N,1);
Ts0 = Tgin * ones(N,1);
Fig0 = Figin * ones(N,1);
Fis0 = 1.5 * ones(N,1);


IC=[Tg0, Ts0, Fig0, Fis0];

t = [0 21];

[t,dydt] = ode15s(@ode_syst2, t, IC);

I get an error message stating that I don't have enough input arguments.

Can you please help me with this one?

Answer 1

[Bill Greene]

The error message is caused by your definition of the ode function:

function [dydt] = ode_syst2(t, Tg, Ts, Fig, Fis)

The definiton should be of the form:

function [dydt] = ode_syst2(t, y)

You have to extract the dependent variable vectors, Tg, etc from the single vector, y.

[Quentin Carboué]

Dear Bill, thank you for your reply, but how can I create a single vector y containing all my dependent variables?

I tried to assign my dependent variablesas following but it still doesn't work.

Tg = y(1:N);
Ts = y(N+1:2*N);
Fig = y(2*N+1:3*N);
Ts = y(3*N+1:4*N);