A non-iterative direct forcing immersed boundary method for strongly-coupled fluid–solid interactions

Abstract

A non-iterative direct forcing immersed boundary method is presented for the strongly-coupled simulations of fluid–solid interactions. While it retains many advantages of the immersed boundary framework by Yang and Stern (2012) [30], especially the simplified field extension strategy for moving boundary treatment and the pointwise integration of hydrodynamic force using the momentum forcing term, the present approach improves upon the previous method in several aspects including optimized computational cost for a strong coupling scheme, reduced algorithm complexity for a straightforward implementation, and enhanced numerical stability for low density ratio problems. Central to these improvements is a simple intermediate step in which the velocity fields around solid bodies are predicted on temporary non-inertial reference frames attached to moving solid bodies in a one-to-one manner. This step enables the explicit inverse of the implicit equations for rigid body dynamics, thus rendering unnecessary the previous predictor–corrector scheme for iteratively adjusting the displacements and velocities of the immersed bodies until reaching a convergence. In addition, a simple, generalized procedure is developed to obtain the interpolation coefficients in a local reconstruction stencil explicitly from the geometric relationship. For verification and validation, the vortex-induced vibration of a circular cylinder and the rotational galloping of a rectangular body are considered first; then several particulate flow problems, including settling and buoyant particles of low density ratios, a settling particle in a small container, and the kissing–drafting–tumbling problem of two settling particles, are studied. The agreement between the present results and the reference data in the literature is excellent. An overall second-order accuracy of the algorithm is verified in two systematic grid convergence tests. The present idea can be easily applied to similar methods for achieving a strong coupling scheme on top of a weak one with a nominal increase in the computational cost. Details of the algorithm are provided to facilitate its implementation in other solvers using non-boundary-conforming grids.

Introduction

The immersed boundary method has long been established as a versatile and cost-effective approach for fluid–structure interaction (FSI) problems since its introduction by Peskin in the 1970s [17], [18]. It is very attractive in that a simple forcing term is added to the momentum equation to represent the effect of a complex immersed boundary on the fluid flow, without modifying the regular finite difference schemes on a fixed Cartesian grid. In the conventional immersed boundary method [19], a constitutive relation for the material elasticity is required to determine the forcing term according to the structure deformations. Structures with no/prescribed deformations can be modeled using very stiff springs; but large stiffness values, which are usually determined ad hoc, have some numerical stability implications. In the so-called direct forcing approach [15], instead of making use of a constitutive relation, the momentum forcing function was derived by directly injecting the velocity boundary condition at the immersed boundary into the discrete-time momentum equation. Essentially it can be viewed as a local solution reconstruction procedure with which the desired boundary conditions is imposed on the immersed boundary, and the explicit appearance of a forcing term in the momentum equation was not even required [6]. This opened up great opportunities for developing algorithms of various features targeted for specific improvements as well as diversified applications (see, for example, [12], [23], [1] among others), since the velocity boundary conditions are ubiquitous in numerical simulations of fluid flows regardless of the boundary properties (rigid, deforming, or elastic) and the motion characteristics (stationary, prescribed, or predicted).

It should be noted that, in Peskin's method, the momentum forcing from an immersed boundary represented by a Lagrangian mesh is formulated as a Dirac delta function, which is usually replaced by a smooth approximation in the actual discretization on the underlying Eulerian grid. For a moving boundary problem, ideally, a discrete delta function can smoothly spread out the boundary force onto the surrounding Eulerian grid points without incurring spurious jumps in the flow field, although the smoothing process results in a blurring fluid–solid interface. On the contrary, in a direct forcing approach [6], the velocity boundary conditions are imposed on grid points adjoining the immersed boundary and satisfied “exactly” as for a boundary-conforming method if the grid points coincide with the immersed boundary. This sharp interface treatment produces temporal and spatial jumps in the momentum forcing when the immersed boundary moves across the fixed grid points. Without an appropriate handling, these jumps can result in spurious oscillations of hydrodynamic forces on the immersed boundary, which could be disastrous for FSI problems sensitive to the variations in the hydrodynamic forces. Actually, this issue was a main motivation in [24] to resort to discrete delta functions for the smooth transfer of velocity and forcing information between the Eulerian and Lagrangian locations. Similar Eulerian–Lagrangian coupling ideas in the direct forcing framework were also presented in [8], [34]. Uhlmann's method has been widely used in particulate flow simulations and inspired quite some follow-up work for improvements or extensions to general FSI problems, such as [25], [33], [20], [10] among others. In particular, a major limitation of [24] restricting the applicable range of density ratios (${\rho }_s/{\rho }_f⩾1.2$ for spherical particles with ${\rho }_s$ and ${\rho }_f$ the solid and fluid densities, respectively) was mitigated to cover cases with ${\rho }_s/{\rho }_f⩾0.3$ in [10]. Unfortunately, most if not all of these improvements or extensions came at the price of additional algorithm complexity and increased computational cost, and hardly addressed several inherent drawbacks of the discrete delta function formulations. For instance, depending on the width of the support of the delta function on the Eulerian grid, the sharpness of the fluid–solid interface is weakened and the computational cost of the Eulerian–Lagrangian coupling process is augmented correspondingly; also, some ad hoc treatments may be necessary if the Eulerian stencils from two different structures overlap [10]. Furthermore, the Lagrangian markers on the immersed boundary are usually required to be distributed evenly with a resolution close to that of the local Eulerian grid (a finer surface mesh will increase the cost, instead of the accuracy) [24], even for stationary geometries or structures with prescribed motions. This could be very difficult, if not impossible, for cases with complex geometries or thousands/millions of particles of arbitrary shapes.

Due to the abovementioned drawbacks, it is still preferred to avoid the dependence on a discrete delta function as the kernel of fluid–structure (i.e., Eulerian–Lagrangian) coupling in many applications by retaining the sharp interface property of the original direct forcing approach [6]. It should be pointed out that the validation in [6] was mainly concerned with stationary immersed boundaries, although a moving boundary problem was demonstrated without force information. Thus, the implications of an immersed boundary moving on a fixed grid in the framework of a fractional-step method were not addressed in [6]. Actually, a close observation in [26] on the role switching of the grid points around the immersed boundary revealed that, when a grid point with a reconstructed solution at the previous time step becomes a normal fluid point at which the governing equations are solved, the non-physical field information from the solid phase may enter the system via the derivatives evaluated at this point. Depending on the deviation from the physical values, these contaminated derivatives may produce spurious force oscillations of various amplitudes when the original method in [6] is applied to moving boundary problems. Therefore, an intuitive remedy of this problem is to exclude the involvement of those points from the solid phase in evaluating the derivatives at such a point. But the penalty will be the loss of the regularity of the discretization stencils for individualizing the treatments in a case-by-case manner. On the contrary, a field extension strategy was proposed in [26] to recover the physical values of the contaminated derivatives by extending the flow field into the solid phase through an extrapolation procedure. The field extension can be performed at the end (or equally, the beginning) of each time step by directly modifying the solution field without tangling with the discretization stencils. Several moving boundary problems ranging from laminar to turbulent flows were carried out to demonstrate the accuracy of this simple approach in [26]. The extension to FSI problems with multiple rigid bodies using a strong coupling predictor–corrector scheme in [27] further verified the effectiveness of the field extension strategy in tackling the issue of spurious force oscillations.

An important feature in [26], [27] is that the hydrodynamic force was calculated through a surface integration on all Lagrangian elements; thus there was still a resolution requirement of the surface mesh, e.g., close to that of the local Eulerian grid, in order to obtain an accurate force evaluation. Recently, Yang and Stern [30] substantially simplified the field extension procedure from [26] and significantly accelerated the predictor–corrector algorithm for fluid–structure coupling from [27]. The setup of extrapolation stencils for the field extension was eliminated to minimize the geometric operations involved in the immersed boundary setup procedure and the fluid solver was moved out of the predictor–corrector iterative loop. The latter was facilitated by a switch from the surface integration to a point integration for the hydrodynamic force evaluation. In the new method the Lagrangian surface discretization is employed for the sole purpose of adequately representing the surface geometries and is not revelant to the accuracy of the force integration any more. Actually, for some simple geometries with analytical representations, such as a circular cylinder or a sphere, a surface discretization is not even required [31]. The new method represents an acceleration of several times or up to one order of magnitude compared with [27], thanks to the fact that the pressure Poisson equation is solved only once in each time step instead of multiple times in [27].

On the other hand, the fluid–structure coupling algorithm in [30] still contains an iterative loop, in which the displacements and velocities of the immersed structures are adjusted and the momentum forcing term is updated until reaching a convergence. For an iterative scheme, the major concerns are usually the numerical stability, computational expense, and algorithm complexity as well. The work in [31] did indicate that the numerical stability property of the predictor–corrector scheme is only slightly better than that of [10] for low density ratio problems. Moreover, the repeated immersed boundary setup within the loop could be very expensive for complex geometries with millions of Lagrangian surface elements even for a robust and efficient procedure like the one presented in [32]. In addition, the predictor–corrector scheme is inevitably quite involving in its implementation for executing different components during different stages of the solution process.

It should be noted that, in general, a weak coupling scheme is optimal in terms of algorithm complexity and computational expense for its non-iterative nature. Unfortunately, the numerical stability of a weak coupling scheme is not satisfactory in many applications. Therefore, strong coupling schemes are usually preferred in FSI problems. Actually, a non-iterative strong coupling scheme was developed for a single rigid body interacting with an incompressible viscous flow in [11] within the framework of a direct forcing immersed boundary method. The salient feature in this approach was the adoption of a conservative formulation of the Navier–Stokes equations in a non-inertial reference frame fixed to the body. Consequently there is no relative motion between the body and the grid, and the immersed boundary only needs to be set up once at the beginning of the simulation. More importantly, with an explicit, linear momentum forcing, a simple linear relationship between the fluid forces/moments and the body velocities was identified in [11], which enables a strongly coupled, direct update of the body velocities without any iterations. In addition, the issue of spurious force oscillations due to a moving immersed boundary on a fixed grid does not exist anymore. However, this approach is limited to the cases of a single rigid body within an infinite domain. If a large-amplitude angular motion is involved, the boundary condition treatment, especially, for the inflow and outflow, could be quite complicated; also, the time step could be excessively small due to the large “grid velocity” from the body rotation near the domain boundary. These limitations are inherent to a method formulated in the non-inertial frame; thus the wider application of the approach in [11] was severely restricted.

In this study, a novel scheme is proposed for the strongly-coupled simulations of fluid–solid interactions on the basis of the direct forcing immersed boundary framework in [30]. The predictor–corrector algorithm in [30] is replaced by a non-iterative scheme similar to the one given in [11]. Therefore, every component in the algorithm is executed only once in each time step; for complex geometries, this represents a huge saving in the immersed boundary setup compared to the predictor–corrector algorithm. Consequently, the design and implementation of the new scheme is considerably simpler with a much reduced algorithm complexity. Along with the simplicity, interestingly, comes an enhanced numerical stability property for low density ratio problems as to be shown later. On the other hand, similar to [30], our new scheme does not have any limitations in the nature of the computational domain and the number of rigid bodies in the domain. In addition, there does not exist any difference in the boundary condition treatments or any extra restriction on the time step. The key component of the new scheme is a simple intermediate step in which the velocity fields around solid bodies are predicted on temporary non-inertial frames. These non-inertial frames are parallel to the inertial frame that the full governing equations of fluid flow is solved on, and are temporarily attached to the moving solid bodies in a one-to-one manner. On a temporary non-inertial frame, the treatment in our algorithm can be regarded as a reformulation of the scheme in [11], thus it enables the direct inversion of the implicit equations for rigid body dynamics as well. After the velocities and hence displacements of all solid bodies are obtained, the full system is still solved in the inertial frame as in [30]. Therefore, those abovementioned limitations of the scheme in [11] can be totally avoided.

Moreover, as a supplementary improvement, an explicit linear relationship between the immersed boundary velocity and the reconstructed flow velocity near the body is developed for the new scheme. In [26] a linear relationship was applied in the local solution reconstruction, but it was implicitly expressed with a small linear system, which is not very convenient and flexible. In this paper, the construction procedure of the interpolation stencils is reformulated to obtain a simple explicit expression and retain the desirable normal-oriented property and the grid point configuration in [26]. Also, a generalized approach is designed to obtain interpolation stencils of various dimensions.