[FMDenaro]

Originally Posted by HumanistEngineer

It was my mistake is that I shall consider the equation as:



Now I focus on deriving the derivative approximation as to space dependent (x) _eff

But the original paper you linked shows Eq.(2) in a non-divergence form...

[HumanistEngineer]

Still I encounter the same problem (abnormal temperature increase) when I include the unsteady thermal diffusivity in the second order derivative.

Here is the Matlab code for the Forward Time Central Space (FTCS) for

function [T] = CDR_FTCS_InletMixing(t_final,n)
%% 1D Convection-Diffusion-Reaction (CDR) equation - (Forward-Time Central-Space)
%   solves a d2u/dx2 + b du/dx + c u = du/dt

%% References
%   pg. 03 - https://www.sfu.ca/~rjones/bus864/notes/notes2.pdf
%   Azad et al - Stability analaysis of finite difference schemes for an advection diffusion equation 

tStart = tic;

%% INPUT
%n  [-]       Number of nodes
%t_final : Total time period for simulation [s]

% Tank Dimensions
height=1.05;        % [m]       Tank height
v=5.2397e-4;        % [m/s]     Water velocity


dx=height/n;        % [m]       mesh size

dt=0.05;            % [s] time step 
maxk=round(t_final/dt,0);    % Number of time steps

NumParam=sprintf('The step size is %f and the time step is %f',dx, dt)

% Constant CDR coefficients | a d2u/dx2 + b du/dx + c u = du/dt
a_const=0.15e-6;       	% [m2/s]Thermal diffusivity 
b=-v;               % Water velocity
c=0;                %!!! Coefficient for heat loss calculation

st1=2*a_const*dt/(dx^2)
st2=v*dt/(2*dx)

%% Inlet Mixing - Decreasing hyperbolic function through tank height
eps_in=20;             % Inlet mixing magnification on diffusion
A_hyp=(eps_in-1)/(1-1/n);
B_hyp=eps_in-A_hyp;

Nsl=1:1:n;

eps_eff=A_hyp./Nsl+B_hyp;

a=a_const.*eps_eff;         % Case 0 
% a=ones(n,1)*a_const;        % Case 1
% a=ones(n,1)*a_const*500;	  % Case 2

%% Initial condition - Tank water at 40 ?C
T=zeros(n,maxk);
T(:,1)=40;

%% Boundary Condition
T(1,:)=90;

%% Tridiagonal Matrix @Right-hand side
subR(1:n-1)=dt*a(1:n-1)/(dx^2)-b*dt/(2*dx);  % Sub diagonal - Coefficient of u_i-1,n
maiR(1:n-0)=1-2*a*dt/(dx^2); % Main diagonal - Coefficient of u_i,n
supR(1:n-1)=a(1:n-1)*dt/(dx^2)+b*dt/(2*dx);  % Super diagonal - Coefficient of u_i+1,n
tdmR=diag(maiR,0)+diag(subR,-1)+diag(supR,1);

%% Boundary Condition - Matrices
tdmR(1,1)=1; tdmR(1,2)=0;
tdmR(end,end-1)=0; tdmR(end,end)=1; 

for j=2:maxk % Time Loop
    Tpre=T(:,j-1);
    T(:,j)=tdmR*Tpre;
    
    if T(end-1,j)>=89.9
        T(:,j+1:end)=[];
        finishedat=j*dt;
        ChargingTime=sprintf('Charging time is %f [s]', finishedat)
      	tElapsed = toc(tStart);
        SimulationTime=sprintf('Simulation time is %f [s]',tElapsed)
        return
    end
end
tElapsed = toc(tStart);
SimulationTime=sprintf('Simulation time is %f [s]',tElapsed)
end

and for:

function [T] = CDR_FTCS_InletMixing_Last(t_final,n)
%% 1D Convection-Diffusion-Reaction (CDR) equation - (Forward-Time Central-Space)
%   solves a d2u/dx2 + b du/dx + c u = du/dt

tStart = tic;

%% INPUT
%n  [-]       Number of nodes
%t_final : Total time period for simulation [s]

% Tank Dimensions
height=1.05;        % [m]       Tank height
v=5.2397e-4;        % [m/s]     Water velocity

dx=height/n;        % [m]       mesh size

dt=0.005;            % [s] time step 
maxk=round(t_final/dt,0);    % Number of time steps

NumParam=sprintf('The step size is %f and the time step is %f',dx, dt)

% Constant CDR coefficients | a d2u/dx2 + b du/dx + c u = du/dt
a_const=0.15e-6;       	% [m2/s]Thermal diffusivity 
b=-v;               % Water velocity
c=0;                %!!! Coefficient for heat loss calculation

st1=2*a_const*dt/(dx^2)
st2=v*dt/(2*dx)

%% Inlet Mixing - Decreasing hyperbolic function through tank height
eps_in=5;             % Inlet mixing magnification on diffusion
A_hyp=(eps_in-1)/(1-1/n);
B_hyp=eps_in-A_hyp;

Nsl=1:1:n;

eps_eff=A_hyp./Nsl+B_hyp;

a(:,1)=a_const.*eps_eff;         % Case 0 
% a=ones(n,1)*a_const;        % Case 1
% a=ones(n,1)*a_const*500;	  % Case 2

%% Initial condition - Tank water at 40 ?C
T=zeros(n,maxk);
T(:,1)=40;

%% Boundary Condition
T(1,:)=90;

%% Tridiagonal Matrix @Right-hand side
sub_subR(1:n-2)=dt/(4*dx^2);
mai_subR(1:n-1)=dt/(dx^2);
sup_subR(1:n-2)=-dt/(4*dx^2);
tdm_subR=diag(mai_subR,0)+diag(sub_subR,-1)+diag(sup_subR,1);
subR(1:n-1)=tdm_subR*a(1:n-1)-b*dt/(2*dx);

maiR(1:n)=1-2*a*dt/(dx^2);

sub_supR(1:n-2)=-dt/(4*dx^2);
mai_supR(1:n-1)=dt/(dx^2);
sup_supR(1:n-2)=dt/(4*dx^2);
tdm_supR=diag(mai_supR,0)+diag(sub_supR,-1)+diag(sup_supR,1);
supR(1:n-1)=tdm_supR*a(1:n-1)+b*dt/(2*dx);

tdmR=diag(maiR,0)+diag(subR,-1)+diag(supR,1);

%% Boundary Condition - Matrices
tdmR(1,1)=1; tdmR(1,2)=0;
tdmR(end,end-1)=0; tdmR(end,end)=1; %%% Something WRONG!!!!!

for j=2:maxk % Time Loop
    Tpre=T(:,j-1);
    T(:,j)=tdmR*Tpre;
    
    if T(end-1,j)>=89.9
        T(:,j+1:end)=[];
        finishedat=j*dt;
        ChargingTime=sprintf('Charging time is %f [s]', finishedat)
      	tElapsed = toc(tStart);
        SimulationTime=sprintf('Simulation time is %f [s]',tElapsed)
        return
    end
end
tElapsed = toc(tStart);
SimulationTime=sprintf('Simulation time is %f [s]',tElapsed)
end

[FMDenaro]

That is written as a function not as a full program? I cannot run without the parameter at line 20

[HumanistEngineer]

Originally Posted by FMDenaro

That is written as a function not as a full program? I cannot run without the parameter at line 20

Yes it is a function. Please make your run as "[Tlast] = CDR_FTCS_InletMixing(3000,1000)" as a simulation for 3000 seconds with 1000 nodes.

Here is the plotter script (not a function this time) but please make sure that function output is Tlast (this plot script based on Tlast variable name!).

Tlast = Tlast(:,any(Tlast,1)); n=length(Tlast(:,1)); h=1:n-1; h=(1.05/n)*h';
slmin=1; slmax=size(Tlast,2); hold on;
plot(h,Tlast(2:end,1)); axis([0 1.2 20 120]);
hsl = uicontrol('Style','slider','Min',slmin,'Max',slmax,...
'SliderStep',[1 1]./(slmax-slmin),'Value',1,...
'Position',[5 5 200 20]);
set(hsl,'Callback',@(hObject,eventdata) plot(h,Tlast(2:end,round(get(hObject,'Value')))) )

[FMDenaro]

I wrote few line in "educational" style, check it. I see reasonable results.








clear, format long e

% Tank Dimensions
height=1.05; % [m] Tank height
v=5.2397e-4; % [m/s] Water velocity
n=1000
dx=height/n; % [m] mesh size
dt=0.1; % [s] time step
n_time_limit=1000
a_const=0.15e-6; % [m2/s]Thermal diffusivity
x(1:n+1)=dx*(1:n+1)-dx;




%% Inlet Mixing - Decreasing hyperbolic function through tank height
eps_in=20; % Inlet mixing magnification on diffusion
A_hyp=(eps_in-1)/(1-1/n);
B_hyp=eps_in-A_hyp;

Nsl=1:1:n+1;

eps_eff=A_hyp./Nsl+B_hyp;

coeff=a_const.*eps_eff; % Case 0

plot(x,coeff)
pause

stab=coeff.*dt/dx/dx;
plot(x,stab)
pause

Told(2:n)=40;
Told(n+1)=Told(n);
Told(1)=90;

n_time_limit=10000;

for k=1:n_time_limit

for i=2:n
Tnew(i)=Told(i)-v*dt*(Told(i+1)-Told(i-1))/(2*dx)+dt*coeff(i)*(Told(i+1)-2*Told(i)+Told(i-1))/(dx*dx);
end
Told(2:n)=Tnew(2:n);
Told(n+1)=Tnew(n);
k
plot(x,Told)
pause

end

[HumanistEngineer]

Originally Posted by FMDenaro

I wrote few line in "educational" style, check it. I see reasonable results.

Thank you Filippo. I could the reach the same reasonable result by basing on the for loop use as yours to calculate the new temperatures. I guess that the tridiagonal matrix formed in my previous codes are misleading. Thanks again.

[FMDenaro]

Ok, Good!