An efficient algorithm of solution for the flow of generalized Newtonian fluid in channels of simple geometries

Abstract

In this paper a problem of stationary flow of generalized Newtonian fluid in a thin channel is considered. An efficient algorithm of solution is proposed that includes a flexible procedure for a continuous approximation of the apparent viscosity by means of elementary functions combined with analytical integration of the governing equations. The algorithm can be easily adapted to circular or elliptic conduits. The accuracy and efficiency of computations are analyzed using an example of the Carreau fluid. The proposed computational scheme proves to be highly efficient and versatile providing excellent accuracy of solution at a very low computational cost.

Appendix A: Numerical solution to the problem of Carreau fluid flow

The numerical solution to the problem with the Carreau law is obtained in the following way. Just as was done in the main body of the paper, we consider only one symmetrical part of the channel y ∈ [0,w/2]. We discretize it by M uniformly spaced nodal points y_i with i = 1,...,M.

When formally integrating equation (1) with respect to y over the span [0,y] under the boundary condition (4)₁ one arrives at the relation:

$$ \eta_{a}\left( \dot \gamma\right)\dot \gamma=\frac{dp}{dx}y. $$

(35)

On the left-hand side of (35) the definition (3) was used. After simple rearrangement of equation (35) and another analytical integration from y to w/2 with the boundary condition (4)₂ we obtain the following formula for the fluid velocity:

$$ v(y)=-\frac{dp}{dx}{\int}_{y}^{w/2}\frac{y}{\eta_{\text{a}}(\dot \gamma)}dy. $$

(36)

For every point of the spatial domain, y_i, equation (35) can be treated as an algebraic equation with respect to $\dot \gamma _{i}=\dot \gamma (y_{i})$:

$$ \eta_{\text{a}}(\dot \gamma_{i})\gamma_{i}=\frac{dp}{dx}y_{i}, \quad i=1,...,M. $$

(37)

This equation is solved iteratively. The iterations are stopped when the relative difference between two consecutive values of $\dot \gamma$ falls below 10^− 12. After obtaining a discrete characteristics $\dot \gamma (y)$, it is used to approximate the integrand from equation (36) by cubic splines. Subsequent integration according to the formula yields the values of velocity at the respective spatial points y_i. The fluid flow rate, Q, is computed by integration of the velocity profile in line with formula (6).

We verify the accuracy of computations by the constructed scheme using two benchmark examples. The first one is an analytical solution obtained for the power-law fluid where the velocity profile and the fluid flow rate are expressed as (Perkowska et al. 2016):

$$ \begin{array}{@{}rcl@{}} v(y)&=&\frac{n}{n+1}\left( -\frac{1}{C}\frac{dp}{dx}\right)^{1/n}\left[\left( \frac{w}{2}\right)^{\frac{n+1}{n}}-y^{\frac{n+1}{n}} \right], \\ && Q=\frac{n}{2n+1}2^{-\frac{n+1}{n}}\left( -\frac{1}{C}\frac{dp}{dx}\right)^{1/n}w^{\frac{2n+1}{n}}. \end{array} $$

(38)

The values of C and n are taken from Table 1. The second benchmark example is the solution to the truncated power-law problem given by the formulae (15)–(17) and (18)–(20). The channel height, w, was set to 10^− 3 m. The accuracy of computations is assessed by the relative error of the fluid flow rate, δQ, and the relative error of the fluid velocity, δv.

In Fig. 8 the error dependence on the mesh density is shown (for δv the maximal value over y is taken) for two values of the pressure gradient: dp/dx = − 5 Pa/m and dp/dx = − 75 Pa/m. As can be seen, the overall accuracy is very good even with only 20 nodal points. The maximal error is of the level of 10^− 8 for the first benchmark and 10^− 6 for the second one. With M = 200 the errors are reduced to the level of 10^− 11 regardless of the considered benchmark.

Fig. 8

The relative errors of the fluid flow rate, δQ, and the maximal relative errors of the fluid velocity, δv, for a the power-law model and b the truncated power-law model

The results show that for a predefined M the accuracy is better for the first benchmark. It stems from the fact that in the second benchmark the integrand on the right-hand side of (36) is not a smooth function of y. Indeed, at the points that correspond to the limiting values of the shear rate ($\dot \gamma _{1}$ and $\dot \gamma _{2}$) η_a is only of C⁰ class. Thus, as cubic spline interpolation of the integrand does not preserve this feature, the overall accuracy is reduced.

To illustrate this problem we show in Fig. 9 the spatial distributions of the velocity errors, δv. For the second benchmark the limiting values of viscosity were reached for $y\sim 9.26 \cdot 10^{-6}$ ($\dot \gamma _{1}$) and $y\sim 1.33 \cdot 10^{-4}$ ($\dot \gamma _{2}$), respectively. Especially in the latter location a sharp error magnification can be observed (Fig. 9b). Naturally, in the first benchmark no such behavior is present as η_a is a smooth function of $\dot \gamma$ for that case. Moreover, the same situation holds for the Carreau problem. Thus, we can conclude that the numerical solution obtained by the proposed numerical scheme in the case of Carreau rheology is more accurate than the one computed for the truncated power-law.

Fig. 9

The spatial distributions of the relative error of fluid velocity, δv, for a the power-law model and b the truncated power-law model. The pressure gradient was dp/dx = − 75 Pa/m

In order to support the last conclusion let us make a comparison of the fluid flow rates, Q, obtained for the Carreau law by the following: (i) the proposed numerical scheme, (ii) the semi-analytical solution delivered in Sochi (2015). Note that the latter assumes solving a respective algebraic equation numerically to find the value of a shear rate at the conduit wall, $\dot \gamma _{\text {w}}$, and subsequent substitution to an analytical formula for Q. Naturally, the error is generated only at the first stage. When looking for $\dot \gamma _{\text {w}}$ we set the relative tolerance for this parameter to 10^− 14, which defines the accuracy of this reference solution.

The relative differences between the fluid flow rates, δQ, computed by respective methods are shown in Fig. 10. Different mesh densities were analyzed for the pressure gradient in the range 1 Pa/m ≤|dp/dx|≤ 150 Pa/m. One can note that the accuracy of the numerical solution delivered by the proposed scheme depends essentially on both, the mesh density and the pressure gradient. The average (over dp/dx) values of fluid flow rate deviation, δQ_av(M), for different M are: δQ_av(20) = 9.9 ⋅ 10^− 7, δQ_av(50) = 2.3 ⋅ 10^− 8, δQ_av(100) = 4.9 ⋅ 10^− 10, δQ_av(200) = 7.0 ⋅ 10^− 12. A comparison of these numbers and the data from Fig. 10 with the characteristics provided in Fig. 8b allows us to support the claim that the solution obtained here for the Carreau law is more accurate than the one computed for the truncated power-law benchmark. Thus, the former can be confidently adopted as a reference example when estimating the accuracy of approximate solution introduced in the section “Approximate solution for the generalized Newtonian fluid.”

Fig. 10

The relative difference between the fluid flow rates, δQ, computed with the Carreau law by means of the proposed numerical scheme and the method presented in Sochi (2015)