Streamline penetration, velocity error, and consequences of the feedback immersed boundary method

The immersed boundary method (IBM) has attracted growing interest in the computational fluid dynamics (CFD) research community due to its simplicity in dealing with fluid–structure interaction (FSI) systems involving complex geometries and large deformations.^1–3 The IBM was first developed by Peskin⁴ for modeling FSI of heart flows. In this method, the fluid dynamics is solved on a fixed Cartesian grid that does not conform to the fluid–structure boundaries, significantly reducing the effort in mesh generation and avoiding the mesh movement for moving boundaries. A boundary force density (also known as Lagrangian force) is distributed to the fluid (also known as Eulerian body force density) in the vicinity of the moving boundaries to account for their interactions.^1,2,4 The Lagrangian force is calculated by considering the principle of virtual work of the elastic body, that is, $Fib=−δE/δX$⁠, where E is the elastic energy function at the boundary configuration X.^3,5 The fluid dynamics is described in the Eulerian form, while the boundary is tracked in the Lagrangian form for the convenience of treating moving boundary cases.^1,4

Since Peskin's pioneer work, great effort has been made in improving the capability of the IBM. For example, Goldstein et al.⁶ extended the IBM by introducing a feedback forcing scheme with two feedback coefficients α and β so that the Eulerian body force density $f=α∫0tu(τ)dτ+βu(t)$⁠, where u is the fluid velocity around the boundary. This scheme is efficient and can be used to simulate complex geometry and turbulent flows.^3,6 From the physical point of view, the α term acts like a spring force on a small volume of fluid around the boundary; and the β term provides a damping effect.^1,6 The two feedback coefficients should be large enough to maintain the accuracy and small enough to avoid a system with a very large stiffness to ensure the numerical stability.^2,3,6–9 A direct-forcing IBM was proposed by Mohd-Yusof¹⁰ and Fadlun et al.¹¹ In this method, the body force density is directly deduced from the momentum equation, which can be taken as a special version of the feedback IBM. As a result, the feedback coefficients are fixed. On the other hand, several sharp-interface IBMs have been developed and applied to many topics such as fish swimming, insect/bird flight, and cardiovascular flows.^12–17 While sharp-interface IBMs have higher-order accuracy, special techniques are required to reduce the numerical oscillations due to the violation of the geometric conservation law near the immersed boundary.^18,19 Mittal and Bhardwaj²⁰ recently reviewed the progress in the development and application of diffused-interface and sharp-interface IBMs to thermofluids problems and discussed their computational challenges.

To improve the computational efficiency, the IBM was combined with the lattice Boltzmann method (LBM) by Feng and Michaelides²¹ to simulate fluid–particle interactions. The immersed boundary-lattice Boltzmann method (IB-LBM) has since been implemented in many applications, including blood flow,^22,23 viscoelastic fluids with complex geometries,^24,25 fish swimming,^7,26,27 flapping wings,^28,29 and flows at moderate and high Reynolds numbers.^30,31 In addition, a direct-forcing IB-LBM was also developed.^7,32,33

Although the feedback IBM (including the direct-forcing version) has the advantage of simplicity, its accuracy and streamline penetration have been widely discussed.^34–42 Several strategies have been proposed to enhance velocity accuracy near the immersed boundary, such as introducing an iterative IBM (also known as multi-direct forcing IBM) into LBM⁴³ and implicit IB-LBM.⁴⁴ The iterative IB-LBM iterates the IBM calculation until the velocity error less than a prescribed threshold. Zhang et al.⁴⁵ showed that usually, five iterations are required to accurately enforce the no-slip boundary condition for flow over a stationary cylinder. In many of the above works, obvious penetration of streamlines through the fluid–solid interface that was generated by their original feedback IBM implementation has been reported. However, our simulations have not observed such severe violations of the no-slip and no-penetration velocity boundary conditions. Therefore, we aim to provide an important guideline for the correct use of the feedback IBMs by comparing the performance of three implementations in modeling a few benchmark cases. In addition, the consequences of the boundary velocity error and the spurious flow penetration have not been discussed, especially for internal flows. These are the motivations for this study.

In this study, we compare the non-iterative IBMs with and without the prediction step and the iterative IB-LBM to clarify the source of streamline penetration and correct the implementation of the IBM. In addition, we discuss the velocity error and the spurious flow penetration and their consequences in external and internal flows. This study will provide an important guideline for the correct use of the feedback IBM.

The remaining parts of this paper are organized as follows: Sec. II introduces three implementations of the feedback IB-LBM. Then, detailed case testings of the numerical method and discussions are provided in Sec. III. Finally, conclusions are given in Sec. IV.

In this study, the lattice Boltzmann method (LBM) with the multi-relaxation-time (MRT) model is adopted for the fluid dynamics. For two-dimensional (2D) FSI cases (e.g., flow-induced vibration of a flexible plate attached behind a stationary cylinder in a channel and a one-sided collapsible channel), the nonlinear dynamics of the structure is treated as a Bernoulli–Euler beam and solved by the finite difference method (FDM). In the three-dimensional (3D) FSI case, the nonlinear dynamics of the structure is solved by the finite element method (FEM). The two-way fluid–structure interaction is implemented by the feedback IBM. The details of these methods are given in Secs. II A–II C. For simplicity, only the 2D numerical method is presented; the extension to 3D is straightforward, which can be also found in our previous publications.^29,46–48

The governing equations for the incompressible fluid flow are

where u is the fluid velocity, ρ is the fluid density, p is the pressure, μ is the dynamic viscosity, and f is the body force density. In the LBM, the computational domain is discretized with a fixed Eulerian grid (lattice grid). The fluid is modeled as a set of fictive particles undergoing streaming and collision over a lattice grid. The macroscopic dynamics of the fluid is the result of the statistical behavior of the particles, which is described by the distribution function $gi(x,t)$ according to^49,50

where $gi(x,t)$ is the distribution function for particles with velocity $ei$ at position x and time t, $Δt$ is the time increment, $Ωi(x,t)$ is the collision operator, and G_i is the forcing term accounting for the body force f. The D2Q9 model is used on a square lattice. “D2” stands for “two dimensions,” while “Q9” stands for “nine particle speeds.” The discrete velocity components of D2Q9 model can be represented as

where $cΔ=Δx/Δt$ with $Δx$ the lattice spacing. Compared with the single relaxation time (SRT) collision model, the multiple relaxation time (MRT) model has been proven to be more numerically stable owing to shear and bulk viscosities are allowed to tune independently.⁴⁹ Increasing the bulk viscosity can often improve the numerical stability by attenuating any spurious pressure waves more quickly.⁵¹ Therefore, the MRT collision model is adopted here and is given by⁴⁹ 

where M is a 9 × 9 transform matrix, and $S=diag(s0,s1,…,s8)$ is a non-negative diagonal 9 × 9 relaxation matrix for D2Q9 model. In matrix S, s₁ is related to the bulk viscosity of the fluid, s₇ and s₈ are related to the shear viscosity of the fluid, and other parameters are free parameters. The physical meaning and values of these parameters be found in Lallemand and Luo⁴⁹ and Luo et al.⁵⁰ The equilibrium distribution function $gieq$ is defined as

where $cs=Δx/(3Δt)$ is the speed of sound, I is the unit tensor, and the weighting factors ω_i are given by $ω0=4/9$⁠, $ω1−4=1/9$⁠, and $ω5−8=1/36$⁠. The mass density ρ, pressure p, and velocity u are, respectively, calculated by

The force scheme proposed by Guo et al.⁵² is adopted to determine G_i,

where $τ=1/s8$ is the non-dimensional relaxation time.

The geometrically nonlinear dynamics of a 2D elastic wall is described by the Bernoulli–Euler beam^48,53

where ρ_s is the linear density of the elastic wall, s is the arc length along the elastic wall, X is the position vector of the Lagrangian point on the elastic wall, T(s) is the longitudinal tension, EI is the bending stiffness (where E is the Young's modulus, $I=h3/12$ is the moment of inertia of the wall cross section, and h is the wall thickness), and F is the hydrodynamic stress exerted by the fluid. The longitudinal tension T(s) is evaluated by the Hooke's law

where Eh is the stretching stiffness. In this study, the FDM is used to discretize the governing equation [i.e., Eq. (10)] of the elastic wall. The mesh size of the solid wall is half of the lattice spacing to avoid “fluid volume leakage.”^3,5,51 The details of the numerical method can be found in Huang⁴⁸ and Ma et al.²⁵ 

The two-way interactions between the fluid and the structure are coupled by the feedback law⁵⁴ 

where $Fib(s,t)$ is the Lagrangian force density, β is the feedback coefficient and $β=2 m/s$ in LBM simulations, $Uib(s,t)$ is the immersed boundary velocity interpolated from the ambient flow, and $U(s,t)$ represents the velocity of the solid wall, $U(s,t)=0$ for a rigid solid wall. In dimensionless form, $β*=β/U0=40$⁠, and $β*$ ranges from 20 to 104.²³ The derivation of the dimensionless form of β can be found in  Appendix A. It is noted that this method is simple, and has been extensively used in both laminar and turbulent flows (see, e.g., Huang and Tian,³ Huang et al.,²³ and Xu et al.³⁰).

A Dirac delta function is used to transfer the interactions between the Lagrangian and Eulerian variables. The velocity interpolation and Lagrangian force distribution to the adjacent Eulerian grid points are calculated according to

where $u(x,t)$ is the fluid velocity, x is the coordinate of the fluid lattice nodes, $δ(x−X(s,t))$ is Dirac delta function, X is the coordinate of the solid wall, f is the body force added in the Navier–Stokes equation to mimic a boundary condition, and ds is the arc length of the immersed boundary.

Because the location of the immersed boundary does not generally coincide with the nodal points of the Cartesian grid, the boundary force needs to be distributed over a band of lattices around each Lagrangian point, resulting in a diffusive boundary.¹ The Dirac delta function is a smooth distribution function, serving as a Gaussian-like function in transferring boundary information between Lagrangian and Eulerian points. The four-point discrete delta function $δh(x−X(s,t))$ developed by Peskin⁵ is used to approximate the Dirac delta function

where $r=(x−X(s,t))/Δx$ for x-component, and $(y−Y(s,t))/Δy$ for y-component.

Here, we consider three implementations of the feedback IB-LBM, summarized in Algorithms 1–3. Algorithm 1 is the conventional feedback IB-LBM without a prediction step. Algorithm 2 is an improved IB-LBM with a prediction step. Algorithm 3 is an iterative feedback IB-LBM. Algorithm 1 could be taken as a special case of Algorithm 3 with one iteration. Algorithm 2 can be taken as a variation of Algorithm 1 by slightly modifying the integration scheme. The computational complexities of Algorithms 1–3 are, respectively, O(N), O(N), and O(mN), where N is the number of the Lagrangian points and m is the total iteration number, see also Zhao et al.⁴² For each algorithm, a partitioned and weakly coupled approach of the fluid and structure solvers is adopted, where each solver is solved sequentially only once at each time step. Consequently, the fluid–solid interface boundary conditions are mismatched by one-half step at the end of each time step. This approach is computationally efficient, while it may cause numerical instability at low structure-to-fluid mass ratios.^15,23 To improve the numerical stability, an iterative feedback IBM, as shown in Algorithm 3, is used at the fluid–solid interface, which allows for local flow reconstruction near the solid boundary. The iteration ensures that the displacement, velocity, and traction boundary conditions at the fluid–solid interface are matched between the fluid and the solid boundary at each time step.²³ Note that in Algorithm 3, the pre-set criterion is an empirical parameter. At low structure-to-fluid mass ratios, the numerical stability can also be improved by implicit/strong coupling.²⁰ 

The flow over a cylinder, flow-induced vibration of a flexible plate attached behind a stationary cylinder in a channel, the flow through a 2D asymmetric stenosis, a one-sided collapsible channel, and a 3D collapsible tube are conducted here to test the proposed methods. Note that two typical groups of cases are particularly covered: internal/confined and external flows including both steady and unsteady flows. As detailed below, the Reynolds numbers considered are not arranged to be uniform across the cases presented, as the focus here is on the streamline penetration, velocity error, and consequences of three implementations of IBM in different types of flows.

A uniform flow past a 2D cylinder is one of the most robust methods used for validating the fluid solvers.⁵⁵ Therefore, this problem is used as a benchmark test for the three algorithms. Figure 1 shows a schematic depiction for the simulation of the flow over the cylinder with boundary conditions and dimensions used. Note that in this work, we only consider flows at low Reynolds numbers and thus do not implement the nonreflecting boundary condition at the outlet.

The Reynolds number Re, based on the cylinder diameter D, undisturbed free-stream velocity U₀, fluid density ρ, and dynamic viscosity μ, determines the flow pattern and drag on a cylinder

The non-dimensional drag C_D and lift C_L coefficients are defined to describe the hydrodynamic forces on the cylinder

where F_D and F_L are the drag and lift forces, respectively. The F_D and F_L are integrated from the forces exerted on the immersed boundary by the ambient fluid

where $Fxi$ and $Fyi$ are the forces acting on the i-th $IB$ point in $x$ and $y$ directions, and $Δs$ is the cylinder mesh size.

A multi-block LBM⁵⁶ is adopted to provide high resolution near the cylinder while remaining low resolution in the far-field to reduce the computational effort. Here, two blocks of the fluid mesh are used. The finest mesh size around the cylinder is $D/100$⁠. The cylinder mesh size is maintained at half of the fluid mesh. The flow over the cylinder has been simulated at $Re=40$⁠, corresponding to a steady flow regime. It, respectively, requires 5.083, 4.633, and 4.25 CPU min for Algorithms 1–3 for calculating 1000 time steps.

We first compare the streamlines of the three implementations of the feedback IB-LBM with the conventional feedback coefficient $β=2 m/s$ and the optimal (in practical meaning) feedback coefficient $β=5.2 m/s$ where the conventional feedback coefficient is amplified 2.6 times. Here, $U0=0.05 m/s$ is used for all tested cases. Therefore, the non-dimensional feedback coefficient $β*=β/U0=40$ for $β=2 m/s$⁠, and $β*=104$ for $β=5.2 m/s$⁠. Figure 2(a) shows that the streamlines exhibit symmetrical vortices behind the cylinder. The vortex center location (a, b) is defined as an indication of the consequences of the streamline penetration, which will be discussed in Table I. Figures 2(b)–2(h) show the streamlines for the Algorithms 1–3, respectively. As is shown in Fig. 2(b), many streamlines penetrate through the cylinder (i.e., spurious flow penetration as described by He et al.⁵⁷) for the conventional feedback IB-LBM (Algorithm 1, $β=2 m/s$⁠), which means the no-penetration boundary condition is not precisely satisfied at the fluid–solid interface. Figure 2(c) shows that the apparent streamline penetration still exists for Algorithm 1 with $β=5.2 m/s$⁠. Similar spurious flow penetration can be also found in the improved (Algorithm 2, $β=2 m/s$⁠) and the one iteration (Algorithm 3, $β=2 m/s$⁠) feedback IB-LBM [i.e., Figs. 2(d) and 2(f)]. By contrast, there is no streamline penetration through the cylinder for the improved (Algorithm 2, $β=5.2 m/s$⁠) and five iterations (Algorithm 3, $β=2 m/s$⁠) feedback IB-LBM [i.e., Figs. 2(e) and 2(g)]. Figure 2(h) shows that there is also no streamline penetration through the cylinder for the one iteration (Algorithm 3) feedback IB-LBM with $β=5.2 m/s$⁠.

The observations suggest that either iterations or a larger feedback coefficient can suppress the spurious flow penetration and improve the no-penetration boundary condition at the fluid–solid interface. These observations have also been reported in previous studies.^35,37,38,58 Compared with Algorithm 1, the velocity prediction step in Algorithms 2 and 3 can help to suppress the spurious flow penetration and reduce the velocity slip on the immersed boundaries.

Table I shows the results from the present simulation and those published in the literature for comparisons. A reasonably good agreement is found between the present results and other numerical simulations and experimental results for all quantities of interest. For the vortex center location (a, b) and drag coefficient C_D computed by present three algorithms, the difference is negligibly small, indicating that iterations do not significantly improve the prediction of the drag and the vortex center location. In the rest of this paper, we only compare Algorithm 3 for different iteration strategies and feedback coefficients.

Here, a moving boundary case, the FSI of a flexible plate behind a stationary cylinder in a channel, is considered. Figure 3 shows the schematic diagram of the geometry and the boundary conditions of this case. This case has been used as a benchmark validation for FSI solvers involving large displacement.^15,64–67

As shown in Fig. 3, a fixed circular rigid body is submerged in an incompressible fluid. A flexible thin beam is attached downstream to the cylinder. The cylinder has a diameter of D and is centered at the origin. The beam has thickness h and length $L=3.5D$⁠. A parabolic velocity profile with averaged velocity U₀ and a constant pressure P₀ is imposed at the inlet. No-slip walls are enforced at the top and bottom sides of the computational domain. The normal and shear stress are set to zero at the outlet. The computational domain is a rectangular box (⁠$x∈[−2D,9D]$ and $y∈[−2.05D,2.05D]$⁠), and the grid size for the fluid and the beam is 0.01D and 0.005D, respectively. The non-dimensional parameters for this case are

where ρ_s is the linear density of the beam, M is the structure-to-fluid mass ratio, and K_b and K_s are the bending stiffness and stretching stiffness of the beam, respectively. E is Young's modulus, and $I=h3/12$ is the moment of inertia of the beam cross section. Here, $Re=100,M=2.0$⁠, and K_b = 1.111. A large stretching stiffness K_s = 500 is chosen to achieve a nearly inextensible beam.⁶⁶ 

Figure 4 compares the time histories of the x- and y-coordinate at the trailing end. It shows that the flexible beam reaches a periodic oscillation at around $tU0/D=40$⁠. A reasonably good agreement is found between the present results and other numerical results of literature.^{15,64,66–68}

Table II shows the comparison of the mean drag coefficient $CD,m=FD,m/(0.5ρU02D)$ (⁠$FD,m$ is the mean drag of the cylinder-beam system), Strouhal number defined as $St=fD/U0$ (f is the oscillation frequency), and vertical oscillation amplitude of the trailing end. Overall, the present results show reasonable agreement with previous results.

The Strouhal number St, which characterizes the flapping shedding frequency, is 0.186 in the present Algorithm 3 with one iteration and $β=2 m/s$ simulation, while that computed by previous studies^64,66,68 is 0.19 (a difference of about 2.1%). In present simulation of Algorithm 3 with five iterations and $β=2 m/s$⁠, St = 0.185 (a difference with previous results of 2.6%), suggesting that the iterations do not affect the prediction of the flapping frequency of the beam. For the one iteration with $β=5.2 m/s$⁠, the prediction is also not affected. Given that three significant digits of St are kept here, and a thin plate is used instead of the plate of finite thickness, the deviations of the frequency predicted by our simulations are reasonable. This case also shows the capability of this solver in modeling FSI involving a thin flexible plate.

Many studies examined the iteration effects of the IBM in the scenarios of external flows (e.g., flow over a cylinder^38,45). However, there is no study on the iteration effects of IBM in internal flows. Therefore, the iteration effects of the IBM are tested here for fluid flow through an asymmetric stenosis with a diameter restriction of 50% at the constriction. Figure 5 shows the schematic diagram of the stenosis, where two red lines represent the upper and lower channel walls. A cosine function dependent on the axial coordinate x is used to describe the upper stenosed channel wall

where D is the diameter of the non-stenosed channel, $a0=0.25$ for the 50% diameter reduction, $x0=6$ is the x coordinate of the center of the stenosis (⁠$x0−L/2≤x≤x0+L/2$⁠), and $L=2D$ is the length of the stenosis. The length and width of the whole computational domain are 16D and $1.2D$⁠, respectively. A steady Poiseuille flow with an averaged velocity U₀ is imposed at the upstream inlet, and a constant pressure p_d is specified at the downstream outlet. The grid size for the fluid and the channel wall is $dx=0.01D$ and $ds=0.005D$⁠, respectively.

The fluid flow through a 2D asymmetric stenosis at Re = 200 is simulated and is validated against the commercial software ANSYS Fluent. In the ANSYS Fluent, a boundary-fitted grid with a total grid size of 498 × 100 quadrilateral elements is generated. A steady Poiseuille flow with averaged velocity U₀ imposed at the inlet using the user-defined function (UDF). The flow is laminar and steady, which is solved by pressure–velocity coupling scheme. This scheme has the advantage of solving the momentum and pressure correction separately, which simplifies computations but can led to slower convergence for poorly constructed grids.⁶⁹ For the space discretization, a least squares cell-based, body force weighted, second-order upwind scheme is chosen for the gradient, pressure, and momentum, respectively. As the flow is simple (i.e., laminar and steady) and at low Reynolds number Re = 200, second-order space discretization is enough for resolving all concerned physical quantities such as the pressure and velocity. The convergence standard $ϵ=1×10−5$ is chosen.

Figure 6 shows the simulated pressure and streamwise velocity contours at $tU0/D=100$ (the flow reaches a steady state) for one iteration (i.e., Algorithm 3 with one iteration and $β=2 m/s$⁠) and five iterations (i.e., Algorithm 3 with five iterations and β = 2 m/s), respectively. As shown in the pressure contours, a low-pressure area is observed at the posterior part of the stenosis. Compared with the pressure contours for one iteration, a much higher pressure region is observed in the upstream channel of the stenosis for five iterations. The velocity contours show that a stable jet flow is formed downstream of the stenosis. The velocity increases significantly due to the constriction of the upper channel wall. The jet flow is stronger for the five iterations. The streamlines penetrate through the upper and lower channel walls in the velocity contours for one iteration, but not for five iterations.

These observations again demonstrate that the iterative IBM can suppress the spurious flow penetration and improve the no-penetration boundary conditions at the walls. As a result, the mass conservation in the channel is more strictly satisfied, which is crucially important for simulating internal flows.

Figure 7 shows the velocity profiles at four different axial positions at x = 5.85 and x = 7. As shown in Fig. 7, the velocity profiles predicted by the present five iterations and one iteration with β = 5.2 m/s agree very well at all the axial positions with the numerical solutions of the commercial software ANSYS Fluent. By comparison, big discrepancies are observed at all the axial positions for the present one iteration with β = 2 m/s. Small discrepancies are observed at all axial positions for the present three iterations with β = 2 m/s. The under-prediction of velocity for the one iteration and three iterations with β = 2 m/s is due to fluid leakage from the channel as a result of spurious flow penetration.

The pressure and wall shear stress (WSS) on the arterial wall are of great interest to the medical community as they play an essential role in the genesis and progression of cardiovascular diseases.^70,71 Therefore, the distributions of pressure and WSS along the upper channel wall are shown in Fig. 8. The pressure and WSS are linearly interpolated based on the corresponding values at 2.5 and 5.0 grid points inward of the channel walls.^23,72 Results of one iteration with β = 2 m/s largely under-predict the pressure and WSS because the no-slip and no-penetration boundary conditions at the channel walls are not exactly satisfied. This issue can be well addressed by five iterations or using β = 5.2 m/s of the IBM. Small discrepancies are observed for three iterations. The iteration allows for local flow reconstruction in the vicinity of the channel walls,²³ and the results agree very well with those of ANSYS. Thus, the observations suggest that the iteration and β = 5.2 m/s both improve the enforcement of the no-slip and no-penetration boundary conditions on the channel walls.

Here, the iteration effects of the IBM are examined for an internal flow with a moving and deformable boundary by considering a 2D incompressible flow in a one-sided collapsible channel. As shown in Fig. 9, a part of the channel wall is replaced by an elastic beam. The elastic beam has length L and is subjected to an external pressure p_e. The rigid channel has a width of D. A steady Poiseuille flow with averaged velocity U₀ is imposed at the upstream inlet, and a constant pressure p_d is specified at the downstream outlet.

The averaged flow velocity at the inlet U₀, channel height D, and fluid density ρ is used to non-dimensionalize this system, giving five non-dimensional parameters: the Reynolds number, the structure-to-fluid mass ratio, the stretching stiffness, the bending stiffness, and the external pressure, which are, respectively, given by

Here, Re = 250, M = 1, $Ks=56.88$⁠, and $Pe=1.95$ are used. A no-slip boundary condition is applied along the channel wall, including the elastic segment. Clamped conditions are used at the two ends of the elastic wall. The remaining parameters are $Lu=5D, Ld=30D,L=5D$⁠, and $Kb/Ks=(h2/12D2)≈10−5$ for a wall thickness h of 1% of the channel height. The nonlinear dynamics of the collapsible channel wall is treated as a Bernoulli–Euler beam with zero initial tension and solved by the finite difference method. The grid size for the fluid and the channel walls is $0.01D$ and $0.005D$⁠, respectively. More computational details can be found in our previous work.^23,46,47,72 Here, the iteration effects of the IBM are examined by considering five cases: one iteration, three iterations, five iterations, dynamic iterations with β = 2 m/s, and one iteration with β = 5.2 m/s. For the dynamic iterations, the iteration is terminated when the maximum velocity error at the immersed boundary is less than a pre-set criterion (i.e., $max(Uerror(s,t))≤5×10−3$⁠)²³ 

For the dynamic iterations with β = 2 m/s, the average iteration number is 6.17. The iteration number is reduced to 1 or 2 using β = 5.2 m/s.

Figure 10(a) shows the comparison of wall shapes from present simulations with the arbitrary Lagrangian-Eulerian (ALE) method of Luo et al.⁷³ The simulations with dynamic iterations and one iteration with β = 5.2 m/s predict the most accurate wall shape, while the one iteration with β = 2 m/s predicts the worst. The prediction improves as the increase in iteration. Figure 10(b) shows the time history of the y-coordinate of the mid-point of the upper elastic wall for different iterations. It shows that the time histories of the y-coordinate are quite different not only in the final value but also in the oscillation trajectory. For the dynamic iterations and one iteration with β = 5.2 m/s, the oscillations are evident during the transition state (⁠$0≤t/T≤50$⁠), while the oscillation amplitude is gradually damped out as time increases. For the other three cases (i.e., one iteration, three iterations, and five iterations), the elastic walls are in a steady state for most of the time. The different stabilities predicted by these five iteration strategies are due to the boundary velocity slip error and the spurious flow penetration on the elastic walls for the one iteration, three iterations, and five iterations. The boundary velocity slip error and the spurious flow penetration on the elastic wall are demonstrated, respectively, in Figs. 11 and 12.

In order to examine the boundary velocity slip error, the velocity error along the elastic wall, as shown in Fig. 11, is calculated at $t/T=80$⁠. It shows that the maximum velocity error is on the collapsible part of the wall (⁠$x/D≈1.4$⁠) for all the iteration strategies. Please note that the pressure inside the channel is generally lower than the external pressure. For the one iteration case, a lower β induces a larger velocity error, which in turn reduces the pressure drop and consequently increases the pressure difference [see Fig. 12(a)], leading to a larger collapse deformation. A larger collapse deformation increases the velocity near the collapsed region leads to a lower pressure enlarging the pressure difference. As discussed in  Appendix A, the physical meaning of the Lagrangian force is to balance the fluid force, which is the pressure difference here. Therefore, according to $Δp*≈ρβΔU*$ (where $Δp*$ is the dimensionless pressure difference), a larger $Δp*$ leads to a larger velocity error, as observed in Fig. 11. In addition, the velocity error profile is similar to the distribution of pressure [not shown for this case but the distribution is close to that shown in Fig. 8(a)].

Serious boundary velocity slip error is on the channel wall from the inlet to the end of the elastic wall, and this error can be significantly (more than an order of magnitude) reduced by five iterations or one iteration with β = 5.2 m/s. For the one iteration with β = 5.2 m/s, the overall velocity error is the smallest (averaged $Uerror≈0.0001$⁠) among all the five iteration strategies. The observations also suggest that either iterations or a larger feedback coefficient can reduce the velocity slip error on the solid wall.

Figure 12 shows the pressure and velocity contours of the one iteration and the dynamic iterations, respectively. As shown in the pressure contours, a high-pressure region is developed upstream of the constriction, and a low-pressure region is observed downstream of the elastic wall. Compared with the pressure contours for the one iteration, a much higher pressure region is observed at the upstream of the constriction for the dynamic iterations, causing the upstream half of the elastic wall to bulge out. A jet flow is developed downstream of the constriction, as shown in the velocity contours. For the dynamic iterations, the high-velocity region at the constriction is more prominent and higher than that of the one iteration, indicating a stronger jet flow. For the one iteration, the streamlines penetrate through the upper and lower channel walls, while no streamline penetrates through the channel walls for the dynamic iterations. Severe streamline penetrations at the upstream elastic wall are observed, and the spurious flow penetrations are less severe downstream of the constriction, especially at $x/D>17$⁠. The one iteration predicts a larger recirculation area due to the narrower constriction of the channel.

Here, the computation of a steady flow in a collapsible tube is conducted to examine the iteration effects of the IBM. The coupled nonlinear system for the FSIs in the collapsible tube is solved by the finite element-immersed boundary-lattice Boltzmann method (FE-IB-LBM) FSIs solver. D3Q19 model has been proven to be accurate enough for low to medium Reynolds number flows,^25,30 and MRT has been proven to be a numerically stable model.^49,51 Therefore, the D3Q19 LBM with MRT model is adopted for the fluid dynamics. The structural equation is solved by an explicit FEM. The strain of the wall is assumed to be small, which enables the tube wall to be treated as a linear elastic structure. The dynamics for the solid is

where ρ_m is the solid density, $us(X,t)$ denotes the solid displacement field, c is the material damping, $σ$ is the Cauchy stress tensor, and b is the body force.

A 3D incompressible flow in a thin-walled collapsible tube of undeformed radius R, length L, and wall thickness h is considered. As illustrated by Fig. 13, the elastic wall is subjected to an external pressure p_e, and the tube wall material has Young's modulus E and Poisson ratio ν_s. The rigid tube has a diameter of D. The averaged flow velocity at the inlet U₀, tube diameter D, and fluid density ρ are used to non-dimensionalize this system, giving four non-dimensional parameters governing this FSI system: the Reynolds number, structure-to-fluid mass ratio, bending stiffness, and external pressure are given, respectively, by

A no-slip boundary condition is applied along the tube wall. The tube is clamped at both ends and is discretized with 26 400 isoparametric 3D solid elements (20 nodes for each element and one element in the thickness direction), giving 185 460 nodal points with three degrees of freedom (displacements) for each structural element. A steady Hagen–Poiseuille velocity profile with average velocity U₀ is imposed at the inlet, and a constant pressure p_d = 0 is specified at the downstream outlet. The initial flow field is initialized as

where u, v, and w are the velocity in $x, y$⁠, and $z$ directions, respectively, and $r=y2+x2$ is the radial distance from the tube centerline.

The axisymmetric solution is remarkably robust, such that a perturbation load is required in order to force the structure collapse into a mode 2 (two-lobed shape) buckling. A perturbation is initially added to the external pressure as in Marzo et al.⁷⁴ 

where X is the horizontal spatial coordinate of the tube as shown in Fig. 13, and the coefficients a₀ and a₁ are

where C is the parameter controlling the degree of perturbation and is chosen as $0.1pe$⁠. p_load was set to zero once the tube had buckled.

The computation parameters are $R=0.5 m, Lu=R=0.5 m,L=Ld=10R=5 m, νs=0.49, E=4559.4 Pa, h=0.05R=0.025 m$⁠, and $ρm=1000 kg/m3$⁠, giving the same non-dimensional parameters used by Hazel and Heil:⁷⁵ Re = 128, M = 25, $Kb=0.030 517, Pe=5.468 75$⁠. To approximate the steady flow, an empirical damping ratio $ζ=c/cc=0.05$ (⁠$cc=2Km$ is the critical damping, $K=E(h/R)3/12(1−νs2)$ is the bending modulus, and m is the mass of the elastic tube) is used here to damp out the oscillations. More computational details can be found in our previous publications.^47,48

Figure 14 shows the comparison of present computations with those of Hazel and Heil⁷⁵ and Zhang et al.⁷⁶ The simulation with dynamic iterations and one iteration with β = 5.2 m/s predicts the most accurate wall shape, while the one iteration with β = 2 m/s predicts the worst. This observation suggests again that either iterations or a larger feedback coefficient can improve the no-slip and no-penetration boundary conditions at the fluid–solid interface.

The non-iterative IBMs with and without the prediction step (Algorithms 1 and 2) and the iterative IBM (Algorithm 3) are examined for external (i.e., flow over a stationary cylinder and flow-induced vibration of a flexible beam attached behind a cylinder) and internal (flow through a 2D asymmetric stenosis, a one-sided collapsible channel, and a 3D collapsible tube) flows. Compared with Algorithm 1, the velocity prediction step in Algorithms 2 and 3 can help to suppress the spurious flow penetration and reduce the velocity slip on the immersed boundaries. The velocity error that occurred in the one iteration feedback IBM with the conventional feedback coefficient β = 2 m/s is of low consequence for external flows as it does not introduce significant deviation in predicting the forces on the structure and the structure deformation. However, it matters for internal flows, as it leads to significant deviation in the force prediction and the structure deformation. In addition, Luo et al.³⁴ showed that reducing the velocity error is also quite important in the multiphase flow simulation (e.g., particle-laden flows).

For external flows, the velocity error on the cylinder does not cause significant deviation in the pressure and shear force on the solid. Therefore, the errors in predicting the drag coefficient C_D and the structure deformation are small, while for flows in a 2D asymmetric stenosis, as shown in Fig. 8(a), the velocity error could cause a large maximum pressure error (35.94% for the case considered here) at the inlet. The boundary velocity slip causes a maximum WSS error of 25.04% at the most constriction region of the upper channel wall [Fig. 8(b)]. Therefore, for internal flows with flexible channel walls, such as the one-sided collapsible channel flow and the 3D collapsible tube flow considered here, the significant velocity error (Fig. 11) on the channel walls causes a significant error in the pressure and shear force. As a result, the error in the force prediction causes a significant error in the wall deformation (Fig. 10), which in turn enlarges the error in velocity.

Results have shown that the iterative IBM or the non-iterative IBM with a larger feedback coefficient can suppress the boundary velocity-slip error and spurious flow penetration on the solid wall. As a result, good force production and structure deformation are achieved. The iteration number and the computational time can be significantly reduced by amplifying the feedback coefficient β. The iterative IBM and the method of amplifying the feedback coefficient could effectively suppress boundary slip and spurious flow penetration. Overall, this strategy is highly effective and efficient and effortless without any extra effort. Based on our extensive testing, β = 5.2 m/s can cause serious numerical instabilities in the cases of low mass ratio and large stretching stiffness, while the iterative IBM with β = 2 m/s works well in this case and is thus more numerically stable. It should be noted that the optimal feedback coefficient will be a bit different for different flow solvers/methods, for example, different force schemes in the LBM or different delta functions.³ On the other hand, the three algorithms presented in this study are diffused-interface IBMs, which requires further effort in order to be applied to high Reynolds number flows.²⁰ 

Just as pointed out by Zhao et al.,⁴² the boundary condition-enforced IB-LBM^77,78 can accurately fulfill the no-slip boundary condition, but it is quite computational challenge when simulating moving boundary problems as a large matrix must be assembled and implicitly resolved at every time step. The feedback IBM (Algorithm 3) presented in our manuscript is an explicit method, in which calculations of the large matrix are avoid and no-slip and no-penetration boundary conditions are well fulfilled either by small feedback coefficient with several iterations (e.g., five iterations with β = 2 m/s) or the optimal feedback coefficient with only one iteration (e.g., one iteration with β = 5.2 m/s). Compared with the explicit technique-based IB-LBM proposed by Zhao et al.,⁴² the feedback IB-LBM (Algorithm 3) presented in our manuscript has an adjustable feedback coefficient β, which is very attractive for complex fluid–structure interaction problems involving large and elastic deformation, large amplitude vibration, low structure-to-fluid mass ratio, and turbulent flows. For example, for numerical instabilities often encountered in low structure-to-fluid mass ratio cases, the numerical stability can be improved by a small feedback coefficient, and the boundary slip and spurious flow penetration can be suppressed by several iterations (e.g., five iterations with β = 2 m/s).

The performance of three implementations of the feedback IBM has been studied. The streamline penetration, velocity error on the immersed boundary, and consequences in the force production and structure deformation are discussed by simulating external (a uniform flow over a cylinder, flow-induced vibration of a flexible beam attached behind a cylinder) and internal flows (flow through a 2D asymmetric stenosis, a one-sided collapsible channel, and a 3D collapsible tube). Results show that the widely reported streamline penetration can be significantly reduced by properly implementing the feedback IBM or employing the iterative IBM. The boundary velocity error does not significantly affect the force production and structure deformation for external flows. However, for internal flows such as the stenosis and the collapsible channel/tube flows, reducing the velocity error by using the iterative IBM or the optimal feedback coefficient substantially improves the prediction of the force distribution and structure deformation. Moreover, the value of the feedback coefficients could be smaller for the iterative IBM, which is very attractive for improving the numerical stability for low structure-to-fluid mass ratio cases but at the expense of more iterations. The results obtained provide an important guideline for the correct use of the feedback IBMs.

See the supplementary material for more results associated with Figs. 6–8, 11, and 12.

The authors would like to thank Dr. Lincheng Xu, Laboratory M2P2, Aix Marseille University for suggestions and discussions on the immersed boundary method. Dr. F.-B. Tian is the recipient of an Australian Research Council Discovery Early Career Researcher Award (Project No. DE160101098). This work was also supported by Asian Office of Aerospace Research and Development (Grant Nos. FA2386-19-1-4066 and FA2386-20-1-4084) and Office of Naval Research Global (Grant No. N62909-20-1-2088). The computation work of this research was performed on the National Computational Infrastructure (NCI) supported by the Australian Government.

The authors have no conflicts to disclose.

Qiuxiang Huang: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead). Zhengliang Liu: Investigation (equal); Methodology (equal); Writing – review and editing (equal). Li Wang: Investigation (equal); Methodology (equal); Writing – review and editing (equal). Sridhar Ravi: Investigation (equal); Writing – review and editing (equal). John Young: Investigation (equal); Writing – review and editing (equal). Joseph Lai: Investigation (equal); Writing – review and editing (equal). Fang-Bao Tian: Conceptualization (lead); Investigation (equal); Methodology (lead); Writing – review and editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

By incorporating the immersed boundary method [see Eqs. (12)–(14)] into the Navier–Stokes momentum equation [Eq. (2)], the momentum equation can be rewritten as

where $ΔU=U−Uib$⁠. The above equation can be non-dimensionalized by the reference ρ, the reference length L, and the characteristic velocity U₀ (i.e., $ρU02/L$⁠) as

where $β*=β/U0$ and $Re=ρU0L/μ$⁠. We may discretize Eq. (A1) with the use of Eq. (15) as

where G includes the convective, viscous, and pressure terms. Equation (A3) can be then non-dimensionalized as

which can also be obtained by discretizing Eq. (A2). From Eqs. (A2) and (A3), we need to consider two aspects in the numerical setup. The first aspect is to keep the dynamical similarity; that is, the $β*$ should be the same for different values of U₀. The other aspect is the numerical coefficient $βΔtΔy$⁠, which acts like a Lagrangian multiplier and should be much larger than one in order to achieve accurate results. Another physical explanation is that the IB force is used to balance the fluid dynamic force, that is, $ρU02=βρΔU$⁠, or in dimensionless form $ΔU*=1/β*$⁠, where $ΔU*$ is the velocity error around the fluid–structure interface. Therefore, feedback coefficients near 20 and 104 would estimate the velocity errors of about 5% and 1%, respectively. For a small coefficient, for example, $β*≤20$⁠, iteration is highly desired to enhance the non-slip boundary condition, while for a large coefficient, for example, $β*≥104$⁠, it may lead to numerical instability, especially for the LBM model used here. It is recommended to tune the non-dimensional feedback coefficient in the range of $[20,104]$ instead of random test. Note that the choice of β is more flexible in other numerical methods such as the finite difference method where $Δt$ is independent of $Δx$ except for the numerical stability requirements. Therefore, researchers can always increase β and reduce $Δt$ to achieve a higher accurate solution at the cost of higher computational expense. It is more limited in the LBM where it requires $Δt=Δx$ for most of the LBM models.

Flexible beam behind a stationary cylinder in a channel, fluid flow through a 2D asymmetric stenosis, and one-sided collapsible channel flow are chosen for the grid refinement study. Three grid sizes (lattice spacings) are conducted: $dx=0.02D, 0.1D$⁠, and $0.005D$⁠. The solid grid size is always maintained at half of the lattice spacing.

Figure 15 shows that the converged solution for the time histories of the x- and y-coordinates at the trailing end produced by $dx=0.005D$ agree well with the computational results of Wang and Tian.⁶⁶ 

Figure 16 shows that the converged solution for velocity profile and pressure distribution produced by $dx=0.005D$ agrees well with the computational results of ANSYS Fluent.

Figure 17 shows the wall shapes and y-coordinate time history for the mid-point of the upper elastic wall for the three grid refinements. The solutions are converged, and the results produced by the finest grid size $dx=0.005D$ agree well with the computational results of Luo et al.⁷³