# In-Situ Combustion Model - MATLAB code

9 views (last 30 days)
I need to solve in MATLAB the following system:
[![][1]][1]
And for it I would like to use these datas and scheme (Crank-Nicolson)
With this initial condition:
I really do not know how I can do this on MATLAB. I tried to put on MATLAB all my datas, as you can see but I do not know as I can finish.
About this code I have got the following error. This error arise in all the pharenteses of the scheme:
. But I don't know exactly where to put these and it seems that these are not about the error that I get.
Many thanks for any help!
tic
clear all
close all
%% INPUT PARAMETERS
Xmax=0.05; % x maximal value
Xmin=0; % x maximal value
tmax=0.008; % x maximal value
tmin=0; % x maximal value
N=101 % number of iterations
dt=10^(-5); % time step value
tol=10^(-6); % precison on the staionnary solution
Pet=1406;
beta=7.44*10^(10);
epsilon=93.8;
theta0=3.67;
u=3.76;
phi=beta*(1-eta)*e^(-epsilon/(theta+theta0));
rho=theta0/(theta+theta0);
V0=[theta;eta];
% discretise the domain
dx= (xmax-xmin)/N;
%%
% Initial conditions
theta(0,t)=0;
eta(0,t)=1;
theta(x,0)=0;
eta(x,0)=0;
%%
V=V0;
VCN=V
%%
% loop through time
nsteps= tmax/dt;
for n=1 : nsteps
% calculate the scheme
for i = 1 : N+1
for j=1 : N+1
theta(x(j),t(i+1))+u*dt/(4*dx)(eta(x(j+1),t(i+1))*theta(x(j+1),t(i+1))-eta(x(j-1),t(i+1))*theta(x(j-1),t(i+1)))-
dt/(2*Pet*dx)*(theta(x(j-1),t(i+1))-theta(x(j),t(i+1))+theta(x(j+1),t(i+1)))-dt/2*phi(x(j),t(i+1))==theta(x(j),t(i))-
u*dt/(4*dx)*(eta(x(j+1),t(i))*theta(x(j+1),t(i))-eta(x(j-1),t(i))*theta(x(j-1),t(i)))+dt/(Pet*2dx)*(theta(x(j-1),t(i))-
2*theta(x(j),t(i))+theta(x(j+1),t(i)))+dt/2*phi(x(j),t(i));
dt/2*phi(x(j),t(i+1))+eta(x(j),t(i+1))==eta(x(j),t(i))-dt/2*phi(x(j),t(i));
end
end
% update t and V
t=t+dt;
V=VCN;
end
function x=x(j);
x(j)=j*dx
end
function t=t(i);
t(i)=i*dt
end
function eta=eta(x,t);
end
function theta=theta(x,t);
end
function phi=phi(x,t);
phi(x,t)=beta*(1-eta(x,t))*e^(-epsilon/(theta(x,t)+theta0))
end
% update t and V
t=t+dt;
V=VCN;
% plot solution
set(fig,'YData',V(1,:));
set(tit,'String',strcat('n =',num2str(n)));
drawnow;
end

John D'Errico on 15 Jun 2019
Edited: John D'Errico on 15 Jun 2019
This is a simple error of syntax, which seems a bit strange for someone who has written such a long code to misunderstand. You should spend a little time learning the basic syntax of MATLAB.
MATLAB uses * to multiply. Or .* if you want to make sure that this is an element-wise matiplication. But it uses an operator to multiply.
The point is though, that
(x-y)(x+y)
is NOT the product of those two numbers, at least not in MATLAB. It may be that when you write it on paper, etc. But in MATLAB? What I wrote there is a syntax error. If you want to multiply numbers, then you genrally need to write
(x-y).*(x+y)
in case they are vectors or matrices. If both x and y are scalars, then it would be sufficient to write:
(x-y)*(x+y)
But in any case, you can see that I put an operator in there, between the pieces. Now go back and read the error message. Look at the point that is indicated. MATLAB sees a situation just like what I wrote, and it does not understand what you intended there, because what you wrote is not proper MATLAB syntax. For example, when MATLAB sees this generic form:
A(B)
it assumes that A is a function, that will be called with an operand B. So, when I write
(A)(B)
it gets confused. Is A a function, and B an argument? Computers need syntax to communicate with you.

Noemi Zeraick Monteiro on 15 Jun 2019
tic
clear all
close all
%% INPUT PARAMETERS
xmax=0.05; % x maximal value
xmin=0; % x maximal value
tmax=0.008; % x maximal value
tmin=0; % x maximal value
dt=10^(-5); % time step value
tol=10^(-6); % precison on the staionnary solution
Pet=1406;
beta=7.44*10^(10);
epsilon=93.8;
theta0=3.67;
u=3.76;
phi==beta*(1-eta)*e^(-epsilon/(theta+theta0));
rho==theta0/(theta+theta0);
V0==[theta;eta];
% discretise the domain
dX=(xmax-xmin)/500;
%%
V=V0;
VCN=V
%%
% Initial conditions
for i=0: 800;
if x==0;
theta(x,t(i))==0;
eta(x,t(i))==1;
end
end
for i=0: 800;
if t==0;
theta(x(i),t)==0;
eta(x(i),t)==0;
end
end
% loop through time
nsteps= tmax/dt;
for n=1 : nsteps
% calculate the scheme
for i = 1 : 800
for j=1 : 800
theta(x(j),t(i+1))+u*dt*(1/(4*dX))*(rho(x(j+1),t(i+1))*theta(x(j+1),t(i+1))-rho(x(j-1),t(i+1))*theta(x(j-1),t(i+1)))-(dt/(2*Pet*dX*dX))*(theta(x(j-1),t(i+1))-theta(x(j),t(i+1))+theta(x(j+1),t(i+1)))-(dt/2)*phi(x(j),t(i+1))== theta(x(j),t(i))-u*dt*(1/(4*dX))*(rho(x(j+1),t(i))*theta(x(j+1),t(i))-rho(x(j-1),t(i))*theta(x(j-1),t(i)))+ (dt/(Pet*2*dX*dX))*(theta(x(j-1),t(i))-2*theta(x(j),t(i))+theta(x(j+1),t(i)))+(dt/2)*phi(x(j),t(i));
(dt/2)*phi(x(j),t(i+1))+eta(x(j),t(i+1))== eta(x(j),t(i))-(dt/2)*phi(x(j),t(i));
end
end
% update V
V=VCN;
end
% plot solution
set(fig,'YData',V(1,:));
set(tit,'String',strcat('n =',num2str(n)));
drawnow;
function x=x(j);
x(j)=j*dx
end
function t=t(i);
t(i)=i*dt
end
function eta=eta(x,t);
end
function theta=theta(x,t);
end
function phi=phi(x,t);
phi(x,t)=beta*(1-eta(x,t))*e^(-epsilon/(theta(x,t)+theta0));
end
Noemi Zeraick Monteiro on 15 Jun 2019
Hi. Sorry, I've got a hard activity and I don't have experience on MATLAB, indeed my deadline is hard, so I'm sorry for my mistakes. I've tried again I got the following error.