A systematic approach for constructing higher-order immersed boundary and ghost fluid methods for fluid–structure interaction problems

https://doi.org/10.1016/j.jcp.2011.12.027Get rights and content

Abstract

A systematic approach is presented for constructing higher-order immersed boundary and ghost fluid methods for CFD in general, and fluid–structure interaction problems in particular. Such methods are gaining popularity because they simplify a number of computational issues. These range from gridding the fluid domain, to designing and implementing Eulerian-based algorithms for challenging fluid–structure applications characterized by large structural motions and deformations or topological changes. However, because they typically operate on non body-fitted grids, immersed boundary and ghost fluid methods also complicate other issues such as the treatment of wall boundary conditions in general, and fluid–structure transmission conditions in particular. These methods also tend to be at best first-order space-accurate at the immersed interfaces. In some cases, they are also provably inconsistent at these locations. A methodology is presented in this paper for addressing this issue. It is developed for inviscid flows and prescribed structural motions. For the sake of clarity, but without any loss of generality, this methodology is described in one and two dimensions. However, its extensions to flow-induced structural motions and three dimensions are straightforward. The proposed methodology leads to a departure from the current practice of populating ghost fluid values independently from the chosen spatial discretization scheme. Instead, it accounts for the pattern and properties of a preferred higher-order discretization scheme, and attributes ghost values as to preserve the formal order of spatial accuracy of this scheme. It is illustrated in this paper by its application to various finite difference and finite volume methods. Its impact is also demonstrated by one- and two-dimensional numerical experiments that confirm its theoretically proven ability to preserve higher-order spatial accuracy, including in the vicinity of the immersed interfaces.

Introduction

Eulerian immersed boundary and ghost fluid methods (for example, see [1], [2], [3], [4] and references cited therein) are usually preferred over alternative approaches based, for example, on the Arbitrary Lagrangian Eulerian framework [5], [6], for the solution of fluid–structure interaction (FSI) problems characterized by large structural motions and/or deformations or by topological changes. Such methods allow a rigid or flexible moving body to penetrate the CFD (computational fluid dynamics) grid. For this reason, they operate on non body-fitted grids. They typically distinguish between “real” (or “active”) fluid grid points lying in the physical fluid domain, and “ghost” (or “inactive” or “fictitious”) fluid grid points lying inside the moving obstacle. This distinction is usually a dynamic one because FSI problems are typically dynamic problems. When at the end of a computational time-instance tn the status of a fluid grid point changes from real to ghost — which is referred to in this paper as the RTG scenario — the value of the fluid state vector at this ghost point must be provided in order to enable the advancement of the flow computation to the next time-instance tn+1. Similarly, when at the end of tn the status of a fluid grid point changes from ghost to real — which is referred to here as the GTR scenario — a new value of the fluid state vector must be provided at this grid point.

For the RTG scenario, most approaches proposed in the literature fall into two categories. The first one consists of methods which project a ghost fluid grid point onto the fluid–structure interface (see Fig. 1.1), interpolate the velocity of the structure at the projection point, and set the value of the fluid velocity vector at the ghost fluid grid point (or the normal component of this vector when the flow is assumed to be inviscid) to this interpolated value (or its normal component for inviscid flows) [7], [2]. This enforces an approximate form of the non-penetration condition at the fluid–structure interface. In this case, the remaining primal variables of the fluid state vector at the ghost fluid grid point (density, pressure, or entropy) are computed by first approximating their counterparts at the projection point using data from nearby real fluid grid points and either interpolation or extrapolation, then extrapolating the obtained values to the location of that ghost fluid grid point. Hence, this category of methods is easy to implement. However, it is only first-order accurate in space. For this reason, it is often equipped with adaptive mesh refinement (AMR) [8], [2] in the vicinity of the fluid–structure interface. This improves its practical accuracy, but at the cost of an increase of its implementational (and computational) complexity, especially in three dimensions (3D). The methods in the second category [4], [9] can be labeled as “mirroring” methods. Not only they project a ghost fluid grid point of interest onto the fluid–structure interface, but they also determine its reflection with respect to the tangential interface — a line in two dimensions (2D) and a plane in three dimensions (3D) — containing the projection point. Since a mirroring point usually falls in the physical fluid domain, the fluid state variables at this point are computed by interpolation. Then, the fluid state vector at the ghost fluid grid point is computed as follows. The fluid velocity vector (or its normal component in the case of an inviscid flow) is obtained by linearly extrapolating the interpolated fluid velocity vector at the mirroring point and the structural velocity vector at the projection point (see Fig. 1.2). The values of the remaining primal fluid variables at the ghost fluid grid point are set to the interpolated values at the mirroring real fluid grid point. In theory, the spatial accuracy of this second category of methods is one order higher than that of the first one described above. However, to the best of knowledge of the authors, second-order spatial accuracy has been reported for such methods only for the case of a stationary fluid–structure interface [9]. For genuine FSI problems with dynamic fluid–structure interfaces, the mirroring procedure is sometimes combined with local grid refinement (for example, see [4]) to enhance the delivered spatial accuracy in the vicinity of the fluid–structure interface.

For the GTR scenario, two different computational strategies can be found in the literature. In the first one, the fluid state vector at a fluid grid point which switches from the ghost to the real status at the end of time-instance tn is determined from the knowledge, at the end of the computational step performed at tn, of its counterparts at nearby real fluid grid points. For example, such a fluid state vector is set in [4] to that of the nearest real fluid point at the end of tn. This approach introduces a spatial discretization error that is proportional to the local grid size. Higher-order spatial approximations of that fluid state vector require data from multiple nearby real fluid grid points and raise the non-trivial issue of the selection of these points, particularly in 3D. In the second strategy which specifically characterizes the ghost fluid method (GFM) [10], [1], all ghost fluid grid points with the potential to switch status at the end of the computational step performed at tn are updated during this computational step as if they were real fluid grid points; consequently, at the end of this computational step, the fluid state vectors at the subset of these grid points which turn into real fluid grid points are identified with their updated ghost values. This approach was successfully incorporated into various computational methods to achieve second- or higher-order spatial accuracy in the solution of fluid–fluid interaction problems with moving material interfaces [10], [11], [12]. For this reason, and because it is simpler than the alternatives, this approach is chosen in this work as a foundation for treating the GTR scenario.

Most approaches overviewed above for populating the fluid state vectors at the grid points that switch status at the end of a computational step are practical and have been successfully applied to the solution of various types of FSI problems. However, many of the fluid–structure computational strategies based on these approaches have also relied on adaptive grid refinement near the fluid–structure interface to achieve higher accuracy, rather than on imposing higher-order discrete transmission conditions. Alternatively, this paper presents a systematic approach for constructing higher-order immersed boundary and ghost fluid methods for FSI problems. This approach equally applies to CFD problems with fixed wall boundaries.

Unlike the projection and mirroring techniques briefly described above, the approach proposed in this paper for populating the state vectors at the ghost fluid grid points takes into account the underlying CFD scheme and focuses directly on preserving its formal order of spatial accuracy. It distinguishes itself from other approaches such as that presented in [13] which is tightly coupled to the fractional step method and incompressible Navier–Stokes problems, by the fact that it is more than a mere higher-order scheme: it is a systematic procedure for constructing higher-order immersed boundary and ghost fluid methods given in principle any preferred CFD scheme for solving a problem of interest. Furthermore, the proposed approach treats both RTG and GTR scenarios in the same manner. More specifically, it constructs at each time-step the state vectors at the ghost fluid grid points so that the extension of the CFD scheme of interest into the ghost fluid domain “matches” an accuracy-preserving modification of this scheme in the neighborhood of the fluid–structure interface that eliminates the need for ghost fluid grid points. This operator matching method is developed here in one and two dimensions, assuming that the flow is inviscid and the dynamic motion of the structure is prescribed. However, its extensions to three dimensions and flow-induced structural motions are straightforward. Its extension to viscous flows follows the same principles, but involves numerous additional technical details; therefore, it will be presented elsewhere.

To this effect, the remainder of this paper is organized as follows. In Section 2, it is shown that for a one-dimensional (1D) advection problem discretized by the Lax–Wendroff [14] scheme, and for the 1D Euler equations discretized in space by the Roe scheme [15] and in time by forward differencing, the basic ghost fluid method for FSI (GFM-FSI) [16], [1] is simply inconsistent at the fluid–structure interface. Since this method relies on constant extrapolation to populate some ghost values, the results presented in Section 2 highlight the importance of the issues raised in this paper. In Section 3, a systematic approach is presented for populating the state vectors at ghost fluid grid points while preserving the formal order of spatial accuracy of the semi-discretization scheme underlying a given immersed boundary or ghost fluid method for CFD. In Sections 4 Application to the second-order Lax–Wendroff scheme, 5 Application to a third-order finite difference scheme, this approach is illustrated with its application to sample second- and third-order finite difference schemes, respectively. In Section 6, its application to sample finite volume schemes is described. In Section 7, the performance of the resulting schemes is assessed. Finally, conclusions are offered in Section 8.

Section snippets

Loss of consistency at a moving wall boundary: case studies

Here and throughout the remainder of this paper, the following notation is used:

  • (x,t) designates the exact value of a variable † at the location x (in Cartesian coordinates) and time t.

  • In contrast, x̃ and ỹ designate the Cartesian coordinates of a point at a fluid–structure interface.

  • The superscript n designates the time-instance tn and the superscript G indicates a ghost value.

  • The subscript i and pair of subscripts i, j designate a grid point in 1D and 2D, respectively.

  • The subscripts t, x

Systematic approach for populating ghost values and achieving higher-order spatial accuracy

The main idea behind the approach proposed in this paper for populating ghost values in a CFD computation with dynamic material interfaces is based on the observation made below.

Consider the boundary value problem (BVP)P(u)=0inΩu=gonΓwhere Ω denotes here a spatial domain of interest. For simplicity, assume that Γ = ∂Ω. Next, suppose that only the values of u on a subdomain Ω  Ω (see Fig. 3.1) are of interest, and that g=u|Γ where Γ = ∂Ω can be somehow obtained. In order to determine the

Application to the second-order Lax–Wendroff scheme

Again, consider the 1D advection equationut+aux=0The internal operator Lio associated with the application of the second-order Lax–Wendroff scheme to the solution of this equation is given byLio(un)=uin-12λui+1n-ui-1n+12λ2ui+1n-2uin+ui-1ni<In

To construct the associated boundary operator Lib, it is recalled that it suffices to achieve first-order spatial accuracy at the interface in order to achieve second-order spatial accuracy globally. Therefore, Lib is chosen here to be the discrete operator

Application to a third-order finite difference scheme

Consider again the 1D advection Eq. (4.1), and its discretization by the following method which combines a third-order upwind-biased finite difference semi-discretization and the third-order TVD Runge–Kutta [23] time-integratorui(1)=uin-16λ-ui+2n+6ui+1n-3uin-2ui-1niIn-2uIn(1)=uInn-12λ3uInn-4uIn-1n+uIn-2nuIn-1(1)=uIn-1n-12λuInn-uIn-2nui(2)=34uin+14ui(1)-16λ-ui+2(1)+6ui+1(1)-3ui(1)-2ui-1(1)iIn-2uIn(2)=34uInn+14uIn(1)-12λ3uIn(1)-4uIn-1(1)+uIn-2(1)uIn-1(2)=34uIn-1n+14uIn-1(1)-12λuIn(1)-uIn-2(1)uin

Application to second-order finite volume schemes

Consider again the 1D Euler Eq. (2.4) and their discretization by a method combining a finite volume approximation based on the Roe flux [15] and MUSCL technique [17] (Roe-MUSCL), and the second-order TVD Runge–Kutta [23] time-integrator. A ghost-free version of this method can be described as followswi(1)=win-ΔtΔxFi+1/2Roe,n-Fi-1/2Roe,niIn-1wIn(1)=wInn-ΔtΔxf(w(tn))-FIn-1/2Roe,nFj+1/2Roe,n=FRoewj+1/2,ln,wj+1/2,rnjIn-1wk+1/2,ln=wkn+12ΔxσknkInwk+1/2,rn=wk+1n-12Δxσk+1nkIn-1wi(2)=wi(1)-ΔtΔxFi+1

Numerical examples and performance assessment

Here, the proposed systematic approach for enhancing the formal order of spatial accuracy of ghost fluid and immersed boundary methods is illustrated, and its performance is demonstrated with the solution of various 1D and 2D transport problems with prescribed moving wall boundary conditions. The errors are computed as the differences between the obtained numerical results and analytical results when these are available, or the obtained numerical solutions and carefully computed reference

Conclusions

A systematic approach for constructing higher-order immersed boundary and ghost fluid methods for CFD in general, and fluid–structure interaction problems in particular is presented in this paper. This method leads to a departure from the current practice of populating ghost fluid values independently from the chosen spatial discretization scheme. At each time-step, it evaluates the ghost data so that the extension of a CFD scheme of interest into the ghost fluid domain “matches” an

Acknowledgments

The authors acknowledge partial support by the Office of Naval Research under Grants N00014-06-1-0505 and N00014-09-C-015, and partial support by the Army Research Laboratory through the Army High Performance Computing Research Center under Cooperative Agreement W911NF-07-2-0027. The first author also acknowledges partial support by a Stanford Graduate Fellowship.

References (24)

Cited by (0)

View full text