# Immersed boundary method for the simulation of lid-driven cavity flow with an embedded cylinder

Shang-Gui Cai, Abdellatif Ouahsine, Julien Favier, Yannick Hoarau, Hassan Smaoui

► To cite this version:

Shang-Gui Cai, Abdellatif Ouahsine, Julien Favier, Yannick Hoarau, Hassan Smaoui. Immersed boundary method for the simulation of lid-driven cavity flow with an embedded cylinder. Coupled Problems in Science and Engineering VI, May 2015, Venice, Italy. pp.1130-1137. hal-01462151

HAL Id: hal-01462151

https://hal.science/hal-01462151v1

Submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
IMMERSED BOUNDARY METHOD FOR THE
SIMULATION OF LID-DRIVEN CAVITY FLOW WITH AN
EMBEDDED CYLINDER

Shang-Gui CAI*, Abdellatif OUAHSINE*, Julien FAVIER**, Yannick HOARAU† AND Hassan SMAOUI*,††

**Key words:** Immersed Boundary Method, Projection Method, Lid-driven cavity flow, Embedded Cylinder

**Abstract.** We present the numerical simulation of lid-driven cavity flow with an embedded cylinder by using a modified implicit immersed boundary method. The immersed boundary method features the non body-conforming mesh merit, which facilitates the simulation with complex geometry in a flow. The paper compares the numerical solutions between the modified implicit immersed boundary method with the traditional body-conforming mesh method as well as the results in the literature, which shows an excellent agreement.

# 1 INTRODUCTION

Immersed boundary (IB) method has emerged in recent years as an alternative to the conventional body-conforming mesh technique [1, 2], demonstrating its advantage to handle complex rigid or flexible boundaries embedded in a flow [3]. Although meshless methods have a successful applications [4], they have difficulties of their own. The establishment of a well body fitted Cartesian mesh generally possess most of the CFD work. IB method significantly simplifies this procedure and speedup the solution by eliminating the re-meshing or mesh deforming process. Through adopting an appropriate boundary force for the presence of immersed structure or solid, the simulation is simply performed

<footnote>* Sorbonne Universités, Université de Technologie de Compiègne, CNRS, Laboratoire Roberval
UMR 7337, Centre de Recherches de Royallieu, CS 60319, 60203 Compiègne Cedex, France
e-mail: shanggui.cai@utc.fr, ouahsine@utc.fr</footnote>

<footnote>** Aix Marseille Université, CNRS, M2P2 UMR 7340, 13451 Marseille, France</footnote>

<footnote>† Université de Strasbourg, CNRS, ICUBE UMR 7357, 67000 Strasbourg, France</footnote>

<footnote>††Direction Technique Eau, Mer et Fleuves (CEREMA)
134, Rue de Beauvais, CS 60039, 60280 Margny-lès-Compiègne, France</footnote>
on a regular Cartesian mesh which can be accomplished by highly efficient solvers. A recent review of the immersed boundary method can be found in [5, 6].

Ever since introduced by Peskin [7] for the simulation of blood flow in a heart with flexible valves, IB method has been extended to a great variety of applications, including civil engineering, aeronautics, hydrodynamics, etc. The original IB method is categorized as “continuous forcing” in [5] as it is implemented into the Navier-Stokes equations before discretized. However, when applied to rigid boundary problems, the original IB method appeals to a spring with an extremely large stiffness to pull the boundary back to its equilibrium position [8], which in general gives rise to numerical difficulties. To overcome this shortness, Mohd-Yusof [9] proposed the direct forcing IB method using a discrete-time derivation of the forcing value to eliminate the spring, which falls in the “discrete forcing” group. Later Uhlamn [10] suggested evaluating the force on the Lagrangian coordinates instead of the Eulerian location, which is much simpler to perform. Nonetheless, those IB methods are all based on an explicit scheme, the no-slip boundary condition on the interface is never exactly satisfied even for the steady flow [11]. However, full implicit scheme is too cumbersome to perform, we need to iterate the whole system to ensure the incompressibility and no-slip boundary condition on the interface in the same time. In the present study we adopt the modified implicit IB method in [12] for accuracy and efficiency.

Almost all IB methods are able to deal with elastic and rigid boundaries. However, in the present work we only focus on the direct forcing IB method with rigid boundaries. Actually the simulation of lid-driven cavity flow with an embedded rigid cylinder is not a difficult task for the traditional body-conforming mesh method. One can mesh the domain with one O-type structured mesh near the cylinder. However the purpose of this paper is mainly to compare the numerical accuracy between the IB method and the traditional one.

In the next section, we discuss the governing equations and its numerical solution by the rotational incremental pressure correction projection method. In the following the implementation of the modified implicit IB method into the projection method is recalled. Then numerical results of lid-driven cavity flow with an embedded cylinder are compared between IB method and traditional method as well as the result in the literature. Conclusions are drawn in the final section.

# 2 METHODOLOGY

## 2.1 Governing equations and projection method

The unsteady viscous incompressible Navier-Stokes equations read

$$ \frac{\partial u_i}{\partial x_i} = 0 \qquad (1) $$

$$ \frac{\partial u_i}{\partial t} + \frac{\partial(u_i u_j)}{\partial x_j} = -\frac{1}{\rho} \frac{\partial p}{\partial x_i} + \nu \frac{\partial}{\partial x_j} \left( \frac{\partial u_i}{\partial x_j} \right) \qquad (2) $$
where $u_i$ are the velocity components, $p$ is the pressure, $\rho$ is the fluid density and $\nu$ is the fluid kinematic viscosity. The equations are integrated in time using the rotational incremental pressure correction projection method of Guermond et al. [13] which is proved to be second order in time for the velocity fields and has a consistent boundary condition for the pressure. The projection method consists of three sub-steps: prediction, projection and correction.

1. Predict an intermediate velocity $u_i^*$ with the pressure value at past time level.

$$ \frac{u_i^* - u_i^n}{\Delta t} + \frac{1}{2}[3C_i^n - C_i^{n-1}] = -\frac{1}{\rho}\frac{\delta p^n}{\delta x_i} + \frac{\nu}{2}(D_i^* + D_i^n) \quad (3) $$

where $C_i = \frac{\delta(u_i u_j)}{\delta x_j}$ and $D_i = \nu \frac{\delta}{\delta x_j} (\frac{\delta u_i}{\delta x_j})$ represent the convective and diffusive terms, respectively. The convective term is treated explicitly in time for simplicity by using a second order Adams-Bashforth scheme, and the diffusive term is discretized with an implicit Crank-Nicolson scheme to eliminate the viscous stability constraint. Therefore, the total scheme is only subjected to the standard CFL condition. The diffusion term is spatially discretized by a second-order central difference operator while the convective term is treated with the TVD scheme described in [14].

2. Project the tentative velocity $u_i^*$ to its divergence-free part $u_i^{n+1}$, by solving the pressure Poisson equation in conjunction with the incompressibility constraint.

$$ \frac{u_i^{n+1} - u_i^*}{\Delta t} = -\frac{1}{\rho} \frac{\delta \phi^{n+1}}{\delta x_i} \quad (4) $$

Taking divergence of Eq. (4) and applying Eq. (1), we obtain

$$ \frac{1}{\rho} \frac{\delta}{\delta x_i} \left( \frac{\delta \phi^{n+1}}{\delta x_i} \right) = \frac{1}{\Delta t} \frac{\delta u_i^*}{\delta x_i} \quad (5) $$

3. Correct or update the velocity and pressure field. Once the pressure correction $\phi^{n+1}$ is known, the final pressure and velocity can be updated as

$$ u_i^{n+1} = u_i^* - \Delta t \frac{1}{\rho} \frac{\delta \phi^{n+1}}{\delta x_i} \quad (6) $$

$$ p^{n+1} = p^n + \phi^{n+1} - \frac{\mu}{2} \frac{\delta u_i^*}{\delta x_i} \quad (7) $$

where $\mu$ is the fluid dynamic viscosity. In Eq. (7) the last term is actually the missing term due to time splitting when implicit scheme is applied to the viscous term in Eq. (3), which is now absorbed into the pressure.
<figure><img src="image_2.png" /><figcaption>Figure 1: Immersed boundary configuration.</figcaption></figure>

## 2.2 Immersed boundary method implementation

The principal idea of immersed boundary method is to represent the structure or solid boundary immersed in a fluid by an artificial force. The force is then added to the momentum equation as a source term to take effect. In fact, how the force is formulated determines each kind of IB method. In the present work, we choose the modified implicit IB method in [12], which can be briefed as

1. After the prediction step Eq. (3), interpolate the tentative velocity from the relying Cartesian mesh (see Fig. 1) to the structure or solid boundary $Γ_b: U_i^* = T(u_i^*)$, where upper case letter represents the quantities on the Lagrangian locations and $T$ is the field transfer operator.

2. Evaluate the force on the immersed boundary: $F_i = (U_i^d - U_i^*)/\Delta t$, where $U_i^d$ is the desired velocity boundary condition we want to impose on $Γ_b$. Here $i$ is the vector component index ($i = 1, 2$ in 2D, $i = 1, 2, 3$ in 3D) and not a specific Lagrangian point.

3. Spread the force to the relying mesh: $f_i^* = S(F_i^*)$, where $S$ is the spreading operator. In Fig. 1 the “•” donates the neighbouring points associated with the Lagrangian variable. The number of points depend on the support width of discrete delta function. In the present work, Peskin’s 4-point version is employed [8].

$$ \delta_h(r) = \begin{cases} \frac{1}{8}(3 - 2|r| + \sqrt{1 + 4|r| - 4r^2}), & |r| \le 1, \\ \frac{1}{8}(5 - 2|r| - \sqrt{-7 + 12|r| - 4r^2}), & 1 \le |r| \le 2, \\ 0, & \text{otherwise.} \end{cases} \quad (8) $$

4. Add the force to Eq. (4) to mimic the presence of structure or solid.

5. Solve the pressure correction Eq. (5) with the artificial force and update the velocity Eq. (6).

Repeat 1-5 until the desired velocity boundary condition is satisfied on these Lagrangian
points. Noting that the pressure correction equation is solved only once at each iteration for efficiency, we propose to perform Eq. (5) after the forcing loop to obtain better convergence and then update the velocity and pressure with Eq. (6), (7). In the present paper, the Krylov iterate solver of conjugate gradient along with the multigrid preconditioner is applied to solving the pressure Poisson equation.

## 3 NUMERICAL RESULTS

The flow of lid-driven cavity with an embedded cylinder is investigated in the present work to assess the accuracy of the improved implicit immersed boundary method. The length and height of cavity are unity ($L = 1$). The top wall is moving with a constant velocity $u_\infty = 1$ while the other three are stationary no slip walls. The cylinder is fixed at the center of the cavity with a diameter $D = 0.4L$, as shown in Fig. 2a for the IB method with uniform mesh arrangement. For comparison, a body-conforming mesh is applied in Fig. 2b. A fine mesh is used here, where the ratio of cylinder diameter to the mesh size is $D/h = 200$ in both cases. The flow pattern is characterized by the Reynolds number $Re = u_\infty L/\nu$, which is set to 1000 in the this study.

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Uniform mesh arrangement for the IB method (the embedded cylinder is represented by the thick dashed line), (b) body-conforming mesh arrangement.</figcaption></figure>

The convergence criteria is set to $||u_i^{n+1} - u_i^n||_2/||u_i^{n+1}||_2 < 1.0E-8$ so that the flow reaches the steady state. Fig. 3 shows the vorticity contour and streamline of the flow at $Re = 1000$ for both IB method and traditional body-conforming mesh method, where the vorticity contour value is varied from -3 to 3 with an increment of 0.4. Fig. 3a, 3c correspond to the IB method for the vorticity contour and streamline, while Fig. 3b, 3d are the results from the traditional method for vorticity and streamline respectively. The numerical simulation presented in this study is very similar to the results of Vanella and Balaras [15]. Apparently, three vortices emerge: one at the upper right near the cylinder and two lower at the left and right corners. It is noteworthy that the upper vortex is
<figure><img src="image_4.png" /><figcaption>Figure 3: Vorticity contour and streamline of lid-driven cavity flow with an embedded cylinder.</figcaption></figure>

generated by the presence of the fixed cylinder. The locations of three vortices center of both methods are listed in Table 1. The two methods yield nearly the same results.

<table><caption>Table 1: Comparison of vortices center positions</caption><thead><tr><th></th><th>upper right vortex</th><th>lower left vortex</th><th>lower right vortex</th></tr></thead><tbody><tr><td>IB method</td><td>(0.6946, 0.6872)</td><td>(0.0790, 0.0721)</td><td>(0.8856, 0.1061)</td></tr><tr><td>Traditional method</td><td>(0.6906, 0.6872)</td><td>(0.0791, 0.0721)</td><td>(0.8849, 0.1063)</td></tr></tbody></table>

The velocity component *u* at the vertical midline and the velocity component *v* at the horizontal midline are plotted in Fig. 4. Since the modified implicit IB method imposes the no-slip boundary condition on the interface exactly, the velocity profiles of both methods match very well.
<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of velocity profiles (a) distribution of velocity component u along vertical midline, (b) distribution of velocity component v along horizontal midline.</figcaption></figure>

## 4 SUMMARY

We have presented the numerical simulation of the lid-driven cavity flow with an embedded cylinder by the improved implicit immersed boundary method. Comparison has been made between the IB method, traditional body-conforming mesh method and the results in the literature. This validation case has demonstrated good agreements for both methods. Even though the fixed cylinder in the present work is not complex, the IB method can be easily extended to the flows with more complicated rigid or deforming, and/or moving geometries. In particular for the fluid-structure interaction problems with multi-structure in artificial shapes, the traditional body-conforming mesh method becomes extremely difficult and time-consuming to ensure the mesh quality especially after re-meshing. This task with IB method will be included in the future work.

## ACKNOWLEDGMENT

The first author acknowledges the financial support of the China Scholarship Council. The calculation is performed on the HPC platform PILCAM2 at the Université de Technologie de Compiègne.

## REFERENCES

[1] Ji, S.C., Ouahsine, A., Smaoui and H., Sergent, P. 3-D numerical simulation of convoy-generated waves in a restricted waterway. *J. Hydrodyn. Ser. B* (2012) **24**:420–429.

[2] Kaidi, S., Rouainia, M. and Ouahsine, A. Stability of breakwaters under hydrodynamic loading using a coupled DDA/FEM approach. *Ocean Eng.* (2012) **55**:62–70.
[3] Favier, J., Revell, A. and Pinelli, A. A lattice-Boltzmann-immersed boundary method to simulate the fluid interaction with moving and slender flexible objects. *J. Comput. Phys.* (2014) **261**:145–161.

[4] Darbani, M., Ouahsine, A., Villon, P., Naceur, H. and Smaoui, H. Meshless method for shallow water equations with free surface flow. *Appl. Math. Comput.* (2011) **217**:5113–5124.

[5] Mittal, R. and Iaccarino, G. Immersed boundary method. *Annu. Rev. Fluid Mech.* (2005) **37**:239–261.

[6] Sotiropoulos, F. and Yang, X. Immersed boundary methods for simulating fluid-structure interaction. *Prog. Aerosp. Sci.* (2014) **65**:1–21.

[7] Peskin, C.S. Flow patterns around heart valves: A numerical method. *J. Comput. Phys.* (1972) **10**:252–271.

[8] Lai, M.-C. and Peskin, C.S. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. *J. Comput. Phys.* (2000) **160**:705–719.

[9] Mohd-Yosuf, J. Combined immersed boundary/b-spline methods for simulation of flow in complex geometries. *Annu. Res. Briefs, Cent. Turbul. Res.* (1997)

[10] Uhlmann, M. An immersed boundary method with direct forcing for the simulation of particulate flows. *J. Comput. Phys.* (2005) **209**:448–476.

[11] Kempe, T. and Fröhlich, J. An improved immersed boundary method with direct forcing for the simulation of particle laden flows. *J. Comput. Phys.* (2012) **231**:3663–3684.

[12] Cai, S.-G., Ouahsine, A., Smaoui, H., Favier, J. and Hoarau, Y. An efficient implicit direct forcing immersed boundary method for incompressible flows. *Proc. of the 3rd International Conference on Mathematical Modelling in Physical Sciences (IC-MSQUARE)*. Madrid, Spain. (2014) **574**:012165.

[13] Guermond, J.L., Minev, P. and Jie, S. An overview of projection methods for incompressible flows. *Comput. Methods Appl. Mech. Engrg.* (2006) **195**:6011–6045.

[14] Smaoui, H., Zouhri, L. and Ouahsine, A. Flux-limiting techniques for simulation of pollutant transport in porous media: Application to groundwater management. *Math. Comput. Model.* (2008) **47**:47–59.

[15] Vanella, M. and Balaras, E. A moving-least-squares reconstruction for embedded-boundary formulations. *J. Comput. Phys.* (2009) **228**:6617–6628.