Sulzer White Paper 1 / 2018

A versatile immersed boundary method for pump design

Mathieu Specklin Ireland
Dipl.-Ing Mathieu Specklin, Research Engineer

"The use of computational fluid dynamics (CFD) calculations for the hydraulic design of conventional centrifugal pump impellers and single-phase liquids is well established. Sulzer initiated a research collaboration with the Dublin City University to work further on a new approach called the immersed boundary method (IB method).

The IB method can be used to solve a single set of conservation equations for several components, like gas, liquids, rigid solids or moving material. There is one main advantage of the IB method over the typical body fitted methods: With the IB method, we can model the interaction of fluids with rigid or moving bodies without the need for a mesh which conforms the interphase boundaries. As a result, it takes less effort and time to set up calculation models and the mesh points for complicated geometries. The accuracy of the predictions with the new method has been confirmed by extensive validation tests. In the future, such models will be able to accelerate the process of wastewater pump design.

The proposed numerical method has also been extended to model the transport of flexible cloth-like structure, which will provide insights into the process of pump clogging." 


Dipl.-Ing Mathieu Specklin, Postgraduate Researcher, Sulzer Pump Solutions Ireland Ltd., Wexford, Ireland


Robert Connolly, Hydraulic Design Engineer, Sulzer Pump Solutions Ireland Ltd.
Ben Breen, Product Development Manager, Sulzer Pump Solutions Ireland Ltd.
Dr. Yan Delaure, Lecturer, Dublin City University, Dublin, Ireland



One of the major challenges in CFD simulation lies in the creation of adequate computational meshes in realistic engineering configurations. The issues arise when dealing with complex, sharp and moving boundaries, removing the possibility of automatic creation of a high quality structured mesh. For moving boundaries or boundary deformations, in addition, the Generalized Grid Interface (GGI) methods used in these contexts are limited to moderate changes, and still require mesh morphing and re-meshing to the detriment of the computational time.

In this framework, Immersed Boundary types of methods (IB methods) are a promising alternative for this sort of problem. The aim of the IB method is to solve a single set of conservation equations for several components (gas, liquids or solids). As a result, this type of method makes it possible to handle multi-phase flow or fluid-structure interaction problems, in a single-phase framework. One of the main interests lies in the simplification of the model, and especially in its suitability for use with Cartesian meshes of perfect orthogonal quality.

The IB method developed in this project has been implemented in the open-source code OpenFOAM. The method is based on a penalty approach and is able to handle rigid moving solids with any kind of geometry, and with sharp and high definition of the solid interface. The formulation of the IB method is introduced along with an overview of second order corrections for velocity and pressure at the interface. The turbulence modelling choices are also presented.

The results of two benchmark test cases are presented, showing the ability of the model to handle both fixed and moving bodies with the expected accuracy. Finally, the results of the IB method for a full engineering case show good agreement with standard CFD techniques and experimental data.

Description of the problem

The mesh generation phase in a Computational Fluid Dynamics study of complex turbomachine applications can be a very significant challenge, which can typically absorb more than 50% of the R&D effort. Achieving a high quality mesh can be particularly troublesome when a lack of periodicity, symmetry and/or clear block structure makes automated structured grid mesh generation by specialised pre-processors very difficult or simply unavailable. Centrifugal pumps with a three-dimensional one-channel impeller, unless based on a simplistic geometry typically fall in this category of difficult geometry.

Fig. 1 illustrates the complexity of this kind of geometry. One can distinguish on this figure the volute and both inlet and outlet pipes of a centrifugal pump, as well as the one-channel impeller, which is driving the flow.

CAD Picture of a pump part
Fig. 1 Section of the CAD representation of a centrifugal pump with its one-channel impeller.

One solution in this case is to rely on hybrid meshing to combine regular hexahedral cells where possible, prism layers near solid surfaces where boundary layer must be accurately measured and other polyhedral cells to facilitate transition at interfaces. It can be difficult to avoid distorted mesh with associated numerical instabilities and loss of accuracy. This is particularly the case when strong gradients in mesh sizes cannot be avoided due to large cell count and high refinements near walls in high Reynolds flow or when interfaces between rotating and static zones in Generalised Grid Interface (GGI) models must be located in narrow passages. Also it should be noted that although polyhedral mesh blocks provide useful flexibility when hexahedral cells are used in the core of the flow domain and can greatly simplify meshing, the resulting hybrid mesh structure can be a problem in turbulence modelling with Large Eddy Simulation (LES) or Hybrid LES and Reynolds Averaged Navier Stokes (RANS) simulations, where the eddy viscosity model typically relies on a coefficient whose calibration is related to the mesh type and size. A uniform cubic mesh as achieved by a ‘Cut Cell’ type in this case provides the ideal solution where the characteristic size Δ can be defined exactly Δ = (Δ x Δ y Δ z)1⁄3.

Numerical solution

The numerical solution proposed in this work is based on an immersed boundary type of approach, which removes all difficulties related to the mesh, by solving the flow problem with moving boundaries on a Cartesian grid of perfect orthogonal quality. Internal boundaries from immersed solid (or flexible) surfaces are simulated using marker functions, which can be configured to identify Eulerian cells as internal or external to the fluid or solid domains. The solution is conceptually simple and varying degrees of sophistication have been introduced to improve the accuracy of boundary flow modelling. The new solution assessed in this article is based on a second order correction method which accounts for both velocity and pressure conditions at immersed surfaces.

The model presented here is implemented in OpenFOAM, which uses a Finite Volume Method for the Navier-Stokes solution. The pressure-velocity coupling is achieved by using the Pressure Implicit with Splitting of Operators (PISO) algorithm [1].

Definition of the immersed solid

The whole domain, including the immersed solid, is discretized with a Cartesian Eulerian grid. In addition, a discrete surface mesh, referred as a Lagrangian grid, is used to represent the internal solid and its surface. This mesh can be simply and easily built from triangular or other polyhedral faces. The surface mesh is used to differentiate the cell centers of the Eulerian grid inside the solid domain, from those inside the fluid domain. In the solution presented in this article, this is achieved using an octree partition of space.

In each Eulerian cell in the vicinity of the surface, the interface normal is computed. In addition, a distance function is defined to estimate the distance to the interface, as an approximation of the distance between the Eulerian cell center to the closest Lagrangian point of the surface mesh. Both normals and distances are used in the IB method algorithm to satisfy the correct boundary condition for velocity and pressure at the interface.

Immersed boundary formulation: A penalty approach

The presence of a solid body immersed in a fluid is taken into account by a momentum source term in the conservation equation of momentum for the fluid, namely the immersed boundary force. The IB method presented in this paper is based on a penalty approach [2], and leads to the following modified momentum equation for incompressible fluid.
Immersed boundary formulation

In this approach, the immersed solid is considered as a porous media with a very small permeability. The characteristic function χ allows localizing the position of the solid. This method makes it possible to impose the fluid velocity inside the solid domain at a specific user-defined value.

For Cartesian grids however, the forcing term of the penalty approach leads to a stair-case definition of the solid geometry, as illustrated in Fig. 2 (left). Furthermore, the cell centers do not necessarily coincide with the Lagrangian solid points representing the immersed body.

Two main errors are introduced with a simple penalty approach. Firstly the interface orientation is not accounted for, and secondly, the Dirichlet condition for the velocity (no-slip condition) is not necessarily satisfied at the exact position of fluid-solid interface. Significant improvement to the method’s accuracy can be achieved by correcting the velocity in the vicinity of the interface. Since the interface position and orientation are easily estimated with a prescribed motion of a rigid body, the no-slip condition can be imposed at the exact interface position even when this does not coincide with the grid cell center or vertices. This requires that the velocity is interpolated taking into account both the surface normal and the cell to surface distance. However, penalization of the velocity only leads to a spurious pressure field inside the solid, which can introduce errors primarily through the pressure gradient term. This issue has rarely been taken into account in the literature for immersed boundary approaches [3], and, to the authors’ knowledge has not yet been addressed with a penalty approach. In this work, two different types of pressure correction are proposed. The first solves the pressure in both fluid and solid domains independently in the pressure equation step. This decoupling is applied on the cells faces, and leads again to a staircase definition of the solid. To avoid losing information on the interface orientation, the second method based on a penalization of the pressure, satisfying a zero-pressure gradient along the normal at the interface has also been tested. Although more stable, the first method has been shown to bring little improvements whereas significant benefits have been linked to the second correction.


The current IB method has been validated against several test cases, covering different types of geometries and conditions for the immersed solid. The results of the IB method were obtained with second order velocity and pressure corrections near the interface. Two test cases, widely studied in the literature, are presented, respectively a fixed cylinder in a cross flow, and an in-line oscillating cylinder in a fluid at rest. In order to assess the accuracy of the model, both integrated and local physical values are compared to numerical and experimental data available in the literature. In addition, for each test case, numerical solutions were also computed on equivalent body-fitted grids.

The IB method solution is obtained on a Cartesian grid. For the standard approach, the body-fitted grid is based on the same mesh size than the IB method grid, as illustrated in Fig. 2 (right). A so-called ‘Cut-Cell’ method is used to fit the Cartesian cells to the solid geometry. An additional layer of identically sized boundary conforming cells is added on the surface of the solid, without any inflations layers. This was implemented in an attempt to minimize differences between the IB method and the body-fitted approaches.

Different grids for CFD calculation
Fig. 2 Sketch of the grids used to resolve the flow around a cylindrical obstacle. (left) Cartesian grid used for the IB method approach, distinguishing both fluid (in blue) and solid (in red) domains. The white line represents the interface of the immersed solid. (right) Equivalent grid in a body-conforming approach.

A third academic test case was considered based on a Wannier flow in the Stokes regime. As the analytical solution is known in this case, a quantitative assessment of the model was possible. Then, the capability of the present IB method to handle sharp geometry was investigated with a simplified static impeller in a channel flow. These last two test cases are not presented here, but have shown that good agreement could be achieved between the IB method results and the equivalent body-fitted ones.

Fixed cylinder in a cross flow

In this test case, a fixed cylinder in a 2D cross flow with a Reynolds number of 40 is considered. The domain is rectangular, with a length of 40D, and a height of 30D, where D is the diameter of the cylinder. The cylinder is placed at 10D from the inlet, and 15D from both top and bottom side walls. A slip-condition is imposed on the latter two boundaries, while a velocity inlet and an outflow boundary conditions are imposed on the side boundaries.

Although integrated values as drag coefficients are not sufficient to assess the accuracy of a model, they are of interest for engineering calculation and have been widely used in the literature to assess new models’ accuracy. One can note that for an immersed boundary formulation, the estimation of the forces applied on the immersed object can require additional and delicate computation, like reconstruction of the solid surface. In addition to the drag, the length of the recirculation zone behind the cylinder has also been investigated.

Both values are listed in table 1, for the body-fitted and the IB method approaches (with velocity and pressure correction), alongside measurements and computations from published studies.

Model CD LS
Coutanceau and Bouard (Exp.) [4] / 2.13
Tritton (Exp.) [5] 1.59 /
Range in Lit. Num. [6, 7, 8] 1.54 to 1.66 2.21 to 2.30
Body-fitted 1.60 2.31
IBM 1.60 2.30
Table 1 Values of the drag coefficient CD and recirculation length LS, obtained with a grid cell size Δ = D ⁄ 40 (refinement equivalent to those found in literature).

The low Reynolds flow considered here means that both the boundary layer over the cylinder and the wake remain laminar. The numerical solution provided here does not include any turbulence model and the mesh resolution is sufficient to resolve the boundary layer with 34 cells before detachment, with the finest grid (Δ = D ⁄ 40). The results obtained with the IB method in this case are in close agreement with the body-fitted results, and both drag coefficient and recirculation lengths are within the ranges found in the literature.

In-line oscillating cylinder in a fluid at rest

A cylinder oscillating in a fluid at rest is considered here. The dimensions of the problem are similar to the previous test case. For the body-fitted simulation, boundary motion is transfered to the cells through Laplacian smoothing which introduce a time dependent stretching.

The equation of motion is defined by x(t) = -Asin(2πft), where the frequency f is set to 0.2Hz and the amplitude A is fixed to 5/2π. This results in a Reynolds number of 100 and a Keulegan-Carpenter number of 5. The present validation case is also frequently used in literature for numerical models handling moving geometries, and laboratory experiments have also been conducted [9].

The oscillations of the cylinder lead to the development of lower and upper boundary layers, which separate and generate two counter-rotating vortices. This generation process stops when the cylinder starts to move backward and finally splits the vortex pair. The velocity contours obtained with our IB method are plotted in Fig. 3 at a time t = 15 s, which corresponds to a phase angle of 0°, when the cylinder just ends its translation and starts to move backward. This velocity contours show that the IB method preserves symmetry, and are in good agreement with the results of Dutsch. Furthermore, the averaging of the drag coefficients over half a period gives similar results for both body-fitted and IB method approaches on the fine grid (40 cells in diameter). Finally, although a detailed assessment of the different stages of corrections is beyond the scope of this article, Fig. 3 illustrates the differences between a first order velocity penalisation applied at cell centers without reference to the relative surface position and orientation and a 2nd order penalisation. The latter method predicts a velocity field inside the solid which is not uniform but interpolated instead.

Comparison of contours of velocity magnitude and streamlines
Fig 3 Contours of velocity magnitude and streamlines (white lines) around the oscillating cylinder at t = 15 s (phase angle of 0°), obtained with a simple penalty approach (left), and a penalty approach with velocity correction (right). The black line represents the cylinder interface, while the vertical red line shows the position of the extracted velocity profiles.

Fig. 4 shows the streamwise and transverse velocity profiles across the center of the cylinder at the same time t = 15 s. IB method and body-fitted approaches are in very good agreement, and match well with the experimental results of Dutsch et al. In general, the differences between both numerical approaches used here are below 15% at the maximum. Previous IB method in the literature do not compare in general to body-fitted simulation, but a similar level of accuracy has been found for our IB method in comparison to the experimental data.

A full engineering case: the sewage pump

A centrifugal sewage pump has been simulated at Best Efficiency point and at part loads. The special feature of this single-stage pump is to drive the flow with a single blade impeller. This type of impeller is able to handle liquids with stringy materials and large solids, without screening while maintaining good hydraulic performance.

Comparison of velocity profiles
Fig. 4 Velocity profiles across the cylinder at t = 15 s (phase angle of 0°).

For the IB method solution, a mesh based on a regular Cartesian grid snapped to the inlet sump and volute boundaries is generated to discretize the whole fluid domain, including the zone swept by the impeller. The presence of the impeller is accounted for as a surface mesh. Removing the impeller boundaries from the grid drastically simplifies the meshing stage, and allows the use of cubic cells everywhere except at wall boundaries. Both Eulerian and Lagrangian meshes for the fluid and impeller are shown in Fig. 5.

The cell size inside the volute is Δ = 2 mm, which is also the average size of the triangular cells for the surface mesh. This setting leads to a Eulerian mesh with about 2 millions cells in total. Here again, the IB method is compared to a standard body-fitted solution. In this case, a grid with similar cells dimensions has been created but in this case the impeller volute has been subtracted from the flow domain and a single prism layer has been inserted at the impeller surface. The resulting mean and max y+ for cells adjacent to the impeller surface in the body fitted mesh are respectively around 230 and 880.

Different meshes for CFD
Fig. 5 (left) Cross-section of the computational domain and mesh. The cells inside the fluid are coloured in cyan, while the cells inside the immersed body representing the impeller are coloured in red; (right) Surface mesh of the one-channel impeller.

This would typically be considered insufficient if a resolution of the boundary layer was sought. It should be noted however that the primary purpose here is to compare body fitted and IB method simulation. In order to model the rotation of the impeller with the body fitted solution, the Generalised Grid Interface (GGI) method is used. With this technique, two different mesh parts are defined, one for the fixed background, and one for the moving zone which includes the impeller surface.

The flow is fully turbulent in all cases and the IB method was extended to include a SA-DES turbulence model. This is part of the family of hybrid RANS-LES models [10] and is based on the Spalart-Allmaras (SA) model for the RANS part. The standard transport equation for the turbulent viscosity has been modified in order to satisfy a zero Dirichlet condition at the immersed interface, as it is assumed in the SA model [11]. This modification has been implemented based on the same penalization approach as for velocity in the momentum equation.

The Detached-Eddy Simulation (DES) model, which uses the distance to the wall as a trigger for the transition between RANS and LES approaches [12], has been combined with our IB method by taking into account the presence of the immersed body to estimate this distance. Pure RANS simulations have been found to provide broadly similar results to the DES models for the pump head in particular at BEP but produce excessive turbulent dissipation failing to reproduce the rich vortical structures that can be expected for this type of high Reynolds flow.

The LES simulations do resolve this but their excessive mesh requirements at wall boundaries make it an impractical solution for transient flow simulations (for example for modelling of multiphase flow). The DES models take advantage of the RANS modelling for wall layers with lower mesh refinement requirements, while switching to LES modelling outside of this region enabling better modelling of turbulence scales up to the filter size. As for any RANS model, relevant wall functions should be used at walls when the grid refinement is not sufficient to resolve the viscous sublayer. This has not been implemented yet and accurate modelling of wall effects would require mesh refinement.

Formula for head calculation

The total losses can be estimated from the hydraulic power delivered by the pump, which is a function of the torque. In this work however, the losses have been subtracted from the hydraulic head, whose calculation does thus only depend on the difference of static pressure and elevation. In addition, the dynamic head has been omitted as the two monitor surfaces used have the same cross section. The evolution of this head for a flow rate of 40 l/s is shown in Fig. 6 (left) for both IB method and body-fitted approaches with DES modelling for turbulence. Although it appears that the IB method slightly overestimates the head in comparison to the body-fitted solution, the difference between the two methods is lower than 5%.

The contours of velocity magnitude predicted by the IB method and body fitted approaches over a horizontal section plane are compared in Fig. 6 (right) in the same conditions. This confirms the similarity in flow structures in particular on the concave part of the impeller. The average size of the vortices is similar for both methods. One can note a slightly higher velocity magnitude in the discharge side for the IB method. In the suction side, although the IB method seems to generate too many small vortices, the flow features are well reproduced. The corrections applied in the IB method solution near the impeller interface were found to improve the comparison, especially when the immersed solid becomes thin (below 4 cells to resolve the thickness), as it is the case at the trailing edge.

Comparison of CFD methods
Fig. 6 Comparison of body-fitted and IB method approaches. (left) Evolution of the pump heads; (right) Velocity fields in a horizontal section plane of the pump (located 8 mm below the outlet centre point). For the IB method approach (right), the impeller interface is shown in white.
Fig. 7 shows the shear strain rate obtained with the IB method in the centre part of the volute, with a DES turbulence model and for a flow rate of 40 l/s. The rich elongated vortical structures are well captured with the DES model. The results are qualitatively similar for both IB method and body-fitted approaches. However, regular rolls along the surfaces of the impeller predicted by the IB method are not observed with standard body fitted computations. This suggests that the IB method boundary treatment, although accurate with laminar flow past a cylinder, may produce local variations in this type of pump flow.
Shear strain rate in a pump
Fig 7 Shear strain rate Q=100000 (coloured by the magnitude of the velocity) obtained with IB method for a flow rate of 40 l/s.
A performance study has been performed on the pump XFP 150E. The hydraulic head results obtained experimentally as well as the numerical heads predicted with the IB method and body-fitted approaches are summarised in Fig. 8 for five different flow rates. The IB method method is providing estimates which are in very good agreement with the experimental data around the best efficiency point (55 l/s). An underestimation of the head tends however to appear with the appearance of large separation and the formation of large trailing vortices. Predictions in this case can be expected to be strongly influenced by accurate modelling of boundary layers and would require much finer wall adjacent mesh. Adaptive meshing could provide a practical solution in this case without seriously compromising computational efficiency.
Comparison of pump heads
Fig. 8 Non-dimensional experimental and numerical heads for different flow rates.


The research presented in this article concerns the validation of a new IB method developed as an alternative to body-fitted simulation based on the GGI method. The method is intended for the simulation of turbomachine applications which are not suitable for automated meshing tools designed to exploit block structures and periodicity in domain geometries and capable of generating a high quality mesh. IB method brings a number of significant additional benefits. It can make use of near perfect cubic mesh in the core of the flow domain removing issues in LES simulations related to mesh size and shape variations. It can model solid bodies with prescribed motion but also fluid structure interaction with flexible bodies. This includes the transport of flexible cloth-like structures as found in urban waste water processing and of particular interest to ongoing research.

Finally, it avoids numerical issues related to the presence of over-hanging nodes in non-conformal mesh interfaces found in GGI methods. This introduces a number of approximations for example when dealing with the transport of multiphase flows. The computational results presented here confirm that the IB method can produce integrated as well as local velocity and pressure values that are in close agreement with experimental data and standard body-fitted CFD simulations. The validation covers a broad range of flow conditions, including low and high Reynolds flows. Both attached and separated flow past fixed and oscillating obstacles have been considered. The IB method is shown to be a promising tool and is undergoing further development to allow for the modelling of flexible and deformable solids and single- and two-phase flow. This will support ongoing research in the study of processes leading to pump clogging.

The financial support of the Irish Research Council (IRC) under its Employment Based Research Scholarship programme (Grant No. EBPPG/2013/63) is acknowledged.


[1] R. Issa, R. I. (1986). Solution of the implicitly discretized fluid flow equations by operatorsplitting. Journal of Computational Physics, 62, 66 - 82.

[2] J. Caltagirone, E. A. (1986). Recirculating flow in pourous media. C. R. Aca. Sci. Paris, 302, 843 - 846.

[3] S. Kang, G. I. (2009). Prediction of wall-pressure fluctuation in turbulent flows with an immersed boundary method. Journal of Computational Physics, 228, 6753 - 6772.

[4] M. Coutanceau, R. B. (1977). Experimental determination of the main features of the viscous flow in the wake of a circlar cylinder in uniform translation. Part 1: Steady flow. J. Fluid. Mech., 79, 231 - 256.

[5] D.J. Tritton, (1959). Experiments on the flow past a circular cylinder at low Reynolds number. J. Fluid. Mech., 6, 547 - 567.

[6] K. Taira, T. C. (2007). The immersed boundary method: A projection approach. Journal of Computational Physics, 225, 2118 - 2137.

[7] L. Parussini, V. P. (2009). Fictitious domain approach with hp-finite element approximation for incompressible fluid flow. Journal of Computational Physics, 228, 3891 - 3910.

[8] S. Xu, Z. J. (2005). An immersed interface method for simulating the interaction of a fluid with moving boundaries. Journal of Computational Physics.

[9] H. Dutsch, F. D. (1998). Low Reynolds number flow around an oscillating circular cylinder at low Carpenter numbers. J. Fluid Mech., 360, 249 - 271.

[10] S-H. Peng, P. D. (2009). Progress in Hybrid RANS-LES Modelling. Berlin: Springer.

[11] S. R. Allmaras, F. T. (2012). Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model. Proceedings of ICCFD7.

[12] P. R. Spalart, W.-H. J. (1997). Comments on the feasibility of LES for wings and on a hybrid RANS/LES approach. Proceedins of the first AFOSR International Conference on DNS/LES.

Sulzer Technical Review


Sulzer Management AG

Neuwiesenstrasse 15

8401 Winterthur