Packed-bed dynamic solution using Euler' s numerical method (explicit).

© 2026 Universidad de Salamanca, Department of Applied Physics, Thermodynamics, Statistical Physics and Energy Optimization Research Group . 
Author : David Pérez Gallego . 
ORCID : https : // orcid . org/0000 - 0003 - 2519 - 6014
Google Scholar : https : // scholar . google . com/citations?user = MuzFipMAAAAJ & hl = es

This work is licensed under the Creative Commons Attribution–NonCommercial–ShareAlike 4.0 International License (CC BY - NC - SA 4.0) . You are free to use, reproduce, and adapt this work for non - commercial purposes, provided that proper attribution is given to the author and the same license is maintained for any derivative works .

Initial Data

ClearAll;
SetDirectory[NotebookDirectory[]];
Off[CreateDirectory::eexist];

Directorio="Directory_Name";   


tfc=300;(*Time*)

materialLecho="Aluminio";(* Aluminio ,   Arena,   PCM_NaNO2  ,    Coarser_GR50 *)
fluid="Air";(*  PCM_60NaNO3_40KNO3 , Air ,   Therminol66*)

L=1;
DiametroLecho=0.7;(*m*)

mdot=0.5;(*Kg/s*) 



Δz=0.01;(*Δz*)
Δt=1;
int=10;(*pasos*)







Tfinlet=400+273;(*Kelvin*)

Tfinitial=25+273;(*Tª inicial del fluido.*)
Tsinitial=25+273;(*Tª inicial de la roca.  *)
Tambiente=25+273;(*Tª ambiente del exterior del tanque*)

Aislante="Fibra de vidrio";
Espesor=0.2;(*de la capa aislante*)




Pambiente=1;(*atm*)(*Presión ambiente*)






DirCarga=Directorio<>"/Carga";

CreateDirectory[{Directorio,DirCarga}];

Physical properties

Clear[T,P,ρs,ρf,Ks,Kf,cpf,cps,h,Us,u,Pr,Reyd,Nu,μf,ν];

If[fluid=="Air",
Clear[Kf];Kf[T_]:=0.05;
Clear[ρf];ρf[T_]:=0.8;
Clear[μf];μf[T_]:=0.0000185;
Clear[cpf];cpf[T_]:=1005;

];




(******************************)

Nz=IntegerPart[L/Δz]+1;(*Número de divisiones espaciales*)




(*********************************************************)
(*********      MATERIALES PARA EL LECHO       **************)
(*********************************************************)


If[materialLecho=="Aluminio",
(*Lecho: Aluminio   Cascetta 2016*)
ρs=3550;(*(kg/(m^3))  densdiad  lecho*)
cps=902;(*(J/kg K)capacidad calorífica  lecho*)
Ks=20;(*W/mK conductividad térmica de las bolas de aluminio de Cascetta*)
ϵ=0.39;(*porosidad*)
dp=0.008;(*(m)  Diámetro de las partículas del lecho*)
];

If[materialLecho=="Arena",
(*Lecho: Arena*)
ρs=3550;(*(kg/(m^3))  densdiad roca lecho:  es arena*)
cps=902;(*(J/kg K)capacidad calorífica roca lecho: arena*)
Ks=0.5;(*conductividad térmica de las bolas de aluminio de Cascetta*)
ϵ=0.39;(*porosidad*)
dp=0.008;(*(m)  Diámetro de las partículas del lecho*)

];




Q[T_]:=(mdot/(ρf[T]));(*flujo volumétrico ((m^3)/s)  *)
Us[T_]:=Q[T]/(π*(((DiametroLecho/2))^2));(*velocidad superficial fluido*)
u[T_]:=Us[T]/ϵ;(*velocidad intersticial fluido*)



ν[T_]:=μf[T]/ρf[T];(*viscosidad cinemática del fluido*)




(*Tubería: es acero*)
ρw=7900;(*(kg/(m^3)) *)
cpw=477;(*(J/kg K)*)
Kw=14.9;(*conductividad térmica*)

(*Aisalente: es fibra de vidrio*)
ρins=32;(*(kg/(m^3)) *)
cpins=835;(*(J/kg K)*)
Kins=0.038;(*conductividad térmica aislante*)



(*********************************************************)
(*********      PARÁMETROS DEPENDIENTES       **************)
(*********************************************************)





ab=((2(DiametroLecho/2+Espesor))/(((DiametroLecho/2))^2));



difusividadTermica[T_]:=Kf[T]/(cpf[T]*ρf[T]);(*difusividad termica fluido*)
Pr[T_]:=ν[T]/difusividadTermica[T];(*Nº de Prandtl*)



Reyd[T_]:=(Us[T]*dp)/(ν[T]);(*Nº de Reynolds*)
(*Nº de Nusslet*)
Nu[T_]:=2.0+1.1*((Reyd[T])^0.6)*((Pr[T])^(1/3));
(*coeficiente transferencia de calor convección entre roca y fluido: *)
h[T_]:=(Kf[T]*Nu[T])/dp;


(*Coeficiente de convección fluido-pared del lecho*)
Nuw[T_]:=2.57*((Reyd[T])^(1/3))*((Pr[T])^(1/3))+0.0936((Reyd[T])^0.8)*((Pr[T])^0.4);(*Nº de Nusslet para pared de acero*)

hw[T_]:=(Nuw[T]*Kf[T])/dp;(*coef. transferencia de calor convección entre pared y fluido*)

(*Coeficiente de convección tubería-aire exterior*)
GrL[T_]:=(9.81*(1/24)*(5)*(L^3))/((ν[T])^2);(*Nº de Grashof; 1/24 es el coef. de expansión térmica*)
RaL[T_]:=GrL[T]*Pr[T];(*Nº de Rayleigh*)
Nuins[T_]:=0.68+(0.67*((RaL[T])^(1/4)))/(((1+(((0.492/Pr[T]))^(9/16))))^(4/9));(*Nº de Nusslet tubería-aire*)
hins[T_]:=(Nuins[T]*Kf[T])/L;(*coef. transferencia de calor entre lecho y aire*)
(*Según Kocak y también Ismail: *)
(*hx[T_]:=1.42(((((T-Tambiente)/L)))^0.25);*)
hx[T_]:=1.42((((75/L)))^0.25);(*Ismail*)


(*Coefiente final de transferenica de calor por radiacción:*)




Clear[U];U[T_]:=1/((1/(π*(((DiametroLecho/2))^2) hw[T]))+((Log[((DiametroLecho+Espesor)/DiametroLecho)])/(2 π L Kins))+1/(π*((((DiametroLecho+Espesor)/2))^2)hx[T]));






as=(6(1-ϵ))/dp(*(Subscript[A, s,total])/(Subscript[V, b])   *);












(*Ficheros de la carga*)
rutaCargaFluido=DirCarga<>"/"<>"Carga_Fluido.txt";
rutaCargaRoca=DirCarga<>"/"<>"Carga_Roca.txt";

Start: computing

AbsoluteTiming[
Clear[Nt];

Nt=IntegerPart[(((tfc))/(Δt))]+1;(*Número de divisiones temporales*)


(*TEMPERATURAS:   T[[tiempo,posición]]*)
Clear[Tf,Ts];Tf=ConstantArray[0,{2,Nz}];Ts=ConstantArray[0,{2,Nz}];


(*Condiciones iniciales*)
Tf[[All,All]]=Tfinitial;(*condición inicial fluido*)
Ts[[All,All]]=Tsinitial;(*condición inicial roca*)
Tf[[All,1]]=Tfinlet;
Clear[Ts100,Tf100];Ts100={};Tf100={};(*Matriz con la solución de las temperaturas *)

Ts100=Append[Ts100,Ts[[1]]];
Tf100=Append[Tf100,Tf[[1]]];
j1=0;(*Para la barra de carga*)
Do[
j1=j1+1;
(*Roca termaliza con frontera del fluido (z=1)*)
Clear[y];

(*
Ts[[2,1]]=y/.Flatten[NSolve[y==Ts[[1,1]]+((Δt)/((1-ϵ)ρs cps))(Ks*((Ts[[1,2]]-2Ts[[1,1]]+Ts[[1,1]])/((Δz)^2))+h[Tf[[1,1]]] *as*(Tf[[2,1]]-y)),y,Reals]];
*)
Ts[[2,1]]=Tf[[2,1]];

For[i=2,i<Nz,i++,

(*PROGRAMA NSOLVE*)
(*Clear[x,y];

{Tf[[2,i]],Ts[[2,i]]}={x,y}/.Flatten[NSolve[x==Tf[[1,i]]+Δt(((cpf[Tf[[1,i]]]*Us[Tf[[1,i]]])/(cpf[Tf[[1,i]]]*Δz* ϵ))(Tf[[2,i-1]]-x)+((Kf[Tf[[1,i]]])/(ϵ*ρf[Tf[[1,i]]]*cpf[Tf[[1,i]]]*((Δz)^2)))*(Tf[[1,i-1]]-2Tf[[1,i]]+Tf[[1,i+1]])+((h[Tf[[1,i]]] as *(y-x))/(ϵ*ρf[Tf[[1,i]]]*cpf[Tf[[1,i]]]))-((U[Tf[[1,i]]]*ab)/(ϵ*ρf[Tf[[1,i]]]*cpf[Tf[[1,i]]]))(x-Tambiente))
&&
y==Ts[[1,i]]+((Δt)/((1-ϵ)ρs cps))(Ks*((Ts[[1,i+1]]-2Ts[[1,i]]+Ts[[1,i-1]])/((Δz)^2))+h[Tf[[1,i]]] *as*(x-y)),{x,y},Reals]];

*)



(*
SOLUCIONES DESPEJADAS:
*)

(*X*)
Tf[[2,i]]=-((Δt (-cps (-1+ϵ) ρs+as Δt h[Tf[[1,i]]]) Kf[Tf[[1,i]]] (Tf[[1,-1+i]]-2 Tf[[1,i]]+Tf[[1,1+i]])+as Δt h[Tf[[1,i]]] (Ks Δt Ts[[1,-1+i]]+(-2 Ks Δt-cps ((Δz)^2) (-1+ϵ) ρs) Ts[[1,i]]+Ks Δt Ts[[1,1+i]])+ab Tambiente Δt ((Δz)^2) (-cps (-1+ϵ) ρs+as Δt h[Tf[[1,i]]]) U[Tf[[1,i]]]+Δz cpf[Tf[[1,i]]] (-cps (-1+ϵ) ρs+as Δt h[Tf[[1,i]]]) (Δz ϵ Tf[[1,i]]+Δt Tf[[2,-1+i]] Us[Tf[[1,i]]]) ρf[Tf[[1,i]]])/(Δz (cps (-1+ϵ) ρs (ab Δt Δz U[Tf[[1,i]]]+cpf[Tf[[1,i]]] (Δz ϵ+Δt Us[Tf[[1,i]]]) ρf[Tf[[1,i]]])-as Δt h[Tf[[1,i]]] (-cps Δz (-1+ϵ) ρs+ab Δt Δz U[Tf[[1,i]]]+cpf[Tf[[1,i]]] (Δz ϵ+Δt Us[Tf[[1,i]]]) ρf[Tf[[1,i]]]))));

(*Y*)
Ts[[2,i]]=((Ks Δt Ts[[1,-1+i]]+(-2 Ks Δt-cps ((Δz)^2) (-1+ϵ) ρs) Ts[[1,i]]+Ks Δt Ts[[1,1+i]]) (ab Δt Δz U[Tf[[1,i]]]+cpf[Tf[[1,i]]] (Δz ϵ+Δt Us[Tf[[1,i]]]) ρf[Tf[[1,i]]])+as Δt Δz h[Tf[[1,i]]] (Δt Kf[Tf[[1,i]]] (Tf[[1,-1+i]]-2 Tf[[1,i]]+Tf[[1,1+i]])+Ks Δt Ts[[1,-1+i]]-2 Ks Δt Ts[[1,i]]+cps ((Δz)^2) ρs Ts[[1,i]]-cps ((Δz)^2) ϵ ρs Ts[[1,i]]+Ks Δt Ts[[1,1+i]]+ab Tambiente Δt ((Δz)^2) U[Tf[[1,i]]]+Δz cpf[Tf[[1,i]]] (Δz ϵ Tf[[1,i]]+Δt Tf[[2,-1+i]] Us[Tf[[1,i]]]) ρf[Tf[[1,i]]]))/(((Δz)^2) (as Δt h[Tf[[1,i]]] (ab Δt Δz U[Tf[[1,i]]]+cpf[Tf[[1,i]]] (Δz ϵ+Δt Us[Tf[[1,i]]]) ρf[Tf[[1,i]]])-cps (-1+ϵ) ρs (as Δt Δz h[Tf[[1,i]]]+ab Δt Δz U[Tf[[1,i]]]+cpf[Tf[[1,i]]] (Δz ϵ+Δt Us[Tf[[1,i]]]) ρf[Tf[[1,i]]])));




];

(*Termalizar Z=L*)
Clear[x,y];


Tf[[2,Nz]]=Tf[[2,Nz-1]];
Ts[[2,Nz]]=Ts[[2,Nz-1]] ; (* Derivadas cero en la frontera del fluido y roca *)


If[Mod[t,int]==0,Ts100=Append[Ts100,Ts[[2,All]]]];
If[Mod[t,int]==0,Tf100=Append[Tf100,Tf[[2,All]]]];

Tf[[1]]=Tf[[2]];
Ts[[1]]=Ts[[2]];


(*barra de carga*)If[Mod[t,10]==0,
NotebookDelete[nota1];
NotebookDelete[nota2];
nota1=PrintTemporary["Carga: ",N[(j1/Nt)*100,4]," %"];
nota2=PrintTemporary[ProgressIndicator[j1,{1,Nt}]];
];

,{t,1,Nt}];

(*Aquí exportamos*)
Export[rutaCargaFluido,Tf100[[All,All]],"Table"];
Export[rutaCargaRoca,Ts100[[All,All]],"Table"];


Print["Tiempo total de cálculo: "];
][[1]]/60

Results: temperature plots

Clear[TemperaturasFluido,TemperaturasRoca];
TemperaturasFluido=Import[rutaCargaFluido,"Table"];
TemperaturasRoca=Import[rutaCargaRoca,"Table"];





Manipulate[
ListPlot[
{
Table[{(i-1)*Δz,TemperaturasFluido[[t,i]]},{i,1,Dimensions[TemperaturasFluido][[2]]}]
,Table[{(i-1)*Δz,TemperaturasRoca[[t,i]]},{i,1,Dimensions[TemperaturasRoca][[2]]}]
}
,PlotRange->{Tfinitial,Tfinlet},Joined->True,AxesLabel->{"z [m]","T [K]"},PlotStyle->{Red,Blue}]


,{t,1,Length[TemperaturasRoca],1}]