Real-time optimization using proper orthogonal decomposition: Free surface shape prediction due to underwater bubble dynamics

Computers & Fluids 36 (2007) 499–512 www.elsevier.com/locate/compfluid

Real-time optimization using proper orthogonal decomposition: Free surface shape prediction due to underwater bubble dynamics D. My-Ha a, K.M. Lim a

a,b

, B.C. Khoo

a,b,*

, K. Willcox

a,c

High Performance Computation for Engineered Systems, Singapore–MIT Alliance, 4 Engineering Drive 3, Singapore 117576, Singapore b Department of Mechanical Engineering, Faculty of Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore c Aerospace Computational Design Laboratory, Massachusetts Institute of Technology, Cambridge, MA 02139, United States Received 8 December 2004; received in revised form 2 August 2005; accepted 16 January 2006 Available online 7 November 2006

Abstract A new data-driven reduced-order modeling approach for real-time optimization applications is presented. The proper orthogonal decomposition (POD) technique is used for the reduced-order model, with the basis functions determined from an ensemble of offline high-fidelity simulations. For optimization in real time, a rapid two-stage approach is taken in the online phase: the POD coefficients are first determined by solving a small optimization problem, and the desired parameters are subsequently obtained by interpolation of the POD coefficients using precomputed information from the simulation ensemble. This method is applied to optimizing parameters for underwater bubble explosions so that a desired free surface shape is generated. For time-critical applications, such as using the water barrier generated to stop sea-skimming objects, the time available for online optimization is limited to about 30 s. Results for two-dimensional simulation, using a personal computer (dual CPU running at 2.8 GHz), show that our new methodology is able to meet such a critical time requirement. For three-dimensional simulations, the time taken for computation increases, and a faster computer or parallel implementation is required.  2006 Elsevier Ltd. All rights reserved.

1. Introduction Model reduction is a powerful tool that allows the systematic generation of cost-efficient representations of large-scale computational systems, such as those resulting from discretization of partial differential equations (PDEs). Reduced-order models are essential for achieving real-time control and optimization of dynamic processes. In many cases, models for such applications yield very large systems that are computationally intensive to solve. A necessary

*

Corresponding author. Address: Department of Mechanical Engineering, Faculty of Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore. Tel.: +65 68742889; fax: +65 67791459. E-mail addresses: [email protected] (D. My-Ha), limkm@nus. edu.sg (K.M. Lim), [email protected] (B.C. Khoo), [email protected] (K. Willcox). 0045-7930/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2006.01.016

element in achieving a real-time simulation capability is the development of accurate, efficient computational models that can be solved sufficiently rapidly to permit decision-making in real time. This paper proposes a datadriven, low-order modeling approach that targets critical real-time optimization and inverse problem applications. Most reduction methods for large-scale systems are based on a projection framework, that is, the large-scale governing equations are projected onto the subspace spanned by a reduced-space basis, yielding a low-order dynamical system. The projection basis can be determined using a number of methods, including Krylov subspace methods [1,2], balanced truncation [3–5], reduced-basis methods [6], and the proper orthogonal decomposition (POD) [7–9]. The POD, also known as Karhunen–Loe´ve decomposition [10] and principal components analysis [11], yields an empirical data-driven basis that is based on a set of solutions or snapshots of the high-fidelity computational model. The POD

500

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

basis functions are optimal in that the error between the original and low-order representation of the snapshot data is minimized for a given basis size. Model reduction methodology has been developed and applied in a wide range of applications, such as fluid dynamics [12] and integrated circuits [13], among others, with focus mainly on approximation of the forward simulation model. Model reduction applications in an optimization context are less common, although success has been demonstrated for cases when the number of parameters in the optimization application is small, particularly in optimal control applications [14,15]. In some cases, an explicit low-order dynamical system may not be required for real-time simulation or decisionmaking needs. Instead, dynamic data obtained from a priori simulation of the large-scale model (i.e., in an offline phase) can be combined with a low-order representation and an interpolation method to provide the required information rapidly in an online phase using a data-driven approach. The latter fast online phase is necessary due to real-time constraints, such as a very short time duration for analysis and decision making, present in the problem of interest. For example, the ‘‘gappy’’ POD is a variation of the POD proposed by Everson and Sirovich [16], which allows missing or incomplete data to be reconstructed in a low-dimensional space. This theory has been used recently to propose fast data-driven methods for optimization of steady flows [17] and real-time reconstruction of unsteady flows [18]. In this paper, we use the POD to derive a loworder basis that describes the system output over a range of parametric variations. By combining the low-order output representation with an interpolation scheme, an optimization problem can be solved very rapidly without the need for online simulation of a dynamical system. POD has been used to describe parametric variations for a number of applications. In Epureanu et al. [19], POD was applied to describe turbomachinery flow with a set of snapshots sampled at both varying time and interblade phase angle but at a single Mach number. The flow at varying Mach number that is close to the value used in snapshots is then predicted using the resulting reduced-order model. Ly and Tran [20] used the basic POD and an interpolation scheme for predicting the steady-state temperature distribution of a flow in a square cavity. In this example, snapshots are the steady temperature distributions taken at varying Rayleigh number. Bui-Thanh et al. [17] used POD to predict the steady-state transonic flow passing an airfoil as the Mach number and the angle of attack are varied. For the approach proposed in this paper, the PODbased interpolation scheme is used directly to represent the system output, rather than being applied to the system state. This modification allows us to target the rapid approximate solution of an optimization problem without needing to solve the underlying governing equations in the online phase. In the following section, a general nonlinear constrained optimization problem is described and the proposed new

methodology is presented. The implementation of the method is then described in detail for a particular application—which seeks to solve an inverse problem to determine the positions, depths, and strengths of a set of bubbles so as to attain a specified water plume free surface shape. Numerical simulation of this system is achieved on the large scale using a boundary integral method. Results for the data-driven low-order modeling method are given in both two and three dimensions. Finally, conclusions are drawn. 2. Optimization problem formulation We consider the following problem, which applies to a wide range of systems, including those governed by PDEs: determine the set of inputs contained in the vector w that best matches a target output y*(x, T, w) at a specified final time t = T. x 2 X is defined over the spatial domain X and t 2 [0, T). One approach to solving this problem is to pose it as a constrained optimization problem, where one must solve the underlying PDEs to determine the state and then the corresponding output y*. This approach has been used with considerable success in a number of applications (see, e.g., [21]). For computations sufficiently fast for real-time applications, it is impractical or near impossible to solve the PDEs (which would require, for example, costly CFD solvers). Instead, we propose a data-driven, low-order modeling approach. Our method is of low complexity, but reduced fidelity—it does not satisfy the underlying governing equations. However, as the results will show, the simple approach can be used in certain applications to obtain sufficiently accurate results for rapid response. 2.1. Reduced-order model To formulate the reduced-order model, the output is represented as a linear combination of q basis vectors q X yðx; t; wÞ ¼ aj ðt; wÞ/j ðxÞ; ð1Þ k¼1

where the time-dependent coefficients aj depend on the parameters in w and the basis vectors /j characterize the spatial variation in output. The basis vectors are determined in an offline phase using high-fidelity simulation data as described in the following subsection. Note that the usual approach in a reduced-basis approximation is to use a lowdimensional representation of the system state; however, in (1), the basis vectors are used to represent the output directly. The coefficients aj and parameters w are determined in the online phase using a combination of interpolation and optimization, as described in more detail below. 2.2. Proper orthogonal decomposition The low-order characterization of the output, given by (1), requires a set of basis vectors /j, j = 1, . . ., q. There are several different methods that could be used to define

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

the basis vectors; in this work, we use the POD, also known as singular value decomposition, Karhunen–Loe´ve decomposition [10] and principal components analysis [11]. The POD basis is a desirable choice since, for a given basis size, it provides the optimal least-squares representation of a given data set. In addition, through a wide range of applications, the POD has been shown to provide an excellent low-dimensional characterization of important flow dynamics. However, it is important to note that the decomposition in (1) applies generally to other relevant choices of basis. The POD basis is computed in the offline phase using the method of snapshots [8], which can be described in general terms as follows. Consider a set of m solutions, or snapshots, yi, i = 1, . . ., m, where yi(x) = y(x, ti, wi) is the output solution of the system for a given time ti and given parameter values wi. The POD basis vectors are chosen to minimize the two-norm of the error due to the representation of the snapshots yi using a truncated number of basis functions. The POD basis can be computed efficiently, as noted by Sirovich [8], since the basis vectors are given by a linear combination of the snapshots m X /j ðxÞ ¼ bji y i ðxÞ; j ¼ 1; . . . ; m; ð2Þ i¼1

where the coefficient bji describes the contribution of the ith snapshot to the jth basis vector. It can be shown [8] that the vector of coefficients bj ¼ T j j ½b1 b2    bjm  satisfies Rbj ¼ kj bj ;

ð3Þ

where R 2 Rmm is the correlation matrix, defined as Rij ¼ 1 ðy i ; y j Þ, where (yi, yj) denotes an inner product over the m spatial domain between yi and yj. The number of basis functions retained in the approximation (1), q, is chosen based on the eigenvalues kj, where only those modes corresponding to large kj are selected. Although the low-order modeling method presented here can be applied over the full time horizon t 2 [0, T), in the context of the stated target problem we will focus on a basis to describe only the output at t = T: q X yðx; T ; wÞ ¼ aj ðwÞ/j ðxÞ: ð4Þ k¼1

In this case, the coefficients aj are a function only of the parameters w. 2.3. Parametric interpolation The POD procedure described above applies to parameter-dependent problems for which the snapshots are collected at a series of varying parameters (see, e.g., Epureanu et al. [19], Ly and Tran [20], Bui-Thanh et al. [17]). As proposed in [20], the solution corresponding to arbitrary parameter value can be constructed by interpolation in POD space as follows.

501

Let yi(x) = y(x, T, wi) denote the snapshot corresponding to final time and parameter values wi, i = 1, . . ., m. The m POD method is applied to the set of snapshots fy i gi¼1 to m obtain the orthonormal basis f/i gi¼1 . The coefficient aij describes the amount of POD vector j needed to repPbasis m resent snapshot i, i.e., y i ¼ j¼1 aij /j , and is given by the projection aij ¼ ð/j ; y i Þ. Then, the POD coefficients aj for intermediate values of w that are in the range of interest but not included in the ensemble can be found by interpolation among the aij . By precomputing and storing the values of aij , this interpolation can be performed rapidly in the online phase. The corresponding approximation of y(x, T, w) is then given by (4). 2.4. Reduced-order optimization formulation The reduced-order optimization problem can now be stated in two phases as follows. First, using the representation (4), determine the set of coefficients aj, j = 1, . . ., m that best matches a target output y* at a specified final time T. Second, using a parametric interpolation model, determine the inputs w that best match the desired coefficients aj. A detailed mathematical formulation of this approach is given in the next section in the context of a specific example for bubble dynamics. It will be demonstrated that this formulation leads to a numerical implementation that can be solved very rapidly. 3. Real-time optimization example: water barrier construction A bubble initiated under high pressure will oscillate and eventually collapse unsymmetrically with a water jet directed away from the initial quiescent free surface, resulting in a water plume in the air. Therefore, a set of underwater gas bubbles created by explosions can lead to the formation of a water plume in the sea surface. A water plume that attains a particular shape can serve certain functionality, such as obstructing aerodynamic objects skimming above the water surface. This shape can be achieved by solving an optimization problem which minimizes the deviation of the simulated free surface from some desired target shape. This application is one example of the general problem described in Section 2. In this case, the target output y* is a particular free surface shape, while the decision variables w are bubble parameters such as lateral positions, depths and strengths—decisions that for practical applications, must be made in real time. For typical applications, such as response to an object travelling at speeds in the range slightly less than unity Mach number, this translates to a decision-making time-frame of about 30 s. The 30 s decision time estimate comes from realistic consideration of the time of flight for an aerodynamic object to the target, starting from the instance of detection, but maintaining sufficient time for the building up of the water barrier due to the underwater explosions. One may note that the created water

502

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

barrier at the desired height largely stays in place for a further few seconds which ensures its effectiveness. A number of authors, including Wang et al. [22,23], Zhang et al. [24,25], Rungsiyaphornrat et al. [26], and Best and Kucera [27], have used the boundary integral method (BIM) for simulation of bubble dynamics and with extensive comparisons to the limited experiments available in the literature. BIM has the distinct advantage of reducing the dimension of problem by one, leading to potential savings of computational time and resources. Nevertheless, it is still computationally very expensive to perform full simulation of the bubbles and free surface interaction for an optimization problem, where such a simulation must be repeated many times. In particular, for situations in which decisions must be effected in real time, the computational cost associated with high-fidelity simulation is prohibitive. 3.1. Optimization problem formulation The problem of water barrier construction is to find the best initial locations, depths and strengths for a set of bubbles beneath the quiescent free surface such that the resultant free surface attained is close to the desired prescribed shape. The geometric features of the problem are shown in Fig. 1. Our objective is to minimize the sum of squares of the difference between the constructed surface S and the desired surface S0 using N0 bubbles. The optimization problem is written as kS  S 0 k22

min

frk0 ;k ;ck gk¼1;2;...;N

s:t:

0

jrk01  rk02 j P dr k 1 ; k 2 ¼ 1; . . . ; N 0 ; k 1 6¼ k 2 min 6 k 6 max k ¼ 1; . . . ; N 0 ; cmin 6 ck 6 cmax k ¼ 1; . . . ; N 0 : ð5Þ

In the above, rk0 , ck, and k are the initial lateral position, depth, and strength of the kth bubble, respectively. The depths and strengths are limited to a specified range defined by lower limits cmin and min, and by upper limits cmax and max, respectively. In addition, the bubbles are required to maintain a minimum lateral separation of dr. As will be discussed in the next section, appropriate selection of dr ensures that the free surface S can be approximated by a linear superposition of the individual free surfaces, S ¼ PN 0 k k k¼1 S , where S is the surface corresponding to the kth bubble.

3.2. Simulation model A bubble is initiated as a tiny point of high pressure at a distance h below the initially quiescent free surface. Assuming that the fluid in the domain X is inviscid, incompressible and irrotational, the governing equation for such potential flow is the Laplace equation $2U = 0, where U is the velocity potential. The integral representation is cðxÞUðxÞ ¼

  oUðyÞ oGðx; yÞ Gðx; yÞ  UðyÞ dS: on on oX

Z

ð6Þ

The boundary of the domain oX comprises the bubble surface and free surface. The normal vector n is directed out of the fluid, and o/on = n Æ $. The Green’s function for a three-dimensional domain is G(x, y) = 1/jx  yj, where x is the field point and y is the source point. The term c(x) is the solid angle at point x. The dimensionless kinematic and dynamic boundary conditions governing the motion of the bubble and the free surface are:

Fig. 1. Problem of water barrier construction.

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

dx ¼ rU; dt  k dU 1 V0 2 2 ¼ jrUj þ d ðz  cÞ   þ1 dt 2 V for the bubble surface; dU 1 ¼ jrUj2 þ d2 z for the free surface; dt 2

ð7Þ

ð8Þ ð9Þ

where x = (x, y, z) is the position of a node on the boundary, and the vertical direction z is positive in an upwards direction. The buoyancy, initial inception position, and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi strength of the bubble are d ¼ qgRm =pref , c = h/Rm, and  = p0/pref, respectively. Here, p0 is the initial pressure inside the bubble and the maximum radius Rm of the Rayleigh bubble is obtained in an infinite fluid domain under the uniform pressure p1 and the hydrostatic pressure pref = p1 + qgh, where g is the gravitational acceleration, p1 is the atmospheric pressure, and q is the density of the surrounding fluid. The scaling factors for the velocity pffiffiffiffiffiffiffiffiffiffiffiffi potential and time are U ¼ R p ref m ref =q and t ref ¼ pffiffiffiffiffiffiffiffiffiffiffiffi Rm q=pref , respectively. In (8), V0 and V are the initial and instantaneous volumes of the gas bubble; and k is the ratio of specific heats. The initial radius of the bubble, R0, is the solution of the Rayleigh–Plesset equation [28]  2  3k d2 R 3 dR R0 ¼  1; ð10Þ R 2 þ dt 2 dt R similar to what has been carried out in Rungsiyaphornrat et al. [26]. The initial velocity of the free surface is assumed to be zero, and that of the bubble surface is negligible. Therefore the initial velocity and potential are equal to zero. If the position of the boundary oX and the velocity potential U on it are known, the velocity v = $U = (oU/ on, oU/os) on the nodes along the boundary can be evaluated. The normal velocity oU/on can be obtained directly by solving the boundary integral equation (6). The tangential velocity oU/os is estimated using the interpolation of potential values at the nodes on the boundary. Then, Eqs. (7)–(9) are numerically integrated in time to find the new position and new velocity potential. The forward time integration is carried out using the second-order predictor– corrector scheme with the timestep Dt defined as Dt ¼ minfDt1 ; Dt2 g; Dt1 ¼

ð11Þ DU

 k 2 max j1 þ 12 jvj þ d2 ðz  cÞ þ  VV0 j

for the bubble surface; DU   for the free surface; Dt2 ¼ 1 2  max 2 jvj þ d2 z

ð12Þ ð13Þ

where DU is chosen to be 0.02–0.03. This procedure ensures that the change in the potential at each node on the boundary is bounded by DU. Details of the overall implementation are given in [22,23,25].

503

The simulation is initiated with a small bubble of initial radius R0 corresponding to a strength  at a depth of c below a quiescent free surface. The evolution of the free surface shape is traced, and the shape with the maximum height is captured for subsequent analysis. 3.3. Reduced-order optimization formulation A reduced-order model will be constructed for the output of interest, in this case the free surface shape at the time of maximum of height, using the expansion (4). The complete set of inputs comprises k and ck, the parameters describing the initial strength and depth of the kth bubble, and the lateral positions rk0 . Due to the translational symmetry associated with the input rk0 , the solution due to a bubble at location rk0 is the same as that due to a bubble at location rk0 þ dr but translated by an amount dr. The input rk0 is therefore not included in w, but rather handled explicitly in the basis function representation. A nominal set of basis vectors, /j, is computed for a reference lateral bubble position (taken to be zero); the basis vectors corresponding to an arbitrary lateral bubble position rk0 are denoted /kj  /j ðrk0 Þ and are computed easily by translation of the nominal case. The input vector w therefore comprises the initial strength and depth of each bubble. The numerical boundary integral method described in the previous subsection is used to generate a set of highfidelity snapshots corresponding to the free surface created by a single bubble at varying values of the initial strength  and depth c. After subtracting the mean response from each snapshot, a set of q POD basis vectors is computed. The reduced-order approximation of a free surface corresponding to the kth bubble initiated at a lateral distance rk0 from the origin with a strength k and depth ck is then given by S k ¼ /ðrk0 Þ þ

q X

akj /kj ;

ð14Þ

j¼1

where / is the mean of the snapshots, and /ðrk0 Þ denotes lateral translation of the mean solution by rk0 . The POD coefficients akj are computed by interpolating the values of aij  aj ði ; ci Þ sampled in the snapshot collection process, as described in Section 2.3. 3.4. Solution method To achieve rapid solution of the reduced-order formulation, the problem can be decomposed into two stages as described in Section 2.4, with a slight variation to handle selection of the lateral bubble positions. We first solve for the POD coefficients fakj gk¼1;...;N 0 ;j¼1;...;q and the lateral positions frk0 gk¼1;...;N 0 of the bubbles. Then, using the determined POD coefficients, the strength k and depth ck of the kth bubble are solved for separately in a second stage. It is noted that the lateral positions frk0 g, which should be part of the input parameters, are handled as described in

504

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

Section 3.3, by the basis function representations and not by the POD coefficients. Hence, these quantities are solved at the first stage. The remaining input parameters, k and ck, are obtained at the second stage, after the POD coefficients are determined. The first-stage problem, which seeks to determine the bubble lateral positions and POD coefficients that best match a desired surface shape, can be written as 2

kS  S 0 k2

min

frk0 ;akj gk¼1;2;...;N

0 ;j¼1;2;...;q

where

S¼

N0 X k¼1

k

/ þ

q X j¼1

jrk01  rk02 j P dr amin 6 akj 6 amax j j

s:t:

! akj /kj

k 1 ; k 2 ¼ 1;.. .;N 0 ; k 1 6¼ k 2 k ¼ 1; ...; N 0 ; j ¼ 1;. ..;q: ð15Þ

The lower and upper bounds on the POD coefficients, amin j and amax , are determined by the corresponding bounds on j the bubble strength and depth, specified in (5). k With the optimal POD coefficients fða Þj gk¼1;...;N 0 ;j¼1;...;q obtained from solution of (15), the second-stage problem to determine the input w can now be solved for each individual bubble in turn. Using the interpolation method described in Section 2.3, the POD coefficients are calculated (in the offline phase) at sampled points over the range of strength and depth variations. A piecewise function is used to approximate the variation of the POD coefficient with strength and depth: akj ð; cÞ ¼ fas;t j;k ð; cÞjt 6  6 tþ1 ; cs 6 c 6 csþ1 g s¼1;...;n1 ;

ð16Þ

t¼1;...;m1

where n and m are the numbers of samples in strength and depth ranges of interest, respectively. The function as;t j;k ð; cÞ depends on the particular interpolation method used. In the present work, linear interpolation is used for calculating the POD coefficients; therefore, as;t j;k ð; cÞ is a linear function. We apply the definition (16) to each bubble individually, yielding the second-stage optimization problem for the kth bubble, which seeks to find the optimal set of strengths and depth to best match the POD coefficients found previously: min k ;ck

s:t:

q X

½ða Þkj  akj 2

j¼1

akj ¼ fas;t j;k ð; cÞjt 6  6 tþ1 ; cs 6 c 6 csþ1 g s¼1;...;n1 ; t¼1;...;m1

j ¼ 1; . . . ; q:

ð17Þ

The problem defined by (17) involves minimization of a function with a piecewise linear constraint. The global minimum can be determined by first minimizing the objective function on each segment as;t j;k , and then choosing the solution with minimal objective value.

The solution to the optimization problem is fully specified when the lateral position, depth and strength of each bubble are determined. The evolution time, s (measured from the time of inception to the time of maximum free surface height), of a bubble at an arbitrary depth and strength is estimated by interpolation among the evolution times for the bubbles with depth and strength considered in the snapshots using a cubic spline interpolation. The relative inception time of each bubble is then calculated based on the largest s of the bubbles employed. The bubble with the largest time, smax, will be initiated at time zero; the initiation for bubble k is delayed by an amount smax  sk. This procedure ensures that the bubbles are synchronized at their time of maximum free surface height. In summary, the procedure for solving the problem of water barrier construction can be described by following steps: In the offline phase: (i) Collect a set of m snapshots corresponding to the free surface created by a single bubble at varying values of the initial strength  and depth c; apply POD on the snapshots to obtain the mean / and the orthonormal basis f/j gqi¼j ; compute the POD coefficients aj corresponding to each snapshot. In the online phase: (ii) For a desired free surface shape, solve (15) to determine the lateral position of the bubbles and the POD coefficients that yield a surface shape that best matches the target. (iii) Solve (17) to determine the depths and strengths of the bubbles that best match the desired POD coefficients. (iv) Find the timings and the sequence of the underwater explosions. In Section 4, this procedure is applied to solve the problem of water barrier construction with several different target free surface shapes. 4. Results and discussion 4.1. Bubble dynamics simulations High-fidelity simulations of bubble dynamics are carried out using the boundary integral method and implementation described in Section 3. Fig. 2 shows the bubble and free surface profiles at selected dimensionless times t for an explosion bubble with the initial parameters c = 1.50,  = 100, d = 0.3. After inception, the bubble starts its growth phase. The free surface is pushed up by the growth of the bubble. A dome is formed on the free surface, attaining its maximum height at a time close to that when the bubble reaches its maximum volume, after which the dome falls with the collapse of the bubble.

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512 1

505

1 4 3 2 1

0.5

4 5 6 7

0.5

0

0 4

4 3

–0.5

5

–0.5

2

z –1

6

z –1

1

7

–1.5

–1.5

–2

–2

–2.5

–2.5

–3 –2

–1.5

–1

–0.5

0

r

0.5

1

1.5

2

–3 –2

–1.5

–1

(a)

–0.5

0

r

0.5

1

1.5

2

(b)

Fig. 2. Evolution of bubble and free surface for c = 1.25,  = 100, d = 0.3 during (a) growth phase at time t (1) 0.051, (2) 0.1496, (3) 0.3091, (4) 0.8072 and (b) collapse phase at time t (4) 0.8072, (5) 1.1317, (6) 1.4195, (7) 1.5454.

As the bubble moves closer to the free surface, the interaction between the two becomes stronger. The bubble will vent into the air during expansion if it comes too close to the free surface; however, this situation is not considered in the present work, since we are interested in the water plume that is formed before the collapse. Similarly, the interaction between two adjacent bubbles is highly nonlinear when they are close to each other, and merging of the two bubbles must be considered. Furthermore, in this case numerical instability may occur in the simulation. For the application considered in this paper, merging of the bubbles is undesirable, since the resultant bubble is unlikely to produce a water plume of reasonable dimension; therefore, the two bubbles should be initiated at a sufficient distance so that they will not be too close to each other during their evolution. As was shown in [29], with a minimum lateral separation of dr = 2.0, nonlinear interaction effects between bubbles are small, and the resultant surface can be accurately constructed using linear superposition. 4.2. Approximation with POD An ensemble of 312 snapshots is generated that contains free surface solutions corresponding to 39 values of initial depth in the range c = [1.25, 5.05] with interval step of 0.1, and eight values of strength in the range  = [100, 800] with interval step of 100. Using just the first three POD modes captures 99.998% of the ‘‘energy’’ of this ensemble, defined as the sum of the two-norm of each P snap3

kj

shot vector, i.e., the POD eigenvalues satisfy Pj¼1 ¼ 312 j¼1

kj

0:99998. In Fig. 3, a free surface corresponding to a single bubble is constructed using three POD modes. The depth and strength parameters for this case were not part of the snapshot ensemble; linear interpolation was used to deter-

0.1

0.08

0.06

z 0.04

0.02

0

–0.02 0

5

r

10

15

Fig. 3. Plot of the predicted free surface using three POD modes and the exact surface (the two curves are almost identical) for c = 2.72,  = 405, which was not a member of the snapshot ensemble.

mine the POD coefficients as described in Section 2.3. It can be seen that the constructed surface matches the actual profile closely. Fig. 4 shows the relative error, e, between the constructed surface and the corresponding actual profile for varying numbers of POD modes. Results are given for both a set of parameters that was included in the snapshot ensemble and a set of parameters that was not sampled. The relative error is defined as e¼

kS act  Sk ; kS act k

ð18Þ

where Sact and S are the actual and approximate free surface profiles, respectively. As the figure shows, accurate representations of the free surface can be obtained with just a handful of POD modes.

506

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512 –1

0.12

10

non member member

0.1

–2

10

Relative error

0.08

0.06

z

–3

10

0.04

0.02 –4

10

0

–5

10

1

2

3

4

5

6

7

8

9

–0.02 15

10

10

5

0

5

10

15

5

10

15

r

No. of POD modes used

(a) Fig. 4. Log-scale plot of the relative error between the predicted and the exact free surface for nonmember case (c = 2.72,  = 405) and member case (c = 1.35,  = 400) versus the number of POD mode used.

0.05

0

–0.05

4.3. Implementing the two stages of optimization –0.1

The two stages of optimization are implemented using a commercial software LOQO. The parameters frk0 g in the first optimization stage appear as indices when the POD basis functions /j are represented by a set of discrete nodal positions (obtained from the BIM simulation), and LOQO is not able to use these parameters frk0 g as inputs for the stage 1 optimization. To alleviate the problem, we generate an approximate closed form, given by f ðr  rk0 Þ ¼

–0.15

z –0.2

–0.25

–0.3

–0.35

k 2

C 1 eC2 ðrr0 Þ , for the mean and first POD mode, to be used in the LOQO software. A comparison of the approximate and exact profiles for the two-dimensional case is shown in Fig. 5. Although the two profiles do not agree fully at the peak and around the base, the approximation functions capture the general variation of a high narrow peak at the axis of symmetry and rapid decay with distance from center. For the stage 1 problem, the lateral positions of the bubbles frk0 g are first obtained using LOQO by considering just the mean and first POD mode using the approximate closed forms. This yields good results for the lateral positions of the bubbles, since the position of the peak of the free surface is accurately represented by the approximate mean mode. With this set of lateral bubble positions, LOQO is used again to obtain the POD coefficients for all modes. In the second run of LOQO, the lateral positions are no longer input parameters, so the representation of POD basis functions by the set of discrete nodal positions can be used. Finally, the second stage optimization problem is also solved using the LOQO software. 4.4. Optimization results Results are presented in both two and three dimensions for the reduced-order optimization methodology. For the

–0.4 15

10

5

0

r

(b) Fig. 5. (a) Mean solution and (b) the first POD mode of the ensemble corresponding to the two-dimensional example. The dashed lines are the exact functions, the solid lines are the approximations using the function 2 zðrÞ ¼ C 1 eC2 r .

application of interest—the generation of a water barrier—we are willing to trade some accuracy for much greater speed of prediction. For practical implementation, a relative error between the predicted surface and the desired shape of order e  10% or less is deemed acceptable. The important outcome is that the predicted surface should reach a certain minimum height and width, and that predictions should be made within 30 s or less. The desired surfaces in two dimensions span the lateral extent from 0 to 14 and are discretized into 141 grid points (i.e., a grid size of 0.1 is employed). The dimensionless maximum height of the surface is set to 0.3. In the first example, we consider a desired free surface that comprises a flat region of height 0.3 above the quiescent free surface, with a width of 10 on top and 14 at the base. Two sinusoidal shaped segments join the elevated flat

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

region on both sides. This target surface is shown as the dash line in Fig. 6(a). The reduced-order modeling algorithm is applied to find the initiation parameters for a set of six bubbles. The optimal solution to the bubble parameters and corresponding prediction of the free surface shape are presented in Fig. 6(b) and (a), respectively. As shown in Fig. 6(a), the main discrepancies lie primarily on the flat part of the surface. The local effects due to each bubble (that decay exponentially from each bubble location) can be seen clearly in the reconstructed surface shape. The reconstructed surface shape is limited by the minimum separation between bubbles and the number of bubbles used. Under these constraints, the predicted surface is considered to match the desired shape reasonably well, with a relative error of e = 8.8%, which is acceptable within the framework of the water barrier application. The

calculation for eight bubbles shows that the result is quite similar since, due to minimum separation requirements, the added two bubbles are located outside the domain of interest and hence their contribution is small. Fig. 6(b) presents the bubble parameters. The lateral position follows a linear distribution, meaning that the bubbles are located almost uniformly along the horizontal direction (r-axis) from r  0 to r  12. From the figure, it can also be seen that the bubbles are all located at approximately the same depth, except for the first one, which is cast deeper. Fig. 6(b) indicates that the first bubble is initiated earliest (at time zero); the others are released almost simultaneously after a delay of approximately 2 · 102 units of dimensionless time. The second example considers a convex sine-shape free surface which is depicted as the dash line in Fig. 7(a) with

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25

0.2

z 0.2

0.15

0.15

0.1

0.1

0.05

0.05

Z

0

0

2

4

6

8

10

12

507

0 0

14

2

4

6

8

r

r

(a)

(a)

10

12

14

16

16 horizontal position depth strength ε/100 inception time 100*t

14 12

12

10

10

8

8

6

6

4

4

2

2

0

0

–2

–2

–4

–4 1

2

3

4

5

horizontal position depth strength ε/100 inception time 100*t

14

6

bubble No.

(b) Fig. 6. Approximate optimal solution to two-dimensional flat shape free surface and bubble parameters for six bubbles. (a) Desired surface (dash) and predicted surface (solid), e = 8.8%. (b) Optimized values of bubble parameters.

1

2

3

4

5

6

bubble No.

(b) Fig. 7. Approximate optimal solution to two-dimensional convex sineshape free surface and bubble parameters for six bubbles. (a) Desired surface (dash) and predicted surface (solid), e = 4.1%. (b) Optimized values of bubble parameters.

508

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

a maximum height of 0.3 at r = 7 and a minimum height of 0.2 at r = 0 and r = 14. The resultant surface corresponding to the approximate optimal solution with six bubbles is shown in Fig. 7(a). Again, fluctuations due to local bubble effects are observed along the surface. With a relative error of e = 4.1%, this solution is considered acceptable. In Fig. 8(a), the predicted surface shows improvement with the addition of two more bubbles. With eight bubbles located in the lateral extent from 0 to 14, the spacing between adjacent bubbles is at the minimum allowable value of dr = 2.0. It can be seen that the agreement between the predicted and desired surface shapes is much better. The relative error is reduced to e = 1.2%. The parametric values for the sets of six and eight bubbles are presented in Figs. 7(b) and 8(b), respectively.

Fig. 9 shows the approximate optimal solution to the third example of a concave sine-shape free surface with a maximum height of 0.3 at the two ends and a minimum height of 0.2 in the center. This solution employs eight bubbles, and has an acceptable relative error of e = 2.8%. From the different kinds of free surface shapes considered, encompassing the broad extremes of flat, convex and concave profiles, the reduced-order approach to solving the optimization problem can be seen to give a viable solution for the surface when a sufficient number of bubbles is employed. Table 1 shows the corresponding computation time with different sizes of the problem for the case of the convex-shaped water barrier. (The timings for the other two water barriers are similar to this case.) This computation time includes those elements of the algorithm that

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25

z 0.2

z 0.2

0.15

0.15

0.1

0.1

0.05

0.05

0 0

2

4

6

8

10

12

0 0

14

2

4

6

8

r

r

(a)

(a)

10

12

14

7

8

16

16 horizontal position depth strength ε/100 inception time 100*t

14

horizontal position depth strength ε/100 inception time 100*t

14

12

12

10

10

8

8

6

6

4

4

2

2

0

0

–2

–2 –4

–4 1

2

3

4

5

6

7

8

bubble No.

(b) Fig. 8. Approximate optimal solution to two-dimensional convex sineshape free surface and bubble parameters for eight bubbles. (a) Desired surface (dash) and predicted surface (solid), e = 1.2%. (b) Optimized values of bubble parameters.

1

2

3

4

5

6

bubble No.

(b) Fig. 9. Approximate optimal solution to two-dimensional concave sineshape free surface and bubble parameters for eight bubbles. (a) Desired surface (dash) and predicted surface (solid), e = 2.8%. (b) Optimized values of bubble parameters.

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512 Table 1 Online computational time for solving the two-dimensional reduced-order optimization problem in the domain of [0, 14] (141 grid points) No. of bubbles

Online computation time (in seconds)

5 bubbles 6 bubbles 7 bubbles 8 bubbles Full simulation of 1 bubble

14 16 20 25 60 s

The solution is obtained by running a LOQO program on an Intel Xeon Dual CPU 2.80 GHz processor with 2.0 Gb RAM.

must be performed in the online phase, but excludes the cost of offline computations (e.g., the computation of the snapshot ensemble). With the current grid, computation time for a two-dimensional problem increases linearly as the number of bubbles increases. For comparison purposes, the computational time associated with a high-fidelity simulation of a single bubble is also given. It can be seen that the approximate optimization approach for prediction

509

of two-dimensional free surface shape is very efficient, and sufficiently rapid to be performed in real time. It is noted that the timing can meet the requirement of about 30 s (decision-making time) as enunciated in Section 3. The desired surface shape for the three-dimensional example is shown in Fig. 10(a). The surface considered is in the region of [0,8] · [0,8] of the x–y plane and is discretized into a mesh with element size of 0.5 in each direction. The method is applied to find the optimal solutions for sets of six, eight and ten bubbles. The predicted surface corresponding to the set of six bubbles is presented in Fig. 10(b). The interpretation of the bubble parameters given in Fig. 11 is similar to the two-dimensional cases except for the lateral positions, which are now presented in the x–y plane. The bubble distribution is symmetric, as is the free surface; the bubbles are essentially placed in two rows, and uniformly distributed within each row. Fig. 10(b) shows that six bubbles are not sufficient for creating the desired surface. There are many large crests and troughs on the predicted surface. The bubbles are located quite far from each other so that the free surface between

0.45

0.45

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25

z

z 0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0 8

0 8 8

6

2

2 0

6

4

4

y

8

6

6

4

4

y

2 0

0

(a)

x

2

x 0

(b)

0.45

0.45

0.4

0.4

0.35

0.35

0.3

0.3

0.25

0.25

z

z 0.2

0.2

0.15

0.15

0.1

0.1

0.05

0.05

0 8

0 8 8

6

y

8

6

6

4 4 2

2 0

0

(c)

y x

6

4 4 2

2 0

x

0

(d)

Fig. 10. (a) Desired surface shape for three-dimensional example, (b) predicted free surface shape for six bubbles, (c) predicted free surface shape for eight bubbles, (d) predicted free surface shape for ten bubbles.

510

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

9

9

8

8 4 3

7

y

6

6

6

5

5 2

4

y

5

3

2

2

1

3

7

2

6

4

3

1

8

7

1

4

1 0

–1 –1

5

0

0

1

2

3

4

5

6

7

8

9

–1 –1

0

1

2

3

4

x

x

(a)

(a)

12

6

8

8

6

6

4

4

2

2

0

0

–2

–2

–4

–4 2

3

4

5

8

9

depth strength ε/100 inception time 100*t

10

1

7

12

depth strength ε/100 inception time 100*t

10

5

6

bubble No.

(b)

1

2

3

4

5

6

7

8

bubble No.

(b)

Fig. 11. Parameter solution for six bubbles: (a) lateral positions, (b) other parameters: depth, strength, and inception time.

Fig. 12. Parameter solution for eight bubbles: (a) lateral positions, (b) other parameters: depth, strength, and inception time.

any two bubbles is depressed. The relative error of the approximation is e = 24.3%, which does not meet the accuracy requirements. Increasing the number of bubbles is expected to improve the solution. Fig. 10(c) shows the free surface corresponding to the set of eight bubbles; the parametric values for the bubbles are given in Fig. 12. Similar to the case of six bubbles, the bubble position and the free surface assume a fairly symmetric distribution; the bubbles sit uniformly in two rows. As shown in Fig. 10(c), the free surface follows reasonably the desired shape, although there are still some small troughs in the middle of the surface and between the crests. This solution, however, shows improvement over the previous one, with the relative error reduced to e = 16.9%. Fig. 10(d) presents the solution for the free surface associated with a set of 10 bubbles. The parametric values for

the bubbles are given in Fig. 13. It is clear that eight of these bubbles are located at similar positions corresponding to the solution for eight bubbles, where the remaining two bubbles are located on the two sides. This is due to the minimum spacing requirement of dr > 2.0, which prevents the additional bubbles from being located more closely. In addition, these two remaining bubbles do not make much contribution since they are located much deeper than the others. The resultant surface is flatter for the most part but the error of e = 17.1% represents a slight increase over the case of eight bubbles. This is due to the fact that, as can be seen in Fig. 10, the error falls mainly on the two sides of the surface (x 6 2 and x P 6). Further computations show that as more bubbles are used, the minimum spacing constraint prevents any further improvement in the solution.

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512

ever, as the number of bubbles increases beyond ten, the computation time increases substantially with the number of bubbles. Based on the current numerical experiments, for applications requiring ten or more bubbles a faster computer or parallelization of the method would be needed for the above three-dimensional optimization to meet the real-time requirements. Nevertheless, it should be noted that a highfidelity simulation of a single bubble in our example requires 280 s of computational time. The solution of the optimization problem, which requires numerous direct simulation runs, would be computationally intractable without the model reduction approach proposed here. Another computationally efficient option would be to calculate the three-dimensional solution as a series of two-dimensional water barriers. This approach could be easily parallelized.

9

8

4

8

3

7

9 7

6

5 y 4

3

2

6

1

5

2

1 10 0

–1 –1

0

1

2

3

4

5

6

7

8

511

9

x

(a)

5. Conclusions

12 depth strength ε/100 inception time 100*t

10

8

6

4

2

0

–2

–4 1

2

3

4

5

6

7

8

9

10

bubble No.

(b)

A new reduced-order optimization methodology is presented. High fidelity simulations are performed offline to generate an ensemble of POD basis functions. The realtime optimization is conducted by a two-stage approach using the POD coefficients as intermediate parameters. The real-time optimization using this new methodology yields approximate results, but reduces the computational time considerably. Results for two-dimensional simulations show that the critical time requirement (of about 30 s) for optimizing the shape of a water plume can be achieved. For three dimensional simulations, the time taken for computation increases, and problems using less than ten bubbles can be solved in the required time. For applications requiring ten or more bubbles computation times are still very fast, but do not meet the 30 s requirement. A viable alternative to the full three-dimensional optimization would be to carry out in parallel a series of two-dimensional optimizations.

Fig. 13. Parameter solution for ten bubbles: (a) lateral positions, (b) other parameters: depth, strength, and inception time.

Acknowledgement Table 2 shows the online computation time of the algorithm for this three-dimensional case. It can be seen that the algorithm performs reasonably well for moderate size problems in three dimensions (such as eight bubbles); howTable 2 Online computational time for solving the three-dimensional reducedorder optimization problem in the domain of [0, 8] · [0, 8] (17 · 17 grid points) No. of bubbles

Online computation time

6 bubbles 8 bubbles 10 bubbles 12 bubbles 14 bubbles Full simulation of 1 bubble

24 s 32 s 50 s 82 s 170 s 280 s

The solution is obtained by running a LOQO program on Intel Xeon Dual CPU 2.80 GHz processor with 2.0 Gb RAM.

This work is funded by the Singapore-MIT Alliance. References [1] Feldmann P, Freund RW. Efficient linear circuit analysis by Pade´ approximation via the Lanczos process. IEEE Trans Comput Aided Des Integrated Circ Syst 1995;14:639–49. [2] Gallivan K, Grimme E, Van Dooren P. Pade´ approximation of largescale dynamic systems with Lanczos methods. In: Proc of the 33rd IEEE conf on decision and control, December, 1994. [3] Moore BC. Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Trans Autom Control 1981;AC-26(1):17–31. [4] Gugercin S, Antoulas A. A survey of model reduction by balanced truncation and some new results. Int J Control 2004;77:748–66. [5] Li J, White J. Low rank solution of Lyapunov equations. SIAM J Matrix Anal Appl 2002;24(1):260–80. [6] Prud’homme C, Rovas D, Veroy K, Maday Y, Patera AT, Turinici G. Reliable real-time solution of parameterized partial differential

512

[7] [8] [9] [10] [11] [12] [13]

[14]

[15]

[16] [17]

[18]

D. My-Ha et al. / Computers & Fluids 36 (2007) 499–512 equations: Reduced-basis output bound methods. J Fluids Eng 2002;124:70–80. Lumley JL. The structures of inhomogeneous turbulent flow. Atmos Turbulence Radio Wave Propagat 1967:166–78. Sirovich L. Turbulence and the dynamics of coherent structures, Part 1: Coherent structures. Q Appl Math 1987;45(3):561–71. Kosambi DD. Statistics in function space. J Indian Math Soc 1943;7:76–88. Loe´ve M. Probability Theory. New York: D. Van Nostrand Company Inc.; 1955. Hotelling H. Analysis of a complex of statistical variables with principal components. J Educ Psychol 1933;24:417–41. p. 498–520. Dowell EH, Hall KC. Modeling of fluid-structure interaction. Annu Rev Fluid Mech 2001;33:445–9. Silveira LM, Kamon M, Elfadel I, White J. A coordinate-transformed Arnoldi algorithm for generating guaranteed stable reduced-order models of RLC circuits. Comput Meth Appl Mech Eng 1999;169(3–4): 377–89. Hinze M, Volkwein S. Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: Error estimates and suboptimal control. Technische Universitat Dresden, SFB609-Preprint-27-2004, 2004. Kunisch K, Volkwein S. Control of Burgers’ equation by reduced order approach using proper orthogonal decomposition. J Optim Theory Appl 1999;102:345–71. Everson R, Sirovich L. The Karhunen–Loe`ve for gappy data. J Opt Soc Am 1995;12:1657–64. Bui-Thanh T, Damodaran M, Willcox K. Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition. AIAA paper 2003, p. 4231. Willcox K. Unsteady flow sensing and estimation via the gappy proper orthogonal decomposition. Comput Fluids 2006;35(2):208–26.

[19] Epureanu BI, Dowell EH, Hall KC. A parametric analysis of reduced-order models of potential flows in turbomachinery using proper orthogonal decomposition. 2001-GT-0434. In: Proc of ASME TURBO EXPO 2001, New Orleans, Louisiana, June 2001. [20] Ly HV, Tran HT. Modeling and control of physical processes using proper orthogonal decomposition. J Math Comput Model. [21] Akcelik V, Bielak J, Biros G, Epanomeritakis I, Fernandez A, Ghattas O, et al., Terascale forward and inverse earthquake modeling. In: Proc of SC2003, 2003. [22] Wang QX, Yeo KS, Khoo BC, Lam KY. Strong interaction between a buoyancy bubble and a free surface. Theoret Comput Fluid Dyn 1996;8:73–88. [23] Wang QX, Yeo KS, Khoo BC, Lam KY. Nonlinear interaction between gas bubble and free surface. Comput Fluid 1996;25(7): 607–28. [24] Zhang YL, Yeo KS, Khoo BC, Chong WK. Simulation of threedimensional bubbles using desingularised boundary integral method. Int J Numer Meth Fluids 1999;31(8):1311–20. [25] Zhang YL, Yeo KS, Khoo BC, Chong WK. Three-dimensional computation of bubbles near a free surface. J Comput Phys 1998;146:105–23. [26] Rungsiyaphornrat S, Klaseboer E, Khoo BC, Yeo KS. The merging of two gaseous bubbles with an application to underwater explosion. Comput Fluid 2003;32:1049–74. [27] Best JP, Kucera A. A numerical investigation of nonspherical rebounding bubbles. J Fluid Mech 1992;245:137–54. [28] Lamb H. The early stages of a submarine explosion. Philos Mag 1923;45:257–65. [29] My-Ha D. Model Order Reduction for Prediction of Free Surface Shape due to Underwater Bubbles Dynamics. ME Thesis, National University of Singapore, Singapore, June 2004.