**2. The general model**

The concept of "shape factor" herein used is applied to a function ܩ of spatial coordinates, such as when ܩൌͲ, a series of closed curves are determined for a range of some parameters contained in ܩ. One typical example is

$$G = 1 - r^2 + \varepsilon r^n \sin n\,\theta \tag{1}$$

In this ሺݎǡ ߠሻ are polar coorfinates, ݊ is an integer number and ߝ is a parameter such as that for ߝൌͲ the curve described by (1) is a circle, and as ߝ increases, the shape evolves to some limiting shape, controlled by ݊. In all cases here considered, the maximum allowable value of ߝ is less than unity, and beyond that value, the curve is no longer a closed one. If ߝൌߝ is the critical, or maximum, allowable value of ߝ, then for the shape factor described by (1), ߝ is found to be


A more general shape factor is

$$G = 1 - r^2 + \varepsilon\_1 \, r^{n\_1} \sin n\_1 \, \theta + \varepsilon\_2 \, r^{n\_2} \sin n\_2 \theta + \dotsb \tag{3}$$

which leads to more complex shapes. Some instances of these shapes are shown in Fig. 1. For the purposes of this presentation, a general shape factor in polar coordinates can be defined as

$$G = h\_0(r) + \varepsilon\_1 h\_1(r, \theta) + \varepsilon\_2 h\_2(r, \theta) + \dotsb \tag{4}$$

in which ݄ଵ, ݄ଶ … are boundary perturbation functions. For the case of channel flow, the structure of (4) may be the same, in which polar coordinated may be substituted by Cartesian coordinates. The specific characteristics of functions ݄ are determined by the nature of the equations of motion and associated boundary conditions.

Two relevant cases can be highlighted, namely, flow with one velocity component, and flow with more than one velocity component.

$$\mathbf{w}(r,\theta,t) = \mathcal{G}(r,\theta)[f\_0(r,t) + \varepsilon f\_1(r,\theta,t) + \varepsilon^2 f\_2(r,\theta,t) + \cdots] \tag{5}$$

$$\Psi = G^2 \{ g\_0(r, t) + \varepsilon \left. g\_1(r, \theta, t) + \varepsilon^2 \left. g\_2(r, \theta, t) + \cdots \right| \right. \tag{6}$$

$$
\Delta \frac{\partial \boldsymbol{w}}{\partial t} - \left(\frac{1}{r} \frac{\partial \boldsymbol{w}}{\partial r} + \frac{\partial^2 \boldsymbol{w}}{\partial r^2} + \frac{1}{r^2} \frac{\partial^2 \boldsymbol{w}}{\partial \theta^2}\right) = \boldsymbol{\phi}(\mathbf{t}) = -\frac{\partial \boldsymbol{p}}{\partial \mathbf{z}}(\mathbf{t})\tag{7}
$$

$$\frac{\partial \mathbf{w}}{\partial \mathbf{z}} = \mathbf{0} \tag{8}$$

$$
\Omega = \frac{\rho a^2}{\mu \mathbf{T\_0}} \tag{9}
$$

$$\mathcal{W} = A\_2(1 - r^2) + A\_4(1 - r^4) + A\_6(1 - r^6) + \dotsb \tag{10}$$

$$\log\_{0} = (1 - r^2) \left\{ A + \frac{\Omega}{4^2} \frac{dA}{dt} (1 + r^2) + \frac{\Omega^2}{4^2 \delta^2} \frac{d^2 A}{dt^2} (1 + r^2 + r^4) + \cdots \right\} \tag{11}$$

$$\frac{d\phi}{dt} = A + \frac{\Omega}{2^2} \frac{dA}{dt} + \frac{\Omega^2}{2^2 \Phi^2} \frac{d^2 A}{dt^2} + \dotsb \tag{12}$$

$$f\_0 = A + \frac{\Omega}{\mathfrak{s}^2} \frac{dA}{dt} (1 + r^2) + \frac{\Omega^2}{\mathfrak{s}^2 \mathfrak{s}^2} \frac{d^2 A}{dt^2} (1 + r^2 + r^4) \cdots \tag{13}$$

$$\mathcal{L} = \Omega \frac{\partial}{\partial \mathbf{t}} - \left( \frac{1}{\mathbf{r}} \frac{\partial}{\partial \mathbf{r}} + \frac{\partial^2}{\partial \mathbf{r}^2} + \frac{1}{\mathbf{r}^2} \frac{\partial^2}{\partial \theta^2} \right) \tag{14}$$

$$L\{f\_1(1-r^2)\} = L\{f\_0 r^n \sin(n\theta)\}\tag{15}$$

$$f\_1 = \sin(n\theta) \,\Sigma\_{l=1}^{\ll} \left\{ \Omega^{\stackrel{l}{}} \frac{\mathrm{d}^{\dagger}\mathrm{A}}{\mathrm{d}t^{\dagger}} \right\} \Sigma\_{s=0}^{2l-2} \left\{ \mathcal{C}\_{ls} r^{l+s} \right\} \tag{16}$$

$$\mathcal{L}\_{10} = \frac{n-3}{16(n+1)}\tag{17}$$

$$\mathcal{C}\_{22} = \frac{72\mathcal{C}\_{10} + n - 2.5}{576(n+2)} \tag{18}$$

$$\mathcal{C}\_{20} = \frac{(n+1)(576\mathcal{C}\_{22} + 1) - 144\mathcal{C}\_{10} - 9}{576(n+1)}\tag{19}$$

### **3.1.2 Steady plastic flow**

In this application it is considered steady flow of a Bingham plastic. Here (5) is also applicable, and the equation of motion, in terms of shear stress, as defined below, is

$$\frac{\mathbf{r}\_{rz}}{r} + \frac{\partial \mathbf{r}\_{rz}}{\partial r} + \frac{1}{r} \frac{\partial \mathbf{r}\_{\theta x}}{\partial \theta} = \boldsymbol{\phi} \tag{20}$$

$$
\tau\_{rz} = -\left(1 - \frac{N}{l}\right) \frac{\partial \nu}{\partial r} \tag{21}
$$

$$
\pi\_{\theta z} = -\left(1 - \frac{N}{l}\right) \frac{1}{r} \frac{\partial \mathcal{W}}{\partial \theta} \tag{22}
$$

and

$$I = \sqrt{\left(\frac{\partial \boldsymbol{\nu}}{\partial r}\right)^2 + \left(\frac{1}{r}\frac{\partial \boldsymbol{\nu}}{\partial \theta}\right)^2} \tag{23}$$

is the second invariant of the rate of deformation tensor. The dimensionless yield stress is

$$N = \frac{\tau\_0 a}{\varkappa\_0 \eta\_o} \tag{24}$$

The momentum equation (20) is the standard one for parallel steady flow. Its structure has been made consistent with (21-22) and with the standard mathematical ordering of terms. The constitutive expressions (21-22) come from the applicable form of the Bingham fluid model. Defining

$$
\omega = \omega\_0 + \varepsilon \,\, w\_1 + \dotsb \tag{25}
$$

then, from (5) it follows

$$\mathbf{w}\_0 = (1 - r^2) f\_0 \tag{26}$$

$$w\_1 = \left(1 - r^2\right) f\_1 + r^n f\_0 \sin n\theta \tag{27}$$

and the following equations are found

$$\tau\_{rz} = N - \frac{\partial \boldsymbol{w}\_o(r)}{\partial r} - \varepsilon \frac{\partial \boldsymbol{w}\_1(r, \theta)}{\partial r} \tag{28}$$

$$\tau\_{\theta z} = \left(\frac{N}{\frac{\partial \mathbf{w}\_0(r)}{\partial r}} - 1\right) \frac{\varepsilon}{r} \frac{\partial \mathbf{w}\_1(r, \theta)}{\partial \theta} \tag{29}$$

$$-\frac{1}{r}\frac{\partial w\_0(r)}{\partial r} - \frac{\partial^2 w\_0(r)}{\partial r^2} = \phi - \frac{N}{r} \tag{30}$$

$$\frac{1}{r}\frac{\partial \boldsymbol{w}\_{1}(\boldsymbol{r},\boldsymbol{\theta})}{\partial r} + \frac{\partial^{2} \boldsymbol{w}\_{1}(\boldsymbol{r},\boldsymbol{\theta})}{\partial r^{2}} + \left(1 - \frac{N}{\frac{\partial \boldsymbol{w}\_{0}(\boldsymbol{r})}{\partial r}}\right) \frac{1}{r^{2}} \frac{\partial^{2} \boldsymbol{w}\_{1}(\boldsymbol{r},\boldsymbol{\theta})}{\partial \boldsymbol{\theta}^{2}} = \mathbf{0} \tag{31}$$

Equations (28-29) are the result of substituting (25) in (23) and of ordering terms in powers of ߝ through a linearization procedure. From (30-31) it is found

$$\mathcal{W}\_0(r) = N(r-1) + \frac{\phi}{4}(1-r^2) \tag{32}$$

$$w\_1(r, \theta) = A\_0(1 + n^2 \sum\_{l=1}^n \frac{(-1)^l \phi^l (n^2 - (l-1)^2)!}{2^l N^l (l^2)!} r^l) \cos(n\theta) \tag{33}$$

$$A\_0 = \frac{2^n [n^2!] N^n}{(-1)^n \phi^n [n^2 (n^2 - 1 \ \ ) (n^2 - 2^2) \dots (n^2 - (n-1)^2)]} \tag{34}$$

$$f\_0 = \frac{w\_0}{1 - r^2} \tag{35}$$

$$f\_1 = \frac{w\_1 - r^n f\_0 \sin n\,\theta}{1 - r^2} \tag{36}$$

$$\frac{1}{r}\frac{\partial (ru)}{\partial r} + \frac{\partial w}{\partial z} = \mathbf{0} \tag{37}$$

$$
\mu \frac{\partial u}{\partial r} + \mathcal{W} \frac{\partial u}{\partial z} = -\frac{2}{Re} \frac{\partial p}{\partial r} + \frac{2}{Re} \left( \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial (ru)}{\partial r} \right) + \frac{\partial^2 u}{\partial z^2} \right) \tag{38}
$$

$$
\mu \frac{\partial \boldsymbol{w}}{\partial r} + \boldsymbol{w} \frac{\partial \boldsymbol{w}}{\partial \boldsymbol{z}} = -\frac{2}{Re} \frac{\partial \boldsymbol{P}}{\partial \boldsymbol{z}} + \frac{2}{Re} \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial \boldsymbol{w}}{\partial r} \right) + \frac{\partial^2 \boldsymbol{w}}{\partial \boldsymbol{z}^2} \right) \tag{39}
$$

$$
\mu = \varepsilon u\_1(r, z) \dots \tag{40}
$$

$$
\omega = \omega\_0(r) + \varepsilon \omega\_1(r, z) \dots \tag{41}
$$

$$P = P\_0(\mathbf{z}) + \varepsilon P\_1(\mathbf{r}, \mathbf{z}) \dots \tag{42}$$

$$
\Psi = \Psi\_0(\mathbf{r}) + \varepsilon \Psi\_1(\mathbf{r}, \mathbf{z}) \dots \tag{43}
$$

$$
\hbar \omega = \frac{1}{r} \frac{\partial \omega}{\partial r} \tag{44}
$$

$$
\mu = -\frac{1}{r} \frac{\partial \Psi}{\partial z} \tag{45}
$$

$$r\_w = 1 + \varepsilon h(r) F(\mathbf{z}) \tag{46}$$

$$h = r; \qquad F(\mathbf{z}) = \sin \omega \mathbf{z} \tag{47}$$

$$\Psi = (r\_W - r)^2 \{ g\_0(r) + \varepsilon \{ g\_{11}(r) \sin(\omega z) + g\_{12}(r) \cos(\omega z) \} + 0(\varepsilon^2) + \cdots \} \tag{48}$$

$$
\Psi\_0 = (1 - r)^2 g\_0(r) = \frac{\Phi}{16}(r^2 - 1)^2 \tag{49}
$$

$$g\_0(r) = \frac{\phi}{16} \tag{50}$$

$$\Psi\_1 = (1 - r)^2 \left( g\_{11}(r) \sin(\omega z) + g\_{22}(r) \cos(\omega z) \right) + 2(1 - r) \operatorname{r} \sin(\omega z) \frac{\varphi}{16} \tag{51}$$

$$\Psi\_1 = \mathcal{H}\_1(\mathbf{r}) \ast \sin(\omega \mathbf{z}) + \mathcal{H}\_2(\mathbf{r}) \ast \cos(\omega \mathbf{z})\tag{52}$$

$$H\_1(r) = a\_0 + a\_2r^2 + a\_4r^4 + a\_6r^6 + a\_6r^9 + a\_{10}r^{10} + a\_{12}r^{12} + a\_{14}r^{14} + a\_{16}r^{16} \tag{53}$$

$$H\_2(r) = b\_0 + b\_2 r^2 + b\_4 r^4 + b\_6 r^6 + b\_8 r^9 + b\_{10} r^{10} + b\_{12} r^{12} + b\_{14} r^{14} + b\_{16} r^{16} \tag{54}$$

$$F(\mathbf{z}) = -0.028\mathbf{z}^4 + 0.434\mathbf{z}^3 - 2.156\mathbf{z}^2 + 3.43\mathbf{z} \tag{55}$$

which was transformed in a Fourier series in the range Ͳݖ ͵in terms of sine and cosine functions that allow a modeling similar, but more complex, to that already described. Examples of typical streamline and isovelocity patterns are shown in figure 7 and 8.

Fig. 7. Streamlines for Re=100 and ε =0.2.

Fig. 8. Isovelocity lines for Re=1 and ε =0.3.

In similar fashion, the a third contour presented is defined by

$$F(\mathbf{z}) = -0.03\mathbf{z}^4 - 0.045\mathbf{z}^3 + 0.405\mathbf{z}^2 + 0.42\mathbf{z} \tag{56}$$

For this case, typical isovelocity and isobaric curves are shown in figures 9 and 10.

Fig. 9. Isovelocity lines for Re=1 and ε =0.3.

$$
\pi\_{\mathfrak{X}\mathfrak{X}} = -\left(1 - \frac{N}{l}\right) \mathcal{Z}\left(\frac{\partial u}{\partial \mathfrak{x}}\right) \tag{57}
$$

$$
\pi\_{\mathcal{YY}} = -\left(1 - \frac{N}{I}\right) 2\left(\frac{\partial v}{\partial y}\right) \tag{58}
$$

$$
\tau\_{\rm yx} = \tau\_{\rm xy} = -\left(1 - \frac{N}{l}\right) \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \tag{59}
$$

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag{60}$$

$$
\mu \frac{\partial u}{\partial x} + \nu \frac{\partial u}{\partial y} = -\frac{1}{Re} \left( \frac{\partial \mathbf{r}\_{xx}}{\partial x} + \frac{\partial \mathbf{r}\_{xy}}{\partial y} \right) - \frac{1}{Re} \frac{\partial P}{\partial x} \tag{61}
$$

$$
\mu \frac{\partial \upsilon}{\partial x} + \upsilon \frac{\partial \upsilon}{\partial y} = -\frac{1}{Re} \left( \frac{\partial \tau\_{xy}}{\partial x} + \frac{\partial \tau\_{yy}}{\partial y} \right) - \frac{1}{Re} \frac{\partial p}{\partial y} \tag{62}
$$

$$
\mu(\mathbf{x}, \mathbf{y}) = u\_0(\mathbf{y}) + \varepsilon u\_1(\mathbf{x}, \mathbf{y}) \dots \tag{63}
$$

$$
\sigma(\mathbf{x}, \mathbf{y}) = \varepsilon \upsilon\_1(\mathbf{x}, \mathbf{y}) \dots \tag{64}
$$

$$P(\mathbf{x}, \mathbf{y}) = P\_0(\mathbf{x}) + \varepsilon P\_1(\mathbf{x}, \mathbf{y}) \dots \tag{65}$$

$$
\psi(\mathbf{x}, \mathbf{y}) = \psi\_0(\mathbf{y}) + \varepsilon \psi\_1(\mathbf{x}, \mathbf{y}) \dots \tag{66}
$$

$$u(\mathbf{x}, \mathbf{y}) = \frac{\partial}{\partial \mathbf{y}} \psi(\mathbf{x}, \mathbf{y}) \tag{67}$$

$$
\Delta \nu (\mathbf{x}, \mathbf{y}) = -\frac{\partial}{\partial \mathbf{x}} \psi (\mathbf{x}, \mathbf{y}) \tag{68}
$$

$$\psi = (\mathbf{y} - \mathbf{y}\_{\text{w}})(g\_0(\mathbf{y}) + \varepsilon g\_1(\mathbf{x}, \mathbf{y}) + \varepsilon^2 g\_2(\mathbf{x}, \mathbf{y}) + \cdots) \tag{69}$$

$$y\_{\mathsf{W}} = 1 + \varepsilon F(\mathsf{y}) \sin(\mathsf{ax})\tag{70}$$

$$
\psi\_0 = (\mathbf{y} - \mathbf{1})^2 g\_0(\mathbf{y}) \tag{71}
$$

$$g\_0(\mathbf{y}) = -\frac{1}{6}(2\mathbf{y} + 4 - 3N) \tag{72}$$

$$\psi\_1 = (1-\chi)^2 \left( g\_1(\mathbf{y}) \sin(\omega \mathbf{x}) + g\_2(\mathbf{y}) \cos(\omega \mathbf{x}) \right) + 2(1-\chi) F(\mathbf{y}) \sin(\omega \mathbf{z}) f\_0(\mathbf{y}) \tag{73}$$

$$\psi\_1 = A(\text{y})\sin(\omega x) + B(\text{y})\cos(\omega x) \tag{74}$$

$$A(\mathbf{y}) = a\_0 + a\_1\mathbf{y} + a\_2\mathbf{y}^2 + a\_3\mathbf{y}^3 + a\_4\mathbf{y}^4 + a\_5\mathbf{y}^5 + a\_6\mathbf{y}^6 + a\_7\mathbf{y}^7 + a\_8\mathbf{y}^8 + a\_9\mathbf{y}^9 \tag{75}$$

$$B(\mathbf{y}) = b\_0 + b\_1 \mathbf{y} + b\_2 \mathbf{y}^2 + b\_3 \mathbf{y}^3 + b\_4 \mathbf{y}^4 + b\_5 \mathbf{y}^5 + b\_6 \mathbf{y}^6 + b\_7 \mathbf{y}^7 + b\_8 \mathbf{y}^8 + b\_9 \mathbf{y}^9 \tag{76}$$

Fig. 15. Isovelocity lines for *Re*=1, *N*=0.1 and ε=0.2.

Fig. 16. Isovelocity lines for *Re*=20, *N*=0.2 and ε=0.3.

Fig. 17. Isovelocity lines for *Re=*100, *N=*0.1 and ε=0.3.

*x*

*x*

*x*

Fig. 15. Isovelocity lines for *Re*=1, *N*=0.1 and ε=0.2.

*y*

*y*

*y*

Fig. 16. Isovelocity lines for *Re*=20, *N*=0.2 and ε=0.3.

Fig. 17. Isovelocity lines for *Re=*100, *N=*0.1 and ε=0.3.

Fig. 18. Axial velocity profiles at x=0.75 for ߝ ൌ ͲǤʹ and *Re =*100.

Fig. 19. Axial velocity profiles at x=2.25 for ߝ ൌ ͲǤʹ and *Re =*100.

Fig. 20. Axial velocity profiles at x=0.75 for ߝ ൌ ͲǤ͵ and *Re =*100.

Fig. 21. Axial velocity profiles at x=2.25 for � � ��3 and *Re =*100.

### **4. Conclusion**

The method here described can lead to very accurate solutions for the velocity field and related variables such as shear stress, rate of flow and pressure in a great variety of flows in tubes and channels. Symbolic software presently available, such Maple and MathCAD make it possible to obtain and compute higher order solutions that, in some cases, may have complex algebraic structures. The fact that for all cases here considered, ie cases where � � 3, � is much less than unity (cf table 2), leads to a regular perturbation scheme that in most cases requires terms up to second order to achieve enough accuracy. The cases when ��� and ��� deserve special mention. For ��� the shape factor (1) describes an excentric circle, and for ��� an ellipse. In this last instance � is not bounded and can take any finite value, which implies that the perturbation scheme would break down if ���. So that, in this particular case, the method is limited to elliptical cross-sections of axes ratio close to unity. The method can be expanded to many more complex flow geometries. This possibility is implicit in the more general shape factor (3), which makes it necessary to develop a compound perturbation scheme, in terms of more than one perturbation parameter. The structure of the shape factor (1) determines that the analysis, especially for � � 3� is more sensible to the perturbation parameter for � � �� ie close to the wall conduit. This requires a careful analysis of series convergency which should define the order of the higher order term considered. On the other hand, in the case of flow in straight tubes, in all cases studied, in a considerable region around the conduit axis, say for � � ���� the flow variables are independent of the boundary geometry and take the values of the corresponding flow in round tubes.

### **5. Acknowledgment**

The authors acknowledge the financial support provided, at different stages of the work here presented, by FONDECYT-CONICYT and DICYT at the University of Santiago of Chile.
