[HumanistEngineer]

Thank you in advance,

I could manage to build a numeric model for a sensible stratified storage tank by means of solving the convection-diffusion equation (1D) by use of Crank-Nicolson scheme.



I try to implement a mixing effect that is changing in a hyperbolic decaying form from the top of the tank (inlet is here) to the bottom. In the above equation, eps_eff is this mixing effect. The details can be found here: Zurigat et al - Stratified thermal storage tank inlet mixing characterization.

Please see the Case 0, 1, and 2 under the section %% Inlet Mixing in the Matlab code given below.

For Case 1 - No mixing effect included so a constant thermal diffusivity through the tank height, here is the result (simulation results are solid lines compared with the experimental result shown with circles):


When Case 0 is active, the thermal diffusivity is replaced from constant to a decaying form through the tank height via multiplication with eps_eff in the code, the results become weird:


If Case 2, the thermal diffusivity is magnified up to a level up at a constant rate at all tank levels, the temperature does not increase above the inlet temperature but Case 0 all the time results in abnormal temperature increase.

Would you please check my codes and/or explain to me why this abnormal temperature increase happens (stability, oscillation,... issues?)? Shall I try upwind or any other scheme in my solution considering the inlet mixing effect? Or any other suggestions?

Here is the Matlab code:

function [T,dh,dt] = CDR_CrankNicolson_Exp_InletMixing(t_final,n)
%% 1D Convection-Diffusion-Reaction (CDR) equation - Crank-Nicolson scheme
%   solves a d2u/dx2 + b du/dx + c u = du/dt

%% References
% pg. 28 - http://www.ehu.eus/aitor/irakas/fin/apuntes/pde.pdf
% pg. 03 - https://www.sfu.ca/~rjones/bus864/notes/notes2.pdf
% http://nptel.ac.in/courses/111107063/24

tStart = tic; % Measuring computational time

%% INPUT
% t_final   : overall simulation time [s]
% n         : Node number [-]

%% Tank Dimensions
h_tank=1;       % [m]       Tank height
d_tank=0.58; 	% [m]       Tank diameter
vFR=11.75;    	% [l/min]   Volumetric Flow Rate
v=((vFR*0.001)/60)/((pi*d_tank^2)/4);    % [m/s]     Water velocity

dh=h_tank/n;   	% [m] mesh size

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

% 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

%% Inlet Mixing - Decreasing hyperbolic function through tank height
eps_in=2;             % 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*20;	% Case 2

%% Initial condition - Tank water at 20 degC
T=zeros(n,maxk);

T(:,1)=20;

%% Boundary Condition - Inlet temperature at 60 degC
T(1,:)=60;

%% Formation of Tridiagonal Matrices
%   Tridiagonal Matrix @Left-hand side
subL(1:n-1)=-(2*dt*a(1:n-1)-dt*dh*b);         	% Sub diagonal - Coefficient of u_i-1,n+1
maiL(1:n-0)=4*dh^2+4*dt*a-2*dh^2*dt*c;          % Main diagonal - Coefficient of u_i,n+1
supL(1:n-1)=-(2*dt*a(1:n-1)+dt*dh*b);         	% Super diagonal - Coefficient of u_i+1,n+1
tdmL=diag(maiL,0)+diag(subL,-1)+diag(supL,1);
%   Tridiagonal Matrix @Right-hand side
subR(1:n-1)=2*dt*a(1:n-1)-dt*dh*b;            	% Sub diagonal - Coefficient of u_i-1,n
maiR(1:n-0)=4*dh^2-4*dt*a-2*dh^2*dt*c;          % Main diagonal - Coefficient of u_i,n
supR(1:n-1)=2*dt*a(1:n-1)+dt*dh*b;           	% Super diagonal - Coefficient of u_i+1,n
tdmR=diag(maiR,0)+diag(subR,-1)+diag(supR,1);

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

tdmR(1,1)=1; tdmR(1,2)=0;
tdmR(end,end-1)=1; tdmR(end,end)=0;

MMtr=tdmL\tdmR;

%% Solution - System of Equations
for j=2:maxk % Time Loop
Tpre=T(:,j-1);

T(:,j)=MMtr*Tpre;

if T(end,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

end

[arjun]

Based on the graphs, your diffusion term is creating a source that should not be there. If i were you i will just set the diffusion term to 0 in the cell connected to boundary and run it.

[FMDenaro]

I see a poor accordance also in the constant thermal diffusivity case...


However, I haven't see the code but if you manage the eps_eff in the Crank-Nicolson scheme have you correctly considered the effect in the time-variying case? The matrix will change at each time step.


Then, how do you discretize the spatial derivatives?

[HumanistEngineer]

Originally Posted by arjun

Based on the graphs, your diffusion term is creating a source that should not be there. If i were you i will just set the diffusion term to 0 in the cell connected to boundary and run it.

Thank you arjun for your quick reply. I tried your suggestion but didn't work. Again I have this abnormal temperature increase.

[HumanistEngineer]

Originally Posted by FMDenaro

I see a poor accordance also in the constant thermal diffusivity case...


However, I haven't see the code but if you manage the eps_eff in the Crank-Nicolson scheme have you correctly considered the effect in the time-variying case? The matrix will change at each time step.


Then, how do you discretize the spatial derivatives?

Thank you FMDenaro. I consider the thermal diffusivity constant (so defined for a mean temperature). So at each time step, the thermal diffusivity or thermal diffusivity times eps_eff stays constant.

Coefficients of approximated temperatures in the Crank-Nicolson scheme applied stay same with time so I calculate it once at the beginning (a fixed matrix). Do you think the problem is this?

[FMDenaro]

Originally Posted by HumanistEngineer

Thank you FMDenaro. I consider the thermal diffusivity constant (so defined for a mean temperature). So at each time step, the thermal diffusivity or thermal diffusivity times eps_eff stays constant.

Coefficients of approximated temperatures in the Crank-Nicolson scheme applied stay same with time so I calculate it once at the beginning (a fixed matrix). Do you think the problem is this?

You wrote “When Case 0 is active, the thermal diffusivity is replaced from constant to a decaying form ...”

Therefore the coefficients of the matrix are time-dependent not constant

[HumanistEngineer]

Originally Posted by FMDenaro

You wrote “When Case 0 is active, the thermal diffusivity is replaced from constant to a decaying form ...”

Therefore the coefficients of the matrix are time-dependent not constant

The coefficients are dependent on the tank height but all are constant through time.

[FMDenaro]

Originally Posted by HumanistEngineer

The coefficients are dependent on the tank height but all are constant through time.

But are you including in the CN integration also the convection part??


Have you checked that a purely explicit FTCS works fine?


Increasing the eps_eff greater to 1 you should see exactly the opposite effect, check also the correct sign in the coefficients, what you het is a sort of anti-diffusion.

[HumanistEngineer]

Originally Posted by FMDenaro

But are you including in the CN integration also the convection part??

Yes. But in another try, I applied a hybrid scheme, considering the upwind scheme for the convection term while the rest with CN scheme simultaneously (the reason is to avoid spurious oscillations in convective-dominant cases). Still same problem occured (abnormal temperature increase) in this hybrid scheme after eps_eff applied.

Originally Posted by FMDenaro

Have you checked that a purely explicit FTCS works fine?

No. What will you recommend of? I am new in this finite difference approximation methods (my understanding in their superiority -i.e. accuracy, stableness, etc. is missing)! Which scheme shall I use? If possible, any source can be rewarding in my solution.

Originally Posted by FMDenaro

Increasing the eps_eff greater to 1 you should see exactly the opposite effect, check also the correct sign in the coefficients, what you het is a sort of anti-diffusion.

I will.

Thank you.

[FMDenaro]

As you are not experienced, start solving the explicit first order integration in time. Use central second order formula for the space derivatives.
At high values of eps_eff the numerical stability constraint is mainly based on

dt*diff_coef/h^2<0.5

For convection dominated case you need a fine grid to avoid numerical oscillations

[HumanistEngineer]

Originally Posted by FMDenaro

As you are not experienced, start solving the explicit first order integration in time. Use central second order formula for the space derivatives.
At high values of eps_eff the numerical stability constraint is mainly based on

dt*diff_coef/h^2<0.5

For convection dominated case you need a fine grid to avoid numerical oscillations

I will try this. Thank you for your time.

[HumanistEngineer]

Originally Posted by FMDenaro

As you are not experienced, start solving the explicit first order integration in time. Use central second order formula for the space derivatives.
At high values of eps_eff the numerical stability constraint is mainly based on

dt*diff_coef/h^2<0.5

For convection dominated case you need a fine grid to avoid numerical oscillations

I changed my code as to Forward Time Central Difference scheme. By changing the mesh size and/or time step, I could have a stable result without mixing effect. However, after I involve of mixing efficiency as a hyperbolic decaying form through the tank height, I again encounter the same problem with abnormal temperature increase as same as in the graph that I present in my first message.

[FMDenaro]

Originally Posted by HumanistEngineer

I changed my code as to Forward Time Central Difference scheme. By changing the mesh size and/or time step, I could have a stable result without mixing effect. However, after I involve of mixing efficiency as a hyperbolic decaying form through the tank height, I again encounter the same problem with abnormal temperature increase as same as in the graph that I present in my first message.

What is the value of alpha*eps_eff? What about the bc at the right?

[HumanistEngineer]

Originally Posted by FMDenaro

What is the value of alpha*eps_eff? What about the bc at the right?

alpha, the thermal diffusivity, is constant as 0.15e-6 m2/s. When magnified with eps_eff, it becomes as an array through the tank height as below (PrintScreen from the Matlab code), which is an example for a node number of 10:



The boundary condition at the right/bottom is dT/dt=0.

[FMDenaro]

Therefore, if I am right this problem has a thermal coefficient that increases going to the right of the diagram.
Just as a test, could you use the same law for eps_eff but changing only the sign so to decrease it?

[HumanistEngineer]

Originally Posted by FMDenaro

Therefore, if I am right this problem has a thermal coefficient that increases going to the right of the diagram.
Just as a test, could you use the same law for eps_eff but changing only the sign so to decrease it?

Here is the result when I reverse the direction of eps_eff (example at 90 °C charge temperature while initial temperature in the tank at 40 °C):

[FMDenaro]

The case now shows a diminishing temperature, right? So, it seems that the numerical solution describes opposite physical behaviors.

[HumanistEngineer]

Originally Posted by FMDenaro

The case now shows a diminishing temperature, right? So, it seems that the numerical solution describes opposite physical behaviors.

Weirdly the hyperbolic decaying form of eps_eff cause this issue, either in the original case or in this reverse direction case (the last having issue on the right boundary). When I multiply a constant of eps_eff with all elements of the thermal diffusivity, alpha, the results are not obtained with this abnormal temperature increase but a high rate of diffusion as expected.

[FMDenaro]

Originally Posted by HumanistEngineer

Weirdly the hyperbolic decaying form of eps_eff cause this issue, either in the original case or in this reverse direction case (the last having issue on the right boundary). When I multiply a constant of eps_eff with all elements of the thermal diffusivity, alpha, the results are not obtained with this abnormal temperature increase but a high rate of diffusion as expected.

Could you post your Matlab code for the simple FTCS case? I want try to have a check

[HumanistEngineer]

Originally Posted by FMDenaro

Could you post your Matlab code for the simple FTCS case? I want try to have a check

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



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