Role of Nonmonotonic Constitutive Curves in Extrusion Instabilities

1College of Computer, National University of Defense Technology, Changsha, Hunan 410073, China 2State Key Laboratory of High Performance Computing, National University of Defense Technology, Changsha, Hunan 410073, China 3School of Chemical Engineering and Analytical Science, Manchester Institute of Biotechnology, University of Manchester, Manchester M13 9PL, UK 4National Supercomputer Centre in Guangzhou and Research Institute on Application of High Performance Computing, Sun Yat-Sen University, Guangzhou, Guangdong 510006, China


Introduction
A great deal of effort has been devoted to investigating instabilities of non-Newtonian fluids over the past decades [1].The instabilities occurring during the extrusion of polymeric fluids are called extrusion instabilities.The existence of extrusion instability would generally reduce the quality of chemical products.Therefore, the researchers have tried hard to characterize the nature of such instabilities.Extrusion instability is firstly observed as small amplitude, shortwavelength periodic distortions on the surface, which is named as sharkskin.With the increase of flow rate, the surface profile of the extrudate might alternate between rough and smooth, which is called spurt or slip-stick.At a higher flow rate, the extrude surface becomes relatively smooth and experiences a long-wavelength distortion.A further increase of the flow rate would result in a gross distortion of extrudate surface.The last feature is usually called wavy flow or gross melt fracture.A typical flow curve about extrusion instabilities with different flow regions can be seen in Figure 1 [2][3][4].Note that flow instabilities are sensitive to the type of polymer materials.For instance, the spurt effect may be observed in linear low density polyethylene (LLDPE), but it would never occur in low density polyethylene (LDPE).Nonetheless, gross melt fracture exists in almost all unstable extrusion cases.
The origin of the extrusion flow instability is widely studied since its first observation; however, it is still controversial at present.Among various theories developed so far, "slip on the surface of the die" [1,3,5,6] and "effects of constitutive instability" [7][8][9] are two well-known explanations.Slip theory attributes the observation that there is a sudden jump of fluid velocity near the wall to the fact that the polymer molecules fail to adhere to the wall once the shear stress near the wall exceeds a critical value.Therefore, special slip-law is integrated into numerical simulations.However, these methods usually fail to simulate asymmetrical oscillations such as wavy flow.Constitutive instability theory regards the instability as intrinsic natures of the polymer itself and describes it with a nonmonotonic constitutive equation.This point of view could be validated by successful numerical simulation of spurt flow or wavy flow without slip boundary conditions imposed on the wall.As a matter of fact, constitutive instabilities are due to the nonexistence of steady state solutions for the flow in the negative slope regime.Several comparisons have been made between these two theories [4,10].The constitutive instability theory is criticized by some researchers [10] for its failure in reproducing some of the experimental phenomena, but probably it is the inaccuracy of the constitutive equation that should be blamed.Besides, it has made great achievements in the simulation of shear banding [11][12][13][14] which is another kind of polymeric instabilities.
When the behaviour of viscoelastic fluids under complex conditions is simulated, the constitutive model plays a vital role in reproducing experimental results.The FENE-CD model is one of the dumbbell models which is widely used in viscoelastic fluid simulations.In order to obtain a better description of the slipping effect, the upper-convected time derivative in the FENE dumbbell model is replaced by the Johnson-Segalman derivative [15] which is named hereafter as FENE-CD-JS.It is generally used to simulate the behaviour of dilute and semidilute polymer solutions.
Previous researches about the extrusion instability were usually based on mathematical analysis of the constitutive equation [16] or relied on simulations without the extrudate region [3].These two approaches successfully explained and verified the possible reason of extrusion instability to some extent.However, their results did not take the extrudate part into account.The challenge in extrusion flow simulations is how to track the severely distorted flow surface.Though some researchers also showed the feature of the extrudate oscillation [17], they revealed neither the flow behaviours in the extrudate part, nor the irregular oscillations as wavy flow observed in experiments.
This paper aims at studying the original cause of extrusion instabilities.The relationship between the wavy extrudate flow outside the die and the oscillation flow inside is analysed.Both symmetrical and asymmetrical oscillations in and out of the die are observed with the modified FENE-CD-JS model, which indicates that the constitutive instability might be the original cause of extrusion instabilities.Moreover, the results imply that the singularity at the die exit may generate or magnify oscillations.
The rest of the paper is organized as follows.Details of the proposed method are presented in Section 2. The numerical results and discussion are reported in Section 3. Conclusions are drawn in Section 4.

Methodology
In this section, the mathematical formulation of FENE-CD-JS model and numerical algorithms used in this paper will be briefly introduced.

Mathematical Model.
The dynamics of the isothermal and incompressible viscoelastic fluids can be governed by the continuity equation and the momentum equation: where  is the fluid density, u is the fluid velocity, w is the velocity of the cell [18,19],  is the pressure, and   is the solvent viscosity.  is the polymer contribution to the stress tensor which can be calculated from a constitutive equation.FENE-CD-JS model is a constitutive equation which considered the behaviour of finite extensible polymer chains.It could reproduce a range of experimental phenomena, but it still remains a relatively simple mathematical form.It can be written as where  is the slip parameter, which reflects the nonaffine response of polymer chains when they are deformed. is the characteristic relaxation time.D is the rate of deformation tensor which can be expressed as D = (∇u + (∇u)  )/2.I is the unit tensor.A = ⟨Q Q⟩ is a conformation tensor in which Q is the end-to-end vector of a polymer molecule and the symbol ⟨⋅⟩ represents an ensemble average.The time evolution of the conformation tensor is given by (tr(A)) is a function of the trace of the tensor A, and it defines a nonlinear coupling between the applied macroscopic flow field and the evolution of the polymer microstructure.By introducing the dumbbell extensibility parameter, , and a new model parameter, , it can be written as The generalised no-linear spring force may be expressed as F = (tr(A))Q, where  is related to the shortest relaxation time  0 and a constant drag coefficient  0 .The term [1−+√(1/2)tr(A)] can be regarded as the conformationaldependent friction coefficient, which is based on the simple scaling friction law by taking hydrodynamic interactions between fluid and isolated chain molecule into account [20,21]. varies from 0 to 1 which is used to tune the spring stiffness.Moreover, the increase of  would result in a decrease of the force F under the same tr(A) [22].The polymeric stress tensor can be written as where  is the number of polymers per unit volume, and   is the polymer viscosity.
For the convenience of comparison, dimensionless variable Reynolds number Re and Weissenberg number Wi are introduced as in which  is the characteristic length (the half width of the die in this work) and  is the characteristic velocity (the average velocity in this work).

Implementation of Free Surface Boundary Condition.
The boundary conditions imposed on the free surface are comprised of the kinematic boundary condition and the dynamic boundary condition.The kinematic boundary condition can be expressed as It means the relative speed between the fluid and the grid perpendicular to the interface is zero.This has been implemented by Tuković and Jasak [23].
In terms of the dynamic boundary condition, it requires that the traction exerted from the fluid to the atmosphere should be equal to the atmospheric forces imposed on the fluid.Supposing the atmosphere pressure is zero, the boundary condition is given by where  is the total stress tensor, namely, the Cauchy stress tensor across the free surface, and T = 2D +   is the extra stress tensor.Unlike in FEM systems, this boundary condition is hard to satisfy directly in FVM systems.Noticing that velocity is stored in the cell center and the usage of the traction boundary condition here is only to calculate the velocity of cells which own the free surface boundary, we could introduce a balance force  * = − on the free surface and zero elsewhere (including the internal field and other boundary fields) to solve this problem.It is known that −∇ + ∇ ⋅   +   Δ = ∇ ⋅ , and then by applying the operation ∇ ⋅  * in Gauss method and adding it into the original momentum equation, the equation is actually solved as where  is the volume of the cell, S is the face area vector, and n is the unit normal vector of face area. represents faces around a cell, and  represents free surface faces.It can be seen that the influence of the total stress imposed on the free surface is offset in this way, which means the traction boundary condition is satisfied equivalently.The momentum equation is actually solved in the form where Γ  indicates the term is only applied on the free surface boundary.A proper pressure boundary condition on the free surface is required when coupling the velocity and the pressure in SIMPLE algorithm.According to the force balance in the normal direction on the interface, it can be derived as

Preserving the Physical Soundness of Molecular Conformation.
In the original FENE-CD-JS model, there only exists a maximum extension of polymer conformation as tr(A) could not exceed the maximum extensional length .However, there should be a minimum length (or volume) as well to prevent tr(A) from becoming negative.Therefore, before solving the constitutive equations, the state of current polymer conformation should be checked.Modification of tr(A) by setting it back to the initial value would be done if numerical underflow occurred.This kind of modification was proposed by Yang [24] and had achieved a great success in the simulation of viscoelastic flows in complex situations.

Other General Numerical Schemes Used in the Simulation.
In addition to the approaches described above, the discrete elastic split stress (DEVSS) numerical strategy, which is proposed by Matallah et al. [25], is used to enforce the stability of numerical simulation.Moreover, underrelaxation algorithm is also used with the relaxation parameters 0.3 and 0.5 for  and , respectively.Second-order backward Euler method is used for time derivatives.The discretized linear systems are solved by PCG (Preconditioned Conjugate Gradient method) with AMG (Algebraic Multigrid method) preconditioning used for pressure, and BiCGstab (Biconjugate Gradient Stabilized method) with Diagonal Incomplete Lower-Upper (DILU) preconditioning for velocity and stress [26,27].

Results and Discussion
3.1.Model Analysis under Rheometric Flow. Figure 2 shows the rheological response in simple shear and uniaxial extensional flows for the FENE-CD-JS model.Note that the viscosity is normalised by the zero-shear viscosity.The upper lines are the uniaxial extensional viscosities, while the lower ones are the shear viscosities.The shear viscosity of the FENE-CD-JS model undergoes a shear thickening or shear thinning.In terms of the elongational viscosity, it is bounded and would reach a plateau value at last.It means the relationship between the strain stress and the strain rate is linear and the coefficient is constant after a specified strain rate.The solid lines and the dashed lines in Figure 2 illustrate the parameter  has a significant influence on both the shear viscosity and the extensional viscosity.Furthermore, larger  would show a higher extensional viscosity.In the case with a large  ( = 1), the shear viscosity experiences a shear thickening at first and then shear thinning.The shear-thickening effect would be eliminated gradually as  decreases, and the curve of  = 0.1 only shows shear thinning.The parameter  is mainly used to control the shear thinning.A larger  would result in a more pronounced shear thinning.From the solid and the dotted lines in Figure 2, it can be seen that the shear-thickening effect is weakened or even disappeared with the increase of .The effect of the viscosity ratio on the extensional viscosity is almost negligible as shown by the solid lines and the dashand-dot lines in Figure 2. The shear viscosity with a larger viscosity ratio is decreased due to relatively smaller solvent viscosity.
The FENE-CD-JS model reproduces various flow regimes as the rheological flow curves shown in Figure 3.The solid lines in Figure 3  setting, that is,  = 1,  = 0.07, and   / = 1/9, the prediction of the first normal stress is almost the same but the absolute value is a little bit higher.The shear stress varies nonmonotonically with the increase of the shear rate.With further decrease of the viscosity ratio,   / = 1/50, the shear stress of the FENE-CD-JS model exhibits even more pronounced nonmonotonical behaviour as shown by the dotted line in Figure 3.Such a nonmonotonical behaviour of constitutive model could give three shear rate states under one shear stress; hence it is capable of generating constitutive instabilities.The constitutive instability in both Poiseuille flow and die swell flow is studied in the following sections so as to shed light on the origin of extrusion instabilities.

Basic Geometric Shape and Boundary Conditions.
A sketch of the geometry and the boundary conditions is shown in Figure 4.The half width of the die  is chosen as the characteristic length.The length of the upstream and the downstream is set as  1 and  2 , respectively.They should be long enough to ensure the flow before the entry of the die exit and further downstream at the outlet is fully developed.In this paper,  = 1,  1 / = 1546, and  2 / = 16.A uniform inflow is imposed at the inlet.No-slip boundary condition is set at the wall and zero derivative is set at the outlet and free surface for velocity.The pressure at the outlet is set as zero, while on the free surface, it is initially set as zero and would change dynamically during the simulation.The conformation tensor (A) at the input is set as the unit tensor and the corresponding polymer tensor is zero.The computational mesh is shown in Figure 5, in which it can be seen that uniform grids were regenerated.Simulations possess 81600 (680 * 120) cells in total with minimum grid relative interval spaces at Δ min / = 0.025 and Δ min / = 0.017.Flow behaviours in a strait pipe flow are also studied with minimum relative interval spaces at Δ min / = 0.025 and Δ min / = 0.0032.

Inlet Outlet
No-slip Free surface   Since flows with monotonic shear stress (the solid line in Figure 3) are extremely stable in both Poiseuille flow and extrusion flow, only the behaviour of flows with nonmonotonic shear stress is analysed in the following sections.

Simulations of the FENE-CD-JS Model in Poiseuille
Flow.The effects of the constitutive instability are firstly investigated in Poiseuille flow.For the model parameters given as the dashed curves in Figure 3, the extrusion flows under nominal flow rates Wi = 0.1, Wi = 0.5, Wi = 1, and Wi = 3 are investigated.The actual local Wi number at the wall, which is evaluated by multiplying the shear rate of the cell most close to the wall and the relaxation time after simulation, equals 0.3, 7.7, 8.9, and 13.6, respectively.Therefore, the first case is located in the smooth region, while the latter three cases are located in the wavy region as shown in Figure 1 (the Wi of local bottom point B is corresponding to 7.2).The flow of the first case remains stable.But under the latter three conditions, they begin to oscillate.For the case of Wi = 0.5, there are only very small oscillations occurring near the wall which almost does not affect the central flow, while for Wi = 1, the flow changes gradually from stable to asymmetrical oscillating.The intensity images of tr(A) in time sequence are shown in Figure 6.It can be seen clearly that the oscillation occurs near the wall where polymer molecules are stretched most.Numerous experimental and numerical simulations have proved that elastic stress tends to destabilize the fluid.For example, instabilities of pipe flow were observed at a lower Re in viscoelastic fluids (about 2230 or even lower) than those in Newtonian fluids (about 5800) [16,28].
Pipe flows under model parameters of the dotted lines in Figure 3 are studied under nominal flow rate of Wi = 0.1, Wi = 0.4, Wi = 0.5, and Wi = 3.By calculating the shear rates after the simulations, the local Wi numbers at the wall are around 0.3, 1.2, 15.8, and 31.8,respectively.The first two cases are located in the smooth region and display completely stable status as expected.In contrast, the latter two cases are located in the spurt and wavy flow regime, and both of them exhibit oscillations.Note that, at Wi = 0.4, the shear rate near the wall has been very close to the local top (point A in Figure 1) of the shear stress curve.Thus, after a slight increase of the nominal flow rate, the flow near the wall would run into the spurt regime and experience a sudden change of the flow rate.Figure 7 illustrates that there is a sudden velocity jump near the wall at Wi = 0.5, in contrast to the case of Wi = 0.4.It is trustworthy that the sudden jump of the shear rate is the direct cause of extrusion instability.Furthermore, the position of the discontinuity is very close to the wall, which has also been discovered in experiments that the polymer could form a thin layer near the wall with large velocity gradient called the spurt layer [7].Similar simulation results are also observed by Fyrillas et al. [8] with Johnson-Segalman model.By drawing all cells in the middle region of the pipe and averaging over the period of last 5 relaxation times of each cell point, a comparison of the relationship between shear rate and shear stress is shown in Figure 8.The error bars in the horizontal direction and vertical direction are the standard deviation of shear rate and shear stress, respectively.The markers represent the average values of each cell during this period.It can be seen clearly that the shear stress linearly increases with respect to the shear rate for the case of Wi = 0.4 as it is located in the smooth region.In contrast, for the case of Wi = 0.5, the shear rate near the wall might depart from the smooth flow regime and jump to a much higher value at nearly the same shear stress.This is due to the existence of negative region in its shear stress flow curve, by which the spurt effect happens.In addition, the length of the error bar for points with higher shear rate is much longer than the ones with lower shear rates, which means the spurt layer is particularly unstable.For the case with Wi = 3, the flow is totally in chaos as shown in Figure 9.This might explain why the polymer would distort severely in gross melt fracture.The figure exhibited by the diamonds in Figure 8 shows that the error bars for Wi = 3 (red bars) are much longer than that for Wi = 0.5 (blue bars) which obviously indicates the oscillation is more serious.

Simulations of the FENE-CD-JS Model in Die Swell.
The effects of constitutive instability in free surface flow are further studied to reveal the possible cause of extrusion instabilities.For flows with nonmonotonic constitutive model (the dashed curves in Figure 3), flow behaviours of die swell process with different flow rates are given in Figure 10.The local Wi number at the wall for the nominal flow rate of Wi = 0.8 is about 9.1.At Wi = 0.1, both flows in the die and in the plug domain are stable.With the increase of flow rate to Wi = 0.5, instabilities emerge and the extrudate oscillates with a small amplitude; however, the flow in the die is still almost stable.By increasing the flow rate further to Wi = 0.8, periodical and symmetry slight oscillations begin to appear in the die, and spurt appears near the wall of the die.In terms of the extrudate part, it oscillates in the same  manner as the flow in the die with a short wavelength.It seems that the wavelength of the periodical oscillation in and out of the die preserves the same order of magnitude, which indicates that there might exist a close relationship between the instabilities in the inside and outside of the die.One possible explanation is that the extrudate instability comes from the unstable flow in the die, and the other one is that the instability is caused by the singularity of the die exit and propagates towards both sides.To increase the flow rate to Wi = 1, more severe distortions are observed.Moreover, the flow is no longer oscillating symmetrically but forms waves.
The shear rates at the wall vary drastically approximately between 70 and 20 when approaching the die exit for both conditions at Wi = 0.8 and Wi = 1, primarily because they are located in the wavy or melt fracture region.Concerning the extrudate domain, it is still smooth at the surface but distorts with a larger amplitude, and the oscillation is also no longer symmetric compared with the previous results.The oscillation of flow in the die impels the movement of the extrudate.Thus, in this case, oscillations in the die could be considered as the driving cause of the extrudate instability.Palza et al. [29] also agreed with the point that spurt began at the entrance of the die and propagated from upstream to downstream along the die by analysis of physical experiments with a long die.Combined with the performance in the straight pipe described in Section 3.3.1,instability in the die which is caused by constitutive instability is more likely to be the origin of wavy flow.The difference between straight pipe flow and extrusion flow is that extrusion flow is more likely to be unstable.This could be attributed to the influence brought by the die exit singularity and the extrudate part.The existence of singularity would enhance the shear rate to spurt region quickly.Moreover, once the extrudate begins to oscillate, it would influence the upstream flow in return and prevent the flow from going back to the stable status.
For the case with a deeper negative region (the dotted curves in Figure 3), a specific simulation of die swell at Wi = 0.5 is shown in Figure 11(c).It is found that the velocity profile is significantly different from cases mentioned above.Although the plug flow region is extremely uniform and stable, the flow in the die oscillates symmetrically.Spurt occurs near the wall periodically which might suggest that spurt effect is the initial cause of instability and is reproduced by constitutive equations.More interestingly, the visual effect    with the time passing by looks like the position of the spurt is moving backward as displayed in Figures 11(c)-11(f).From these time serial figures, it could be seen clearly that the blue region near the wall and the red region in the middle of the die are both moving backward from the die exit to the input region, which is a counterintuitive finding.In order to demonstrate that this is not an illusion caused by frame rate, the generation of the first spurt appearance at the wall and first ellipsoid-like phenomena of the central flow are given in Figures 11(a) and 11(b), respectively.From these two figures, it could be seen clearly that the instability firstly occurred near the wall at the die exit and then affected the running of the central flow.Instabilities were apparently propagating backward.It is reasonable that spurt firstly happened near the die exit since polymers are always stretched most severely and the shear rate is also the highest near the singularities.Time and mesh density convergence tests were implemented, and the results were the same.This counterintuitive phenomenon might suggest that one kind of instability is related with the singularity at the die exit.Considering the fact that slight oscillation was also observed in pipe flow, it is reasonable to believe that the singularity generated or magnified the instability in extrusion flows.This point may be supported by the experimental effect that periodical oscillations close to the exit were found [30] and it owned higher relative pressure fluctuations than other upstream positions of the die [29].Although spurt phenomenon is seen in the simulation, no obvious evidence has been proved that slip effects in experiments are also an intrinsic property of constitutive instability and could be integrated into a constitutive equation.Moreover, sharkskin was not seen in the extrudate region in the simulation.Therefore, this result could only be used as an implication that sudden jump of velocity at wall could be reproduced by a proper constitutive equation, but not an obvious evidence that sharkskin is also originally caused by constitutive instabilities.

Conclusion
In this work, the original cause of extrusion instability was discussed systematically with the modified FENE-CD-JS model.By carefully studying different flow regimes, the conclusion may be drawn that nonmonotonic constitutive model reproduced wavy flow, which supported the theory that the constitutive instability may be the origin of the extrusion instabilities in spurt and wavy flow.Some welldocumented experimental phenomena (symmetrical and asymmetrical oscillations in and out of the die) were reproduced in this paper, which provided numerical evidences to support the constitutive instability theory.Detailed analysis showed that wavy flow may be attributed to the instability emerging in the die.Moreover, experiments have shown that spurt instabilities may be generated or magnified at the die exit and propagate to both sides.In future work, more constitutive equations with worm-like rheological flow curves can be tested so as to further investigate the influence of nonmonotonic shear stress on extrusion instabilities.

Figure 3 :
Figure 3: The normalised first normal stress and shear stress for the FENE-CD-JS model in a simple shear flow with various parameters.

3. 3 .
Results and Discussion.In this section, the relationship between wavy flow and constitutive instability is investigated by using the FENE-CD-JS model.The rheometric behaviour of FENE-CD-JS model with different parameters has been illustrated in Section 3.1.The model is firstly studied under Poiseuille flow to reveal possible flow instabilities.Afterwards, die swell simulations are implemented to study the relationship between wavy flow and constitutive instability.