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

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.

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.

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.

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.

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/M→0. 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.

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

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.

For this work and with the setup described above, the iterative convergence criterion was set to L∞<10−8 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

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.

Fig. 9.

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.

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.

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.

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.

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,

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

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).

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