MATLAB: Ode15s PDE equation system strange solution

method of linesode15s

Dear matlabbers,
I've been stuck with a problem the last few days, and I do not know if I should either re-think completely the problem I want to solve, or it is maybe a bug on the ode15s solver. I think it is the first one :). I have a system of two equations in partial derivatives: the heat equation and the mass diffusion equation that I'm trying to solve using the method of lines. This means that the variable I'm trying to find (temperature or concentration) depends on time and space. In other words, the solution consists on the evolution of the variable with time for each of the nodes I have divided the space into. The system consists on only one dimension. I need to solve both of them with the same call to the ode solver since in the original problem, there is an extra term that couples both equations. Unfortunately, ode15s is only able to give me a consistent solution (but not accurate) for one of the variables: the one that is called first. The solution I get for the other is an oscillating thing that does not make any sense. I do not think it is a bug on my code, since I can solve each of the equations separately with the same scripts, if I place comments on the lines that refer to the other variable.
The script that calls the function is:
clear all
close all
R = 8.31; 
MH2O = 18; 
Mair = 29.97; 
Cp_food = 2273; 
rho_food = 1299; 
k_food = 0.34; 
alpha_food = k_food / (rho_food * Cp_food); 
Tair = 200; 
ms = 0.95; 
Yinf = 0.08 * MH2O / Mair;
c0 = 0.05 * ones(50,1); %kg H2O kg d.s.-1
T0 = 25 * ones(50,1); 
y0 = [T0; c0];
cair = Yinf;
hav = 60; %W m-2 K-1
DHevap = 2264.76 * 1000;
L_food = 0.03; 
t = [0 3600];
D   = 1e-6;      % Diffusivity (m2/s)
kmav = 0.08;    % mass transfer coefficient (m/s)
nx=length(c0);
delta_x = L_food/(nx-1);
pos = [1:1:nx];
save parameters
options=odeset('RelTol',1e-9,'AbsTol',1e-9);
[t,y]=ode15s(@coupled_transfer,t,y0,options);
The function containing the discretized differential equations is:
function dydt = coupled_transfer(t,y)
load parameters
% size(y)
T = y(1:nx);
c = y(nx:end);
dTdx2 (1) = (T(3) - T(1))/delta_x^2; %probably the same as dTdx2 (1) = 0;
dCdx2 (1) = (c(3) - c(1))/delta_x^2; %probably the same as dCdx2 (1) = 0;
for j = 2:nx-1;
dTdx2 (j) = (T(j+1) + T(j-1) - 2*T(j))/delta_x^2;
dCdx2 (j) = (c(j+1) + c(j-1) - 2*c(j))/delta_x^2;
end
dTdx2(nx) =  - hav * (T(nx)-Tair)/ k_food;
dCdx2(nx) =  - kmav   * (c(nx)-cair)/D;
dCdt = D * dCdx2';
dTdt = alpha_food * dTdx2';
dydt = [dTdt; dCdt];
I suspect that the issue is due to the different times required by each of the equations to reach equilibrium. Actually, the equation for the concentration reaches equilibrium in every node much sooner that the one for the temperature. I admit that how the solvers work is fairly unknown to me. You can find attached what the solutions look like for each variable. I get these two warnings (I might not get them if instead of one hour I simulate 20 s).
    Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
    > In ode15s (line 589)
      In parameters (line 51) 
    Warning: Failure at t=1.891510e+02.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.547474e-13) at time t. 
    > In ode15s (line 668)
      In parameters (line 51)
I hope I have exposed my problem in a clear way.
Your comments are very appreciated.

Best Answer

c=y(nx+1:end);
Best wishes
Torsten.
Related Question