Solving system of ODE and PDE in 1D

I am trying to solve a system of ODE and PDE (3 ODEs and 1PDE). The governing equations, initial and boundary conditions used are listed in the attached Model_Eqns - Copy.pdf file. Here is the code so far I have written. Apparently, I can't find any bugs but if I run the code I encounter numerical errors. I am not sure where is the problem. Can anyone please help me?

% Discretization parameters
L = 4.0; % length of tube
N=2000; % no. of discretization in X direction
X = linspace(0,L,N); % discretise the domain 
n = numel(X);

% time interval
ti=0;
tend=300;   % end time 
array_length=2*tend;
tspan = linspace (ti,tend,array_length); % given time span
%-------------------------------------------------------------

% Initial conditions for ODE
Tgas0      = 293.0*ones(n,1);        
Tbed0     = 293.0     *ones(n,1);     
Tw0     = 293.0 *ones(n,1);     
Q0      = 2.5*ones(n,1);     
% X0      = 0     *ones(n,1);     % Conversion rate
%-------------------------------------------------------------

% setting up the solver
y0 = [Tgas0; Tbed0; Tw0; Q0];
%  y0 = [Tgas0; Q0];
[T,Y] = ode23(@(t,y) solution(t,y,X,n),tspan,y0);
Y=Y.';
%-------------------------------------------------------------

% saving results
Tgas = Y(1  :  n,:);
Tbed = Y(n+1:2*n,:);
Tw = Y(2*n+1:3*n,:);
Q  = Y(3*n+1:4*n,:);
% X  = Y(4*n+1:5*n,:);
%-------------------------------------------------------------

figure
subplot(2, 1, 1);
%plot (T, Tw(N,1:tend), 'g', T, Tbed (N*0.5,1:tend), 'red',   T, Tgas (N*0.5,1:tend),'b');
plot ( T,  Tw (N,1:600), 'red');
%plot ( T,  Tbed (N,1:600), 'red');
% plot ( T,  Tgas (1:600), 'red');
xlabel('time')
ylabel('Temperature')   % how to define water temperature at the outlet 
legend('Water at outlet');


subplot(2, 1, 2);
plot(T, Q (N, 1:600));
title('Loading');
xlabel('Time');
ylabel('Loading at middle [mol/kg_m_o_f]');
%----------------------------------------------------------------


function DyDt = solution(~,y,X,n)

% Define Constants 

%water
V_w=0.197635; %velocity of water m/s
m_w= 0.156938; % mass of water kg
Cp_w= 4184; %J/kg.K

% Bed 
% alpha_bed = 3.33942E-6; % m2/s thermal diffusivity 
H_ads  = 26000;   % J/mol Heat of sublimation of CO2
m_bed= 0.15; %kg =mass of MOF 
Cp_bed= 620; % J/kg.K Cp of MIL 101
h_w= 600; % HTC bed to water W/m2.K fitted
A_tube = 0.091985755; % m2 tube surface area

% Gas 
h_gas= 25; %W/m2.K HTC bed to gas -fitted
A_fin= 0.606602; %m2 
m_gas=0.007286; %kg 
Cp_gas= 846; % J/kg.K
m_MOF=0.15; %kg 
Cv_gas= 28.93; % J/mol.K


% Adsorption parameters 
K_LDF  = 6.3211E-2; % s-1
Q_eq   = 30;    % mol/kg


% Define matrix

%Change in variables with time
DTgasDt = zeros(n,1); 
DTbedDt = zeros(n,1); 
DTwDt = zeros(n,1); 
DQDt = zeros(n,1); 


% from saving results
Tgas = y(1:n);            
Tbed = y(n+1:2*n);        
Tw = y(2*n+1:3*n);      
Q  = y(3*n+1:4*n);  

%     X  = y(4*n+1:5*n);      % Conversion rate


% Constants for Governing equations 
GC_A=(h_gas*A_fin);  % GC=Gas Constant
GC_B=(m_gas*Cp_gas);
GC1= GC_A/GC_B; 
GC_2= (m_MOF*K_LDF*Cv_gas)/GC_B;

BC_A=(H_ads*m_MOF*K_LDF);  % BC= Bed Constants 
BC_B=(m_bed*Cp_bed);
BC_C=(h_w*A_tube);
BC1= BC_A/BC_B; 
BC2= GC_A/BC_B; 
BC3=BC_C/BC_B; 

WC_A=(m_w*Cp_w); % WC= water constants 
WC1=BC_C/WC_A; 


% Governing equations 
% Tgas, Tbed and Q  => ODE
for i = 1:n
  DTgasDt(i) = (GC1*(Tbed(i)-Tgas(i)))- (GC_2*(Q_eq - Q(i))*Tgas(i)) ;
  
  DTbedDt(i)= (BC1*(Q_eq - Q(i))) + (BC2*(Tgas(i)-Tbed(i))) +(BC3* (Tw(i)-Tbed(i))); 
  
  DQDt(i) = K_LDF*(Q_eq - Q(i));

end

%water
% BC at x = 0
DTwDt(1) = 293.0;   

% interior and last grid point
for i = 2:n
  DTwDt(i) = V_w* ((Tw(i)-Tw(i-1))/(X(i)-X(i-1))) + (WC1*(Tbed(i) - Tw(i))) ;     
end

%--------------------------------------------------------------------------

DyDt = [DTgasDt; DTbedDt;DTwDt; DQDt]; 

end

Answer 1

[Torsten]

Except for

DTwDt(1) = 293.0;

which must of couse read

DTwDt(1) = 0.0;

it seems the equations as written are correctly implemented.

So my guess is that there must be something wrong with the equations themselves.

[JUBAIR AHMED SHAMIM]

Hi Torsten. Thanks for your input again. Actually the correct expression for inlet BC would be

Tw(1)=293.0; [not DTwDt] as our inlet water temperature is 20 deg C. If I supply BC using DTwDt, it is not supplied at all the time steps at X=1.

However my conclusion as follows: I tried to use very small CFL numbers using small time steps and del X. Also tried stiff ODE solvers like ODE15s. But results are same. After some time steps solution becomes unstable.

I think Matlab ODE solvers not capable of handling discretized equations of our system. More robust numerical solvers may be necessary like gauss Seidel, Jacobi or TDMA.

[JUBAIR AHMED SHAMIM]

Hi Torsten, I have found the error. Vw in my equation for water should be negative. Now it woks.

DTwDt(i) =  (WC1*(Tbed(i) - Tw(i))) - V_w*((Tw(i)-Tw(i-1))/(X(i)-X(i-1));

and the initial part I have modified slightly to make CFL smaller and using ODE15s for the stiff equations:

% Discretization parameters
L = 3.7; % length of tube
step=0.01;
N=L/step; % no. of discretization in X direction
X = linspace(0,L,N); % discretise the domain 
n = numel(X);

% time interval
ti=0;
tend=300;
delt=0.001;% end time 
Nt=tend/delt;
tspan = linspace (ti,tend,Nt); % given time span
%-------------------------------------------------------------

CFL= 0.197*(delt/step); % water flow velocity 0.197 m/s
display (CFL)

and BC for water at inlet should be supplied as:

Tw(1) = 293.0;   

[Torsten]

Your initial condition for Tw was Tw(1) = 293 K. If you set DTwDt(1) = 0, Tw(1) will remain 293 K throughout the simulation - no need to reset Tw(1).

The number of output times you choose is irrelevant for the integration - they don't fix the internal stepsize of the ode integrator. Thus choosing delt = 0.1 or delt = 0.001 will only change the number of output times, but doesn't have an effect that CFL is better satisfied.

[JUBAIR AHMED SHAMIM]

Thanks for the comment. I agree choosing output time step may not have effect on the integration. But if I don't choose this, the solver outputs different step size at different times if I change something in the code. So plotting of results becomes difficult. Selecting a specific time step at least this problem of plotting results can be solved.

[Torsten]

What I meant to say was: you can of course prescribe output times, but the vector you prescribe should not influence the result you get.

Thus definining 600 output times or 60 output times will influence the length of T and Y, but not the solution curves themselves.