Numerical simulation of a stabilizing Poiseuille-type polymer fluid flow in the channel with elliptical cross-section

This work is devoted to the numerical analysis of stabilization of the incompressible viscoelastic polymer fluid flow in the channel with elliptical cross-section. To describe the flow, mesoscopic rheological relations are used, and resolving non-stationary equations are derived. For solving them a special pseudo-spectral method is developed and implemented. As time increases, under certain conditions on the parameters of flow the solution to the non-stationary problem stabilizes and converges to the one of three branches of the solution to the corresponding stationary problem. It is shown that the variation of the parameters describing polymer microstructure leads to the switch of stabilized solution between these branches. The work provides the results of simulation of the flow stabilization and the analysis of the threshold values of parameters at which the switching occurs.


Introduction
Polymer fluid flows demonstrate complex dependencies of the viscosity on the strain rates, the effects of relaxation and anisotropic behavior, which is due to the complex structure of macromolecular coils of polymer. Even the most basic types of the flow, like the Poiseuille and Couette ones, discover new effects that are obviously not inherent to the classical Newtonian fluids.
To take into account such effects in our study, the mesoscopic rheological relations of the Pokrovskii-Vinogradov model are used, see [1]- [4]. After transforming them, the equations describing non-stationary flows in the channel with elliptical cross-section are derived under the assumption that a constant pressure drop acts along the channel. The goal of the work is to study the process of stabilization of these flows and to compare the stabilized solutions of the non-stationary equations with the solutions of stationary formulation, which were recently obtained in [3].
In order to simulate the flows, a special pseudo-spectral method based on the Fourier and Chebyshev approximations was developed. The method uses the ideas of the algorithms without saturation from [5] and allows one to perform fast and highly-accurate computations. The latter turns out to be crucial for the careful comparison of the solutions to the stationary and nonstationary formulations. Firstly, it allowed us to study the sensitivity of the stabilization process to the parameters of the model. Secondly, the accurate computations were the key ingredient in discovering and description of the process of switching the stabilized solution to the nonstationary problem from one branch of the solution to the stationary problem to another.

Formulation of the problem 2.1. Governing equations
Following the works [1], [2], [7]- [9] and the paper [6] we can write the equations of the rheological mesoscopic Pokrovskii-Vinogradov model describing the polymer fluid flows. These equations in the dimensionless form in Cartesian coordinate system (x 1 , x 2 , x 3 ) take the following form divu = 0, Here t is the time; u 1 , u 2 , u 3 are the components of the velocity vector u.
is the Reynolds number ; ρ(= const) is the density of medium; η * 0 , τ * 0 are the initial values of the share viscosity and of the relaxation time as T = T 0 , where T 0 is the temperature of environment; l and u H are the characteristic length and velocity, respectively.
Further we shall use the notations x, y, z for the spatial coordinates (x = x 1 , y = x 2 , z = x 3 ).
Remark 1 In the works [6,10] the cases k = β, k = 1.2β have been considered. In this work the more general case is studied: k = c k β, where c k > 0 is new parameter of the microstructure.
To pose the boundary value problem in the domain Ω, we pass to the elliptic coordinates ζ, γ and make the corresponding change of variables in (5)- (7). We set the distance between foci of the ellipse Ω equal to 2δ, the length of the large semi-axis equal to 1, and the length of the small semi-axis equal to r < 1. Therefore, the change of the variables can be written as In the order to obtain the resolving system of equations in the elliptic coordinate system, the derivatives with respect to the variables y, z in (5)-(7) were replaced by the derivatives with respect to ζ, γ.
Let us supplement the resolving equations with the no-slip condition on the boundary of Ω and with the free boundary condition on the line between foci where ζ = Z r is the equation of the boundary of ellipse. Assuming that the fluid is at the rest state in zero moment of time we set the initial data for the resolving equations as follows:

Stationary solutions
It is worth noting that the setup (4) in the stationary case allows for the considerable reduction of the resolving equations (see the details in [6], [3]). If the unknowns of the problem do not depend on the time variable: u = u(ζ, γ), a ij = a ij (ζ, γ), i, j = 1, 3, P = P(ζ, γ), A = const, then system (1)-(3) reduces to the one resolving equation for the component u of the velocity vector.
However, as mentioned in [3], in order to derive this equation, one needs to solve the nonlinear cubic equation relating the diagonal and non-diagonal components of the anisotropy tensor or, more precisely, the quantities σ = a 22 + a 33 and µ = a 2 12 + a 2 13 . The analysis of this cubic equation showed that for small enough value of the pressure drop it has three real solutions that can be expressed using Cardano's formula. Let us denote these solutions by σ 1 (µ), σ 2 (µ), σ 3 (µ) and the corresponding branches of the solution of stationary problem -by 1,2,3.
A natural question arises: which of the three possible stationary flows is implemented in real life? Non-stationary formulation derived above is aimed at studying this issue. To do this, we shall assume that the pressure drop in this formulation is constant and run the simulations till the solution become stable at some moment of time t = t * . Comparison of the solution u(t * , y, z) of non-stationary problem with three different branches of the solution of stationary problem corresponding to σ 1,2,3 (µ) will show, which stationary flow implements in reality.

Approximation of the unknown functions and their derivatives
We shall search for the solutions of problem (5)-(7), (9)-(11) written in elliptic coordinates in the class of sufficiently smooth functions that are 2π-periodic with respect to the coordinate γ. To design the computational algorithm, it is necessary to approximate the values of functions u(t, ζ, γ) and α ij (t, ζ, γ) at the time moments t = t n by elements of finite-dimensional function spaces. The derivatives of u and α ij with respect to spatial variables ζ and γ must be approximated too. To this end, we use the spaces of algebraic and trigonometric polynomials. To make the projection onto such spaces, we use the modified Lagrange polynomials with Chebyshev nodes and the polynomials with Dirichlet kernel. More precisely, the functions u(t n , ζ, γ) and α ij (t n , ζ, γ) are interpolated by tensor products of such polynomials. Some important properties of such approximations including the absence of saturation, the slow growth of the Lebesgue constants, and the error estimates are discussed in [12], chapter 3.
The domain Ω can be represented as follows: where are the cardinal functions of the interpolation, T K (ζ) is the Chebyshev polynomial of the power K; D M is the Dirichlet kernel. The modified cardinal function l b kK (ζ) in (14) accounts for the boundary conditions (9), (10).
For approximating equations (5)- (7) in Ω we shall use the collocation method with the nodes (ζ k , γ m ). To this end, we differentiate expressions (13), (14) with respect to ζ and γ and pass to the limit (ζ, γ) → (ζ l , γ q ), l = 0, K, q = 1, M. Application of the l'Hopital's rule and the basic properties of the cardinal functions allows one to represent the first derivatives of polynomial (13) at each point (ζ l , γ q ) as a linear combination of the values of (α n ij ) km -the first derivative with respect to ζ is the linear combination with variable index k = 0, K and fixed m = q; the first derivative with respect to γ is the linear combination with variable index m = 1, M and fixed k = l. The same is valid for the first and the second derivatives of (14).
We collect the coefficients of these linear combinations into matrices and denote by A 1 , B 1 the (K + 1) × (K + 1)-and M × M-matrices of the coefficients for the first order derivatives of α n ij (ζ, γ) with respect to ζ and γ correspondingly. By analogy, A 1 , B 1 are the matrices obtained for the first order derivatives of u n (ζ, γ); and A 2 , B 2 are the matrices obtained for the second order derivatives. For more details on this approximation of the derivatives, which is usually associated with pseudo-spectral methods, we refer the reader to the works [5,11]. Let us denote by Υ n ij , (Υ n ij ) ζ -K × M-matrices with the elements α n ij (ζ k , γ m ) and ∂α n ij ∂ζ (ζ k , γ m ), K = K + 1; by U n , U n ζ , U n ζζ -K × M-matrices with the elements u n (ζ k , γ m ), ∂u n ∂ζ (ζ k , γ m ), ∂ 2 u n ∂ζ 2 (ζ k , γ m ).
For approximation of the derivatives with respect to the variable ζ one obtains the following matrix formulas: For approximation of the derivatives with respect to the variable γ one has ζγ , U n ζγ are determined by analogy with (Υ n ij ) ζ , U n ζ . For constructing the algorithm we apply the spectral decomposition of the matrices A 2 , B 2 approximating the second derivatives where R A , R B are the matrices of eigenvectors of A 2 and B 2 ; and D A , D B are the diagonal matrices of eigenvalues d k A , d m B of the matrices A 2 and B 2 . To construct decomposition (17), we used the fact that the matrix B 2 is symmetric, which follows from the explicit expressions for its elements. Therefore, R −1 B = R T B . By using the interval analysis, it was shown that the eigenvalues of the matrix A 2 are all distinct and real numbers. It justifies decomposition (17) for A 2 .

Problem of linear algebra
After linearization of the equations (5), (6) written in elliptic coordinates with respect to the unknown functions α ij we used (12) together with collocation method and approximations (15), (16) over the spatial variables. Dropping-off the error terms we passed to the linear matrix equations that express the matrices Υ n+1 ij through Υ n ij and U n . Let us write for instance one of these equations, namely, the one for Υ n+1 11 , which follows from the first equation of system (6): Here, the dot denotes the element-wise product of matrices, E I is the K×M-matrix filled up with unit elements, R, Q are the K × M-matrices containing the values of R(ζ, γ) = sinh ζ sin γ δ(sinh 2 ζ+cos 2 γ) , Q(ζ, γ) = cosh ζ cos γ δ(sinh 2 ζ+cos 2 γ) written in the nodes ζ k , γ m , respectively (they appear after passage to elliptic coordinates). The initial condition for (18) is posed in accordance with (11): Υ 0 11 = (0) K×M is zeroes matrix. Approximation of other five equations for α ij can be written by analogy.
The approximation of (7) is as follows: we passed to the following matrix equation: It is easy to see that the elements v n+1 km of the unknown matrix V n+1 express as follows: where g n km are the elements of the matrix G n . Here, the time step τ must satisfy the condition , ∀k, m. Taking into account that ae 2 = (W Re) −1 , we can write , ∀k, m.
After obtaining the elements of the matrix V n+1 , the solution on the (n + 1)st time step can be computed as U n+1 = R A V n+1 R B . Since the problem is non-linear, the additional iterations on each time step have been applied to resolve non-linearity. It should be noticed that in comparison to the classical scheme of the collocation method the proposed algorithm essentially reduces the number of operations. Indeed, for computing of Υ n+1 ij in accordance with (18) it is necessary to do KM operations (they are componentwise divisions). For computing U n+1 the most consuming operations are the multiplications of matrices:

Results of simulations and discussion
The designed algorithm was used for solving non-stationary problem (5)-(7), (9)- (11). The pressure drop in all simulations was constant: A(t) ≡ A 0 . The non-linear stationary problem strictly corresponding to the considered non-stationary one was solved using the relaxation method (for details see [3] and the algorithm from [5]). It should be noted that the relaxation method in the stationary problem converges only when certain (not necessary all) branches σ 1,2,3 (µ) are used in the resolving equation for stationary velocity field (see section 2.3). The indexes of branches, for which the relaxation process converges, as well as, the index of soughtfor branch yielding the stabilized solution of non-stationary problem (5)-(7), (9)-(11) depend upon the parameters of the model.
In simulations we set ε S = 10 −7 and checked that t * practically does not depend on the number of collocation nodes K, M and on τ . For the results presented below we set K, M in the diapason 21-41, τ ∈ [0.01, 0.25]. Let us denote by D u (ζ, γ) the relative deviation of the non-stationary velocity at the moment t * from the solution of the stationary equation for the velocity. Simulations show that in the case c k = 1 and β < 0.5 (this values have been introduced in [1,10]) the relaxation method used for searching the solution of the stationary problem converges only when using two branches of the function σ(µ): σ 1,3 (the expressions for them are given in [3]). Moreover, the stabilized solution of the non-stationary problem coincides with the solution obtained using the branch σ = σ 3 . The values of t * and ∥D u (ζ, γ)∥ obtained in this case for various values of the parameters W , r and c k are shown in table 1. It can be seen that the dependence of the stabilization time t * on the value of W is strongly non-linear. The time t * is also very sensitive to the variation of the channel's width r. While decreasing r, the flow stabilizes much faster. This also remains valid for much larger values of W and A 0 than those given in the table.
An interesting picture can be observed when varying c k : the deviations ∥D u (ζ, γ)∥ remain small for practically all values of c k except for the small vicinity of a certain point (in the case shown in table 1 this point is c k ≈ 1.25). In this vicinity the deviations increase very fast. Such behavior is related to switching the stabilized solution of the non-stationary problem from the branch 3 of the solution of stationary problem to the branch 1. It means that the solution in this vicinity has strong singularity. Table 1. Values of t * and ∥D u (ζ, γ)∥ obtained in tests with Re=10 and variation of the parameters W , r, c k

Variation of W for
Variation of r for Variation of c k for Simulations showed that the convergence of the relaxation process in stationary problem takes place not for all possible branches of σ(µ). The convergence for each branch can appear and disappear when changing the parameters of the problem. The most essential influence is exerted by β and c k . Moreover, a certain relation between these parameters provides the switching of stabilized solution of the non-stationary problem from one branch of the solution of the stationary problem to another.
• c k < 1.232 -there are two stable solutions of the stationary problem (branches σ 1 , σ 3 ); the stabilization goes to the solution of σ 3 . • c k ≈ 1.232 -stationary solution of the branch σ 3 disappears; stabilization switches to the solution corresponding to the branch σ 1 . • 1.232 < c k < 1.377 -there is only one solution of the stationary problem corresponding to the branch σ 1 , and the non-stationary solution converges to it. • c k ≈ 1.377 -three branches of the stable solution appear (branches σ 1 , σ 2 , σ 3 ).
• c k > 1.377 -three branches of the stationary solution are preserved; stabilization goes to the solution of σ 1 .
In this test the threshold values of the parameter c k have been obtained in numerical simulations and they are approximate. The described scenario conserves when changing the parameters of the model (of course the threshold values of c k also change).
Let us denote by c 3→1 k the value of c k , at which the stabilized solution switches from the branch σ 3 to the branch σ 1 , figure 2 shows the dependence of c 3→1 k on β for different values of W . It can be seen that for small W the curves c 3→1 k (β) (the gray and the dashed lines) practically coincide. However, while W increases the curve changes and goes below the previous ones. The same picture can be observed while increasing A 0 .
Let us denote by c 1,2,3 k the minimal value of the parameter c k that provides the convergence of the relaxation method for all three branches of function σ(µ) in the stationary problem.