Solving a system a coupled ODE and PDE

10 次查看（过去 30 天）
Quentin Carboué
Quentin Carboué 2021-2-18
评论： Quentin Carboué 2021-2-20
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?

请先登录，再回答此问题。

回答（1 个）

Bill Greene
Bill Greene 2021-2-19
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.
  1 个评论
Quentin Carboué
Quentin Carboué 2021-2-20
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by