An immersed boundary method for incompressible flows in complex domains
Introduction
The presence of moving fluid–solid interfaces in a fluid flow is a common feature in many natural phenomena and engineering processes, such as the flows within mechanical pumps and cardiovascular systems, and those around aeroplanes, fish and birds. To approximate the moving fluid–solid interface in computational fluid dynamics simulations, much research in recent years has focused on the immersed boundary method (IBM), which has been applied for diverse classes and scales of problems e.g. flow-structure interactions around a heart valve [1], flows inside the cylinder of an IC engine [2], particulate flows, where the flow interactions with the surface of the particle are resolved [3], and flows around marine and aerial organisms [4]. Using an appropriate forcing scheme, the IBM enforces the no-slip condition on the immersed fluid–solid interface that is located inside a fixed domain, and the Lagrangian and Eulerian frameworks are used to model the interface and the flow domain, respectively. This circumvents the necessity to modify or reconstruct a boundary-conformal Eulerian mesh whenever the fluid domain deforms following the moving interface. Hence the implementation of a simulation framework using the IBM is substantially simpler and computationally more efficient than remeshing every time-step.
IBMs can generally be classified into two types: smooth- and sharp-interface types. The IBM was originally proposed as a smooth-interface type by Peskin [1]. In Peskin's work the solid–fluid interface, typically named the immersed boundary (IB), is represented by a set of Lagrangian markers which are spaced closer to each other than the size of the local Eulerian mesh cell. Fundamental steps of Peskin's IBM are the interpolation of flow velocities from the fluid nodes onto the marker points on the IB, the evaluation of the forces at the marker points, followed by the spreading of the forces back onto the fluid nodes, where the forces are treated as source terms in the momentum equations. Note that the term “spreading” is used interchangeably with the term “convolution” [5] and “regularisation” [6] in the literature. The terms “fluid node” and “IB marker point” are also used interchangeably with the terms “Eulerian point” and “Lagrangian point”, respectively. Both interpolation and spreading operations often employ distributed weights which are calculated from a smoothed approximation of the Dirac delta function, typically named the discrete delta function (DDF), centred at the Lagrangian points. The discrete delta functions that are typically used are derived for uniform Cartesian meshes, and their construction ensures the satisfaction of conditions, such as those presented by Peskin [7], for the continuity and compactness of the DDF, its conversation properties and translational invariance. Frequently used discrete delta functions are the compact 3-point [8] and 4-point [7] regularised delta functions. Using these distributed weights, the spreading of the forces at the Lagrangian points along the IB towards more than one layer of neighbouring Eulerian points results in a smeared representation of the boundary, hence the “smooth-interface” terminology. For the approximation of a boundary with zero thickness, the smooth-interface IBM is typically first-order accurate due to the jump of the normal derivative across the boundary [7], which occurs in typical applications of the IBM.
In contrast to the smooth-interface approach, a sharp-interface formulation imposes a velocity corresponding to the no-slip condition at the Lagrangian points, by constraining the nodal variables of the nearest neighbouring Eulerian cells using an interpolation function, typically a linear one [9], [10], [11], [4]. This results in the sharp representation of the IB, which can show a second-order spatial convergence behaviour on some specific features of the flow. Although this is favourable for some applications, an important drawback of the sharp-interface formulation is that with transient flows, the behaviour can become unphysical with spurious flow patterns for cases with moving boundaries, especially when a coarse mesh and/or a small time-step are used, due to geometric inconsistencies in the discretisation of the continuity constraint [12], and to the sudden transfer of cells from the solid to fluid phase (and vice-versa) between consecutive time-steps [13]. Yang and Balaras [14] have proposed a robust field-extension procedure that extends the pressure and velocity fields into the solid body in order to implicitly treat problematic grid cells, and which is shown to significantly reduce the occurrence of spurious oscillation. Seo [12] have also proposed an algorithm to substantially reduce the oscillations by adopting a conservative cut-cell approach, which complexify the implementation of the method. In addition to their numerical complexity, sharp-interface formulations require extra work and approximations in order to integrate the hydrodynamic forces on the IB.
The original smooth-interface scheme, as put forward by Peskin [1], has been developed for an IB that represents a thin elastic fibre, regarded as a soft solid structure. As outlined by Peskin [7], the interacting force field between the fluid and the solid phase is evaluated from a constitutive equation based on the structural deformation of the fibre. The calculated boundary force based on the deformation is then imposed in the solution for the following iteration or the new time level. Such a forcing scheme can be extended to approximate rigid interfaces by applying very large values of fibre stiffness [15], [16]. However, this results in a stiff system of equations to be solved, leading to a high computational cost. In order to approximate rigid boundaries, Goldstein et al. [17] proposed a feedback-forcing scheme that adopts principles from control theory, where the forcing for a new time-step is determined based on the differences between the interpolated fluid velocities at the previous time-step and the desired velocities at the Lagrangian points. The scheme involves a number of ad hoc parameters which need to be carefully selected in order to ensure that the no-slip condition is enforced satisfactorily at every time level, which often requires very small time-steps [16], [18], [19].
A direct evaluation of the forces required to enforce the IB, typically characterised as a direct-forcing method, was considered by Mohd-Yusof [20] in their pseudo-spectral code to approximate rigid boundaries. In their work, the direct forcing terms are derived by directly considering the discretised momentum equations governing the flow, including the IB. This formulation has also been adapted to a finite-difference framework, by Fadlun et al. [2] and is shown to be stable for larger time-steps. Moreover, this method does not rely on ad hoc parameters.
Lima E Silva et al. [21] proposed a variation of the direct forcing formulation, named “Physical Virtual Model” (PVM), which algebraically evaluates the required forcing at each of the Lagrangian points from the local flow equations, instead of doing it at the Eulerian points. To calculate the forces at the Lagrangian points from the flow equations, the derivatives for velocity and pressure that appear in the flow equations are approximated using second-order Lagrange polynomials. The computed forces at the Lagrangian points are subsequently distributed towards the nearby Eulerian points, analogously to the spreading operation of smooth-interface formulations. A drawback of this approach is that the polynomial approximations are not necessarily consistent with the discretisation of the momentum equations from which the forces are derived.
The family of IBMs which use direct forcing to calculate the boundary forces in conjunction with the spreading of these forces are classified as smooth-interface direct forcing formulations. Uhlmann [3] proposed a straightforward smooth-interface approach which follows the direct forcing formulation. An appropriate splitting of the direct forcing equation results in a simplified approach, where the interpolation to the Lagrangian points is only required from an intermediate velocity field to evaluate the IB forces. The intermediate velocity field results from a preliminary step that solves the equations that also arise from the splitting, where the equations are similar to the momentum equations without the presence of IB forces. Such a preliminary step can be incorporated into a framework which is based on a fractional-step time advancement [22], without adding a significant computational cost. Since the smooth-interface formulation is used to impose the boundary forces, the temporal variations of the drag on a moving IB are smooth, unlike those with sharp-interface formulations. However, more recently [23], [24] it was shown that the implementation of Uhlmann [3] results in a violation of the no-slip condition in simulations where the Reynolds number is low or the CFL number is large. For example, Kempe and Fröhlich [23] show that with the formulation of Uhlmann [3], the simulation of the flow past a stationary cylinder with a CFL number around unity leads to an error in the no-slip enforcement of 13% around the stagnation point. To improve the enforcement of the no-slip condition at every discrete time level along the IB, more recent works apply the concept of multi-direct forcing [25], [26], [23], [27], which iteratively corrects the Lagrangian forces within the predictor step of the fractional-step framework, before solving the equation for pressure correction.
Most of the smooth-interface direct forcing IBMs, which treat the forces as momentum source terms in the flow equations rely on the fractional-step framework, e.g. [3], [28], [29], [30], [24], since this framework permits an efficient implementation of the IB forcing as discussed above. However, the IBMs using this framework typically enforce the no-slip condition onto the intermediate velocity field, before a correction to satisfy the continuity equation is applied. Given that the temporal change of the pressure gradient across the IB may be large and the continuity equation is approximated by a pressure Poisson equation, the corrected velocity field can significantly violate the no-slip condition. This leads to a permeability through the IB. Even with the use of the multi-direct forcing procedure, the approximation of a flow with a low Reynolds number, a large CFL number or transient changes of the flow field results in a significant permeability through the IB [31]. The adaptation of the above schemes towards pressure-correction-based frameworks, which iteratively correct both the velocity and pressure fields until they converge at every time-step, without a significant increase in computational effort is not straightforward [32]. For frameworks that use pressure-correction schemes, Blais et al. [33] and Yi et al. [34] have put forward an algorithm which iteratively corrects both the IB forces alongside the pressure-velocity coupling until both the coupling and no-slip condition are satisfied within a specified tolerance. However, because Blais et al. [33] apply the correction of the forces in a similar form as the transient term, the convergence rate of the non-linear iterations for a time-step strongly depends on the significance of the transient term with respect to the other terms in the momentum equations.
Most IBMs are developed with a staggered variable arrangement for the pressure and velocity [34]. It is generally accepted that adopting a staggered variable arrangement is the optimal choice for Cartesian meshes [35]. However, for general mesh elements a colocated variable arrangement has significant advantages [36], [37]. To prevent the decoupling of the pressure and velocity fields which can be produced by Navier–Stokes solvers using colocated variable arrangements, momentum-weighted interpolation (MWI) is typically required to determine the fluid fluxes at cell faces [38]. However, the application of MWI in conjunction with the forces resulting from the IB may lead to complications, as these forces typically change rapidly in space and time. This is illustrated, for instance, in the work of Pinelli et al. [5], where the velocity field is shown to exhibit unphysical oscillation in the close vicinity of the IB, even though the no-slip condition is weakly enforced. Most IBMs that use a colocated variable arrangement are based on the sharp-interface approach, which typically replaces the momentum equations with interpolation equations at the forced Eulerian points [39], [4], [11], thus avoid the treatment of the IB forces as source terms. In order to remedy the decoupling for a smooth-interface IBM, where the forces are necessarily treated as source terms, Yi et al. [34] proposed the use of linear interpolation of the velocity at cell faces instead of using the MWI for cells with non-zero IB forces. Such an implementation works best for the standard 4-point discrete delta function which satisfies the even–odd condition.
There are many computational flow studies that benefit from the combination of an unstructured and fixed mesh, and a moving and possibly deforming IB. This is for instance demonstrated in the work of Roman et al. [40], where a more efficient simulation is achieved when a boundary-fitted curvilinear mesh is used to discretise the flow domain inside an S-channel, in which the IBM is used to approximate the boundary of an internal valve. As mentioned above, the discrete delta functions, which are widely used in smooth-interface IBM, have been derived for their application on uniform and Cartesian meshes. The application on other types of mesh requires modification of the interpolation weights. Pinelli et al. [5] adopted Reproducing Kernel Particle Methods (RKPM) [41] to redistribute the interpolation/spreading weights so that they satisfy the moment conditions and force conservation when the discrete delta functions are applied to curvilinear meshes. This has been successfully applied to an unstructured mesh framework [42].
Based on the recent developments in the literature discussed above, there is a need to further improve the formulation and the implementation of the IBM in order for the method to be relevant for many types of applications and more general numerical frameworks. To model a moving solid–fluid interface, IBM frameworks that use the smooth-interface approaches have been shown to be robust as well as computationally efficient, in comparison with sharp-interface approaches. Almost all IBM frameworks have been developed for Cartesian meshes and therefore use a staggered variable arrangement. When simulations on a non-Cartesian mesh are considered, a staggered variable arrangement becomes rather cumbersome [43] and a colocated variable arrangement is more favourable [36], [37]. However, the application of the smooth-interface IBM for a colocated variable arrangement is still in its infancy.
This paper presents a smooth-interface direct forcing scheme that directly uses the discretised momentum equations to evaluate the IB forces at the Lagrangian points. In the context of the finite volume discretisation, this paper also addresses the applicability of the direct forcing IBM as part of a fully-coupled pressure-velocity framework, where mass and momentum conservation are solved simultaneously, with a colocated variable arrangement and arbitrary mesh elements [36]. The main features of the IBM proposed in this paper can be summarised as follows: The forces resulting from the IB are determined directly at the Lagrangian points from the interpolated Eulerian flow variables and their derivatives that appear in the momentum equations without further approximations in time and space. The estimated Lagrangian forces are communicated to the Eulerian mesh in the form of a regularised source field, via a spreading step. The application of this force, which is consistent with the discretisation of the momentum equations, results in an immediate and accurate satisfaction of the no-slip condition on the IB, regardless of the CFL or Reynolds numbers. Additionally, the present formulation is not limited to Cartesian meshes but is applicable to any types of meshes.
Alongside this framework inspired from widely used features of the fields of IB modelling and – more generally – CFD, this manuscript presents the following novel contributions: Firstly, a method that is based on the work of Pinelli et al. [5] is proposed to compute the volumes associated with Lagrangian markers. The approach is shown to be stable and accurate for a wide range of ratios between the distance separating the Lagrangian markers and the Eulerian mesh spacing. Additionally, it allows to account for the arbitrary nature of the Lagrangian forces acting on the IB, as well as the specificities of the present implementation of the Eulerian momentum source terms. Secondly, a rigorous approach for the treatment of the source terms arising from the presence of the IB, based on the work of Denner and van Wachem [36] and Bartholomew et al. [44], is proposed. It allows pressure gradient and IB source terms to numerically balance each other, and is shown to prevent the development of unphysical flow patterns in the close vicinity of the IB.
In order to validate and assess the accuracy of the present formulation, a number of flow configurations are considered and the results obtained for these are validated with findings from the literature. The flow past a stationary cylinder as well as an oscillating cylinder is studied. These two configurations are similar to the typical benchmark configurations applied in the literature. A third case concerns an initially uniform flow, which is suddenly confronted with a stationary cylinder. A case where a cylinder oscillates within a wall-bounded flow domain is also considered in order to demonstrate the ability of the present formulation to conserve the mass of the displaced fluid as the IB moves. Finally, the unsteady flow past a sphere is considered. The above cases are studied on both Cartesian and unstructured meshes.
The present paper is organised as follows: Section 2 describes the framework of the fully-coupled flow solver with colocated variable arrangement and a generic description of the smooth-interface IBM. This is followed by the formulation of the present IBM. In Section 3, the computed results for cases with stationary and moving IB cases are presented, from which the validity, robustness and accuracy of the present method are demonstrated. Finally, conclusions are presented in Section 4.
Section snippets
Numerical framework
The flow of an incompressible Newtonian fluid, with constant fluid properties, is governed by the Navier–Stokes equations: where is the vector of fluid velocity, p is the pressure, is the source term that represents the force due to the immersed boundary (IB), and ρ and μ are the density and dynamic viscosity of the fluid, respectively. The continuity equation is given by
Using a finite-volume spatial discretisation and the first-order backward Euler
Validation and results
In this section, the presented smooth-interface direct-forcing IBM is validated and analysed in terms of accuracy, convergence and robustness. The present IBM is compared with different implementations from previous studies for a variety of flow cases. The rate of convergence of the method with regard to the refinement of the Eulerian mesh is also established. The validation of the framework is done by considering a number of two- and three-dimensional flow cases past a stationary or
Conclusions
This paper introduces a new formulation of the Immersed Boundary Method (IBM), named the Regularised Discrete Momentum Forcing (RDMF). The RDMF evaluates the boundary forces at the Lagrangian control points by directly using the interpolated algebraic discretised terms of the momentum equations. These interpolated discretised terms are applied in the direct forcing equation. The force evaluation is formulated without further approximation in time or space, other than using a discretised
Acknowledgements
M.H.A.A. acknowledges the financial support from the Government of Malaysia through Ministry of Higher Education Malaysia and Universiti Teknologi Malaysia. Advice and support from Fabian Denner on working with the main numerical framework (www.multiflow.org) adopted in this work are gratefully acknowledged. The computations were performed using the facilities of the Imperial College High Performance Computing Service.
References (62)
Flow patterns around heart valves: a numerical method
J. Comput. Phys.
(1972)- et al.
Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations
J. Comput. Phys.
(2000) An immersed boundary method with direct forcing for the simulation of particulate flows
J. Comput. Phys.
(2005)- et al.
A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries
J. Comput. Phys.
(2008) - et al.
Immersed-boundary methods for general finite-difference and finite-volume Navier–Stokes solvers
J. Comput. Phys.
(2010) - et al.
The immersed boundary method: a projection approach
J. Comput. Phys.
(2007) - et al.
An adaptive version of the immersed boundary method
J. Comput. Phys.
(1999) Modeling complex boundaries using an external force field on fixed Cartesian grids in large-eddy simulations
Comput. Fluids
(2004)- et al.
Derivation and validation of a novel implicit second-order accurate immersed boundary method
J. Comput. Phys.
(2008) - et al.
A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations
J. Comput. Phys.
(2011)
Sources of spurious force oscillations from an immersed boundary method for moving-body problems
J. Comput. Phys.
An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries
J. Comput. Phys.
An immersed boundary method with formal second-order accuracy and reduced numerical viscosity
J. Comput. Phys.
Modeling a no-slip flow boundary with an external force field
J. Comput. Phys.
Stability characteristics of the virtual boundary method in three-dimensional applications
J. Comput. Phys.
An immersed-boundary finite-volume method for simulations of flow in complex geometries
J. Comput. Phys.
Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method
J. Comput. Phys.
Application of a fractional-step method to incompressible Navier–Stokes equations
J. Comput. Phys.
An improved immersed boundary method with direct forcing for the simulation of particle laden flows
J. Comput. Phys.
A pre-conditioned implicit direct forcing based immersed boundary method for incompressible viscous flows
J. Comput. Phys.
A direct-forcing fictitious domain method for particulate flows
J. Comput. Phys.
A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows
J. Comput. Phys.
An immersed boundary technique for simulating complex flows with rigid boundary
Comput. Fluids
On the application of immersed boundary, fictitious domain and body-conformal mesh methods to many particle multiphase flows
Int. J. Multiph. Flow
A semi-implicit immersed boundary method and its application to viscous mixing
Comput. Chem. Eng.
Fully-coupled pressure-based finite-volume framework for the simulation of fluid flows at all speeds in complex geometries
J. Comput. Phys.
A ghost-cell immersed boundary method for flow in complex geometry
J. Comput. Phys.
An improved immersed boundary method for curvilinear grids
Comput. Fluids
An unstructured solver for simulations of deformable particles in flows at arbitrary Reynolds numbers
J. Comput. Phys.
Conservation properties of a new unstructured staggered scheme
Comput. Fluids
Unified formulation of the momentum-weighted interpolation for collocated variable arrangements
J. Comput. Phys.
Cited by (24)
A versatile sharp boundary ghost-node method for moving rigid boundary fluid flow with meshless nodes distribution
2024, Engineering Analysis with Boundary ElementsMicrostructure-based prediction of hydrodynamic forces in stationary particle assemblies
2024, International Journal of Multiphase FlowDrag, lift and torque correlations for axi-symmetric rod-like non-spherical particles in locally linear shear flows
2024, International Journal of Multiphase FlowA hybrid immersed boundary method for dense particle-laden flows
2023, Computers and Fluids