# Solving a system a coupled ODE and PDE

38 views (last 30 days)
Quentin Carboué on 18 Feb 2021
Commented: Quentin Carboué on 20 Feb 2021
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.

Bill Greene on 19 Feb 2021
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é on 20 Feb 2021
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);

### Categories

Find more on Electromagnetics in Help Center and File Exchange

R2020a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!