1+ component (Propagation = blocks) lumpedTES
2+ % TES:2.5
3+ % This component models a thermal energy storage (TES) or a sand battery which
4+ % stores thermal energy in form of sensible heat of the sand.
5+ % This battery is charged by electrical nodes, and discharged by moist air nodes.
6+ % The sand is heated when electric current flows through the resistive heater.
7+ % For removing the heat out from the sand battery, air must be blown through the metallic pipe terminals at node A and B.
8+
9+ % Copyright 2023-24 The MathWorks, Inc.
10+
11+ nodes
12+ A = foundation.moist_air.moist_air;
13+ n = foundation.electrical.electrical;%-:left
14+ B = foundation.moist_air.moist_air;
15+ p = foundation.electrical.electrical;%+:right
16+ H = foundation.thermal.thermal;
17+ end
18+ annotations
19+ A:Side = left; % Air flow inlet port
20+ B:Side = right;% Air flow outlet port
21+ H:Side = bottom;% Thermal port
22+ end
23+
24+ annotations
25+ UILayout = [UIGroup('Layout',coilShape,numCoilX,numCoilY,cellRefT)...
26+ UIGroup('Heater Element',coilWidthHeater,coilDiameterHeater,coilDepthHeater,nturnsHeater,heaterTubeDia,heaterTimeconst,heaterRes,heaterEmissivity,heaterMaxT)...
27+ UIGroup('Air Duct',coilWidthCooler,coilDiameterCooler,coilDepthCooler,nturnsCooler,pipeHydraulicDia)...
28+ UIGroup('Sand Properties',sandDensity,sandThCond,sandSpHeat,sandInitT)];
29+ Icon = 'ThermalBattery.svg';
30+ end
31+
32+ parameters
33+ coilShape = layoutShape.HelicalCoil; % Coil shape
34+ numCoilX = {1,'1'}; % Number of coil rows
35+ numCoilY = {1,'1'}; % Number of coils in each row
36+ coilDepthHeater = {0.1, 'm'}; % Depth of heating coil
37+ nturnsHeater = {2,'1'}; % Number of heating coil turns
38+ coilDepthCooler = {0.1, 'm'}; % Depth of cooling coil
39+ nturnsCooler = {2,'1'}; % Number of cooling coil turns
40+ sandDensity = {1500, 'kg/m^3'}; % Density of sand
41+ sandThCond = {1.2, 'W/(K*m)'}; % Thermal conductivity of sand
42+ sandSpHeat = {850, 'J/(K*kg)'}; % Specific heat of the sand
43+ cellRefT = {298.15, 'K'}; % Reference temperature for heating potential
44+ sandInitT = {298.15, 'K'}; % Initial sand mass temperature
45+ pipeHydraulicDia = {0.01, 'm'}; % Hydraulic diameter of air pipe
46+ heaterTubeDia = {0.01, 'm'}; % Diameter of heater tube
47+ heaterEmissivity = {0.8 , '1'}; % Heater tube emissivity
48+ heaterTimeconst = {0.2, 's'}; % Heater filament thermal time constant
49+ heaterRes = {1, 'Ohm'}; % Heater filament electrical resistance
50+ heaterMaxT = {750, 'degC'}; % Heater filament temperature upper limit
51+ end
52+ parameters (ExternalAccess = none)
53+ stfBtz = {5.670374419e-8, 'W/(K^4*m^2)'};
54+ numTube = numCoilX * numCoilY; % total no of cooling or heating tube
55+ hotSandMass = numTube * hotSandVol * sandDensity; % mass of the sand associated with heating element
56+ coldSandMass = numTube * coldSandVol * sandDensity; % mass of the sand associated with cooling element
57+ pipeArea = numTube * pi * pipeHydraulicDia * coolerLen; % area for heat conduction B/W sand and air pipe
58+ eqHydrlcD = numTube * pipeHydraulicDia; % equivalent pipe hydraulic diameter
59+ heaterArea = numTube * pi * heaterTubeDia * heaterLen; % area for heat conduction B/W sand and electrical resistance
60+ coldSandEquivalentDia = (biotFactor * pipeHydraulicDia + coilWidthCooler) / 2; % equivalent OD of cold sand volume
61+ radiationCoefficient = stfBtz * heaterEmissivity;
62+ coilWidthHeater = {0.1, 'm'}; % Width of heating coil
63+ coilDiameterHeater = {0.1, 'm'}; % Diameter of heating coil
64+ coilWidthCooler = {0.1, 'm'}; % Width of cooling coil
65+ coilDiameterCooler = {0.1, 'm'}; % Diameter of cooling coil
66+ end
67+
68+ parameters (ExternalAccess = none)
69+ biotFactor = {7.04, '1'}; % ratio of sand column diameter to cooling pipe diameter to maintain Biot Number = 0.1
70+ end
71+
72+ if coilShape == layoutShape.UshapedCoil
73+ annotations
74+ [coilWidthHeater] : ExternalAccess = modify
75+ [coilWidthCooler] : ExternalAccess = modify
76+ end
77+ parameters (ExternalAccess = none)
78+ heaterLen = nturnsHeater * coilWidthHeater + (nturnsHeater - 1) * (pi / 2) * coilDepthHeater / nturnsHeater; % length of heater element for S-type option
79+ coolerLen = nturnsCooler * coilWidthCooler + (nturnsCooler - 1) * (pi / 2) * coilDepthCooler / nturnsCooler; % length of cooler element for S-type option
80+ sandConductionAr = (biotFactor / 4) * (coilDepthCooler + coilDepthHeater) * (pipeHydraulicDia + heaterTubeDia); % heat conduction area B/W hot and cold sand
81+ sandConductionL = (coilWidthCooler + coilWidthHeater) / 2;
82+ hotSandVol = (biotFactor * coilWidthHeater * heaterTubeDia * coilDepthHeater) - (heaterLen * pi * (heaterTubeDia ^ 2) / 4);
83+ coldSandVol = (biotFactor * coilWidthCooler * pipeHydraulicDia * coilDepthCooler) - (coolerLen * pi * (pipeHydraulicDia ^ 2) / 4);
84+ end
85+ else
86+ annotations
87+ [coilDiameterHeater] : ExternalAccess = modify
88+ [coilDiameterCooler] : ExternalAccess = modify
89+ end
90+ parameters (ExternalAccess = none)
91+ heaterLen = sqrt((2 * pi * coilDiameterHeater) ^ 2 + (coilDepthHeater / nturnsHeater) ^ 2); % length of heater element for Helix option
92+ coolerLen = sqrt((2 * pi * coilDiameterCooler) ^ 2 + (coilDepthCooler / nturnsCooler) ^ 2); % length of cooler element for Helix option
93+ sandConductionAr = (coilDepthCooler + coilDepthHeater) * (coilDiameterCooler + coilDiameterHeater) / 4; % heat conduction area B/W hot and cold sand
94+ sandConductionL = (coilDiameterCooler + coilWidthHeater) / 2;
95+ hotSandVol = (coilWidthHeater ^ 2 * coilDepthHeater) - (heaterLen * pi * (heaterTubeDia ^ 2) / 4);
96+ coldSandVol = (coilDiameterCooler ^ 2 * coilDepthCooler) - (coolerLen * pi * (pipeHydraulicDia ^ 2) / 4);
97+ end
98+ end
99+
100+
101+
102+ equations
103+ assert(heaterLen > 0, 'Heating element length must be greater than 0')
104+ assert(coolerLen > 0, 'Air pipe length must be greater than 0')
105+ assert(sandConductionAr > 0)
106+ assert(sandConductionL > 0)
107+ assert(hotSandVol > 0, 'Calculated sand volume is negative or zero for given inputs')
108+ assert(coldSandVol > 0, 'Calculated sand volume is negative or zero for given inputs')
109+ assert(hotSandMass > 0, 'Calculated sand mass is negative or zero for given inputs')
110+ assert(coldSandMass > 0, 'Calculated sand mass is negative or zero for given inputs')
111+ assert(coldSandEquivalentDia > 0)
112+ assert(numCoilX >= 1, 'Number of coil rows must be a positive integer')
113+ assert(numCoilY >= 1, 'Number of coils must be a positive integer')
114+ end
115+
116+ variables (ExternalAccess = observe)
117+ hPotential = {value = sandSpHeat * hotSandMass * (sandInitT - cellRefT) + sandSpHeat * coldSandMass * (sandInitT - cellRefT), priority = priority.high}; % battery sensible heat potential
118+ end
119+
120+ variables
121+ filamentTemp = {298.15, 'K'}; % Heater filament temperature
122+ end
123+
124+ components(ExternalAccess = observe)
125+ hotSand = foundation.thermal.elements.mass(mass = hotSandMass, num_ports = int32(1), sp_heat = sandSpHeat, T = {value = sandInitT, priority = priority.high});
126+ coldSand = foundation.thermal.elements.mass(mass = coldSandMass, num_ports = int32(1), sp_heat = sandSpHeat, T = {value = sandInitT, priority = priority.high});
127+ heater = foundation.electrical.elements.thermal_resistor(K_d = {800, 'W/K'}, R0 = heaterRes / numTube, T0 = {298.15, 'K'}, alpha = {5e-05, '1/K'}, tc = heaterTimeconst, T = {value = sandInitT, priority = priority.high});
128+ HTradiate = foundation.thermal.elements.radiation(area = numTube * heaterArea, rad_tr_coeff = radiationCoefficient);
129+ airDuct = foundation.moist_air.elements.pipe(Dh = eqHydrlcD, Nu_lam = 3.66, RH_ws = 1, Re_lam = 2000, Re_tur = 4000, area = pipeArea, length = coolerLen, length_add = coolerLen / 5,...
130+ moisture_trace_gas_source = int32(0), roughness = {1.5e-05, 'm'}, shape_factor = 64, wall_condensation = int32(0));
131+ HTcond1 = foundation.thermal.elements.conduction(L = coolerLen, d_in = pipeHydraulicDia, d_out = coldSandEquivalentDia, th_cond = sandThCond, thermal_type = int32(1), wall_geometry = int32(2));
132+ HTcond2 = foundation.thermal.elements.conduction(area = sandConductionAr, th_cond = sandThCond, thermal_type = int32(1), thickness = sandConductionL, wall_geometry = int32(1));
133+ end
134+
135+ intermediates (ExternalAccess = none)
136+ hPot = sandSpHeat * hotSandMass * (hotSand.M.T - cellRefT) + sandSpHeat * coldSandMass * (coldSand.M.T - cellRefT);
137+ heaterCoilT = heater.H.T;
138+ end
139+
140+ equations
141+ hPotential == hPot;
142+ filamentTemp == heaterCoilT;
143+ end
144+
145+ connections
146+ connect(n, heater.n);
147+ connect(HTradiate.A, heater.H);
148+ connect(B, airDuct.B);
149+ connect(A, airDuct.A);
150+ connect(HTcond1.B, airDuct.H);
151+ connect(H, HTcond2.B);
152+ connect(H, HTcond1.A);
153+ connect(H, coldSand.M);
154+ connect(p, heater.p);
155+ connect(HTcond2.A, HTradiate.B);
156+ connect(HTcond2.A, hotSand.M);
157+ end
158+ end
0 commit comments