Particle Settling in a Non-Newtonian Fluid Medium Processed by Using the CEF Model

Abstract

In this study, the settling of small particles in a non-Newtonian fluid medium is considered. The simulation of this problem according to the fluid mechanics principles may be realized by the flow of a non-Newtonian fluid around a sphere falling along the centerline of a cylindrical tube. The knowledge of the rate of settling of particles in practice is particularly significant in determining the shelf life of materials such as foodstuffs, cleaning materials and many others. Thus, this problem has great importance in many natural and physical processes and in a large number of industrial applications such as chemical, genetic and biomedical engineering operations. The majority of the theoretical, experimental and numerical studies available in the literature deal with Newtonian fluids. Conversely, for non-Newtonian fluids the problem is considerably more complex. It is well-recognized that extensional behaviour in non-Newtonian fluids plays a major role in complex flows. Most non-Newtonian fluids such as polymeric solutions and melts exhibit shear-thinning behaviour. In this study it is aimed to determine the equations governing this process and some important conclusions about the properties of polymeric liquids related to their viscoelastic constitution are drawn. Effectively, it is found that for polymeric liquids, the elastic behaviour characterized by the normal stress coefficients, implies relatively increased normal stresses with respect to the generalized Newtonian fluids, whereas the shear stresses tend to decrease, thus changing somewhat the category of the flow from shear-flow into extensional flow in a small rate. Hence, the viscoelastic property of the polymeric liquids must be stressed by their constitutive equation choice, which led us to the CEF model.

1. Introduction

As an example to the application of Criminale-Ericksen-Filbey (CEF) fluid, we can consider the settling of small particles in a non-Newtonian fluid medium. The geometry of the problem: schematic diagram of a sphere falling through a fluid in a cylinder is shown in (Figure 1). Due to axisymmetry the problem can be considered as 2-D.

Figure 1.

Schematic diagram of a sphere falling through a fluid in a cylinder.

This problem has great importance in a large number of industrial applications. Because these materials are rather polymeric and consequently viscoelastic, it is obligatory to include the elastic behaviour in the analysis. The notations are the usual ones.

The quantities of interest are the stresses: we must consider the effects of the normal stress in addition to those of the shear stresses and to use a constitutive equation which includes normal stress coefficients, ${\upsilon }_1$ and ${\upsilon }_2$ as well as viscosity coefficient $\eta$ , while taking into account their shear-thinning variation [1, 2, 3].

1.1. The non-Newtonian fluid chosen

Many constitutive equations developed from the continuum mechanics or microstructural viewpoints, are used to describe the behaviour of non-Newtonian fluids. Among them the second-order Rivlin-Ericksen model is generally preferred because it describes the real behaviour of the fluids with sufficient accuracy and also because its application is not very cumbersome. On the other hand, the material coefficients used being constant, it is not in good agreement with experimental results in case the shear strain-rate is not very small. The CEF constitutive equation removes this draw-back by taking these coefficients variable and dependent of the shear strain-rate. That is why this model is used in the study of the settling of small particles in a non-Newtonian fluid medium.

1.1.1. The CEF equation

The constitutive equation of the CEF fluid is:

$$\mathit{\tau }=–p\mathit{I}+\eta {\mathit{A}}_{\mathit{1}}+\left({\upsilon }_{\mathit{1}}+{\upsilon }_{\mathit{2}}\right){{\mathit{A}}_{\mathbf{1}}}^{\mathbf{2}}–\frac{1}{2}{\upsilon }_1{\mathit{A}}_2$$ E1

The merit of CEF constitutive equation is that it stresses the dependence of the viscosity coefficient on shear strain-rate, that is, it takes into account the shear thinning (or shear-thickening) effects, and those of normal stresses which are also dependent on the shear strain-rate. In steady-state shear flow an extremely wide class of viscoelastic constitutive equations simplifies to CEF equation. The first term of the CEF equation for $\mathrm{\tau }$ is just $\mathrm{\eta }\left(\dot{\gamma }\right)\dot{\gamma }$ ; the other two terms, containing ${\upsilon }_1$ and ${\upsilon }_2$ , describe the elastic effects associated with the normal stresses [4].

The Rivlin-Ericksen tensors involved in the CEF equation are [5]:

$${\mathit{A}}_{\mathit{1}}=2\mathit{d}=\nabla \mathit{V}+\nabla {\mathit{V}}^{\mathit{T}}$$

$${{\mathit{A}}_{\mathit{1}}}^{\mathit{2}}=4{\mathit{d}}^{\mathit{2}}$$ E2

$${\mathit{A}}_{\mathit{2}}=\nabla \mathit{a}+\nabla {\mathit{a}}^{\mathit{T}}+\mathit{2}\nabla \mathit{V}\nabla {\mathit{V}}^{\mathit{T}}$$

The first invariant of d is null

$$I_d=0$$ E3

The shear strain-rate in term of the second invariant II_d is [6]

$$\dot{\gamma }=2\sqrt{{\mathrm{II}}_{\mathrm{d}}}=\sqrt{\frac{1}{2}\mathrm{tr}{\mathit{A}}_1^2}=\sqrt{\frac{1}{2}\mathrm{tr}{\mathit{A}}_2}$$ E4

The velocity gradients of A₁, ${\mathrm{A}}_1^2$ , A₂ appearing in Eq. (1) are given in extenso.

1.1.1.1. Material functions in steady-state shear flows

At high shear rates, the viscosity of most polymeric liquids decreases with increasing shear rate (Figure 2). For many engineering applications this is the most important property of polymeric fluids.

Figure 2.

Shear-thinning in a typical non-Newtonian fluid.

An especially useful form has been described by Carreau [4] for the viscosity coefficient (n = 0.364).

$$\eta ={\eta }_0{\left(1+32.32\cdot \mathit{tr}\cdot {\mathit{A}}_1^2\right)}^{-0.318}$$ E5

The normal stress coefficients may be handled as below:

If η = η_o and υ₁ = υ₁₀ for $\dot{\gamma }$ = 0, according to $\lambda \dot{\gamma }$ Weissenberg number, η/η_o and $\frac{N_1}{{\dot{\gamma }}^2{\upsilon }_{10}}=\frac{{\upsilon }_1}{{\upsilon }_{10}}$ curves may be displayed as in the Figure 3 [7]. Elastic effects are observable in a steady simple-shear flow through normal stress effects. This is demonstrated in Figure 4 [8].

Figure 3.

Non-linear results.

Figure 4.

The relaxation time (defined as $\frac{N_1}{\tau \dot{\gamma }}=\frac{{\upsilon }_1}{\eta }$ ) is plotted against the shear-rate.

Treating the curve in the Figure 4 by the least square method, the formula below can be found for the first normal stress coefficient υ₁:

$$\frac{{\mathrm{\upsilon }}_1}{\mathrm{\eta }}\approx {10}^{\left[-0.169{\left({\text{log}}_{10}\dot{\gamma }\right)}^2-0.76{\text{log}}_{10}\dot{\gamma }-0.821\right]}$$ E6

The most important points to note about υ₂ are that its magnitude is much smaller than υ₁, usually about 10–20% of υ₁, and that it is negative [4].

Hence, we shall take

$${\upsilon }_2=–0.15{\upsilon }_1$$ E7

Consequently

$${\upsilon }_1+{\upsilon }_2=0.85{\upsilon }_1$$ E8

The choice of Carreau formula is justified by the behaviour similarity of inelastic and viscoelastic fluids concerning the viscosity, and that of the formula about the normal stress coefficients by the fact that the particle settling problem has characteristics close to dilute suspensions.

2. An example to CEF fluid application

As an example to the application of the CEF fluid, we can consider the settling of small particles in a non-Newtonian fluid medium. The simulation of this problem according to the fluid mechanics principles may be realized by the flow of a non-Newtonian fluid around a sphere falling along the centerline of a cylindrical tube [9, 10, 11, 12, 13, 14].

The knowledge of the rate of settling of particles in practice is particularly significant in determining the shelf life of materials such as foodstuffs, cleaning materials and many others. Also, in oil and gas drilling it is important to understand the distribution of loose material, removed by the drill bit and carried to the surface by the drilling mud. Thus, this problem has great importance in many natural and physical processes and in a large number of industrial applications such as chemical, genetic and biomedical engineering operations. The cylindrical tube is considered as stationary. The drag coefficient must be calculated. The equations determining the motion of a sphere in an incompressible fluid under isothermal conditions will be given for a non-Newtonian viscous fluid exhibiting shear-thinning and using cylindrical coordinates [15]. As this study is the simulation of the slow motion of small particles, we can suppose the flow irrotational. Furthermore, due to the data of the considered problem, the motion may be assumed steady, axisymmetric and the fluid incompressible. The global cylindrical coordinate system (r, $\mathrm{\theta }$ , z) is shown in Figure 1, cited firstly in the introduction section showing schematic diagram of the problem studied. The local area coordinates are L₁, L₂, and L₃.

Due to the lack of analytical solutions, one will have to resort to numerical methods. The finite element method (FEM) will be used for this purpose.

2.1. Governing equations

Continuity equation

$$\mathit{\nabla }\mathit{.}\mathit{V}=\frac{\partial v_r}{\partial r}+\frac{\partial v_z}{\partial z}+\frac{v_r}{r}=0$$ E9

Motion equations

$$\rho \frac{\mathit{D}\mathit{V}}{\mathit{Dt}}=\mathit{\nabla }\mathit{.}\mathit{\tau }$$ E10

They have two projections in the axisymmetric problem considered. Hence, there are overall three equations. Unknowns are v_r, v_z, p. Supposing a creeping flow, the left-hand side may be taken null, and the motion equations reduce to

$$\mathit{\nabla }\mathit{.}\mathit{\tau }=\mathit{0}$$ E11

According to the assumptions of steady, axisymmetric, incompressible and irrotational flow properties above mentioned, we can write

$$\frac{\partial }{\partial t}=0\frac{\partial }{\partial \theta }=0\mathit{\nabla }\mathbf{.}\mathit{V}=0{\mathrm{\nu }}_{\mathrm{\theta }}=0$$ E12

2.2. Nondimensionalization

Nondimensional variables and scales are introduced as

$$r^{\ast }=\frac{r}{a}\hspace{1em}{z}^{\ast }=\frac{z}{a}\hspace{1em}{v_r}^{\ast }=\frac{v_r}{V_s}\hspace{1em}{v_{\theta }}^{\ast }=\frac{v_{\theta }}{V_s}\hspace{1em}{v_z}^{\ast }=\frac{v_z}{V_s}$$

$$p^{\ast }=p\frac{a}{V_s{\eta }_o}{t}^{\ast }==t\frac{{\eta }_o}{a^2\rho }{\eta }^{\ast }=\frac{\eta }{{\eta }_o}{{\upsilon }_1}^{\ast }=\frac{{\upsilon }_1}{{\upsilon }_{10}}{{\upsilon }_2}^{\ast }=\frac{{\upsilon }_2}{{\upsilon }_{20}}{D}^{\ast }=\frac{D}{{\eta }_0{\mathit{aV}}_s}$$

$${\left({\mathit{A}}_{\mathit{1}}\right)}^{\ast }=\left({\mathit{A}}_{\mathit{1}}\right)\frac{a}{V_s}{\left({\mathit{A}}_{\mathit{1}}^{\mathit{2}}\right)}^{\ast }=\left({\mathit{A}}_{\mathit{1}}^{\mathit{2}}\right){\left(\frac{a}{V_s}\right)}^2{\left({\mathit{A}}_{\mathit{2}}\right)}^{\ast }=\left({\mathit{A}}_{\mathit{2}}\right){\left(\frac{a}{V_s}\right)}^2{\psi }^{\ast }=\frac{\psi }{a_2{V}_s}$$ E13

For brevity * notation will be elected by implication.

2.3. Explicit expressions of the dimensionless governing equations

Continuity equation

$$\frac{\partial v_r}{\partial r}+\frac{\partial v_z}{\partial z}+\frac{v_r}{r}=0$$ E14

Projection of the motion equation on the r axis

$$\begin{array}{l}-\frac{\partial p}{\partial r}+\frac{\partial \eta }{\partial r}{\left(A_1\right)}_{\mathit{rr}}+\frac{\partial \eta }{\partial z}{\left(A_1\right)}_{\mathit{zr}}+\eta \left[\frac{\partial {\left(A_1\right)}_{\mathit{rr}}}{\partial r}+\frac{\partial {\left(A_1\right)}_{\mathit{zr}}}{\partial z}+\frac{{\left(A_1\right)}_{\mathit{rr}}-{\left(A_1\right)}_{\mathit{\theta \theta }}}{r}\right]+\\ +0.85K\left\{\frac{\partial {\upsilon }_1}{\partial r}{\left(A_1^2\right)}_{\mathit{rr}}+\frac{\partial {\upsilon }_1}{\partial z}{\left(A_1^2\right)}_{\mathit{zr}}+{\upsilon }_1\left[\frac{\partial {\left(A_1^2\right)}_{\mathit{rr}}}{\partial r}+\frac{\partial {\left(A_1^2\right)}_{\mathit{zr}}}{\partial z}+\frac{{\left(A_1^2\right)}_{\mathit{rr}}-{\left(A_1^2\right)}_{\mathit{\theta \theta }}}{r}\right]\right\}-\end{array}$$ E15

$$-\frac{1}{2}K\left\{\frac{\partial {\upsilon }_1}{\partial r}{\left(A_2\right)}_{\mathit{rr}}+\frac{\partial {\upsilon }_1}{\partial z}{\left(A_2\right)}_{\mathit{zr}}+{\upsilon }_1\left[\frac{\partial {\left(A_2\right)}_{\mathit{rr}}}{\partial r}+\frac{\partial {\left(A_2\right)}_{\mathit{zr}}}{\partial z}+\frac{{\left(A_2\right)}_{\mathit{rr}}-{\left(A_2\right)}_{\mathit{\theta \theta }}}{r}\right]\right\}=0$$

Projection of the motion equation on the z axis

$$\begin{array}{l}-\frac{\partial p}{\partial z}+\frac{\partial \eta }{\partial r}{\left(A_1\right)}_{\mathit{rz}}+\frac{\partial \eta }{\partial z}{\left(A_1\right)}_{\mathit{zz}}+\eta \left[\frac{\partial {\left(A_1\right)}_{\mathit{rz}}}{\partial r}+\frac{\partial {\left(A_1\right)}_{\mathit{zz}}}{\partial z}+\frac{{\left(A_1\right)}_{\mathit{rz}}}{r}\right]+\\ +0.85K\left\{\frac{\partial {\upsilon }_1}{\partial r}{\left(A_1^2\right)}_{\mathit{rz}}+\frac{\partial {\upsilon }_1}{\partial z}{\left(A_1^2\right)}_{\mathit{zz}}+{\upsilon }_1\left[\frac{\partial {\left(A_1^2\right)}_{\mathit{rz}}}{\partial r}+\frac{\partial {\left(A_1^2\right)}_{\mathit{zz}}}{\partial z}+\frac{{\left(A_1^2\right)}_{\mathit{rz}}}{r}\right]\right\}-\end{array}$$ E16

$$-\frac{1}{2}K\left\{\frac{\partial {\upsilon }_1}{\partial r}{\left(A_2\right)}_{\mathit{rz}}+\frac{\partial {\upsilon }_1}{\partial z}{\left(A_2\right)}_{\mathit{zz}}+{\upsilon }_1\left[\frac{\partial {\left(A_2\right)}_{\mathit{rz}}}{\partial r}+\frac{\partial {\left(A_2\right)}_{\mathit{zz}}}{\partial z}+\frac{{\left(A_2\right)}_{\mathit{rz}}}{r}\right]\right\}=0$$

where K is the normalization coefficient.

Taking as numerical example

a/R = 0.2, a = 0.05 m, and V_s = 0.016 m/s, and because for $\dot{\gamma }\to 0$ , we have $\frac{{\upsilon }_1}{\eta }\to \frac{{\upsilon }_{10}}{{\eta }_0}=1$

which ensues

$$K=\frac{{\upsilon }_{10}}{{\eta }_o}\frac{V_s}{a}=0.32$$ E17

The coefficient K is the special case for $\dot{\gamma }$ = 0 of the Weissenberg number. It is worthwhile to notice that for K = 0 we go back to the generalized Newtonian fluid.

2.3.1. Stream function

$$\psi =\text{constant}-{\int }_0^r{\mathit{rv}}_z\mathit{dr}\approx \pi -{\int }_0^r{\mathit{rv}}_z\mathit{dr}$$ E18

2.3.2. Boundary conditions (dimensionless)

At the inlet of the flow: v_r = 0 v_z = 1.

Along the cylindrical tube (r = R): v_r = 0 v_z = 1.

Along the centerline of the cylindrical tube (r = 0): v_r = 0 $\frac{\partial v_z}{\partial r}=0$

On the surface of the sphere: v_r = v_z = 0.

At the outlet of the flow: v_r = 0 $\frac{\partial v_z}{\partial r}=0$ .

$$p=0\left(\text{atmospheric pressure}\right)$$ E19

2.4. Dimensionless stress components

Introducing $K=\frac{{\upsilon }_{10}}{{\eta }_0}\frac{V_s}{a}$ in the CEF constitutive equation Eq. (1), using Eq.(13) nondimensional formulas and Eq. (7), ${\tau }^{\ast }=\tau \frac{a}{V_s{\eta }_o}$ stress tensor components can be written electing the * notation by implication as follows:

$$\begin{array}{l}{\tau }_{\mathit{rr}}=-p+\eta {\left(A_1\right)}_{\mathit{rr}}+K\left[0.85{\upsilon }_1{\left(A_1^2\right)}_{\mathit{rr}}-\frac{1}{2}{\upsilon }_1{\left(A_2\right)}_{\mathit{rr}}\right]\\ {\tau }_{\mathit{rz}}=\eta {\left(A_1\right)}_{\mathit{rz}}+K\left[0.85{\upsilon }_1{\left(A_1^2\right)}_{\mathit{rz}}-\frac{1}{2}{\upsilon }_1{\left(A_2\right)}_{\mathit{rz}}\right]\\ {\tau }_{\mathit{\theta \theta }}=-p+\eta {\left(A_1\right)}_{\mathit{\theta \theta }}+K\left[0.85{\upsilon }_1{\left(A_1^2\right)}_{\mathit{\theta \theta }}-\frac{1}{2}{\upsilon }_1{\left(A_2\right)}_{\mathit{\theta \theta }}\right]\\ {\tau }_{\mathit{zz}}=-p+\eta {\left(A_1\right)}_{\mathit{zz}}+K\left[0.85{\upsilon }_1{\left(A_1^2\right)}_{\mathit{zz}}-\frac{1}{2}{\upsilon }_1{\left(A_2\right)}_{\mathit{zz}}\right]\end{array}$$ E20

3. Application of the finite element method to fluids

Due to the assumptions Eq. (12), only v_r, v_z velocity components and pressure p are chosen as primitive variables.

All stress variables are expressed in terms of the velocities and the pressures and the stresses are then back-calculated after finding the velocities v_r, v_z and the pressure p.

For stability the pressure field must be interpolated with a polynomial one order lower than the velocity terms. Here we have chosen linear pressure and quadratic velocity fields over the element.

A quadratic triangle for velocities with six nodes, three on the vertices and the three others on mid-sides, and a linear triangle for pressures, both processed with area coordinates (Figure 5).

Figure 5.

(a) A linear triangular element; and (b) a quadratic triangular element.

In our problem gravitational and inertial forces are neglected and the motion is obtained as in the first Stokes problem by drawing the cylinder opposite to the sphere displacement with a constant, upward V_s velocity. To solve the nonlinear equations iterative optimization techniques are used.

3.1. Linear and quadratic triangles

The $\left\{N^e\right\}$ shape functions for velocities corresponding to the three vertice nodes and three mid-side nodes (1, 2, …, 6) are respectively using the well-known area coordinates L₁, L₂, L₃.

$$L_{\mathit{1}}\left(\mathit{2}{L}_{\mathit{1}}\right);{\mathit{\text{L}}}_{\mathit{2}}\left(\mathit{2}{L}_{\mathit{2}}\right);L_{\mathit{3}}\left(\mathit{2}{L}_{\mathit{3}}\right);{\mathit{\text{4L}}}_{\mathit{1}}{\mathit{\text{L}}}_{\mathit{2}};{\mathit{\text{4L}}}_{\mathit{2}}{\mathit{\text{L}}}_{\mathit{3}};{\mathit{\text{4L}}}_{\mathit{3}}{\mathit{\text{L}}}_{\mathit{1}}$$ E21

and those concerning the pressure for nodes 1, 2, 3 are L₁, L₂, L₃, (Figure 6).

Figure 6.

Area coordinates.

The element shown satisfies the LBB condition and thus gives reliable and convergent solutions for velocity and pressure fields [16].

The derivatives of the interpolation functions with respect to the global coordinates are of the form

$$\left\{\begin{array}{l}\frac{\partial N_i^e}{\partial r}\\ \frac{\partial N_i^e}{\partial z}\end{array}\right\}=J^{-1}\left\{\begin{array}{l}\frac{\partial N_i^e}{\partial L_1}\\ \frac{\partial N_i^e}{\partial L_2}\end{array}\right\}$$ E22

with the Jacobian matrix of transformation

$$J=\left|\begin{array}{l}\frac{\partial r}{\partial L_1}\frac{\partial z}{\partial L_1}\\ \frac{\partial r}{\partial L_2}\frac{\partial z}{\partial L_2}\end{array}\right|$$

The element area

$$\mathit{\text{drdz}}=\frac{1}{2}\text{det}{\mathit{JdL}}_1{\mathit{dL}}_2$$ E23

The integration must be carried out over the elemental volume of the axisymmetric geometry rdrdθdz. Since the solution is independent of the θ coordinate, the integration with respect to θ yields a multiplicative constant 2π. (Table 1) [17].

Table 1.

Quadrature weights and points for triangular elements.

3.2. Interpolation formulas

According to the constraint

$$L_1+L_2+L_3=1$$

we eliminate the dependent variable L3 for the sake of simplicity. Thus, we have

$$\begin{array}{l}{v}_r=\left(-L_1+2L_1^2\right)v_{r1}+\left(-L_2+2L_2^2\right)v_{r2}+(1-3L_1-\\ -3L_2+2L_1^2+4L_1{L}_2+2L_2^2)v_{r3}+4L_1{L}_2{v}_{r4}+\\ +4\left(L_2-L_1{L}_2-L_2^2\right)v_{r5}+4\left(L_1-L_1{L}_2-L_1^2\right)v_{r6}\end{array}$$

$$\begin{array}{l}{v}_z=\left(-L_1+2L_1^2\right)v_{z1}+\left(-L_2+2L_2^2\right)v_{z2}+(1-3L_1-\\ -3L_2+2L_1^2+4L_1{L}_2+2L_2^2)v_{z3}+4L_1{L}_2{v}_{z4}+\\ +4\left(L_2-L_1{L}_2-L_2^2\right)v_{z5}+4\left(L_1-L_1{L}_2-L_1^2\right)v_{z6}\end{array}$$

$$p=L_1\left(p_1-p_3\right)+L_2\left(p_2-p_3\right)+p_3$$ E24

3.3. Basis of application

The flow domain is meshed using linear and quadratic triangular elements [9]. Three unstructured meshes are generated by an adaptive mesh generator, as given in Table 2 and shown in Figure 7. Independent unknown variables are, due to axisymmetry vr, vz and p. The continuity equation being interpreted as an additional relation among the velocity components, i.e., a constraint satisfied in an approximate least squares sense [18], it may be omitted in the effective area consideration around the mid-side nodes, in order to keep the balance between the equation and unknown variable numbers. Resolving the global equation system obtained while taking into account the boundary conditions and integrating numerically by means of Gauss quadrature over the effective area around each node, the values of the variables are found [19, 20].

Mesh	Sphere surface nodes	Elements	Nodes
AM1	21	644	1401
AM2	31	948	2035
AM3	41	1242	2645

Table 2.

Summary of finite element meshes.

Figure 7.

Mesh patterns around sphere, a/R = 0.2.

As a typical numerical example after [9], the values below are processed:

$$a=0.05\mathrm{m} a/R=0.2 \mathit{Vs}=0.016\mathrm{m}/\mathrm{s}$$

3.3.1. Comparison and contour patterns

In order to provide a basis for a comparison, and to pose the behavioral difference between inelastic and viscoelastic fluids, an example taken from the literature [9], for the simpler case ${\mathrm{\upsilon }}_1$ = 0, gives contours of streamfunction, velocity, stress, pressure and viscosity at parameter settings of n = 0.5 and We = Weissenberg number = 2.5 (Figure 8). As it can be seen in the Figure 8, the flow accelerates as the fluid approaches the sphere.

Figure 8.

Contours of streamfunction (ψ), radial velocity (v_r), axial velocity (v_z), pressure (p), stress components (τ_rr, τ_zz, τ_θθ, τ_rz), and viscosity (η); a/R = 0.2, We = 2.5, and n = 0.5.

The extrema contour levels for the inelastic fluid are again shown in a table (Table 3).

Flow domain	Minimum	Maximum
ψ	0.0000	3.1250
νr	−0.1933	0.2171
νθ	—	—
νz	0.0000	1.1249
P	−4.9296	7.8010
τrr	−2.6545	3.1418
τθθ	−1.4139	1.8356
τrz	−1.3276	2.9703
τzz	−3.6693	2.8291
η	0.244	1.4672

Table 3.

Extrema contour levels.

3.3.2. The attempt of getting ${\mathrm{\upsilon }}_1={\mathrm{\upsilon }}_2=0$ condition as a special case

The reasons why the results obtained by Townsend [9] cannot be achieved for the generalized fluid by replacing ${\upsilon }_1={\upsilon }_2=0$ in our equations given above for the generalized condition are explained by Kaloni [21] in the literature.

Nonetheless, we think that is beneficial to make a comparison with Townsend [9] solutions.

3.3.3. Mesh-independent results

The model used in this study, which is called as Mesh#1, was obtained as a result of a three-stage optimization. In the “Initmesh” stage, an “adaptive mesh” was formed. In the second stage number of triangles was increased with the “refinemesh” function to obtain an improvement in this first mesh. Finally, at the third stage, an optimization was achieved with the “jigglemesh” function. Therefore, Mesh#1 is an optimal mesh, and it can be considered sufficient in terms of “mesh independence” alone. Mesh#2 is an adaptive mesh, and the number of triangles in this mesh, was increased considerably with respect to Mesh#1 (Figure 9 and Table 4).

Figure 9.

Mesh#1 and Mesh#2.

Definition	Matrices for points; P	Matrices for edges; E	Matrices for triangles; T	Vector of nodal values; V
Mesh#1	2 / 205	7 / 62	4/346	1/1715
Mesh#2	2 / 311	7 /84	4/536	1/2625

Table 4.

Mesh comparison.

Because the subject of this study is based on the fact that the elastic effect in polymeric fluids cannot be neglected, Mesh#1 and Mesh#2 results in Table 5 were presented as comparison of the normal stresses. In this table, the average difference between each mesh is approximately 2%. Because the results of two very different meshes mostly overlap, the results obtained are proven to be independent from mesh configuration, i.e., “mesh independence” is achieved. The generated mesh is improved according to the initmesh, refine mesh and jigglemesh MATLAB programs and the mesh independence is tested and ensured.

Stresses	τrr		τθθ		τzz
Extreme values	Min	Max	Min	Max	Min	Max
Mesh#1	−6.0019	3.3792	−6.0115	3.4000	−6.1339	3.7165
Mesh#2	−5.9810	3.2325	−5.8722	3.4022	−6.3396	3.8825

Table 5.

Comparison of normal stresses.

3.4. Explanatory flow characteristics remarks and results of our the theoretical problem solved numerically, cited before the sections findings and discussion and conclusions to better find out the problem studied

3.4.1. Explanatory remarks

The starting example of “Settling of Small Particles in Fluid Medium” problem in the literature which has important applications in many fields such as the shelf life of foodstuffs etc., is given for “generalized Newtonian fluid” medium [9]. As is known, the motion equation is Navier-Stokes in this case, and viscosity coefficient $\mathrm{\eta }$ is function of the shear-rate $\dot{\gamma }$ . Although this paper contains useful information, especially for velocity and stress contours, the elastic effects should not be neglected because in practice the fluids encountered in this field are non-Newtonian. This study is to cover this difference. Thus, CEF model was taken as fluid, and for emphasizing the elastic effect, the normal stress coefficients ${\mathrm{\upsilon }}_1$ and ${\mathrm{\upsilon }}_2$ was included. The simulation model is in the form of the flow of a non-Newtonian fluid around a sphere falling along the centerline of a cylindrical tube. The problem was tried to be solved numerically by writing continuity and motion equations and using FEM. The partial differential equation system obtained was transformed into a set of simultaneous, nonlinear, algebraic equations. Here, we chose linear triangular elements for pressures and quadratic triangular elements for velocities.

Because of the boundary conditions, the equations are “overdetermined.” The problem is axisymmetric, the global coordinates r, z were transformed into local coordinates L1, L2, L3; and Gauss quadrature was used for integrating numerically over the effective area around each node. In the solution of the algebraic equation system, naturally optimized formulation with the aid of (MATLAB-Optim set f-solve) function was used instead of the methods like “weak form.”

It is useful to give following explanatory remarks:

Although the differential equation is third order, the shape functions are chosen as linear for pressures and quadratic (2nd degree) polynomials for velocities. The degree difference between pressures and velocities is compulsory for stability. These are thought to be sufficient in the literature examples. Using third order polynomial would extend the problem extremely, and make the solution unreachable. Considering the elastic effect ( ${\mathrm{\upsilon }}_1$ , ${\mathrm{\upsilon }}_2$ coefficients) made the problem already harder compared to the generalized Newtonian fluid example [9]. Zienkiewicz especially stresses that FEM application of very refined and sophisticated models does not always yield better solutions [19].

Because the motion is relative, the sphere was taken constant in the problem; however, a motion from bottom to top with a constant Vs velocity was introduced to the cylinder, as done in the 1st Stokes problem. In this case, the inlet is through the bottom surface, while the outlet is through the top surface of the cylinder.

As an example to viscoelastic fluid, Carreau type Separan 30 was chosen [7]. The coefficient $\lambda$ = 8.04 s was taken from the table given by Tanner [7]. This fluid was used frequently in experiments, and detailed information about its properties is given in the literature [4].

The boundary conditions are generally the conditions at the inlet and outlet, and—for velocities—the no slip conditions at the sphere surface and cylinder wall. It was assumed that vr = 0 vz = 0 at sphere surface and vr = 0 vz = 1 (i.e., Vs) at cylinder wall. These values are observed clearly on the contours (Figures 10 and 11).

Figure 10.

Radial velocity (v_r) [−0.1772 −0.1333 −0.0893 −0.0454 −0.0014 0.0425 0.0865 0.1304 0.1743 0.2183].

Figure 11.

Axial velocity (v_z) [0.0000 0.1802 0.2737 0.3671 0.4606 0.5541 0.6476 0.7411 0.8345 0.9280 1.0215].

Contour values were found by dividing the space between the max and min contours by ten, and were written on the figures without covering them up. That is why contour labels are generally fractional. If an integer contour (e.g., zero) is to be drawn, it can be shown with linear interpolation on the figure.

As seen in pressure (p) contours (Figure 12), similarly to the example Rameshwaran, Townsend [9], max and min values are on the “stagnation” points on the sphere (at the vicinity of the vertical axis). Here vr = vz = 0; since the total energy is constant (Bernoulli at inviscid flow), if the kinetic energy is zero, potential energy is an extremum, i.e., the pressure is at its extremum value. This is what is seen perspectively in the 3-D surf (surface) drawing (Figure 13).

Figure 12.

Pressure (p) [−3.4962 −2.4411 −1.3860 −0.3310 0.7241 1.7792 2.8343 3.8893 4.9444 5.9995].

Figure 13.

3-D surf drawing of pressure.

The extreme values concerning the contours are gathered in Table 6. The flow is not uniform: the existence of the sphere at the center, the wall effect of cylinder (drag) and the boundary conditions prevent the flow from being uniform. Thus, max and min values are obtained.

Flow domain	Minimum	Maximum
ψ	0.0000	3.1416
νr	−0.1772	0.2183
νz	0.0000	1.0215
p	−3.4962	5.9995
τrr	−6.0019	3.3792
τzz	−6.1339	3.7165
τθθ	−6.0115	3.4000
τrz	−0.3701	0.8390
η	0.0845	1.0000
υ1	0.0029	0.9981

Table 6.

Extrema contour levels.

In the adaptive “mesh” drawn, there are 755 nodes. At these nodes vr, vz, p values are taken as unknowns and they were determined as a result of calculations. Because the contours followed the nodes, zigzag course was obtained instead of a smooth curve. If the number of nodes had been increased extremely, we could have got curves instead of the zigzag course. This is not a drawback in terms of the results. 10 contours were drawn by dividing space between the max and min values into 10 fragments. Contours are the geometric milieu of the nodes which have the same value of the related independent variable. Horizontal projections of the curves (level curves) that are obtained by intersecting the surface shown in 3D in the surf drawing with the horizontal planes yield the contours. Contours were drawn with the pdecont, pdemesh, pdesurf commands at MATLAB. These 10 values, below the contours are given separately.

3.4.2. Technical details for the solution of the equation system

We transformed the partial equation system with the aid of FEM to nonlinear overdetermined algebraic equation system. As we have numerous simultaneous and highly nonlinear equations, which moreover are overdetermined due to the existence of the boundary conditions, we were compelled to resort to the optimization techniques to resolve the equation system [22].

The end of the iterative optimization process is determined according to the fulfillment of the Error computation, Standard deviation, Histogram and minimum and maximum values of the minimized equation vector F, Euclidean norm and Infinite norm (Tables 7–9). We determined that the results are in the acceptable order. We used the (Optim set f-solve) function in MATLAB for iterative optimization. The basis of the method is based on “least squares method.”

Average error	0.1004
Standard Deviation	0.1820
η_avg	0.1104
C	1.6502

Table 7.

Error computation, standard deviation, the average of η values at sphere surface (η_avg) and drag coefficient (C).

Error spaces	0–0.5	0.5–1.0	1.0–1.5	1.5–2.0	2.0–2.5
Equation number	1637	66	11	1	0
Percentage (%)	95.4519	3.8484	0.6414	0.0583	0

Table 8.

Histogram results concerning equation errors; 1715 equations were used in total.

Min., Max. and Norms:	Min(F)	Max(F)	||F||2	||F||∞
	−1.3449	1.6251	8.6075	1.6251

Table 9.

Euclidean and infinite norms and minimum and maximum values of the minimized equation vector F.

3.4.3. Numerical solution

Obtained mesh configuration contains (linear and quadratic) 755 nodes in total (Figure 14). Totally, 1715 equations are used. Extrema contour levels are given in Table 6. Error computation, Standard deviation, Drag coefficient in Table 7 and Histogram results about equation errors are given in Table 8. Minimum and maximum values of the minimized equation vector F, Euclidean and Infinite norms are shown in Table 9.

Figure 14.

Mesh scheme for a/R = 0.2.

Stream function $\psi$ , radial velocity vr, axial velocity vz, pressure p, normal stress ${\tau }_{\mathit{rr}}$ , axial stress ${\tau }_{\mathit{zz}}$ , axisymmetric stress ${\tau }_{\mathit{\theta \theta }}$ , shear stress ${\tau }_{\mathit{rz}}$ , viscosity coefficient $\mathrm{\eta }$ , normal stress coefficient v1 contours and error histogram are given in Figures 10–13 and in Figures 15–27.

Figure 15.

Stream function (Ψ) [0.00 0.01 0.05 0.20 0.44 0.77 1.25 1.88 2.67 3.14].

Figure 16.

3-D mesh drawing of radial velocity.

Figure 17.

3-D surf drawing of radial velocity.

Figure 18.

3-D mesh drawing of axial velocity.

Figure 19.

3-D surf drawing of axial velocity.

Figure 20.

3-D mesh drawing of pressure.

Figure 21.

Normal stress (τ_rr) [−6.0019 −4.9595 −3.9172 −2.8748 −1.8325 −0.7902 0.2522 1.2945 2.3369 3.3792].

Figure 22.

Axial stress (τ_zz) [−6.1339 −5.0394 −3.9449 −2.8504 −1.7560 −0.6615 0.4330 1.5275 2.6220 3.7165].

Figure 23.

Azimuthal stress (τ_θθ) [−6.0115 −4.9658 −3.9200 −2.8743 −1.8286 −0.7829 0.2628 1.3086 2.3543 3.4000].

Figure 24.

Shear stress (τ_rz) [−0.3701 −0.2358 −0.1014 0.0329 0.1673 0.3016 0.4360 0.5703 0.7047 0.8390].

Figure 25.

Viscosity coefficient (η) [0.0845 0.1862 0.2879 0.3896 0.4914 0.5931 0.6948 0.7965 0.8982 1.0000].

Figure 26.

Normal stress coefficient (υ₁) [0.0029 0.1135 0.2241 0.3346 0.4452 0.5558 0.6663 0.7769 0.8875 0.9981].

Figure 27.

Equation numbers producing error are shown on the vertical axis.

4. Findings and discussion

In order to make a comparison and to compare the findings given in the literature in connection with settling of small particles for inelastic fluids with the behaviour of viscoelastic non-Newtonian fluids having variable material coefficient that show shear-thinning due to being mostly polymeric and viscoelastic nature of aforementioned materials; based on an example given by Rameshwaran et al. [9] for generalized Newtonian fluid, results of a typical example (CEF fluid) discussed in this study are compared visually with stream function, velocities, pressures and stresses with contour curves, and numerically with extrema tables. For both fluids the contour curves are shown in Figure 8 and in Figures 10–13 and in Figures 15–27, and the contour levels are shown numerically in Table 3 and 6. The contours are drawn for the $\psi$ stream function, vr, vz velocities, p pressure, ${\tau }_{\mathit{rr}}$ , ${\tau }_{\mathit{\theta \theta }}$ , ${\tau }_{\mathit{zz}}$ , and ${\tau }_{\mathit{rz}}$ stresses, η viscosity coefficient, υ1 normal stress coefficient (Figures 10–13, 15–27). These contours and their extrema are compared with the generalized Newtonian fluid given in the literature [9].

As explained above, although the transition from system we used to a simpler example in the literature taking υ1 as equal to 0 is not completely possible, this comparison was very useful in general terms and allowed for new findings.

This comparison results can be explained as follows:

$\psi$ stream function contours are completely similar for each two fluids. The extreme values match up with each other. vr radial velocity contours are also in complete harmony to each other. Although max values being 0.2171 and 0.2183 are approximately equal, there is small difference between minimums, −0.1933 and −0.1772. At vz contours, harmony is perfect, extreme values are respectively 0.00/0.00 and 1.1249/1.0215. Because vz = 1 must be at cylinder wall as a boundary condition, the value we found is more realistic. The courses of p contours are similar. The contour given in the literature is not very detailed, however our contour is completely detailed and also it gives a better idea as it is colorful. What attracts attention is that max and min values of each contour are at the same place at source and sink points of the sphere surface. Extreme values are −4.9296/−3.4962 and 7.8010/5.9995. The pressure drops for the fluids between sphere’s front and back faces caused by sphere motion, which settles along the cylinder axis are 7.8 − 4.93 = 2.87 and 5.9995 − 3.4962 ≈ 2.50, respectively. The small difference between the pressure drops can be explained by the elastic effect.

The viscosity coefficient (η) contours are harmonic. The contour given by us is colorful and detailed, and it continues throughout the cylinder. Shear-thinning in this contour obviously manifests itself. In the vicinity of the sphere where the shear-rate $\dot{\gamma }$ is very large, η coefficient becomes smaller and approaches to zero. Extreme values are respectively 0.244/0.0845 and 1.4672/1.0000. For both fluids Carreau equation and shear-thinning curve is used, and according to this extreme values given for the generalized Newtonian fluid seem to be exaggerated. The contour of the first normal stress coefficient υ1 properly illustrates the formula it was derived from and shows a course parallel to contour η. In conformity to the relaxation time curve $\frac{{\upsilon }_1}{\eta }$ given by Barnes et al. [8], this ratio approaches to 1 for small values of $\dot{\gamma }$ and goes to zero for large values. Extrema values are 0.0029/0.9981. Likewise ratio of the υ1 and η extreme values is 0.9981/1.0000 ≈ 1.0 for small $\dot{\gamma }$ and 0.0029/0.0845 ≈ 0.034 for large $\dot{\gamma }$ .

The main differences are at the stresses. These support our expectations. The magnitudes of normal stresses ${\tau }_{\mathit{rr}},{\tau }_{\mathit{\theta \theta }},{\tau }_{\mathit{zz}}$ in absolute values are increased considerably. Thus, elastic effect appeared by inclusion of υ1 in our equations. In spite of this, the shear stress ${\tau }_{\mathit{rz}}$ in absolute value decreased.

So, the shear effect decreased relatively, while the extensional effect increased. The category of the flow slightly changed from shear flow toward extensional flow (Table 10).

Stresses	τrr		τθθ		τzz		τrz
Extreme values	Min	Max	Min	Max	Min	Max	Min	Max
Inelastic	−2.6545	3.1418	−1.4139	1.8356	−3.6693	2.8291	−1.3276	2.9703
Viscoelastic	−6.0019	3.3792	−6.0115	3.400	−6.1339	3.7165	−0.3701	0.8390

Table 10.

Comparison of the normal and shear stresses.

5. Conclusions

Since the materials under consideration are rather polymeric and consequently viscoelastic, it is compulsory to include the elastic behaviour in the analysis by taking the effects of the normal stress into account in addition to those of the shear stresses and to use a constitutive equation appropriate to this. In this study, this point of view is tried to be ensured, and in order not to neglect the viscoelastic effect, a constitutive equation such as CEF, which includes normal stress coefficients, υ1, υ2, is used.

The results of the analysis, as expressed above in the findings and discussion section, confirms this idea. Comparison of the normal and shear stresses given in Table 10 exposed an increase which reaches 3 times in absolute value in the normal stresses and in parallel a decrease which reaches 4 times in the shear stresses. This table reveals that elastic effect in polymeric fluids cannot be neglected. In this manner a fact has been emphasized, and an important point in fluid literature has been clarified.