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 <sup>N</sup><sup>e</sup> f g shape functions for velocities corresponding to the three vertice nodes and three midside nodes (1, 2, …, 6) are respectively using the well-known area coordinates L1, L2, L3.

Particle Settling in a Non-Newtonian Fluid Medium Processed by Using the CEF Model http://dx.doi.org/10.5772/intechopen.75977 73

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

$$\rm L\_1(2\ L\_{1-1});\rm L\_2\ (2\ L\_{2-1});\rm L\_3(2\ L\_{3-1});\rm 4L\_1\newline L\_2;4\ L\_2\ 4\ L\_3;4\ L\_3\ \rm L\_1\tag{21}$$

and those concerning the pressure for nodes 1, 2, 3 are L1, L2, L3, (Figure 6).

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

Figure 6. Area coordinates.

Along the centerline of the cylindrical tube (r = 0): vr = 0 <sup>∂</sup>vz

<sup>∂</sup><sup>r</sup> ¼ 0.

Vsη<sup>o</sup>

<sup>τ</sup>rr ¼ �<sup>p</sup> <sup>þ</sup> <sup>η</sup>ð Þ <sup>A</sup><sup>1</sup> rr <sup>þ</sup> <sup>K</sup> <sup>0</sup>:85υ<sup>1</sup> <sup>A</sup><sup>2</sup>

τθθ ¼ �<sup>p</sup> <sup>þ</sup> <sup>η</sup>ð Þ <sup>A</sup><sup>1</sup> θθ <sup>þ</sup> <sup>K</sup> <sup>0</sup>:85υ<sup>1</sup> <sup>A</sup><sup>2</sup>

<sup>τ</sup>zz ¼ �<sup>p</sup> <sup>þ</sup> <sup>η</sup>ð Þ <sup>A</sup><sup>1</sup> zz <sup>þ</sup> <sup>K</sup> <sup>0</sup>:85υ<sup>1</sup> <sup>A</sup><sup>2</sup>

are then back-calculated after finding the velocities vr, vz and the pressure p.

<sup>τ</sup>rz <sup>¼</sup> <sup>η</sup>ð Þ <sup>A</sup><sup>1</sup> rz <sup>þ</sup> <sup>K</sup> <sup>0</sup>:85υ<sup>1</sup> <sup>A</sup><sup>2</sup>

3. Application of the finite element method to fluids

p ¼ 0 atmospheric pressure

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

2.4. Dimensionless stress components

η0 Vs

formulas and Eq. (7), <sup>τ</sup><sup>∗</sup> <sup>¼</sup> <sup>τ</sup> <sup>a</sup>

notation by implication as follows:

At the outlet of the flow: vr = 0 <sup>∂</sup>vz

Introducing <sup>K</sup> <sup>¼</sup> <sup>υ</sup><sup>10</sup>

72 Polymer Rheology

primitive variables.

element.

niques are used.

3.1. Linear and quadratic triangles

<sup>∂</sup><sup>r</sup> ¼ 0

<sup>a</sup> in the CEF constitutive equation Eq. (1), using Eq.(13) nondimensional

rr � <sup>1</sup> 2

υ1ð Þ A<sup>2</sup> rz

θθ � <sup>1</sup> 2

zz � <sup>1</sup> 2

1 

> 1

1 

rz � <sup>1</sup> 2

1 

Due to the assumptions Eq. (12), only vr, vz velocity components and pressure p are chosen as

All stress variables are expressed in terms of the velocities and the pressures and the stresses

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

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

The <sup>N</sup><sup>e</sup> f g shape functions for velocities corresponding to the three vertice nodes and three midside nodes (1, 2, …, 6) are respectively using the well-known area coordinates L1, L2, L3.

(19)

stress tensor components can be written electing the \*

υ1ð Þ A<sup>2</sup> rr

υ1ð Þ A<sup>2</sup> θθ

(20)

υ1ð Þ A<sup>2</sup> zz

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

$$\begin{Bmatrix} \frac{\partial N\_i^\epsilon}{\partial r} \\\\ \frac{\partial N\_i^\epsilon}{\partial z} \end{Bmatrix} = J^{-1} \begin{Bmatrix} \frac{\partial N\_i^\epsilon}{\partial L\_1} \\\\ \frac{\partial N\_i^\epsilon}{\partial L\_2} \end{Bmatrix} \tag{22}$$

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

<sup>1</sup> <sup>þ</sup> <sup>4</sup>L1L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

2 vr<sup>5</sup> <sup>þ</sup> <sup>4</sup> <sup>L</sup><sup>1</sup> � <sup>L</sup>1L<sup>2</sup> � <sup>L</sup><sup>2</sup>

<sup>1</sup> <sup>þ</sup> <sup>4</sup>L1L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

2 vz<sup>5</sup> <sup>þ</sup> <sup>4</sup> <sup>L</sup><sup>1</sup> � <sup>L</sup>1L<sup>2</sup> � <sup>L</sup><sup>2</sup>

<sup>þ</sup> <sup>L</sup><sup>2</sup> <sup>p</sup><sup>2</sup> � <sup>p</sup><sup>3</sup>

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

a ¼ 0:05 m a=R ¼ 0:2 Vs ¼ 0:016 m=s

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 υ<sup>1</sup> = 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

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

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

2 vr<sup>2</sup> þ ð<sup>1</sup> � <sup>3</sup>L1�

2 vz<sup>2</sup> þ ð<sup>1</sup> � <sup>3</sup>L1�

<sup>2</sup>Þvr<sup>3</sup> þ 4L1L2vr4þ

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

<sup>2</sup>Þvz<sup>3</sup> þ 4L1L2vz4þ

vr<sup>6</sup>

vz<sup>6</sup>

1

1

<sup>þ</sup> <sup>p</sup><sup>3</sup> (24)

http://dx.doi.org/10.5772/intechopen.75977

75

1 vr<sup>1</sup> þ �L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

<sup>þ</sup> <sup>4</sup> <sup>L</sup><sup>2</sup> � <sup>L</sup>1L<sup>2</sup> � <sup>L</sup><sup>2</sup>

1 vz<sup>1</sup> þ �L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

<sup>þ</sup> <sup>4</sup> <sup>L</sup><sup>2</sup> � <sup>L</sup>1L<sup>2</sup> � <sup>L</sup><sup>2</sup>

p ¼ L<sup>1</sup> p<sup>1</sup> � p<sup>3</sup>

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

Figure 8, the flow accelerates as the fluid approaches the sphere.

vr ¼ �L<sup>1</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

vz ¼ �L<sup>1</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

3.3. Basis of application

found [19, 20].

3.3.1. Comparison and contour patterns

Table 2. Summary of finite element meshes.

� <sup>3</sup>L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

� <sup>3</sup>L<sup>2</sup> <sup>þ</sup> <sup>2</sup>L<sup>2</sup>

with the Jacobian matrix of transformation

$$J = \begin{vmatrix} \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{vmatrix}$$

The element area

$$drdz = \frac{1}{2} \text{det} \mathbf{J} dL\_1 dL\_2 \tag{23}$$

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

$$L\_1 + L\_2 + L\_3 = 1$$



Table 1. Quadrature weights and points for triangular elements.

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

$$\begin{aligned} v\_r &= \left(-L\_1 + 2L\_1^2\right)v\_{r1} + \left(-L\_2 + 2L\_2^2\right)v\_{r2} + (1 - 3L\_1 - 2L\_2)v\_{r3} \\ &- 3L\_2 + 2L\_1^2 + 4L\_1L\_2 + 2L\_2^2)v\_{r3} + 4L\_1L\_2v\_{r4} + \\ &+ 4\left(L\_2 - L\_1L\_2 - L\_2^2\right)v\_{r5} + 4\left(L\_1 - L\_1L\_2 - L\_1^2\right)v\_{r6} \\\ 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\_1L\_2 + 2L\_2^2)v\_{z3} + 4L\_1L\_2v\_{z4} + \\ &+ 4\left(L\_2 - L\_1L\_2 - L\_2^2\right)v\_{z5} + 4\left(L\_1 - L\_1L\_2 - L\_1^2\right)v\_{z6} \\\ p &= L\_1(p\_1 - p\_3) + L\_2\left(p\_2 - p\_3\right) + p\_3 \end{aligned} \tag{24}$$

#### 3.3. Basis of application

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

∂N<sup>e</sup> i ∂L<sup>1</sup> ∂N<sup>e</sup> i ∂L<sup>2</sup>

∂z ∂L<sup>1</sup> � � � � � � � � �

detJdL1dL<sup>2</sup> (23)

L1 L2 L3 W Geometric locations

1/3 1/3 1/3 1

1/2 0 1/2 0 1/2 1/2 1/2 1/2 0

1/3

∂z ∂L<sup>2</sup> 9 >>>=

>>>;

(22)

8 >>><

>>>:

∂N<sup>e</sup> i ∂r ∂N<sup>e</sup> i ∂z

9 >>=

>>; ¼ J �1

∂r ∂L<sup>1</sup>

� � � � � � � � �

drdz <sup>¼</sup> <sup>1</sup> 2

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 θ

L<sup>1</sup> þ L<sup>2</sup> þ L<sup>3</sup> ¼ 1

Number of integration points Degree of polynomial Order of the residual Location of integration points

∂r ∂L<sup>2</sup>

8 >><

>>:

J ¼

with the Jacobian matrix of transformation

yields a multiplicative constant 2π. (Table 1) [17].

1 O (h2 )

2 O (h3 )

Table 1. Quadrature weights and points for triangular elements.

The element area

74 Polymer Rheology

3.2. Interpolation formulas

According to the constraint

1 (Pressures)

3 (Velocities) 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].

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

$$a = 0.05 \text{ m} \quad a/R = 0.2 \quad \text{Vs} = 0.016 \text{ m/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 υ<sup>1</sup> = 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.

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


Table 2. Summary of finite element meshes.

Figure 7. Mesh patterns around sphere, a/R = 0.2.
