Linear Driving Force Model for Gas Adsorption in Fixed Bed

Dear Community

I am currently trying to solve a system of ODE containing energy balance equation during gas adsorption in a fixed bed. I am atatching my initial script herewith. There are four equstions. Eq, 1 to Eq. 3 respectively are energy balance for water passing though the bed (i.e., taking away adsorption heat), adsorption heat generation in the bed and heat transferred to the gas domain from bed during the adsoption. In Eq. 2 there is a heat generation term H_ads which should be multiplied with differential amount of gas adsorbed at each timestep (dX_ads/dt). Eq. 4 represents in Linear Driving Force Model as: dX_ads/dt = K_LDF*(X_eq-X_ads), where X_eq and K_LDF are constant.

so, basically, say at time step t=2 sec, H_ads should be multiplied by the differential value of dX_ads/dt at time step t=1sec, and X_ads in eq. 4 should be updated to the cumulative value of X_ads for all previous time steps. Currently, I am unable to program eq, 2 and eq. 4 in that way in my mode. If someone can help me to overcome this difficulty, I shall remain grateful.

I have attached the matlab file herewith too.

Shamim, UTokyo

clear all
close all
clc


% Constants_water
hw= 800;
A_tube=0.091985;
rho_water=1000;
Cpw=4182.0;
Mw=85.086;
C_w_num= (hw*A_tube)/(Mw*Cpw);
display (C_w_num)
%---------------------------
% Constants_gas
hgas= 4.11;
A_fin=0.6018;
rho_gas=1.98;
Cpgas=37.35;
Mgas=0.0072864;
C_gas_num= (hgas*A_fin)/(Mgas*Cpgas);
display (C_gas_num)
%---------------------------
% Constants_bed
Cpmof=658.9;
Mmof=0.15;
C_bedA_num=(hgas*A_fin)/(Mmof*Cpmof);
C_bedB_num=(hw*A_tube)/(Mmof*Cpmof);
display (C_bedA_num)
display (C_bedB_num)
%---------------------------------------------
% Adsorption parameters 
H_ads_num= 150;   % J/mol
K_LDF_num= 6.3211E-3; % s-1
MOF_num=0.15;
X_eq_num= 30; % mol/kg

%C_ads_num=K_LDF*Mmof;



% Equations 
syms T_w(t) T_bed(t)  T_gas(t) X_ads(t) C_w C_gas C_bedA C_bedB X_eq H_ads K_LDF MOF
display (syms)

eq1 = diff(T_w,t) == (C_w*(T_bed - T_w));   % water
eq2 = diff(T_bed,t) == ((C_bedA*(T_gas - T_bed))+ (C_bedB*(T_w-T_bed)))+ H_ads*MOF*X_ads; %bed
eq3 = diff(T_gas,t) == (C_gas*(T_bed - T_gas)); %gas
eq4 = diff(X_ads,t)== (K_LDF*(X_eq-X_ads)); %LDF

eqns = [eq1 eq2 eq3 eq4];
[V, vars] = odeToVectorField(eqns);


%Note----------------------------------------------------------------------------
% In eq4, i get the cumulative value of X_ads which is multiplied with
% H_ads in eq.2. I need the differential value of X_ads at each time step
% to be multiplied with H_ads in eq.2.  
%--------------------------------------------------------------------------------------------------------------------

%time step
ti=0; tf=300; tspan = [ti tf];

%Initial values 
Tw_initial =20;
Tbed_initial =20;
Tgas_initial=20;
X_ads_initial=2.5;

yInit = [Tw_initial Tbed_initial Tgas_initial X_ads_initial];


%Solve ODE45
Fun = matlabFunction(V, 'vars' ,{'t','Y', 'C_w', 'C_gas', 'C_bedA', 'C_bedB', 'X_eq', 'H_ads', 'K_LDF', 'MOF' })

ODE = @(t,Y)Fun(t,Y, C_w_num, C_gas_num, C_bedA_num, C_bedB_num, X_eq_num, H_ads_num, K_LDF_num, MOF_num)

[t,res] = ode45(ODE,tspan,yInit);


% plot results
subplot(2,2,1);
plot(t, res(:,1))
xlabel('time')
ylabel('T-water')


subplot(2,2,2);
plot(t, res(:,2))
xlabel('time')
ylabel('T-bed')


subplot(2,2,3);
plot(t, res(:,3))
xlabel('time')
ylabel('T-gas')

subplot(2,2,4);
plot(t, res(:,4))
xlabel('time')
ylabel('Adsorbed Amount')

[Torsten]

Please include the equations you try to solve in a mathematical notation.

[JUBAIR AHMED SHAMIM]

Hi Torsten Many thanks for your message. I have included the equations in mathematical notations in attached JPG file. Also, so far, I have made some progress. Now my gas temperature and water temperature follows the expected trend I want. But bed temperature becomes constant. But bed temperature should also have the same trend as gas and water temperature. Not sure where is the problem. These are all first order ODE. But I used "ODE to vector field'' from an example code. This command usually used for transforming the second order equations into first order equations. Is that can be a problem? if so, how to use ODE 45 for a system of first order equations for my case I am not sure yet. Also, attaching my latest code using some changes made by Sam Chak.

[Torsten]

Xeq must depend on the concentration of the adsorptive. And you don't have a balance equation for the adsorptive. In my opinion, your model as written is at least strongly incomplete.

% In eq4, i get the cumulative value of X_ads which is multiplied with

% H_ads in eq.2. I need the differential value of X_ads at each time step

% to be multiplied with H_ads in eq.2.

Just insert the differential ; it is K_LDF*(X_eq-X_ads).

[JUBAIR AHMED SHAMIM]

Hi Torsten, Thanks for your reply. I have done exactly same thing you said. just multiply with K_LDF*(X_eq-X_ads). Regarding model, yes, it's not complete yet. I will keep adding auxiliary equations from now on. I just made the framework to solve for governing equations at first.

Answer 1

[Sam Chak]

Hi @JUBAIR AHMED SHAMIM

Editted to run your supplied code. The state responses (Tbed, Twater, Tgas) converge to the same value around 20.07.

% Water parameters
hw         = 800;
A_tube     = 0.091985;
rho_water  = 1000;
Cpw        = 4182.0;
Mw         = 85.086;
C_w_num    = (hw*A_tube)/(Mw*Cpw);
display (C_w_num)

C_w_num = 2.0681e-04

% Gas parameters
hgas       = 4.11;
A_fin      = 0.6018;
rho_gas    = 1.98;
Cpgas      = 37.35;
Mgas       = 0.0072864;
C_gas_num  = (hgas*A_fin)/(Mgas*Cpgas);
display (C_gas_num)

C_gas_num = 9.0885

% Bed parameters
Cpmof      = 658.9;
Mmof       = 0.15;
C_bedA_num = (hgas*A_fin)/(Mmof*Cpmof);
C_bedB_num = (hw*A_tube)/(Mmof*Cpmof);
display (C_bedA_num)

C_bedA_num = 0.0250

display (C_bedB_num)

C_bedB_num = 0.7446

% Adsorption parameters 
H_ads     = 6000;       % J/mol
K_LDF_num = 6.3211E-2;  % s-1
Mass_MOF  = 0.15;
X_eq_num  = 30;         % mol/kg
C_ads_num = (H_ads*Mass_MOF)/(Mmof*Cpmof); 
display (C_ads_num)

C_ads_num = 9.1061

% Equations 
syms T_w(t) T_bed(t)  T_gas(t) X_ads(t) C_w C_gas C_bedA C_bedB X_eq H_ads K_LDF MOF C_ads
display (syms)

  14×1 cell array

    {'C_ads' }
    {'C_bedA'}
    {'C_bedB'}
    {'C_gas' }
    {'C_w'   }
    {'H_ads' }
    {'K_LDF' }
    {'MOF'   }
    {'T_bed' }
    {'T_gas' }
    {'T_w'   }
    {'X_ads' }
    {'X_eq'  }
    {'t'     }

eq1 = diff(T_w,t)   ==    C_w*(T_bed - T_w);                                                        % water
eq2 = diff(T_bed,t) == C_bedA*(T_gas - T_bed) + C_bedB*(T_w - T_bed) + C_ads*K_LDF*(X_eq - X_ads);  % bed 
eq3 = diff(T_gas,t) ==  C_gas*(T_bed - T_gas);                                                      % gas
eq4 = diff(X_ads,t) ==  K_LDF*(X_eq  - X_ads);                                                      % LDF

eqns      = [eq1 eq2 eq3 eq4];
[V, vars] = odeToVectorField(eqns)  % check the order of state equations

V = 

vars = 

% time interval
ti    = 0; 
tf    = 300; 
tspan = [ti tf];

% initial values 
Tbed_initial  = 20;
Tw_initial    = 20;
Tgas_initial  = 20;
X_ads_initial = 2.5;
yInit         = [Tbed_initial Tw_initial Tgas_initial X_ads_initial];

% ODE45 solver
odefun   = matlabFunction(V, 'vars' ,{'t', 'Y', 'C_w', 'C_gas', 'C_bedA', 'C_bedB', 'X_eq', 'K_LDF', 'C_ads'});
ODE      = @(t, Y) odefun(t, Y, C_w_num, C_gas_num, C_bedA_num, C_bedB_num, X_eq_num, K_LDF_num, C_ads_num);
[t, y] = ode45(ODE, tspan, yInit);

Tbed = y(:,1);
Tw   = y(:,2);
Tgas = y(:,3);
Xads = y(:,4);

% plot results according to the order of equations [eq1 eq2 eq3 eq4]
subplot(2,2,1);
plot(t, Tw), grid on, % ylim([20.06 20.08])
xlabel('Time')
ylabel('T-water')

subplot(2,2,2);
plot(t, Tbed), grid on, % ylim([20.06 20.08])
xlabel('Time')
ylabel('T-bed')

subplot(2,2,3);
plot(t, Tgas), grid on, % ylim([20.06 20.08])
xlabel('Time')
ylabel('T-gas')

subplot(2,2,4);
plot(t, Xads), grid on
xlabel('Time')
ylabel('Adsorbed Amount')

[Sam Chak]

Hi @JUBAIR AHMED SHAMIM,

I have edited and correct your code in my Answer. It seems that Tbed and Tgas will converge to the initial value of Twater, regardless of the initial values of Tbed and Tgas. I'm unfamiliar with adsorption systems. Are there partial derivative in the system of equations?

[JUBAIR AHMED SHAMIM]

HI Sam Chak, Many thanks for your reply. This is the other version of the code (attached) I just made without using odeToVectorField command. You are right. All equations here are first order equation. So, we don't have to use odeToVectorField which is used to convert second order equations into first order equations. It gives similar trend to your results. Thanks.

[Sam Chak]

Hi @JUBAIR AHMED SHAMIM

Good to hear that another version your code produces the same results. If you find the corrected code helpful in guiding you making the framework to solve for governing equations, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Thanks a bunch!

[JUBAIR AHMED SHAMIM]

Hi Sam I have one more question. If I want to add two more terms in Eq. (1) and Eq. (2) of the code you corrected yesterday, then how do i do it? I have attached the equations in the jpg file herewith. But now it becomes a system of ODE and PDE. Can it be solved using pdepe commend. I found example where a system of PDE was solved using this command. But not sure it can be used for system of ODE and PDE. Any input will be very helpful.

[Torsten]

You must have a 5th equation for the concentration of the adsorptive depending on x. Otherwise you cannot compute X_eq(x).

"pdepe" is not suited to solve your PDE system because it is not designed to solve ODEs or equations that only contain first-order derivatives in space like your equation for T_water.

You will have to discretize the spatial derivatives and solve the resulting system of ordinary differential equations for the functions in the grid points using ode15s. Look up "method-of-lines" for more details.

[JUBAIR AHMED SHAMIM]

Hi Torsten, I also thought so. Discretization will be necessary. I will try to understand coding of that method from now using method of lines. Regarding X_eq(x), I will add the fifth equation later. I think handling all the equations same time is more complex. So I will build the code block by block. For now, I am assuming pressure and Temperature in my system is constant. So using the adsorption isotherm from a paper for CO2 adsorption into MIL-101, I can get a X_eq at defined pressure and Temperature.

[Torsten]

Especially the decrease of adsorptive concentration in x-direction will change X_eq. I don't know whether there are big changes in P and T for your process.

[JUBAIR AHMED SHAMIM]

Hi Torsten, I am sorry to ask you again. Today I finally tried to solve the 1D model using method of lines, ODE45 and central difference approximation as you suggested last time. I am attaching the equations I tried to solve and also the matlab code. I tired to build this code from an example that used backward difference scheme. When I changed to central difference scheme, I am getting error. I understand somewhere I made mistake inm defining matrix but I can not solve it.

I am also attaching an image for expected trends for temperatures. But unfortunately, I can not get the expected trends too. If you can take a look and tell me how to solve, I shall be grateful. Thank you for your time and cooperation agaian.

[Torsten]

You need to incorporate boundary conditions for Tw at x = 0 and for Tbed at both ends (x=0 and x=3).

This way, errors like trying to access Tbed(n+1) in the expression

DTbed2DX2(i) =  (Tbed(i+1)-2*Tbed(i)+Tbed(i-1))/(dX(i)^2);  % second derivative

will vanish.

Take care that you really give values to the time derivatives for 1 <= i <= n, not only 2 <= i <= n like in your code.

Since the water flows from left to right, its gradient should be taken upwind:

dTwdx(i) = (Tw(i)-Tw(i-1))/(X(i)-X(i-1))

[JUBAIR AHMED SHAMIM]

HI Torsten, I changed the formula for dTwdx(i) as you suggested. Also, I found the error for central difference approximation. Now the code is running but I don't get expected results. Boundary conditions are implemented at both ends I can see from the matrix of Tbed. I am attaching a screenshot too.But for some reasons results are erroneous and getting negative. Not sure where did I make the mistake. Could you please take a look into my code please one more time:

Another thing i have realized that, for Tbed i have to use Neumann BC i.e., dTbed/dX =0 at both ends of the bed. however, coding is getting complex for that and i am not sure how to implement it. Aslo, the equations for Q and Tgas are not PDE, so only IC should be ok. but currently, i used dirichlet BC for these two. further, for Twater i want to use dirichlet at inlet (constant inlet temperature) but Neumann at outlet dTw/dx=0.

clear all
close all
clc

% Discretization parameters
L = 3; % length of tube
N=10; % no. of discretization in X direction
X = linspace(0,L,N); % discretise the domain 
n = numel(X);

% time interval
ti=0;
tend=50;   % end time 
array_length=tend;
tspan = linspace (ti,tend,array_length); % given time span
%-------------------------------------------------------------

% Dirichlet boundary condition at x=0 (no of BC = no of eqns)
Tgas0      = 20*ones(n,1);        
Tbed0     = 20     *ones(n,1);     
Tw0     = 20 *ones(n,1);     
Q0      = 2.5*ones(n,1);     
% X0      = 0     *ones(n,1);     % Conversion rate
%-------------------------------------------------------------

% setting up the solver
y0 = [Tgas0; Tbed0; Tw0; Q0];
[T,Y] = ode45(@(t,y) solution(t,y,X,n),tspan,y0);
Y=Y.';
%-------------------------------------------------------------

% saving results
Tgas = Y(    1 :  n,:);
Tbed = Y(  n+1:2*n,:);
Tw = Y(2*n+1:3*n,:);
Q  = Y(3*n+1:4*n,:);
% X  = Y(4*n+1:5*n,:);
%-------------------------------------------------------------

figure
subplot(2, 1, 1);
plot (T, Tw(N,1:tend), 'g', T, Tbed (N*0.5,1:tend), 'red',   T, Tgas (N*0.5,1:tend),'b');
%plot ( T, Tw (N,1:2*tend), 'red');
xlabel('time')
ylabel('Temperature')   % how to define water temperature at the outlet 
legend('Water-outlet', 'Bed_middle','Gas_middle');


subplot(2, 1, 2);
plot(T, Q (N, 1:tend));
title('Loading');
xlabel('Time');
ylabel('Loading at middle [mol/kg_m_o_f]');
%----------------------------------------------------------------


function DyDt = solution(t,y,X,n)

% Define Constants 
% h=0.1; % step-size
%water
rho_w =1000; % Kg/m3
Cp_w = 4200; % J/kg.K
V_w = 0.000157411; %m3 from excel calc. sheet 
Flow_rate_water = 0.5; % LPM 
H_w =1500; % W/m2.K 1500 as fitting /Heat Transf. Coeff.
A_tube= 0.091985; % m2

% Gas
hgas       = 411;
A_fin      = 0.6018;
rho_gas    = 1.98;
Cpgas      = 37.35;
Mgas       = 0.0072864;


%Bed
Cpmof      = 658.9;
Mmof       = 0.15;
conductance=0.002804812;

% Adsorption parameters 
H_ads  = 1500;   % J/mol
K_LDF  = 6.3211E-2; % s-1
Mass_MOF    = 0.15;
Q_eq   = 30;    % mol/kg


% Dependent variables for water energy balance 
V_dot_w = Flow_rate_water*1.66667E-5; %m3/sec.
C1_water = (V_dot_w/V_w);
C2_water = (H_w*A_tube)/(rho_w*Cp_w*V_w);
C_gas  = (hgas*A_fin)/(Mgas*Cpgas);
C_bedA = (hgas*A_fin)/(Mmof*Cpmof);
C_bedB = (H_w*A_tube)/(Mmof*Cpmof);
C_ads= ((H_ads*Mass_MOF)/(Mmof*Cpmof));
C_cond= (conductance*A_tube)/(Mmof*Cpmof);
%-------------------------------------------------------------

% Define matrix

%Change in variables with time
DTgasDt = zeros(n,1); 
DTbedDt = zeros(n,1); 
DTwDt = zeros(n,1); 
DQDt = zeros(n,1); 


% Change in variables with space
    DTwDX = zeros(n,1); % Change in temperature with space
    DTbedsqDXsq= zeros(n,1);

% Dependent variables
Tgas    = zeros(n,1);
Tbed    = zeros(n,1);
Tw    = zeros(n,1);
Q    = zeros(n,1);

% from saving results
Tgas = y(1:n);            
Tbed = y(n+1:2*n);        
Tw = y(2*n+1:3*n);      
Q  = y(3*n+1:4*n);  

%     X  = y(4*n+1:5*n);      % Conversion rate


% Define finite (central) difference approximation
for i=2:n-1
       
    %dX(i)= X(i+1)-X(i);
    DTwDX(i) =    (Tw(i)-Tw(i-1))/(X(i)-X(i-1));  % first derivate backward difference
    DTbed2DX2(i) =  (Tbed(i+1)-2*Tbed(i)+Tbed(i-1))/((X(i+1)-X(i))^2);  % second derivative central difference
end
 
% Define governing equations 

for i=2:n-1
    
    DTgasDt(i)= C_gas* (Tbed(i)-Tgas(i));

    %DTbedDt(i)= C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i));   %without second derivative      
    
    DTbedDt(i)= C_cond*DTbed2DX2(i)+ C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i)); % with second derivative              ;
    
    DTwDt(i) = C1_water* DTwDX(i) + C2_water*[Tbed(i) - Tw(i)] ;

    DQDt(i) = K_LDF*(Q_eq - Q(i));
end
 
% Return function
DyDt = [DTgasDt; DTbedDt;DTwDt; DQDt]; 

end

[Sam Chak]

Hi @JUBAIR AHMED SHAMIM

I recently viewed a video clip featuring an experiment conducted by Abhishek Aravalli, titled "Fixed Bed Gas Adsorption Column." I was intrigued by the experiment and had a question regarding the manipulable parameter that governs or influences the adsorption dynamical behavior of gases on solid materials. Could you kindly identify this parameter within the equations or images presented below?

https://www.youtube.com/watch?v=IJ2Az0aci8I

[JUBAIR AHMED SHAMIM]

Hi Sam chak thanks for your question. The parameter that controls the adsorption behavior in my model is K_LDF and Q_eq. Here, I just showed the governing equations. But in the complete model there will be several auxiliary equations too: one for K_LDF (this is the gas transport parameter within the pores and there are several ways to model it) and one for Q_eq (this is called adsorption isotherm). Here I am assuming constant values for now to setup the code for governing equations first. Then, I shall include auxiliary equations in the next step. Now I am facing problem to setup BCs in my model.

[Torsten]

for i=2:n-1
    
    DTgasDt(i)= C_gas* (Tbed(i)-Tgas(i));
    %DTbedDt(i)= C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i));   %without second derivative      
    
    DTbedDt(i)= C_cond*DTbed2DX2(i)+ C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i)); % with second derivative              ;
    
    DTwDt(i) = C1_water* DTwDX(i) + C2_water*[Tbed(i) - Tw(i)] ;
    DQDt(i) = K_LDF*(Q_eq - Q(i));
end

Your loop still runs from 2 to n-1. So you still don't supply boundary conditions for Tbed and Tw. And supplying DTgasDt amd DQDt for i = 1 and i = n should be not problem since the equations for Tgas and Q are ordinary differential equations.

[Sam Chak]

Hi @JUBAIR AHMED SHAMIM,

Firstly, thank you for your information. I have a question regarding the setup of boundary conditions. In the YouTube video, the individual named Abhishek Aravalli initiated the experiment by choosing the Adsorbent Material, configuring the Bed Height and Particle Size, and regulating the steady flow rate at which the gas flows through the column. Could you please explain the purpose or significance of setting up these boundary conditions?

Your insights would be greatly appreciated.

[Torsten]

@JUBAIR AHMED SHAMIM

I suggest you use the discretized equations below.

I don't know which boundary condition you want to use for Tbed at x = 0. If you tell me, I will try to deduce DTbedDT(1)

C1_water has a wrong unit ; it must be m/s instead of 1/s in your code.

Q_eq mustn't be constant as in your code - it must depend on T and Q.

I didn't check the other parameters of your code for correctness.

% Tgas and Q
for i = 1:n
  DTgasDT(i) = C_gas* (Tbed(i)-Tgas(i));
  DQDt(i) = K_LDF*(Q_eq - Q(i));
end
% Tbed
DTbedDT(1) = ?  % Which boundary condition do you want to set at x =0 ?
for i = 2:n-1
  DTbedDt(i)= C_cond*((Tbed(i+1)-Tbed(i))/(X(i+1)-X(i)) - (Tbed(i)-Tbed(i-1))/(X(i)-X(i-1)))/...
              ((X(i+1)+X(i))/2 - (X(i)+X(i-1))/2) + ...
              C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i)); % interior grid points
end
DTbedDt(n) = -C_cond*(Tbed(n)-Tbed(n-1))/((X(n)-X(n-1))*(X(n)-(X(n)+X(n-1))/2))+...
             C_bedA*(Tgas(i)-Tbed(i)) + C_bedB*(Tw(i)-Tbed(i))+ C_ads*K_LDF*(Q_eq-Q(i));  % dTw/dx = 0 at x = L
% Tw
DTwDt(1) = 0.0;   % Constant inlet temperature for water over time
for i = 2:n
  DTwDt(i) = -C1_water* (Tw(i)-Tw(i-1))/(X(i)-X(i-1)) + C2_water*(Tbed(i) - Tw(i)) ;   % interior and last grid point
end

[JUBAIR AHMED SHAMIM]

Hi Torsten, this is really great. You have saved my weeks of work. Thanks a lot. Now I understand how to set the BCs. For Tbed at x = 0 I want to use the same BC as at x=L, i.e., dTw/dx = 0 at x = 0. To do this, is it okay if I just replace n with 1 in the expression of DTbedDt(n)?

And thanks for the comment on C1 and Q_eq. Let me double check the unit of C1 and I will incorporate Langmuir isotherm for Q_eq, so that it will depend of temperature and pressure.

[Torsten]

DTbedDt(1) = C_cond*(Tbed(2)-Tbed(1))/((X(2)-X(1))*((X(1)+X(2))/2-X(1))) + ...
             C_bedA*(Tgas(1)-Tbed(1)) + C_bedB*(Tw(1)-Tbed(1))+ C_ads*K_LDF*(Q_eq-Q(1));  % dTw/dx = 0 at x = 0

And replace the "(i)" in the equation for DTbedDt(n) above by "(n)":

DTbedDt(n) = -C_cond*(Tbed(n)-Tbed(n-1))/((X(n)-X(n-1))*(X(n)-(X(n)+X(n-1))/2))+...
             C_bedA*(Tgas(n)-Tbed(n)) + C_bedB*(Tw(n)-Tbed(n))+ C_ads*K_LDF*(Q_eq-Q(n));  % dTw/dx = 0 at x = L

[JUBAIR AHMED SHAMIM]

Thanks a lot for all the efforts you have made. I will implement this BC and let you know the results in this week.

[JUBAIR]

Hi Toresten

I have one more question. Coould you please explain why C1_water is negative in the following code? I will highly appreciate your answer.

% Tw

DTwDt(1) = 0.0; % Constant inlet temperature for water over time

for i = 2:n

DTwDt(i) = -C1_water* (Tw(i)-Tw(i-1))/(X(i)-X(i-1)) + C2_water*(Tbed(i) - Tw(i)) ; % interior and last grid point

end