Experimental Analysis and Numerical Modeling of a Shell and Tube

Subscriber access provided by La Trobe University Library

Article

Experimental analysis and numerical modeling of a shell and tube heat storage unit with phase change materials Tilman Johannes Barz, Christoph Zauner, Daniel Lager, Diana C. López C., Florian Hengstberger, Mariano Nicolás Cruz Bournazou, and Klemens Marx Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b01080 • Publication Date (Web): 01 Jul 2016 Downloaded from http://pubs.acs.org on July 6, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Experimental analysis and numerical modeling of a shell and tube heat storage unit with phase change materials Tilman Barz,

∗, †

Christoph Zauner,

Hengstberger,

†

†

Daniel Lager,

†

Diana C. López C.,

Mariano Nicolás Cruz Bournazou,

¶

‡

Florian

and Klemens Marx

†

†AIT Austrian Institute of Technology GmbH, Giengasse 2, 1210 Vienna, Austria ‡Technische Universität Berlin, Chair of Process Dynamics and Operation, Sekr. KWT-9,

Str. des 17. Juni 135, 10623 Berlin, Germany ¶Technische Universität Berlin, Institute of Biotechnology, Department of Bioprocess

Engineering, Ackerstr. 71-76, D-13355 Berlin, Germany E-mail: [email protected]

Abstract Thermal storages are part of highly integrated energy systems. The development of accurate and reduced models is critical for ecient simulations on a system-level and the analysis of the storage design, control and integration. We present the experimental analysis and numerical modeling of a lab-scale shell and tube latent heat thermal energy storage (LHTES) unit with a (latent) storage capacity of about 10-15 kWh. The phase change material (PCM) is a high density polyethylene (HD-PE) with phase change temperatures between 120◦ C and 135◦ C. An ecient 2D numeric storage model is derived which accounts for design and material parameters of PCM, storage and heat transfer uid (HTF). Dierent probability distribution functions are used to model 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the PCM apparent specic heat capacity. From these functions the state of charge (SOC) can be predicted, which indicates the extend to which a LHTES is charged relative to storeable latent heat. Model predictions are tted to experimental data from thermophysical measurements and from LHTES operation with partial and full charging/ discharging. The storage model agrees with experimental results. However, thermophysical material analysis and storage operation revealed that phase transition is noticeably aected by applied heating and cooling rates. It is shown how the neglection of these phenomena leads to limitations in the predictive quality of storage internal PCM temperatures.

Introduction Thermal energy storages (TES) using phase change materials (PCMs) have been subject to extensive research in the last three decades. They oer a signicantly higher energy density compared to sensible heat storage systems. A multiplicity of applications is reported, e.g. 5

in systems for ice storage, conservation and transport of temperature sensitive materials, building applications, solar heating and cooling, thermal control applications in buildings and spacecrafts, o-peak electricity storage systems, waste heat recovery systems, and others. 13 However, most of the phase change problems consider temperature ranges between 0 ◦ C and 60◦ C suitable for domestic heating applications. 4 A review of applicable materials can be

10

found elsewhere. 1,3 Most PCMs suer from low thermal conductivity, hence heat transfer enhancement techniques are required to eectively charge and discharge LHTES. 1,3,4 An overview of numerical and experimental studies to assess the eects of dierent geometries and congurations of PCM containers, e.g. at plate, cylinder, sphere, nned tube systems, dispersion of PCM

15

with high conductivity particles, and other can be found elsewhere. 1,4 Besides these studies, the inuence of operating parameters, e.g. the inlet temperature and the mass ow rate of the heat transfer uid (HTF), 4 as well as the inuence of specic thermophysical and trans2


Page 2 of 55

Page 3 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


port phenomena, e.g. natural convection in the liquid PCM or temperature dependence of liquid and solid density, have been subject of interest. 1

20

Modeling of LHTES

Dierent reviews on the modeling and numerical solution of heat

transfer problems in LHTES have been published. They are organized according to the problem geometry and specic conguration or application, 5 the mathematical problem formulation as well as analytical and numerical solution method 6 or the availability and use of experimental data for validation. 7 Besides analytical methods, various numerical methods 25

have been developed for a better handling of the multidimensional phenomena involved in PCM. They involve either the description of the pure heat conduction problem with a phase change, or a combined description of the convection/diusion phase change problem. 6,8 The models also dier on how to deal with the moving boundary between phases during melting and solidication. 6 Some of the models aim for calculation of the exact location of the moving

30

solid/liquid interface as a function of time and space (strong formulation). Others apply a set of equations which applies over the entire region of the domain where the phase change takes place (weak formulation). 6,8,9 For the rst approach (strong formulation) dierent numeric implementations exist applying either a xed or variable space or time grid. Popular examples using the second approach (weak formulation) are the enthalpy method, and specically

35

the so called apparent heat capacity method. 10 Here the governing equation is similar to a single phase equation without explicit treatment of the moving interface. Instead a mushy transition zone between the two phases is considered where a single enthalpy- or apparent heat capacity-temperature curve is applied. 3,6,11,12 This is of interest for commercial-grade materials and mixtures with phase change temperature ranges which can reach an order of

40

10K or more. 12 Signicant density changes during phase change are common for most PCMs. While there is only a minor eect on the heat transfer, 12 a volume decrease may lead to creation of voids and separation of PCM from the heat transfer surface. These voids aect the

3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

heat transfer and thus the performance of the LHTES considerably. 1,12,13 Inclusion of the 45

temperature-density dependence leads to a problem with moving boundaries and an internal material velocity due to density changes. 14 A solution can be found by dening a new space coordinate xed to the material. For simplicity a constant density is often assumed. 8 Besides of the density, also the thermal conductivity and heat capacity depend on phase state and temperature. The weak problem formulation allows for direct implementation of

50

this dependency 11,15 where the Kirchho transformation can be applied to obtain a reduced simple form. 14 Another phenomenon is convection within the PCM arising temperature and/or concentration gradients in the liquid phase under the action of buoyancy forces. 6 This ow, transporting energy, can aect the melting considerably and change the shape of the melt-

55

ing front. 12 The Boussinesq approximation for laminar ows can be applied to treat this convective ow. 8 However, convection inside PCM modules is often disregarded to simplify the model. The resulting heat conduction controlled phase change problem is valid for PCM with a high viscosity which applies to most PCM materials discussed for LHTES systems. Hysteresis loops of the enthalpy/temperature in thermally induced phase transformation

60

have been experimentally identied using dierential scanning calorimetry (DSC). 16,17 Hysteresis is normally not considered in the modeling of LHTES. 12 A simple hysteresis model uses a shift of the enthalpy curve according to a dierential temperature resulting in a delay of the phase change. 15 Although practically relevant, kinetics of the phase transition are normally not considered

65

in the modeling of LHTES. Studies on the isothermal melt-and cold-crystallization kinetics of syndiotactic polypropylene using DSC and macrokinetic modeling is described in Supaphol and Spruiell 18 . Lastly, supercooling (and less common superheating) is a phenomenon which is normally not considered in the modeling of LHTES. 12 Supercooling is caused mainly by primary

70

nucleation and crystal growth mechanisms. Subcooling produces the same eects in the

4


Page 4 of 55

Page 5 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperature as supercooling. Subcooling can be caused by kinetics in the solidication, i.e. when cooling rates and heat transfer in the PCM are high compared to the crystallization kinetics. For the operation of LHTES subcooling can be highly problematic when it is in the order of magnitude of the driving temperature dierence while discharging the storage, 75

i.e. solidication. 12 A simple model considering subcooling is presented elsewhere. 15

Optimization and control of TES systems

The inherent transient behavior of TES

systems requires eective design and operation strategies considering the process dynamics. 1921 A review with focus on applications in combined heat and power systems, building systems, and solar thermal power systems can be found elsewhere. 19,22 However, TES sys80

tems are frequently designed in view of quasi-static energy duties, adopting widely used eciency metrics commonly used for heat exchangers, i.e. NTU method. 23,24 Limitations exist regarding the extrapolation of the correlations to dierent storage congurations, geometries and operating parameters as well as the accuracy of these methods. 25 Dynamic optimization and advanced control approaches based on detailed LHTES models have been

85

used recently in applications of thermal energy management. 20,26 However, most models are not validated and base on elaborate computational uid dynamics and multi-physics tools leading to extremely time-consuming optimization procedures. 25 Accordingly, they are suitable for studying the detailed unit operation, rather than for studying devices on a system level.

90

It is important to develop system-level models which allow for the analysis of the interactions of TES with all other components of the system. 19 This contribution deals with the development and experimental validation of a reduced shell and tube LHTES model. The model can be eciently solved yielding suciently accurate predictions. It considers storage geometry, thermophysical data and corresponding correlations and can therefore be adapted

95

to dierent designs and PCMs. The paper is organized as follows: The experimental set-up of a shell and tube laboratory scale LHTES and methods for the determination of material

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

properties are briey presented in section 'Experimental'. A reduced two-dimensional model for the LHTES is derived in section 'Modeling'. Numeric solution methods for the simulation and model validation are discussed in section 'Numeric solution'. The tting and selection 100

of storage models and applicable thermophysical property correlations as well as an error analysis of estimated parameters are presented in section 'Fitting models and parameter estimation'. The results for the predicted stored energy and state of charge are discussed in section 'Predicted stored energy and state of charge'. Finally, a summary is given.

Experimental 105

Latent heat thermal energy storage Several designs of LHTES systems have been proposed. 4 The shell and tube design, with the PCM at the shell side and the HTF owing through the pipes in the center, is claimed to have a benecial performance with respect to charging and discharging of such systems. 2730 It is reported that a two-dimensional heat transfer problem is dominant in both, the pipe

110

and the shell where the heat transfer in the PCM is limiting the performance of an LHTES. 4 Improvements can include ns (circular or longitudinal) and the introduction of a matrix into the PCM with a high thermal conductivity, e.g. metal matrix 31 or carbon bers. 28 Due to its simplicity mostly nned systems are used. This evens out temperature gradients in radial direction while temperature gradients in axial direction increases indicating an improved

115

heat transfer. 32 A shell and tube LHTES with circular ns is considered, see Fig. 1. The system is composed of in total 72 steel tubes of each 2.5 m in length. To improve the heat transfer there are 312 aluminum ns in perpendicular direction to each tube installed.

6


Page 6 of 55

Page 7 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


port at he nsfer tra d i flu

fins pipes 4

port

3

cross-sec�on

2

at he nsfer tra d i flu

1

Figure 1: Design sketch of the shell and tube latent heat thermal energy storage (LHTES). Temperatures of the PCM are measured at axial positions 1, 2, 3, 4. R HD6070EA, available from INEOS Polyolens. This is a The storage material is Rigidex 120

high density polyethylene (HD-PE) with a narrow molecular weight distribution, suitable for a wide range of injection and compression molding applications. A comprehensive summary of design parameters is given in Table 1. Thermal energy for heating and cooling is supplied by a controllable laboratory heat transfer system ITHW 350 / 50kW from LAUDA DR. R. WOBSER. The measurement setup includes ow and temperature measurements. The ow

125

rate of the HTF is measured by ultrasonic ow measurement using the FLEXIM UMFLUXUS F6V. A temperature compensation is applied to account for the change in physical properties of the HTF during the ow measurement. The temperature of the HTF is determined at the inlet and at the outlet of the LHTES using platinum thermometers. Temperatures of the PCM are analyzed over the length of the storage at the axial positions 0.10, 0.87, 1.64

130

and 2.40 m. At each axial position 4 individual thermocouple sensors are installed exactly between two ns. Mean values are considered for further analysis. In the following we refer to these measurements as: HTF inlet, HTF outlet and PCM 1, 2, 3, 4.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 55

Table 1: Overview of LHTES design parameters Item storage casing material casing length inner cross-section inner cross-section area HTF type approximate HTF hold-up tube material number of tubes length inner diameter outer diameter ns material orientation number of ns thickness total surface area of the ns PCM material material hold-up specic latent heat storage capacity (latent heat)

Symbol

Value

Unit

− lap aap × bap Aap − VH

steel 2.7 320 × 360 0.1152 Marlotherm SH 170

− m mm m2 − kg

− nT lT din = 2rin dout = 2rout

P235GH 72 2.5 13.5 16.5

− m m mm mm

− − nF sF AF _total

AlMg2.5 perpendicular 312 0.25 59.8

− − − mm m2

− mP ∆hP ∆EP

R HD6070EA Rigidex 170 183 − 260 10 − 15

− kg kJ/kg kW h

8


Page 9 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Material properties Dierential Scanning Calorimetry (DSC) 135

Specic heat capacity cp (T ) and PCM

phase transition enthalpy ∆hP were determined by Heat Flow Dierential Scanning Calorime-

try (hf-DSC) using a NETZSCH DSC 404C hf-DSC. The PCM measurement procedure developed in the IEA SHC TASK 42/ECES Task 29 was used. 33 Aluminum crucibles with pierced lids were used and a Helium atmosphere was applied.

Laser Flash Analysis 140

In a laser ash analysis a laser pulse is applied to a coplanar

discoid sample and the time dependent temperature curve at the rear side of the specimen is detected. A NETSCH LFA 427 allowing for determination of the thermal diusivity is used. The thermal conductivity λ(T ) is determined from ρ(T ), cp (T ) and the thermal diusivity considering the error propagation. 34

Applied measurement procedures 145

The HTF thermophysical properties have been t-

ted to measured data taken from product information (Rev: 12/06), see Eqs. (S1), (S2), (S3), (S4). Tubes and ns material properties have been tted to measured data from hfDSC and Laser Flash Analysis, see Eqs. (S5), (S6). PCM material properties have been tted to measured data from hf-DSC and laser ash analysis, see Eqs. (S7), (S8).

Modeling 150

The shell and tube LHTES model considers convective heat transport in the HTF, transient forced convection heat transfer between HTF and the tube wall, heat conduction in axial direction of the tube wall, heat exchange with the PCM and conduction in radial direction in the PCM. The phase change phenomena are modeled using the apparent heat capacity method.

155

For all tubes of the LHTES unit, equal ow of the HTF and equal temperature distribution on both the shell and tube side is assumed. Thus, in the model only one tube is 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 55

considered and boundary eects near to the limit of the storage device, e.g. energy losses to the surrounding, are ignored. The ns are considered only indirectly by an increased heat conduction in the PCM. The thermal capacity of the ns is lumped into the thermal capacity 160

of the tube wall. A schematic view of the radial modeling domain is derived in Figure 2 B from the cross-section of the LHTES in Figure 2 A.

Figure 2: A: schematic of the LHTES; B: radial modeling domain for a single tube.

Heat transfer uid (HTF) The forced convective ow of the HTF in the tube is modeled in 1D in axial direction x. Heat transfer to the tube inner wall is considered. Radial uid ow, axial heat conduction in the uid, viscous dissipation, external forces and compressibility are neglected. 35 The internal energy equation in Eq. (1) is then expressed in terms of uid temperature TH with

du = cp (T )dT using specic heat capacity for incompressible ows cp ≈ cv .

ρH cp,H

∂TH ∂TH 2 = −vH ρH cp,H − q˙H ∂t ∂x rin

(1)

The initial and boundary conditions for Eq. (1) read:

TH (t = 0, x) = TH0 (x) ,

10

TH (t, x = 0) = THin (t)


(2)

Page 11 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Moreover, the uid velocity vH depends on the mass ow m ˙ H,total :

vH =

m ˙ H,total nT π(rin )2 ρH

(3)

In Eqs. (1) and (3), rin is the tube inner radius and q˙H represents the heat ux density, i.e. the heat transfer from uid to the tube. q˙H is calculated with the heat transfer coecient α to (4)

q˙H = α(TH − TW )

In Eq. (4), TH and TW are the uid and inner tube wall temperature at axial position x, respectively. The heat transfer coecient α is obtained from the following dimensionless numbers and correlations for forced convective heat transfer inside a tube. 36 Nu =

α din , λH

Re =

vH din , νH

Pr =

νH ρH cp,H λH

(5)

In Eq. (5), din = 2rin , λH is the thermal conductivity coecient and νH is the kinematic viscosity. All HTF properties are temperature dependent, i.e. ρH (T ), cp,H (T ), λH (T ), νH (T ), 165

with correlation given in Eqs. (S1),(S2),(S3),(S4). The calculation for local Nu as function of Re and Pr for laminar and turbulent ow as well as the transition region is given in Eqs. (S9),(S10),(S11).

Tubes and ns The wall temperature TW is modeled in 1D in axial direction x assuming a constant temperature in radial direction. Heat conduction in axial direction and heat transfer at the inner and outer tube wall are considered. The inuence of density changes on the internal energy is neglected. The internal energy equation, Eq. (6), is expressed in terms of the wall

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 55

temperature TW and the inner energy du = cp (T )dT assuming cp ≈ cv .

(Aρcp )W

∂TW ∂ = AT ∂t ∂x

∂TW λT ∂x

+ 2π (rin q˙H − rout q˙P )

(6)

The initial and boundary conditions read:

TW (t = 0, x) =

0 (x) , TW

λW

∂TW = 0, ∂x t,x=0

λW

∂TW =0 ∂x t,x=L

(7)

The thermal capacity of the ns is added to the thermal capacity of the tube wall. Accordingly, (Aρcp )W = (AT ρT cp,T + AF ρF cp,F ), with the subscripts T and F indicating the tube and n, respectively. AT and AF are the corresponding cross-sectional areas. In Eq. (6), λT is the thermal conductivity coecient of the tube, q˙H and q˙P represent heat ux densities from the HTF to the tube inner wall (see Eq. (4)) and from the wall to the PCM (see Eq. (8)), respectively. Finally, rin is the inner, rout is the outer tube radius.

∂TP q˙P = λP ∂r r=rout

(8)

All properties are temperature dependent, i.e. ρT (T ), cp,T (T ), ρF (T ), cp,F (T ), with correla170

tions given in Eqs. (S5), (S6).

Phase change material (PCM) The control volume of the PCM around the tube is simplied. A cylindrical shell domain is considered ranging from the outer tube radius r = rout to r = rend , see Figure 2 B. rend is p 2 . determined from the cross-section of the LHTES lled with PCM, rend = Aap /nT /π + rout 175

It is assumed that the transport phenomena in the cylindrical shell are controlled by heat conduction in radial direction only. Natural convection in the liquid phase, as well as diusion and radiation are ignored. Moreover, the ns reaching into the PCM are considered only indirectly by prescribing an improved thermal conductivity in radial direction of the PCM 12


Page 13 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


applying an eective thermal conductivity λP . By using the apparent heat capacity method, 180

an explicit treatment of the phase change moving boundary is avoided. A temperature interval for the phase change is assumed and a corresponding mushy transition zone is dened. Mixed material properties with a smooth transition between phases (solid, mushy and liquid) are then applied over the entire region of the PCM. The inuence of density changes on the internal energy is neglected. The internal energy equation for the PCM is expressed in terms of the PCM temperature TP applying the apparent specic heat capacity in Eq. (11) with du = c˜P (T )dT . For cylindrical coordinates and with explicit consideration of the temperature dependency of λP and c˜P we get: 37,38

1 ∂ ∂TP = ρP c˜P ∂t r ∂r

∂TP r λP ∂r

(9)

Eq. (9) is solved with the following initial and boundary conditions:

TP (t = 0, r) =

TP0 (r) ,

TP (t, r = rout ) = TW (t) ,

∂TP =0 λP ∂r t,r=rend

(10)

The correlations for the temperature dependent PCM properties, i.e. ρP (T ), λP (T ), c˜P (T ), 185

are given in Eqs. (S7), (S8) and (11). Experimental data of ρP (T ) and λP (T ) and corresponding predictions from these correlations are shown in Fig. 7. The correlation for the apparent specic heat capacity c˜P (TP ) reads:

c˜P = 1000.0 (cp,P + b1 φ) ,

with

cp,P = a0 + a1 TP

(11)

c˜P (TP ) is modeled using a combination of a linear term for the specic caloric heat cp,P (TP ) and a probability distribution function (PDF) φ(TP ) for the specic latent heat. Dierent PDFs have been considered which are: 190

1. Normal distribution (Eq. (A1)), where µ is the location parameter and σ is the scale parameter, 13



2. Gumbel Minimum (Eq. (A2)), also known as Extreme Value Type I distribution, where

µ is the location parameter and β is the scale parameter, 3. Adaptation of the Weibull distribution (Eq. (A3)), where γ >1.1 is the shape param195

eter, µ is the location parameter and α is the scale parameter, 4. Adaptation of the Lognormal distribution (Eq. (A4)), where σ is the shape parameter,

µ is the location parameter and m is the scale parameter. Together with the parameters a0 , a1 and b1 in Eq. (11) the PDF parameters have been used to t calculated c˜P (TP ) (prediction) to DSC measurement data. The best results in terms 200

of the residual value (quality of t) are obtained for the adapted Lognormal distribution (combination of Eqs. (11) and (A4). The results are shown in Fig. 3 where it can be seen that the experimental data is tted reasonably well. Noticable deviations exist only for the solidication curve with 0.5 K/min and around 120◦ C. It is concluded that a model based on the Lognormal distribution is the best option to represent individual DSC data. 120 measurement prediction

solidification A

100

0.5 K/min T_ = "hP = 260.4 kJ/kg

80

c~P [kJ/kg/K]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

60

! melting

0.5 K/min 237 kJ/kg

2 K/min 231.5 kJ/kg

2 K/min 224.3 kJ/kg 10 K/min 207.8 kJ/kg

40

10 K/min 182.9 kJ/kg

20

0 105

110

115

120

125 T [°C]

130

135

140

145

Figure 3: Experimental data from DSC for the apparent specic heat capacity c˜P (caloric and R HD6070EA PCM (black dotted line) and model predictions using latent heat) of Rigidex the relations in Eqs. (11) and (A4) (grey solid line). Note that six individually tted models are used to predict measurement lines for melting and solidication at dierent temperature rates. 14


Page 15 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

205


In Fig. 3, the values for the specic latent heat ∆hP are obtained from the integration of the peak area of the measurement data. It can be seen that higher heating and cooling rates

T˙ lead to peak broadening and increased hysteresis of the phase change. Both phenomena could be attributed to the eect of crystallization kinetics in the HD-PE. It is found that the relatively large variations in the values of ∆hP result from the initial condition of the 210

melting or solidication experiment, e.g. the initial crystallinity of the sample. 33

Stored energy The total stored energy E absorbed or released by the storage device between time [t0 , t] is obtained from the solution of the dierential equation:

dE =m ˙ H,total dt

TH (x=L)

Z

cp,H (T ) dT ,

E(t = 0) = E 0

(12)

TH (x=0)

with E 0 being the initial absorbed energy at time t = 0.

State of charge The state of charge, SOC is a parameter which indicates the extend to which a LHTES is charged relative to storable latent heat. SOC can have values ranging from 0 − 1 with 1 meaning molten and 0 meaning solid PCM. While for materials with a distinct phase change temperature SOC can be either 0 or 1, for materials discussed here SOC can reach any value in between and is temperature dependent. The local soc(TP ) is obtained from the integration of the latent heat φ(TP ) in Eq. (11)

Z

TP

soc(TP ) =

Z

+∞

φ(T )dT , −∞

φ(T )dT = 1 −∞

15


(13)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

In time domain we get:

∂soc(TP ) ∂soc(TP ) ∂TP ∂TP = = φ(TP ) ∂t ∂T ∂t ∂t | {zP }

(14)

φ(TP )

With TP and soc being functions of r, x the global SOC of the complete storage is computed as mean value by integration in radial and axial direction

dSOC = dt

R L R rend 0

∂soc(TP ) r dr dx ∂t rout R L R rend r dr dx rout 0

,

SOC(t = 0) = SOC 0

(15)

with SOC 0 being the initial state of charge at time t = 0.

215

Numeric solution Simulation The numeric solution is based on a xed grid technique. Partial dierential equations are transformed into a sparse system of ordinary dierential equations (ODEs) by applying a discretization in space. A nite dierence method, namely central dierence scheme, is applied to the energy equations for the HTF, the tube and ns and the PCM domain (for details see the Supporting Information). Central dierence and upwind schemes could prove very useful for the modelling of uid ows. 39 The equations for the HTF (Eq. (1)) and the tube and ns (Eq. (6)) are discretized in axial direction for i = 1, · · · , Ne mesh points. Moreover, at each discrete axial position i the equation for the PCM (Eq. (9)) is discretized in radial direction for j = 1, · · · , Nr mesh points. This yields a system of nonlinear ODEs with dimension Nx = Ne (Nr + 2). The state vector is given in Eq. (16).

x = [TH,1 , · · · , TH,Ne , TW,1 , · · · , TW,Ne , TP,1 , · · · , TP,(Ne ·Nr ) ]T

16


(16)

Page 17 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


For the computation of E and SOC in Eqs. (12) and (15), two additional dierential equations need to be solved. These two equations can be treated separately from the rest of the nonlinear ODE system, as they are pure quadrature equations. In Eq. (15), trapezoidal 220

numerical integration is applied based on the solutions of TP at discrete mesh points (solution of Eq. (9)). The ODEs can be eciently solved using algorithms for problems with an arbitrarily or banded sparse Jacobian. 40 For the performed simulations Ne = 50 and Nr = 7. A further increase of these values showed no inuence on the predicted temperatures. This yields

225

Nx = 450 and a sparse unstructured Jacobian with 867 elements. Applying Matlab's sparse reverse Cuthill-McKee ordering algorithm symrcm a banded Jacobian with bandwidth equal

9 is generated. The CVODE solver with banded Jacobian option from the SundialsTB Matlab Toolbox is used for the solution of the ODEs. 40 For the experiments shown in this contribution, simulation times did not exceed 0.7 sec (average 0.25 sec) using a Intel(R) 230

Core(TM) i7-4600U CPU @ 2.10GHz.

Parameter estimation Temperature measurements are considered for model validation, see section 'Experimental', with errors assumed to be independent-additive and normally distributed with a standard deviation σ = 0.1 K. Measurements are taken at discrete time instances (data sampling) and 235

for several experiments, each with a unique input signal (HTF ow and temperatures). The tting of model predictions to measurements and the estimation of parameter values is done by solving the unconstrained parameter estimation problem. The parameter output sensitivities can be generated using ODE solvers with sensitivity analysis capabilities. 4042 In this case study the sensitivities are generated using the CVODES solver. 40 The parameter

240

estimation problem is solved using Matlab's trust-region-reective algorithm lsqnonlin and user dened derivatives.

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 55

Fitting models and parameter estimation In the following the results of the parameter estimation are presented. The tting of model predictions to experimental data and the model validity are analyzed. The precision of the 245

parameter estimates is assessed by a statistic analysis. In the lab-scale experiments, measurements of the inlet and outlet HTF temperatures and internal PCM temperatures have been collected during transient operation with varying mass ows of the HTF. These measurements are: HTF inlet, HTF outlet, PCM 1, 2, 3, 4, see section 'Experimental'. After a careful analysis of the model and the measurements

250

the model parameters shown in Table 2 have been selected for parameter estimation. Four dierent models are considered (Mod1, Mod2, Mod3 and Mod4) each with a dierent set of parameters θ. They dier in the correlation used to model c˜P . Table 2: Vector of estimated parameters θ and respective reference values. Four dierent models (Mod1 , · · · , Mod4) are considered where the parameters θ5 , · · · , θ7 are taken from dierent equations. parameter

θ1

θ2

θ3

θ4

reference value Mod1 symbol origin symbol Mod2 origin symbol Mod3 origin Mod4 symbol origin

96.162 c2 Eq. c2 Eq. c2 Eq. c2 Eq.

0.1541 c3 (S8) c3 (S8) c3 (S8) c3 (S8)

1.7658 a0 Eq. a0 Eq. a0 Eq. a0 Eq.

211.61 b1 (11) b1 (11) b1 (11) b1 (11)

θ5

θ6

θ7

100.0 4.2661 1.00 µ σ Eq. (A1), Normal distribution µ β Eq. (A2), Gumbel distribution µ α γ Eq. (A3), adpt. Weibull distr. µ m σ Eq. (A4), adpt. Lognormal distr.

Each model allows for the adaptation c˜P and λP . Parameters not listed in Table 2 are kept constant using values given in Eqs. in (S1) - (S8), in Eq. (11) a1 = 6.5655 × 10−3 . 255

Moreover, a correction by a factor of 2.88 has been applied to the heat transfer coecient

α in Eq. (4). This is necessary since the used local heat transfer correlations (Eqs. (5), (S9), (S10) and (S11)) are very conservative and since the model does not account for the additional heat transfer at the front and rear tube plates of the LHTES. A second corrective 18


Page 19 of 55

coecient has been introduced to increase (Aρcp )W in Eq. (6) by multiplication with 2.13. 260

This is necessary to account for the additional thermal capacity of the tube plate, front and rear header and outer wall of the LHTES which are not considered. Both corrections have been determined by tting the model to empty LHTES measurements, i.e. before lling the LHTES with PCM.

Individual tting of experimental data 265

Two experiments, one for melting and one for solidication have been considered. In both, the HTF ow rate is kept constant with low ows and laminar ow conditions in the pipes with Reynolds Numbers of 2100 − 2220. In the melting experiment the temperature is increased from approximately 105 to 155◦ C and in the solidication experiment reduced from 155 to 105◦ C. The results are given in Table S1. 160 Mod1 Mod2 Mod3 Mod4

150

140

T [°C]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


PCM 1 PCM 2 PCM 3 PCM 4

130

HTF inlet HTF outlet

120

1.02 kg/sec

110

measurement prediction (Mod3)

100 0

10

20 c~P [kJ/kg/K]

30

0

20

40

60 time [min]

80

100

120

Figure 4: Individual tting of experimental data from melting at m ˙ H,total = 1.02 kg/sec. Left: Predicted apparent specic heat capacity c˜P . Right: Experimental data and model predictions for Mod3 (using adpt. Weibull distribution). For parameter values see Table S1.

270

The temperature proles for the melting experiments are shown in Fig. 4. The predicted internal PCM and external HTF temperatures t the measurements very well. The average error can be as low as 0.8 K when Mod3 (adpt. Weibull distr.) is applied. With the other

c˜P models slightly higher average temperature errors are obtained, see Table S1. Analysis of 19



the experimental data from melting experiment show that the temperature gradient in the 275

phase change region 120 − 140

◦

C is about 2 − 4 K/min. The comparison of the predicted

c˜P (Fig. 4, left) with DSC measurements with similar temperature gradients (Fig. 3) show very good agreement. 160 Mod1 Mod2 Mod3 Mod4

150

measurement prediction (Mod4)

PCM 1 PCM 2 PCM 3 PCM 4

140

HTF outlet T [°C]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

HTF inlet

130

120 1.02 kg/sec 110

100 0

20

40 60 80 c~P [kJ/kg/K]

100

0

20

40

60 time [min]

80

100

120

Figure 5: Individual tting of experimental data from solidication at m ˙ H,total = 1.02 kg/sec. Left: Predicted apparent specic heat capacity c˜P . Right: Experimental data and model predictions for Mod4 (using adpt. Lognormal distribution). For parameter values see Table S1. The temperature proles for the solidication experiments are shown in Fig. 5. In comparison to the other models, Mod4 (adpt. Lognormal distr.) gives much better tting of 280

internal PCM and external HTF temperatures with an average error of only 0.53 K, Table S1. In the solidication temperature region the temperature gradients are 4.5 − 20 K/min. The comparison of the predicted c˜P (Fig. 5, left) with DSC measurements with similar temperature gradients (Fig. 3) show a large discrepancy indicating structural errors in the model. In the parameter estimation procedure the optimizer compensates for these errors

285

leading to physically inapporiate predictions of c˜P . For the solidication (particularly at high temperature gradients), the most important structural model error certainly is the neglection of subcooling. This is indicated by deviations of the internal PCM temperatures between experiment and model. Another important structural model error is the neglection of the temperature hysteresis. This can be seen in 20


Page 21 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

290


the temperature shift between the predicted c˜P in Fig. 4 and Fig. 5. This hysteresis behavior has been conrmed in the DSC measurements (Fig. 3).

Simultaneous tting of experimental data All experiments were carried out with laminar ow conditions in the pipes with Reynolds Numbers of 2100 − 2220. 295

For the individual tting of data from melting or solidication experiments, the predictions are in very good agreement with internal PCM and external HTF temperature measurements. However, for the simultaneous tting of data from multiple (melting and solidication) experiments no equally good match can be found. In particular the internal PCM temperatures show considerable deviations caused by structural model errors. How-

300

ever, from an operational viewpoint the key variables are the HTF inlet/ outlet temperatures. Therefore, internal PCM temperatures are neglected for the simultaneous tting. The results are shown in Table 3, and Figs. 6, 7, and 8 (left). Table 3: Estimated parameter values for the simultaneous tting of data from eight experiments in Fig. 6. Dierent models Mod1 - Mod4 are considered (see Table 2). All parameter values are normalized with reference values given in Table 2.

θ1 Mod1 Mod2 Mod3 Mod4

2.219 2.762 1.851 0.518

θ2

θ3

θ4

θ5

θ6

θ7

3.715 3.514 4.070 6.539

0.832 0.476 0.772 0.673

0.956 1.236 1.045 1.095

1.237 1.265 1.328 1.328

1.619 2.201 2.724 1.917

1.132 0.947

21


error max. 3.87 3.47 3.94 3.74

[K] avg. 0.66 0.63 0.65 0.67

∆hP [kJ/kg] 201.2 207.5 188.5 190.0


T [°C]

160 measurement predicton (Mod2)

140 120

0.415 kg/sec

100

T [°C]

160 140

HTF inlet 120

HTF outlet

0.6 kg/sec

100

T [°C]

160 140 120

1.03 kg/sec

100 160 T [°C]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 55

140 120

1.98 kg/sec

100 0

50 time [min]

100

0

50 time [min]

100

Figure 6: Simultaneous tting of experimental data from melting (left) and solidication (right) at dierent HTF ow rates m ˙ H,total . Model predictions are shown for Mod2 (Gumbel Minimum distribution). Parameter values are given in Table 3).

22


Page 23 of 55

1.2

6P [W/m/K]

1 0.8

measurement prediction (Mod1-Mod4)

0.6 0.4 0.2 0 950

;P [kg/m3 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


900 850 800 750 80

100

120

140

160

T [°C]

Figure 7: Experimental data of PCM density ρP and thermal conductivity λP (uncertainties indicated as error bar) and corresponding tting (black lines) from the correlations in the Supporting Information - thermophysical data. Predictions for Mod1 - Mod4 using parameter values in Table 3. Generally for all c˜P models a very good agreement of external HTF (outlet) temperatures is obtained with a minor average error of 0.63 − 0.67 K. Comparing predicted c˜P in Fig. 8 305

(left), with data from DSC measurements in Fig. 3, very high deviations are found indicating structural errors of the model. Note that the predicted c˜P in Fig. 8 (left) stretches over a wide temperature range of 105 − 140◦ C. Despite these signicant deviations between model predictions and measurements for c˜P the temperature of the HTF at the storage outlet shows high agreement. Moreover, dierences in the prediction of c˜P for Mod1 - Mod4 (see Fig.

310

8, left) do not lead to signicant deviations in the error in predicted external HTF (outlet) temperatures, see Table 3. Fig. 7 shows experimental data and predicted ρP and the eective thermal conductivity coecient λP . The predicted λP is about three times higher than the measured values. This is the direct result of the conductivity enhancement, i.e. the ns installed around the tube.

23



315

Cross validation The quality of the model obtained from the simultaneous tting of experimental data has been assessed by a cross validation (Fig. 8). As expected, external HTF temperatures are predicted with very good accuracy, whereas internal PCM temperatures show relatively large deviations. The reasons are discussed before. 160 measurement prediction (Mod2)

140

T [°C]

120

Mod1 Mod2 Mod3 Mod4

100 80 0

10 20 c~P [kJ/kg/K]

30

0

50

100

150 time [min]

200

250

300

Figure 8: Cross validation using experimental data for m ˙ H,total = 1.09 kg/sec. Left: Predicted apparent specic heat capacity c˜P . Right: Experimental data and model predictions for Mod2 (using Gumbel Minimum distribution). For parameter values see Table 3.

Predicted stored energy and state of charge The stored energy E and state of charge SOC are calculated from Eqs. (12) and (15). Fig. 9 illustrates exemplarly the obtained curves for the cross validation experiment (Fig. 8) E0-20 1

E0-15

A

E0-10

!

0.5

E0-5

0

state of charge SOC [-]

320

stored energy E [kWh]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 55

E0 0

50

100

150 time [min]

200

250

300

Figure 9: Calculated values of stored energy E and state of charge SOC for the cross validation experiment (Fig. 8). E is multiplied by −1 and E 0 is the initial absorbed energy at time t = 0, see also Eq. (12).

24


Page 25 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


It is shown that E and SOC show a similar trend. The main dierences result from the fact that SOC indicates changes in the latent heat of the PCM only, whereas E indicates 325

changes in thermal energy of the whole storage, i.e. latent and sensible heat of the PCM and sensible heat of shell and tube materials. Accordingly, variations of SOC are limited to the the phase change temperature region. Variations of E do always occur with absolute values depending on the reference E 0 .

Summary and discussion 330

Thermal storages are part of highly integrated energy systems. It is therefore important to develop reduced models for ecient and robust numeric solutions. These system-level models are necessary for analysis of the design, of the interactions of the storage with other components and for the development of advanced control strategies. Publications with experimental validated LHTES models are rather scarce. Especially, studies using experimental data from

335

a realistic LHTES operation, including partial and full charging/discharging, variable HTF ow rates and inlet temperatures are not available. LHTES models for thermal networks or detailed ow sheet simulations need to accurately predict the following key variables: the HTF outlet temperatures and the state of charge which indicates the extend to which a LHTES is charged relative to storable latent heat.

340

The presented 2D shell and tube LHTES model has been successfully validated with experimental data from a realistic LHTES operation. The model includes detailed information about storage geometry and measured thermophysical data. It also includes a relation for the calculation of the state of charge. It easily can be adapted to several storage dimensions and dierent PCMs. For the studied shell and tube design a set of key model parameters for

345

this adjustment is determined based on parameter identiability analysis (see Appendix B for details). The model can be used for simulations on a system-level where, e.g. a complete industrial process is included.

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 55

Dierent distribution functions are used to directly consider the measured apparent heat capacity data of the PCM (using DSC) and it is shown how these functions can be imple350

mented. For the used PCM (HD-PE), the phase transition is noticeably aected by the heating and cooling rates. In the DSC measurements peak broadening and hysteresis eects can be observed, which is in line with results for many dierent PCM materials. 33 Storage internal PCM measurements show clearly the relevance of these eects which lead to subcooling and temperature hysteresis. For the studied lab-scale storage, it turns out that modeling

355

of these complex phase transition phenomena can be neglected, while still achieving satisfactory results in the prediction of LHTES outlet temperatures. However, from the specic results in this contribution only, no conclusive statement is permitted. Therefore additional studies on non-equilibrium modeling of phase transition of PCM, e.g. using macrokinetic models, 18 and their implementation in models for storage devices is highly desirable.

360

Acknowledgement This work was partly funded by the Austrian Klima- und Energiefonds within the programme e!MISSION in the project "StoreITup-IF" (number 848914) and by the Austrian Research Funding Association (FFG) within the programme Bridge in the project "modELTES" (number 851262).

365

Appendix A: Relations for the apparent specic heat capacity Four dierent PDF are taken from NIST/SEMATECH 43 . Note that, the selection is restricted to PDF with a closed (analytic) form. Symbol τ is used to indicate temperatures in

26


Page 27 of 55

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[◦ C]. Normal distribution: 1 (τ − µ)2 φ (τ ; µ, σ) = √ exp − 2σ 2 σ 2π

(A1)

N

Gumbel Minimum (Extreme Value Type I) distribution:

1 φ (τ ; µ, β) = exp β G

τ −µ β

exp − exp

τ −µ β

(A2)

Adapted Weibull distribution (adaptation, mirrored around µ):

φW (τ ; γ, µ, α) =

  

γ α

− τ −µ α

γ−1

exp

τ −µ γ α

, τ

Experimental Analysis and Numerical Modeling of a Shell and Tube

Recommend Documents