Data-based model reduction for phase change problems with convective heat transfer

Latent heat thermal energy storages (LHTES) exploit the high energy density of phase change material (PCM). The typically low thermal conductivity of PCM limits the charging and discharging rates and poses considerable challenges for dynamic storage operation. To operate LHTES efficiently and to exploit their full potential, new methods are required to obtain accurate and fast models for state of charge estimation and control tasks. In LHTES the heat transfer in low viscosity PCM is driven by conduction and also significantly by convective transport. In previous works, various high-precision models have been developed which employ finite element, difference and volume methods to solve the coupled Navier–Stokes and energy equations, but they incur large computational effort. In the present work, a novel, high-fidelity model reduction technique is proposed to achieve real-time capability while preserving high model accuracy. The idea is to short-cut the laborious solution of the Navier–Stokes equations by an efficiently parametrized, data-based model which approximates the stream function of the typical convection flow pattern by singular value decomposition. To account for the complexity of the solution-dependent flow domain, a suitable transformation method is proposed. The efficiency and accuracy of the proposed reduction method is demonstrated in typical operating modes.


Motivation
Due to the volatile nature of renewable energy sources, their efficient use requires the application of thermal energy storage systems (TES). TES are essential to decouple energy supply and consumption, which in turn improves the performance and reliability of energy systems as a whole [28]. In particular, the energy-intensive industries rely on balancing energy supply and demand to achieve the necessary flexibility for an effective use of resources and to incorporate renewable energy sources.
There are three main types of storage media in TES systems: sensible, latent and chemical, see Gil et al. [24]. In sensible TES, the energy is stored with increasing/decreasing temperature of the storage media, e. g. in the packed rock bed thermal storage. Latent heat thermal energy storage (LHTES) mainly use the phase transition to store energy in a small temperature range. The third storage mechanism is based on completely reversible chemical reactions to store and fully recover thermal energy, e.g. in oxidation-reduction reactions.
LHTES consist of phase change material (PCM) to store energy with high energy density and at an almost constant temperature level, see Agyenim et al. [1] and Zalba et al. [62]. A drawback of many PCMs is their low thermal conductivity which limits the charging and discharging rates and poses significant challenges for dynamic operation of the storage. Some methods to increase the thermal conductivity are summarized by Tao and He [56] and Ibrahim et al. [28]. Bondareva et al. [9] investigated the addition of nanoparticles in PCM to increase conduction. Typical LHTES configurations incorporate aluminum fins in PCM for thermal conductivity enhancement, such as the setup shown in Fig. 1. Due to the complex dependency of the temperature distribution on the total energy content and the low thermal conductivity of PCM, the thermodynamic state is distributed over the domain and thus the state of charge cannot be measured directly. For an efficient implementation of LHTES in industrial energy systems, knowledge of the distributed thermodynamic state and its dynamic behavior is of crucial importance. In order to realize state of charge estimation with an observer such as the extended Kalman filter as well as to enable modelbased control, an accurate and real-time capable model is required that contains all relevant effects. The two most relevant effects that occur when modeling the heat transfer of low-viscosity PCMs are conduction and convection. The computation of these effects with finite element, finite difference or finite volume methods for the coupled energy equation and the Navier-Stokes equations requires high computational effort and leads to simulation models which typically cannot be computed in or faster than real-time, see Kasper [29].

PCM modeling approach
Relevant heat transfer mechanisms in PCM are conduction and natural convection. Dutil et al. [21] and Liu et al. [35] review options for mathematical modeling and simulation of PCM. Analytical solutions only exist for a limited number of phase change problems, as for example for the one-dimensional Stefan-problem, see Radhakrishnan and Balakrishnan [44].
Fortunato et al. [23] state that the effect of natural convection is widely neglected in modeling PCM thermal storage systems due to its complexity. However, the legitimacy of such simplification strongly primary flow liquid ratio x switch liquid ratio which accounts for a switching flow regime depends on the type of PCM encapsulation and is potentially violated for PCM with low viscosity in the liquid phase. Bondareva and Sheremet [10] investigated the influence of geometric parameters of fins such as length and width on convection modes and thus on the melting rate. In a numerical study with an experimentally validated model, Vogel et al. [59] analysed the impact of natural convection on melting in so-called flat plate LHTES filled with the common high-temperature PCM sodium nitrate and potassium nitrate (KNO3-NaNO3) and found a significant heat transfer enhancement depending on the geometry of PCM enclosures due to convective transport. This finding was confirmed by Kasper [29], who numerically investigated the impact of natural convection on melting of PCM in an encapsulation depending on its orientation.
When considering heat transfer due to conduction in PCM in more than one dimension, numerical methods, as for example the enthalpy method and the effective heat capacity method, which were summarized by Liu et al. [35], have to be applied. Nedjar [40] state that especially finite element methods are able to handle complex coupled thermomechanical problems with various and complex boundary conditions. The effective heat capacity method was, for example, applied by Tenchev et al. [57] in a moving-mesh finite element model considering both conduction and natural convection. A similar approach was pursued by Kasper [29], using the effective heat capacity method and an adaptation of the finite difference code published by Seibold [51] to solve the twodimensional Navier-Stokes equations arising in convection modeling. This provided an experimentally validated coupled finite element/ finite difference model for PCM cavity simulations and is adopted in the present work as a reference model for the model reduction approach.

Innovation
The modeling approaches for PCMs described in the above section are computationally highly demanding and not suitable for real-time models. Fig. 2 outlines the main contribution of this work -a model reduction method aiming to cut down computational requirements strongly while retaining high accuracy. The energy and the Navier--Stokes equations are coupled via the velocity field u. A first analysis of the method proposed by Kasper [29] showed that the Navier-Stokes equations require about 80% of the computation time. Therefore, in this work a novel reduction method is developed to model the relevant convection effects in PCM in a simplified form without fully solving the Navier-Stokes equations. The model reduction approach is data-based and relies on the stream function derived from the velocity field obtained in a reference simulation. A suitable transformation is introduced to map the active flow domain (whose shape depends on the solution itself) to a unit domain to account for the changing flow domain of the PCM. The stream function snapshots are decomposed into modes of space and time using singular value decomposition (SVD). It is found that the temporal behavior of the dominant modes of the velocity field can be obtained from thermal properties of the domain, such as liquid ratio and temperature gradient. The resulting reduced model is able to reconstruct the stream function and thus the velocity field without solving the Navier-Stokes equations in the real-time simulation. The new model reduction method is evaluated based on simulations of typical storage operating modes and shows highly accurate results achieved with considerably lower computation times.

Overview
Model reduction methods have been a major research topic in recent decades and have been addressed in numerous reviews. Benner et al. [5] conducted a survey of model reduction methods for linear systems and Reis and Stykel [45] provided an overview of approaches in coupled systems. Antoulas et al. [2] divided the reduction methods into two main categories: moment matching based methods and SVD-based methods. The former group can be implemented iteratively but does not have global error bounds, see Antoulas et al. [2] for an overview of linear implementations and Astolfi [4] for the nonlinear enhancement of this method. The approach developed in this paper is based on the latter group. An example of SVD-based reduction methods in linear systems is balanced truncation, where the system is transformed via principal component analysis to a basis in which the hard-to-reach states are simultaneously difficult to observe and are simply truncated for the reduced model. This method was first introduced by Moore [38], and a short overview of implementations is given by Gugercin and Antoulas   Fig. 2. Coupling of the dominant effects in the detailed and reduced model. [26]. Other SVD-based methods for linear systems are the Hankel-norm approximation and the singular pertubation, see Glover [25] and Liu and Anderson [36], respectively. Proper orthogonal decomposition (POD), first introduced by Lumley [37], is a popular method for nonlinear complexity reduction problems. The basic idea is to extract the most significant characteristics of a system's behavior and represent them in a set of orthonormal basis vectors. For this purpose, the POD basis vectors are determined empirically by examining sample data collected over a range of relevant system dynamics and identifying the most energetic modes. The system's governing equations are projected onto the reduced-order subspace defined by the POD basis vectors. This low-dimensional description of the system's dynamics is obtained directly from the Galerkin projection of the governing equations on the POD modes. This yields an explicit POD reduced model that can be solved instead of the original system, see Chinesta et al. [14] and Volkwein [60].
More recent methods for model reduction of nonlinear systems are the dynamic mode decomposition (DMD) developed by Schmid [50] and the proper generalized method (PGD), see Chinesta et al. [15] and Néron and Ladeveze [41]. In contrast to the POD method, the PGD technique is an "a priori" iterative model reduction method which does not rely on particular bases. DMD is data-driven and, in contrast to POD, the DMD modes attempt to describe the dynamics in the timeseries rather than best reconstructing a dataset, see Tu et al. [58]. Brunton et al. [13] proposed the data-based regression method SINDY, combining DMD and sparse regression and exemplified it in applications of uniform flow vortex shedding and discovering the sparse equations structure in a noisy Lorentz system. Other model reduction techniques are based on the combination of machine learning and SVD, which demands a lot of training data, see e.g. Prasad and Bequette [43]. Furthermore, Swischuk et al. [55] demonstrate the importance of embedding physical constraints as well as of incorporating knowledge such as conservation laws within the learned models.

Model reduction in fluid dynamics
Particularly in the field of fluid dynamics, modeling often reaches the limits of computing power and is therefore the subject of current research regarding model reduction methods. The coupled partial differential equations to be solved in this area represent a particular challenge. Rowley and Dawson [48] and Lassila et al. [32] summarized model reduction approaches in flow analysis. Bistrian and Navon [7] compared the DMD and POD approaches to derive a reduced model for the shallow water equation. Rowley et al. [47] presented a framework for applying POD to compressible fluids with small temperature gradients and moderate Mach numbers. Rowley [46] successfully applied a balanced POD method to a channel flow. Dumon et al. [20] used the PGD method to solve the Navier-Stokes equations in the case of a liddriven cavity for different Reynolds numbers. Liberge and Hamdouni [34] applied the POD technique to a flow around an oscillating cylinder using a fictitious domain. Stabile and Rozza [54] compared two different pressure stabilisation strategies for POD methods in a lid-driven cavity problem and in a flow around a circular cylinder for moderate Reynolds number. Kutz [30] emphasized the potential of deep learning in fluid dynamics, but also pointed out unresolved issues such as the number of layers/nodes and the amount of training data required. Applicatons of reduction techniques using neural nets in fluid dynamics can be found in Wang et al. [61] for forced convection and in San and Maulik [49] for natural convection.
Most of the model reduction approaches for fluid dynamics presented are highly complex and laborious to adapt. One research gap is the lack of methods to efficiently reduce dominant flow patterns to a simple model. In the present work, an easy-to-implement data-based approach is developed to address this issue.

Model reduction with varying domains
A large number of problems is found on partial differential equations with time-varying spatial domains, such as melting/solidification processes, crystal growth, and hydraulic fracturing. However, only a few model reduction approaches have been developed specifically considering this aspect so far. Armaou and Christofides [3] expressed the partial differential equation system of a diffusion-reaction process with respect to an appropriate time-invariant coordinate derived from an analytical expression. A POD method was applied to the resulting transformed time-invariant system. Fogleman et al. [22] applied the POD to flows of an engine combustion process within a time-varying domain. They transformed the velocity fields to a fixed domain in such a way that the divergence-free property remained preserved, which is relevant also in the present context. Dauvergne and del Barrio [18] presented a simulation-free POD model reduction method for multidimensional heat conduction problems with phase change. They decomposed the problem into a heat conduction and source-term problem for the phase change. Narasingam and Kwon [39] partitioned the domain into multiple subdomains and then applied POD to each local domain. Recently Sidhu et al. [52] developed a model reduction method for the fracture propagation in a hydraulic fracturing process with moving boundaries. They could consider the time-varying spatial domain as a time-invariant one by assuming the fraction width zero where the fraction has not propagated into the domain.
Our newly developed model reduction method solves the problem of varying flow domains by first converting the velocity field into the stream function which preserves the continuity of the flow per definition. Then the stream function is mapped to a unit domain in which the dominant modes can be efficiently identified.

Main contributions
To the best of the authors' knowledge, no model reduction technique for convection problems with phase change and the associated problem of a varying flow domain has been presented so far. The main contributions of this paper are as follows: • A novel and efficient high-fidelity model reduction concept for convective heat transfer in phase change problems is developed, which takes changing flow domains into account. • The new approach considerably reduces the computational effort and realizes the real-time capability of the model. • Simulation studies are conducted and discussed, showing the effectiveness and accuracy of the proposed model reduction in typical operation modes of an LHTES application. • In contrast to the solution of the complete Navier-Stokes equations, the accuracy of the reduced model is less dependent on the time step and mesh size.

Paper structure
The paper is organised as follows: In Section 2, the detailed model and its fundamental equations are presented. Section 3 describes the model reduction approach through stream functions, transformation and SVD as well as the assumptions made. In Section 4, simulation studies to demonstrate the effectiveness of the proposed model reduction are performed. In Section 5 the results of the model reduction are compared and different aspects are discussed.

Detailed model
The reduction approach is demonstrated using a two-dimensional rectangular PCM domain with aluminum fins as heat conducting structures. For reasons of symmetry, modelling of a LHTES configuration as illustrated in Fig. 3 is reduced to a single, so called PCM cell, consisting only of one fin section of the storage.

Fundamental Equations
A detailed multi-phase thermodynamic model including melting/ solidification, heat conduction, and natural convection in two dimensions is considered for the PCM cell, developed and thoroughly documented in Kasper [29]. It underlies the basic conservation laws of mass, momentum and energy. Multiple PCM cell models were successfully coupled with a Ruths steam storage in a co-simulation in Pernsteiner et al. [42]. In order to simplify modeling of a PCM cell to the significant occurring phenomena, the following assumptions are made therein: • Heat transfer is driven by conduction and natural convection, both being relevant phenomena. • Material properties, apart from density, are constant in each phase but they can take on different values for the liquid and solid phases. • Density is assumed constant and the same for both phases except for the buoyancy term in the liquid domain (Boussinesq approximation). Therein, the only body force arises due to the gravitational acceleration in vertical direction. This is a common and widely used approach, see Vogel et al. [59]. • Phase change is modeled by means of the apparent heat capacity method (see [8,17]), meaning that the phase transition takes place in a small temperature region T ∈ [T m − ε, T m + ε] defined by the melting temperature T m and mushy region parameter ε.
• The phase ratio only depends on T (no temperature hysteresis, no dynamics, no rate dependency, no subcooling). • No relevant surface effects nor dentrite growth occur.
• Dissipation and pressure terms are neglected in the energy Eq. (1) due to small Eckert numbers, Ec≪1. • The depth of the PCM enclosure in z-direction, is assumed large enough for wall boundary layer effects to be negligible, hence the problem is reduced to two dimensions (no heat flow or convection in z-direction).
The governing equations within the framework of these assumptions are the energy Eq. (1), continuity Eq. (2) and Navier-Stokes Eqs. (3) as follows: Therein, the temperature field T = T(x, y, t) is treated as dependent variable in the energy Eq. (1), and the velocity field u = [u, v] T with spatial components u(x, y, t) and v(x, y, t) is treated as dependent variable in the Navier-Stokes Eq. (3), see Chorin [16]. The variable p PCM stands for the pressure in the PCM. The symbols ρ PCM , c, k and μ denote the parameters density, apparent heat capacity heat conductivity and dynamic viscosity, respectively. The force density f describes the buoyancy force which is calculated via the Boussinesq approximation [11] as given in Huang et al. [27] by the volumetric thermal expansion coefficient β, the constant (reference) density ρ 0 , a reference temperature T ref , and the gravitational standard acceleration vector g, g = g with g = 9.81 m/s 2 .

Discretization
The PCM cell conduction/convection model developed by Kasper [29] uses a finite element discretization to model the time-dependent energy Eq. (1) and applies a finite difference scheme to solve the incrompressible Navier-Stokes Eqs. (2) and (3) in each time step. The obtained velocity field is used for the next time step of the energy equation.
The finite element method uses a standard Galerkin finite element approach with four-noded bilinear rectangular elements for the energy Eq. (1). The set of incompressible two-dimensional Navier-Stokes Eqs. (2) and (3) is treated separately from the heat transfer part of the model via a finite difference discretization based on a MATLAB® code available on the course homepage(http://math.mit.edu/~gs/cse/) of "Computational Science and Engineering" at the Massachusetts Institute of Technology. This open source code, created and documented by Seibold [51], was implemented and extended to meet specific requirements in the diploma thesis of Kasper [29]. The main numerical concept used therein is the fractional step method [16,19], which is applied to split the Navier-Stokes system into equations that are significantly simpler to work with.

Boundary conditions
The geometry of the considered PCM cell suggests symmetry reduction of the simulation to single aluminum fin sections. Consequently, the symmetry boundaries, in this case chosen as upper and lower boundaries of the computational domain, are treated as adiabatic: For the left and right boundaries type-3 boundary conditions, also known as Robin or heat-transfer boundary conditions, are prescribed following Newton's law of cooling, where α in , α out are heat transfer coefficients, q is the specific heat flux across the boundary, and T in , T out present boundary temperatures at the left and right wall surfaces, respectively. Regarding the velocity field u = [u, v] T , no-slip boundary conditions are set for the domain boundaries, Furthermore, the velocity is set to zero if the PCM is not completely liquefied: This is also done for any solid material (here, aluminum) in the computational domain, see Kasper [29].

Features of the model
The model presented above is able to model heat transfer by convection and conduction in PCM. It fully solves the Navier-Stokes equations, thus providing a highly accurate velocity field in the liquid domain. It was validated in Kasper [29] against experimental data from Brent et al. [12], resulting in an error of 3% in the liquid fraction. However, for realistic problem sizes it is not real-time capable by a factor of 5 or more and it is found that the computation of the Navier-Stokes equations consumes about 80% of the time. To achieve real-time capability, a model for high accuracy heat transfer and good accuracy of the velocity field is developed, utilized to determine the convection terms in the energy equation and hence achieve the overall solution quickly. This enables a model suitable for estimating the state of charge and controlling an LHTES.

Reduced model
The basic idea is to replace the solution of the Navier-Stokes equations by a simplified data-driven model, that explains the velocity field from few selected properties of the energy equation, such as the size of the relevant domain in which convection acts and the temperature distribution driving the convection. The reduction relies on snapshots of data of the velocity field derived from high-fidelity simulation studies. First, the assumptions under which the proposed model reduction is admissible are discussed. Then the pre-processing of the sampled data is described: flow domain identification, representation of the velocity field via the stream function and the transformation to a unit domain. The data is decomposed into dominant modes of space and time through SVD. A model for the magnitude of the dominant modes is developed and the reduced model is presented, see Fig. 4.

Assumptions
The proposed model reduction approach is based on the following assumptions: • The flow phenomena consist of one or only a few dominant flow patterns. • The fluid is incompressible and its motion can be fully described in two dimensions, so a scalar stream function exists which defines the velocity field. • The flow domain varies slowly compared to the time step size. This assumption is justified due to the low thermal conductivity and the high latent heat of PCMs. • Convection can be neglected during solidification, see Pernsteiner et al. [42]. During discharging, the PCM solidifies on the heatconducting structure where it acts as an insulator due to its low thermal conductivity. Therefore the temperature gradient in the remaining liquid domain is low and convection becomes negligible.

Flow domain identification
Representative simulations in Kasper [29] and Pernsteiner et al. [42] show that relevant convection phenomena only occur in a specific region within the entire domain which will be called "primary flow domain" in the following. This primary flow domain flow = flow (t) is identified in the overall domain according to the following properties: The primary flow domain flow (t)⊆ is thus defined as follows: The primary flow domain is the definition area of the stream function and the starting point of the transformation into a unit domain. The proportion of the liquid PCM in the primary flow domain to the entire domain x liquid , with x liquid ∈ [0, 1], is introduced as primary flow liquid ratio:

Stream function
The stream function is a well-known scalar representation to describe a planar velocity field of an incompressible fluid and was first introduced by Lagrange [31]. Partial derivatives of the stream function yield the velocity field, which in turn is automatically divergence-free and thus satisfies the continuity Eq. (2), see Bestehorn [6]. Therefore, the stream function is chosen as a representation of the velocity field and is derived from snapshots of simulation data in this work. The stream function Ψ = Ψ (x, y, t) : flow →R is defined through Therein, u and v are the spatial components of the velocity field u in the spatial x and y-directions.

Transformation
To account for the varying liquid area of the PCM, the stream function in the time-dependent primary flow domain Ψ = Ψ (x, y, t) : flow (t) →R is mapped to a stream function in a time-invariant unit domain Ψ * = Ψ * (ξ, η, t) : unit →R via the coordinate transformation Ω : flow (t)→ unit , see Fig. 5. The transformation itself depends on the shape of flow and hence ultimately on the temperature field of the solution and is evaluated using equally spaced nodes in x-direction and linear interpolation. This is done with the aim to simplify the flow patterns and ease their treatment as Ψ * is now defined on a domain of constant shape.

Singular value decomposition
Sirovich [53] introduced the method of snapshots to efficiently determine POD modes. Therefore, the data set is selected as time snapshots containing the spatial distribution and reflecting the system dynamics. Indicating by Ψ * i the vector of values of Ψ * at the time t i , the data matrix Z, consists of n snapshots in time (columns) from the spatial distribution of the stream function values (rows). Dominant modes in time and space are identified from the snapshots. The economy-size SVD is utilized to decompose the snapshot matrix, by projecting it onto an orthonormal basis given by the left singular vectors U. The matrix U contains column-wise the spatial modes, and V represents the temporal behavior of each mode row-wise. The diagonal matrix Σ comprises the singular values which indicate the signal energy content of each mode. Identified by the largest singular values, typically only a few modes, represent the major contributions. Their dominant modes can be used as a reduced description of the relevant behavior of the system. The reduced model (U red , Σ red and V red ) is obtained by selecting the corresponding rows/columns from U, ΣandV, respectively.

Temporal behavior of the spatial modes
To establish the temporal behavior of the spatial modes, it is common to apply the Galerkin method to project the partial differential equation system onto the spatial basic functions, yielding a system of ordinary differential equations for the modal coordinates. As these projections (numeric integration over the domain) are cumbersome, a different, data-based approach is chosen in this work. The temporal behavior of the dominant modes given by the matrix V red as well as the magnitude of the velocity field itself seems to depend on quantities of the energy equation, such as the temperature distribution in the PCM or the size of the primary flow domain. Therefore, it is attempted to explain the temporal behavior of the spatial modes V red by the primary flow liquid ratio x liquid , the maximum temperature in the flow domain T max , the maximum temperature difference in the domain T spread , with T min as minimum temperature in the domain, and the temperature near the heated wall (left side) in the horizontal symmetry plane of the PCM cell T corner , see Fig. 3. These four measurements, form the basis of a cubic polynomial approach for a regressor matrix R.
In order to consider a switching flow regime in the regressor R at a certain liquid ratio x switch , an activation function γ, is introduced. The regressor matrix R consists of 62 regressors (columns) at the n snapshots in time (rows) and is defined as follows: r 2 2 r 4 r 2 3 r 1 r 2 3 r 2 r 2 3 r 4 r 2 4 r 1 r 2 4 r 2 r 2 4 r 3 ].
A matrix Θ containing the parameter vectors is identified for an ARX model using a least square algorithm: The estimation of the current magnitude of the spatial modes is then given by with r as regressor using the current measurement values, r = [ r con r lin r qu r cub γr con γr lin γr qu γr cub ].

Reduced model architecture
The reduced model described above consists of an offline preprocessing and an online real-time capable part.
During pre-processing, dominant modes are identified by the SVD of stream function data to obtain a reduced system, U red , Σ red and V red . An ARX model approximates the modal coordinates v red of the dominant spatial modes U red based on four characteristic measurements (20). Fig. 6 shows the online architecture of the reduced model developed in the section above. It consists of the finite element code of the detailed model for the energy equation. The finite difference approach for the Navier-Stokes equations is replaced by the reduced stream function model, After back-transforming the stream function Ψ * within the unit domain to the stream function Ψ within the primary flow domain, the velocity field is obtained by numeric differentiation and inserted into the energy Eq. (1).

Simulation studies
As described in the previous section, our reduction method consists of two parts: the pre-processing/model creation and the real-time model. To demonstrate the novel method, a suitable geometry is first defined and the material parameters are determined. Then simulations are performed to generate training data for the creation of the reduced model. Two additional load profiles are defined to validate the reduced model against the detailed high-fidelity model. Finally it is shown that the reduced model can be calculated on a coarser grid and larger time step without significantly compromising the result.

Simulation setup
The simulation setup is based on Pernsteiner et al. [42]. The size and shape of a PCM cell has been chosen to allow load profiles on an industrially relevant scale.
In order to demonstrate the capability of the model reduction approach, five different model instances are simulated. The high-fidelity model, which completely solves the Navier-Stokes and energy Eqs. (1)-(3), relies on the code developed by Kasper [29] and is described in Section 2. In addition, it is used to generate training and validation data for the model reduction method. The regular-reduced model fully solves the energy equation, but uses a reduced stream function model instead of the Navier-Stokes equations. The coarse-reduced model is computed on a coarser mesh and larger time step to demonstrate mesh and time step independency of the model reduction approach and is compared to a coarse-high-fidelity model. A conduction-only model, which completely neglects natural convection, acts as a reference simulation.

Geometry
The geometry of the PCM cell is defined in Fig. 3 and Table 1. An insulation layer is attached to the right side of the PCM cell to reduce heat loss to the environment. The left side of the PCM cell has good heat transfer properties for charging and discharging via heat flows. The assumed heat transfer coefficients are listed in Table 2.

Material parameters
The PCM is an eutectic mixture of potassium nitrate and sodium nitrate KNO 3 -NaNO 3 (see [59]), enclosed in an aluminum encapsulation. The PCM is either liquid (L) or solid (S). The material properties of the aluminum encapsulation and the PCM are listed in Table 3. The density of solid PCM is used for both phases. A mushy region parameter of ε = ±0.5K was chosen. Influence of the variation of this parameter on the result of the high-fidelity simulation was studied in Kasper [29].

Load profiles
The PCM cell is charged/discharged by a heat flow resulting from the temperature difference on its left side (7). The input temperature on the left side, T in , is set as shown in Fig. 7 to generate training and validation data. Five different load profiles are used to train the reduced model: Four load profiles with a constant temperature for 4-h operation time and one load profile with a varying cycle for 5-h operation time. The performance of the reduced model is demonstrated using two 5-h validation load cycles. The right side of the PCM cell is insulated and the temperature responsible for the heat loss, T out , is fixed at 20 • C for all load profiles (8). The initial temperature in the PCM cell is T 0 = 218 • C.

Mesh and time step size
A mesh of square elements with a side length of Δx = Δy = 0.5 mm discretizes the geometry of the PCM cell of the high-fidelity model, the conduction-only model and the regular-reduced model. A time step of Δt = 0.1s is applied for the temporal discretization. These settings resulted as optimal values from an independence study of mesh size and time step. In order to demonstrate the capability of calculating the reduced model on a coarser mesh and larger time step, the side length of the elements and the time step is doubled for the coarse-reduced model. The properties of the different models are summarized in Table 4.

Settings of the reduced model
The reduced model was created as described in Section 3. The data from the high-fidelity simulations are used and only the charging intervals for the model training are considered (since natural convection during discharging is neglected). The SVD of the stream function data results in a single dominant spatial mode, see Fig. 8a. The temporal behavior of this dominant mode is characterized by the ARX model (28)-(21) using properties of the domain, e.g. temperature spread or liquid ratio in the PCM cell. The parameter for a switching flow regime in (21) is set to x switch = 0.5 and the excellent fit of the ARX model is shown in Fig. 8b.

Results
The results are evaluated under two aspects: accuracy and computational effort. These criteria are compared between the five different models listed in Table 4.

Solution accuracy
The stored enthalpy in the PCM cell is essential in order to assess the state of charge and to control the LHTES. Therefore, Fig. 11 compares the stored normalized enthalpy, of the different models for the training and validation load cases, respectively. In (31) H 0 is the initial enthalpy and H latent the latent heat of the PCM cell. The overall error is listed in Table 5. As a second criterion, the shape and progression of the melting fronts in the PCM are evaluated qualitatively, see Fig. 12.

Computational effort
The computational effort is assessed through the real-time factor, which is listed for the different models in Table 5. Eq. (32) consists of the computation time t comp required to simulate a PCM cell until the end time t end . The computation was performed using MATLAB® on the processor Intel® Core™; i7-8550U with a base frequency of 1.8 GHz.  Table 3 Material properties of the LHTES.

Discussion
At the end, the stored enthalpy of the regular-reduced model and the coarse-reduced model is 7.9% and 7.1% lower for the validation load cycle A as well as 3.3% and 5.8% lower for the validation load cycle B than of the high-fidelity model, respectively. The coarse-high-fidelity model differs by 6.3% and 3.7% and the conduction-only model shows an error of 35.4% and 26.5% for the validation load cycles A and B, respectively, see Fig. 11 and Table 5. The coarse-reduced model is as good as the regular-reduced model, but more than 8 times faster and even 44 times faster than the high-fidelity model. Therefore, the coarsereduced model requires only 0.65 h instead of 28.6 h computation time for 5 h simulation time. Fig. 10 shows a comparison between accuracy and computation time. A coarsening of the grid and an enlargement of the time step causes a significant decrease in accuracy in the highfidelity model, but not in the reduced model. The shape and progression of the melting fronts is slightly different between the reduced and the high-fidelity models, but still in good agreement after multiple charging/discharging cycles, see Fig. 12 for the validation load cycle A. Therefore, the reduced models are able to accurately estimate the state of charge, determine the location of the melting fronts, as well as serve as a basis to control the LHTES in real-time.

Mesh coarsening and time step enlargement
Solving the Navier-Stokes equations is often sensitive with respect to the chosen spatial discretization Δx = Δy of the high-fidelity model. The spatial discretization, in turn, determines the discrete time step Δt of the  simulation. The relationship between Δt, Δx and the velocity u is given by the Courant-Friedrich-Lewy number, which should be less than 1 to correctly resolve convective heat transfer, see Lewy [33]. The CFL number indicates the maximum number of cells which a considered quantity passes through per time step. In the novel model reduction approach, the computation of the Navier-Stokes equations is replaced by a stream function model. Mesh and time step sizes are therefore only limited by the solution of the energy equation through the finite element model. As a result it is seen that the mesh can be significantly coarsened and the time step can be enlarged. This reduces the computational effort enormously while maintaining high accuracy.

Further improvements
As seen above, the results of the proposed model reduction method are highly satisfactory. In order to further improve the model reduction approach in future work, two aspects are highlighted below: the transformation of the stream function into a unit domain and the decomposition of the stream function using SVD.

Transformation
The transformation of the stream function into a unit domain utilizes a simple and robust linear interpolation. This distributes the arising distortions evenly, but also leads to a distorted mapping of the flow's boundary layers. Additionally, the authors studied piecewise-linear mapping approaches that are designed to keep the boundary layers undistorted, but the similarity of the transformed stream function could not be further improved yet. Still, more complex transformations that seek to accurately represent the boundary layer structure while compressing the wake interior could yield a simpler stream function representation via the SVD method and should thus be studied deeper in the future.

Radial basis functions
The dominant modes are extracted from the stream function via SVD. The dominant modes approximate the grid values of the stream function optimally. For the proposed model reduction approach, however, the velocities, i.e., the spatial partial derivatives of the stream function are the actual quantities of interest. Selecting more than one mode for the model reduction approach improves the grid values but not necessarily their derivatives. In an improvement of the presented method, the SVD weighting could be adjusted to focus on the precision of representing the partial derivatives of the stream function and radial basis functions could be utilized to represent the stream function modes and allow to evaluate accurate derivatives of the stream function.

Conclusions
The model reduction method for dominant flow patterns developed in this work replaces the Navier-Stokes equations with a reduced stream function model. The stream function model is data-based and parametrized from simulations of a high-fidelity model. In order to consider the solution-dependent flow region (melting and solidification processes), the stream function is mapped from the original flow domain to a unit domain. While the shape and size of the stream function differ greatly in the original domain at different times, the transformed stream functions show similar shape. The reduction of the stream function model is SVDbased and the velocity field of the flow domain can be reconstructed solely from properties of the flow domain, e.g. the temperature distribution. The reduced model can be computed on a coarser grid and therefore on a larger time step without significantlydecreasing its accuracy.
The novel contribution of this work is the easily adaptable reduction technique for problems with varying domains. Its accuracy and computational speed is demonstrated in simulation studies. The reduced model can be computed up to 44 times faster than the high-fidelity model. In the validation cases, the reduced models yielded a maximum error of 7.9% instead of 35.4% error caused by a conductiononly model that completely neglects natural convection. Thus the reduced model is sufficient for the intended applications of state of charge estimation and model-based control of an LHTES. The model captures the dominant dynamics and the remaining uncertainties, which are in the range of the actual accuracy of the numerical model of 3%, can be compensated by an observer or a controller, respectively. To further improve the reduction efficiency, radial basis functions and an advanced transformation of the stream function are emphasized.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.