Multiscale Numerical Simulations of Branched Polymer Melt Viscoelastic Flow Based on Double-Equation XPP Model

The double-equation extended Pom-Pom (DXPP) constitutive model is used to study the macro and micro thermorheological behaviors of branched polymermelt.The energy equation is deduced based on a slip tensor.Theflowmodel is constructed based on a weakly-compressible viscoelastic flowmodel combined with DXPPmodel, energy equation, and Tait state equation. A hybrid finite element method and finite volume method (FEM/FVM) are introduced to solve the above-mentioned model. The distributions of viscoelastic stress, temperature, backbone orientation, and backbone stretch are given in 4 : 1 planar contraction viscoelastic flows. The effect of Pom-Pom molecular parameters and a slip parameter on thermorheological behaviors is discussed. The numerical results show that the backbones are oriented along the direction of fluid flow in most areas and are spin-oriented state near the wall area with stronger shear of downstream channel. And the temperature along y = −1 is little higher in entropy elastic case than one in energy elastic case. Results demonstrate good agreement with those given in the literatures.


Introduction
Branched polymer becomes more and more concerned because of its unique structural characteristics and properties, now its development is one of the fastest in macromolecular materials.Branched polymer has more complex thermorheological behavior compared with other polymers, and its rheological behavior depends on its topological structure of branched molecules [1,2].As compared to the linear polymer, when the main chain of branched polymer introduces a certain number and length of branched chains, the viscoelasticity is significantly different.Branched polymer in shear flow shows the similar strain softening but has a longer relaxation time at the end of branched molecular chains because of the limitation of branched chain.Moreover, branched polymer in elongation flow has entirely different strain softening.Therefore, branched has a great influence on polymer viscoelastic properties.
In recent decades, some researchers have developed many viscoelastic constitutive models for describing the rheological behavior of polymer based on different theories [3].
Among them, the model based on molecular theory can more truly reflect the rheological properties of fluid and can more fully reflect the flow of the fluid [4].And for all we know, a branched polymer melt can be considered as a melt in which a certain concentration of branched molecules is embedded in a viscous melt.In this, Mcleish and Larson [1] proposed a Pom-Pom model based on Doi-Edwards's peristaltic tubes theory.In the Pom-Pom model, they simplified each branched molecule to a molecule with only two branched points at each end, and with a certain number of arms at each branched points.This model is not completely consistent with the topological structure of branched molecules, but it is an important breakthrough in the field of viscoelastic constitutive models.This model introduces the important branched information and distinguishes the orientation relaxation time and extension relaxation time of backbone.It can also study the relaxation time of branched molecules and their effects on the above two relaxation time.Subsequently, Verbeeten et al. [5] improved the Pom-Pom model and proposed an extended Pom-Pom (XPP) model by using the slip tensor.This model overcomes some defects 2 Advances in Mathematical Physics of Pom-Pom model, such as the discontinuity of steady-state stretching, the unrestricted orientation under the high strain, and the unpredicted second normal stress difference.In addition, Clemeur et al. [6,7] proposed a Double Convected Pom-Pom (DCPP) model in order to solve the problem that the solution of XPP model is not unique.However, the DCPP model suffered from numerical instability in the numerical simulation.On the basis of this, Clemeur and Debbaut [8] proposed a modified DCPP model and Wang et al. [9] given the Simplified Modified Double Convected Pom-Pom (S-MDCPP) model with good numerical stability and easy programmable ability.
Generally, there are two kinds of Pom-Pom molecule constitutive models: single-equation model and double-equation model.Due to the simple solution and easy programming of the single-equation model, many studies have used the singleequation XPP model to simulate the viscoelastic flows [10][11][12][13][14][15][16][17], but it cannot describe some micro information.Doubleequation XPP (DXPP) model can describe the microscopic orientation and stretch of backbone and study the influence of microscopic molecular parameters on the rheological behavior of branched polymers.However, due to the complexity of the DXPP model, there are few reports on the numerical simulation of this model.Therefore, DXPP model is used to study the microscopic information of the orientation and stretch of branched molecules in this paper.
In the past twenty or thirty years, the numerical simulation of viscoelastic flow has been developing rapidly and the main numerical methods are finite element method, finite volume method, and meshless method.Although there are many numerical methods, they each have their own advantages and disadvantages.There is no certain method to "dominate the world."There is only one method to solve a problem when appropriate or not.Therefore, the combination of the merits of various methods to form a hybrid algorithm will be a trend of numerical simulation [16,18,19].In this paper, the hybrid finite element method and finite volume method (FEM/FVM) are proposed based on the advantages of finite element method and finite volume method and the characteristics of the solved problem.
In addition, since the actual polymer processing is often a nonisothermal viscoelastic flow problem, the effects of temperature are also considered.The slip tensor of viscoelastic fluid actually affects the energy equation; that is, the energy equation is also different for different slip tensor [20,21].Therefore, we will give the derivation of the energy equation based on the slip tensor and study the influence of the slip parameter on the temperature.
Above all, the DXPP model is used to study the macroand micro-rheological information of branched polymer melt.The energy equation based on the slip tensor is deduced and used to study the influence of slip parameters on the temperature.Subsequently, based on the characteristics of weakly-compressibility and high specific heat capacity of the polymer melt, the hybrid FEM/FVM method is used to solve the above model, and the macro and micro thermorheological properties of the branched polymer are discussed according to the numerical simulation results.

Mathematical Models
where ∇ S ≡ S/ + u ⋅ ∇S − S ⋅ ∇u − (∇u)  ⋅ S denotes the upper convected time derivative of orientation tensor S; D is the rate of deformation tensor; the slip tensor B is defined as where  is a material parameter, defining the amount of anisotropy;  0 is the relaxation time of the backbone tube orientation; the exponential stretch relaxation time   =  0  −](Λ−1) ensures the stretch relaxes very fast and stays bounded for high strains;  0 is the relaxation time for the stretch, and ] = 2/, where  is the amount of arms at the end of a backbone; tr(⋅) is the trace; Λ is the backbone tube stretch and its material derivative Substituting (2) into (1) gives the orientation equation Viscoelastic stress equation is where  0 is the plateau modulus; I is the unit tensor.
In conclusion, ( 3) and ( 4) constitute a DXPP model describing the backbone tube stretch and orientation using two decoupled equations; (5) denotes the viscoelastic stress.Here, the model is extended with a second normal stress difference when  ̸ = 0.By defining   =  0  0 as the viscosity of polymer, We =  0 / as the Weissenberg number, and  =  0 / 0 as the relaxation time ratio, dimensionless DXPP model can be written as We where  is the ratio of Newtonian viscosity to total viscosity and  and  are the velocity and length of the dimensionless parameters, respectively.

Governing Equations.
In the polymer processing, the weakly-compressibility of polymer melt cannot be ignored.Therefore, the weakly-compressible flow conservation equation is used to describe the polymer melt flow.For weaklycompressible viscoelastic flows, the conservation equations for mass and momentum can be expressed as follows, respectively, where Re =   /  denotes the Reynolds number;   and   are the density and viscosity of the dimensionless parameters, respectively.
The energy equations of different viscoelastic fluids also vary due to the different slip tensors.The derivation of the energy equation based on the XPP fluid slip tensor is given below.
The general form of the energy equation based on the slip tensor is as follows: where   is the specific heat,  is the temperature,  is the heat flux,  is the Cauchy stress tensor,   is the material parameter, and their expressions are Substituting ( 2) and ( 12) into (11) gives the energy equation The second term in the right-hand side of ( 13) reflects the contribution of entropy elasticity.The last one reflects the contribution of energy elasticity and  ∈ [0,1].The dimensionless energy equation is where Pe =    0 / 0 is the Peclet number, Br =    2 /( 0   ) is the Brinkman number, and   ,  0 , and  0 are the temperature, specific heat, and coefficient of heat transfer of the nondimensional parameters, respectively, where In addition, a P-V-T equation of state is necessary to satisfy the completeness of governing equations because of considering the compressibility of the polymer melt.Tait state equation [18] is usually considered as the classical empirical equation and is capable of describing both the liquid and solid regions.So Tait state equation is used in this paper.

Numerical Methods
The flow of brand polymer melts is governed by the conservation of mass, momentum, and energy equations, Tait state equation, together with a DXPP constitutive model.The numerical simulation of the above model employs hybrid FEM/FVM [18] method.The momentum equations are solved by the FEM, in which a discrete elastic viscous stress split (DEVSS) scheme is used to overcome the elastic stress instability, and an implicit scheme of iterative weaklycompressible Crank-Nicolson-based split scheme (WCNBS) is used to avoid the Ladyzhenskaya-Babuška-Brezzi (LBB) condition.The energy and DXPP equations are solved by the FVM, in which an upwind scheme is used for the strongly convection-dominated problem of the energy equation.
To analyse the accuracy of the algorithm mentioned above, we construct the DEVSS scheme based on (10) and consider its discretization in the time domain within a typical time subinterval [  ,  +1 ], which give us the form of the Wilson- method as follows: where D is an added variable for constructing DEVSS scheme, 0 ≤ ,  1 ,  2 ,   ,   ≤ 1.
The truncation error of ( 16) is Based on formula ( 17), ( 16) has the second-order accuracy when  =  1 =  2 =   =   = 0.5, which is adopted in this study for the Crank-Nicolson scheme.
Since the energy equation is deduced based on the slip tensor and the viscoelastic stress is calculated using a DXPP model that can describe the backbone orientation and stretch of the polymer molecules, we will detail the solution of the energy equation and the DXPP model based on the nonstaggered grid under the framework of the FVM.The energy equation and the DXPP model can be normalized as follows: where , Γ are constants;  and   are the physical quantities and source term which are defined in Table 1.The terms from left to right in (18) represent the time, convective, diffusive, and source contributions, respectively.The discretization of the energy equation ( 13), orientation equation ( 6), and stretch equation ( 7) can be written as the following form by a generalized quantity ; that is, where   is the source term after discretization in ( 13), (6), and In this paper, the hybrid FEM/FVM method described above is used to solve the weakly-compressible flow model based on the DXPP constitutive model.Details are as follows.
Step 2. Solve the momentum and mass conservation equations to calculate u, , and  under the framework of the FEM.
Step 3. Solve the energy equation and the DXPP constitutive equation to obtain , S, and Λ under the framework of the FVM.Step 4. Use expression (8) to calculate the polymer stress ().
Step 5. Substitute  into the momentum equation (10) to ensure the calculation of the coupling of the physical quantities.

Numerical Simulation and Analysis
The 4 : 1 planar contraction flow is a benchmark test example and has been widely studied [9][10][11].In the 4 : 1 planar contraction flow, the fluid flows into the narrower channel from the wider channel with a simple shear flow far from the contraction region, a pure elongation flow along the central axis, a complex strong shear flow near the wall, and a mixture of shear and elongation flows near the reentrant corners.In fact, the contraction flow widely exists in the processing of polymer materials, such as the polymer extrusion and injection molding.Therefore, this example not only verifies the correctness of the proposed algorithm and model, but also provides the basis for the processing of polymer materials.
The sketches of the lower half of the 4 : 1 planar contraction geometry and computational mesh are shown in Figure 1.The lengths of upstream and downstream channel are both 16, in which  denotes the height of downstream channel.The structured triangular mesh and rectangular mesh are used in FEM and FVM, respectively.It is noted that the mesh is refined in the near of the reentrant corner.
The initial and boundary conditions are as follows.
We choose the polymer High-Density Polyethylene (HDPE) Sclair 2714 made by Nova Chemicals Inc. as fluid.The material parameters of HDPE, which are obtained from the materials database of Moldflow software, are shown in Tables 2 and 3, respectively.A series of meshes is used for the FEM/FVM method to ensure spatial convergence based on salient-corner vortex cell size.Mesh characteristics, detailing numbers of elements (FEM/FVM) or volumes (SLFV) [22] and smallest mesh spacing employed, and salient-corner vortex cell size are   provided in Table 4.Moreover, the quantitative information regarding mesh convergence of the salient-corner vortex cell size is compared to the results in literature [22] for the Re = 0, covering the range 1 ≤ We ≤ 10.The information in Table 4 demonstrates that convergence with mesh refinement has been achieved for the range of parameters considered.Figure 3 shows the numerical results for the stress   along the axis of symmetry with different values of the amount of arms  and relaxation time ratio .The values of stress   increase with the values of  increase and they have obvious changes near the reentrant corner.However, the values of   are almost no change when the  increases to a certain extent.In addition, the values of   under different values of  tend to be consistent when the flow is fully developed.For different values of , the values of   have similar change trends with different values of  except they decrease with the values of  increase.

Numerical Solutions of Temperature.
Figure 4 shows the distribution of the temperature near the reentrant corner at different Peclet numbers when the Weissenberg number is fixed to 1.0.It is observed from Figure 4 that the temperature is lower near the wall, and low-temperature region is getting smaller and smaller with the increasing of Pe number.This is because the effect of heat convection gradually increases with the increase of Pe number.Then it causes the heat transport ratio of heat convection to heat dissipation changes.
The temperature along  = −1 with different slip parameters  is shown in Figure 5.It is seen that the temperature values are slightly higher when  = 1 than those when  = 0; that is, the temperature of the entropy elasticity is slightly higher than that of the energy elasticity.This is consistent with the result of literature [20].The most intuitive way to describe the backbone orientation of branched polymer molecules on the molecular scale is to use the specific information of the second-order orientation tensor S. The ellipse method is adopted to obtain the backbone orientation state for two-dimensional cases.For ellipse method, the eigenvalues and eigenvectors are first obtained by computing the second-order matrix corresponding to S, and then the eigenvectors and eigenvalues represent the major axis's direction and length of the orientation ellipse, respectively.
The backbone orientation in 4 : 1 planar contraction flow is shown in Figure 6.As can be seen in Figure 6, the backbones are oriented along the direction of fluid flow in most areas; the backbones are spin-oriented state near the wall area with stronger shear of downstream channel.This is because abrupt contraction flow area and rapid increase fluid velocity lead to backbone quickly spin near the reentrant corner.Near the wall and away from the area of the reentrant corner, the shear stress is the largest, and the backbone orientates first along the wall and then spins with the flow.Along symmetry axis, the backbone near the reentrant corner is first stretched and exhibits uniaxial tension state due to the velocity gradient increases; then the other backbones are in turn oriented in the horizontal axis.x Λ q = 2 q = 4 q = 8 q = 12 q = 15 the size of angle vortex, the fluid stretching thickening [17], and shear rate distribution.The appearance of this fact will help to understand the phenomenon of polymer wall slip and extrusion instability.Figure 8 shows the numerical results for backbone stretch along the axis of symmetry with different values of the amount of arms  and relaxation time ratio  with We = 10.It is observed that the value of Λ is increasing with the increasing value of  and there are obvious changes near the reentrant corner.However, the value of Λ is almost no change when  increases to a certain extent.In addition, the values of Λ under different values of  tend to be consistent when the flow is fully developed.For different values of , backbone stretches have similar change trends with different values of  except they are decreasing with the increasing values of .

Conclusions
In this paper, the DXPP constitutive model, which describes backbone orientation and stretch, is used to study the thermorheological behaviors of branched polymer melt.The hybrid FEM/FVM method is used to solve the nonisothermal weakly-compressible viscoelastic flow model coupled with a DXPP model.The distribution of viscoelastic stress, temperature, and backbone orientation and stretch are given.The effect of Pom-Pom molecular parameters and a slip parameter on thermorheological behaviors is studied.All numerical results can prove the models and numerical methods mentioned above are valid.
For 4 : 1 planar contraction flow, the backbone orientation is along the flow direction in most of contraction area and is spin in the downstream stronger shear near wall area.Stress increases with the increase of the value of  and decreases with the increase of the value of .The backbone stretch increases with the increase of the values of We number and , and it decreases with the increase of the values of .The variable trend of stress and backbone stretch for the different values of  is the same with the full development of polymer melt flow.In addition, the temperature along the centre line is little higher in entropy elastic case than one in energy elastic case.The macroscopic thermal rheological behavior in the flow field is a true reflection of the microscopic topological structure of the polymer melt.The characterization of the microscopic information helps to further study the flowinduced residual stress and other complicated behaviors in the process of polymer melt processing and provide a theoretical basis to improve the polymer product performance.

4. 3 .
Numerical Solutions of Backbone Orientation.The Pom-Pom molecular model describes the relaxation time of branched macromolecules separately on two different time scales and introduces backbone stretch parameters to describe the tensile behavior.Using the DXPP model allows one to investigate the complexity rheological behavior on the molecular scale.

Figure 3 :Figure 4 :
Figure 3: The influence of different parameters on   : (a)  and (b) .
. The coefficients   ,   ,   ,   , and   can be expressed as the combination of the convection term and the diffusion term; that is, The form of (| Δ |) can be different under different discretization schemes for the convection term.In order to solve the convection-dominated problem caused by the high specific heat capacity and high Weissenberg number, (| Δ |) equals 1 for the upwind scheme in this paper.All the coefficients are formulated as follows.

Table 3 :
Property parameters of HDPE.The numerical results for the first normal stress   , first normal stress difference  1 , and second normal stress difference  2 near the reentrant corner are illustrated in Figures2(a), 2(b), and 2(c), respectively.It is seen that the stress contours near the reentrant corner are smooth and  2 is not zero.This proves that the given FEM/FVM method is feasible.