Packed bed storage modeling probelm

Hello, I have to make a numerical model of a packed bed storage full of spheres where I have to predict both fluid temperature and solid temperature across the bed length. Its a transient 1D coupled heat transfer problem where I know only the fluid inlet temperature. I am a newbie at matlab and I tried to develop code with the help of chatgpt for charging phase but it doesn't work. I also have to derive code for storage and discharging phase. Can you kindly help me with it?

rhoS = 2750; % Solid density [kg/m^3]

cpS = 850; % Solid specific heat capacity [J/(kg*K)]

ks = 1.28; % Effective thermal conductivity of solid [W/(m*K)]

epsilon = 0.36; % Void fraction

has = 47223.11473; % Heat transfer coefficient between solid and fluid [W/(m^2*K)]

Uw = 4.5118; % External heat transfer coefficient [W/(m^2*K)]

Tinfinity = 21 + 273.15; % Ambient temperature [K]

TinF = 83 + 273.15; % Inlet fluid temperature [K]

TinS = 21 + 273.15; % Initial solid temperature [K]

L = .14605; % Height of the packed bed [m]

u = 5.94; % Specific HTF mass flow rate [kg/(m^2*s)]

D = 0.1016; % Diameter of the bed [m]

dp = 0.01407; % Particle diameter [m]

rhoF = 1.08; % Fluid density [kg/m^3]

cpF = 1008; % Fluid specific heat [J/(kg*K)]

% Set up computational domain

Nx = 50; % Number of spatial steps

dx = L / Nx; % Spatial step size

dt = 1; % Time step size in seconds

Nt = 720; % Number of time steps

% Initialize temperature profiles

TS = ones(Nx, 1) * TinS; % Initial solid temperature profile

TF = ones(Nx, 1) * TinF; % Initial fluid temperature profile

TF(1) = Tinfinity

% Simulation loop

for n = 1:Nt

TS_new = TS;

TF_new = TF;

for i = 2:Nx-1

% Fluid temperature equation

TF_new(i) = TF(i) + dt * (-u * (TF(i) - TF(i-1)) / dx + (has/(rhoF * cpF * epsilon)) * (TS(i) - TF(i)));

% Solid temperature equation

d2TSdx2 = (TS(i+1) - 2*TS(i) + TS(i-1)) / dx^2;

TS_new(i) = TS(i) + dt * ((ks/(rhoS * cpS * (1 - epsilon))) * d2TSdx2 ...

+ (has/(rhoS * cpS * (1 - epsilon))) * (TF(i) - TS(i)) ...

+ (Uw * D * pi / (rhoS * cpS * (1 - epsilon))) * (Tinfinity - TS(i)));

end

% Update temperatures

TS = TS_new;

TF = TF_new;

% Apply boundary conditions

TF(1) = TinF; % Constant inlet temperature for fluid

TS(1) = TS(2); % Adiabatic boundary for solid at inlet

TS(Nx) = TS(Nx-1); %Adiabatic boundary for solid at outlet

end

% Plotting the results

x = linspace(0, L, Nx);

plot(x, TS-273.15, 'r', x, TF-273.15, 'b');

xlabel('Bed Length (m)');

ylabel('Temperature (°C)');

legend('Solid Temperature', 'Fluid Temperature');

title('Temperature Distribution Along the Bed Length');

grid on;

Answer 1

[Torsten]

rhoS = 2750;         % Solid density [kg/m^3]
cpS = 850;        % Solid specific heat capacity [J/(kg*K)]
ks = 1.28;           % Effective thermal conductivity of solid [W/(m*K)]
epsilon = 0.36;       % Void fraction
has = 47223.11473;             % Heat transfer coefficient between solid and fluid [W/(m^2*K)]
Uw = 4.5118;          % External heat transfer coefficient [W/(m^2*K)]
Tinfinity = 21 + 273.15; % Ambient temperature [K]
TinF = 83 + 273.15;  % Inlet fluid temperature [K]
TinS = 21 + 273.15;  % Initial solid temperature [K]
L = .14605;             % Height of the packed bed [m]
u = 5.94;           % Specific HTF mass flow rate [kg/(m^2*s)]
D = 0.1016;           % Diameter of the bed [m]
dp = 0.01407;           % Particle diameter [m]
rhoF = 1.08;         % Fluid density [kg/m^3]
cpF = 1008;          % Fluid specific heat [J/(kg*K)]

% Set up computational domain
Nx = 50;            % Number of spatial steps
dx = L / Nx;        % Spatial step size
dt = 0.0001;             % Time step size in seconds
Nt = 7200000;          % Number of time steps

% Initialize temperature profiles
TS = ones(Nx, 1) * TinS; % Initial solid temperature profile
TF = ones(Nx, 1) * TinF; % Initial fluid temperature profile
TF(1) = Tinfinity;

% Simulation loop
for n = 1:Nt
    TS_new = TS;
    TF_new = TF;
    for i = 2:Nx-1
        % Fluid temperature equation
        TF_new(i) = TF(i) + dt * (-u * (TF(i) - TF(i-1)) / dx + (has/(rhoF * cpF * epsilon)) * (TS(i) - TF(i)));
        
        % Solid temperature equation
        d2TSdx2 = (TS(i+1) - 2*TS(i) + TS(i-1)) / dx^2;
        TS_new(i) = TS(i) + dt * ((ks/(rhoS * cpS * (1 - epsilon))) * d2TSdx2 ...
            + (has/(rhoS * cpS * (1 - epsilon))) * (TF(i) - TS(i)) ...
            + (Uw * D * pi / (rhoS * cpS * (1 - epsilon))) * (Tinfinity - TS(i)));
    end
    % Fluid temperature equation (last point)
    TF_new(Nx) = TF(Nx) + dt * (-u * (TF(Nx) - TF(Nx-1)) / dx + (has/(rhoF * cpF * epsilon)) * (TS(Nx) - TF(Nx)));
    
    % Update temperatures
    TS = TS_new;
    TF = TF_new;
    
    % Apply boundary conditions
    TF(1) = TinF;      % Constant inlet temperature for fluid
    TS(1) = TS(2);     % Adiabatic boundary for solid at inlet
    TS(Nx) = TS(Nx-1); %Adiabatic boundary for solid at outlet
end 

% Plotting the results
x = linspace(0, L, Nx);
plot(x, TS-273.15, 'r', x, TF-273.15, 'b');
xlabel('Bed Length (m)');
ylabel('Temperature (°C)');
legend('Solid Temperature', 'Fluid Temperature');
title('Temperature Distribution Along the Bed Length');
grid on;

[soulhunter]

Thank you for the answer. But the thing is the result doesn’t reflect the real scenario. The fluid temperature should decrease by around 10-20 degrees and the solid will gain heat from the incoming fluid to rise from its initial temperature. I have checked the energy equation and it’s okay. I can’t find out why it is not reflecting the physical behaviour.

[Torsten]

Changing Nt to

Nt = 7200000;          % Number of time steps

gives you the above result (which should be almost steady-state).

[soulhunter]

Hi, the solid temperature should be rising from TinS = 21 + 273.15 due to hot fluid flow. Can you kindly explain why the solid temperature is begining from almost 83 degrees? Is there something wrong in my physics?

[Torsten]

This is a time-dependent simulation. Fluid and solid temperature equalize over time, and I plotted the temperatures after an "infinitly" long time.

If you let water of 83 degrees flow over a metal of 21 degrees, the metal will approach 83 degrees in the long term.