A Dissipative Particle Dynamics Study of Flow Behaviors in Ultra High Molecular Weight Polyethylene/Polyamide 6 Blends Based on Souza-Martins Method

This paper presents our study on the use of dissipative particle dynamics (DPD) simulations to discover the flow behavior in ultra high molecular weight polyethylene/polyamide 6 (UHMWPE/PA6) blends associated with extensional-shear coupled flow, based on the Souza-Martins method, for the first time. By way of simulations, we aimed at investigating the mesoscopic morphology and alignment behavior in response to extensional-shear coupled flow, in comparison with simple shear flow and simple extensional flow. Our results reveal that the aggregation of polymers is noticeable under zero flow, as expected. Within the considered range of extensional-shear coupled rates, the morphology transforms from micelle-like clusters to a chain-like network structure by increasing coupled rates from 0.01 to 2.0. Furthermore, it shows a linear distribution along the flow direction at a high coupled rate. It can be concluded that the flow behaviors in UHMWPE/PA6 blends are significantly impacted by extensional-shear coupled rates. The orientation behavior induced by extensional-shear coupled flow is more obvious than shear flow, even though flow variations and mass fractions yield less effects on the distribution behaviors of UHMWPE/PA6 blends. The DPD results are verified by mean square displacement (MSD) as a function of simulation time and relative concentration distribution along Z direction.


Introduction
It is well-established that morphology and interfacial behavior play key roles in the rheological properties, when blending or mixing immiscible polymer melts, and thus affects the final mechanical properties of the blends [1,2]. Hence, it is vital to elucidate the connection between applied flow fields and morphology evolutions for optimization of the processing, and therefore, the resulting properties of blends, taking interfacial tension into account [3,4].Various thermodynamic variables change continuously over this interfacial region and the thickness of this interfacial region is associated with the finite range of molecular interaction. The finite interaction range is characterized by the free energy of the system, which is dependent not only on the local composition but also on the composition of the immediate environment [5]. For immiscible polymer pairs with infinite molecular weights, interfacial tension can be calculated according to the Flory-Huggins interaction parameter (χ) between the corresponding monomers; e.g., A and B, given the assumption of complete immiscibility [4,6]. Mixing performance is related to the energy to mix, reflected in overall pressure drop for all designs [7]. Small chains cluster at the interface, resulting in the reduction of the interfacial tension and the Gibbs free energy of the system [4,8].
Factually, for ultrahigh molecular weight polyethylene (UHMWPE) samples prepared by eccentric rotor extrusion [9][10][11], the oriented structure has been experimentally observed and evidenced by Raman spectra [12]. However, it is plagued by experimental difficulty in studying the dynamic process of chain conformation. Compared with experimental study, the unique advantage of computational simulation, e.g., coarse-grained molecular dynamics (CGMD) and dissipative particle dynamics (DPD) simulation, lies in directly recording, tracking and observing the dynamic evolution process in real time, offering a complement to experimental studies.CGMD and DPD are coarse-grained simulation techniques based on molecular dynamics (MD), where each beadrepresents several atoms or repeat units in MD [13,14]. Even though they are employed for studying systemswithlarger lengths and time scales than classical MD dynamics [15][16][17], only in DPD simulation the interfacial properties can be taken into consideration. For engineering applications, mixing UHMWPE with polyamide 6 (PA6) is an effective and convenient route to reduce the melt viscosity.In our earlier report [18], great efforts have been made to compare the complex phase morphology and alignment evolution of the immiscible UHMWPE/PA6 blends in response to extensional and shear flow by means ofDPD simulation, in which the extensional flow and shear flow were separately imposed. During the actual process of eccentric rotor extrusion, the polymer melts encounter an extensional and a shear deformation at the same time; that is, they undergo extensional-shear coupled flow. The extensional-shear coupled flow can be realized by adjusting the magnitude of extensional loading in Z direction and X-Y shear loading by utilizing finite element simulation [19]. Finite element method is a typical success of continuum mechanics in predicting materials' response and failure at macroscopic level. To extend the concepts of continuum mechanics, i.e., stress and strain, to the nano-scale for exploring the microscopic behavior, attracts increasing attention [20,21].Inspired by this, we make an attempt to construct the extensional-shear coupled flow in DPD simulation, for the first time. In Reference [18], based on the Souza-Martins method [22], the extensional flow was driven by an external pressure along the length L z of the simulation box. Simultaneously, the Lees-Edwards boundary condition [23] was first proposed in 1972 for achieving a linear velocity profile, and then it was used for producing a shear flow in Y direction with the velocity gradient in X direction [24]. Likewise, we can establish shear flow utilizing the Souza-Martins method, similar to the extensional flow, in which external pressure along the lengths L x and L y of the simulation box are imposed on each bead with identical magnitude, respectively. In this case, the extensional flow and shear flow was simultaneously imposed on each bead, being consistent with the actual process of eccentric rotor extrusion. The Souza-Martins method is found to function more effectively under an extensional-shear coupled flow field, than Lees-Edwards boundary conditions. Hence, in this contribution, we present a study on the use of DPD simulations to investigate the flow behavior in UHMWPE/PA6 blends subjected to extensional-shear coupled flow by varying extensional-shear coupled rates, in comparison with simple shear flow and simple extensional flow, based on the Souza-Martins method.
where a is the statistical segment (monomer) length, k B is the Boltzmann constant and T is the absolute temperature. The χ of UHMWPE and PA6 pairs can be estimated from the change in energy during mixing per unit volume using Equation (2) [25] and Equation (3) [26].
where φ i is the volume fraction and an implicit condition relating to this equation is that the lattice is filled completely, namely, φ A + φ B = 1. V m and V are the average molar volume of the beads and unit volume of polymer system, respectively. E coh is the calculated cohesive energy, which can be directly calculated by MD simulation as depicted in detail in our earlier work [18]. Then the repulsive parameter α ij between different types of beads is calculated according toχ, see Equation (4). When i equals j, the repulsive parameter α ii between two identical type beads can be obtained using Equation (5): where the bead number-density, ρ, is chosen as 3 and k B T is the conservative interaction potential chosen as 1.

Construction Basis in Flow Field
As mentioned above, in the actual process of eccentric rotor extrusion, the polymer melts undergo a combination of extensional deformation and shear deformation; see Figure 1. The χ of UHMWPE and PA6 pairs can be estimated from the change in energy during mixing per unit volume using Equation (2) [25]and Equation (3) [26].
where i φ is the volume fraction and an implicit condition relating to this equation is that the lattice is filled completely, namely, A V and V are the average molar volume of the beads and unit volume of polymer system, respectively. coh E is the calculated cohesive energy, which can be directly calculated by MD simulation as depicted in detail in our earlier work [18]. Then the repulsive parameter ij a between different types of beads is calculated according toχ, see Equation (4). When i equals j, the repulsive parameter ii a between two identical type beads can be obtained using Equation (5): where the bead number-density, ρ, is chosen as 3 and kBT is the conservative interaction potential chosen as 1.

Construction Basis in Flow Field
As mentioned above, in the actual process of eccentric rotor extrusion, the polymer melts undergo a combination of extensional deformation and shear deformation; see Figure 1. For simple shear flow, the velocity gradient ( ) t γ  is given by Equation (6) [27]: For simple extensional flow, the velocity gradient ( ) t ε  can be expressed as [28]: .
Polymers 2019, 11, 1275 4 of 12 The velocity gradient ∇u for a combination of shear flow and extensional flowhas the following equation [29]: In Materials Studio, pressure is another basic thermodynamic variable that defined as the force per unit area and the general form is shown in Equation (9): Each element of the tensor is the force applied on the surface of an infinitesimal cubic volume with edges parallel to the x, y, and z axes. The first subscript represents the direction of the normal to the plane and the second subscript refers to the component of the force.
In an isotropic simulation, the off-diagonal elements are zero and the diagonal elements are equal, where the forces are the same in all directions and there is no viscous force, see Equation (10): where the scalar quantity, p, is the equivalent hydrostatic pressure. Sometimes, especially in materials science studies, the diagonal elements are acknowledged as the tensile or normal stress and the off-diagonal elements are known as the shear stress. Both pressure and stress are input and reported in GPa in Materials Studio.
Consequently, in Materials Studio, simple shear flow (P S ), simple extensional flow (P E ) and extensional-shear coupled flow (P E−S ) can be described as Equations (11)-(13), respectively, in which the equivalent hydrostatic pressure p is zero.
P xx P xy P xz P yx P yy P yz P zx P zy P zz P xx P xy P xz P yx P yy P yz P zx P zy P zz P xx P xy P xz P yx P yy P yz P zx P zy P zz

Simulation Details
On the basis of Equations (2)-(5), the repulsive parameter α ij between UHMPWPE and PA6 is 38.18, 52.08 and 57.83 in UHMWPE/PA6 blend with mass fractions of 70/30, 50/50, 30/70, respectively. More details can be found in our earlier work [18].The simulation was carried out using the Mesocite module in Accelrys Materials Studio software. The extensional-shear coupled flow is generated by imposing a force on each bead in three directions, in which the pressure and stress are controlled by Souza-Martins method, according to Equations (11)- (13).Initially, the simulation was accomplished in a cubic box with a bead mass of 685 amu and length scale of 12.8 Å in reduced DPD unit. The coarse-grained model consists of 42,521 beads (11,596 beads for every UHMWPE chain and 51beads for every PA6 chain).Initial condition for the simulation: DPD simulation of UHMWPE/PA6 blends were performed under zero flow after 1.2 × 10 5 steps geometry optimization and 1.0 × 10 5 steps dynamic simulation, with a time step of 1fs in reduce DPD unit to ensure that the system reached a complete equilibration, judging from the individual stable stage in the time evolution of pressure and total potential energy of the system (see Figure S1).

Results and Discussion
In this work, Souza-Martins method is performed for driving the shear flow field, instead of Lees-Edwards boundary conditions. The Lees-Edwards boundary conditions have been commonly applied in MD to generate a shear flow with a constant shear rate and extended into DPD to eliminate the wall effect, in which the periodic boundary conditions are applied for the whole domain, including the top and bottom boundaries. A constant force was applied on all beads to generate simple shear flow, and the shear rates varied from 0.01 to 2, where . γ = 0.01-2, . ε = 0 (P yx = P xy = 1.0 × 10 −5 -2.0 × 10 −3 GPa and P xx = P yy = −1/2P zz = 0 GPa). In this subsection, the extensional flow also has been taken into consideration and it brings an significant effect on the morphological behaviors of UHMWPE/PA6 blends, under a range of extensional rates, in which . γ = 0, . ε = 0.01-2 (P yx = P xy = 0 GPa and P xx = P yy = −1/2P zz = −1.0 × 10 −5 -2.0 × 10 −3 GPa). Figure 3 presents the typical equilibrium structures as a function of extensional rates . ε. It shows almost the same structural evolutions as the condition of shear flow when the extensional rate is low (i.e., . ε = 0.05). With an increased . ε of 1.0, the PE and PA6 molecules in the blended system are elongated and exhibit strong tendency to parallel align along the flow direction. By increasing . ε to 2.0, the micelle-like structures finally show linear distribution and transforms to chain-like network structures.
(d)γ = 0.5, ε = 0 (e)γ = 1.0, ε = 0 (f)γ = 2.0, ε = 0   To further investigate the effect of extensional-shear coupled rate on the morphology and distribution in the blended UHMWPE and PA6, a combination of the above simple shear flow and extensional flow was then applied. ε, the molecular chains are highly oriented. The distribution of PE and PA molecules are uniform by changing coupled rates, maintaining a good compatibility between PE and PA6. Based on the analysis above, it can be concluded that the coupled rate applied to the UHMWPE/PA6 blends more significantly affect the morphologies than the distribution of the polymers, which should be taken into account in the preparation of UHMWPE/PA6 blends. It is worth noting thatfrom the enlarged view of Figures 2f and 4f, the orientation of UHMWPE and PA chains induced by extensional-shear coupled flow is more obvious than shear flow, which is in accordance withexperimental Raman spectra results in Reference [12]. results above, along with the consideration of the shear rates and extensional rates, we chose the shear rate and extensional rate as 1.0, in our simulations as presented below. In response to simple extensional flow, the micelle-like cluster is elongated along the flow direction. By applying simple shear flow and extensional-shear coupled flow, they enable chain-like network morphologies and uniform linear distribution, even though mass fraction varies. These results agree well with the systems of UHMWPE/PA6 blends at 50/50 mass fractions.     To shed more light onto the morphology and alignment behavior corresponding to extensionalshear coupled flow, we plot the mean square displacement (MSD) of PA6 as a function of simulation time for different coupled rates in Figure 7, as well as shear rates and extensional rates. MSD can be acknowledged as more effective parameter to characterize the average motion velocity of particles compared with the diffusion coefficient and it is calculated by Equation (14) [30,31]:

Extensional-shear coupled flow
where ( ) r t  and (0) r  refer to the coordinate of a particle at time t and 0, respectively. The slope of To shed more light onto the morphology and alignment behavior corresponding to extensional-shear coupled flow, we plot the mean square displacement (MSD) of PA6 as a function of  Figure 7, as well as shear rates and extensional rates. MSD can be acknowledged as more effective parameter to characterize the average motion velocity of particles compared with the diffusion coefficient and it is calculated by Equation (14) [30,31]: where → r (t) and → r (0) refer to the coordinate of a particle at time t and 0, respectively. The slope of the MSD-time curve represents the intensity of activity in polymer chains (or beads). We noticed that the slopes of MSD curve in Figure 7a exhibit slight increases upon increasing coupled rates because increasing coupled rates can enhance the random motion velocity of polymer chains, consequently raising the diffusivity of polymer chains, which provides clear proof that the alignment in UHMWPE/PA6 blends are significantly impacted by extensional-shear coupled rates, similarly in the case of shear flow and extensional flow (see Figure 7b,c). Figure 8 displays the relative concentration distribution of PA6 along Z direction. From the number and magnitude of the peak values, the distribution behavior can be evaluated.The fewer peaks and smaller the peaks values, the more uniform the distribution presents in the blends. An almost similar trend of the number and magnitude of peak values in relative concentration distribution curves is observed by comparingextensional-shear coupled flow, simple shear flow and simple extensional flow. This indicates that flow variations yield less effect on the distribution behaviors, in agreement with previous DPD results.

Conclusions
In this contribution, the Souza-Martins method was adopted in DPD simulation for predicting the flow behaviors in UHMWPE/PA6 blends, in which the effects of simple shear flow, simple extensional flow and extensional-shear coupled flow by varying shear rates, extensional rates and extensional-shear coupled rates, respectively, on the morphologies and distributionwere investigated for the first time. It was found that the flow variations and mass fractionsyielded lesser effects on the distribution behaviors. Nevertheless, the orientation behavior induced by extensional-shear coupled flow was more obvious than shear flow, in accordance with experimental Raman spectra results.The morphologies were more significantly affected by extensional-shear coupled rates in the case of extensional-shear coupled flow, as well as shear rates in the case of simple shear flow and extensional rates in the case of simple extensional flow. Within the considered range of extensional-shearcoupled rates, the morphologies evolved from micelle-like structure to chain-like structure as the coupled rates increased from 0.01 to 2.0, showing a linear distribution along the flow direction. It was evident by the increase of MSD value upon increasing rates and an almost similar trend of relative concentration distribution while subjecting the system to the same rates in relation to different flow.