Quadrature weights and points for triangular elements.

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

### Keywords

- non-Newtonian fluids
- CEF equation
- shear strain-rate
- particle settling
- FEM
- optimization
- optimization techniques
- elastic behaviour

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

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,

### 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:

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

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

The first invariant of **d** is null

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

The velocity gradients of **A**_{1}, **A**_{2} 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.

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

The normal stress coefficients may be handled as below:

If *η = η*_{o} and υ_{1} = υ_{10} for *η/ηo* and

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

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

Hence, we shall take

Consequently

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*, *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*_{1}, *L*_{2}, and *L*_{3}.

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

Motion equations

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

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

### 2.2. Nondimensionalization

Nondimensional variables and scales are introduced as

For brevity * notation will be elected by implication.

### 2.3. Explicit expressions of the dimensionless governing equations

Continuity equation

Projection of the motion equation on the *r* axis

Projection of the motion equation on the *z* axis

where *K* is the normalization coefficient.

Taking as numerical example

*a/R* = 0.2, *a* = 0.05 m, and *Vs =* 0.016 m/s, and because for

which ensues

The coefficient *K* is the special case for *K* = 0 we go back to the generalized Newtonian fluid.

#### 2.3.1. Stream function

#### 2.3.2. Boundary conditions (dimensionless)

At the inlet of the flow: *vr =* 0 *vz =* 1.

Along the cylindrical tube (*r* = *R*): *vr =* 0 *vz =* 1.

Along the centerline of the cylindrical tube (*r* = 0): *vr =* 0

On the surface of the sphere: *vr = vz* = 0.

At the outlet of the flow: *vr =* 0

### 2.4. Dimensionless stress components

Introducing

## 3. Application of the finite element method to fluids

Due to the assumptions Eq. (12), only *vr*, *vz* 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 *vr*, *vz* 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).

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 *Vs* velocity. To solve the nonlinear equations iterative optimization techniques are used.

### 3.1. Linear and quadratic triangles

The *1*, *2*, …, *6*) are respectively using the well-known area coordinates *L*_{1}, *L*_{2}, *L*_{3}.

and those concerning the pressure for nodes *1, 2, 3* are *L*_{1}, *L*_{2}, *L*_{3}, (Figure 6).

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

with the Jacobian matrix of transformation

The element area

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

### 3.2. Interpolation formulas

According to the constraint

we eliminate the dependent variable *L*3 for the sake of simplicity. Thus, we have

### 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 |

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

#### 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 *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.

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 |

#### 3.3.2. The attempt of getting υ 1 = υ 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

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

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 |

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 |

### 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

Because of the boundary conditions, the equations are “overdetermined.” The problem is axisymmetric, the global coordinates *r*, *z* were transformed into local coordinates *L*1, *L*2, *L*3; 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 (

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

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

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

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 |

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 |

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 |

Min., Max. and Norms: | Min(F) | Max(F) | ||F||_{2} | ||F||_{∞} |
---|---|---|---|---|

−1.3449 | 1.6251 | 8.6075 | 1.6251 |

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

Stream function *vr*, axial velocity *vz*, pressure *p*, normal stress *v*1 contours and error histogram are given in Figures 10–13 and in Figures 15–27.

## 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 *vr*, *vz* velocities, *p* pressure, *η* 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:

*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 *η* 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 *η* extreme values is 0.9981/1.0000 ≈ 1.0 for small

The main differences are at the stresses. These support our expectations. The magnitudes of normal stresses

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 |

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