NUMERICAL PREDICTIONS OF STEADY AND UNSTEADY FLOWS OF FENE-CR VISCOELASTIC FLUIDS WITH A NEWTONIAN SOLVENT

In this work, we presente numerical results of the flows of FENE-CR viscoelastic fluids with a Newtonian solvent in a channel and in a cross-slot geometry. The governing equations for transient incompressible flows using the FENE-CR fluid are solved by the finite difference method on a staggered grid. The fluid is modelled by a Marker-and-Cell type method and an accurate representation of the fluid surface will be employed. The full free surface stress conditions are considered. Indeed, an implicit treatment to solve the momentum equation is employed. Predictions for fully developed pipe flow are compared with the analytic solution to validate the numerical technique. Other problem perform in this work is the low Reynolds number viscoelastic flows in a cross-slot geometry. In particular, Rocha et. al. [12] carried out a systematic study in a planar cross-slot geometry using a FENE-CR fluid. In order to demonstrate that the numerical method developed in this work can capture the asymmetry predicted by Rocha et. al., we present results for flows of FENE-CR fluids in the planar cross-slot geometry. About all simulations presented in this work, our results are in quantitative or qualitative agreement with numerical predictions or analytical solutions.


INTRODUCTION
Our research team have studied a wide range of problems involving viscoelastic free surface flows.The technique employed is well known by many researchers from this area, which is a finite difference method on a staggered grid and the fluid is modelled by a Markerand-Cell approach ( [7]).This metodology have been capable of simulating many problems involving unsteady and steady viscoelastic free surface flows, such as, jet buckling, dieswell, impacting drop and fountain flow (see, as example, [8,14,15,16,17]).
In this work, we are specially concerned with confined flows of viscoelastic fluids in a planar cross-slot geometry.This problem is very important to investigate purely elastic instabilities of viscoelastic flows because Newtonian fluids at the same condition remain prefectly steady.Several researchers have been studing this problem (or similar) in 2D and 3D dimensions ( [1,10,11,12]).Poole, Alves and Oliveira [11] successfully used a numerical technique with a simple viscoelastic constitutive equation (UCM) to model this steady asymmetry under creeping-flow conditions (Re = 0.0).In that work, they probed that the inertia decreases the degree of asymmetry and stabilizes the flow.The effects of inertia to Re = 0.01 do not influence the magnitude of the asymmetry, the datas to Re = 0.0 and Re = 0.01 are absolutly the same.However for Re = 1, 2, and 5 the influence of inertia is stronger.Arratia, Thomas, Diorio and Gollub [3] have reported two flow instabilities in a cross-slot geometry.In the first instability the flow becomes deformed and asymmetric, but remains steady.The second instability occurs at higher strain rates and the flow becomes unsteady and fluctuates non-periodically in time.Recently, Rocha, Poole, Alves and Oliveira [12] extended previus studies in infinitely extensionable models (UCM and Oldroyd-B fluids) to finite extendable non-linear elastic models (FENE-P and FENE-CR models).They provided quantitative data of benchmark quality for flows in two geometries, a planar cross-slot geometry and roundedcorner geometry, by exploring the effects of finite extensibility (L 2 ), Weissenberg number (W i) and concentration of polymer solution (β).
In the present work, we claim to get a good agreement between our results with Rocha et al. [12] ones.Moreover, predictions for fully developed flow of FENE-CR fluid with Newtonian solvent into a channel is also presented in this work.We compare the exact solution with numerical results to validate our metodology for FENE-CR fluids.In Section 2, 3 and 4, we present the governing equations, computational domain, boundary conditions and method of solution.The validation of approach by simulating viscoelastic free surface flows into a channel is presented in Section 5.The Section 6 is dedicated to presente and discuss numerical results of the flow of FENE-CR fluids in a planar cross-slot geometry.

GOVERNING EQUATIONS
Let consider incompressible and isothermal flows described by the mass conservation equation and the equation of motion where t is the time, u is the velocity vector, p is the pressure, ρ is the fluid density, g is the gravitational field and τ is the extra-stress tensor.For considering viscoelastic fluid with Newtonian solvent, the extra-stress tensor is splited by where 2η S D and T represent the Newtonian solvent contribution and the viscoelastic contribution tensors, respectively.D = 1 2 ∇u + (∇u) T is known as the rate of deformation tensor and η S represents the solvent viscosity.
The viscoelastic contribution tensor T is linked to configuration tensor A of the FENE-CR model ( [5]) by the Kramers' form as, where the dimensional number λ is the fluid relaxation time, η P represents the polymeric viscosity and the dimensionless function f (tr(A)) is given by The evolution equation for the configuration tensor of the FENE-CR model is given by where the symbol The operator D/Dt = ∂/∂t + u • ∇ and (∇u) T is the transpose of the velocity gradient.
Equations ( 1)-( 5) are solved for the dependent variables u, p, A and T by setting the following non-dimensionalization, where H, U and ρ denote typical length, average velocity and density scalings, respectively.The amount of Newtonian solvent contribution is controlled by the constant non-dimensional β = η S η 0 , with η 0 = η S + η P representing the total viscosity.Therefore, the system of equations become (the bars are omitted for clarity) ) where H are the Froude, the Reynolds and the Weissenberg numbers, respectively.

COMPUTATIONAL DOMAIN AND BOUNDARY CONDITIONS
The governing equations are solved by the finite difference method on a staggered grid with cell dimensions δx × δy.The velocities are calculated at the cell faces and all other quantities are computed at the cell center.A detailed description of this arrangement can be found in [7,18].A fundamental concept behind the MAC method is the classification of the grid cells according to their position relative to the fluid as, cells that do not contain fluid (Empty (E) ); cells full of fluid and which are not contiguous with any empty cells (Full (F)); cells that contain fluid and have at least one face contiguous with an empty cell face (Surface (S)); cells that define a rigid boundary (Boundary (B)); cells that define an inflow boundary (Inflow (I)) and cells that define an outflow boundary (Outflow (O)).
To solve equations ( 7) - (10) we need to specify appropriate boundary conditions for u.For the momentum equation, we employ the no-slip condition (u = 0) on rigid boundaries (B cells).At fluid entrances (I cells) the velocity is specified by u n = U and u m = 0, while at the fluid exits (O cells), we impose homogeneous Neumann conditions, namely ∂u n ∂n = 0, and ∂u m ∂n = 0, where the subscripts n and m denote the normal and tangential directions to the fluid inlet (and to the fluid exit), respectively.
The S cells contain the free surface.On these free boundaries, it is necessary to impose special boundary conditions for the velocity and pressure fields.We shall consider unsteady free-surface flows of a viscous fluid moving into a passive atmosphere (which we may take to be at zero pressure).In the absence of surface tension the normal and tangential components of the stress must be continuous across any free surface, so that on such a surface (see [4]) where the stress tensor is given by In equations (11), n = (n x , n y ) represents a unit vector normal and external to the surface, and m = (m x , m y ) is a unit vector tangent to the free surface (Fig. 1).By using two-dimensional Cartesian co-ordinates and taking m = (n y , −n x ), the stress conditions (11) can be written as Eqs. ( 12)-( 13) represent the appropriate boundary conditions on the fluid free surface in the absence of surface tension.

Computation of the configuration tensor A on mesh boundaries
When solving the convective terms of the FENE-CR constitutive equation one needs to employ an effective method for approximating the first order derivatives contained within the material derivative; otherwise the numerical solution can deteriorate and produce unphysical results.It is known that the use of a second order central difference scheme can lead to oscillatory solutions while first order upwind schemes tend to produce solutions that contain excessive numerical diffusion.In this work we employ a high order stabilized upwind scheme known as CUBISTA (Convergent Universally Bounded Interpolation Scheme for the Treatment of Advection) [2].This method requires the values of a generic variable, say φ, that can be positioned upstream (φ U ), downstream (φ D ) or remote-upstream (φ R ) with respect to a reference point (see Fig. 2) at which the variable is being approximated.
Reference points used for the CUBISTA upwind scheme.
It is clear then that, when calculating the convective terms of the components of the non-Newtonian stress tensor in zones that are near to the boundary domain, the values of the components of the stress on mesh boundaries are required.On O cells (outflows), it is considered ∂A ∂n = 0, on I cells (inflows), it is setted A = 0, on B cells (rigid boundaries), it is solved by equations ( 9) as follows.
Let consider the cases: rigid boundaries parallel to the x−axis and rigid boundaries parallel to the y−axis.On wall parallel to x-axis, as exemplo, from the no-slip condition we have ∂ ∂x = 0 ⇒ ∂v ∂y = 0. Thus, only the derivatives ∂u ∂y are not zero.In this case, equations (9) (in cartesian coordinate) reduce to Note that, the boundary conditions on the fluid free surface requires values to non-Newtonian contribuition T (see Eqs. ( 11)).Simply, the tensor T is updating by equations ( 10), except on outflows, which is employed ∂T ∂n = 0.

METHOD OF SOLUTION
It will be presented the cross-slot problem that requires a small Reynolds number to study purely-elastic instabilities ocurring for the viscoelastic flows in such geomety.Therefore, to avoid the parabolic stability restriction inherent in explicit schemes, the momentum equation ( 7) is integrated implicitly in time by the Crank-Nicolson method, following the ideas of Oishi et.al. [9].In this case, the Navier-Stokes equations ( 7) and ( 8) may be rewritten as Following the ideas of the projection method, first proposed by Chorin [6] e Temam [13], every smooth vector field can be decomposed as a sum of a gradient of a potential function and a divergence-free vector field, The gradient of the potential ψ is the minimum perturbation that will make the final velocity field satisfy the incompressibility constraint at the next time level.
Taking the divergence of Eq. ( 19) and imposing mass conservation on u, one obtains the following Poisson equation for ψ We suppose that at time t n the velocity field u(x, t n ) = u (n) and the configuration tensor A(x, t n ) = A (n) are known and the values of u and A on the boundary are given.The velocity field u(x, t n+1 ) = u (n+1) and the configuration tensor A(x, t n+1 ) = A (n+1) , where t n+1 = t n + δt, are calculated by the following steps: Step 1: The tensor A is discretized using the second-order accuracy Runge-Kutta (RK21) method.To describe it, the Eq. ( 9) is re-written as where From Eqs. ( 21) and ( 22), calculate an intermediate configuration tensor A (n+1) by explicit forward Euler discretization, Thus, update the non-Newtonian tensor T (n+1) by approximating the Eq. ( 10), Step 2: Calculate a provisional velocity field u (n+1) using the Eq. ( 17), where the boundary conditions for u (n+1) are the same as those for u (n+1) and p (n) is an approximation to p (n+1) .
Step 3: Solve the Poisson equation (20) for the potential function ψ, The boundary conditions for Poisson's equation are ψ (n+1) = 0 on outflows while the homogeneous Neumann condition is used for fixed boundaries and inflows.Details about the calculation of ψ at free surfaces can be found in [9].
Step 4: Compute the final velocity from decoupling equation ( 19) Step 5: Compute the pressure field from Note that, to obtain the final pressure field, we substitute equation ( 27) into equation ( 25) and comparing it with equation ( 17) to obtain the correction equation (28).
Finaly, update the non-Newtonian tensor T (n+1) from Eq. (10), Step 7: Update the markers positions: The last step in the calculation is to move the markers to their new positions.This is performed by solving for each particle.In this work, we follow the ideas of Oishi et.al. [8] using the second-order RK21 scheme to solve equation (31): where the velocities are calculated at the required positions using a bilinear interpolation from the nearest values in the mesh.The x is an intermediate position for the more accurate calculus using equation (33).Details about the calculation of the pressure at free surfaces are presented in Oishi et.al. [9,8]

VALIDATION OF THE APPROACH
To validate the methodology presented in the Section 4 we considered two-dimensional Cartesian flows in a channel formed by two parallel plates as shown in Fig. 3.
At the channel entrance we imposed the analytic values of the velocity field u(y) as well the analytic values of the components of the configuration tensor, A xx , A xy and A yy .Recall that A and T are related through the Kramers' form as Eq.(10).At the channel walls the velocity satisfied the no-slip condition and the configuration tensor A was calculated as already mencioned in this work.Indeed, at the channel exit the velocity and the configuration tensor A was computed by the homogeneous Neumann condition.
The simulation commenced with the channel empty.The fluid was then injected at the channel entrance so that the channel was gradually filled.Initially, there was a free surface within the channel.On the free surface, the boundary conditions were the free surface stress conditions presented in this work.The solution was calculated using the numerical method presented in Section 4.
At time t = 0s the channel was empty and the fluid had been injected in at the channel entrance until the channel is completely full and steady state has been achieved.Figure 4 shows the contours of the velocity and the components of the non-Newtonian tensor T at the time t = 3.5s while Fig. 5 displays the contours at the later time t = 50s.Both computations were performed on mesh M50.From Fig. 4, we point out that at time t = 3.5s the channel is still being filled and the contours are not parallel which points to the transient state of the flow.
However, at the time t = 50s, Fig. 5 shows that the channel is completely filled and the contour lines are all parallel indicating that steady state has essentially been reached.The solution within the channel should therefore be the analytic solution which was imposed at the channel entrance.To verify this fact, we calculated the values of the variables u, A xx , A xy , A yy , T xx and T xy at the middle of the channel (x = 5) using the results obtained on the five  meshes and compared these with the analytic solution.These values are shown in Fig. 6 for three meshes (M 10, M 30, M 50) for clarity the graphics.Indeed, we can observe in Fig. 6 that the numerical solutions obtained using the three meshes are in good agreement with the analytic solution.To demonstrate the convergence of the numerical methodoly we computed the Euclidean norms of the errors between the exact solution (ExSol) and the numerical solutions (NumSol) for each simulation by Eq. (34); the results are displayed in Table 1.  1 shows that as the mesh is refined the numerical solution converges to the analytic solution.

NUMERICAL RESULTS OF THE FLOW OF FENE-CR FLUIDS IN A PLANAR CROSS-SLOT GEOMETRY
In this Section, we will show that methodology presented in this work is capable of simulating the flow of FENE-CR fluids with a Newtonian solvent in a planar cross-slot geometry in comparison with the numerical results obtained by Rocha et al. [12].
The cross-slot problem requires a small Reynolds number to study purely-elastic instabilities ocurring for the flows in such geomety.We point out that in this Section we consider the Reynolds number of Re = 0.01 instead of zero as it was considered in [12].However, Poole, Alves and Oliveira probed in [11] that the effects of inertia to Re = 0.01 do not influence the magnitude of the asymmetry studied in that work.The datas to Re = 0.0 (creepingflow) and Re = 0.01 (low-Re) are absolutly the same but for Re = 1, 2, and 5 the influence of inertia is stronger (see [11], Fig. 3(a), pg.3).
As Rocha et al., we consider the domain, displayed in the Fig. (7), completely full of fluid.There are two vertical inflows, one on the left and other on the right hand.In the first, the fluid is injected from West to East and on the other one, the fluid flows in the opposite way.The outflows are put on horizontal place, on top and bottom.On the top outflow, the fluid leaves domain from South to North, and on the bottom outflow, the fluid goes away to South direction.
We set the analytical solution of the fully developed flow on the inflows for velocity and stress, and on the outflows and rigid boundaries we set the conditions already presented in this paper.The inflow length is H and the channel width is 5H.The central square of the cross-slot is H × H.In the paper [12], the authors considered 10H for the top and bottom channel lengths.They have confirmed to be sufficiently long to eliminate end effects.However, they imposed just the average value U on the inflow.For the simulations presented in this Section, it was not noticed any influence about the difference between these lengths.By our numerical method, we have to perform a wide range of this problem to confirm whether 5H is enough to eliminate end effects.The mesh adopted by Rocha al. is generated by five blocks: two horizontal channels with 50 × 51 cells (blocks I and V); two vertical channels with 51 × 50 cells (blocks III and IV); and the central square of the cross-slot with 51 × 51 cells (Block II).In our method, we can not do a local refinement, then we have to consider the same δx and δy along the domain.This fact justifies our choice about the length 5H.
The meshes, used in this Section, are provided in Table 2.The length H is setted by H = 1, so that the blocks are characterized as: Block I: −5.5 ≤ x ≤ −0.5, −0.5 ≤ y ≤ 0.5, Block II: −0.5 ≤ x, y ≤ 0.5, Block III: −0.5 ≤ x ≤ 0.5, −5.5 ≤ y ≤ −0.5, Block IV: −0.5 ≤ x ≤ 0.5, 0.5 ≤ y ≤ 5.5 and Block V: 0.5 ≤ x ≤ 5.5, −0.5 ≤ y ≤ 0.5.First of all, we simulated a flow of Newtonian fluid into a cross-slot geometry and we compared these results (by Mesh 2) with the numerical results from Rocha et al. (Fig. 10 from the paper [12], pg.66) to validate our technique.Figure 8 shows the profile of strain rate • ε= ∂u/∂x along the central line y = 0.As it is shown in Fig. 8, there is a good agreement between our numerical solution and the numerical solution from [12]. Figure 8 reveals that our numerical technique is capable of simulating flows of Newtonian fluid into a cross-slot geometry.To verify the convergence of the velocity u for the flows of FENE-CR fluid into this kind of geometry, we used both Mesh 1 and Mesh 2. Figure 9 a) shows the velocity profile along the central line (y = 0) for a Weissenberg number of W i = 0.2, L 2 = 200.0 and β = 0.1.Figure 9 a) reveals that the numerical results present a good level of agreement between the two meshes, even Mesh 1 to be very coarse for this problem.Indeed, we obtained the numerical strain rate results for flows of FENE-CR fluid to compare them   with the same results on the refined mesh to verify the convergence of the present numerical technique for viscoelastic FENE-CR fluids.Figure 9 b) also displays the Newtonian profile just for emphasizing the difference between the • ε= ∂u/∂x on the Newtonian and viscoelastic model.From Figure 9 b), the strain rate of FENE-CR flow decreases even after to reach the central square (−0.5 ≤ x, y ≤ 0.5).However, near the central point ((0, 0)) this rate increases slightly.
In all these simulations we have considered Re = 0.01, L 2 = 200.0 and β = 0.1.The Weissenberg number is a dimensionless parameter used to measure the effect of elasticity of the fluid on the flow.As predicted by Rocha et. al. ([12], Table 2., thirdy collumn -pg 61), the flow does not present asymmetry for Newtonian flow and for W i = 0.3; for W i = 0.45 the flow starts to get asymmetryc and it is totally asymmetryc for W i = 0.7.Look at the contours plots of velocity u and the normal tensor T yy at the Figures 10 and 11, respectively.Our results are in qualitative agreement with the results in [12].

CONCLUSIONS
This paper presente numerical results of free surface flows in a 2D-channel and confined flows in a planar cross-slot geometry.Both of problems, FENE-CR viscoelastic fluids with a Newtonian solvent are considered.The numerical solutions obtained from the free surface flows in a 2D-channel are in good agreement with the analytical solution and the convergence of the numerical solution to analytical solution can be noticed with a refinement of meshes.Other problem performed in this work is the low Reynolds number viscoelastic flows in a cross-slot geometry.This problem is very important to investigate purely elastic instabilities of viscoelastic fluids.In particular, Rocha et al. [12] carried out a systematic study in a planar cross-slot geometry using a FENE-CR fluid.Our results are in quantitative and qualitative agreement with numerical predictions from [12] so that this work demonstrate that the numerical method developed is capable of capturing the asymmetry studied by Rocha et al..

Figure 1 .
Figure 1.Flow in a channel: unit normal (n) and tangential (m) vectors to the surface of fluid.

Figure 3 .
Figure 3. Definition of the computational domain for 2D channel flow.

Figure 4 .
Figure 4. Numerical results of profiles of the components: (a) u, (b) T xx , (c) T xy , (d) T yy , at time t = 3.5s on the mesh M 50.

Figure 5 .
Figure 5. Numerical results of profiles of the components: (a) u, (b) T xx , (c) T xy , (d) T yy at time t = 50s on the mesh M 50.

Figure 6 .
Figure 6.Numerical and analytic solutions of components: (a) u, (b) A xx , (c) A xy , (d) A yy , (e) T xx and (d) T xy , at time t = 50s at position x = 5.The solid line represents the exact solutions while the numerical solution on the three meshes are represented by symbols.

Figure 8 .
Figure 8. Profile of strain rate for Newtonian flow along the central line y = 0. Comparison between our numerical solution and the numerical solution from [12] a) b)

Figure 9 .
Figure 9. a) Velocity profile |u|/U along the central line y = 0. b) Mesh refinement studies of the profile of strain rate for FENE-CR fluid along the central line y = 0.For both, W i = 0.2, L 2 = 200.0 and β = 0.1.

Table 2 .
Meshes: number of cells for each block.