Reduction of Nonlinear Embedded Boundary Models for Problems with Evolving Interfaces

Reduction of Nonlinear Embedded Boundary Models for Problems with Evolving Interfaces

Maciej Balajewicz111Postdoctoral Fellow Charbel Farhat222Vivian Church Hoff Professor of Aircraft Structures Department of Aeronautics and Astronautics Department of Mechanical Engineering Institute for Computational and Mathematical Engineering
Stanford University, Stanford, CA 94305-4035, U.S.A
Abstract

Embedded boundary methods alleviate many computational challenges, including those associated with meshing complex geometries and solving problems with evolving domains and interfaces. Developing model reduction methods for computational frameworks based on such methods seems however to be challenging. Indeed, most popular model reduction techniques are projection-based, and rely on basis functions obtained from the compression of simulation snapshots. In a traditional interface-fitted computational framework, the computation of such basis functions is straightforward, primarily because the computational domain does not contain in this case a fictitious region. This is not the case however for an embedded computational framework because the computational domain typically contains in this case both real and ghost regions whose definitions complicate the collection and compression of simulation snapshots. The problem is exacerbated when the interface separating both regions evolves in time. This paper addresses this issue by formulating the snapshot compression problem as a weighted low-rank approximation problem where the binary weighting identifies the evolving component of the individual simulation snapshots. The proposed approach is application independent and therefore comprehensive. It is successfully demonstrated for the model reduction of several two-dimensional, vortex-dominated, fluid-structure interaction problems.

keywords:
data compression, embedded boundary method, evolving domains, immersed boundary method, interfaces, model reduction, singular value decomposition
journal: Journal of Computational Physics

1 Introduction

Two common computational frameworks for the solution of problems with evolving domains and interfaces are the Arbitrary Lagrangian Eulerian (ALE) Donea (); CF1 (); CF2 (); ALERezoning () and Embedded Boundary Method (EBM) frameworks. In an ALE framework, the computations are performed on an interface-fitted mesh which is deformed using a mesh motion algorithm MM1 (); MM2 (); MM3 () to follow the evolution of the interface. In an EBM framework, the interface is embedded in a background mesh --- also called an embedding mesh --- and allowed to intersect it as it evolves (Peskin_AN_2002, ; Mittal_ARFM_2005, ; Wang_IJNMF_2011, ; FIVER, ). In the latter case, the interface separates the computational domain in two regions: a ‘‘ghost’’ (or fictitious) region usually associated with the volume delimited by the embedded interface (boundary, or body), and a ‘‘real’’ region usually associated with the remainder of the computational domain. Both computational frameworks have advantages and shortcomings. An ALE framework is particularly advantageous for problems where the evolution of the interface is characterized by small-amplitude motions and deformations. An EBM framework is often indispensable for problems where the interface undergoes large-amplitude motions and deformations, and/or topological changes.

Irrespectively of the chosen framework however, the computational cost of high-fidelity simulations of problems with evolving domains and interfaces can be prohibitive. For this reason, interest in Model Order Reduction (MOR) methods continues to grow (Rewienski_LAA_2006, ; Thuan1, ; Amsallem1, ; Astrid_AC_2008, ; Chaturantabut_JSC_2010, ; Carlberg1, ; Balajewicz_JFM_2013, ; Carlberg_JCP_2013, ; Navon1, ; Navon2, ). Most of these methods are projection-based. Therefore, they rely on the computation of suitable basis functions. Typically, these are constructed by computing solution snapshots of some instances of the family of problems of interest and compressing them using the Proper Orthogonal Decomposition (POD) or Singular Value Decomposition (SVD) method. In an interface-fitted framework, this computation is straightforward because the computational domain does not contain a ghost region. In an EBM framework however, this computation is complicated by the presence of a ghost region where the solution snapshots may or may not be populated, and in the latter case, the populated values may or may not have a physical meaning. The situation is further exacerbated when the interface evolves in time, as in this case the partitioning into meaningful and meaningless components of the collected snapshots also evolves in time, which complicates further the compression of the collected snapshots.

Therefore, the main objective of this paper is to propose an effective method for the compression of solution snaphsots collected during the simulation of problems with evolving domains and interfaces using an EBM, in view of enabling the projection-based reduction of high-fidelity embedded boundary computational models. This method consists of formulating the data compression problem as a weighted low-rank matrix approximation problem, where the weighting is binary and identifies the ghost/real partition of the individual snapshots. Therefore, the basis functions derived using this approach are optimal for the real component of the partition of the computational domain.

The weighted low-rank approximation problem formulated in this paper is pervasive. It can be found, for example, in:

• Computer vision applications with occlusion.

• Signal processing applications with irregular measurements in time and space.

• Control problems with malfunctioning measuring devices.

This problem is a generalization of the popular compressed sensing problem (Donoho_2006, ; Candes_2008, ). It also shares many similarities with the matrix completion problem (Candes_2009, ; Cai_2010, ; Jain_ACM_2013, ; Vandereycken_SIAM_2013, ). For all these reasons, many algorithms for the solution of this problem are readily available (for example, see (Srebro_ICML_2003, ; Buchanan_2005, ; Chen_2008, )).

Whereas the approach outlined above for compressing solution snaphots computed by an EBM in view of constructing a reduced-order basis suitable for model reduction is application independent and therefore quite general, it is developed and demonstrated in the remainder for this paper for Fluid-Structure Interaction (FSI) problems. To this effect, the remainder of this paper is organized as follows.

In Section 2, the main notation used in this paper is laid out to facilitate its reading. In Section 3, a general EBM framework for Computational Fluid Dynamics (CFD) is overviewed, and the standard projection-based MOR approach is recapitulated. In Section 4, the proposed method for compressing CFD snapshots computed by an EBM is introduced. In Section 5, this method is applied to three different FSI problems. Finally in Section 6, conclusions are offered and prospects for future work are summarized.

2 Notations and definitions

Throughout this paper, matrices are denoted by bold capitals (ex. ), vectors by bold lower case (ex. ), and subscripts identify rows and columns (for example, is the entry at the -th row and -th column of matrix ). The ‘‘colon’’ notation is used to specify columns or rows of a matrix (for example denotes the -th column and denotes the first rows of matrix ). and denote the identity and null matrices, respectively.

For two matrices and of equal dimension , the Hadamard product is the matrix of the same dimension whose elements are given by

 (A⊙B)i,j=Ai,j⋅Bi,j (1)

If is an matrix and is a matrix, the Kronecker matrix product is the block matrix

 A⊗B=⎡⎢ ⎢ ⎢⎣A1,1B⋯A1,nB⋮⋱⋮Am,1B⋯Am,nB⎤⎥ ⎥ ⎥⎦ (2)

The operator denotes the vectorization of the matrix formed by stacking the columns of into a single column vector

 vec(A)=⎡⎢ ⎢⎣A:,1⋮A:,m⎤⎥ ⎥⎦ (3)

For a vector of dimension , the operator returns the diagonal matrix

 diag(x)=⎡⎢ ⎢⎣x1⋯0⋮⋱⋮0⋯xn⎤⎥ ⎥⎦ (4)

For a matrix that is , the operator returns the diagonal matrix

 diag(A)=⎡⎢ ⎢ ⎢⎣diag(A:,1)⋯0n⋮⋱⋮0n⋯diag(A:,m)⎤⎥ ⎥ ⎥⎦ (5)

The standard Euclidean norm of a vector and Frobenius norm of a matrix are denoted by and , respectively, and defined as follows

 ∥x∥2=(n∑i=1x2i)12,∥A∥F=(n∑i=1m∑j=1A2i,j)12 (6)

Two of the applications considered in this paper involve compressible flows governed by the Navier-Stokes equations. In this case, the fluid state vector is denoted by , and the fluid pressure, density, and velocity magnitude are denoted by , , and , respectively. The free-stream conditions are designated by the subscript .

3 Reduction of embedded boundary models

In this work, FSI problems are considered as a vehicle to explain the snapshot compression method proposed for constructing a Reduced-Order Basis (ROB) suitable for the nonlinear reduction of embedded boundary models. It is stressed however that this method is general. In principle, it applies to any EBM framework designed for solving any computational problem. Before presenting this method in Section 4 however, the stage is set here by describing the generic problem of interest, and its reduction process.

3.1 Generic nonlinear fluid-structure interaction problem

Consider the problem of computing an unsteady compressible flow around a rigidly moving body contained in a fixed region . A representative geometry of this problem is illustrated in the left part of Figure 1. Assume that the boundary and fluid/body interface are equipped with the appropriate boundary and fluid-structure transmission conditions, respectively. In an EBM framework, a background Eulerian mesh is typically used to discretize the global domain . The nodes of this mesh that are occluded by at time are usually referred to as ‘‘ghost’’ nodes; the others are referred to as ‘‘real’’ nodes. A semi-discretization on this mesh of the Navier-Stokes equations governing the fluid subsystem yields a set of nonlinear ordinary differential equations (ODEs) which can be written as

 dwdt+f(w)=0 (7)

where denotes time, , is the number of fluid variables per mesh node, denotes the total number of mesh nodes, and denotes a nonlinear function of that represents the convective and diffusive fluxes.

Without any loss of generality, it is assumed throughout the remainder of this paper that Eq. (7) is discretized in time using an implicit linear multistep scheme. Hence, if denotes a discretization of the time-interval and , the discrete counterpart of Eq. (7) at time-step is

 R(wn)=s∑j=0αjwn−j+s∑j=0βjf(wn−j)=0 (8)

where is the order of accuracy of the chosen time-integrator and and are two constants that characterize it.

In an EBM framework for CFD, the fluid-structure transmission conditions require special attention because the fluid/body interface does not coincide in general with the nodes of the background mesh. Over the years, a large number of different boundary treatment schemes has been developed for addressing this issue (for example, see  Mittal_ARFM_2005 () and Zeng12 () for a review of this and related issues). This is noteworthy because the method proposed in this paper for compressing solution snapshots computed by EBMs based on these schemes is actually independent of them.

3.2 Nonlinear model reduction

In a projection-based MOR method, the state vector is approximated in a global affine trial subspace as follows

 wn≈˜w:=w0+Uan (9)

where is a matrix whose columns contain the basis of this subspace, is the reduced dimension, denotes the generalized coordinates in this basis at time , and is the initial condition. Substituting Eq. (9) into Eq. (8) yields , which represents an overdetermined system of equations with unknowns. Consequently, the residual is projected onto a test subspace represented by , yielding the square system . In a Galerkin projection, . In a least-square projection,  Carlberg1 (), where is the Jacobian of  (8) with respect to . In either case, the vector of generalized coordinates is then obtained by solving a square system of nonlinear equations using Newton’s method or a preferred variant.

Computing the vector of generalized coordinates using a least-squares projection approach is equivalent to solving at each time-step the minimization problem Carlberg1 ()

 minan∥∥R(w0+Uan)∥∥2 (10)

For nonlinear, non self-adjoint problems such as those represented in this case by the set of ODEs (7), this approach is more robust than its Galerkin counterpart. Hence, its application in the context of an EBM framework is given here further attention.

Since only a portion of the Eulerian computational domain corresponds to the real flow region, minimizing the norm of the global residual vector --- that is, including its ghost component --- is neither necessary nor convenient. Instead, let denote a binary vector identifying the fluid/body --- or more generally, the ghost/real --- partition of the computational domain at time

 mni={1,if i∈Ω∖¯¯¯¯B(tn)0,if i∈B(tn) (11)

The idea here is that at each time-step, the norm of only the real component of --- that is, that associated with the real flow region --- needs to be minimized to obtain the vector of generalized coordinates needed to determine the Reduced-Order Model (ROM) approximation (9). In other words, the idea here is that in the context of an EBM, the minimization problem (10) can be reformulated as

 minan∥∥mn⊙R(w0+Uan)∥∥2 (12)

Equation (12) above is a nonlinear equation. Therefore, it must be solved iteratively using, for example, the Gauss-Newton method, which can be summarized as

for solve

 minΔa(i)∥mn⊙(J(i)UΔa(i)+R(Uan,(i)))∥2 (13)

where the superscript designates the -th iteration, is the increment of the sought-after solution at the -th Gauss-Newton iteration, and is determined by a convergence criterion.

Equation (13) constitutes a -dimensional Gauss-Newton ROM of the High-Dimensional Model (HDM) --- that is, the high-fidelity CFD model --- represented here by Eq. (7). The normal equations associated with Eq. (13) can be written as

 (J(i)U)Tdiag(mn)J(i)UΔa(i)=−(J(i)U)Tdiag(mn)R(Uan,(i)) (14)

where the superscript denotes the transpose.

In any case, whether a Galerkin or least-squares approach is chosen for determining the vector of generalized coordinates , the time-invariant ROB must be first determined. In MOR applications, basis functions are usually constructed from solution snapshots. These snapshots are collected, assembled into a matrix, and compressed using, for example, SVD. In an interface-fitted computational framework, this computation is straightforward because the computational domain does not contain a ghost region. In an EBM framework however, this compuation is complicated by the presence and evolution of a ghost/real partition of the computational fluid domain. In the following section, a method for constructing optimal fluid basis functions from snapshots computed by an EBM is presented.

4 Construction of a reduced-order basis for an embedded boundary model

To this effect, this section is organized in three parts. First, the concepts of solution snapshots and a snapshot matrix are recalled. Then, the problem of constructing an optimal ROB for embedded boundary models is formulated as a weighted low-rank matrix approximation problem. Finally, an efficient iterative scheme for solving this problem is specified.

4.1 Eulerian snapshots

A solution snapshot, or simply a snapshot, is defined here as a state vector computed as the solution of Eq. (8) for some instance of its parameters --- that is, for some specific time or specific value of the set of flow parameters or boundary/initial conditions underlying this governing equation --- on the background Eulerian mesh discretizing the domain . A snapshot matrix is defined as a matrix () whose columns are individual snapshots. For all practical purposes, the main focus of this paper is on unsteady flows and on snapshots associated with different time-instances . Hence, for .

Because the proposed method for computing a ROB for an embedded boundary model does not utilize information from the occluded region of the computational domain, the ghost components of these snapshots can take arbitrary values. Consequently, the proposed method, which is described below, is independent of the specifics of the EBM framework used to generate the snapshots.

4.2 Weighted low-rank matrix approximation problem

Constructing basis functions from snapshots in the spirit of the POD method --- that is, using data compression --- can be formulated mathematically as a low-rank matrix approximation problem as follows.

For a given snapshot matrix , find a lower rank matrix that solves the minimization problem

 minrank(˜X)=k∥X−˜X∥F (15)

where . In this problem, the rank constraint can be taken care of by representing the unknown matrix as , where and , so that problem (15) becomes

 minU∈RcN×k,V∈Rk×K∥X−UV∥F (16)

It is well-known that the solution of the above low-rank approximation problem is given by the SVD of : specifically, and where .

Unfortunately, when the snapshot matrix is generated using an EBM framework, the solution outlined above cannot be expected to yield an optimal ROB. This is because: (a) the snapshots computed using an EBM contain information from both the flow (or real) region of the computational domain and its occluded (or ghost) region, and (b) this data inconsistency is not accounted for in the standard low-rank matrix approximation problem (15). Hence, this issue is resolved here by proposing the alternative weighted low-rank matrix approximation problem

 minU∈RcN×k,V∈Rk×K∥M⊙(X−UV)∥F (17)

where is a binary mask matrix identifying the ghost/real partition of the snapshot matrix

 M:=[m0m1⋯mK] (18)

where has already been defined in (11). This binary weighting ensures that the values of at the ghost nodes do not play a role in the construction of the basis functions, which in turn implies that the derived basis functions are optimal for the real components of the snapshots associated with .

In the general case (i.e. ), the weighted low-rank approximation problem formulated above is not reducible to the un-weighted problem. Therefore, its solution is not given by the SVD factorization of the snapshots 333The weighted low-rank matrix approximation problem is reducible to the un-weighted problem only for the special case where ; see B for details.. Furthermore, no closed form solution of this problem is known. Hence, it must be solved by numerical iterations. In this work, the Alternating Least Squares (ALS) algorithm is chosen for this purpose. This algorithm is the most widely applicable and empirically successful approach for solving this and related problems (Chen_2008, ; Okatani_2007, ; Mitra_2010, ).

4.3 Alternating least squares algorithm

The ALS algorithm takes advantage of the bi-linearity of the representation . Using the connection between the Frobenius and Eucledian norms and the notation for the Kronecker matrix product (2), problem (17) can be re-written as

 minU∈RcN×k,V∈Rk×K∥M⊙(X−UV)∥F ⇔minU∈RcN×k,V∈Rk×K∥diag(M)(vec(X)−vec(UV))∥2 (19a) ⇔minU∈RcN×k,V∈Rk×K∥diag(M)(vec(X)−(IK⊗U)vec(V))∥2 (19b) ⇔minU∈RcN×k,V∈Rk×K∥∥diag(M)(vec(X)−(VT⊗IcN)vec(U))∥∥2 (19c)

Problem (19b) is a linear least-squares problem for if is known. Problem (19c) is a linear least-squares problem for if is known. This suggests the following simple iterative solution procedure:

1. Guess (for example, initialize as the first left-singular vectors of )

2. Repeat until convergence

1. Solve problem (19b) for given

2. Solve problem (19c) for given

A more detailed outline of the above algorithm is given in Algorithm 1, where the stopping criterion has been omitted for the sake of brevity (see (Markovsky_book_2012, ; Markovsky_JCAM_2013, ) for a discussion of this topic).

The superscripts of and in Algorithm 1 designate an ALS iteration. The computational cost of this algorithm is roughly times the computational cost of a single thin SVD of the snapshot matrix times the number of iterations for convergence --- that is, , where denotes the performed number of iterations. In Section 5, it is shown that in general, suffices to obtain good results. The convergence of this algorithm is monitored using the normalized weighted projection error

 e(i):=∥M⊙(X−U(i)V(i))∥2F/∥M⊙X∥2F (20)

This error is guaranteed to decrease monotonically to a local minimum (Markovsky_book_2012, ; Markovsky_JCAM_2013, ).

A MATLAB implementation of Algorithm 1 is provided in A.

5 Applications

Three FSI problems are considered here to illustrate the method proposed in this paper for constructing basis functions suitable for the nonlinear model reduction of embedded boundary models, and demonstrate its features and performance. The first one is a simple one-dimensional model problem based on Burgers’ equation. It has the merit of illustrating the basic idea and being easily reproducible by the interested reader. The second problem focuses on the prediction of a two-dimensional unsteady viscous flow around a cylindrical body in large-amplitude heaving motion. Because this problem is formulated in an unbounded fluid domain, it is also solvable using an interface- (or body-) fitted ALE framework. Hence, it offers a venue not only for assessing the intrinsic performance of the proposed method for constructing a suitable ROB for an embedded boundary model that is more representative than that of the previous example, but also assessing the performance of the resulting nonlinear ROM relative to that of a counterpart ROM constructed using a body-fitted computational framework. The third problem focuses on the solution of a two-dimensional unsteady turbulent flow inside a square cavity containing a rotating ellipsoidal body. In this case, the large-amplitude motion of challenges the robustness if not feasibility of a body-fitted ALE framework and calls instead for an EBM framework. Hence, this third problem highlights the need for and demonstrates the performance of the computational methodology proposed in this paper.

For each FSI problem outlined above, the governing (fluid) equations are semi-discretized using the central finite difference method. They are discretized in time using the backward Euler implicit scheme with a constant time-step chosen so that does not travel more than one layer of nodes during this time-step. The Newton method is used to solve all nonlinear equations arising from the implicit time-discretization.

Whenever the flow problem of interest is solved using a body-fitted ALE framework, a ROB for this problem is constructed by computing snapshots of the solution at different time instances , then compressing them using the SVD. In this case, the nonlinear ROM of interest is generated using the GNAT method Carlberg1 (); Carlberg_JCP_2013 () based on a least-squares projection.

On the other hand, whenever an EBM framework is used to solve the flow problem of interest, a ROB for this problem is constructed by computing snapshots of the solution at different time instances , and compressing them using the ALS algorithm described in Algorithm 1. In this case, the nonlinear ROM of interest is constructed using a variant of the same GNAT method where at each time-step, the vector of generalized coordinates is computed by solving (12).

5.1 One-dimensional fluid-structure model problem based on the viscous Burgers equation

First, the following one-dimensional FSI problem based on the viscous Burgers equation, a periodic boundary condition for the fluid subsystem, and non-homogeneous Dirichlet boundary conditions for the body subsystem in lieu of the fluid-structure transmission conditions is considered

 (21)

where is a Reynolds-like number and is set to , and

 ξ(t)=(1−0.1)/2+0.1sin(2πt) (22)

The above initial boundary value problem models a flow problem in an unbounded, one-dimensional, fluid domain in the middle of which a rigid, linear body of length is placed and set into the oscillatory harmonic motion characterized by the amplitude 0.1, frequency , and velocity . At each time , defines the position of the left extremity of , and defines the position of its right extremity. Hence, this problem is a one-dimensional instance of the generic FSI problem discussed in Section 3.1.

The simplicity of problem (21) is such that it can be solved using both an ALE framework on a moving mesh, and an EBM framework on a fixed mesh. Hence, both frameworks are considered here, primarily for the purpose of comparisons. In both cases, Eq. (21) is discretized on a uniform Cartesian mesh with ( nodes).

Specifically, due to the unbounded fluid domain assumption and the periodic fluid boundary condition, the body-fitted ALE framework considered here for the solution of problem (21) is equipped with a rigid motion of the mesh that is identical to that of . As for the considered EBM framework, it is equipped with a first-order ghost fluid scheme where the ghost values of are populated at each time-instance only at the two ghost fluid nodes that are the nearest to the left and right boundaries of , using the same value . This EBM framework is illustrated in Figure 2, where the circles filled in black represent the real nodes, the empty ones represent the ghost fluid nodes, and the empty squares designate the subset of ghost fluid nodes whose ghost fluid values are populated. Both computational frameworks deliver the same HDM solution which is graphically depicted in Figure 3.

For each of the body-fitted ALE and EBM models, three ROBs of dimension , 20, and 40 are constructed using snapshots of the solution of problem (21), and three corresponding nonlinear ROMs are generated.

Figure 4(a) reports on the convergence of the ALS algorithm using the normalized weighted projection error (20). It illustrates a well-known property of this algorithm, namely, its guarantee to deliver a monotonic convergence to a local minimum. Figure 4(b) reports on the variation with the ROM dimension of the EBM ROM error defined as

 e(i)ROM:=∥M⊙(X−U(i)A)∥2F/∥M⊙X∥2F (23)

where is the matrix of snapshot solutions of problem (21), is the ROB computed at the -the iteration of the ALS algorithm, and is the matrix of generalized coordinates of the approximate solutions associated with and can be written as

 A:=[a0a1⋯aK] (24)

--- that is, each -th column of the matrix product is the EBM ROM solution of problem (21) at time , based on the ROB computed at the -th iteration of Algorithm 1.

At the initial (-th) ALS iteration, the ROB is constructed by performing the SVD directly on EBM solution snapshots where the ghost values of the solution --- that is, the values of this solution at the occluded mesh nodes --- are populated with zeros. As shown in Figure 4(b), this approach --- referred to here as the ‘‘naive approach’’ --- delivers a poor performance, which illustrates the need for an alternative method for constructing EBM ROBs. Figure 4(b) also reveals that a single ALS iteration significantly improves the computed EBM ROM solution. The results reported in Figure 4(b) also show that after a single ALS iteration, the constructed EBM ROMs deliver already a comparable accuracy to that of the body-fitted ALE ROMs of same dimensions.

Figure 5 displays the convergence of the EBM ROM solutions at for , and contrasts these solutions with their counterparts obtained using the body-fitted ALE ROMs of same dimensions. The reader can observe that both types of ROMs exhibit comparable convergence behavior and accuracy, thereby demonstrating the effectiveness of the proposed method for constructing suitable EBM ROBs.

5.2 Two-dimensional fluid-structure interaction problem in an unbounded fluid domain

Next, the simulation of a compressible viscous flow around a heaving rigid cylindrical body with a circular cross section of radius is considered. The fluid is assumed to be a perfect gas with the ratio of specific heats and is modeled using the two-dimensional compressible Navier-Stokes equations. The cylinder is assumed to be infinitely long, and therefore the flow is effectively modeled as a two-dimensional problem around a disk. The physical fluid domain is assumed to be unbounded, but the computational fluid domain is bounded by a square of edge length equal to as shown in Figure 6. The free-stream Mach number is set to , the Reynolds number based on the cylinder diameter to , and the Prandtl number to .

The center of the disk representing , , is set into the vertical sinusoidal motion

 yO(t)=2rsin(t)=sin(t) (25)

where denotes the abscissae of . Hence, the total center-to-center displacement of is , which is a relatively large displacement. Yet, because the fluid domain is also unbounded in this case, and because is rigid, the FSI problem outlined above can be reliably solved by a body-fitted ALE framework where the CFD mesh is rigidly moved to follow the vertical oscillations of .

Hence, two similar but strictly speaking different meshes are constructed for discretizing the spatial domain shown in Figure 6. The first mesh is designed for the EBM framework. It is non-uniform, has elements whose size in the vicinity of the cylindrical body is uniform and at its finest is characterized by , and embeds the circular body . The second mesh is similar except for the fact that it is a body-fitted mesh.

The EBM for CFD used for solving this FSI problem is illustrated in Figure 7, where the real fluid nodes are shown as solid circles and the ghost fluid nodes as empty circles. The kinematic transmission condition (or adherence condition) on the fluid/body interface is enforced using a first-order ghost-fluid scheme. At each time-step, the ghost fluid values at the ghost fluid nodes that are the nearest to the fluid/body interface (these are identified by squares in Figure 7), are populated as follows. The ghost fluid velocities at these nodes are set to the velocity of the immersed body . The ghost fluid pressures and densities at these same nodes are computed by constant extrapolation of their counterparts at the nearest neighboring real fluid nodes. On the other hand, the ghost fluid values at all other ghost nodes, which are labelled ‘‘G’’ in Figure 7, are not populated. At the end of each computational time-step, the fluid state vector at the ‘‘transitional’’ nodes (labelled ‘‘T’’ in Figure 7) --- that is, the nodes whose status switches from ghost to real --- are populated using the ghost fluid values at these nodes from the previous time-step. Snapshots of the vorticity and pressure fluctuations computed using this EBM and the HDM are shown in Figure 8.

First, a pair of unsteady simulations are performed for using the body-fitted ALE and EBM frameworks outlined above, in order to generate in each case snapshots for model reduction. Then, ROBs and nonlinear ROMs are constructed as outlined above.

Figure 9(a) reports on the convergence of the ALS algorithm using the projection error  (20). As expected, this error decreases monotonically.

For this problem, the ROM error is assessed using the relative error in the lift coefficient measured in the norm

 e(i)L:=∥CCFDL−CROML∥2/∥CCFDL∥2 (26)

where and . This error is reported in Figure 9(b). As before, the ROB obtained at the -th iteration of the ALS algorithm is computed using the so-called naive approach. Figure 9(b) shows that as expected, the nonlinear ROM constructed using this ROB performs poorly. However, after only one ALS interation, the ROM constructed using the ROB delivers a good performance. The results reported in Figure 9(b) also demonstrate that after a single ALS iteration, the generated EBM ROMs deliver an accuracy that is comparable to that of their ALE ROM counterparts of the same dimensions.

Figure 10 displays the evolution of the instantaneous coefficients of lift and drag, and , where and . The solid curves correspond to computations performed using the EBM HDM, whereas the dashed curves correspond to counterparts performed using the constructed nonlinear ROMs (and a single iteration of the ALS algorithm in the case of the EBM ROMs). The reader can observe that, as expected, the solutions predicted using both sets of ROMs exhibit comparable accuracies and converge to their counterparts obtained using the EBM HDM.

5.3 Two-dimensional turbulent fluid-structure interaction problem in a cavity

Finally, the two-dimensional computation of a viscous flow inside a square cavity with a rotating ellipsoidal body is considered here. Specifically, the objective is to predict the flow past the ellipsoidal body during the time-interval in which it rotates a full degrees. This FSI problem with large-amplitude displacements and rotations is chosen because unlike the previous two problems, it cannot be solved using a body-fitted ALE framework without some form of repeated remeshing to avoid mesh entangling. Hence, it is representative of a large family of FSI problems for which the EBM framework is preferred, if not essential. The geometry of this problem is illustrated in Figure 11.

Again, the fluid is modeled as a perfect gas with the ratio of specific heats . The Reynolds number based on the ellipse tip is set to and the Prandtl number to . The flow is assumed to be initially quiescent.

The square domain is discretized using a uniform Cartesian grid with elements. The same EBM outlined for the previous FSI problem with a heaving cylindrical body is applied to its solution. Computational snapshots of the vorticity contours obtained using the EBM HDM and three different EBM ROMs of dimensions varying between and are illustrated in Figure 12. All EBM ROMs are constructed with a ROB obtained using a single iteration of the ALS algorithm.

As in the previous two cases, the EBM ROM solution is shown to converge to its EBM HDM counterpart when the dimension of the ROM is increased. In particular, the reader can observe in Figure 12(c) that the vorticity field predicted by the EBM ROM of dimension reproduces remarkably well that obtained using the EBM HDM. For , the vorticity fields predicted by the EBM ROM and EBM HDM are almost indistinguishable, thereby demonstrating the effectiveness of the proposed method for constructing ROBs and ROMs for embedded boundary models.

6 Conclusions

Embedded boundary models are popular for the solution of nonlinear problems with evolving domains and/or moving interfaces. The application of projection-based model order reduction methods to the construction of reduced-order versions of such models calls for the construction first of their underlying reduced-order bases. The collection and compression of solution snapshots using various instances of a computational model of interest is both a popular and effective approach for generating reduced bases. For a traditional interface-fitted computational model, it is a straightforward procedure that has been popularized by the Proper Orthogonal Decomposition and Singular Value Decomposition methods. For an embedded boundary computational framework, this approach faces the problem of dealing with the partitioning of the overall computational domain into a real subdomain and a ghost one corresponding to the region of the overall computational domain that is occluded by the embedded boundary. For time-dependent problems, this complication is further exacerbated by the evolution in time of this partitioning as it implies an evolution in time of the partitioning between meaningful (real) and meaningless (ghost) entries of the collected solution snapshots. This makes the problem of compressing these snapshots in view of constructing a reduced-order basis a difficult task. To address this issue, this paper formulates the snapshot compression problem as a weighted low-rank approximation problem where the binary weighting identifies the evolving component of the individual simulation snapshots, and proposes to solve this problem using the Alternating Least Squares (ALS) algorithm. This approach is applicable in principle to any embedded boundary value problem. It is successfully demonstrated in this paper for three different fluid-structure interaction problems of increasing complexity. In all considered cases, it is shown to deliver reduced-order bases and models that effectively reproduce the high-dimensional solutions even when the flow is vortex-dominated, and the immersed body underoges large displacements and rotations.

Appendix A MATLAB implementation of the ALS algorithm

The following is a simple MATLAB implementation of Algorithm 1 presented in Section 4.3. The inputs of the function ALS are the snapshot matrix X where the ghost values are set to zero, a binary weighting matrix M, the desired dimension of the subspace of approximation k, and the maximum number of iterations it_max. The outputs of this function are the low-rank approximation matrices U and V.

1 function [U,V] = ALS(X,M,k,it_max)
2    [cN, K] = size(X);
3    [U_star, ~, ~] = svd(X,'econ');
4    U = U_star(:,1:k);
5    clear U_star;
6    V = zeros(k,K);
7    for p = 1:it_max
8        for j = 1:K
9            J = find(s(:,j));
10            x = full(X(J,j));
11            u = U(J,:);
12            V(:,j) = u\x;
13        end
14        for i = 1:cN
15            I = find(M(i,:));
16            x = full(X(i,I));
17            v = V(:,I);
18            U(i,:) = x/v;
19        end
20    end

Appendix B Solution of the weighted low-rank matrix approximation problem when rank(M)=1

If is a rank 1 weighting matrix, it can be written as for some and . In this case,

 ∥M⊙(X−˜X)∥2F =∑i,jsitj(X−UV)2i,j =∑i,j(√sitjXi,j−(√siUi,:)(√tjV:,j))2 (27a)

Defining as , as , and as , the above result can be re-written as

 ∥M⊙(X−˜X)∥2F =∥(X′−˜X′)∥2F (28a)

where . This shows that when is a rank 1 weighting matrix, the weighted low-rank approximation problem (17) reduces to its un-weighted counterpart (15).

An example of a rank 1 weighting matrix is the weighting matrix appropriate for an embedded boundary model where the embedded body is stationary --- that is, --- but the flow is unsteady. In this case,

 M:=[m∗m∗⋯m∗] (29)

Acknowledgments

The authors acknowledge partial support by the Office of Naval Research under Grant N00014-11-1-0707, and partial support by the Army Research Laboratory through the Army High Performance Computing Research Center under Cooperative Agreement W911NF-07-2-0027. The content of this publication does not necessarily reflect the position or policy of any of the aforementioned institutions, and no official endorsement should be inferred.

References

• (1) J. Donea, An arbitrary Lagrangian-Eulerian finite element method for transient fluid-structure interactions, Computer Methods in Applied Mechanics and Engineering 33 (1982) 689--723.
• (2) C. Farhat, M. Lesoinne, N. Maman, Mixed explicit/implicit time integration of coupled aeroelastic problems: three-field formulation, geometric conservation and distributed solution, International Journal for Numerical Methods in Fluids 21 (1995) 807--835.
• (3) C. Farhat, A. Rallu, K. Wang, T. Belytschko, Robust and provably second-order explicit-explicit and implicit-explicit staggered time-integrators for highly nonlinear fluid-structure interaction problems, International Journal for Numerical Methods in Engineering 84 (2010) 73--107.
• (4) R. Loubère, P.H. Maire, M. Shashkov, J. Breil, S. Galéra, ReALE: A reconnection-based arbitrary-Lagrangian-Eulerian method, Journal of Computational Physics 229 (2010) 4724--4761.
• (5) C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for two-dimensional dynamic unstructured fluid meshes. Computer Methods in Applied Mechanics and Engineering 163 (1998) 231--245.
• (6) C. Degand, C. Farhat, A three-dimensional torsional spring analogy method for unstructured dynamic meshes. Computers and Structures 80 (2002) 305--316.
• (7) M. Fossati, R.A. Khurram, W.G. Habashi, An ALE mesh movement scheme for long-term in-flight ice accretion, International Journal for Numerical Methods in Fluids 68 (2012) 958--976.
• (8) C. S. Peskin, The immersed boundary method, Acta Numerica 11 (0) (2002) 479--517.
• (9) R. Mittal, G. Iaccarino, Immersed boundary methods, Annual Review of Fluid Mechanics 37 (2005) 239--261.
• (10) K. Wang, A. Rallu, J. F. Gerbeau, C. Farhat, Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluid--structure interaction problems, International Journal for Numerical Methods in Fluids 67 (9) (2011) 1175--1206.
• (11) C. Farhat, J.-F. Gerbeau, A. Rallu, FIVER: a finite volume method based on exact two-phase riemann problems and sparse grids for multi-material flows with large density jumps, Journal of Computational Physics, 231 (2012) 6360--6379.
• (12) M. Rewieński, J. White, Model order reduction for nonlinear dynamical systems based on trajectory piecewise-linear approximations, Linear algebra and its applications 415 (2) (2006) 426--454.
• (13) T. Lieu, C. Farhat, M. Lesoinne, Reduced-order fluid/structure modeling of a complete aircraft configuration, Computer Methods in Applied Mechanics and Engineering 195 (2006) 5730--5742.
• (14) D. Amsallem, C. Farhat, An interpolation method for adapting reduced-order models and application to aeroelasticity AIAA Journal, 46 (2008) 1803--1813.
• (15) P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing point estimation in models described by proper orthogonal decomposition, Automatic Control, IEEE Transactions on 53 (10) (2008) 2237--2251.
• (16) S. Chaturantabut, D. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing 32 (5) (2010) 2737--2764.
• (17) K. Carlberg, C. Bou-Mosleh and C. Farhat, International Journal for Numerical Methods in Engineering, 86 (2) (2011) 155--181.
• (18) M. J. Balajewicz, E. H. Dowell, B. R. Noack, Low-dimensional modelling of high-reynolds-number shear flows incorporating constraints from the Navier--Stokes equation, Journal of Fluid Mechanics 729 (2013) 285--308.
• (19) K. Carlberg, C. Farhat, J. Cortial, D. Amsallem, The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows, Journal of Computational Physics 242 (2013) 623--647.
• (20) J. Du, F. Fang, C.C. Pain, I.M. Navon, J. Zhu, D.A. Ham, POD reduced-order unstructured mesh modeling applied to 2D and 3D fluid flow, Computers and Mathematics with Applications, 65 (2013) 362--379.
• (21) D. Xiao, F. Fang, A. G. Buchan, C.C. Pain, I.M. Navon, J. Du, G. Hu Non-linear model reduction for the Navier-Stokes equations using residual DEIM method, Journal of Computational Physics, 10.1016/j.jcp.2014.01.011.
• (22) N. Srebro, T. Jaakkola, Weighted low--rank approximations, in: ICML 3 (2003) 720--727.
• (23) A. M. Buchanan, A. W. Fitzgibbon, Damped newton algorithms for matrix factorization with missing data, in: Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, IEEE 2 (2005) 316--322.
• (24) P. Chen, Optimization algorithms on subspaces: Revisiting missing data problem in low-rank matrix, International Journal of Computer Vision 80 (1) (2008) 125--142.
• (25) D. L. Donoho, Compressed sensing, Information Theory, IEEE Transactions on 52 (4) (2006) 1289--1306.
• (26) E. J. Candès, The restricted isometry property and its implications for compressed sensing, Comptes Rendus Mathematique 346 (9) (2008) 589--592.
• (27) E. J. Candès, B. Recht, Exact matrix completion via convex optimization, Foundations of Computational Mathematics 9 (6) (2009) 717--772.
• (28) J. F. Cai, E. J. Candès, Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM Journal on Optimization 20 (4) (2010) 1956--1982.
• (29) P. Jain, P. Netrapalli, S. Sanghavi, Low-rank matrix completion using alternating minimization, in: Proceedings of the 45th annual ACM symposium on Symposium on theory of computing, ACM (2013) 665--674.
• (30) B. Vandereycken, Low-rank matrix completion by riemannian optimization, SIAM Journal on Optimization 23 (2) (2013) 1214--1236.
• (31) T. Okatani, K. Deguchi, On the wiberg algorithm for matrix factorization in the presence of missing components, International Journal of Computer Vision 72 (3) (2007) 329--337.
• (32) K. Mitra, S. Sheorey, R. Chellappa, Large-scale matrix factorization with missing data under additional constraints, in: Advances in Neural Information Processing Systems (2010) 1651--1659.
• (33) I. Markovsky, Low rank approximation: algorithms, implementation, applications, communications and control engineering, Springer (2012).
• (34) I. Markovsky, K. Usevich, Software for weighted structured low-rank approximation, Journal of Computational and Applied Mathematics 256 (2013) 278--292.
• (35) X. Zeng, C. Farhat, A systematic approach for constructing higher-order immersed boundary and ghost fluid methods for fluid and fluid-structure interaction problems, Journal of Computational Physics 231 (2012) 2892--2923.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters