Immersed Boundary Method

Zongxian Liang, Developer and Adwaith Gupta, CFD Engineer at Flow Science, contributed this blog.

The accurate prediction of forces and energy losses is crucial in the success of modeling many engineering problems such as discharge from an orifice plate, flow past an obstacle and flow in a sudden contraction pipe. The upcoming release of FLOW3D v12.0, features a new numerical option, the Immersed Boundary Method, that accurately predicts flows for such problems.

Immersed Boundary Method, or precisely, ghost-cell based Immersed Boundary Method, improves the accuracy of numerical flux calculations at the solid-fluid interfaces. In this blog, Zongxian Liang, developer at Flow Science provides validation examples of a sudden contraction pipe and a ship hull where the Immersed Boundary Method improves the accuracy of head loss and resistance force. Brief mathematical details of the Immersed Boundary Method and the ghost cell approach are discussed at the end of the blog.

Sudden contraction pipe

One of the fluid problems that demonstrate the accuracy of the immersed boundary solver is estimating the head loss of water in a contraction pipe, as shown in Figure 1. The pipe experiences a sudden contraction from a large section of 3m in diameter to a small section of 1m. Flow rate at the inlet of the large pipe is 4m3/sec. Velocity heads are measured 3.5m upstream and downstream with respect to the contraction location. Based on the different assumptions used in the literature [1,2], the theoretical value of the head loss should be in the range of 0.494m to 0.711m. In this simulation, one Cartesian mesh block is used for the whole geometry and a pressure boundary is specified with zero stagnation pressure at the exit. The two-equation k-ω model is enabled with maximum turbulent mixing length set to 0.07m. Gravity is active in the negative x direction. The finish time is set to 15s at which the flow is considered to be steady and fully developed.

Contraction pipe immersed boundary method
Figure 1. Schematic describing the shape of contraction pipe, location of flux surfaces for measuring velocity head, and flow direction.

A two-variable parametric study was conducted to investigate the effect of cell size and mesh alignment with respect to the sudden contraction location. The cell sizes of the mesh are 0.1m and 0.05m. Four mesh alignments were tested, as listed in Table 1: “Aligned” denotes the baseline case where the mesh plane is aligned with the sudden contraction location, and “Z-shifted: X%” indicates that the mesh plane is shifted an X% of the cell size away from the contraction location in the z-direction.

The predictions given by the immersed boundary solver all fall in the theoretical range of 0.494 to 0.711. Especially for the cases with a smaller cell size of 0.05m, the head losses are within 4.6% of the median of the theoretical values, 0.603.

    Grid Alignment
  Aligned Z-shifted: 25% Z-shifted: 50% Z-shifted: 75%
Cell size (m) 0.1 0.706 0.710 0.647 0.666
0.05 0.607 0.620 0.575 0.589

Table 1: Head loss (m) predicted by the immersed boundary solvers

Resistance force on a model of ship hull

The accurate prediction of forces in external flow dynamics problems is typically of interest in design phases. An example of an external flow problem with a free surface is the calculation of the total resistance force on a NAVY vessel model hull. In this case, the hull length is 5.72m, and the draft is 0.248m. Based on the hull length and mean flow velocity of 2.10m/s, the Reynolds number in the simulation is approximately 12×106. Because the problem is symmetric, only half of the hull is modeled. The computational domain is 30m long by 8m wide and 5.5m deep, and has three nested mesh blocks with the smallest cell size equal to 0.02m. Velocity and outflow conditions are used for the inlet and outlet boundaries, respectively. RNG turbulence model with dynamic turbulent length scale calculation is used to model the turbulent flow, and the second-order, monotonicity preserving scheme for the momentum advection.

Animation showing the pressure contours near the free surface and resistance force history.

The total resistance force is 22.62N for half of the hull model, based on a total drag coefficient of 0.0423 in experiments [3]. The force from the simulation is the sum of the shear and pressure forces in the x direction, averaged from 35 to 50s. The immersed boundary solver gives 22.43N, 0.8% lower than the experimental result.


We see from these two validations that the Immersed Boundary Method approach estimates the head loss and the force correctly within the theoretical range provided in the literature. Immersed Boundary Method improves the accuracy of estimating fluxes at a very fundamental level and is expected to improve simulation results for many, if not most, applications. To know more about the ghost-cell based Immersed Boundary Method development, continue reading. Otherwise, stay tuned for our next blog post!

Ghost-Cell Based Immersed Boundary Method

In FLOW-3D, the free-slip boundary condition is applied to the advection of velocities to eliminate numerical boundary layers caused by fractional cell areas and volumes near solid boundaries. In order to estimate a rational flux for the control volume, the immersed boundary solver computes a fluid velocity in the solid that implicitly satisfies the boundary condition, as illustrated in Figure 2. The fluid cell in the solid is called ghost cell and this method is generally referred to as the ghost-cell approach.

Immersed boundary method
Figure 2. Flux at the left face of control volume (enclosed by blue dash-dot lines) is computed using the velocity of the ghost cell u_(i-1) inside the solid.

To enforce the boundary condition, an image-point (open diamond denoted by IP) of the ghost cell is created in the fluid region by extending a line segment from the ghost cell (red diamond, denoted by GC) to the wall in the normal direction and that intersects the wall at the midpoint (open circle denoted by BI for boundary-intercept point) between GC and IP. For a non-penetrating boundary condition at BI and assuming the tangential velocities at GC and IP are the same as the wall surface velocity, the velocity of the ghost cell is estimated by the following equation:

(1)    \displaystyle {{\vec{u}}_{{GC}}}={{\vec{u}}_{{IP}}}-2\left[ {\left( {{{{\vec{u}}}_{{IP}}}-{{{\vec{u}}}_{{BI}}}} \right)\cdot \vec{n}} \right]\vec{n}

where \displaystyle {{\vec{u}}_{{GC}}}\displaystyle {{\vec{u}}_{{IP}}}  and  \displaystyle {{{{\vec{u}}}_{{BI}}}} are the fluid velocity at the ghost cell, the image-point and the boundary-intercept point, respectively, and   is the unit normal vector at the boundary. We use trilinear interpolation to approximate the value of velocity  using fluid velocities in the neighboring cells:

(2)    \displaystyle {{\vec{u}}_{{IP}}}=\sum\limits_{{i=1}}^{8}{{\left( {{{\alpha }_{i}}{{u}_{i}}} \right)}}\hat{i}+\left( {{{\beta }_{i}}{{v}_{i}}} \right)\hat{j}+\left( {{{\gamma }_{i}}{{w}_{i}}} \right)\hat{k}

where u1, v1 and w1 are the velocity of interpolation nodes surrounding the image-point and α1, β1 and γ1 are interpolation coefficients. The details of computing the interpolation coefficients can be found in Ref. 4.

The interpolant invokes the velocity value of other ghost cells, and a coupled system of the ghost cells is thus generated. We use a Jacobi based iteration method in conjunction with convergence acceleration techniques to obtain a fast solution of the resulting system.


  1. White, F. M., Fluid Mechanics (McGraw-Hill Book Company, 2003).
  2. Saleh, J., Fluid Flow Handbook (McGraw-Hill Professional, 2002).
  3. Larsson, L., Stern, F. & Bertram, V., Benchmarking of computational fluid dynamics for ship flows: the Gothenburg 2000 workshop. Journal of Ship Research 47 (1), 63–81 (2003).
  4. Mittal, R. et al., A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. Journal of computational physics 227 (10), 4825-4852 (2008).