Packed bed storage modeling probelm

13 views (last 30 days)
soulhunter
soulhunter on 20 Apr 2024
Edited: Torsten on 20 Apr 2024
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;

Accepted Answer

Torsten
Torsten on 20 Apr 2024
Edited: Torsten on 20 Apr 2024
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;
  4 Comments
soulhunter
soulhunter on 20 Apr 2024
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
Torsten on 20 Apr 2024
Edited: Torsten on 20 Apr 2024
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Fluid Dynamics 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!