**1. Introduction**

Nonlinear differential rheological equations of state are being used increasingly often to describe the rheological properties of viscoelastic fluids and to solve fluid mechanics problems related to polymer melts and solutions. First of all, the features of viscoelastic polymer-based materials behavior include stress-relaxation phenomena, the presence of elastic properties, nonlinear dependence of effective viscosity on shear rate, the occurrence of normal stresses in shear flows, and the ability to swell.

The rheological models based on relaxation equations of state and the structure of polymer molecules, which describe all the above phenomena that occur in flows of viscoelastic polymer materials in the channels of production equipment, are well known in the literature. Analytical solutions can be obtained using simple equations or under rather rough simplifying conditions. As a consequence, most of the studies deal with unimodal viscoelastic rheological models, e.g., Gruz and Pinho [1,2] (Poiseuille-Couette flows of PTT fluids), Oliveira [3] (Fene-P and Giesekus fluid flows in round pipes and flat channels), Schleiniger and Weinacht [4], or a number of works by Oliveira et al. [5] and Coelho et al. [6,7] (isothermal and non-isothermal Fene-P and PTT fluid flows in round pipes and flat channels [5,8]), Hashemabadi [9,10]. Nevertheless, the rheological properties of polymeric melts and concentrated polymer solutions are very often more complex than predicted by such unimodal models, and hence multimode models providing more adequate description are required. In [11], an analytical solution was obtained for a simplified multimode rheological model of a PTT fluid flow. In most cases, calculations of multimode viscoelastic flows are carried out by numerical methods, e.g., [12,13]. The disadvantage of the numerical approach to solution of the considered problems is that for each particular case, it is necessary to perform a full complex of numerical studies, which, due to strong nonlinearity of the considered equations, requires significant computational resources. In contrast, the parametric representation of solution yields distributions of velocities and stresses for a wide variety of problems with any degree of accuracy.

Analytical methods for viscoelastic fluid flows described in [4–6] employed simplified rheological models which led to relations connecting the nonzero components of the

**Citation:** Vachagina, E.; Dushin, N.; Kutuzova, E.; Kadyirov, A. Exact Solution for Viscoelastic Flow in Pipe and Experimental Validation. *Polymers* **2022**, *14*, 334. https:// doi.org/10.3390/polym14020334

Academic Editors: Célio Bruno Pinto Fernandes, Salah Aldin Faroughi, Luís L. Ferrás and Alexandre M. Afonso

Received: 10 December 2021 Accepted: 12 January 2022 Published: 15 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

elastic strain tensor and the shear-rate gradient. Using the solution to the momentum transfer equation in projection onto the direction of a viscoelastic medium motion for such simplified models, the authors [4–6] then obtained a correlation between a specific selected unknown function and a transverse coordinate. In this case, the remaining unknown functions are expressed in terms of the selected one using the obtained relations. Our solution method does not have the drawback of using simplified rheological models and can be adapted to an arbitrary number of the rheological model modes. Some experimental studies were carried out to visualize the viscoelastic medium flow in a round pipe to check the obtained results of theoretical studies. The experimental data were compared with the results obtained using the Giesekus model and Extended Pom-Pom.

#### **2. Problem Statement**

#### *2.1. Mathematical Model*

A steady laminar isothermal flow of viscoelastic incompressible fluid in around pipes was considered. We assumed that the velocity vector had a single axial velocity component *Vz*, which is a function of the only variable *r* of a cylindrical coordinate system *r*, *ϕ*, *z* with the z-axis directed along the pipe axis. Under the assumptions made, the system of equations for the momentum transfer and continuity in the selected coordinate system can be written as

$$-\frac{\partial P}{\partial z} + \frac{1}{r}\frac{d(r\sigma\_{rz})}{dr} = 0,\\ -\frac{\partial P}{\partial r} + \frac{1}{r}\frac{d(r\sigma\_{rr})}{dr} - \frac{\sigma\_{\theta\phi}}{r} = 0,\\ -\frac{\partial P}{\partial \phi} = 0,\\ \frac{dV\_z}{dz} = 0 \tag{1}$$

where *P* is the pressure; *σrr*, *σϕϕ*, *σrz* are physical components of the stress tensor in the cylindrical coordinate system. Let us assume that the liquid adheres to the pipe wall. It follows from (1) that

$$
\sigma\_{rz} = -|\mathbb{C}\_0|r/2,\ \hat{\text{OP}}/\partial z = \mathbb{C}\_0 = const = -|\mathbb{C}\_0|\tag{2}
$$

The rheological equation of state for a multimode viscoelastic fluid can be written as

$$\boldsymbol{\sigma} = \sum\_{k=1}^{n} \boldsymbol{\sigma}\_k + \boldsymbol{\sigma}\_{N\prime} \boldsymbol{\sigma}\_N = 2\eta\_N \boldsymbol{\mathcal{D}}$$

where *n* is the total number of modes; σ is the extra stress tensor; σ*<sup>N</sup>* is the Newtonian part of extra stress tensor; *η<sup>N</sup>* is the solvent viscosity; σ*<sup>k</sup>* is the elastic part of extra stress tensor corresponding to each mode.

Two rheological models of viscoelastic behavior were used to determine the elastic components of each mode:

Giesekus model [14]:

$$
\boldsymbol{\sigma}\_{k} + \lambda\_{k} \boldsymbol{\overline{\sigma}}\_{k} + \frac{\boldsymbol{\alpha}\_{k} \lambda\_{k}}{\eta\_{k}} \boldsymbol{\sigma}\_{k} \cdot \boldsymbol{\sigma}\_{k} = 2 \eta\_{k} \mathbf{D}\_{\prime} \left( k = 1, \ldots, n \right) \tag{3}
$$

Single-equation eXtended Pom-Pom [15] (hereinafter mentioned as eXt. Pom-Pom):

$$f(\boldsymbol{\sigma}\_k)\boldsymbol{\sigma}\_k + \lambda\_k \binom{\boldsymbol{\nabla}}{\boldsymbol{\sigma}\_k} + \frac{\lambda\_k \boldsymbol{a}\_k}{\eta\_k} (\boldsymbol{\sigma}\_k \cdot \boldsymbol{\sigma}\_k) + \frac{\eta\_k}{\lambda\_k} (f(\boldsymbol{\sigma}\_k) - 1)\mathbf{I} = 2\eta\_k \mathbf{D} \tag{4}$$

where *f*(σ*k*) = <sup>2</sup> *εk* exp<sup>h</sup> 2(Λ*k*−1) *qk* i<sup>1</sup> <sup>−</sup> <sup>1</sup> Λ*k* + <sup>1</sup> Λ*k* 2 1 − *αkλ<sup>k</sup>* 2 *tr*(σ*<sup>k</sup>* ·σ*<sup>k</sup>* ) 3*η<sup>k</sup>* 2 , Λ*<sup>k</sup>* = q 1 + *λk tr*σ*<sup>k</sup>* 3*η<sup>k</sup>* is the backbone stretch, *ε<sup>k</sup>* = *λs*0*<sup>k</sup> λk* , *k* = 1, . . . , *n*; *k* is the current mode number, <sup>∇</sup> σ = *<sup>d</sup>*<sup>σ</sup> *dt* − σ · <sup>∇</sup>**V***<sup>T</sup>* − ∇**<sup>V</sup>** · <sup>σ</sup> is the upper convective derivative of the tensor <sup>σ</sup>; **<sup>I</sup>** is the unit tensor; p is the pressure; **D** = ∇**V** + (∇**V**) *T* /2 is the rate of deformation tensor; **V** is the velocity; *λ<sup>k</sup>* is the relaxation time; *η<sup>k</sup>* is the polymeric viscosity; *q<sup>k</sup>* is the number of arms at the backbone extremity of the Pom-Pom molecule; *α<sup>k</sup>* and *ε<sup>k</sup>* are the rheological parameters.

#### *2.2. Parametric Method*

In the present paper, we suggest that a solution method based on the search for an unknown functional dependency *Vz*(*r*) should be used in a parametric form. Let us consider a solution method for multimode Giesekus fluid flows in round pipes. For this purpose, we rewrite (3) in a cylindrical coordinate system according to [16]. After some transformations and introduction of the parameter

$$
\rho\_k = \frac{\lambda\_k}{\eta\_k} \left| \sigma\_{rr(k)} \right| \tag{5}
$$

we obtain the following expression for nonzero components of the elastic stress tensor and the shear rate gradient:

$$
\sigma\_{rz(k)} = -\frac{\eta\_k}{\lambda\_k} A\_{k\prime} \,\sigma\_{zz(k)} = \frac{\eta\_k \rho\_k}{\mathfrak{a}\_k \lambda\_k} \left( \frac{2 - \mathfrak{a}\_k (\rho\_k + 1)}{(1 - \rho\_k)} \right) \,\,\sigma\_{rr(k)} = -\frac{\eta\_k}{\lambda\_k} \rho\_k \tag{6}
$$

$$\frac{dV\_z}{dr} = -\frac{1 + \rho\_k (1 - 2a\_k)}{\left(1 - \rho\_k\right)^2} \frac{A\_k}{\lambda\_k} (k = 1, \dots, n) \tag{7}$$

Hereinafter *A<sup>k</sup>* = √ *ρk* (1−*αkρ<sup>k</sup>* ) <sup>√</sup>*α<sup>k</sup>* . Introducing the parameter *ρ<sup>k</sup>* into (6), we get

$$\sigma\_{rz} = -\frac{\eta\_N (1 - 2\alpha\_k \rho\_k + \rho\_k)}{\lambda\_k (1 - \rho\_k)^2} A\_k - \sum\_{j=1}^n \frac{\eta\_j}{\lambda\_j} A\_j \ (k = 1, \dots, n) \tag{8}$$

$$\widetilde{r} = \frac{2\mathcal{R}}{\left|\mathbb{C}\_{0}\right|} \left( \frac{\eta\_{N}\left(1 - 2a\_{k}\rho\_{k} + \rho\_{k}\right)}{\lambda\_{k}\left(1 - \rho\_{k}\right)^{2}} A\_{k} + \sum\_{j=1}^{n} \frac{\eta\_{j}}{\lambda\_{j}} A\_{j} \right), \left(\widetilde{r} = r/R\right) \tag{9}$$

The remaining components of the stress tensor can be obtained similarly: *σzz*, *σrr*.

The parameter *ρ<sup>k</sup>* was introduced for each of the modes; therefore, to obtain the parametric dependences *dVz*(*r*)/*dr*, *Vz*(*r*), and *σij*(*r*) we select the parameter *ρ*1, corresponding to the first mode. The remaining parameters *ρk*(*k* 6= 1) must be expressed in terms of the main one *ρ*<sup>1</sup> using equations (7).

The parametric dependence of the dimensionless shear rate gradient will have the form:

$$\begin{cases} \begin{array}{c} \frac{d\boldsymbol{v}\_{z}(\rho\_{1})}{d\tilde{r}} = -\frac{1+\rho\_{1}(1-2\boldsymbol{a}\_{1})}{\left(1-\rho\_{1}\right)^{2}} \frac{A\_{1}}{\tilde{\lambda}\_{1}\mathcal{W}i} \\\ \tilde{r}(\rho\_{1}) = \frac{1}{\mathcal{K}r} \left( \frac{\tilde{\eta}\_{\mathcal{N}}(1-2\boldsymbol{a}\_{1}\rho\_{1}+\rho\_{1})A\_{1}}{\tilde{\lambda}\_{1}\left(1-\rho\_{1}\right)^{2}} + \sum\_{j=1}^{n} \frac{\tilde{\eta}\_{j}\sqrt{\rho\_{j}\left(\rho\_{1}\right)\left(1-\boldsymbol{a}\_{j}\rho\_{j}\left(\rho\_{1}\right)\right)}}{\sqrt{\pi}\_{j}\tilde{\lambda}\_{j}} \right) \end{array} \tag{10}$$

where *ρj*(*ρ*1) is determined with any degree of accuracy from (7) using numerical methods, for example, by dichotomy method; <sup>e</sup>*λ<sup>j</sup>* <sup>=</sup> *<sup>λ</sup>j*/*λa*, *<sup>η</sup>*e*<sup>j</sup>* <sup>=</sup> *<sup>η</sup>j*/*η*<sup>0</sup> and *<sup>η</sup>*e*<sup>N</sup>* <sup>=</sup> *<sup>η</sup>N*/*η*<sup>0</sup> are dimensionless simplexes; *λ<sup>a</sup>* = *n* ∑ *i*=1 *λiηi*/ *n* ∑ *i*=1 *ηi* is the average relaxation time; *η*<sup>0</sup> = *η<sup>N</sup>* + *n* ∑ *i*=1 *ηi* is the maximum possible fluid viscosity; *v<sup>z</sup>* = *Vz*/*V<sup>a</sup>* is the dimensionless velocity; *V<sup>a</sup>* = *Q*/ *πR* 2 is the average bulk velocity; *Q* is the fluid flow rate through the cross-section of the circular pipe; *R* is the radius of the circular pipe. Here *Wi* = (*λaR*/*Va*) is the Weissenberg number, *Kr* = *Wi* · *Eu* · Re<sup>0</sup> · Γ is the dimensionless complex; *Eu* = |*C*0*l*|/ *ρV<sup>a</sup>* 2 is the Euler number; Re<sup>0</sup> = (*ρVa*(2*R*))/*η*<sup>0</sup> is the Reynolds number; *l* is the pipe length; ∆*P* = |*C*0*l*| is the pressure drop in a round pipe across the length *l*.

To obtain the dependence of the dimensionless velocity *v<sup>z</sup>* on the dimensionless variable in a parametric form, it is necessary to use the condition

$$\int\_{0}^{1} v\_{z}\widetilde{r}d\widetilde{r} = \frac{1}{2}\tag{11}$$

The parametric dependence of the axial component of the dimensionless velocity on the dimensionless coordinate can be obtained from (10) by integrating these relations

$$\begin{cases} \begin{aligned} v\_{z}(\rho\_{1}) &= \frac{\tilde{\eta}\_{\text{N}}}{2Kr\tilde{W}\tilde{\lambda}\_{1}^{-2}\alpha\_{1}} \frac{\left(1+(1-2a\_{1})\rho\_{1}^{w}\right)^{2}\rho\_{1}^{w}\left(1-a\_{1}\rho\_{1}^{w}\right)}{\left(1-\rho\_{1}^{w}\right)^{4}} \\ &- \frac{\tilde{\eta}\_{\text{N}}}{2Kr\tilde{W}\tilde{\lambda}\_{1}^{-2}\alpha\_{1}} \frac{\left(1+\left(1-2a\_{1}\right)\rho\_{1}\right)^{2}\rho\_{1}\left(1-a\_{1}\rho\_{1}\right)}{\left(1-\rho\_{1}\right)^{4}} + \frac{1}{2Kr\tilde{W}l} \sum\_{j=1}^{n} \frac{\tilde{\eta}\_{j}}{\tilde{\lambda}\_{j}^{-2}a\_{j}} B\left(\rho\_{j}\right) \\ \tilde{r}(\rho\_{1}) &= \frac{1}{Kr} \left(\frac{\tilde{\eta}\_{\text{N}}\left(1-2a\_{1}\rho\_{1}+\rho\_{1}\right)A\_{1}}{\lambda\_{1}\left(1-\rho\_{1}\right)^{2}} + \sum\_{j=1}^{n} \frac{\tilde{\eta}\_{j}\sqrt{\rho\_{j}\left(\rho\_{1}\right)\left(1-\kappa\_{j}\rho\_{j}\left(\rho\_{1}\right)\right)}}{\sqrt{\pi\_{j}\lambda\_{j}}}\right) \end{aligned} \tag{12}$$

where *B ρj* = *ρ w* R1 *ρj* (1+(1−2*αj*)*ρj*)(1−2*αjρj*) (1−*ρj*) <sup>2</sup> *dρ<sup>j</sup>* is calculated analytically; *ρ w* 1 is the parameter

value on the pipe wall.

The direct use of analytical expressions (12) is hampered by the fact that there is a functional dependence between the dimensionless complexes *Wi* and *Kr*, which can be determined using (11).

Similarly, a parametric problem solution for multimode fluid flows can be obtained using eXt. Pom-Pom can be obtained. In this case, it is convenient to use the following expression as a parameter (*ρ<sup>k</sup>* )

$$
\sigma\_{\varphi\varphi(k)} = \frac{\widetilde{\eta}\_k}{\widetilde{\lambda}\_k} \frac{\eta\_0}{\lambda\_d r^2} (-\rho\_k)\_\prime \,\,\widetilde{\eta}\_k = \eta\_k / \eta\_{0\prime} \,\, \widetilde{\lambda}\_k = \lambda\_k / \lambda\_d . \, (0 \le \rho\_k \le 1) \tag{13}
$$

Then

$$\begin{cases} \begin{aligned} v\_z(\rho\_1) &= \int\_1^y \dot{\gamma}\_k(\rho\_1) (d\tilde{r}/d\rho\_1) / \tilde{\lambda}\_k d\rho\_1 \\ \tilde{r}(\rho\_1) &= -\left(\frac{\tilde{\eta}\_N}{\tilde{\lambda}\_m} \dot{\gamma}\_m(\rho\_m) + \sum\_{k=1}^n \frac{\tilde{\eta}\_k}{\tilde{\lambda}\_k} \frac{\lambda\_k}{\eta\_k} \sigma\_{rz(k)}(\rho\_k)\right) (Kr)^{-1} \end{aligned} \end{cases} \tag{14}$$

where dependences *ρk*(*ρ*1) can be obtained from the definition of the shear rate gradient

$$\dot{\gamma} = \frac{dv\_z}{d\tilde{r}} = \frac{\dot{\gamma}\_1}{\dot{W}i\tilde{\lambda}\_1} = \dots = \frac{\dot{\gamma}\_n}{\dot{W}i\tilde{\lambda}\_n} \tag{15}$$

A more detailed description of the method is provided in [17].
