Model reduction by the Modal Identification Method in forced convection: Application to a heated flow over a backward-facing step

International Journal of Thermal Sciences 49 (2010) 1354e1368

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

Model reduction by the Modal Identification Method in forced convection: Application to a heated flow over a backward-facing step Y. Rouizi a, M. Girault a, *, Y. Favennec b, D. Petit a a

Institut Pprime, CNRS, ENSMA, Université de Poitiers, Département Fluides, Thermique, Combustion, ENSMA - Téléport 2, 1 avenue Clément Ader, BP 40109, 86961 Futuroscope Chasseneuil Cedex, France b LTN, UMR CNRS 6607, Université de Nantes. La Chantrerie, Rue Christian Pauc, BP 50609, 44306 Nantes Cedex 03, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 29 June 2009 Received in revised form 22 February 2010 Accepted 23 February 2010 Available online 20 April 2010

This numerical study focuses on the use of the Modal Identification Method to build reduced models for problems of heat convection and diffusion. The principle is to minimize a cost function based on the difference between the outputs (velocity and/or temperature) of a detailed model and the outputs of a reduced one. The reduced model structure is defined from the partial differential equations governing fluid mechanics and heat transfer in the physical system. In this paper, an advectionediffusion problem is studied: forced heat convection is considered with an incompressible, stationary, laminar 2D flow. Physical properties of the fluid are temperature independent, hence velocity is independent of temperature. The system under consideration is a channel flow over a backward-facing step with a timevarying heat flux density applied upstream of the step. Three types of reduced models have been investigated: steady fluid mechanics only, unsteady heat transfer for a given constant Reynolds number, and unsteady heat transfer for any constant Reynolds number within the range [100:800]. In this last case, the thermal reduced model is weakly coupled to the fluid reduced one. Results show that reduced models fit very well with detailed ones, and allow a large decrease of computing time.
2010 Elsevier Masson SAS. All rights reserved.

Keywords: Model reduction Modal identification method Optimization Particle swarm optimization Forced convection

1. Introduction In the last decades, the detailed computation of complex physical systems such as flows in channels, combined with heat transfer, has received a growing attention. The applications may concern the direct conception of cavities or pipes along with the optimization of such systems and even the application of on-line control of e.g. flows or temperature. The numerical simulation of such flows, and consequently the optimization and further the on-line control, require high computational resource even though the power of the computers and the improvements of detailed simulation algorithms have been growing exponentially. To face these computational drawbacks, the use of reduced order models presents undeniable interest to the engineer, especially when the solving of large matrix systems requires significant computing times and memory storage. Reduced models are composed of a restricted number of equations obtained through different mathematical approaches, such as solving eigenvalue problems and optimization problems. Many reduction methods have been developed, for linear or non-linear problems. Although we will cite a few ones in the * Corresponding author. Tel.: þ33 5 49 49 81 34. E-mail address: [email protected] (M. Girault). 1290-0729/$ e see front matter
2010 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.ijthermalsci.2010.02.011

following, the present paper is dedicated to the latest extensions of the Modal Identification Method. The reader interested in reviews on reduction techniques can refer to [1] or [2,3] for instance. When considering reduction techniques applied to linear systems, one may cite methods using a change of basis, like the modal methods and the balanced truncation method [4]. In the frame of heat transfer by advection and diffusion, modal methods consist in solving an eigenvalue problem relative to the advection-diffusion operator, and to make a selection of dominant modes according to energetic or temporal criteria. Several selection methods exist [5], the most simple one being a truncation based on the largest time constants [6]. The balanced realization [4] is based on truncation after transforming the dynamical system to a balanced form whose controllability and observability Grammians become diagonal and equal. The balanced truncation is characterized by two important properties: the asymptotic stability is preserved in the reduced model, and an error bound is available. The difficulty of the balanced truncation is that two matrix Lyapunov equations have to be solved, which is very expensive for large systems. The Branch Eigenmodes Reduction Method also uses a change of basis. It consists of solving a specific spectral problem called Branch problem, which concerns the advection-diffusion operator along

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

Nomenclature A B C Cp H h J M N n NR Ns P p q Re T T t u1, u2 UN V ! V x x1, x2 Y z

state matrix input matrix output matrix heat capacity J/kg K output matrix step height m cost function modal matrix DM order RM order number of Reynolds numbers number of time steps vector containing discrete pressure values pressure Pa dimension of output vector Reynolds number vector containing discrete temperature values temperature K time s velocity components m/s mean velocity m/s vector containing discrete velocity components local velocity vector m/s state vector for TRM and TCRM spatial coordinates m output vector state vector for FRM

dof FRM MIM POD PSO RM TCRM TRM

1355

degrees of freedom Fluid Reduced Model Modal Identification Method Proper Orthogonal Decomposition Particle Swarm Optimization Reduced Model Thermal Coupled Reduced Model Thermal Reduced Model

Greek symbols 3 maximum error l thermal conductivity W/m K m dynamic viscosity kg/m s n kinematic viscosity m2/s 4 heat flux density W/m2 P vector associated with heat advection term J vector associated with momentum advection term r density kg/m3 s mean quadratic error q vector of RM parameters Superscripts f fluid T Transposition sign t thermal Subscripts d relative to heat diffusion f fluid t thermal

Abbreviations DM Detailed Model

with a Steklov boundary condition. This particular boundary condition allows to obtain a basis adequate to handle nonlinearities [7]. The reduced basis is then obtained by selecting or amalgaming the most dominant modes according to a particular criterion (temporal, energetic, .). The Branch reduction method has notably been used in non-linear heat diffusion, especially for solving inverse heat transfer problems [8,9], but also for building low-order models for advection-diffusion problems [10]. Another approach uses reduced basis (RB) approximation based on Galerkin projection onto a Lagrange space [11]. The method relies on an underlying large sized (N dof) Finite Element model which is an approximation to the considered infinite-dimensional PDEs, and on samples of solutions of this reference model for different occurrences of parameters of interest (such as boundary conditions, loads, physical properties, etc). A Galerkin projection onto a n-dimensional ðn  NÞ Lagrange space associated with sampling parameters leads to a reduced model with n dof. Rigorous upper bound for the error between the reduced basis approximation and the reference solution is provided in the case of linear problems and for problems with polynomial nonlinearities. The method also provides an algorithm for determining the adequate data sampling to ensure that for some chosen outputs, this error remains bounded within a prescribed tolerance. Under the assumption of “affine in parameter” operators (or operators that can be written under approximate affine form), Offline computations involving N-sized elements (for the reduced basis construction) are decoupled from on-line computations involving only n-sized elements (calculus of outputs using the reduced model for any value of the sampling parameter vector). Such decoupling allows to use the method in the frame of multi-query (numerous

appeals to the solver, in the case of parametric studies for instance) or real-time (control-command for instance) processes. The RBGalerkin method has been applied to several types of problems, among which forced and natural heat convection problems [12,13]. A commonly-used approach in fluid mechanics is the Proper Orthogonal Decomposition method coupled with a Galerkin projection of the partial differential equations (POD-G). POD is a powerful method also known as Karhunen-Loève decomposition [14,15], and was introduced in fluid mechanics by Lumley in 1967 for the identification of coherent structures in turbulent flows [16]. So far, the method has been widely used in fluid mechanics and heat transfer (see for instance recent papers [17e20] to cite but a few). The principle of the method consists in obtaining some orthogonal modes from the analysis of covariance matrices built with some data covering space and time domains. The data signal is then expanded on a truncation of the basis formed by these spatial POD modes. The time-varying coefficients of this decomposition carry the system dynamics. Classically, POD requires solving an eigenvalue problem relative to a two-point spatial correlation matrix, with time averaging (under assumption of ergodic processes). However, for applications involving a large number of degrees of freedom, the spatial correlation matrix may become very large, and hence the computation of POD modes may become infeasible. To avoid this difficulty, Sirovich introduced the method of “snapshots” which has proved to be a powerful tool for the computation of the eigenmodes. The Snapshot POD consists to solve an eigenvalue problem relative to a two-point temporal correlation matrix, with space averaging, before computing the spatial eigenmodes. For more details the reader can see [21e23].

1356

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

Other approaches exist like the balanced POD [24,25] which combine the POD and concepts from balanced realization. One of the main characteristic of POD is its optimality in the sense of the minimization of the reconstruction error built on the signal energy. Indeed, very few modes are sometimes enough for capturing the main part of the energy of the data. We can quote the example of Ravindran [26] where 9 modes capture 99.9% of the energy in the case of a laminar flow over a backward-facing step channel. Liberge [27] in its application of flow around a cylinder at Re ¼ 200, notices that 7 modes are sufficient to reconstruct the solution with a 1% error. It should be noticed that POD is not a model reduction method, but a data compression method. To obtain a dynamic reduced model, the method has to be completed by a projection of the obtained basis on the equations that govern the system. This is usually a Galerkin-type projection. As recalled by Bui-Thanh et al. [28], it should be noted that the POD optimality property is only relative to the reconstructed data and does not concern the solution of a POD-based reduced order model, for which accuracy is therefore not guaranteed by a small POD reconstruction error. Moreover, standard POD-Galerkin based methods are not designed to target specific outputs. In addition, no information about the governing equations is included in the determination of the basis functions (such information appears in the Galerkin projection). Bui-Thanh et al. have proposed a goal-oriented, model-constrained optimization approach for model reduction that is able to target a particular output functional, span an applicable range of dynamic and parametric inputs, and respect the underlying governing equations of the system [28]. Their approach has therefore similarities with the method developed since many years at our lab, and which is presented in the following. The approach used in this paper is the Modal Identification Method (MIM). MIM uses concepts stemming from the automatics community and has been developed in heat transfer to identify low-order models in linear thermal diffusion processes from experimental data [29], and then has been extended to non-linear heat transfer [30,31]. An extension of the method to fluid flows has been proposed [32,33], the identification method being reformulated using the adjoint state method to compute the gradient of the cost function to be minimized [34]. The Modal Identification Method relies on two main steps:  The first stage is to define a structure of equations for the reduced model. Starting from the local partial differential equations describing the physics, a state space representation, which can arise from any type of space discretization, is considered as a general detailed model structure. Eventually, the reduced model equations are similar to those of the state space representation under modal form.  Once the structure is chosen, the next step is to identify the parameters that define the reduced model through minimization of a cost function that measures the difference between the outputs of the reduced model and the outputs of the system (detailed model or real process). The identification process therefore leans on the solution of an inverse problem of parameter estimation, where the data to be fitted are either coming from simulations performed with the detailed model or from in-situ measurements. It should be noted that no geometrical information and no physical properties appear explicitly in the MIM reduced models. All that information is in fact contained in the data, and therefore implicitly transmitted to the reduced model. When performing MIM, one obtains a measurement of the loss of information in the reduced model in comparison to the reference

one, but only for the input-output data used in the identification process. For linear problems, the superposition principle implies that the reduced model will be as good for any input as it is for the input used in the identification procedure (for time-dependent problems, the data input signal has to cover the whole frequency spectrum). In the case of non-linear problems, if the loss of information obtained for a “sufficiently rich” input signal is considered as satisfying by the user, it gives hope that the reduced model solutions will also be of good quality for other input values, but it does not give any assurance. This is similar to the approach proposed in [28] and to POD-Galerkin for instance. For these three methods, the optimality is relative to the used data. But for the first two approaches, the reconstruction error is measured under the constraint of the reduced model equations, while for the PODGalerkin no information about the model equations is included in the POD stage itself but only in the Galerkin projection stage. Of course, for all these reduction methods, one can always analyze the reduced model quality for any input by computing a posteriori the error induced by the reduced model compared to the reference one, if the solution of the reference model is available. An approach for computing estimates and bounds for the approximation error of POD reduced models has been proposed in [35]. Note that a formal comparison between the POD-Galerkin method and the MIM has shown that both can be formulated equivalently although the general ideas behind both methods are completely different [36]. Let us insist here on another important point: when performing MIM, as the reduced model equations are indirectly derived from PDEs, we ensure that we can use as data for the reduced model identification, fields coming from numerical codes, as long as those fields are good approximations of the solution of the related partial differential equations. It implies that the formulation of the reduced model does not depend on the numerical code which provides the data to be used. It does not even depend on the class of numerical discretization method. It is worth underlining that the Branch Eigenmodes Reduction Method [10], the Reduced Basis e Galerkin method [11] and the Goal-oriented, model constrained, optimization approach [28], all rely on a reference model or detailed model (DM) approximating the governing PDEs. In addition, it should be noted that while the first technique does not require DM simulations, the two other methods use numerical data (DM solutions). Whereas those three approaches require a reference model, POD-Galerkin based reduction methods and Modal Identification Method rely only on the governing PDEs and on numerical or measured data. It means that for these two approaches:  In the case of numerical data, only the solutions of a DM are used, the knowledge of this DM itself is not needed.  In the case of in-situ measured data, no DM and no solution of a DM are required. In our mind, this last point is crucial because in most cases, building a model that fits with the actual process is a very difficult task e sometimes impossible as regards to the desired accuracy e which usually requires having recourse to “trial and error” procedures and/or resolution of inverse problems to estimate properly some quantities such as thermophysical parameters or heat exchange coefficients for instance. Obtaining a reduced model that mimics the behavior of a reference model is desirable, but if the reference model fails at simulating the actual process, “experimental modeling”, that is building a low-order model from in-situ measured data (and possibly assumptions about mathematical description of physical phenomena), may become an interesting option. Such an approach allows to avoid dealing with approximations specific to classical discrete models (geometry, boundary

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

conditions, heat exchange coefficients, thermophysical parameters, .) and to take into account some minor nonlinearities whose impact would be difficult to assess. The Modal Identification Method, as well as the approach developed by Bui-Thanh et al. [28], shows features such as the use of optimization techniques, the goal-oriented framework (the reduced model is built for outputs of interest) and the construction of basis functions under the constraint of the reduced model equations. Furthermore, MIM does not require the knowledge of a Detailed Model, and is therefore well suited for experimental modeling. Such an approach has already been conducted in heat diffusion, for both linear [37] and non-linear problems [38], and also in turbulent forced convection [39]. In [39], the low-order model was relative to the advectionediffusion equation with constant mean velocity field, and was linking pollutant concentrations at a set of control locations to emission rates of pollutant sources. The aim was not to build a reduced model for the velocity field. Keeping in mind that our objective is to apply our approach to actual heated flows in order to perform experimental low-order modeling, the technique has been extended to the case of reduced models for steady NaviereStokes equations and tested on academic laminar flow configurations in a numerical framework [32,33]. The purpose of this paper is to evaluate the ability of the Modal Identification Method to reproduce and to predict the results obtained by a numerical code solving the steady NaviereStokes equations, coupled with the unsteady energy equation. For reasons of simplicity, and also because it is a very well known code among engineers, we chose to use the finite-volume code Fluent 6.3.26 [40] to provide us with the flow and temperature data. After giving the structure of the reduced model and introducing the method of identification, the approach is illustrated with a flow over a backward-facing step. Although this case has a very simple geometry, it features rich flow patterns like separation and recirculation bubbles. The literature offers many numerical and experimental studies on 2D steady incompressible flows over the backward-facing step. Two parameters characterize this flow: the Reynolds number and the channel expansion ratio. In the present study, the flow can be heated upstream before being convected. A similar backward-facing step problem, where the flow was heated by prescribed wall temperature, has been treated using a POD method [19]. Three types of reduced models (RMs) are presented, built and tested:  A Thermal Reduced Model (TRM), which allows to compute the temperature field for unsteady thermal boundary conditions and for a given constant Reynolds number.  A Fluid Reduced Model (FRM), related to steady fluid mechanics only, dedicated to the computation of the velocity field, for any constant Reynolds number in the range [100:800].  A Thermal Coupled Reduced Model (TCRM), which is weakly coupled to the FRM and allows to compute the unsteady temperature field for any constant Reynolds number in the same range and for any time-varying applied heat flux density.

1357

The MIM approach related to the TRM has been presented in [41], and the way to identify the FRM has been developed in [32,33]. The formulation and the identification of the TCRM are the main originality of the present paper, which is a logical continuation of our previous works. Concerning the organisation of the work presented hereafter, Section 2 presents the studied system and introduces the governing equations along with the corresponding boundary conditions. It also gives some elements of literature for detailed modeling comparison and validation. In section 3 the discrete forms of the governing equations are given with the reduced order model formulations. This section also gives some optimization elements for the identification of the reduced models. Next, section 4 proposes some numerical results of model reduction for the backward-facing step problem and compares those results to those of the reference model. Section 5 is then dedicated to conclusions and to future work. 2. The studied system This section describes the 2D incompressible laminar flow over a backward-facing step with an expansion ratio of 2.0, along with some numerical results of simulations. A schematic diagram of the considered geometry is shown in Fig. 1. It consists of a backwardfacing step in a duct where the step height is h ¼ 1 cm. The upstream height is also h ¼ 1 cm, hence the downstream height is 2h. The chosen fluid, assumed to be Newtonian and incompressible, is air with the following constant properties: dynamic viscosity m ¼ 1.81  10e5 kg/(m s), density r ¼ 1.205 kg/m3, heat capacity Cp ¼ 1005 J/(kg K), and thermal conductivity l ¼ 0.0262 W/(m K). The coordinate system is defined as shown schematically in this
!
! figure, where the e1 and e2 coordinate directions denote respectively the streamwise and transverse directions. The flow geometry was considered by Armaly et al. [42] who chose to define the Reynolds number Re, as follows:

Re ¼

UN 2h

n

(1)

where UN is the inlet mean velocity and n ¼ m=r is the kinematic viscosity of the fluid. The use of 2h as the characteristic length scale is attributed to the fact that the Reynolds number is based upon the hydraulic diameter of the inlet channel. The flow entering the channel is assumed to be fully developed and is described using a u1-velocity parabola for laminar flow. The inlet velocity profile is prescribed for x1 ¼ 4h; x2 ˛½h : 2h:

(

Ren u1 ð4h; x2 Þ ¼ 3 3 ðx2  hÞð2h  x2 Þ h u2 ð4h; x2 Þ ¼ 0:

(2)

! where (u1, u2) are the components of the velocity vector V in the
!
! ð e1 ; e2 Þ coordinate directions. In the present work, the outflow boundary condition that assumes that the flow is fully developed at the outlet section is used. This can be done when the outlet is

Fig. 1. The backward-facing step scheme.

1358

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

! ! located far enough downstream, assuming that v V =vx1 ¼ 0 . On all other boundaries a null velocity is prescribed. At the entrance of the pipe, the fluid temperature is TN ¼ 300 K and the outflow boundary condition leads to vT=vx1 ¼ 0. Next, a uniform heat flux varying in time 4(t) is located at x1 ˛½2h : 0 and x2 ¼ h, while the other walls are treated as adiabatic surfaces. The thermal initial condition is T(t ¼ 0) ¼ 300 K.

12

Armally Zang Erturk Present

10

8

2.1. The flow description

25 6

Barton [43] stated that when using an inlet channel upstream of the step, significant differences occur for low Reynolds numbers, however, they are localized in the sudden expansion region. In this study, to minimize its possible effect on the numerical solution, it has been decided to use an inlet length of 4h. Erturk [44] believed that the location of the outflow boundary is also very important for the accuracy of the numerical solution. The outflow boundary condition has been applied at the exit channel x1 ¼ 30h; x2 ˛½0 : 2h. The location of the outflow boundary has been chosen sufficiently far from the step to not affect the position of the recirculation zones. The Reynolds number at the inlet is chosen (through the inlet mean velocity UN) in order to lead to a stable flow. Many authors like Gresho et al. [45], Gartling [46] concluded that it is possible to obtain a steady solution for this flow until Re ¼ 800 but have not established when the flow becomes unstable. Other authors like Barkley et al. [47] continued this stability analysis up to Re ¼ 1, 500 and stated that the flow remains stable. As far as we are concerned, the range for the Reynolds number between 100 and 800 has been chosen to make sure the flow remains stable. 2.2. The reference model The computational fluid dynamics (CFD) software (Fluent 6.3.26) [40] has been used to carry out computations. Grid independence tests have been performed using several grid densities, and the reattachment location on the stepped wall has been used as the criterion. In the present study, 144, 247 nodes have been used for the computations. The experimental study provided by Armaly et al. [42] and numerical studies provided by Erturk [44] and Zang et al. [48] have been used to validate our model. For all Reynolds numbers used, there is a main recirculation region, whose length X1 increases with the Reynolds number (see Fig. 2). For a Reynolds number equal to 400, a second recirculation bubble appears attached to the upper wall of the channel between X2 and X3. Fig. 2 presents the evolution of adimensioned lengths X1/h, X2/h and X3/h with respect to the Reynolds number from 100 up to 800. The results follow the trend as observed in the literature. It has been observed that several benchmark solutions are available in the literature when the heat flux is applied in the stepped wall. However, to our knowledge, there is no benchmark solutions available when the heat flux is applied in the stepped wall just before the sudden expansion region. Once the fluid simulation methodology is validated, the thermal simulation in forced convection problem with the CFD code should be accurate. Note that the aim of the present paper is to illustrate the reduction methodology with its advantages, and not to give a complete study of the test problem. The required CPU time for flow simulations was about 2 h on a dual-core bi-processor AMD Opteron 2.2 GHz with 3 GB of RAM on a HP DL 145G2 data processing server. 3. Model reduction The aim is first to find a general structure for the reduced model equations. In Section 3.1 we present the governing equations and in Section 3.2 we consider a spatial discretization of these equations.

X2/h Erturk X2/h Present X3/h Erturk X3/h Present

20 15

4

10 5 400

2 100

200

300

400

500

500

600

700

600

700

800

Fig. 2. Size of the main recirculation region length X1/h as a function of the Reynolds number Re. The locations of the detachment point X2/h and the reattachment point X3/h as a function of the Reynolds number Re are shown in the inset.

Then, in Section 3.3, we develop a modal formulation from which is derived the reduced model structure. Eventually, the parameters of these equations are identified through an optimization problem in Section 3.4. 3.1. The governing equations The governing equations of the problem are the momentum equations, the continuity equation, and the energy equation:

! ! 1 ! ! vV þ ð V $VÞ V  nD V þ Vp ¼ 0 r vt

(3)

! V$ V ¼ 0

(4)

l vT ! DT ¼ 0 þ V $VT  r Cp vt

(5)

! where V is the velocity vector, p is the pressure and T is the temperature. 3.2. The discretized form of NaviereStokes and energy equations in the frame of MIM Rather than proposing hereafter a numerical scheme for solving the governing equations along with boundary conditions, this section consists in presenting a general matrix formulation of the discrete problem arising from equations (3)e(5). This formulation is a system of coupled ordinary differential equations of first order in time, also called “state space representation”. This state space representation is the starting point for writing down the reduced model structure. With this goal in mind, let us consider spatial meshes for velocity components, temperature and pressure, and let us call:  V ¼ ½ðVi Þ; i ¼ 1; .; 2Nf T ¼ ½ðu1 ; u2 Þi ; i ¼ 1; .; Nf T the vector of size 2Nf, containing the velocity components u1 and u2 at the discretization points;  T ¼ ½ðTi Þ; i ¼ 1; .; Nt T the vector of size Nt, containing the temperature value at the discretization points;  P ¼ ½ðPi Þ; i ¼ 1; .; Np T the vector of size Np, containing the pressure value at the discretization points.

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

 the Thermal Reduced Model (TRM) for the computation of the temperature field when 4(t) is varying, and for a constant velocity field;  the Fluid Reduced Model (FRM) relative to the velocity field of a stationary flow, for a given Reynolds number Re in the range [100:800];  the Thermal Coupled Reduced Model (TCRM) designed to calculate temperature for a varying flux 4(t) and for different values of Re (in the considered range). The TCRM is therefore weakly coupled to the FRM.

where Nf, Nt and Np depend on the chosen discretization method and on the mesh used. The state space equation relative to (3), (4) is (see [33] for details):

V_ ¼ Af V þ Q f JðVÞ þ Bf Re

(6)

where: !  Af V is the discrete form of the linear diffusion term D V ;  Q f JðVÞcontains quadratic terms of velocity components. One part corresponds to the discrete form of the non-linear ! ! advection term ð V $VÞ V while the other part comes from manipulation on the pressure gradient term Vp [33];  Bf Re is the discrete form of the flow boundary condition. Vector JðVÞ of size NJ ¼ Nf ð2Nf þ 1Þcontaining quadratic terms is written as:

2 2 ; V2Nf 1 V2Nf ; V2N V2N f 1 f

iT

(8)

where:  Ad T corresponds to the discrete form of the diffusion term DT;  Q t PðV; TÞ contains the discrete form of the transport term ! V $VT;  Bt is the input vector which applies the contribution of 4(t). The temperature field is coupled to the velocity field through the vector PðV; TÞ of size NP ¼ 2Nf Nt that represents the advection term:

h PðV; TÞ ¼ V1 T1 ; V2 T1 ; .; V2Nf T1 ; V1 T2 ; .; V2Nf T2 ; .; iT V1 TNt ; .; V2Nf TNt

(10)

where matrix EðVÞ depends only on V. Equation (8) then becomes:

(7)

Q f is a ð2Nf ; NJ Þ matrix. Note that though the components of matrices Af , Q f and Bf depend on the chosen discretization method and mesh, only the structure of equation (6) is of interest for our low-order model identification method. In a same manner, the space discretization of the energy equation (5) with associated boundary conditions leads to the following state space equation:

T_ ¼ Ad T þ Q t PðV; TÞ þ Bt 4ðtÞ

3.3.1. The form of the thermal reduced model for constant velocity The formulation of the TRM is deduced from equation (8). The temperature field is coupled to the velocity one through the vector PðV; TÞ that represents the advection term. Equation (9) can be modified by writing:

PðV; TÞ ¼ EðVÞT

h

JðVÞ ¼ V12 ; V1 V2 ; .V1 V2Nf ; V22 ; V2 V3 ; .; V2 V2Nf ; .;

1359

(9)

Q t is a ðNt ; NP Þ matrix. The temperature field is also dependent on thermal boundary conditions through the term Bt 4ðtÞ. 3.3. The form of the reduced models Our goal is now to define a reduced model structure for state equations (6) and (8). In classical modal reduction methods, a spectral problem associated with the state matrix (Af of equation (6) or Ad of equation (8)) is solved, and a selection of dominant eigenmodes is performed according to some temporal or energetic criteria. The Modal Identification Method also uses a reduced model expressed in a modal form with few degrees of freedom, but leans on a different approach. In fact, instead of solving a spectral problem and then selecting the dominant eigenmodes, the MIM aims at identifying a reduced set of eigenmodes from numerical or experimental data. Three kinds of reduced models are developed further:

T_ ¼ Adc T þ Bt 4ðtÞ

(11)

where matrix Adc ¼ Ad þ Q t EðVÞ takes into account both diffusion and advection phenomena. Adc is fixed for a constant velocity field. Equation (11) therefore underlines how the temperature field is only dependent on the thermal boundary condition when the velocity field remains constant. Let us also introduce the vector of outputs Y t of dimension q which enables to select a part of the computed field (q  Nt):

Y t ¼ CtT

(12)

Let us now consider the eigenvalue problem associated with matrix Adc on equation (11), and let M be the modal matrix of Adc . In order to find a modal formulation of equation (11) one introduces the change of variable:

T ¼ MX

(13)

that leads to:



X_ ¼ F 0dc X þ G0t 4ðtÞ Y t ¼ H 0t X

(14)

F 0dc ¼ M 1 Adc M is the diagonal matrix of the eigenvalues of matrix Adc . G0t ¼ M 1 Bt and H 0t ¼ C t M. The principle of model reduction using the Modal Identification Method is to build a reduced model having the same structure but b t the with a state vector x of size nt where nt  Nt and with Y corresponding approximation of the detailed model outputs:

(

x_ ¼ F dc x þ Gt 4ðtÞ b t ¼ Ht x Y

(15)

Note that because of the non-reciprocity of convective heat transfer (i.e. the operator Adc involved in (11) is not self-adjoint), not all eigenvalues and eigenvectors of matrix Adc are real numbers. Some form pairs of conjugate complex numbers [49]. Hence, the identified components of F dc , Gt and H t should also be either real or conjugate complex numbers. However, TRMs with real valued parameters have been identified in the present work. 3.3.2. The form of the fluid reduced model In a similar way than previously, we describe here how to obtain a form of reduced model from the velocity field equation (6). Let us consider q points among the Nf ones of the velocity mesh. An output

1360

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

Table 1 Vector q, cost function J ðqÞ, mean quadratic error s and maximum error 3 for each kind of reduced model. q is the number considered outputs. Ns and NR are respectively the number of time steps and Reynolds numbers used in the identification.

q

TRM

FRM

F dc , H t , Gt

Gf , H f , G f

q X Ns  X

J

Y tij 

2 bt Y ij

i¼1 j¼1

TCRM

2q X NR  X i¼1 k¼1

qffiffiffiffiffiffiffiffiffi

s

st ¼

sfu[ ¼

3

t 3t ¼ sup Y tij  Yb ij

J qNs

i;j

3fu[

F d , Gt , H t , Gt

Y fik 

qffiffiffiffiffiffiffiffiffi J u[ qNR

2 bf Y ik

u[

¼ J u1 þ J u2

q X Ns X NR  2 X bt Y tijk  Y ijk i¼1 j¼1 k¼1

[ ¼ 1; 2

b f ¼ sup Y fik  Y ik i;k

st ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J qNs NR



u[



3t ¼ sup Y tijk  Yb ijk t

i;j;k

vector Y f of dimension 2q is introduced, allowing to observe the velocity components u1 and u2 at each one of the q selected points:

mechanics problem (6) and the heat transfer problem (8) are weakly coupled.

Y f ¼ Cf V

V_ ¼ Af V þ Q f JðVÞ þ Bf Re

(20a)

where the matrix C f is the output matrix. In order to find a modal formulation, equation (6) is written by using the change of variable:

T_ ¼ Ad T þ Q t PðV; TÞ þ Bt 4ðtÞ

(20b)

V ¼ Mf Z

The observation equation (12) is also considered, allowing the selection of a set of temperature components:

(16)

(17)

where M f is the modal matrix of the state matrix Af . Therefore equations (6) and (16) become:

_ Z ¼ F 0f Z þ G0f JðZÞ þ G0f Re

(18)

Y f ¼ H 0f Z

where F 0f ¼ M 1 f Af M f is the diagonal matrix of the eigenvalues of matrix Af , and G0f is a ð2Nf ; NJ Þ matrix that applies the non-linear terms JðZÞ. In fact, JðVÞ ¼ JðM f ZÞ ¼ Lf JðZÞ where Lf is Next, a ðNJ ; NJ Þ matrix and therefore G0f ¼ M 1 f Q f Lf . G0f ¼ M 1 f Bf is the vector that applies Re, and the output equation is now expressed as Y f ¼ C f M f Z ¼ H 0f Z. Similarly as for the TRM, the reduction procedure is relative to a new state vector z, whose dimension nf is much lower than the one in the original model ðnf  2Nf Þ. Hence, the corresponding FRM has the following structure, where the steady flow is considered:

(

0 ¼ z þ Gf JðzÞ þ Gf Re b Y

f

(19)

¼ Hf z

JðzÞ is a vector of size nJ ¼ nf ðnf þ 1Þ=2. f

b the approximation of the This is the FRM structure, with Y detailed model output Y f . The components of the reduced matrices Gf , Gf and H f are now the parameters that describe the fluid reduced model. 3.3.3. The form of the thermal coupled reduced model The aim of this part is to find a model able to give the unsteady temperature outputs whatever the value of Re. To do so, the fluid

Y t ¼ Ct T We recall that using the change of variable (17) in equations (20a) and (16), one obtains the equations (18). Let us now define F 0d the diagonal matrix of eigenvalues of matrix Ad , and M t the matrix of corresponding eigenvectors. One has therefore F 0d ¼ M 1 t Ad M t . The change of variable:

T ¼ Mt X

(21)

is then considered. Both equations (17) and (21) are introduced in equation (20b). The advection term Q t PðV; TÞ now writes:

  Q t PðV; TÞ ¼ Q t P M f Z; M t X ¼ Q t Lt PðZ; XÞ where Lt is a ðNP ; NP Þ matrix. Let us call 0 1 0 G0t ¼ M 1 t Bt ; Gt ¼ M t Q t Lt and H t ¼ C t M t

The modal form of (20a,b), considering the steady state for the fluid mechanics point of view, is therefore written as:

8 0 ¼ Z þ G0f JðZÞ þ G0f Re > > < f 0 Y

¼ Hf Z

0 0 0 _ > > : X t ¼ F d X0 þ Gt PðZ; XÞ þ Gt 4ðtÞ

(22)

Y ¼ Ht X

The reduced model equations are going to be similar as equations (22) but with a fluid state vector z of dimension size nf  2Nf and a thermal state vector x of dimension size nt  Nt .

Fig. 3. The studied system.

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

1361

Table 2 Thermal reduced model identification results for the three test cases. Evolution of the cost function and mean quadratic error with respect to the TRM order nt. case 1 (q ¼ 135)

nt

1 2 3 4 5 6 7 8 9

case 2 (q ¼ 213)

st ðKÞ

J ðK2 Þ

st ðKÞ

J ðK2 Þ

st ðKÞ

2.23  10þ2 2.55  10þ1 1.14  10þ1 1.22  10þ0 1.19  101 3.70  102 2.47  102 1.67 3 10L2 1.67  102

0.11  10þ0 3.55  102 2.37  102 7.73  103 2.42  103 1.35  103 1.10  103 9.07 3 10L4 9.07  104

3.83  10þ3 1.23  10þ2 4.07  10þ1 5.08  10þ0 1.58  10þ0 5.57  101 1.49  101 3.02 3 10L2 2.98  102

3.46  101 6.18  102 3.56  102 1.26  102 7.03  103 4.17  103 2.16  103 9.71 3 10L4 9.63  104

7.11  10þ5 6.73  10þ4 1.99  10þ4 4.30  10þ3 9.56  10þ2 1.85  10þ2 5.10  10þ1 1.74 3 10D1 1.72  10þ1

1.81  101 5.57  102 3.03  102 1.41  102 6.64  103 2.92  103 1.53  103 8.95 3 10L4 8.90  104

The FRM for the whole velocity field is then expressed as:

(

0 ¼ z þ Gf JðzÞ þ Gf Re

3.4. Principle of identification

(23)

bf ¼ H z Y f

and the TCRM is written as:

(

x_ ¼ F d x þ Gt Pðz; xÞ þ Gt 4ðtÞ b t ¼ Ht x Y

(24)

where JðzÞ and Pðz; xÞ previously defined in (7) and (9) are here of size nJ ¼ nf (nf þ 1)/2 and nP ¼ nf nt respectively. b t are respectively the approximations of the velocity b f and Y Y

and temperature output vectors Y f and Y t . Because of the one-way coupling between temperature and velocity, the identification procedure of the TCRM is performed through two stages:  One first identifies the fluid reduced model equation (23) defined by the matrices Gf , Gf and H f . This enables to compute the reduced fluid state vector z for a given Reynolds number.  In the second stage, one identifies matrices F d , Gt , Gt and H t of the TCRM which also depends on the reduced fluid state vector z.

NB: only the fluid reduced state z of size nf appears in the thermal state equation of the TCRM, and not the velocity field of size 2Nf.

The previous developments enabled us to establish the structure of the reduced models. The identification approach does not require the solution of any eigenvalue problem. Once the reduced model structures are established, the aim is to identify the parameters q (i.e. the components of involved vectors and matrices) for each reduced model. The procedure leans on the minimization of the mean squared discrepancy J ðqÞ based on the difference between the response of the detailed model Y and b ¼ Y b ðqÞ given by the reduced model when the same the outputs Y input signal is applied. The vector q of parameters to be identified for each kind of reduced model is defined in Table 1. The mean quadratic error s associated with cost function J ðqÞ, and the maximum error 3, allow us to assess the quality of the reduced model identification. The expressions of J ðqÞ, s and 3 are also given for each kind of reduced model in Table 1. Note that the errors are defined either on the temperature variable, or on separated velocity components. For the TRM defined by equation (15), the parameters q to be identified are the matrices F dc , H t and Gt . This dynamic reduced model is associated to a time-varying input 4(t). The cost function to be minimized integrates the squared differences for the temperature at Ns time steps and at the q selected locations. For the static FRM defined by equation (19), the parameters q to be identified are the matrices Gf , H f and Gf . The cost function to be minimized integrates the squared differences for both components of the velocity at the q selected locations and for the NR Reynolds numbers used.

2

DM TRM

304

1.5 Position x2/h

303

2

300 Flux W/m

Temperature (K)

case 3 (q ¼ 144, 247)

J ðK2 Þ

302

250 200

0.5

150 100 0

301 0

1

5

10

10

15

20

20

25

30

30

Time (s) Fig. 4. Temperature evolution for case 2 at the location of the maximum error. The test heat flux (W/m2) is shown in the insert.

0 300

300.5

301

301.5 302 Temperature (K)

302.5

303

Fig. 5. Temperature profile (K) for case 1 at time t ¼ 15.3 s and x1/h ¼ 6 for both the detailed model (solid line) and the thermal reduced model of order 8 (dashed line).

1362

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

Fig. 6. Temperature fields (K) for case 3 for both the detailed model (top) and the reduced model (bottom) at t ¼ 16.6 s.

Fig. 7. The difference (K) between temperature fields obtained with the detailed model and the reduced model of order 8 at t ¼ 16.6 s.

Next, for the TCRM defined by equation (24), the parameters to be identified are the matrices F d , Gt , H t and Gt . This dynamic reduced model is associated to a time-varying input 4(t) and to the Reynolds number. The cost function to be minimized integrates the squared differences for the temperature at Ns time steps, at the q selected locations, and for the NR Reynolds numbers used. The optimization problem consists in finding the optimal vector q such that:

q ¼ arg min½J ðqÞ

(25)

For both TRM and FRM, the minimization has been performed through a BFGS quasi-Newton [50] approach. In spite of the (a priori) non-convexity of the cost function to be minimized, satisfying results have been achieved with such a deterministic optimization method. For the TCRM, it has been noticed that results were very sensitive to the initial guessed values of unknown parameters introduced in the minimization algorithm. Consequently, one had to try several initializations to get a cost function that decreased significantly. This is the reason why we turned to stochastic methods and decided to use the Particle Swarm Optimization (PSO) [51]. 3.4.1. Principle of the quasi-Newton type method 0 From a first guessed set of components q , a series defined by:

qi ¼ qi1 þ ai di

(26)

is built, where i is the iteration index, a is the step size and d is the search direction computed from the cost function gradient. This gradient VJ ðqÞ is computed through the adjoint state method as described in [32,34]. As the relations between the outputs and the states are linear (see equations (15) and (19)), matrices H t and H f are not included in the parameters vector q. For a given iteration, q is fixed, the state i

i

Table 3 Validation of thermal reduced model with an input test signal 4(t), for 3 different data sets. x1/h and x2/h are the coordinates of the maximum error location. Case

Order nt st ðKÞ

1 (q ¼ 135) 8 2 (q ¼ 213) 8 3 (q ¼ 144, 247) 8

vector (x or z) is computed and matrix H t or H f is calculated by the linear least squares method. The minimization of the cost function J is first computed for order n ¼ 1. Using the identified RM of order 1 as initial guess, the RM of order 2 is then identified. The order n is then increased until a prescribed accuracy between the DM and the RM is obtained or until J cannot be decreased in a significant manner. The general algorithm for the reduced model identification can be found in [33].

3ðKÞ

x1/h

x2/h Time (s)

8.06  104 3.20  103 6 0 2 0 8.99  104 7.14  103 8.14  104 1.20  102 1.88 1

15.3 17 16.6

3.4.2. The principle of the Particle Swarm Optimization method (PSO) This algorithm was first described in 1995 by James Kennedy and Russell C. Eberhart. The PSO algorithm is based on a metaphor of social interaction. Let D ¼ dimðqÞ the number of the unknown parameters. PSO uses a swarm of S particles searching in the unknown parameters space of dimension D. Each one of the S particles is therefore localized in the search space by D coordinates forming the vector q of parameters to be estimated. Each particle has also semi-stochastic velocities dq. Particles communicate in order to collectively reach the optimal position q. In addition to their tendency to move in each direction according to a fraction of their current velocity, the individual particles are stochastically drawn toward the positions of their own previous best performance pi and the best previous performance of their neighbors pg . PSO is an order 0 optimization method. Such methods do not require computation of the cost function gradient, only evaluations of the cost function are needed. 0 A population of particles is initialized with random positions q and a function J ðqÞ is evaluated for each particle. The algorithm in pseudo-code follows:

Table 4 Evolution of the cost function value J , the mean quadratic errors sui , i ¼ 1, 2 and the maximum errors 3ui , i ¼ 1, 2 with respect to the reduced model order nf. Order nf J (m2/s2) 1 2 3 4 5 6 7

4.96  10þ2 7.22  10þ1 1.33  10þ1 2.18  10þ0 4.09  101 5.51  102 4.80 3 10L3

su1 (m/s)

3u1 (m/s)

su2 (m/s)

3u2 (m/s)

2.92  102 1.09  102 4.63  103 1.84  103 7.80  104 2.79  104 8.15 3 10L5

1.31  101 4.95  102 2.36  102 1.28  102 3.71  103 1.44  103 4.91 3 10L4

1.62  104 9.92  105 6.44  105 4.06  105 2.64  105 1.58  105 8.54 3 10L6

4.22  102 1.89  102 9.02  103 5.52  103 2.41  103 1.03  103 3.39 3 10L4

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

(1) In the first test case a specific set of outputs at q ¼ 135 points located at the x1/h ¼ 6 cross-section is used; (2) The second test case involves 9 cross-sections located between x1/h ¼ e2 and the x1/h ¼ 6 (q ¼ 213 points indicated by dashed lines in Fig. 3); (3) The last case involves the whole temperature field (q ¼ 144, 247).

Initialize population For i ¼ 1 to S i

if J ðq Þ < J ðpi Þ then i

pi ¼ q end if pg ¼ minðneighborsðpi ÞÞ For d ¼ 1 to D i

i

i

i

dq ¼ cdq þ l rand1 ðpg  q Þ þ l rand2 ðpi  q Þ i

i

1363

i

q ¼ q þ dq Next d Next i

where rand1 and rand2 are 2 random numbers taken from a uniform distribution in [0:1]. The stopping criterion depends on the type of problem to be solved. Usually the algorithm is run for a fixed number of function evaluations or until a specified error bound is reached. The parameters c and l can act on the compromise of exploration of the space of the unknown parameters and of exploitation of an optimal position. The choice of parameters remains largely empirical. A complete analysis of the algorithm has been made by Clerc and Kennedy [52]. Standard PSO is easy to code, we have therefore written our own PSO algorithm. This home-made code uses parameters c ¼ 0.729 and l ¼ 1.494. This set has been determined by Clerc [53] and tested in [54]. A 20 particles swarm with a circular neighborhood of size 3 has been used. These values are recommended by Clerc [51]. 4. Model reduction on the studied system To illustrate the methodology presented in Section 3, we present in this section some computational results of model reduction for the flow over the backward-facing step. This section is divided into three parts:  The first one is relative to the thermal reduced model: the velocity field within the pipe is given through a fixed Reynolds number at the inlet, and a time-varying heat flux is applied upstream the step;  The second part is relative to the fluid reduced model: only the velocity components are of interest, and so for different Reynolds numbers;  Next, the third part of this section is relative to the thermal coupled reduced model: the inputs are the Re value and the time-varying heat flux, and the outputs are a set of chosen resulting temperatures.

4.1. The thermal reduced model (TRM) for a constant velocity field 4.1.1. The thermal reduced model identification In this part, the velocity field is fixed at Re ¼ 500 and a heat flux density 4(t) is applied on the part of the inlet duct just upstream the step on a length of 2h, as shown in Fig. 3. It has been verified that for 4 ¼ 300 W/m2, the maximum temperature in the channel leads to a Richardson number Ri ¼ Gr/Re2 about 0.1. This allows us to make the assumption of forced convection. Temperatures are computed with the DM (144,247 nodes) when an increasing step of 300 W/m2 is applied during Ns ¼ 300 time steps of 0.1 s. The obtained timevarying temperatures, from which is subtracted the inlet temperature TN ¼ 300 K, are then used as data for the model identification. Outputs of TRM are therefore temperature variations with respect to TN. Three cases, each one involving a specific set of outputs, are considered:

Table 2 summarizes the identification results for the TRMs of order nt ¼ 1 to 9, and gives the evolution of the cost function J and mean quadratic error st as defined in Table 1.

Table 5 Mean and maximal velocities of both components u1 and u2 for data associated with Re ¼ 100, 200, ., 800.

u1 u2

Max (m/s)

Mean (m/s)

8.76  101 1.04  101

1.73  101 4.95  103

Table 6 Evolution of errors sui and 3ui , i ¼ 1, 2 for the seven considered validation tests, i.e. for Reynolds numbers from 150 to 750 by steps of 100. Re

su1 (m/s)

150 250 350 450 550 650 750

3.98 2.06 1.25 5.81 5.77 1.35 2.57

      

su2 (m/s)

104 104 104 105 105 104 104

2.23 1.20 7.47 3.55 2.74 6.78 1.26

      

104 104 105 105 105 105 104

3u1 (m/s) 2.27 1.05 5.34 2.82 2.56 5.00 1.00

      

103 103 104 104 104 104 103

3u2 (m/s) 1.44 7.15 3.82 1.81 1.47 3.19 6.09

      

103 104 104 104 104 104 104

Table 7 Length of the main recirculating region X1/h obtained with on the one hand the detailed model and, on the other hand, the identified fluid reduced model of order 7. Re

(X1/h)DM

(X1/h)FRM

Error %

150 250 350 450 550 650 750

3.962 5.872 7.507 8.843 9.885 10.71 11.44

4.038 5.841 7.450 8.770 9.809 10.67 11.39

1.917 0.514 0.748 0.819 0.767 0.438 0.459

Table 8 Detachment locations of the recirculation bubble on the upper wall, X2/h. Re

(X2/h)DM

(X2/h)FRM

Error %

450 550 650 750

7.769 8.190 8.663 9.145

8.215 8.503 8.982 9.423

5.733 3.821 3.682 3.035

Table 9 Reattachment locations of the recirculation bubble on the upper wall, X3/h. Re

(X3/h)DM

(X3/h)FRM

Error %

450 550 650 750

11.633 14.478 17.011 19.361

11.290 14.281 16.727 19.069

2.951 1.359 1.668 1.510

1364

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

Fig. 8. Horizontal velocity (m/s) field for both the detailed model (top) and the reduced model of order 7 (bottom), in the Re ¼ 550 case.

Fig. 9. Vertical velocity (m/s) field for both the detailed model (top) and the reduced model of order 7 (bottom), in the Re ¼ 550 case.

According to Table 2, TRM of order 8 is chosen because the cost function J does not decrease significantly for order 9. Note that for the three cases, st is less than 103 K.

4.2. The fluid reduced model (FRM) 4.2.1. The fluid reduced model identification In the present work, as we consider having no prior information on the effect of the Reynolds number on flow patterns, the intuitive choice is to use data velocity fields corresponding to a regular partition of the interval [100:800]. Therefore, 8 velocity fields are computed for Reynolds numbers from 100 to 800 by steps of 100 using the detailed model involving 144, 247 mesh

2

2

1.5

1.5

x2/h

x2/h

4.1.2. The thermal reduced model validation The purpose of validation is to verify that the TRM is able to reproduce the behavior of the detailed model when it is submitted to another heat flux density than the step used for the reduced model identification. As the heat transfer problem is linear with respect to 4 and because we have used a 4 step function for the identification, the TRM is indeed expected to work with any other 4 input signal. The validation is performed with the test heat flux density 4(t) shown in insert of Fig. 4. Fig. 5 gives for case 1 the temperature profile along the x1/h ¼ 6 cross-section at t ¼ 15.3 s, when the error is maximum, and Fig. 4 shows the temperature evolution in case 2, at the point where the error is maximum. For case 3, Fig. 6 shows the

temperature fields obtained with both detailed and reduced models, at t ¼ 16.6 s. In order to better see the difference between both temperature fields, the difference in absolute value is plotted on Fig. 7. Table 3 summarizes all the validation results and shows that the temperatures computed with each TRM are very close to the DM ones, with a low st . The maximum error 3, with its position in space and time, is also given.

1 x1/h=6 x1/h=6 x1/h=14 x1/h=14 x1/h=30 x1/h=30

0.5

0

0

0.1

0.2 0.3 0.4 Horizontal velocity (m/s)

DM FRM DM FRM DM FRM

0.5

Fig. 10. Comparison of the detailed model with the fluid reduced model of order 7: horizontal velocity u1 for Re ¼ 550 at locations x1/h ¼ 6, x1/h ¼ 14 and x1/h ¼ 30.

1

0.5

0 -0.04

x1/h=6 x1/h=6 x1/h=14 x1/h=14 x1/h=30 x1/h=30 -0.03

DM FRM DM FRM DM FRM

-0.02 -0.01 0 Vertical velocity (m/s)

0.01

0.02

Fig. 11. Comparison of the detailed model with the fluid reduced model of order 7: vertical velocity u2 for Re ¼ 550 at locations x1/h ¼ 6, x1/h ¼ 14 and x1/h ¼ 30.

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

1365

Fig. 12. The studied system.

4.2.2. The fluid reduced model validation The aim here is to validate the FRM finding out if it is able to reproduce with accuracy the outputs of the DM. Note that these tests are performed with Reynolds numbers 150, 250, ., 750 that have not been used in the identification procedure. Table 6 presents the validation results. We can see that the FRM of order 7 can reproduce the DM results with satisfying accuracy. Tables 7, 8 and 9 confirm the validity of the FRM of order 7 by showing respectively that it is able to reproduce the length X1/h of the main recirculating region and the locations of detachment X2/h and reattachment X3/h of the recirculation bubble on the upper wall for Re  400. Figs. 8 and 9 show respectively the velocity components u1 and u2, computed by the DM and the FRM of order 7. We note a good agreement on the velocity fields given by both models. For the Re ¼ 550 case, we have outlined respectively in Figs. 10 and 11, the profiles of the horizontal velocity u1 and the vertical velocity u2 at x1/h ¼ 6, x1/h ¼ 14 and x1/h ¼ 30. We can see that the profiles computed with the FRM of order 7 are in very good agreement with those calculated with the DM.

4.3. The thermal coupled reduced model (TCRM) 4.3.1. The thermal coupled reduced model identification In this section the aim is to identify a TCRM able to give the temperatures whatever the value of the Reynolds number and whatever the heat flux density 4(t) applied at the heater. The reduced fluid state z which depends on the Reynolds number is solution of the state equation of the order 7 FRM identified in Section 4.2. The TCRM is weakly coupled to the FRM through this reduced fluid state z. In order to form a set of representative data of the dynamics of our system, the temperature field is calculated with the DM (q ¼ 144, 247) for 6 Reynolds numbers from 300 to 800 by steps of

304 DM TCRM DM TCRM

Point Point Point Point

A A B B

303 Temperature (K)

points. All nodes whose x1 coordinate is included in the range [e2h:30h] are included in the reduction process. This gives a total of 139, 677 nodes. Using these fields as data, a series of seven FRM of orders nf ¼ 1 to 7 can be obtained through the identification procedure. Corresponding results are summarized in Table 4. This table gives the evolution of the cost function J with respect to the reduced model order, as well as the mean quadratic errors su1 and su2 between the detailed model and the fluid reduced model solutions for u1 and u2 components respectively, and also the maximum errors 3u1 and 3u2 . We note that the cost function value decreases with respect to the reduced model order nf, and that the best minimization is obtained for order 7, for which the cost function J and associated errors sui and 3ui are the lowest, and actually low when compared to the mean and maximal velocities presented in Table 5.

302

301

300

0

5

10

15 Time (s)

20

25

30

Fig. 13. Temperature evolution at 2 points A and B in the x1/h ¼ 6 cross-section. Comparison detailed model/thermal coupled reduced model of order 7.

2

DM TCRM

Table 10 TCRM identification results using FRM of order nf ¼ 7. Order nt

Case 1 (q ¼ 135) J (K )

1 2 3 4 5 6 7 8 9 10

9.90 5.49 2.67 7.67 3.56 1.48 6.41 5.20 4.99 3.03

         

Case 2 (q ¼ 213)

st (K)

2

10þ3 10þ2 10þ2 10þ1 10þ1 10þ1 10þ0 10þ0 10þ0 10þ0

2.85 6.71 4.68 2.51 1.71 1.10 7.25 6.53 6.40 4.99

         

st (K)

2

J (K ) 101 102 102 102 102 102 103 103 103 103

3.78 9.55 1.05 4.50 2.35 1.21 6.42 3.63 2.23 8.40

         

10þ4 10þ3 10þ3 10þ2 10þ2 10þ2 10þ1 10þ1 10þ1 10þ0

0.44 0.22 7.37 4.84 3.50 2.51 1.83 1.37 1.08 8.40

         

10þ0 10þ0 102 102 102 102 102 102 102 103

Position x 2/h

1.5

1

0.5

0 300

300.5

301

301.5 302 Temperature (K)

302.5

303

Fig. 14. Temperature profile along the x1/h ¼ 6 cross-section at t ¼ 30 s. Comparison detailed model/thermal coupled reduced model of order 7.

1366

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

DM TCRM DM TCRM DM TCRM

Temperature (K)

320

Point Point Point Point Point Point

C C D D E E

315

Fig. 17. Comparison of the detailed model and the thermal coupled reduced model of order 10: temperature profile at x1/h ¼ [1, 2, 3, 4, 5, 6] and t ¼ 30 s. 310

Table 11 Validation of the thermal coupled reduced model with the input test signal 4(t), for 2 different data sets. x1/h and x2/h are the coordinates of the maximum error location.

305

300

0

5

10

15 Time (s)

20

25

30

Case

Order nt

st (K)

3(K)

x1/h

x2/h

Time (s)

case 1 (q ¼ 135) case 2 (q ¼ 213)

7 10

6.16  103 7.27  103

1.83  102 3.85  102

6 1

1.075 0

16.9 15.4

Fig. 15. Evolution of the temperature computed with the detailed model and the thermal coupled reduced model of order 10 at 3 points C,D,E.

100. A heat flux density 4(t) ¼ 300 W/m2 is applied during 300 time steps of 0.1 s for each Reynolds number. The obtained temperature fields, from which is subtracted the inlet temperature TN ¼ 300 K, are then used as data for the identification. Two cases, each one with a specific set of outputs, are considered: (1) the cross-section of 135 points located at x1/h ¼ 6. (2) the zone of 213 points indicated by dashed lines in Fig. 12 (e2  x1/h  6). Table 10 summarizes the identification results for TCRMs of order nt ¼ 1 to 10. The order of the TCRM has been chosen such as st ¼ 102 K. Table 10 shows that for case 1 the TCRM of order 7, and for case 2 the TCRM of order 10, give good approximations of the DM data in both cases. 4.3.2. The thermal coupled reduced model validation The purpose of the validation is to verify that the TCRM is able to reproduce the behavior of the DM when it is subjected to the test flux 4(t) (see insert of Fig. 4) and for a value of Reynolds number different than those used in the identification process. Here this test has been performed for Re ¼ 550.  Validation of the thermal coupled reduced model in case 1 Fig. 13 shows the temperature evolution at two points A and B (defined in Fig. 12), for the DM and the TCRM of order 7. Fig. 14 shows the temperature profile at x1/h ¼ 6 and t ¼ 30 s computed by both models.  Validation of the thermal coupled reduced model in case 2 Fig. 15 shows the temperature evolution computed with the DM and the TCRM of order 10, for 3 points C, D and E located within the 2  x1/h  6 zone (defined in Fig. 12). In Figs. 16 and 17 are respectively plotted the temperature profiles at x1/h ¼ 2, 1, 0 and at x1/h ¼ 1, 2, 3, 4, 5, 6, for time t ¼ 30 s.

Fig. 16. Comparison of the detailed model and the thermal coupled reduced model of order 10: temperature profile at x1/h ¼ [e2, e1, 0] and t ¼ 30 s.

It can be seen that the profiles computed with both the detailed and the reduced model are in very good agreement. Table 11 summarizes the validation results and shows that the temperatures computed with each TCRM are very close to the DM ones, with low st . The maximum error 3, with its position in space and time, is also given.

5. Conclusion This study is an extension of the Modal Identification Method for forced convection flows, and proposes a formulation of the reduced model that includes a one-way coupling between velocity and temperature. The case of 2D, laminar, incompressible, steady flow of a Newtonian fluid, with transient heat transfer, has been considered to illustrate the method. Starting from the structure of the partial differential equations of the physical system, a structure of reduced model in a modal basis is defined. Three reduced model formulations have been derived from local governing equations:  One is for the transient energy equation, for a given steady velocity field corresponding to a fixed Reynolds number. In this case, the Thermal Reduced Model input is a time-varying heat flux density and the outputs are a set of temperatures (the whole field or a part of it). The identified reduced order model has been validated with a test heat flux input different from the one used in the identification process. Three very distinct cases have been tested. The first and the second cases were more original since we used specific sets of outputs located on a single cross-section or at specific nodes in 9 cross-sections. The third was more classical since we considered the whole temperature field. All these results are in very good agreement with those computed with the detailed model.  Another one is for the NaviereStokes and continuity equations in steady state. In this case, the Fluid Reduced Model can be used to compute the velocity field whatever the value of the Reynolds number within the range [100:800] for which it has been built. The identified reduced model has been validated with Reynolds numbers that were not used in the reduction process. We have found that the Fluid Reduced Model was able to predict with accuracy the velocity field. It has also been validated comparing the results for the length of the first recirculation bubble, and the locations of the detachment and reattachment of the second recirculation bubble on the upper wall.  The last one is for the transient energy equation coupled with the momentum equation, whatever the value of Reynolds

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368

number governing the steady velocity field and whatever the heat flux density. This Thermal Coupled Reduced Model uses the reduced state vector computed with the steady Fluid Reduced Model. The results show that the temperatures computed with the TCRM are very close to the reference model one, for two distinct cases of outputs. For each reduced model, the identification algorithm consists in minimizing a quadratic function based on the discrepancy between the outputs given by the reduced model and the ones given by the original detailed model. For the fluid reduced model and the thermal reduced model, a gradient method can be used. For the thermal coupled reduced model, due to the sensitivity of results to the initial guessed unknown parameters, we used the stochastic zero-order Particle Swarm Optimization method. While the CPU time required for the detailed simulations was about 2 h on the data processing server, the CPU time required for running the reduced order models were about 0.01 s on the same server, where most of this time was dedicated in reading the data (matrices involved in the reduced order models). This ratio of speed-up computation of about 106 allows us to plan the use of such reduced models for on-line control purposes. Our current works deal with the numerical implementation of such on-line control algorithms based on thermal reduced models. These algorithms include Kalman filters to reconstruct on-line the state from few measurements. Next work will concern the use of the presented thermal coupled reduced models in order to control, for instance, some temperatures by acting on both the thermal and the fluid boundary conditions. Moreover, although numerical data are used in the present work, it should be noted that the structure of reduced model equations developed in this paper is also valid for in-situ measured data, of course under the assumption that the process is well described by the governing PDEs from which the reduced model equations are derived. Our long-term objective is to control experimental thermal systems through such low-order models.

References [1] W.H.A. Schilders, H.A. Van der Vorst, J. Rommes, Model Order Reduction: Theory, Research Aspects and Applications, Volume 13 of Mathematics in Industry. Springer, 2008. [2] D.J. Lucia, P.S. Beran, W.A. Silva, Reduced-order modeling: new approaches for computational physics. Progress in Aerospace Sciences 40 (1e2) (2004) 51e117. [3] A.C. Antoulas, An overview of approximation methods for large-scale dynamical systems. Annual Reviews in Control 29 (2) (2005) 181e190. [4] B.C. Moore, Principal component analysis in linear systems: controllability, observability and model reduction. IEEE Transactions on Automatic Control AC-26 17 (1981) 32. [5] M.T. Ben Jaafar, R. Pasquetti, D. Petit, Model reduction for thermal diffusion: application of the Eitelberg, Marshall and aggregation methods to a heat transmission tube model. International Journal of Numerical Methods in Engineering 29 (1990) 599e617. [6] S.A. Marshall, Control. An approximation method for reducing the order of linear system (1966) 642e643. [7] O. Quéméner, A. Neveu, E. Videcoq, A specific reduction method for the branch modal formulation: application to a highly non-linear configuration. International Journal of Thermal Sciences 46 (9) (2007) 890e907. [8] E. Videcoq, O. Quéméner, M. Lazard, A. Neveu, Heat source identification and on-line temperature control by a branch eigenmodes reduced model. International Journal of Heat and Mass Transfer 51 (19e20) (2008) 4743e4752. [9] E. Videcoq, O. Quéméner, W. Nehme, A. Neveu, Real time heat sources identification by a branch eigenmodes reduced model. In 6th International Conference on Inverse Problems in Engineering, Theory and Practice, Dourdan (France), and in Journal of Physics: Conference Series 135 (freely Available on Line, Paper N 012101), 2008. [10] F. Joly, O. Quéméner, A. Neveu, Modal reduction of an advection-diffusion model using a branch basis. Numerical Heat Transfer, Part B: Fundamentals 53 (2008) 1e20. [11] G. Rozza, D.B.P. Huynh, A.T. Patera, Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial

[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]

1367

differential equations: application to transport and continuum mechanics. Archives of Computational Methods in Engineering 15 (3) (2008) 229e275. N.C. Nguyen, K. Veroy, A.T. Patera, Handbook of Materials Modeling, Chapter Certified Real-time Solution of Parametrized Partial Differential Equations pp. 1529e1564 (2005). http://www.springerlink.com/content/rq1824h2h1q6444n. G. Rozza, N.C. Nguyen, A.T. Patera, S. Deparis, Reduced basis methods and a posteriori error estimators for heat transfer problems, in: Proceedings of the ASME HT 2009 Summer Conference, Heat Transfer Confererence. ASME, New York, 2009. K. Karhunen, Zur spektraltheorie stochastischer prozesse. Annales Academiae Scientiarum Fennicae, Series A1, Mathematica-Physica 34 (1) (1946) 7. M. Loève, Probability Theory. Van Nostrand, Princeton, NJ, USA, 1955. J.L. Lumley, The structure of inhomogeneous turbulent flows, in: A.M. Yaglom, V.I. Tartarski (Eds.), Atmospheric and Radio wave Propagation. Nauka, Moscow, 1967, pp. 166e178. M. Bergmann, L. Cordier, Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models. Journal of Computational Physics 227 (16) (2008) 7813e7840. M. Bergmann, C.-H. Bruneau, A. Iollo, Enablers for robust pod models. Journal of Computational Physics 228 (2009) 516e538. D. Alonso, A. Velazquez, J.M. Vega, A method to generate computationally efficient reduced order models. Computer Methods in Applied Mechanics and Engineering 198 (33e36) (2009) 2683e2691. D. Alonso, A. Velazquez, J.M. Vega, Robust reduced order modeling of heat transfer in a back step flow. International Journal of Heat and Mass Transfer 52 (5e6) (2009) 1149e1157. L. Sirovich, Turbulence and the dynamics of coherent structures part i: coherent structures. Quartely Applied Mathematics 45 (October 1987) 561e571. L. Sirovich, Turbulence and the dynamics of coherent structures part ii: symetries and transformations. Quartely Applied Mathematics 45 (3) (October 1987) 573e582. L. Sirovich, Turbulence and the dynamics of coherent structures part iii: dynamics and scaling. Quartely Applied Mathematics 45 (1987) 583e590. C.W. Rowley, Model reduction for fluids using balanced proper orthogonal decomposition. International Journal of Bifurcation and Chaos 15 (3) (March 2005) 997e1013. K. Willcox, J. Pereire, Balanced model reduction via the proper orthogonal decomposition. AIAA Journal 40 (11) (2002) 2323e2330. S.S. Ravindran, A reduced-order approach for optimal control of fluids using proper orthogonal decomposition. International Journal for Numerical Methods in Fluids 34 (2000) 425e448. E. Liberge, Modèles réduits obtenus par la méthode de POD-Galerkin pour les problèmes d’interaction fluide structure. PhD thesis, Université de Poitiers, 2008. T. Bui-Thanh, K. Willcox, O. Ghattas, B. van Bloemen Waanders, Goal-oriented, model-constrained optimization for reduction of large-scale systems. Journal of Computational Physics 224 (2) (2007) 880e896. R. Pasquetti, D. Petit, Analyse modale d’un processus de diffusion thermique: identification par thermographie infrarouge. International Journal of Heat and Mass Transfer 31 (3) (1988) 487e496. M. Girault, D. Petit, Identification methods in nonlinear heat conduction. part i: model reduction. International Journal of Heat and Mass Transfer 48 (January 2005) 105e118. O. Balima, Y. Favennec, D. Petit, Model reduction for heat conduction with radiative boundary conditions using the modal identification method. Numerical Heat Transfer, Part B 52 (2007) 107e130. O. Balima, Y. Rouizi, Y. Favennec, D. Petit, Reduced modelling through identification on 2D incompressible laminar flows. Inverse Problems in Science and Engineering 17 (3) (2009) 303e319. Y. Rouizi, Y. Favennec, J. Ventura, D. Petit, Numerical model reduction of 2D steady incompressible laminar flows: application on the flow over a backwardfacing step. Journal of Computational Physics 228 (6) (2009) 2239e2255. Y. Favennec, M. Girault, D. Petit, The adjoint method coupled with the modal identification method for nonlinear model reduction. Inverse Problems in Science and Engineering 14 (2006) 153e170. C. Homescu, L.R. Petzold, R. Serban, Error estimation for reduced-order models of dynamical systems. SIAM Journal on Numerical Analysis 43 (4) (2005) 1693e1714. O. Balima, Y. Favennec, M. Girault, D. Petit, Comparison between the modal identification method and the POD-galerkin method for model reduction in nonlinear diffusive systems. International Journal of Numerical Methods in Engineering 67 (2006) 895e915. E. Videcoq, D. Petit, A. Piteau, Experimental modelling and estimation of time varying thermal sources. International Journal of Thermal Sciences 42 (3) (2003) 255e265. M. Girault, E. Videcoq, D. Petit, Estimation of time-varying heat sources through inversion of a low order model built with the modal identification method from in-situ temperature measurements. International Journal of Heat and Mass Transfer 53 (1e3) (2010) 206e219. M. Girault, D. Maillet, F. Bonthoux, B. Galland, P. Martin, R. Braconnier, J.-R. Fontaine, Estimation of time-varying pollutant emission rates in a ventilated enclosure: inversion of a reduced model obtained by experimental application of the modal identification method (22pp). Inverse Problems 24 (1) (2008) 015021. Fluent. http://www.fluent.com/. M. Girault, D. Petit, Resolution of linear inverse forced convection problems using model reduction by the modal identification method: application to

1368

[42]

[43]

[44]

[45]

[46]

Y. Rouizi et al. / International Journal of Thermal Sciences 49 (2010) 1354e1368 turbulent flow in parallel-plate duct. International Journal of Heat and Mass Transfer 47 (17e18) (2004) 3909e3925. B.F. Armaly, F. Durst, J.C.F. Pereira, B. Schonung, Experimental and theoretical investigation of backward-facing step flow. Journal of Fluid Mechanics 127 (1983) 473e496. I.E. Barton, The entrance effect of laminar flow over a backward-facing step geometry. International Journal of Numerical Methods in Fluids 25 (1997) 633e644. E. Erturk, Numerical solutions of 2D steady incompressible flow over a backward-facing step, part i: high Reynolds number solutions. Computers & Fluids 37 (6) (2008) 633e655. P.M. Gresho, D.K. Gartling, J.R. Torczynski, K.A. Cliffe, K.H. Winters, T.J. Garrat, A. Spence, J.W. Goodrich, Is the steady viscous incompressible two-dimensional flow over a backward facing step at Re ¼ 800 stable? International Journal of Numerical Methods in Fluids 17 (1993) 501e541. D.K. Gartling, A test problem for outflow boundary conditions e flow over a backward-facing step. International Journal of Numerical Methods in Fluids 11 (1990) 953e967.

[47] D. Barkley, M.G.M. Gomes, R.D. Henderson, Three-dimensional instability in flow over a backward-facing step. Journal of Fluid Mechanics 473 (2002) 167e190. [48] Y. Zang, R.L. Street, J.R. Koseff, A non-staggered grid, fractional step method for time-dependent incompressible NaviereStokes equations in curvilinear coordinates. Journal of Computational Physics 114 (1994) 18e33. [49] K. El Khoury, A. Neveu, Analyse modale des systèmes thermiques en présence de transferts non-réciproques. International Journal of Heat and Mass Transfer 32 (1989) 213e226. [50] P.E. Gill, W. Murray, M.H. Wright, Practical Optimization. Academic Press, London, 1992. [51] M. Clerc, L’optimisation par essaims particulaires. Hermes-Lavoisier, 2005. [52] M. Clerc, M. Kennedy, The particle swarm e explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation 6 (1) (2002) 58e73. [53] M. Clerc, The swarm and the queen: towards a deterministic and adaptive particle swarm optimization, in: Proceedings ICEC, vol. 3, 1999, pp. 1951e1957. [54] R.C. Eberhart, Y. Shi, Comparing inertia weights and constriction factors in particle swarm optimization. vol. 1, 2000, pp. 84e88.