ASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 Published online 13 July 2010 in Wiley Online Library (wileyonlinelibrary.com) DOI:10.1002/apj.481 Special Theme Research Article Proper orthogonal decomposition and its applications Sanjeev Sanghi1 * and Nadeem Hasan2 1 2 Department of Applied Mechanics, Indian Institute of Technology Delhi, New Delhi, India Department of Mechanical Engineering, Aligarh Muslim University, Aligarh, India Received 18 August 2009; Revised 21 May 2010; Accepted 31 May 2010 ABSTRACT: The proper orthogonal decomposition (POD) has become a very useful tool in the analysis and lowdimensional modelling of flows. It provides an objective way of identifying the ‘coherent’ structures in a turbulent flow. The application of POD to the case of a thermally driven two-dimensional flow of air in a horizontal rotating cylinder is presented. The data for the POD analysis are obtained by numerical integrations of the governing equations of mass, momentum and energy. The decomposition based on POD modes or eigenfunctions is shown to converge to within 5% deviation of the computational data for a maximum of 15 modes for the different cases. The presence of degenerate eigenvalues is an indicator of travelling waves in the flow, and this is confirmed by symmetry in both space and time for the corresponding eigenfunctions. Wave speeds are also determined for these travelling waves. Furthermore, low-dimensional models are constructed employing a Galerkin procedure. The low-dimensional models yield accurate qualitative as well as quantitative behaviour of the system. Not more than 20 modes are required in the low-dimensional models to accurately model the system dynamics. The ability of low-dimensional models to accurately predict the system behaviour for the set of parameters different from the one they were constructed from is also examined. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. KEYWORDS: natural convection; rotational buoyancy; proper orthogonal decomposition INTRODUCTION The proper orthogonal or Karhunen–Loève decomposition is an established technique for extracting the most ‘energetic’ or dominant spatial structures in the average sense from an ensemble of functions representing an infinite-dimensional process. It provides a spatial basis of a vector space onto which the average square projection of the ensemble is maximum. The ensemble may comprise data over a series of experiments or data obtained at different time instants during a single experiment. Once the spatial structure of the ensemble is captured in the spatial basis functions or modes, the members of the ensemble can be represented as linear superpositions of these modes. The advantage of using the proper orthogonal or K-L basis is that finite truncations of the modal expansion of the ensemble members yield the least error in the mean square sense compared with the other possible orthogonal bases having the same dimension. In the context of fluid mechanics, this means that the structure of a complex spatiotemporal field such as the velocity field, representing a *Correspondence to: Sanjeev Sanghi, Department of Applied Mechanics, Indian Institute of Technology Delhi, New Delhi, India. E-mail: sanghi@am.iitd.ac.in 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Curtin University is a trademark of Curtin University of Technology flow, could be captured in an average sense by a relatively few number of proper orthogonal decomposition (POD) modes. Thus, a data ensemble comprising the information of the velocity field at different instants of time could be represented by linear superposition of few POD modes weighted with time varying coefficients. If the partial differential equations, such as the Navier–Stokes equations, governing the fluid flow are projected onto these spatial POD modes, a system of ordinary differential equations for the temporal coefficients is obtained. Hence, the time evolution of the system could be studied in a state space involving few degrees of freedom making it easier to investigate the underlying physics. The earliest attempts to apply the above ideas in studying fluid motions came in the context of identifying and studying the role and dynamics of large-scale average spatial structures of the fluctuating velocity field of a turbulent flow. These large-scale structures are the dominant ones in terms of the average energy of the fluctuating flow and are commonly referred to as coherent structures. However, due to lack of a precise definition of a coherent structure, there is a considerable amount of subjectivity in its identification in a turbulent flow. The POD provides an objective way of identifying such structures. Asia-Pacific Journal of Chemical Engineering The ideas involving extraction of dominant spatial structures using POD and the related low-dimensional modelling have been mostly applied to turbulent flows. Low-dimensional models based on POD modes have been used to study the dynamics of wall turbulence in a turbulent boundary layer. These works were based on experimental data.[1,2] Moin and Moser[3] used direct numerical simulation (DNS) data for studying a turbulent channel flow. They employed Fourier modes in the streamwise and spanwise directions and POD modes in the wall normal direction. Sirovich et al .[4] also employed the DNS data for channel flow and found the intermittency associated with the burst-sweep process as observed in experiments. Rempfer and Fasel[5] employed the data generated from DNS calculations of transition of the laminar boundary layers. They found that the empirical eigenvalues and eigenfunctions almost occur in pairs, implying that there is an approximate translational invariance of relatively slowly growing structures. In an attempt to understand the fundamental mechanisms responsible for self-sustenance of turbulence, most of the studies involving POD and low-dimensional modelling have been performed in the domain of turbulent flows. There have been very few studies involving the application of POD to unsteady laminar flows where the unsteadiness is either due to inherent instabilities or due to some unsteady forcing. More recently, however, applications of POD to problems such as the driven cavity flows, flow past bluff bodies and flows with heat transfer have appeared in the literature. Cazemier et al .[6] applied POD and low-dimensional modelling to investigate the flow in a square lid-driven cavity. They obtained the first 80 POD modes from a dataset comprising of 700 snapshots taken from DNS calculations at Re = 22,000. Jørgensen[7] applied the ideas of POD and low-dimensional modelling to study the axisymmetric lid-driven flow in a cylindrical cavity with a rotating rod. The problem was investigated by fixing the aspect ratio of the cylindrical cavity and the ratio of rod radius to the cylinder radius. The bifurcations exhibited by the full DNS calculations were captured faithfully by the POD-Galerkin models. Ma and Karniadakis[8] investigated the stability and dynamics of three-dimensional periodic states in flow past a circular cylinder using low-dimensional modelling. They have shown that the limit cycle is reproduced very accurately with only 20 three-dimensional modes. Galletti et al . have considered a POD-Galerkin model for two-dimensional vortex shedding past a confined square cylinder.[9] They investigated the validity of such a model for Reynolds numbers and blockage ratios that are different from those from which the lowdimensional model was derived. Their work demonstrates that reliable results can be obtained over a short interval of time. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. PROPER ORTHOGONAL DECOMPOSITION There are some studies involving the application of POD and POD-Galerkin modelling of flows with heat transfer. Turbulent Rayleigh-Bénard convection problems were considered by Sirovich and co-workers.[10 – 12] The change in the structure of the flowfield with increase in Ra was studied via empirical eigenfunctions obtained by application of the Karhunen–Loève procedure to the numerical data. Sahan et al .[13] identified the organized spatiotemporal structures for transitional flow with heat transfer in a periodically grooved channel. In this study, the buoyancy effects were neglected and the eigenfunctions for velocity and thermal fields were calculated separately. Low-dimensional models based on four modes yielded accurate results. Podvin and Le Quéré[14] investigated the two-dimensional buoyancy-driven flow of air in a differentially heated tall cavity. They employed a coupled type of decomposition, where composite vector eigenfunctions having components representing the velocity and thermal fields were obtained. In their low-dimensional models, the effect of neglected modes is modelled through a Heisenberg type of model. For a slightly supercritical Rayleigh number, they have shown that a two-mode model captures the dynamics reasonably accurately. At some distance from the critical point, the flow becomes chaotic and a ten-dimensional model successfully captures the dynamics. Thus, the POD analysis and the associated low-dimensional modelling is an effective reduction tool not only for turbulent flows but also for transitional laminar flows. Hasan and Sanghi[15] considered a POD analysis and low-dimensional modelling of thermally driven two-dimensional flow of air in a horizontal rotating cylinder, subject to Boussinesq approximation. The problem is unsteady due to the harmonic nature of the gravitational buoyancy force with respect to the rotating observer. The data for the POD analysis are obtained by numerical integrations of the governing equations of mass, momentum and energy. The method of snapshots is applied to the DNS snapshots of the fluctuating flow field. The translational symmetry in both space and time of the pair of modes having degenerate (equal) eigenvalues confirms the presence of travelling waves in the flow for several cases of Ra . The shifts in space and time of the structure of the degenerate modes are utilized in the estimation of wave speeds in a given direction. The governing equations for the fluctuations are derived and low-dimensional models are constructed employing a Galerkin procedure. For each of the five cases of Ra , the low-dimensional models yield accurate qualitative as well as quantitative behaviour of the system. Sufficient modes are included in the low-dimensional models so that the modelling of the unresolved scales of motion is not needed to stabilize their solution. Not more than 20 modes are required in the low-dimensional models to accurately model the system dynamics. The ability of low-dimensional Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj 121 122 S. SANGHI AND N. HASAN Asia-Pacific Journal of Chemical Engineering models to accurately predict the system behaviour for the set of parameters different from the one they were constructed from is also examined for the system under consideration. PROPER ORTHOGONAL DECOMPOSITION: THEORY Proper orthogonal decomposition or Karhunen–Loève expansion is a statistical technique that has been known for quite some time but has gained popularity in the subject of fluid mechanics over the last decade. The technique provides an optimal spatial basis for capturing the spatial structure of the most energetic fluctuations in the flowfield in an average sense. Let Um = U(x, τm ) (m varying from 1 to M ) be M realizations or snapshots of the fluctuating flowfield obtained through experimentation or numerical simulations. In this case, the snapshot Um comprises both the fluctuating velocity components and the fluctuating temperature field with the three values (U1 , U2 , U3 ) ≡ (u , v , θ ) assigned at each point. In order to decompose such an ensemble of snapshots, the K-L basis functions φ are taken to be vector functions with three components at each point (φ1 , φ2 , φ3 ) associated with the two velocity components and the thermal field. The vector space in which the decomposition is sought has an inner product defined as According to Lumley and Poje,[16] the introduction of a scaling factor γ is necessary to balance the velocity and temperature fluctuation energies so that the K-L basis captures the structure of the most energetic fluctuations of both the temperature and the velocity fields in a composite manner. They have shown that the proper value of γ that maximizes the average of the square of the projection of the data onto the basis function φ is given as γ = < (u u + v v ) > dA < θ θ > dA (1) The ‘< >’ denotes the ensemble average and in this scenario is taken to be the time average as a result of time stationarity of the flowfield. The K-L basis function φ is governed by the Fredholm Integral of the first kind given as Equation (2) is an eigenvalue problem with λ as the eigenvalue and φ ≡ (φ1 , φ2 , φ3 ) as the eigenfunction or mode. It is known that the eigenvalue problem expressed via (2) has countably infinite, nonnegative eigenvalues and corresponding eigenfunctions. The eigenfunctions {φ}∞ 1 or modes are orthogonal and form a complete basis of an infinite-dimensional linear space spanned by them. The K-L basis is utilized for carrying out the decomposition of the fluctuating flowfield as U(x, τ ) = Rij (x , x )φj (x )dA = λφi (x ), (i , j ) ∈ {1, 3} (2) 2010 Curtin University of Technology and John Wiley & Sons, Ltd. ∞ a (p) (τ )φ (p) (x), (3) 1 where the superscripts in parentheses represent the POD quantum numbers. Once the eigenfunctions are known, the temporal coefficients in (3) can be readily obtained as (4) a (p) (τ ) = (U, φ (p) ) The total energy captured by the expansion in the average sense is given as E = <(U, U) > = < a (p) a (r) > (φ (p) , φ (r) ) = D Rij (x , x ) = <Ui (x , τ )Uj (x , τ ) > p (f1 g1 + f2 g2 + γ f3 g3 )dA (f, g) = The kernel of the integral in (2) is the two-point correlation tensor defined as r λ(p) (5) (1) If the eigenvalues {λ}∞ > 1 are ordered such that λ (2) (3) λ > λ . . ., the sequence on the right of (3) converges more rapidly than with any other basis and only a few modes are needed to obtain a good enough reconstruction of the original data ensemble Ui . The most commonly employed criteria for truncating the infinite sequence in (3) is the retention of the number of modes which capture more than 90% of the average energy in the ensemble Ui with the condition that none of the neglected modes has more than 1% of energy of the most energetic mode. The discrete solution of the eigenvalue problem in (2) involves handling a fairly large algebraic eigenvalue problem in terms of the values of the eigenfunctions at the discrete mesh points. To overcome this difficulty, Sirovich[17] proposed the method of snapshots. In this method, the eigenfunction φ is taken as a linear combination of the snapshots given as φ= M αm U m (6) 1 Utilizing (6) and the definition of the two-point correlation tensor, it can be readily shown that the eigenvalue Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering PROPER ORTHOGONAL DECOMPOSITION problem in (2) can be transformed to an equivalent M × M algebraic eigenvalue problem, Aα = λα (7) 1 (Um , Un ). where Amn = M Since relatively few snapshots are needed to capture the temporal structure of non-turbulent flows, the method of snapshots is preferred in most cases. This method yields at most M modes. GOVERNING EQUATIONS This work applies POD and POD-Galerkin models to study the problem of two-dimensional thermally driven flow in a horizontal steadily rotating cylinder. Similar application in the near wall turbulent boundary layer was first presented by Aubry et al .[1] The problem has been investigated numerically by Hasan and Sanghi,[18] hereafter referred to as HS. A Cartesian frame of reference attached to the rotating cylinder has been employed in the study. To an observer attached to the cylinder, the gravity vector rotates resulting in a time-periodic gravitational buoyancy (GB) force driving the flow. Thus, the problem involves unsteady forcing and is therefore inherently unsteady in nature. It is assumed that the cylinder has been rotating steadily for a sufficient length of time so that the fluid is in a state of solid-body rotation under isothermal conditions. For a figure of the problem, please refer to HS. A temperature perturbation in the form of a periodic cosine function of the angular location of the spatial point is imposed on the wall of the cylinder. The fluid motion relative to the rotating cylinder is induced by the combined action of GB and centrifugal or rotational buoyancy (RB). The dimensionless governing equations of mass, momentum and energy subject to Boussinesq approximation in a Cartesian frame rotating with the cylinder are given as Mass : ∇ · V = 0 (8) dV = −∇pm − Rag Prθ n̂g − Ra Prθ r Momentum : dτ − 2Ta1/2 Pr(k̂ x V ) + Pr∇ 2 V (9) Energy : dθ = ∇ 2θ dτ (10) The vector n̂g in (9) is the unit vector indicating the instantaneous direction of the rotating gravity vector and is defined as n̂g = −î sin(Ta1/2 Pr τ ) − ĵ cos(Ta1/2 Pr τ ) The dimensionless system (8)–(10) shows four dimensionless parameters on which the solution depends. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. These dimensionless parameters are[1] gravitational Rayleigh number (Rag ),[2] rotational Rayleigh number (Ra ),[3] Taylor number (Ta) and[4] Prandtl number (Pr). These are defined as g βTR 3 2 βTR 4 , Ra = , νκ νκ 2 R 4 ν Ta = , Pr = 2 κ ν Rag = (11) The computational study of HS has correlated the spatiotemporal dynamics and the heat transfer characteristics to the changes in Ra and Ta over a wide range at a fixed Rag = 105 and Pr = 0.71. The computational study revealed that for Ra < 105 , the spatial structure of the flow is quite sensitive to the changing orientation of the gravity vector in the rotating frame. In the time domain, for Ra < 105 , the flow exhibits a periodic behaviour with frequency locked onto the rotation frequency of the gravity vector. It was concluded that the temporal variations in the flow for Ra < 105 are caused by harmonic forcing of the GB force. For Ra > 105 , the large-scale spatial structure of the flow is quite insensitive to the harmonic time varying GB force. Examination of the temporal structure via fast Fourier transform revealed bifurcations from periodic to quasi-periodic states for Ra ∈ (106 , 107 ). The study on variation of time mean heat transfer over the hot portion of the wall of the cylinder with Ra revealed that significant control of flow via rotation could be achieved. The mean heat transfer characteristic is highly non-monotonic, indicating the aiding and mitigating effects of rotation on convection in different ranges of Ra . It was shown that at low rotation rates (Ra < 2 × 103 ), the flow is governed essentially by the GB force with the centrifugal effects being negligible. For Ra > 105 , the role of gravity progressively diminishes and the flow for Ra in the range (5 × 105 , 107 ) is driven essentially by the centrifugal or RB forces. In this range, the large-scale flow structure consists of a twocell convection with the interface of the cells aligned along the horizontal diameter and an increase in Ra brings about an increase in convection. The two cells become increasingly distorted and multiple rolls appear as Ra approaches 107 . In the range given by 2 × 103 ≤ Ra ≤ 5 × 105 , the flow is controlled by both the body forces. Under the competing influence of the two buoyancy forces, it is shown that the convection is in general suppressed with very low flow velocities and associated heat transfer. The fluctuation levels are large at low rotation rates (Ra ∼ 102 ) and decrease progressively as the influence of the rotating gravity vector decreases with increase in Ra . The fluctuation levels become extremely small to render the flow to be almost steady at Ra ∼ = 105 . At 5 higher rotation rates (Ra > 10 ), the fluctuation levels Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj 123 124 S. SANGHI AND N. HASAN again start to increase as an outcome of bifurcations to quasi-periodic states. As summarized above, there is a wide variation in the spatiotemporal structure of the flow as Ra is varied. This study is motivated by the fact that for a complex unsteady flow, such as the one under consideration, the POD analysis can provide useful information regarding the flow physics by representing the flow in terms of few modes. Furthermore, it is of interest to examine the feasibility of constructing lowdimensional models and to assess the performance of such models for the class of problem under study. The procedure for carrying out the decomposition and the results of the decomposition are presented in the next section. POD ANALYSIS From this point onwards, unless stated otherwise, it is understood that snapshot refers to the instantaneous spatial distribution for the fluctuating flow field. The fast convergence of POD modes is observed as a small number of modes are needed to capture almost 99% energy of the snapshots in the average sense. The distribution of the spectrum is shown in Fig. 1a–e where the respective eigenvalues are plotted. The spectrum of eigenvalues in each of these cases decays rapidly. At Ra = 102 , seven modes capture more than 99% of the total energy. For Ra = 103 , eight modes are needed to capture more than 99% of the total energy. For Ra = 104 and Ra = 105 , only two modes are sufficient. For Ra = 106 , only six modes are needed. From a qualitative viewpoint, it is shown that the spatial structure of the instantaneous flow field as obtained by DNS calculations can be faithfully captured using few modes. For the case of Ra = 102 , ten modes appear to capture the instantaneous flow structure very well. Another observation from the POD reconstructions is that while 10–15 modes capture the relevant scales for the case of Ra = 102 and 103 , only two modes are needed for the cases of Ra = 104 and 105 . This can be correlated to the temporal structure of the flow field for these cases, investigated via Fast Fourier Transform (FFT) of the time histories, in the earlier work of the authors. As shown in the computational study of HS, the flows for Ra = 102 –105 are periodic, and the frequency spectrum exhibits a fundamental frequency equal to the rotation frequency of the gravity vector with high-frequency low-power harmonics. These harmonics represent less energetic small scales in the flow. As Ra is increased from 102 to 105 , the effects of gravity are reduced and the centrifugal buoyancy starts to affect the flow. The number of harmonics progressively decreases, indicating the reduction in the scales of the fluctuating motion. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pacific Journal of Chemical Engineering The POD reconstructions also reflect the phenomenon of reduction in scales as Ra is increased from 102 to 105 , with lesser number of modes needed to capture the structure of the fluctuating flow in the average sense. It is observed that for Ra = 103 and 106 , a few number of eigenvalues are degenerate, i.e. occur in pairs. The former case is discussed in detail in this work. The eigenvalues for modes[3,4] ,[7,8] and[9,10] are some of the more energetic degenerate pairs. The degeneracy of the eigenspectrum has its roots in the symmetries of the flow solution as demonstrated in the works of Aubry and co-workers.[1,19,20] To observe such symmetries, time histories of the temporal coefficients of a degenerate pair of modes (a (q) (t), a (q+1) (t)) and the contributions of the individual modes of the selected degenerate pair to the flow field at an instant are compared. It is observed that the contour patterns are shifted in the circumferential direction. This clearly suggests a translational shift in the circumferential direction in the spatial structure of the degenerate pair of modes. Some of the structures experience some distortion along with a shift in the circumferential direction. The distortion is expected due to nonlinear convective interactions and diffusive effects. It is also seen that the pair of coefficients a (3) and a (4) as well as other pairs are almost identical except for a phase shift in time. It is observed that the spatial structures of φ3(3) and φ3(4) appear to exhibit a translational symmetry in the circumferential direction. Furthermore, while mode 3 is shifted in the counter clockwise direction relative to mode 4 in space, mode 3 also leads mode 4 in time. This clearly suggests that the degenerate pair[3,4] represents a propagating structure or a travelling wave in the clockwise direction with the flow. Rempfer and Fasel[5] also exploited the space-time symmetry of a propagating flow structure to demonstrate the existence of travelling waves in a transitional boundary-layer flow over a flat plate. Similar space-time symmetries are exhibited by the other degenerate pairs. The deviations from the perfect translational symmetry in space and time can be attributed to the fact that the propagating structures represent modulated travelling waves. Thus, it can be argued that at Ra = 103 , the flow exhibits travelling waves propagating in the circumferential direction. In the numerical study of HS, it was argued on the basis of spatial flow structure that at Ra = 103 , the centrifugal buoyancy force is too weak to exert any significant influence on the flow. The flow is essentially controlled by the periodic GB force. Thus, the travelling waves found are generated due to the combined action of gravity, the fluid inertia and the viscous forces. Such waves have not been reported in the earlier studies on rotating rectangular containers.[21,22] This also demonstrates the power of POD analysis in detecting organized structures from a computational database. Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering PROPER ORTHOGONAL DECOMPOSITION (b) 104 104 103 103 102 102 101 101 λ λ (a) 100 100 10-1 10-1 10-2 10-2 10 20 10 30 20 (c) (d) 104 103 103 102 102 101 101 100 100 10-1 λ λ 30 p p 10-1 10-2 10-2 10-3 10-4 10-3 10-5 10-4 5 10 p (e) 15 20 10-6 5 10 p 15 20 104 103 102 101 λ 100 10-1 10-2 10-3 10-4 10 20 30 p Figure 1. pth eigenvalue, λ, as a function of p at (a) Ra = 102 , (b) Ra = 103 , (c) Ra = 104 , (d) Ra = 105 and (e) Ra = 106 . 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj 125 126 S. SANGHI AND N. HASAN Asia-Pacific Journal of Chemical Engineering LOW-DIMENSIONAL MODELS The low-dimensional model refers to the system of Ordinary Differential Equation (ODE) governing the evolution of temporal coefficients in (9). The standard procedure for obtaining such system of ODEs is Galerkin projection. In this procedure, the governing equations for the flow field are projected onto the individual POD modes or basis functions. The Galerkin projection requires the basis functions to fulfil certain conditions. In the context of a Galerkin procedure for an incompressible flow, the basis functions should be capable of representing all solenoidal velocity fields that satisfy the boundary conditions. In this regard, the POD eigenfunctions provide an excellent choice of basis functions. As a consequence of the fact that the POD eigenfunctions can be represented as linear combinations of the flowfield snapshots, any property or boundary condition of the flow, which is expressed via linear homogenous equations, is passed on to these individual basis functions. The incompressibility constraint, no-slip boundary conditions and periodic boundary conditions are some of the examples of such properties. In this study, since the POD eigenfunctions have been obtained for the fluctuating component of the flow field, the homogeneous boundary conditions for the fluctuating or unsteady flow (u = v = θ = 0) are passed onto each eigenfunction making it vanish at the domain boundary. Since each individual eigenfunction is capable of satisfying the boundary condition, no compatibility constraints are imposed on these basis functions. The Galerkin projection leads to an infinite system of ODEs for the temporal modal coefficients of the form da (p) = F (p) (a (1) , a (2) , . . . , τ ), dτ p = 1, 2, . . . , ∞ (12) The infinite system in (12) is truncated at some level of the quantum number p to yield the low-dimensional model with the hope that the effect of truncation on the system dynamics is not significant. More would be said regarding the truncation of the infinite system in (12) when the validation and performance assessment of the low-dimensional model is presented. In order to carry out the Galerkin projection for the present problem, the governing equations for the fluctuating or unsteady flow have to be obtained from the full set of Eqns (8)–(10). Following the Reynolds decomposition approach, the instantaneous flow variables in Eqns (8)–(10) are considered to comprise a steady or mean part and a fluctuating component. Symbolically, this is expressed as vi (x , y, t) = Umi (x , y) + ui (x , y, t) θ (x , y, t) = m (x , y) + θ (x , y, t) p(x , y, t) = Pm (x , y) + p (x , y, t) (13) 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Substituting in Eqns (8)–(10) and time averaging, the terms lead to equations governing the mean flow. The mean flow equations are then subtracted from Eqns (8)–(10) to yield the equations for the fluctuations. The governing equations for the fluctuating flowfield in Cartesian tensor notation are ∂ uj ∂ xj =0 (14) 1/2 { (m + θ ) sin(Ta Pr τ ) >}δi 1 1/2 ∂ Pr τ ) = Rag Pr − < θ sin(Ta 1/2 ∂ τ +{ (m + θ ) cos(Ta Pr τ ) >}δi 2 −<θ cos(Ta1/2 Pr τ ) ∂p ∂Umi − uj − 2Ta1/2 Pr ei 3j uj − Ra Prθ xi − ∂xi ∂xj ui ∂ < ui uj > ∂ 2 ui ∂ui ∂ui − Umj − uj + + Pr ∂xj ∂xj ∂xj ∂xj ∂xj (15) θ uj > ∂< ∂m ∂θ ∂θ ∂θ = −uj − Umj − uj + ∂τ ∂xj ∂xj ∂xj ∂xj + ∂ 2θ ∂xj ∂xj (16) In equations (14)–(16), (i , j ) ∈ (1, 2). The POD expansions for the fluctuating variables (u1 , u2 , θ ) are substituted in (15) and (16), and then the Galerkin projection is carried out to obtain the system of ODEs involving the temporal coefficients. For notational convenience, the mean and the fluctuating thermal field are denoted as Um3 and u3 . The various correlations, i.e. <ui uj > needed to close the system of equations have been expressed in terms of the POD modes in the following manner: < a (r) a (s) > φi(r) φj(s) <ui uj >= = λ(r) φi(r) φj(r) , i = 1, . . . , 3, j = 1, 2 r (17) The Galerkin projection yields the following system of equations: da (p) = Lpr (τ )a (r) − Qprs a (r) a (s) dt + Fp (τ ), p, r, s = 1, 2, . . . , ∞ (18) Truncating the above system by retaining the first N modes in the system in (18) yields an N -dimensional model.[15] The term Fp (τ ) in (18) is an harmonic forcing function associated with the rotating gravity vector. Thus, the system in (18) has a character of a nonlinear forced oscillator. Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering In order to validate and assess the performance of the low-dimensional models, the five cases of Ra which have been analyzed using POD are considered. Comparison between the temporal coefficients obtained through projection of the POD modes onto the computational data and the solutions of the POD-Galerkin models is taken to be the criterion for valid and accurate solutions of these low-dimensional models. The truncation of the POD-Galerkin models for accurate capturing of the system dynamics is guided by the convergence characteristics of the POD mode based reconstructions. Another issue related to the truncation of the PODGalerkin models is the effect of neglected modes on the long-term behaviour of the low-dimensional model. It is known that if sufficient number of modes are not included then the required amount of energy of the system is not dissipated and accumulates in the first few modes or the large scales of the motion. Thus, in the large time limit, the system becomes unstable. The remedy is to include appropriate number of modes or model the effect of neglected modes. The modelling approach introduces additional parameters in the model which have to be suitably calibrated. In contrast, a model-free approach is employed in this study where enough number of modes has been utilized not only to stabilize the low-dimensional system but to accurately reproduce the system dynamics. For Ra = 102 , a five-dimensional model is found to be sufficient to carry out integrations of the lowdimensional model for very long times (∼400 dimensionless units) without any blowing up of the solution. Thus, sufficient amount of dissipation is achieved even for a five-mode low-dimensional model. The lowdimensional models for the different cases are constructed from a few modes greater than the modes needed to achieve reconstructions of the flow field within 5% of the DNS values. The temporal evolution of the modal coefficients from the low-dimensional models is compared with that of the modal coefficients obtained via projection of the eigenfunctions onto the DNS snaps. The initial conditions for the system of first-order ODEs is taken to be the values of the temporal coefficients obtained by projecting the corresponding eigenfunctions onto the first DNS snap. At this stage, it is worth discussing an important aspect of the definition of the inner product defined for the purpose of a combined decomposition of the velocity and the thermal fields. The inner product involves the scaling factor γ . The choice of the scaling factor that maximizes the projection of the data ensemble onto a fixed set of modes is given as in (1). In the work of Podvin and Le Quéré,[14] it has been mentioned that an ad hoc value of unity for the scaling factor also works reasonably well and that the value based on the criterion given in (1) only serves to provide slightly better performance of the low-dimensional model. However, 2010 Curtin University of Technology and John Wiley & Sons, Ltd. PROPER ORTHOGONAL DECOMPOSITION when a value of unity was tried in this study, the eigenfunctions failed to capture the dominant structures of the fluctuating velocity and thermal fields accurately. This is not surprising when one considers the values of γ obtained using (1) for the five cases of Ra . For the five cases in this study, the value of γ was found to vary between 104 and 107 . Such values exist primarily due to the choice of the velocity and the temperature scales in this study. In the study by Podvin and Le Quéré[14] , the scalings employed were such that both the velocity and temperature are O[1] and therefore the value of γ was reported to be 4.68. Evidently, a value of unity being of comparable order still yielded reasonable results from their low-dimensional models. However, in this case low-dimensional models provide proper results only when the appropriate value of γ is used. At Ra = 102 , the low-dimensional model obtained by retaining the first 20 modes in (18) captures the temporal evolution of the modal coefficients quite faithfully. In fact, the time evolution of the modal coefficients as predicted by the model and the corresponding DNS values for the first five modes is virtually indistinguishable. For the 10th and 15th modes, the agreement is still quite good in the qualitative sense. The time period of the various temporal coefficients is found to be 0.088 dimensionless units. The period of rotation of the gravity vector at this Ra is 0.0885 dimensionless units. Therefore, the phenomenon of frequency locking is captured quite well. Another feature that can be observed in the solution of the low-dimensional model is the fact that while the more energetic structures have a simple periodic structure in time, the less energetic structures (10th and 15th modes) have a greater amount of complexity in their temporal structure. This supports the proposition made earlier that the low-power high-frequency harmonics revealed by the FFT in the temporal structure of the flow are reflections of the dynamic evolution of the less energetic structures or smaller scales of flow. As shown in the study of HS, the phenomenon of frequency locking persists for Ra ≤ 105 . This is reflected in the solutions of the low-dimensional models for Ra = 103 , 104 and 105 . The time periods of the various temporal coefficients for Ra = 103 , 104 and 105 are found to be 0.0278, 0.00886 and 0.0028 dimensionless units, respectively. These agree quite well with the corresponding periods of the rotating gravity vector given as 0.02798, 0.00885 and 0.002798, respectively. As mentioned earlier, for Ra ≤ 105 , the fluctuation levels of the flow decrease with increase in Ra . This feature is also readily captured by the low-dimensional models. The decreasing complexity in the flow structure, as confirmed by the POD analysis, is also confirmed by the fact that the size of the low-dimensional model needed to accurately model the system dynamics progressively decreases from 20 to 4 with increase in Ra upto 105 . For Ra = 106 , the flow is quasi-periodic and the level Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj 127 128 S. SANGHI AND N. HASAN of fluctuation again begins to rise. The trajectories of the temporal coefficients in the time domain obtained from a ten-dimensional model reproduces the dynamics quite nicely as the trajectories follow their DNS paths. CONCLUSIONS This study demonstrates the effectiveness of POD as a tool for the analysis of a complex fluid flow phenomenon involving unsteady flow with heat transfer. It has been demonstrated that relatively few modes are needed to capture the structure of the unsteady flowfield within a specified level of accuracy. The presence of travelling waves in the flow has been detected by employing POD decomposition to the computational data at Ra = 103 . At Ra = 103 , the waves propagate predominantly in the circumferential direction in the clockwise manner and are believed to be caused by the interaction between the GB, the fluid inertia and the viscous forces as the centrifugal or RB is shown to have a negligible effect on the flow at this value of Ra . These travelling waves are actually found to exist for values of Ra in the range 103 –1.9 × 103 . The space-time symmetries of the degenerate pairs of modes is utilized to estimate the average circumferential wave speeds for these cases. The low-dimensional models based on POD modes have been constructed and issues involved in their closure arising out of the specific nature of the problem under consideration have been successfully resolved. An important aspect of POD analysis that has not received due importance for non-isothermal flows is the definition of the inner product space in which the combined decomposition of velocity and temperature field is desired. It has been shown that for situations where the order of magnitude of velocity and the temperature fields are not comparable, the choice of 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Asia-Pacific Journal of Chemical Engineering the proper value of the scaling factor γ is critical to the success of the low-dimensional model. Having included sufficient number of modes in the low-dimensional models so that sufficient amount of dissipation is present to stabilize their solution for several hundred rotation cycles of the gravity vector, not more than 20 mode models were required to accurately capture the temporal dynamics of the system for the five cases under consideration. REFERENCES [1] N. Aubry, P. Holmes, J.L. Lumley, E. Stone. J. Fluid Mech., 1988; 192, 115–173. [2] S. Sanghi, N. Aubry. J. Fluid Mech., 1993; 247, 455–488. [3] P. Moin, R.D. Moser. J. Fluid Mech., 1989; 200, 471–509. [4] L. Sirovich, K.S. Balle, L.R. Keefe. Phys. Fluid A, 1990; 2, 2217–2226. [5] D. Rempfer, H. Fasel. J. Fluid Mech., 1994; 260, 351–375. [6] W. Cazemier, R.W.C.P. Verstappen, A.E.P. Veldman. Phys. Fluid A, 1998; 10, 1685–1699. [7] B.H. Jørgensen, PhD thesis, Department of Energy Engineering, Technical University of Denmark, 2000. [8] X. Ma, G.E. Karniadakis. J. Fluid Mech., 2002; 458, 181–190. [9] B. Galletti, C.H. Bruneau, L. Zannetti, A. Iollo. J. Fluid Mech., 2004; 503, 161–170. [10] L. Sirovich, H. Park. Phys. Fluid A, 1990; 2, 1649–1658. [11] L. Sirovich, A.E. Deane. J. Fluid Mech., 1991; 222, 251–265. [12] A.E. Deane, L. Sirovich. J. Fluid Mech., 1991; 222, 231–250. [13] R.A. Sahan, A. Liakopoulos, H. Gunes. Phys. Fluid, 1997; 9, 551–565. [14] B. Podvin, P.L. Le Quéré. Phys. Fluid, 2001; 13, 3204–3214. [15] N. Hasan, S. Sanghi. J. Fluid Mech., 2007; 573, 265–295. [16] J.L. Lumley, A. Poje. Phys. Fluid, 1997; 9, 2023–2031. [17] L. Sirovich. Quart. Appl. Mech., 1987; 45, 561. [18] N. Hasan, S. Sanghi. J. Heat Transfer Trans. ASME, 2004; 126, 963–984. [19] N. Aubry, R. Guyonnet, R. Lima. J. Nonlinear Sci., 1992; 2, 183–215. [20] N. Aubry, W.Y. Lian. Lectures Appl. Math., 1993; 29, 71–85. [21] F.J. Hamady, J.R. Lloyd, K.T. Yang, H.Q. Yang. J. Heat Transfer, 1994; 116, 136–143. [22] Y.T. Ker, T.F. Lin. Int. J. Heat Mass Transfer, 1996; 39, 3193–3210. Asia-Pac. J. Chem. Eng. 2011; 6: 120–128 DOI: 10.1002/apj

1/--страниц