Molecular dynamics simulation of injection of polyethylene fluid in a variable cross-section nano-channel

Numerical simulation of injection of polyethylene fluid in a variable cross-section nano-channel was carried out using the molecular dynamics method. The effects of the nano-channel cross-section and the external force on the rheological behavior and structural properties of the polyethylene fluid were investigated. It was found that an absorbed layer appeared near the wall and the thickness of the absorbed layer increased with increasing cone angle of the nano-channel. The injection distance of the polyethylene fluid decreased with increasing cone angle and decreasing external force. In the nano-channel with cone angle 45 ° , polyethylene particles uniformly filled the whole channel and were stretched along the flow direction. Uniaxial stretching of particles was enhanced when the external force was strengthened, which facilitates injection of the polyethylene fluid into the nano-channel.

With the development of nano-manufacturing, the technologies of precise distribution and delivery have been applied to a multitude of areas such as medical and biochemical engineering. Nano-injection systems, e.g. the nano-injector with an electromagnetic micro pump [1], nano-injection for capillary gas chromatography [2] and injection of ceramic bodies in nano-channels [3], have received much attention [4,5]. Injection molding of polymers is also utilized in many areas [6,7]. A polymer fluid exhibits unusual viscoelastic behavior in a nano-channel because of the topological constraints [8,9], so the process of polymer injection should be precisely controlled. The variation of any operating condition, such as temperature, viscosity, injection speed and pressure, can change the fluid properties of a polymer and destroy the injection process. Thus the complex *Corresponding author (email: lqzx@cqu.edu.cn) behavior of polymer flow in nano-channels should be fully understood.
Molecular dynamics simulations avoid the difficulties of experimental investigation, and have become an effective tool in analyzing the rheological behavior of polymer fluids. Several researchers [10][11][12] have employed non-equilibrium molecular dynamics simulations to understand the properties of polymer fluids in Couette shear flow. One can distinguish from literatures three models, namely the dumbbell model [13], the freely jointed chain model [14] and the Rouse model [15], to simplify the structures of the polymer chains. In the dumbbell model two beads that represent the molecular segments of several monomers are connected by a Hookean spring, which denotes the entropic effects on the end-to-end vector of the polymer. The freely jointed chain model depicts the polymer as a chain consisting of segments; the orientation of a segment is totally independent of the other segments in the chain. In the Rouse model, a Gaussian chain is considered to be physically equivalent to a string of beads connected by harmonic springs; each bead is regarded as a Brownian particle. The results obtained using these models are sometimes not in agreement with experimental observations, and assigning values to the representative parameters of the polymer is crucial to the accuracy of molecular dynamics simulations. Yashiro et al. recently reported characteristic parameters for polyethylene [16], which has been used as a typical polymer in many fields. In the present study the molecular dynamics method coupled with the unit atom model [17] was developed to simulate injection of polyethylene chains in a nano-channel with variable cross-section, and to analyze the impacts of channel structure and external force.

Simulation method and model
Polyethylene is formed by polymerization of ethene (CH 2 =CH 2 ), i.e. it consists of repeated methylene (-CH 2 -) modules. According to the unit atom model, polyethylene can be considered as a system consisting of many particles in which each particle represents a methylene unit, as shown in Figure 1 where the bond length, r, represents the distance between two particles, bond angle, θ, the angle between three particles, and torsion angle, φ, the angle between two planes consisting of four particles.
The potential energy of the system, E tot , can be expressed as the sum of the chemical bonding potential (bond stretch, bending and torsion) and the non-bonding potential (van der Waals) between the united atoms: where E bs (r) is the 2-body potential for bond length r, E bs (θ) the 3-body potential for bond angle θ, E to (φ) the 4-body potential for torsion angle φ, and ( ) vw E r another 2-body potential for non-bonded distance r within the cutoff radius of 8 nm. The potential functions adopted here were expressed as follows with the parameters listed in Table 1 [16]: where r 0 and θ 0 represent the equilibrium bond length and bond angle, respectively. The torsion potential was expressed by a cosine polynomial that has a stable point at φ=180° and two metastable points at φ=±67.5°. Τhe van der Waals potential was the 12-6 Lennard-Jones type with equilibrium distance about 4.5 nm. Figure 2 shows the shape and dimension of the nanochannel into which the polyethylene was injected. The channel consisted of a cubic chamber with size 300 Å×300 Å×300 Å connected to a duct with a square cross section of 100 Å×100 Å through a frustum of a rectangular pyramid. The cubic chamber was defined as zone A, and the other parts as zone B. The channel length was assumed to be constant, namely 800 Å, while the length of the rectangular frustum was considered as a variable, as shown in Figure 3, so as to investigate the impacts of channel structure. Thus the corresponding cone angle of the frustum was assigned the values α=45°, 71.6°, 78.7° and 90° in the simulations.
Initially 10000 polyethylene chains were randomly distributed in zone A, and each chain had 30 particles, corresponding to a density of 0.2588 g cm -3 . The system was relaxed by a molecular simulation for 20000 fs to reach an C (kJ mol -1 nm -6 ) 6.907×10 9  equilibrium state prior to injection. The injection process was conducted by moving the left movable wall in zone A along the x-direction with velocity 100 m s -1 . During this period the pressure gradient was simulated by imposing an external force on each particle in the system so as to keep the polymer chains homogeneous in the longitudinal direction [18]. The channel walls were considered as continuum surfaces that had no interactions with the fluid particles; the system temperature was kept at 450 K by the velocity scaling method; and the Verlet scheme was adopted for the numerical integration with a time step of 1 fs. The radius of gyration is an important parameter for characterizing the rheological properties of polymers. The mean square radius of gyration in the x-direction, which is a measure of the elongation and constriction of the polyethylene in the x-direction, was calculated as follows: where r pix is the x-position of the ith particle in the pth chain (i=1,···,N, and p=1,···,M) and r pgx =Σ i r pix /N represents the position of the center of mass of the pth chain, M is the number of chains, and N is the number of particles. The mean square radius of gyration in the y-and z-directions was obtained using a method similar to eq. (6).

Stability of polyethylene chains
Prior to injection, the polyethylene chains in zone A were relaxed until the equilibrium state was reached. Figure 4 shows the dependence of the particle number distributions on bond length, bond angle and torsion angle after 20000 fs relaxation. The particle number showed sharp peaks at the equilibrium bond length r 0 =1.533 Å and equilibrium bond angle θ 0 =113.3°, and two peaks at torsion angles φ=67.5 and 179°. This indicates that most of the polymer chains reached the equilibrium state. Because of the constraints between the chains, some particles with high energy were also observed in the vicinity of the equilibrium state, and this prevented the dihedral angle from reaching absolute stability, i.e. φ=180°. A visualization of the polyethylene chains in zone A after relaxation is presented in Figure 5, where some very small, uniformly distributed void spaces can be observed (the gray dots represent the methylene particles).

Rheological properties of polyethylene chains during injection
After the relaxation, the left wall of the cubic chamber was pushed to move with a velocity of 100 m s -1 , and an external  force was applied to each particle in the system to push them towards zone B. The injection process was stopped when the wall reached the interface between zones A and B because of the reduced cross section. Thus, the simulation time for the injection was 300000 fs.
(i) Effects of channel structure. Figure 6 shows the distributions of polyethylene chains in different nano-channels after injection when the external force was f=20×10 -23 Å -1 . With the help of the external force, all of the polymer chains entered zone B for all four nano-channels, and the flow distance increased with decreasing cone angle of the channel.
The simulation results reveal that the largest position that the polyethylene fluid reached increased from 556.56 to 576.23, 581.37 and 587.78 Å as the cone angle of the channel decreased from 90° to 78.7°, 71.6° and 45°, respectively. The reason for that trend is that the volume of the square duct is obviously large at large cone angle, hence decreasing the flow distance and particle density, although the resistance to liquid flow decreases with increasing cone angle. Figure 6 also indicates that after injection the void spaces in the polymer chains disappeared for the channel with cone angle 45°, whereas there were numerous heterogeneous void spaces in the straight nano-channel.
To further understand the distribution and movement of the polymer chains in the different channels, the nanochannel was divided into 30 slices from the wall to the central axis so as to obtain profiles of average density and average velocity of the methylene particles along the y-and z-directions using a statistical method. Figure 7 shows the density profiles of the polyethylene chains along the different directions. In all the cases the average chain density was much lower at the y-and z-walls but increased rapidly to a peak value as the slice number increased. The peak value and its position relative to the wall decreased with decreasing cone angle of the channel. The region between the wall and the position corresponding to the peak value is referred to as the absorbed layer. For the four channels that were investigated, the thinnest adsorbed layer was found in the channel with α=45°, the average density therein was the smallest, and the polyethylene chains were distributed uniformly in a high   proportion of the channel cross-section. Furthermore, the density distributions of the methylene particles along the yand z-directions presented in Figure 8 show that similar to the chain density distribution, while in all cases the largest particle density was near the wall, the channel with cone angle 45° showed the largest particle density at the wall. The particle density decreased with increasing distance from the wall in all cases, and it achieved a uniform distribution only in the centre of the channel with α=45°. With increasing cone angle, the thickness of the absorbed layer as well as the particle number therein increased, leading to reduction in the density of particles in the bulk regions. The particle densities in the bulk regions for the channels with α=71.6 and 78.7° were smaller than in the case α=45°. For the channel with α=90°, the absorbed layer was the thickest, and a serious fluctuation of the particle density and nonuniform particle distribution occurred in the bulk region. Figures 7 and 8 show that the methylene particles uniformly filled the whole channel with α=45° without any void space; while in the other cases the particle density was lower in the bulk region and exhibited larger fluctuations, and some void spaces were present.
The average velocity profiles of the methylene particles along the y-and z-directions are displayed in Figure 9. Because of the resistance from the tapered channels with α<90°, the average particle velocity was smallest at the wall and gradually increased towards the center of the channels. The average velocities in the different channels became closer as the slice moved from the wall towards the center, and approached a constant value. It is inferred from these results that the cone angle and the length of the pyramid frustum channel have a significant influence on the particle velocity profile. Small cone angle and long frustum channel produce significant resistance to fluid flow, leading to decrease in the flow velocity. Among the simulation cases, the lowest particle velocity was found in the channel with α=71.6° because of the combined effects of these two factors.
The chain configuration has an obvious influence on the polymer flow, thus the transformation of the chain configuration can be used to illustrate the behavior of the polyethylene fluid during the injection process. Figure 10 shows the changes of the configuration of a polyethylene chain that was near the wall during the injection process, in the four nano-channels. Initially, the chain curled up near the wall in zone A (Figure 10(a)). The configurations of the same polyethylene chain in the four channels at the end of the injection process are shown in Figure 10(b)-(e); it is apparent that the polymer chain was stretched to different extents in the four cases. The variation of the structure of another polyethylene chain located at the channel center is presented in Figure 11. The initial state of the chain with "I"  shape is shown in Figure 11(a), while the structures of the chain in the different channels at the end of the injection process are displayed in Figure 11 [19]. The distributions of the polyethylene chains at 320000 fs under different external forces, given in Figure 13, show that the polyethylene chains entered zone B pushed by the external forces, and the flow distance increased for strong force (the increment is more pronounced at stronger external force). The simulation results indicate that the initial  void spaces vanished after injection, and the polyethylene chain flow ceased at x=462.8, 473.9 and 587.8 Å for f 1 , f 2 and f 3 , respectively. Figure 14 presents the velocity distributions of the methylene particles along the y-and z-directions for different external forces. In all cases, the velocity had a minimum value at the wall and became larger towards the center line. For larger external forces the average velocity of methylene particles decreased.
The configuration changes of the polyethylene chains near the wall and at the channel center are shown in Figures  15 and 16, respectively, for different external forces after the injection process. Not surprisingly, the polyethylene chains stretched along the flow direction under the effects of external forces. Under external force f 3 , in particular, most of the polyethylene chains in the channel were parallel to the x-direction, which is beneficial to the polyethylene chains entering the small channel, and hence facilitates the injection process. The external forces not only pushed the polyethylene chains forward but also changed their structures and hence their resistance to flow. For example, the stretch in the y-and z-directions increased the flow resistance and reduced the occupation in the flow direction, whereas the converse results were obtained for stretch along the flow direction. This helps to explain why the simulation results in Figures 13-16 reveal that the largest external force (f 3 ) led to the smallest average velocity but the largest flow distance. Figure 17 shows the evolution of the radii of gyration along the x-, y-and z-directions under different external forces. Similar behavior of the radii of gyration was observed under external forces f 1 and f 2 , for which 2 gx R decreased at first because of the resistance from the pyramid frustum, then increased when the polyethylene chains entered the square duct, to recover the initial value. Opposite variations in 2 gy R and 2 gz R were found for the two cases. By contrast, the changes of the radii of gyration for case f 3 were not significant except for the increase in 2 gx R after 170000 fs and its final decrease. The simulation results presented in Figure 17 reveal that external force f 3 led to stretch of polyethylene chains only along the flow direction and thereby facilitated polyethylene injection.

Conclusions
Molecular dynamics simulations were performed to investigate injection flow of polyethylene in nano-channels with variable cross-sections. Impacts of channel structure (cone angle) and external force were explored. The following conclusions were drawn: (1) adsorption occurs near the channel wall, and the thickness of the adsorbed layer increases as the cone angle increases; (2) polyethylene chains are uniformly distributed in the channel with cone angle 45°, and stretch along the flow direction is observed during injection, thereby reducing flow resistance and increasing occupation in the flow direction; (3) the flow distance increases with decreasing cone angle and increasing external force; (4) applying an external force of 20×10 -23 J Å -1 in the channel with α=45° facilitates accomplishment of the injection process.