You are viewing a javascript disabled version of the site. Please enable Javascript for this site to function properly.
Go to headerGo to navigationGo to searchGo to contentsGo to footer
In content section. Select this link to jump to navigation

Application of a maritime CFD code to a benchmark problem for non-Newtonian fluids: the flow around a sphere


The ship’s resistance and manoeuvrability in shallow waters can be adversely influenced by the presence of fluid mud layers on the seabed of ports and waterways. Fluid mud exhibits a complex non-Newtonian rheology that is often described using the Herschel–Bulkley model. The latter has been recently implemented in a maritime finite-volume CFD code to study the manoeuvrability of ships in the presence of muddy seabeds. In this paper, we explore the accuracy and robustness of the CFD code in simulating the flow of Herschel–Bulkley fluids, including power-law, Bingham and Newtonian fluids as particular cases. As a stepping stone towards the final maritime applications, the study is carried out on a classic benchmark problem in non-Newtonian fluid mechanics: the laminar flow around a sphere. The aim is to test the performance of the non-Newtonian solver before applying it to the more complex scenarios. Present results could also be used as reference data for future testing. Flow simulations are carried out at low Reynolds numbers in order to compare our results with an extensive collection of data from the literature. Results agree both qualitatively and quantitatively with literature. Difficulties in the convergence of the iterative solver emerged when simulating Bingham and Herschel–Bulkley flows. A simple change in the interpolation of the apparent viscosity has mitigated such difficulties. The results of this work, combined with our previous code verification exercises, suggest that the non-Newtonian solver works as intended and it can be thus employed on more complex applications.


The ship’s resistance and manoeuvrability in shallow waters can be adversely influenced by the presence of fluid mud layers on the seabed of ports and waterways. With today’s computers it is possible to investigate this problem using Computational Fluid Dynamics (CFD) codes [21,25]. However, modelling the complex non-Newtonian rheology of mud is a challenging task. For engineering purposes, on the other hand, mud can be reasonably well modelled as a Herschel–Bulkley fluid [12,46] (see e.g. the textbook of Irgens [24] for an extensive overview of non-Newtonian models). Hence, with the aim to investigate the effects of muddy seabeds on marine vessels, the Herschel–Bulkley model has been recently implemented in the CFD solver ReFRESCO [45], which is developed by the Maritime Research Institute of the Netherlands (MARIN) in collaboration with several non-profit organisations around the world.

Our previous code verifications [30,31] ensured that the Herschel–Bulkley model was correctly implemented. However, even if the code is correct, fully-converged solutions for realistic non-Newtonian problems may still be difficult to obtain. In fact, ReFRESCO is optimised and verified exclusively for maritime applications, which are typically concerned with Newtonian fluids such as air and water. The code also presents features that are common in general commercial CFD codes, such as the finite-volume method (FVM) and SIMPLE-type solution algorithms. While these features are standard for maritime applications, they are less common to simulate non-Newtonian flows. The latter are characterised by a shear-dependent viscosity that makes the equations stiffer and thus more difficult to solve.

In this paper, the accuracy and robustness of the non-Newtonian solver of ReFRESCO are tested on the laminar flows of Herschel–Bulkley fluids around a sphere. Simulations are carried out also for power-law and Bingham fluids, which are just particular cases of Herschel–Bulkley fluids.

The flow of power-law, Bingham and Herschel–Bulkley fluids around a sphere is a classic problem in non-Newtonian fluid mechanics. While it is probably the simplest three-dimensional flow, it also exhibits features that are typical of the flow around ships, such as boundary layer development and flow separation. These reasons make the flow around a sphere a useful benchmark problem as a stepping stone towards numerical simulations of ships navigating through fluid mud.

The non-Newtonian flow around a sphere has been extensively studied in the past decades, both experimentally (e.g [1,2,4,5,9,27,42,44]) and numerically (e.g [68,13,14,22,23,3335,43]), because of its interest for both industrial and natural processes such as sediment and fluid transport, sedimentation and fluidisation (the reader is also referred to the book of Chhabra [11] for an exhaustive survey on the topic). These processes are typically found in applications concerning offshore, dredging and ocean engineering.

The present work differs from the previous numerical studies in a number of ways. It is remarked that, although the flow field will be discussed, the main focus is on the comparison with literature data. Therefore, more attention is dedicated to the estimation of the uncertainties, which are of utmost importance for validation. The numerical uncertainties were estimated with a robust method and they are provided for the pressure, frictional and total drag for all the considered test cases. The uncertainties due to the regularisation of the Bingham and Herschel–Bulkley models have been quantified as well. The estimated uncertainties, together with the computed pressure and frictional coefficients reported herein, will also allow detailed and rigorous validations of alternative numerical methods in the future. Finally, this article also covers some issues related to the iterative convergence, which can often be problematic for non-Newtonian flow simulations.

The rest of the paper is divided as follows. The governing equations and test cases are outlined in Sections 2 and 3, respectively. The numerical methods and setup are explained in Section 4, whereas results are discussed in Section 5. Finally, the conclusions are summarised in Section 6.

2.Problem formulation and governing equations

The problem of a sphere d moving in an infinite medium is modelled as a sphere fixed at the centre of a cylindrical tube with a uniform inflow U (Fig. 1). The tube has both diameter and length equal to D. The origin of the Cartesian reference frame is placed at the sphere centre with the z-axis aligned with flow direction.

Fig. 1.

Schematic representation of the problem.

Schematic representation of the problem.

The laminar, steady, isothermal and single-phase flow of an incompressible fluid is governed by the following continuity and momentum equations (in Cartesian coordinates):

where u=(ux,uy,uz) is the velocity vector, p is the pressure, ρ is the density and τ is the deviatoric stress tensor. For a class of non-Newtonian fluids called generalised Newtonian fluids [24], the stress tensor can be expressed in a Newtonian-like form as
where γ˙=2SijSij is the shear rate (s1) and Sij=12(u+uT) is the deformation rate tensor. For Herschel–Bulkley fluids, the apparent viscosity, μ(γ˙), is modelled as
where τ0 is the yield stress (Pa), τ=τijτij/2, n is the flow index, k (Pa sn) is the consistency parameter which has dimensions of a viscosity when n=1. The infinite viscosity means that the fluid cannot deform (Sij=0) when the stress level is below the yield stress. This rheological characteristic is called viscoplasticity. The Herschel–Bulkley model reduces to the Bingham model when n=1, to the power-law model when τ0=0, and to the Newtonian model (μ=k) when both n=1 and τ0=0.

The issue related to the infinite viscosity in equation Eq. (3) when τ0>0 is avoided using regularisation methods (see e.g. [37] for an overview on the subject). These methods consist in approximating the original rheological model with a smooth function that produces large but finite values of the viscosity. In this work, the Bingham model has been modified with the popular Papanastasiou regularisation [36]:

where m (s) is the regularisation parameter.

For Herschel–Bulkley fluids, the Papanastasiou regularisation would still produce an infinite viscosity for γ˙0 when n<1. This is avoided using the modification proposed by Souza Mendes and Dutra [38]:

In summary, the apparent viscosity of Bingham fluids is calculated using Eq. (4), whereas Eq. (5) is used for Herschel–Bulkley fluids. In the limit of m, both regularised models tend to their respective ideal (non-regularised) models, as shown in Fig. 2 for Herschel–Bulkley fluids.

Fig. 2.

Effect of the regularisation parameter on the Herschel–Bulkley model.

Effect of the regularisation parameter on the Herschel–Bulkley model.

3.Non-dimensional numbers and test cases

The Herschel–Bulkley flow around a sphere is characterised by two non-dimensional numbers: the generalised Reynolds number

(6)Re=ρU2k(U/d)n(inertiaviscous stress),
and the Bingham number
(7)Bn=τ0k(U/d)n(yield stressviscous stress).
For power-law and Newtonian fluids the Bingham number is zero, whereas when n=1 the Reynolds number reduces to the canonical Re=ρUd/k.

The non-dimensional friction, pressure and drag coefficients are calculated respectively as

where S is the sphere surface having n as its outward normal vector, (Fp)z and (Ff)z are the z-component of the pressure and frictional force vectors, respectively.

Furthermore, throughout the paper, the regularisation parameter m will be expressed in its non-dimensional form,


Numerical simulations are performed for all the combinations of n=1,0.8 and 0.6, Bn=10 and 100 and Re=10 and 100. This choice is made to maximise the number of considered cases for which there is published CFD data in the literature while also running a feasible number of simulations.

In addition, simulations are performed for creep flow (Re<1) of Bingham fluids, for two reasons. First, under the latter flow regime, obtaining a fully converged solution with a SIMPLE-type solver can be challenging (if not impossible in some cases). Second, literature data for Bingham creep flow are in excellent agreement with each other, contrary to data for Re=10 and 100. For these two reasons, the Bingham creep flow is an interesting case to test the accuracy and robustness of the solution approach, even though it is far from the typical flow conditions encountered in maritime applications.

For Bingham creep flow, the drag coefficient as defined in Eq. (8) becomes extremely large. The reason is that inertia is virtually zero for creep flows, hence 1/2ρU2 is no longer a suitable factor to non-dimensionalise the forces. Since viscous effects are dominant, a better alternative could be to use kU/d. Nonetheless, we have adopted the more common practice of using the Stokes coefficient [7], defined as

The Stokes coefficient is thus a measure of how large is the drag force compared to the exact drag coefficient for Newtonian creep flow. The latter is equal to 24/Re according to the Stokes’ law [39].

4.Numerical methods and setup

4.1.Flow solver

The CFD code used for the present work is ReFRESCO, in which the governing equations are discretised with a second-order finite-volume method for unstructured mesh with cell-centered co-located variables. Mass conservation is ensured with a pressure-correction equation based on a SIMPLE-like algorithm [26] and a pressure-weighted interpolation technique to avoid spurious oscillations [32]. The convective term in the momentum equation is linearised with the Picard method and discretised with a central difference scheme. Other features such as cavitation and turbulence models (the latter also for Herschel–Bulkley fluids [29]) are also implemented but not used in the present work. These and other numerical techniques are described in detail elsewhere (e.g. [19]).

4.2.Boundary conditions

The boundary conditions are set as follows. A uniform inflow uz=U is imposed at the inlet (z=D/2, see also Fig. 1), whereas outflow conditions (uz/z=0) are imposed at the outlet boundary (z=D/2). The no-slip/no-permeability conditions (u=0) are set at the sphere surface while the free-slip conditions (un=0, uz/n=0) are applied to the outer cylindrical tube. For the pressure field, Neumann conditions (p/n=0) are applied to all the boundaries, thus a pressure reference is imposed in one point at the inlet boundary.

4.3.Domain size

In order to mimic a sphere settling in an infinite domain, the outer cylindrical tube diameter must be sufficiently large compared to the sphere diameter. The influence of the tube-to-sphere diameter ratio, D/d, was assessed by computing the drag coefficients on four grids having size D=25d, 50d, 100d and 200d respectively, with the grid 50d having the second finest refinement shown in Table 1 (Section 4.4) and the regularisation parameters in Table 3 (Section 4.5).

The uncertainty in the solution due to the domain size has been estimated with the method of Eça and Hoekstra [17] by replacing the grid size with the tube diameter D. Since this method was designed to estimate the discretisation uncertainties, we have checked the validity of this procedure by performing additional calculations for Newtonian creep flow, whose exact solution for the unbounded domain yields CS=1, by virtue of Eq. (11). With the adopted procedure the extrapolated CS for d/D=0 was found to be 0.99995 (Fig. 3, left), thus the procedure is deemed reliable for the present calculations. Note that the 4% uncertainty shown in Fig. 3 is irrelevant in this work as the Newtonian creep flow will not be considered further.

Fig. 3.

Convergence of the Stokes/drag coefficients with D1/Di (Di=200d,100d,50d,25d) for Newtonian creep flow (Re=0.01) (left) and Newtonian and power-law flows (Bn=0) at Re=10 (right). The percentages indicate the domain uncertainty for the case with D1=50d.

Convergence of the Stokes/drag coefficients with D1/Di (Di=200d,100d,50d,25d) for Newtonian creep flow (Re=0.01) (left) and Newtonian and power-law flows (Bn=0) at Re=10 (right). The percentages indicate the domain uncertainty for the case with D1=50d.

The uncertainties were then estimated for all the test cases, and the two largest values obtained with D=50d were 0.12% and 0.09%, which corresponded to the Newtonian and power-law cases (Bn=0) for Re=10, respectively (Fig. 3, right). This is not surprising since the strongest disturbance of the sphere to the flow field occurs for such test cases (see also Fig. 8 in Section 5.1). For Bn>0, the domain size uncertainties were found to be less than 0.005%. Therefore all subsequent calculations were performed with D=50d.

In conclusion, the maximum domain-size uncertainty for Newtonian and power-law fluid is about 0.1%, whereas for all the other test cases the influence of the domain size can be safely neglected.

4.4.Domain discretisation

The domain was discretised with two series of three-dimensional multi-block structured grids, with the layout as sketched in Fig. 4. Four geometrically similar grids were generated for each series in order to estimate the discretisation uncertainties. The number of grid cells and the size of the first cells away from the sphere surface are reported in Table 1. One series was generated for Re=0.01 (creep flow) of Bingham fluids and the other for the higher Reynolds numbers (Re=10 and 100). Note that the series used for creep flow calculations has higher resolution in order to better capture the steep velocity gradient near the sphere walls occurring for large Bn.

Fig. 4.

Illustration of the coarsest grid of Series 1.

Illustration of the coarsest grid of Series 1.
Table 1

Grid size ratio, number of cells Ni and thickness hw of the first cell away from the sphere surface

iSeries 1Series 2

Re=10 and 100Re=0.01

42.59129 0640.008002.60156 8080.00200
31.90324 0700.005711.91397 0720.00143
21.38857 6000.004081.381 054 2080.00102
11.002 238 3680.002921.002 761 0880.00073
Table 2

Discretisation uncertainties in percentage of the corresponding drag component for the pressure, frictional and total drag on the finest grid



The discretisation uncertainties on the force coefficients are given in Table 2 and they were estimated with the method of Eça and Hoekstra [17]. All the calculations were carried out using the regularisation parameters discussed in the following section. Except for the two cases with Bn=100 and n=0.6, all the uncertainties are below 1%. The larger uncertainties observed in CDf for Bn=100 and n=0.6 can be explained by the decrease of the apparent viscosity at the sphere surface (high-shear region) when n is decreased. In fact, the already steep velocity profile at the sphere walls for Bn=100 (see also Fig. 11 in Section 5.1) becomes even steeper when the viscosity is reduced, which would thus require a finer grid to be more accurately captured.

4.5.Regularisation methods

The use of regularisation methods produces regularisation errors, which are defined as the difference between the solution obtained with the regularised and non-regularised rheological models. Large regularisation parameters are needed to minimise these errors, but this may lead to very slow or even stagnating convergence of the residuals in the iterative solver, as will be shown in Section 4.6.

For practical applications, however, low regularisation parameters may be acceptable since the rheology of many real fluids is better captured by regularised models (e.g. [15,18]). On the other hand, for the purpose of comparison with literature data it is important to minimise these errors, also to avoid possible cancellation of regularisation and discretisation errors that can cause a spurious agreement with literature data. Figure 5 shows an example in which CS on an a very coarse grid with M=20 is nearly identical to the CS on a fine grid with M=200. Despite the latter case is numerically more accurate, both CS are very close to literature data (cf. with Table 7) as a result of errors cancellation.

Fig. 5.

Grid convergence of the Stokes coefficient for Bn=197.5 using M=20 and 200. The difference δ is between CS on an a very coarse grid with M=20 and on a fine grid with M=200.

Grid convergence of the Stokes coefficient for Bn=197.5 using M=20 and 200. The difference δ is between CS on an a very coarse grid with M=20 and on a fine grid with M=200.
Fig. 6.

Convergence of the Stokes coefficient with M1/Mi (Mi=200,100,50,20) for Re=0.01 with Bn=59.59 (left) and Bn=544.6 (right). The percentages indicate the regularisation uncertainty for M=200.

Convergence of the Stokes coefficient with M1/Mi (Mi=200,100,50,20) for Re=0.01 with Bn=59.59 (left) and Bn=544.6 (right). The percentages indicate the regularisation uncertainty for M=200.

The choice of the proper regularisation parameter is found to be problem-dependent. With regard to the laminar flow around a sphere, previous numerical studies used very different values, both dimensional and non-dimensional. For instance, Blackery and Mitsoulis [8] used m=200 s, whereas Beaulne and Mitsoulis [6] suggest to keep the product MBn equal to 103. On the other hand, Gavrilov et al. [22] and Nirmalkar et al. [34,35] have used M=103 and m=106 s, respectively.

For the current work, the regularisation parameter for each test case has been gradually increased until the regularisation uncertainty became less than 1%. The question, however, is how to estimate the regularisation uncertainty. Frigaard and Nouar [20] showed that, for typical Bingham shear flows, the Papanastasiou regularisation errors tend to zero with first order as 1/M0. As a first approximation, it is thus reasonable to estimate the regularisation uncertainty with the same method used for the discretisation uncertainty, with the grid spacing replaced by 1/M (Fig. 6). The final selected values of M that ensured a regularisation uncertainty below 1% are reported in Table 3. We found that the maximum and average regularisation uncertainties in percentage of the drag coefficient are 0.77% and 0.2%, respectively.

Table 3

Selected non-dimensional regularisation parameter M=mU/d for each test case

Re=10Re=100Creep flow, Re<1


Finally, we remark that the values reported in Table 3 may not be sufficient to accurately capture the so-called yielded surface, i.e. the locus of points where τ=τ0 (see e.g. [10,28]). This aspect, however, is both beyond the scope of this work and unimportant for practical maritime applications, hence it is no further discussed.

4.6.Iterative convergence and viscosity interpolation scheme

As mentioned in the previous section, iterative convergence can become difficult when using large regularisation parameters. This issue seems rather common when using SIMPLE-like algorithms (see e.g. [40,41]). In this work, we found that the rate of convergence of the residuals is significantly influenced by the choice of the interpolation scheme for the viscosity.

Within the finite-volume method, the diffusive term in Eq. (1) requires the evaluation of the apparent viscosity, μ, at the cell faces. Since in ReFRESCO the computational node coincides with the cells’ centroid, the face value must be obtained by interpolation. Assuming for simplicity that a cell face e is halfway between two Cartesian grid cells P and E, the simplest second-order interpolation scheme to calculate μe reads

Another simple alternative is to first evaluate γ˙e from linear interpolation of γ˙P and γ˙E and then to calculate the viscosity at the face e, i.e.

The two schemes have the same computational cost, thus, in principle, there is no reason to prefer one scheme to the other. Furthermore, when we performed a code verification exercise (not shown here) similar to that in [30], the two schemes exhibited the same accuracy and rate of convergence of the residuals. However, on the present problem, the first scheme (Eq. (12)) turned out to be more robust as it was possible to obtain a fully converged solution using relatively large regularisation parameters,11 as shown in Fig. 7. With the second scheme (Eq. (13)), on the other hand, the residuals were stagnating already with regularisation parameters that were one order of magnitude lower than those used to obtain a fully-converged solution with the first scheme. Therefore, an interpolation scheme that uses the cell centre values of the apparent viscosity was eventually adopted in ReFRESCO.

Fig. 7.

Effect of the two viscosity interpolation schemes on the iterative convergence of the pressure residuals. Test case: Re=0.01, Bn=59.59.

Effect of the two viscosity interpolation schemes on the iterative convergence of the pressure residuals. Test case: Re=0.01, Bn=59.59.

For this work and with the setup described above, the iterative convergence criterion was set to L<108 and the iterative uncertainties were estimated with the method of Eça and Hoekstra [16]. For all the test cases, we found that the iterative uncertainties were virtually zero (<0.0005%).

4.7.Total numerical uncertainty

The final numerical uncertainty in the computed drag coefficient was calculated assuming that all the uncertainty components are dependent on each other, hence

where Udom, Udis, Ureg and Uit are the uncertainties produced by the finite domain size, the domain discretisation, the regularisation parameter and the iterative method, respectively. The total numerical uncertainties are reported in Table 4 and, for all the test cases, they do not exceed 1%.

Table 4

Total numerical uncertainties, Unum, in percentage of the computed drag coefficient



5.Results and discussion

5.1.Flow field

Some features of the computed flow field are now discussed. More detailed descriptions have been already discussed by other authors [7,22,34,35], thus they will not be completely repeated here. The contours of the velocity, shear rate and viscosity for Re=10 and 100 are shown in Figs 8 to 10, respectively. Note that flow field for Bingham creep flow (Re=0.01) is not shown as it is very similar to the the case with Re=10 and Bn=100.

Fig. 8.

Velocity contour and streamlines. The flow is from bottom to top.

Velocity contour and streamlines. The flow is from bottom to top.
Fig. 9.

Contour plot of the shear rate γ˙. The black isoline corresponds to γ˙=1 s1 for Bn=0, whereas it corresponds to τ=τ0 for Bn=10 and 100.

Contour plot of the shear rate γ˙. The black isoline corresponds to γ˙=1 s−1 for Bn=0, whereas it corresponds to τ=τ0 for Bn=10 and 100.
Fig. 10.

Contour plot of μ(γ˙)/k, which is equal to 1 for Newtonian fluids.

Contour plot of μ(γ˙)/k, which is equal to 1 for Newtonian fluids.

For Newtonian and power-law fluids (Bn=0), the flow is attached to the sphere for Re=10, whereas at Re=100 the flow appears separated, exhibiting the characteristic toroidal eddy behind the sphere (Fig. 8, top). The effect of decreasing n is less evident. When n<1, the apparent viscosity becomes lower than the Newtonian viscosity in the region where γ˙>1 (i.e. the region within the black isoline in Fig. 9, top). This leads to a thinner boundary layer for n<1 compared to Newtonian fluids. Another effect of reducing n is the smaller wake eddy. The latter observation agrees qualitatively with the results of Dhole et al. [14] for power-law fluids.

For yield stress fluids (Bn>0), the viscosity increases significantly, especially in the undisturbed flow region away from the sphere (Fig. 10, middle and bottom panels). The large viscosity damps advection, leading to the disappearance of the toroidal eddy behind the sphere. The fore-aft symmetry typical of creep flow is thus restored. For Bn=100, the flow appears symmetrical with respect to the equatorial plane (z=0), both for Re=10 and 100 (Fig. 8, bottom).

When Bn>0, the fluid far upstream is undeformed with very high viscosity, i.e. it behaves as a solid-like material (τ<τ0). The fluid is then deformed as it encounters the sphere and the viscosity is reduced (τ>τ0). The black isoline in the middle and bottom plots of Fig. 9 identifies the yielded surface, i.e. the locus of points where τ=τ0. This surface shrinks with increasing Bn and decreasing n, whereas it grows with higher Re. These observations are in line with those of Nirmalkar et al. [35]. Furthermore, small unyielded regions are observed near the stagnation points. These regions are usually referred to as ‘polar caps’ [7,8].

As Bn increases, the solid-like region (uniform undisturbed flow) tends to expand towards the sphere, thus ‘squeezing’ the fluid close to the sphere. As a result, the fluid has to increase its velocity in order to keep the flow rate constant along the tube cross-sections. Another way to see this is that the streamlines become denser near the sphere.

Furthermore, for yield stress fluids, the velocity in the equatorial plane exhibits a local maximum, which leads to a local viscosity maximum22 (Fig. 10), as also observed in [22,34]. Such steep variation of viscosity over a relative short distance was a major cause for the difficult iterative convergence. We found in fact that the largest residuals were always near the local viscosity maximum.

For ideal (i.e. non-regularised) Bingham and Herschel–Bulkley fluids, the fluid region far from the sphere would have an infinite viscosity and zero shear rate. The effect of the regularisation is evident in Figs 9 and 10: the shear rate is very low (but not zero) and the viscosity is very large (but not infinite).

Fig. 11.

Axial velocity profiles for Bingham fluids (n=1) at Re=100; (a) along y on the plane z=0; (b) along the centreline (x=y=0) near the rear stagnation point.

Axial velocity profiles for Bingham fluids (n=1) at Re=100; (a) along y on the plane z=0; (b) along the centreline (x=y=0) near the rear stagnation point.

Finally, the velocity profiles for Bingham fluids at Re=100 are plotted in Fig. 11 and compared with available literature data [22,34]. The above mentioned disappearance of the recirculation region with increasing Bn is clearly visible in Fig. 11(b). Overall, our results agree fairly well with literature data, especially with those of Gavrilov et al. Some visible discrepancies are observed in Fig. 11(a) for Bn=100, which probably stem from the different settings (e.g. grids, regularisation parameters, domain size, etc.) used in this work and in the work of Nirmalkar et al. [34].

In conclusion, the flow field is qualitatively in line with previous studies from the literature, and it quantitatively agrees with the results of Gavrilov et al. and Nirmalkar et al.

5.2.Drag coefficients and comparison with literature

The drag coefficient is often the only quantity of interest in many engineering applications. It is therefore a meaningful quantity to assess whether the present code can reproduce the results from the literature.

Table 5

Pressure, frictional and total drag coefficients for all the considered test cases



The drag coefficients and its components are reported in Table 5. The ratio CDp/CDf is observed to increase with the non-Newtonian character of the fluid (i.e. with higher Bn or lower n). This observation is consistent with previous numerical studies [14,22,34]. Furthermore, CD appears to decrease with n, which is due to the shear-thinning effect the reduces the viscosity in the high-shear region near the sphere.

Table 6

Drag coefficient, CD, from the present calculations and the literature. All the literature data is from numerical simulations and not from the correlations proposed in the respective articles

Data fromRe=10Re=100

Bn=0present work4.3084.1784.0211.0890.9400.790
Dhole [14]4.2814.0863.7691.0620.9210.759
Tripathi [43]4.314.173.761.020.9200.780
Gavrilov [22]4.3344.1673.9691.0950.9410.787
Bn=10present work43.3940.2237.314.7044.3343.983
Gavrilov [22]42.1538.9035.914.5874.2173.880
Nirmalkar [34,35]43.6340.4037.434.7434.3634.015
Bn=100present work320.8309.2300.632.1831.0230.16
Gavrilov [22]307.9295.4284.730.9529.6928.61

A direct comparison of CD with the literature data is given in Tables 6 and 7. Note that making a rigorous validation is not possible because the numerical uncertainties for the literature data are not known. Nevertheless, we have applied the validation procedure proposed by ASME [3] by replacing experimental data with numerical data from the literature. According the procedure, the modelling error, δmodel, is estimated by comparing two quantities: the (expanded) validation uncertainty,

and the comparison error,
where S is our numerical solution value and D is the literature data. Uinput is the uncertainty due to the input parameters, which is zero for the present work. Unum is the uncertainty in our numerical data, and it is reported in Table 4. Ulit is the uncertainty in the literature data. Since this is unknown, it was assumed Ulit=1%.

Table 7

Stokes coefficients, CS, for Bingham creep flow from the present calculations and the literature

Data fromRe=0.01

present work6.39515.2782.78253.6428.0674.5
Beris [7]6.3915.2482.77252.2426.9669.7
Liu [28]6.3815.2182.67253.6426.0671.9
Nirmalkar [34]15.2582.83427.5673.5

E and Uval define an interval within which δmodel falls, i.e.

It is however remarked that, since the comparison is made with literature instead of experiments, the modelling errors are expected to be zero. A successful validation with literature data must thus yield |E|Uval.

Fig. 12.

Comparison error (± validation uncertainty) for the drag coefficient in percentage of our data for Re=10 and 100 (top) and for Re=0.01 (bottom).

Comparison error (± validation uncertainty) for the drag coefficient in percentage of our data for Re=10 and 100 (top) and for Re=0.01 (bottom).

The comparison errors and the validation uncertainties are plotted in Fig. 12, from which the following observations are made:

  • For power-law fluids (Bn=0), the maximum |E| with Dhole et al. [14] and Tripathi et al. [43] is about 7%, which is well outside the uncertainty range. We found that the difference with respect to Dhole et al. lies in the pressure component, where |E| reaches about 20% relative to our data. The reasons for such difference are however not known. It is possible that the uncertainty in the literature data is actually greater than 1%, therefore assuming Ulit=1% might have led to an underestimation of the validation uncertainty. Nevertheless, the agreement with the more recent results of Gavrilov et al. [22] is better, with comparison errors are either within or very close to Uval.

  • For Bingham and Herschel–Bulkley fluids (Bn>0) at Re=10 and 100, the larger discrepancies are found with respect to the data of Gavrilov et al., with |E|>Uval for all cases. On the other hand, the agreement with Nirmalkar et al. [34,35] is excellent, with differences that are well within the validation uncertainties.

  • The best agreement with literature data is achieved for Bingham creep flow (Re=0.01), with all the comparison errors being within the validation uncertainties.

  • Our results tend to be on the over-predicting side, i.e. E leans towards the right side of Fig. 12. We observed that CD tends to decrease with grid refinement and to increase with higher regularisation parameters. It is thus possible that the other authors have used finer grids and/or lower regularisation parameters. Another possible reason for a systematic increase of the sphere drag may be expected from the effect of the tube-to-sphere diameter ratio [11]. However, it was found in Section 4.3 that the uncertainty in CD due to the finite domain size was between 0.05% and 0.12% for power-law fluids and below 0.001% for the other non-Newtonian cases. Thus, the reasons for the over-predicting trend remain unclear.


We have tested the viscous-flow solver ReFRESCO on the laminar non-Newtonian flow around a sphere as a stepping stone towards more practical maritime applications.

Some difficulties were encountered when using large regularisation parameters, which led to stagnating residuals and thus large iterative errors. A determining factor turned out to be the choice of the viscosity interpolation scheme. Although obtaining a fully converged solutions remained challenging, the improved iterative convergence provided a significant reduction of the iterative and regularisation errors. This allowed a more compelling comparison with literature.

The flow field around the sphere exhibited the expected behaviour for the considered test cases and it is both qualitatively and quantitatively consistent with previous numerical studies. The maximum difference in drag coefficient between the calculated values and the values from the literature reached 7% (i.e. outside the validation uncertainty range) for power-law fluids, whereas it was insignificant for Bingham creep flow. Overall, the agreement with previous studies is good, with discrepancies that are, in most cases, close or within the validation uncertainties.

Combining the evidence of our previous code verification exercises [30,31] and the results of this work, it is concluded that the power-law, Bingham and Herschel–Bulkley models are implemented correctly and that the code is capable of reproducing literature data with good accuracy despite the high non-linearity introduced by the non-Newtonian viscosity. This provides confidence to employ ReFRESCO for more complex applications such as a ship sailing through fluid mud. Finally, present data can be used for future verifications and to derive new correlations for the drag coefficient.

Credit author statement

Stefano Lovato: Conceptualization, Software, Investigation, Validation, Visualization, Writing – Original Draft. Serge Toxopeus: Software, Validation, Supervision, Writing - Review & Editing. Just Settels: Validation, Supervision, Writing – Review & Editing. Geert Keetels: Supervision, Funding acquisition, Writing – Review & Editing.

Conflict of interest

Serge Toxopeus is an Editorial Board Member of this journal, but was not involved in the peer-review process nor had access to any information regarding its peer-review.


1 Sufficiently large, at least, to keep the regularisation uncertainty below 1%.

2 Note that, in three dimensions, this local maximum is actually a circle around the sphere.


This research is funded by the Division for Earth and Life Sciences (ALW) with financial aid from the Netherlands Organization for Scientific Research (NWO) with Grant No.ALWTW.2016.029. Calculations were performed on the Reynolds (TU Delft) and Marclus4 (MARIN) clusters. The research support from MARIN is partly funded by the Dutch Ministry of Economic Affairs. The authors wish to thank Andrey Gavrilov and Neelkanth Nirmalkar for kindly providing their data. This study was carried out within the framework of the MUDNET academic network:



R.W. Ansley and T.N. Smith, Motion of spherical particles in a Bingham plastic, AIChE Journal 13: (6) ((1967) ), 1193–1196. doi:10.1002/aic.690130629.


A.S. Arabi and R.S. Sanders, Particle terminal settling velocities in non-Newtonian viscoplastic fluids, Canadian Journal of Chemical Engineering 94: (6) ((2016) ), 1092–1101. doi:10.1002/cjce.22496.


ASME PTC Committee, Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer: ASME V&V 20, The American Society of Mechanical Engineers (ASME) (2009).


D.D. Atapattu, R.P. Chhabra and P.H.T. Uhlherr, Wall effect for spheres falling at small Reynolds number in a viscoplastic medium, Journal of Non-Newtonian Fluid Mechanics 38: (1) ((1990) ), 31–42. doi:10.1016/0377-0257(90)85031-S.


D.D. Atapattu, R.P. Chhabra and P.H.T. Uhlherr, Creeping sphere motion in Herschel–Bulkley fluids: flow field and drag, Journal of Non-Newtonian Fluid Mechanics 59: (2–3) ((1995) ), 245–265. doi:10.1016/0377-0257(95)01373-4.


M. Beaulne and E. Mitsoulis, Creeping motion of a sphere in tubes filled with Herschel–Bulkley fluids, Journal of Non-Newtonian Fluid Mechanics 72: (1) ((1997) ), 55–71. doi:10.1016/S0377-0257(97)00024-4.


A.N. Beris, J.A. Tsamopoulos, R.C. Armstrong and R.A. Brown, Creeping motion of a sphere through a Bingham plastic, Journal of Fluid Mechanics 158: ((1985) ), 219–244. doi:10.1017/S0022112085002622.


J. Blackery and E. Mitsoulis, Creeping motion of a sphere in tubes filled with a Bingham plastic material, Journal of Non-Newtonian Fluid Mechanics 70: (1–2) ((1997) ), 59–77. doi:10.1016/S0377-0257(96)01536-4.


B.J. Briscoe, M. Glaese, P.F. Luckham and S. Ren, The falling of spheres through Bingham fluids, Colloids and Surfaces 65: (1) ((1992) ), 69–75. doi:10.1016/0166-6622(92)80176-3.


G.R. Burgos, A.N. Alexandrou and V. Entov, On the determination of yield surfaces in Herschel–Bulkley fluids, Journal of Rheology 43: (3) ((1999) ), 463–483. doi:10.1122/1.550992.


R.P. Chhabra, Bubbles, Drops, and Particles in Non-Newtonian Fluids, CRC Press, (2006) . doi:10.1201/9781420015386.


P. Coussot and J.M. Piau, On the behavior of fine mud suspensions, Rheologica Acta 33: (3) ((1994) ), 175–184. doi:10.1007/BF00437302.


G. Dazhi and R.I. Tanner, The drag on a sphere in a power-law fluid, Journal of Non-Newtonian Fluid Mechanics 17: (1) ((1985) ), 1–12. doi:10.1016/0377-0257(85)80001-X.


S.D. Dhole, R.P. Chhabra and V. Eswaran, Flow of power-law fluids past a sphere at intermediate Reynolds numbers, Industrial & Engineering Chemistry Research 45: (13) ((2006) ), 4773–4781. doi:10.1021/ie0512744.


N.Q. Dzuy and D.V. Boger, Direct yield stress measurement with the Vane method, Journal of Rheology 29: (3) ((1985) ), 335–347. doi:10.1122/1.549794.


L. Eça and M. Hoekstra, Evaluation of numerical error estimation based on grid refinement studies with the method of the manufactured solutions, Computers & Fluids 38: (8) ((2009) ), 1580–1591. doi:10.1016/j.compfluid.2009.01.003.


L. Eça and M. Hoekstra, A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies, Journal of Computational Physics 262: ((2014) ), 104–130. doi:10.1016/


K.R.J. Ellwood, G.C. Georgiou, T.C. Papanastasiou and J.O. Wilkes, Laminar jets of Bingham-plastic liquids, Journal of Rheology 34: (6) ((1990) ), 787–812. doi:10.1122/1.550144.


J.H. Ferziger and M. Perić, Computational Methods for Fluid Dynamics, Springer, Berlin Heidelberg, Berlin, Heidelberg, (1996) . doi:10.1007/978-3-642-97651-3.


I.A.A. Frigaard and C. Nouar, On the usage of viscosity regularisation methods for visco-plastic fluid flow computation, Journal of Non-Newtonian Fluid Mechanics 127: (1) ((2005) ), 1–26. doi:10.1016/j.jnnfm.2005.01.003.


Z. Gao, H. Yang and M. Xie, Computation of flow around Wigley hull in shallow water with muddy seabed, Journal of Coastal Research 73: ((2015) ), 490–495. doi:10.2112/SI73-086.1.


A.A. Gavrilov, K.A. Finnikov and E.V. Podryabinkin, Modeling of steady Herschel–Bulkley fluid flow over a sphere, Journal of Engineering Thermophysics 26: (2) ((2017) ), 197–215. doi:10.1134/S1810232817020060.


D.I. Graham and T.E.R. Jones, Settling and transport of spherical particles in power-law fluids at finite Reynolds number, Journal of Non-Newtonian Fluid Mechanics 54: (C) ((1994) ), 465–488. doi:10.1016/0377-0257(94)80037-5.


F. Irgens, Rheology and Non-Newtonian Fluids, Springer International Publishing, Cham, (2014) , pp. 1–190. doi:10.1007/978-3-319-01053-3.


S. Kaidi, E. Lefrançois and H. Smaoui, Numerical modelling of the muddy layer effect on ship’s resistance and squat, Ocean Engineering 199: ((2020) ), 106939. doi:10.1016/j.oceaneng.2020.106939.


C.M. Klaij and C. Vuik, SIMPLE-type preconditioners for cell-centered, colocated finite volume discretization of incompressible Reynolds-averaged Navier–Stokes equations, International Journal for Numerical Methods in Fluids 71: (7) ((2013) ), 830–849. doi:10.1002/fld.3686.


K. Koziol and P. Glowacki, Determination of the free settling parameters of spherical particles in power law fluids, Chemical Engineering and Processing: Process Intensification 24: (4) ((1988) ), 183–188. doi:10.1016/0255-2701(88)85001-3.


B.T. Liu, S.J. Muller and M.M. Denn, Convergence of a regularization method for creeping flow of a Bingham material about a rigid sphere, Journal of Non-Newtonian Fluid Mechanics 102: (2) ((2002) ), 179–191. doi:10.1016/S0377-0257(01)00177-X.


S. Lovato, G.H. Keetels, S.L. Toxopeus and J.W. Settels, An Eddy-viscosity model for turbulent flows of Herschel–Bulkley fluids, Journal of Non-Newtonian Fluid Mechanics 301: ((2022) ), 104729. doi:10.1016/j.jnnfm.2021.104729.


S. Lovato, S.L. Toxopeus, J.W. Settels, G.H. Keetels and G. Vaz, Code verification of non-Newtonian fluid solvers for single- and two-phase laminar flows, Journal of Verification, Validation and Uncertainty Quantification 6: (2) ((2021) ). doi:10.1115/1.4050131.


S. Lovato, G. Vaz, S.L. Toxopeus, G.H. Keetels and J.W. Settels, Code verification exercise for 2D Poiseuille flow with non-Newtonian fluid, in: Numerical Towing Tank Symposium (NuTTS), (2018) .


T.F. Miller and F.W. Schmidt, Use of a pressure-weighted interpolation method for the solution of the incompressible Navier–Stokes equations on a nonstaggered grid system, Numerical Heat Transfer 14: (2) ((1988) ), 213–233. doi:10.1080/10407788808913641.


K.A. Missirlis, D. Assimacopoulos, E. Mitsoulis and R.P. Chhabra, Wall effects for motion of spheres in power-law fluids, Journal of Non-Newtonian Fluid Mechanics 96: (3) ((2001) ), 459–471. doi:10.1016/S0377-0257(00)00189-0.


N. Nirmalkar, R.P. Chhabra and R.J. Poole, Numerical predictions of momentum and heat transfer characteristics from a heated sphere in yield-stress fluids, Industrial & Engineering Chemistry Research 52: (20) ((2013) ), 6848–6861. doi:10.1021/ie400703t.


N. Nirmalkar, R.P. Chhabra and R.J. Poole, Effect of shear-thinning behavior on heat transfer from a heated sphere in yield-stress fluids, Industrial & Engineering Chemistry Research 52: (37) ((2013) ), 13490–13504. doi:10.1021/ie402109k.


T.C. Papanastasiou, Flows of materials with yield, Journal of Rheology 31: (5) ((1987) ), 385–404. doi:10.1122/1.549926.


P. Saramito and A. Wachs, Progress in numerical simulation of yield stress fluid flows, Rheologica Acta 56: (3) ((2017) ), 211–230. doi:10.1007/s00397-016-0985-9.


P.R. Souza Mendes and E.S.S. Dutra, Viscosity function for yield-stress liquids, Applied Rheology 14: (6) ((2004) ), 296–302. doi:10.1515/arh-2004-0016.


G.G. Stokes, On the effect of the internal friction of fluids on the motion of pendulums, in: Mathematical and Physical Papers, Cambridge University Press, Cambridge, (2010) , pp. 1–10. doi:10.1017/CBO9780511702266.002.


A. Syrakos, G.C. Georgiou and A.N. Alexandrou, Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method, Journal of Non-Newtonian Fluid Mechanics 195: ((2013) ), 19–31. doi:10.1016/j.jnnfm.2012.12.008.


A. Syrakos, G.C. Georgiou and A.N. Alexandrou, Performance of the finite volume method in solving regularised Bingham flows: Inertia effects in the lid-driven cavity flow, Journal of Non-Newtonian Fluid Mechanics 208–209: ((2014) ), 88–107. doi:10.1016/j.jnnfm.2014.03.004.


H. Tabuteau, P. Coussot and J.R. de Bruyn, Drag force on a sphere in steady motion through a yield-stress fluid, Journal of Rheology 51: (1) ((2007) ), 125–137. doi:10.1122/1.2401614.


A. Tripathi, R.P. Chhabra and T. Sundararajan, Power law fluid flow over spheroidal particles, Industrial & Engineering Chemistry Research 33: (2) ((1994) ), 403–410. doi:10.1021/ie00026a035.


L. Valentik and R.L. Whitmore, The terminal velocity of spheres in Bingham plastics, British Journal of Applied Physics 16: (8) ((1965) ), 1197–1203. doi:10.1088/0508-3443/16/8/320.


G. Vaz, F. Jaouen and M. Hoekstra, Free-surface viscous flow computations: Validation of URANS code FRESCO, in: Proceedings of OMAE2009, Honolulu, Hawaii, USA, (2009) , pp. 425–437. doi:10.1115/OMAE2009-79398.


R.W. Wurpts, 15 years experience with fluid mud: Definition of the nautical bottom with rheological parameters, Terra et Aqua (2005).