In-Situ Combustion Model - MATLAB code

5 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:
Follow the code. I think I need to add some informartion about approximations, for instance, for
. 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

Accepted Answer

John D'Errico
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.
  2 Comments
Noemi Zeraick Monteiro
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
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.
kzAyr.png

Sign in to comment.

More Answers (0)

Categories

Find more on MATLAB in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!