Quadrature weights and points for triangular elements.
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.
- non-Newtonian fluids
- CEF equation
- shear strain-rate
- particle settling
- optimization techniques
- elastic behaviour
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, and as well as viscosity coefficient , 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:
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 is just ; the other two terms, containing and , describe the elastic effects associated with the normal stresses .
The Rivlin-Ericksen tensors involved in the CEF equation are :
The first invariant of
The shear strain-rate in term of the second invariant
The velocity gradients of
18.104.22.168. 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  for the viscosity coefficient (n = 0.364).
The normal stress coefficients may be handled as below:
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 .
Hence, we shall take
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 . 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 (
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
They have two projections in the axisymmetric problem considered. Hence, there are overall three equations. Unknowns are
According to the assumptions of steady, axisymmetric, incompressible and irrotational flow properties above mentioned, we can write
Nondimensional variables and scales are introduced as
For brevity * notation will be elected by implication.
2.3. Explicit expressions of the dimensionless governing equations
Projection of the motion equation on the
Projection of the motion equation on the
Taking as numerical example
2.3.1. Stream function
2.3.2. Boundary conditions (dimensionless)
At the inlet of the flow:
Along the cylindrical tube (
Along the centerline of the cylindrical tube (
On the surface of the sphere:
At the outlet of the flow:
3. Application of the finite element method to fluids
Due to the assumptions Eq. (12), only
All stress variables are expressed in terms of the velocities and the pressures and the stresses are then back-calculated after finding the velocities
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
3.1. Linear and quadratic triangles
shape functions for velocities corresponding to the three vertice nodes and three mid-side nodes (
and those concerning the pressure for nodes
The element shown satisfies the LBB condition and thus gives reliable and convergent solutions for velocity and pressure fields .
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
3.2. Interpolation formulas
According to the constraint
we eliminate the dependent variable
3.3. Basis of application
The flow domain is meshed using linear and quadratic triangular elements . 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
|Mesh||Sphere surface nodes||Elements||Nodes|
As a typical numerical example after , 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 , for the simpler case
= 0, gives contours of streamfunction, velocity, stress, pressure and viscosity at parameter settings of
The extrema contour levels for the inelastic fluid are again shown in a table (Table 3).
3.3.2. The attempt of getting condition as a special case
The reasons why the results obtained by Townsend  cannot be achieved for the generalized fluid by replacing in our equations given above for the generalized condition are explained by Kaloni  in the literature.
Nonetheless, we think that is beneficial to make a comparison with Townsend  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
|Matrices for edges
|Matrices for triangles
|Vector of nodal values
|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.
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 . As is known, the motion equation is Navier-Stokes in this case, and viscosity coefficient is function of the shear-rate . 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 and 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
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 ( , coefficients) made the problem already harder compared to the generalized Newtonian fluid example . Zienkiewicz especially stresses that FEM application of very refined and sophisticated models does not always yield better solutions .
Because the motion is relative, the sphere was taken constant in the problem; however, a motion from bottom to top with a constant
As an example to viscoelastic fluid, Carreau type Separan 30 was chosen . The coefficient = 8.04 s was taken from the table given by Tanner . This fluid was used frequently in experiments, and detailed information about its properties is given in the literature .
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
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 (
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.
In the adaptive “mesh” drawn, there are 755 nodes. At these nodes
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 .
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
|Min., Max. and Norms:||Min(F)||Max(F)||||F||2||||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.
, radial velocity
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.  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
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:
stream function contours are completely similar for each two fluids. The extreme values match up with each other.
The viscosity coefficient (
The main differences are at the stresses. These support our expectations. The magnitudes of normal stresses in absolute values are increased considerably. Thus, elastic effect appeared by inclusion of υ1 in our equations. In spite of this, the shear stress 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).
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.
Celasun Ş, Öztürk Y. CEF model in the industrial application of non-Newtonian fluids. Proceedings of the Parallel CFD Conference, Kyoto, Japan. May 2002. 20-22
Celasun Ş, Öztürk Y. An example to the numerical processing of nonisothermal flow of non-Newtonian fluids. Proceedings of the Parallel CFD Conference, Moscow, Russia. May 2003
Celasun Ş. Flow of Newtonian and non-Newtonian fluids. European International Journal of Science and Technology. 2014; 3:3
Bird RB, Armstrong RC, Hassager O. Dynamics of Polymeric Liquids. Vol. 1, Fluid Mechanics. U.S.A: John Wiley; 1987
Accurate solution of steady and rotational motion of Rivlin-Ericksen fluid in between two coaxial, porous, rotating cylinders [thesis]. Öztürk Yılmaz: I.T.U. Department of Mechanical Engineering, 1975
Huilgol RR. Continuum Mechanics of Viscoelastic Liquids. New York: John Wiley & Sons Inc; 1975
Tanner RI. Engineering Rheology. Oxford: Clarendon Press; 1988
Barnes HA, Hutton JF, Walters K. An Introduction to Rheology. New York: Elsevier; 1989
Rameshwaran P, Townsend P, Webster MF. Simulation of particle settling in rotating and non-rotating flows of non-Newtonian fluids. International Journal of Numerical Method Fluids. 1998; 26, 851-874, Swansea, U.K.
Bohlin T. On the drag on a rigid sphere moving in a viscous fluid inside a cylindrical tube. Transactions of Royal,Institute of Technology (Stockholm). 1960:155
Dennis SC, Ingham R, Singh SN. The slow translation of a sphere in a rotating viscous fluid. Journal of Fluid Mechanics. 1982; 117:251-267
Haberman WL, Sayre RM. Motion of rigid and fluid spheres in stationary and moving liquids inside cylindrical tubes. Hydromechanics Laboratory, Research and Development Report, No.1143, Department of the Navy, David Taylor Model Basın, 1958
Maxworthy T. An experimental determination of the slow motion of a sphere in a rotating viscous fluid. Journal of Fluid Mechanics. 1965; 23:373-384
Manogg GJ, Townsend P, Webster MF. Numerical simulation of multilayer injection moulding. Journal of Non-Newtonian Fluid Mechanics. 1997; 68:153-167
Yuan SW. Foundations of Fluid Mechanics. New Delhi: Prentice – Hall; 1969
Oden JT, Carey GF. Finite Elements, Mathematical Aspects. Vol. Vol. IV. New Jersey: Prentice Hall; 1983
Cowper GR. Gaussian Quadrature Formulas for Triangles. Ottowa: National Research Council of Canada; 2007
Reddy JN, Gartling DK. The Finite Element Method in Heat Transfer and Fluid Dynamics (Second Edition). New York: CRC Press; 2001
Zienkiewicz OC, Taylor RL. The Finite Element Method (Fourth Edition). Vol. 1. London: Mc Graw – Hill; 1994
Celasun Ş. Three-dimensional lubrication theory in viscoelastic short-bearing. Industrial Lubrication and Tribology. 2013; 65(6):pp.357-368. Emerald Group Publishing Limited
Kaloni PN. Some remarks on useful theorems for the second order fluid. Journal of Non-Newtonian Fluid Mechanics. 1989; 31:115-120
Adolph T, Schönauer W, Koch R, Knoll G. High Performance Computing in Science and Engineering. Münich: Garching; 2009