Application of a Projection Method for Simulating Flow of a Shear-Thinning Fluid

: In this paper, a ﬁrst-order projection method is used to solve the Navier–Stokes equations numerically for a time-dependent incompressible ﬂuid inside a three-dimensional (3-D) lid-driven cavity. The ﬂow structure in a cavity of aspect ratio δ = 1 and Reynolds numbers ( 100, 400, 1000 ) is compared with existing results to validate the code. We then apply the developed code to ﬂow of a generalised Newtonian ﬂuid with the well-known Ostwald–de Waele power-law model. Results show that, by decreasing n (further deviation from Newtonian behaviour) from 1 to 0.9, the peak values of the velocity decrease while the centre of the main vortex moves towards the upper right corner of the cavity. However, for n = 0.5, the behaviour is reversed and the main vortex shifts back towards the centre of the cavity. We moreover demonstrate that, for the deeper cavities, δ = 2, 4, as the shear-thinning parameter n decreased the top-main vortex expands towards the bottom surface, and correspondingly the secondary ﬂow becomes less pronounced in the plane perpendicular to the cavity lid.

Steady two-dimensional (2-D) LDC flow has been studied in literature at different Reynolds number [12,16]-Re = U lid L/ν, where U lid is the lid velocity, L the cavity length and ν the kinematic viscosity of the fluid. The 2-D LDC, however, deviates from actual experiments in 3-D mainly due to the no-slip conditions imposed by end walls in the cross-flow direction span [9][10][11]17]. Albensoeder et al. [17] reported that the fluid bifurcates from steady to periodic, and eventually to turbulent flow when it goes above the critical Reynolds number, Re cr ∼ O 10 3 . This is, however, dependent on cavity aspect ratios, i.e., length-to-width, λ = L/W, and depth-to-width, δ = H/L [18,19]. Most fluids in nature exhibit Newtonian behaviour, in which the dynamic viscosity, µ, at constant temperature and pressure is equal to τ/γ-τ is the shear stress, andγ is the shear rate [20,21]. For an incompressible fluid, the equation governing the tensorial shear stress (in Cartesian coordinates) is [22,23] where µ is the dynamic viscosity, τ is the shear stress tensor, andγ the shear rate tensor. Here,γ is the rate of strain tensor,γ = ∇u + ∇u , which is given by [24][25][26] |γ| = 1 2 I Iγ = 1 2 γ :γ where I Iγ is the second invariant ofγ, ∇u is the velocity gradient tensor and ∇u is its transpose. Similarly, for the stress tensor, τ, we have [27][28][29] |τ| = 1 2 where I I τ is the second invariant of τ. The correlation between shear stress and shear rate defines the flow behaviour of a liquid. For a Newtonian fluid, there is a linear correlation between the shear stress and shear rate, and the slope of this line is the dynamic viscosity, µ. Recognition that stress in a fluid can have a nonlinear or temporal (or both) dependence on the rate of deformation developed in the late nineteenth and early twentieth centuries; we now refer to such materials as non-Newtonian fluids. For a non-Newtonian fluid, a.k.a nonlinear or complex fluid, the slope of shear stress versus shear rate curve is not constant, and hence is a function of τ orγ. The most common category of the complex fluids is the pure viscous fluid in whichγ is determined only by the current value of τ at any particular point in the flow domain. Shear-thinning is perhaps the most widely encountered type of time-dependent non-Newtonian fluid behaviour in engineering practice for ceramics [30][31][32][33]. It is characterised by an apparent viscosity µ a which gradually decreases with increasing shear rate: or, in terms of the apparent viscosity, Equation (5) is known as the Ostwald-de Waele power-law model. Varying n between 0 < n < 1 will yield dµ/dγ < 0; i.e., shear-thinning fluids are characterised by a value of n (power-law index) smaller than unity. Many polymer melts and solutions exhibit a value of n in the range 0.3-0.7 depending upon the concentration and molecular weight, etc., of the polymer. Even smaller values of power-law index (n∼0.1-0.15) are encountered with fine particle suspensions such as, e.g., kaolin-in-water and bentonite-in-water. Naturally, the smaller the value of n, the more shear-thinning the material is. The other constant, k, in Equation (5) (consistency index) is a measure of consistency of the substance, which is well described by Jabbari and Hattel [31].
A computational fluid dynamics (CFD) approach to problems showing non-Newtonian behaviour is generally more complicated than that for a corresponding Newtonian case due to the complexity introduced by the constitutive equation of the employed model in the diffusion terms of the Navier-Stokes equations. In particular, it is clear from Equation (5) that the proposed model for viscosity introduces nonlinearities in the diffusion terms-the highest-order differential operators. Different numerical methods have been employed in the past in order to simulate viscous and viscoplastic fluid flows; the most widely known simulations have employed power-law and Bingham models. However, for the LDC problem, most available studies have focused on modelling elastic and viscoelastic fluid flows: see, for example, [34][35][36]. There are also some studies related to the LDC problem for viscoplastic fluids (power-law and Bingham models) in which rather different approaches were used for solving the NS equations, viz., the smoothed particle hydrodynamics (SPH) method [37], and the PIM (point interpolation method) meshfree method [38]. Despite remarks by these authors, there seems to be no clear-cut overall advantage to non-standard methods in the present context.
In this work, we present results from simulation of flow inside a 3-D LDC. We begin with the mathematical formulation and brief discussion of the numerical analytic procedures employed, all of which are fairly well known. We validate a particular Fortran code in a range of Re numbers corresponding to steady flow by using existent, highly-accurate, direct numerical simulations for comparisons with our results. We then proceed to a problem involving a generalised Newtonian fluid with the shear-thinning power-law behaviour. The influence on the flow pattern in the 3-D cavity with several depth-to-width aspect ratios, δ, as well as Reynolds number, are investigated. These values of Reynolds number and aspect ratios are relevant to processing and manufacturing of ceramics and polymers.

Mathematical Model
Projection methods such as proposed by Chorin [39], Bell et al. [40] and Gresho [41], e.g., are used in problems involving time-dependent incompressible fluids; the diffusion-convection equations are first solved to predict intermediate velocities which are then projected onto the space of divergence-free vector fields [41]. The incompressible NS equations (the differential form of momentum equation combined with the continuity equation) may be written as where ρ is density, and ν is kinematic viscosity. These equations are solved in the interior of a domain Ω consisting of the unit cube, (0, 1) × (0, 1) × (0, 1), shown in Figure 1, with initial conditions u ≡ 0 and no-slip boundary conditions. In the projection 1 method analyzed in [41], and employed here, an intermediate velocity, u * , is computed using the momentum equations with pressure gradient terms absent-given here in semi-discrete form: where u n is the velocity at the n th time level (u = (u, v, w) T ). Equation (7) comprises an easily-solved multi-dimensional system of Burgers' equations. In the present work, the system is coded using a Fortran 77/90 research code based on Douglas and Gunn time splitting [42] with quasilinearisation (a function-space Newton's method) of nonlinear advective terms.
In the next (projection) step, the velocity at the (n + 1)th time level is computed as follows: The right-hand side of the above equation requires knowledge of the pressure, p, at (n + 1) level. This is obtained by taking the divergence and requiring that ∇ · u n+1 = 0, which is the divergence-free (continuity) condition, thereby leading to the following Poisson equation for p n+1 : It is important to emphasize, however, that this is not the true, physical pressure although it differs from physical pressure by ∼ O(∆t/Re).
Since the boundary condition for u * in Equation (7) is no slip, this leads to a boundary condition for the approximate pressure in the form (see, e.g., Foias et al. [43]) It is important to note that boundary conditions must be imposed on physical boundaries, and, in the context of a staggered grid (see Figure 2), this implies that Dirichlet conditions on velocity components can be directly applied (on the boundary) only for the normal-direction component on each boundary. The tangential components must be averaged across the boundary by employing an "image" cell outside the solution domain. Neumann conditions used here for pressure are also implemented by means of image cells, as must also be done in the context of unstaggered gridding. This is done in all cases to maintain second-order accuracy of interior-point discretizations. While using transient simulation, it is important to make sure that the Courant-Friedrichs-Lewy (CFL) condition is preserved. For a given resolution size, the CFL condition enforces an upper limit on the local time step [45]. In general, the CFL condition for a three-dimensional, transient problem is as follows: where C max is normally 1 for explicit solvers. For implicit or semi-implicit, this value can increase up to a bigger value (e.g., 30), and, hence, choosing larger time step size will not harm the accuracy of the solution. However, in the presence of nonlinear phenomena, i.e., very strong shockwaves and non-Newtonian viscosity, an implicit scheme linearising the governing equations will exhibit a restraint on C max . This is due to a variation in time of the linearisation Jacobian from one iteration to the next. In the light of above a relatively small time step size, ∆t = 0.003 s is used in this study. Both Equations (6a) and (7) are valid in the case of constant viscosity. However, inserting a non-Newtonian viscosity will increase the nonlinearity of the system of equations. In order to account for the non-Newtonian behaviour of the fluid and reduce the nonlinearity of the system, one possibility is to replace the generalised apparent kinematic viscosity, Equation (6a), by the constant kinematic viscosity derived from Equation (5) [46]. The kinematic viscosity, however, is a function of ∆u, which makes the viscous term nonlinear. According to Taylor's series of expansion, the viscous terms become: in which D (∆u) = 1 2 ∆u + ∆u , and Inserting Equation (13) into Equation (12) and neglecting the higher order term O ∆u − ∆u old 2 will result in the Newton type iteration. Neglecting the second term on the right-hand side of Equation (13), the viscous term becomes linearised and fixed-point iterations will be achieved as follows: The rheological parameters (k and n) of the fluid are those previously studied by Jabbari and Hattel [31], as a common shear-thinning non-Newtonian fluid used in producing thin substrates for fuel cells and magnetic refrigeration applications. Results from the rheology experiment conducted by Jabbari et al. [30] showed that the La 0.85 Sr 0.15 MnO 3 (LSM) ceramic slurry follows the Ostwald-de Waele power-law behaviour with k = 3.31 and n = 0.9. In the present study, while keeping k constant, the n parameter is chosen to have different values of 0.9 and 0.5 to investigate the influence of greater deviations from Newtonian fluids (shown in Figure 3).

Results and Discussion
We first present a fairly extensive set of results for a Newtonian fluid (applying n = 1 in Equation (5)) in a cubical cavity, for which much is available in the existing literature [7,11,17,47,48]. This provides thorough validation of the new code. Then, we conduct similar studies for a generalised shear-thinning case with different aspect ratios as well as the Re numbers.

Code Validation for Newtonian Fluid
Prior to proceeding with validation in the context of shear-thinning fluids, it is essential to ensure grid independence of the computations for a simpler Newtonian fluid. The aim is to achieve a maximum change of no more than 1% for a certain variable from one grid level to a more refined one. For this purpose, three different uniform grids with the same number of divisions N x = N y = N z in each coordinate direction were used: namely, 128 3 , 256 3 , and 384 3 grids.
Results are shown in Figure 4 for the normalised x-velocity component u = u/U lid along the vertical centreline for a Newtonian fluid. Selecting as a criterion the value of u min , it can be seen that grid independence is achieved with the 256 3 grid (refining further to 384 3 showed a change of less than 1%), and, therefore, based on this grid refinement study, this grid is used for all later computations.    To further validate the developed code, we present a comparison of velocity profiles obtained from the present simulations with the Navier-Stokes solutions published in the literature by Albensoeder and Kuhlmann [47]. These cases correspond to flow in a lid-driven cubical cavity (δ = 1) at different Reynolds numbers. Shown in Figure 5 are the normalised xand y-components of velocity profiles along the vertical (CL y ) and horizontal (CL x ) centreline planes, respectively, in Ω (see Figure 1). Solid lines represent results obtained from the present simulations, whereas the discrete points denote the 3-D pseudo-spectral solution of Navier-Stokes equations in [47]. It is seen that the comparison is excellent at all Reynolds numbers considered here.

Shear-Thinning Fluid
Results of simulation for generalised Newtonian power-law behaviour in the LDC are presented here for three different Reynolds numbers (Re = 100, 400, 1000). For the power-law fluid, the Re is calculated by Re = ρL n /k · U n−2 lid . The shear rate of the flow is updated in every time step through the generalised kinematic viscosity. In order to validate the present code for generalised Newtonian power-law behaviour, we have compared the results for δ = 1 and n = 0.5 at Re = 100 with those published by Neofytou [48]. Figure 6 depicts the normalised x-components of velocity profiles in the vertical (CL y ) centreline plane, showing a fairly good agreement. Apparent discrepancies are due to the fact that the referenced work was a 2-D calculation. The critical values for Reynolds number, Re cr , for the transition from steady to periodic flow (a Hopf bifurcation), and for different values of k and n (shear-thinning parameters) are shown in Figure 7. It is seen that, for the current study, k = 3.31, for all n-values, the critical Reynolds number is above 1000. Results for δ = 1 together with three different n values (see Equation (5)) for the shear-thinning parameter are illustrated in Figure 8. It should be noted that the more n is decreased the greater the deviation from Newtonian behaviour. The results show that for all cases increasing Re causes higher gradients in addition to increasing peak values for u and v; this can be attributed to the gradual reduction of viscous effects that occurs with increasing Re.
The corresponding streamlines coloured by the velocity magnitude inside the cavity for the three different Reynolds numbers as well as three different n values are illustrated in Figure 9. At low Reynolds number, the top row, the region occupied by the primary vortex is relatively small due to the highly diffusive transport of vorticity, and considerable thickness of the boundary layers near the walls. On the other hand, at high Reynolds numbers, the middle and bottom rows, the primary vortex grows in size and occupies most of the central region of the cavity cross-section; the boundary layer is seen to be limited to very thin regions near the walls, as should be expected. In these cases, the vorticity transport is dominated by advection created due to motion of the top lid. The plot of streamlines revealed that a pair of small secondary vortices appear near the bottom corners of the cavity.
As Re is increased, the centre of the primary vortex (x vor , y vor ) shifts towards the downstream wall, and simultaneously it moves down towards the bottom wall of the cavity (as seen in [47]). In 2-D simulations, the centre of the primary vortex has been observed to shift towards the downstream side of the cavity until the Reynolds number increases to about 100, but there onwards it asymptotically moves towards the geometric centre of the cavity [16] with further increase in the Reynolds number. It has been reported [14,16] that the results of 2-D and 3-D simulations agree very well for Re = 0.01 and Re = 100. However, for Re = 400 and Re = 1000, the predicted locations (x vor , y vor ) from 2-D and 3-D flow models are significantly different [14,16]. We believe this difference is due to enhanced motion in the transverse direction. In particular, the entire flow field becomes unsteady at much lower Re in 3-D, and the vortex centre location, itself, does not remain in a fixed location. In particular, although results are still steady for Re = 1000, they transition to periodic behaviour at slightly higher Re, suggesting very long transients at Re = 1000. It is also noted that by increasing the shear-thinning effect, decreasing n, the regions with high velocity magnitudes become narrower (while shifting to the top surface), and hence showing boundary layer effect. Moreover, the secondary vortex tends to diminish with decreasing n at all three Reynolds numbers.  The effect on the corresponding displacement of the primary vortex core (physical point in a 3-D domain) is summarized in Table 1 which presents centre-plane x-y coordinates. It can be seen that the centre of the main vortex shifts slightly towards the lid with decreasing Re in all cases, whereas, for shear-thinning fluids, the displacement is more marked. It is interesting to note that, by reducing the shear-thinning parameter, n, from 1 to 0.9, the main vortex moves towards the lid; however, by further decreasing n to 0.5, the main vortex shifts back to the centre of the cavity. This can be easily interpreted from the velocity profiles shown in Figure 8, where the velocity gradients show less variation along the CL y and CL x for n = 0.5, especially for higher Re. Results of simulations for deeper (δ > 1) cavities are presented with a spanwise aspect ratio of λ = 1. Figure 10 presents the velocity streamlines coloured by the vorticity magnitude for Re = 1000 with δ = 2 (top row) and δ = 4 (bottom row), as well as for different n values (n = 1, 0.9, 0.5). It can be seen that for a Newtonian fluid (top-left) a deep cavity with δ = 2 possesses another main vortex in the bottom region, which has an opposite span-direction compared to the main-top-vortex. The phenomenon occurs for a Newtonian fluid with δ = 4 (bottom-left) leading to a third main vortex in the bottom region which has a different span-direction compared to the vortex above it, but the same as that of the main-top vortex.
For a shear-thinning fluid with a slight deviation from Newtonian (top-mid), n = 0.9, the secondary main vortex (in the bottom region) has expanded somewhat resulting in smaller area being occupied by the main-top vortex. This is much more pronounced for the same fluid with a deeper cavity (bottom-mid), δ = 4, where the secondary main vortex has grown and diffused to the third main-vortex region. It is interesting to see that, for a fluid with even lower n (top-left), the main-top-vortex occupies most of the cavity volume. This phenomenon is due to more shear-thinning effects of the fluid with lower n close to the top surface.
Impact of Re on the main-top vortex for deeper (δ = 2, 4) cavities with n = 0.5 are illustrated in Figure 11. As it can be seen, by increasing the Re (from left to right column), the main-top vortex expands towards the bottom surface of the cavity. It should be highlighted that, for δ > 1, there exists a critical Reynolds number above that the flow near the bottom wall will lose symmetry about the line CL y . The flow symmetry will be regained below the critical Reynolds number or by increasing the aspect ratio [19]. Velocity vector in the x = 0.5 plane (Φ in Figure 1) inside the cavity for δ = 2, together with different Reynolds numbers (Re = 100, 400, 1000) and different values of n = 1, 0.9, 0.5, are given in Figure 12. It is noticed that, by increasing the shear-thinning effect, going from left column to the right, the secondary flow becomes weaker. This is due to the shear-thinning effect, for which the fluid tends to flow more in the in-plane direction (x) rather than in the y-z plane. It can be also seen that, in the plane Φ, a pair of secondary-flow vortices is found at small Reynolds number (Re = 100), and these gradually shift towards the lower corner of the cavity side-walls with increasing Reynolds numbers (going from the top row to the bottom row of the figure). In addition, two pairs of longitudinal vortices begin to appear at the upper corner of the side-walls (in fact, a closer examination revealed the presence of a few additional weak vortices, thus making the flow field quite complex). These secondary vortices gain in strength with increasing Reynolds number. The above mentioned phenomena are somewhat less pronounced for fluid with a high shear-thinning effect. Our examination also showed that the same trend exists for the deeper cavity (δ = 4).

Conclusions
In this paper, the 3-D LDC flow has been studied in rectangular cavities using a first-order projection method. The developed model has been validated against already published literature, and is further extended to simulate flows with generalised Newtonian power-law behaviour as a shear-thinning fluid.
The results show that the vortex structure in the cubed cavity (δ = 1) changes considerably with Re, where increasing Re results in a second primary vortex appearing beside the first primary vortex. For deep cavities (δ > 1), the results show the existence of a series of counter-rotating primary vortices that are placed vertically along the cavity-height. The size and number of these vortices depend on the cavity aspect ratio (δ) and the Reynolds number.
The present study, moreover, shows that the centre of the main vortex shifts towards the lid with decreasing Re in all cases, whereas, for a shear-thinning fluid, the displacement is more marked. It is, moreover, noted that by reducing the shear-thinning parameter, n, from 1 to 0.9, the main vortex shifts towards the lid, however, by further decreasing n to 0.5, the main vortex shifts back to the centre of the box.
Complex flow problems, such as transition to turbulence, and turbulent flow of non-Newtonian power-law fluids in the cavities with different aspect ratios may be of interest for future studies. Migration of a secondary phase, i.e., particles, in the cavity flow of non-Newtonian fluids could be another interesting phenomenon to investigate. This will be beneficial to applications in processing and manufacturing of ceramics and polymers where particles are suspended in non-Newtonian fluid.