Low order dynamical system for fluid-rigid body interaction problem using POD method

ABSTRACT This paper describes the Reduced Order Modeling (ROM) for fluid rigid body interaction problem and discusses Proper Orthogonal Decomposition (POD) utilisation. The principal difficulty for using POD being the moving domains, a referenced fixed domain is introduced. The POD is applied for the velocity field obtained on the fixed domain. Then a method to reduce dynamical system in rigid body fluid interaction is developed. This method uses the fictitious domain method approach, which consists in treating the entire fluid-solid domain as a fluid. The rigid body is considered as a fluid, by using a high viscosity which can play the role of a penalisation factor of the rigidity constraint, and by adding a distributed Lagrange multiplied associated to this constraint in the weak formulation. The fluid flow problem is then formulated on the reference domain and POD modes are used in the weak formulation. The results are compared with computational solution and discussed.


INTRODUCTION
Classical problems in fluid structure interaction (FSI) require a long computational time.Indeed FSI simulation needs to couple fluid and structure solvers, using the smaller time step, and an important mesh near the fluid solid interface, which have to be remeshed with the interface movement.In active control context, where compute in real time is asked, current methods are not adapted.In a way to obtain short computational time, reduced-order modelisation with low order dynamical system has been developed.Different methods had been proposed for FSI context, the most significant are referred by Dowell [1].
There are two possible ways to construct ROM.The most famous one uses the notion of eigenmodes of the flow fluid.This approach characterises a field in terms of a relatively small number of global modes.By modes we mean a distribution of variables that characterises a gross motion of the studied physical system.There are several techniques to find these modes.One of these techniques concentrates on extracting the eigenmodes from the used model [2] (finite element for example).However, in case of a very high-dimensional system, extracting eigenmodes can be computationally very expensive.Thus, we use another method such as the method of balanced modes [3] or Proper Orthogonal Decomposition (POD) which will be explained in details in section 2. The second technique to determine ROM is the input/output model.This method uses a transfer function, that typically receives input in structure modes and gives as output generalised forces weighted by structural modes [4].
We have chosen to study POD capacities in fluid structure interaction.Indeed, this method was introduced in 1967 by Lumley [5] in fluid mechanics in order to extract coherent structures in a turbulent fluid flow.It has been intensively used since the 90's in many applications, such as flows in a driven cavity [6], in boundary layer [7], or particle dispersion [8].
In structure mechanics, POD is a recent investigation domain similar to modal analysis [9] in case of structure vibration.There are only few studies available in Fluid Structure Interaction.
First, this paper recalls the well-known POD method and its application in fluid mechanics.Next, the classical ALE formulation for Navier Stockes equations for fluid rigid bodies interaction is recalled.In the following section, a POD application method to a general case of moving boundary problems is presented.Then, a low order dynamical system in case of fluid-rigid body interaction problem is proposed and finally tested on differents cases.

THE PROPER ORTHOGONAL DECOMPOSITION (POD) 2.1. THE POD FORMULATION
In this section, the POD method is briefly introduced.A detailed methodology is already stipulated in literature [5,10], for bibliography see [11,12].
Considering a space Ω ʚ R n , n = 1, 2 or 3, T ʚ R, x ∈ Ω, t ∈T.The POD consists in finding a determinist function ψ, in a Hilbert space H, which gives the optimum representation of a velocity field v (x, t) ∈H (Ω, T ), by solving the following maximisation problem: (1) where denotes a statistic average operator, (•, •) denotes the inner product of H and •  2 H the associated norm.In the case of H = L 2 (Ω), the maximisation of problem (1) leads to solve the following eigenvalue problem: ( where R is the symmetric spatial correlation tensor, defined non-negative: ( Moreover, if R is continuous, the following operator is compact.Then the Hilbert-Schmidt theorem assures that there exists a set of positive eigenvalues (λ i ) i ≥1 which decrease to 0 λ λ λ λ and a set of eigenmodes (Φi) i ≥1 which is a Hilbertien basis for H. Thus, v can be decomposed according the eigenmodes as: (6) where a i are the temporal coefficients.(Φ i ) i ≥1 are named modes.

POD MODES PROPERTIES •
The spatial modes, (Φ i ) are orthogonal and can be normalised: (7) and they fullfill the boundary conditions.In case of an incompressible fluid, the velocity field fullfills the free divergence, ie div Φ i = 0.

•
The temporal coefficients a i (t) are obtained from the projection of v onto the (Φ i ) basis: • Moreover, they are not corellated and the eigenvalues are the temporal average: (without summation on the repeated indices) • The eigenvalue λ i is the energy captured by the mode Φ i .For a given N, the POD decomposition is the best energy decomposition which can be obtained.

THE SNAPSHOT POD
Solving equation ( 2) can be computationally intensive in higher dimensional problem.In order to minimise the computational times, the so-called snapshot method has been introduced by Sirovitch [13].
Let N m be the node number, n c the component number and Φ a mode POD.If sampling of M realisations (M << N m n c ) of the flow is sufficient to describe the problem, then we search the temporal coefficients a k such as: (10) Introducing the temporal average , and using the inner product of L 2 (Ω), we have to solve the following eigenvalue problem: (11) Hence, the spatial modes Φ i are obtained from equation (10) and the temporal coefficients a i are found by solving equation (8).
The use of the classical or snapshot method depends on the data type.In case of experimental data, the classical method is used.However, in case of computational simulation with a significant grid and time limit, the snapshot method is preferred.Note that when a non-stationary problem is considered, another solution would be the bi-orthogonal decomposition [14].

APPLICATION IN FLUID MECHANIC
In this section some general results obtained in the last decade in reduced order modeling in fluid mechanics are presented.
Considering an incompressible fluid in a rigid domain Ω f , with the density ρ f and the viscosity µ f , the fluid velocity v and the pressure p follow the dimensionless momentum equation of the flow: where R e is the Reynolds number.v is decomposed on the POD basis truncated at N modes: (13) The velocity decomposition (13) is introduced in the equation (12) and projected onto the POD vector Φ n .Thanks to the orthogonality and the free divergence of the POD basis Φ, the following system has been obtained: (14) where (15) (16) (17) (18) with n f being the outward normal on the domain Ω f considered of the boundary ∂Ω f .
The coefficient D n which contains the pressure term p can be avoided.Indeed, for more cases (for example the driven cavity) the velocity field is equal to zero on the boundary, and as the POD vector complies with homegeneous boundary conditions, D n is null.In other cases, some particular methods have been developed.We can mention Rempfer [15], who uses a vorticity formulation, Aubry [16] modelled this term for the study of the viscous sublayer in a channel flow and Allery [11] uses a penalisation method.This method will be explained later in section 6.
Thus, an N order dynamical system is obtained (equation 14).It consists in solving a ordinary differential equation system in time with the coefficients B, C, D, which are non-time dependent and computed only once.In the practice, N is in order ten, that is why this system is denoted low order dynamical system.Aubry [16] constructed the first model for studying the motion of structures in the wall region of a flat plate turbulent.Next, various configurations have been studied, for example near wall boundary layer ( [15], [17], channel flow [18], [19], driven cavity [6].Further examples deal with [12], [20].
One exemple of this method efficiency is given by Allery [21], [8], who applied a method to study the Coanda effect [21] and to model the fluid flow for computation of particle dispersion in a two-dimensional ventilated cavity [8].He underlines that for the Coanda effect only six modes are sufficient to completely capture the spatial structure of the flow and to obtain a good reconstruction with a low order dynamical system.In the second article, he indicates that only four modes are necessary to obtain a residual between the reconstructed velocity and the snapshot velocity less than 1.8 10 −2 .

ALE DESCRIPTION OF THE NAVIER-STOKES EQUATIONS
In this section, the ALE method for fluid rigid body interaction problems is exposed.This method is used to obtain the snapshots of the velocity field.
Consider a rigid body immersed in an incompressible fluid.Figure 1 shows a schematic description of the problem domain of interest, where Ω s (t) is the domain occupied by the moving rigid body, whose the center of gravity is denoted by G; Ω f (t) is the moving spatial domain upon which the fluid motion is described; Γ I (t) is the interface between Ω s (t) and Ω f (t), and n the outward normal.The interface Γ I (t) moves according to the position of the rigid body Ω s (t).
The motion of the fluid is governed by the incompressible Navier-Stokes equations which are given as follows in the ALE description [22,23,24]: where v is the fluid velocity, w the fluid mesh velocity, p is the fluid pressure, µ is the dynamic viscosity.The boundary, ∂Ω f (t)\Γ I (t), is divided into two parts on which the following boundary conditions are specified: (20) (21) where n r is the unit outward normal vector to Γ r and: (22) is the Cauchy stress tensor.
In this study Γ l and Γ r are assumed to be fixed in space.This assumption leads to (23) The velocity on the moving interface Γ I (t) is denoted by v I .This velocity is unknown, but: because of the non-slip condition on Γ I (t).In other words, are dealing with the Lagrangian description on Γ I (t).

THE RIGID BODY MOTION EQUATIONS
In our application we consider a planar motion of a system of rigid bodies.In this case the equation of motion is solved by the finite elements method and it can be written: (25) where u is the vector of structural displacement, and M, C and K are the mass, damping and stiffness matrices, F and b contain the fluid force F f and gravity force respectively.The fluid force is given by the following equation: (26)

MESH MOTION DESCRIPTION
The mesh velocity vector w may be arbitrarily specified though it has to satisfy the following boundary conditions, where v s is the velocity vector of the rigid body nodes on the interface.
In this study, the mesh velocity in the fluid domain is determined by solving the following equation (28) where λ(x) is judiciously chosen to control the mesh deformation.In this case, λ(x)is taken equal to 1 for all x ∈ Ω f (t).

FLUID-RIGID BODY INTERACTION ALGORITHM
Implicit Euler method is used to the time discretisation of fluid equations (19) and finite elements method to the spatial discretisation.
To solve the coupling equations, the following explicit scheme is used [25].Supposing that time t = t n , the fluid velocity and pressure fields, the rigid body displacement and position are known.The time step of the Navier-Stokes equations solver is the same as in the rigid body equations. 1.
The body equations (25) are solved in order to compute the displacement velocity v S at time t n+1 .Then, the position of The mesh velocity equation ( 28) is solved and the velocity w (n+1) of the fluid nodes displacement at time t n+1 is determined.

3.
The rigid body and the fluid nodes are moved at the predicted position by solving the equation x (n+1) = w (n+1) ᭝t + x n for all the mesh nodes.Then the fluid domain Ω f (t n+1 ) is defined.4.
The fluid equations (19) are solved in the domain Ω f (t n+1 ).

5.
The fluid forces acting on the rigid body are computed using the equations ( 26).This explicit algorithm is easy to implement, but it is only of order one and requires a small time step to its stability.An implicit scheme can also be used.

NUMERICAL APPLICATION 4.4.1. Rigid cylinder in an annular fluide space
The first studied case considers a two-dimensional rigid cylinder immerged in an annular cavity (Figure 2(a)).At the beginning, the fluid is at rest and the rigid cylinder is removed from its equilibrium position.Due to the spring effect the solid moves, starts to oscillate and generates a fluid flow.The fluid flow generates fluide forces, which damp the fluid oscillations.
The motion of the fluid is governed by the incompressible Navier Stokes equations in the ALE formulation and computed using Castem code (CEA, 2005) during 6.28 s using the followings: R The initial coordinates of the rigid body center are x 0 =(0, 005, 0), and "at rest" the length of spring is equal to 0.1 m.
For the spatial discretisation of the Navier Stokes equations the finite element Crouzeix-Raviart (Q 2 -P 1 ) has been used.For the velocity-pression coupling, a projection method has been applicated, and a SUPG method for the convection term stabilisation has been employed.
The rigid body displacement follows the following equation: where The damping parameter ξ can be computed on one pseudo-oscillation period, which allows to evaluate the numerical solution obtained by Castem.Indeed, where δ denotes the decrement logarithmic curve of rigid body displacement.On Figure 4.4.1, the analytical solution and that obtained by Castem show a good fit, which validates the solution.

Rigid body immerged in a channel
The second case studied is concerned with a cylindrical solid rigid immerged in a channel fluid flow.The rigid body is attached to the bottom wall by a spring.In this case the solid movement is imposed, i.e. there is no displacement along the x axes.For the fluid parameters, we consider the fluid density ρ f = 1000 kg⋅m −1 , the viscosity µ f = 0.001 kg/m⋅s, the inlet velocity where f denotes the vortex shedding frequency.According to Koopman [26], the cylinder oscillation frequency is chosed to be equal to vortex shedding frequency.In fact, at this frequency the lock-in phenomenon is achieved, i.e. the vortex shedding frequency is the same as when the rigid body is fixed.Thus, the solid displacement follows the equation: The results have been validated by compared lift coefficient frequency.The lift coefficient curve is shown on Figure 4.4.2.This adimensionless frequency is equal to the adimensionless vortex shedding frequency.This results validate the computational velocity field obtained.100 snapshots have been taken all the 241.68 s.
Figures 6 and 7 show differents snapshots taken on one rigid oscillation period.

POD APPLICATION TO MOVING BOUNDARY PROBLEMS
The solution obtained by this method is computed on a time-variant grid.POD method cannot be applied in this case: indeed a fixed domain is necessary to apply POD.For example, the snapshot method leads to solve eigenvalue problem (11), thus compute the L 2 inner product of velocity field at two different instants.But the fluid velocity fields (or structure velocity fields) are not define on the same domain for two different time steps.Snapshot method cannot be applied.One solution consists in considering a global domain Ω fixed in time which contains all the time variant domains and in searching POD vectors for the global velocity fluid on this domain.Thus in case of a fluid structure interaction problem the global fluid velocity is: where v f (respectively v s ) denotes the fluid velocity (respectively solid velocity), I the unit and I Ω f (x, t) the characteristic function of fluid domain is defined as the following: (34) The main objective is to obtain low order dynamical systems using POD basis truncated at N modes.In fluid mechanics, on a fixed domain, the Navier Stokes equations are projected on the POD vectors, and thanks to the orthonormality of POD modes, a low order dynamical system whose size is equal to the number of POD vector selected is obtained.For Fluid Structure Interaction, the ALE formulation of the Navier Stokes equations cannot be used, because POD basis have to be compute on a fixe domain.

LOW ORDER DYNAMICAL SYSTEM FOR FLUID RIGID BODY INTERACTION PROBLEMS
Next the low order dynamical system for fluid rigid body interaction problems is studied.This case is chosen due to previous works about eulerian formulation of rigid bodies [27,28].
This approach uses domain fictitious method, which consists in treating the entire fluidsolid rigid domain (the fictitious domain) as a fluid, by using Navier-Stokes equations for solid rigid domain and adding a rigid constraint to the solid domain: where ρ s is the solid density, D(v) = 0 denotes the rigid contrainst, v Γ the velocity fluid at fluid-structure interface, D(•) denotes the strain operator: x t 70 Low order dynamical system for fluid-rigid body interaction problem using POD method (35) Thus the Navier Stokes equations are extended to the rigid domain.To apply the rigidity contrainst a lagrange multiplier D(λ), which correspond to a stress field, is used [27].The weak formulation using POD basis leads to the following: Find a field v such as div v = 0 and for all Φ a virtual velocity field, div Φ = 0: where ρ f is the fluid density, µ f the fluid viscosity and the solid viscosity µ s is the penalisation factor of the rigidity constraint and Γ f = ∂Ω f \Γ I .
ρ and ν are defined on the global domain Ω: Thus the weak formulation is obtained for the global domain Ω with information about fluid and solid domains which are contained in density ρ and viscosity µ functions.
The basis obtained by solving ( 2) is truncated at N modes.For example N is searched as > 0.9999, where λ i denotes the i th eigenvalue of POD problem, and M the snapshot number.Thus the velocity field v is evaluated by using this truncated basis as equation ( 13): (38 This decomposition is introduced in (36) and the following dynamical system is obtained There are some differences with low order dynamical system obtained using POD basis in classic fluid mechanics exposed in section 3.In fact coefficients A, B, … are time dependent and must be computed at each time step.

The term
pΦ n ndx is avoided by a stress formulation in conjunction with a penalisation of the non-homogeneous Dirichlet boundary conditions [11].This method consists in considering that on the non zero velocity fluid boundary Γ σ , the strain boundary condition is the following: where u Γ σ denotes the velocity computed on Γ σ , u BC the velocity imposed and G a constant.Choosing G relatively larger than F  leads to consider that Finally, the initial problem is transformed into a more simple system of ordinary differential equation in a i (t) with few degrees of freedom.Indeed, in practice the POD method gives a basis which is maximal in terms of energy, with only a few functions.

APPLICATION 7.1. ONE DIMENSIONAL EXAMPLE
The first case study considers a one dimensional problem of fluid modelled with the Burgers equation coupled with spring mass equation (Figure 8).The Burgers equation is solved with the ALE method [22,23,24] with an initial fluid velocity v (x, 0) = sin (2πx) on Ω f (0), and the initial interface position Γ I (0) =1.
∫ Γf Low order dynamical system for fluid-rigid body interaction problem using POD method y 0 x

Energy contribution of POD modes
We apply POD method developped in Section 2.3.The first six POD modes have more than 99.99% of the kinetic energy (Table 1).The first mode gives information about the most energetical structure of fluid, and is constant on the solid domain.The followings modes capture fluctuations at the fluidstructure interface and have constant value on rigid domain.
Next, the reduced system using the 6 first modes is solved with N = 6.Each POD mode having a energy contribution, truncated the basis at the first N more significant modes, means fail the smaller energy scale, which are the most dissipatives.Thus, the low-order truncation of the POD basis inhibits generally all the transfers between the large and the small (unresolved) scales of the fluid flow.Consequently, to recover the effects of the truncated modes, that is generally of the small scales, we use an "eddy viscosity" [16].The viscosity for the mode i is multiplied by (1 + i0.001).
Reducing the system with only six modes, gives a good result with an error less 4%, illustrated in Figure 9(a).The figure 9(b) shows that the error is relatively not very significant.
The reconstructed temporal modes resulting from the reduced system and those resulting from the snapshot method are compared on figure 10.The difference is very small and comes from computational approximation of modes and derivatives at the interface.Indeed, the   interface position is between nodes and we need to know for example Φ i (Γ I ) by interpolating on the first two nearest modes.
The first mode which is consistent with the mean velocities properties, captures at least 99% of the kinetic energy in all tested cases.Other modes are used to keep velocity variations in the fluid.We can see that where the domain is always structured, the POD modes are constant.In fact we supposed that the body is a rigid structure.First, the POD vectors have been searched and the reconstructed velocity has been defined as and is compared to the initial velocity.On Figure 12, we can see that we have a good reconstruction with 3 modes and a maximal error near to the fluid-solid interface.
In fact the first POD eigenvalue contains 99.2% of the total kinetic energy and with three vectors almost 99.99% of the energy is captured.That is why with only three POD vectors the reconstructed velocity is a good approximation of the initial velocity.
Next, the low order dynamical system with three modes is constructed and the temporal coefficients obtained are compared with those obtained by computing the POD vector (at each snapshot .There is a good conformity between them, for example for the first temporal coefficient a 1 (Figure 13(a)).
In the following the low order dynamical system during a longer period than the snapshot period will be presented.The solution obtained has not been compared to a numerical solution, but the gravity center displacement can be predicted by an analytical solution and compared to those computed.In fact, in this case the analytical solution is being evaluated on the snapshot period.The solution for a period longer than 10 times the snapshot period gives good prediction for the gravity center position (Figure 13(b)).It is an adequate criteria to conclude that the reduced system obtained gives a good result with a few degrees of freedom and this case gives a good prediction for simulate a period longer than the snapshot period.

EXAMPLE 2
In what follows, the ROM method is applied to the second case presented in section 4.4.2.First, the POD vectors have been computed, and the results indicate that almost ten modes are needed to obtain almost total energy.Results are identical in comparison with the reconstructed velocity on the truncated basis (Figure 14).
Figure 15 shows the two isovalues components of the two first modes.We can see that the first vector contains the most important information about fluid flow compared to velocity field at different time steps (Figure 6).
The first POD vector does not capture very well the fluid fluctuation near the moving interface.Those are the following POD vectors which contain this information.That is why, although Reynolds number is very low (200), more POD vectors are needed compared to the POD application in fluid mechanics (see section 3).
The low order dynamical system using the first ten modes is constructed and computed with a good adequation to the initial velocity field.Indeed the maximal error in time is about 5% in L 2 norm.The difference between the isovalues of the reconstructed velocity field using reduced system solution and snapshot velocity field at the maximal error time are showed on Figure 16.The error is located on a very small area and does not disturb the rest of the flow.The result presented makes it possible to obtain the same Strouhal number.

CONCLUSION
This paper introduces an application method of the proper orthogonal decomposition for moving boundary problems, more specifically in the case of fluid rigid body interaction problem.The method uses a global domain fixed in time for computing the POD vectors of a global velocity fluid.To construct low order dynamical system for a fluid rigid body interaction problem, a fictitious domain method has been used and a variational formulation using POD bases is obtained.The method gives good results, but precision seems to be lost at the fluid solid interface.This method can be applied to 3D problems without difficulties.In the case of a deformable structure the current formulation cannot represent the deformations of the solid domain.A tensor deformation for the solid domain has to be introduced.It will be studied in a further work.

Figure 2
Figure 2 Schematic description of the problem domain.

Figure 3
Figure 3 Solid gravity center displacement.

Figure 4
Figure 4 Schematic description of problem.

Figure 6
Figure 6 Velocity field snapshots at different time.

Figure 8
Figure 8 Scheme of the fluid and structure domain.

Figure 9
Figure 9 Reconstruction error evaluation.

Figure 12 20 −
Figure 12  Isovalues of the difference between the reconstructed solution with 3 modes and the initial at the 70 snapshot.

Figure 13
Figure 13  Low order dynamical system result.

Figure 14 Figure 15
Figure 14 Reconstructed error in L 2 norm.

Table 1
Kinetic energy contribution for the five first modes.