Simulation of Fountain Flow based on Molecular Dynamics Method

In this paper, numerical simulations of the fountain effect on the flow front of polymer dilute solutions are carried out based on molecular dynamics principles and the FENE dumbbell model. The Euler method is used to solve the constitutive equations of the mechanical model for the simple shear flow field and the dislocation equations of the dumbbell molecules, and subsequently obtain the tracer flow lines and the dumbbell distribution of the flow front. The results are used to calculate the stress field, analyze the rheological evolution law, study the effect of temperature and shear rate and other parameters on the model and analyze the effectiveness of the FENE dumbbell molecular model for the fountain effect. The study shows that due to the fountain effect the dumbbell molecules stretch more with increasing shear rate, the polymer stress at the flow front will show a complex change of stress overshoot and then stabilize, and the dumbbell molecules are oriented along the fountain flow line.


Introduction
Most polymer solutions fall under the category of viscoelastic fluids.The crux of studying these viscoelastic fluids is to ascertain the microscopic characteristics inherent in the polymer.Unfortunately, the conventional methods of continuous medium mechanics fall short in providing insights into the microscopic movements of polymer molecules.To bridge this gap, we propose the use of polymer molecular mechanical models.These models have the potential to investigate the rheological properties of high polymer systems, thereby exposing the intrinsic relationship between the molecules' microstructure and morphology, and their macroscopic rheological performance.In the fluid filling process, two distinctive regions are observable in the flow of high polymer solutions: the mainstream region and the leading edge region.Each of these regions exhibits significantly different polymer flow behavior.The mainstream region nearly shows unidimensional flow, whereas the flow in the leading edge region is much more complex.Owing to the influence of shear flow, high polymer solutions align along the flow direction, continuously compelling the fluid to move from the center towards the mold wall, thus displaying a characteristic "fountain flow" pattern.This fountain flow significantly influences the microstructure, as well as the mechanical and thermophysical properties of the end products.Traditional methodologies such as control volume, Lagrange-Euler methods, and Level Set methods fall short in simulating the fountain flow phenomena occurring at the leading edge; this necessitates a specific description for the leading edge fluid.Rose [1] was the pioneer in defining and simulating the interfacial evolution within a capillary by leveraging the fluid displacement method.E. Mitsoulis et al. [2] conducted a parameterized investigation on Newtonian fluids within planar and axisymmetric geometries, delivering a comprehensive analysis of the influence of various fluid mechanic parameters on fountain flow, thereby significantly enhancing the standard of fountain flow research.Wen Yan et al. [3] constructed a three-dimensional model describing the fountain effect at the leading edge during the mold-filling process of short fiber-reinforced polymer melts.This model dynamically tracked fiber motion and orientation, effectively reflecting the three-dimensional fountain effect at the leading edge of the melt.Bechara et al. [4] simulated the impact of fountain flow on fiber orientation and distribution by combining finite volume models and mechanical models.This innovative approach extended the conventional methodology by integrating mechanical methods.Papathanasiou et al. [5] analyzed the influence of fountain flow occurring near the advancing free surface on fiber orientation during the filling flow, providing a systematic overview and discussion of experimental and computational results related to this topic.Liu et al. [6], leveraging the Rolie-Poly model, explored the effects of inlet velocity, polymer molecular weight, and temperature on the evolution of the leadingedge fountain effect of polystyrene melts, effectively simulating the viscoelastic behavior of polystyrene.Borzenko et al. [7][8][9], through the Ostwald-de Waele rheological equation, examined the impact of pseudo-plastic and dilatational characteristics of the liquid medium on the fountain flow motion of non-Newtonian fluids during the tube filling process, thereby explaining the elimination of contact line singularity in traditional methods.Smit et al. [10] predicted the instability of fountain flow using the PTT-XPP model, compared the elastic effects in the flow, and validated the impact of reduced critical Weissenberg numbers on fountain flow characteristics.Jong et al. [11] employed a visualization mold design to study the influence of gas back pressure and mold temperature on fountain flow, thus advancing visualization research concerning fountain flow.Velmaat et al. [12] simulated the fountain effect using a meshless method.Lastly, Ren et al. [13] implemented an improved SPH method to track particle motion within the mold cavity, verifying the fountain flow effect at the melt's leading edge.In comparison to traditional methodologies for simulating fountain flow, which may involve certain oversights or approximations, molecular dynamics methods provide an accurate simulation of molecular movement and interactions.These methods also operate on a smaller time scale, thus offering more detailed data.The FENE (Finitely Extendable Nonlinear Elastic) dumbbell molecular model is a molecular dynamics method used for simulating polymer molecular motion.Established on the foundations of the dumbbell molecular model, it introduces additional elastic energy to simulate the interactions among polymer molecules.This model postulates that intramolecular forces can be depicted by an array of spring forces or constraint conditions [14], offering high efficiency, accuracy, and adaptability.In this paper, a FENE dumbbell molecular model, based on the principles of molecular Brownian dynamics and tailored to the peculiarities of the rheological process of polymer solutions, is proposed.This model enables simulation analysis of the fountain effect occurring during the mold flow process in polymer dilute solutions.By harnessing the resolved configurational and positional information, we compute the variance and evolutionary trends of stress at the leading edge and assess the impact of various parameters on the model.The model's validity is affirmed through an analysis of fountain flow lines.

Theoretical Foundation
The Finitely Extensible Nonlinear Elastic (FENE) dumbbell model, composed of two small beads connected by a spring of a specified length as depicted in Figure 1, is a molecular dynamics model utilized to simulate the interaction between polymer molecular chains.This model is ideally suited for numerical simulation studies of polymer dilute solutions.In this context, the beads simulate interacting particles within the solution, while the massless spring exemplifies the constraints of the molecular chain.The position vectors of the beads are represented by   ,   , whereas Q signifies the dumbbell's connecting vector between the beads, and c r stands for the center of mass of the dumbbell molecular chain.The spring constraint force between the rigid beads signifies the influence of entropy elasticity and is directly associated with the configuration vector of the dumbbell.Within the FENE dumbbell molecular model, according to Warner's law of elasticity, the effective constraint force s i F of the finitely extendable nonlinear spring is denoted as [15] 1 ( / ) Wherein (2) The four forces represented in the above equation are fluid resistance (including surface tension), spring constraint force, Brownian force, and shear stress.As the bead travels through the solvent, liquid resistance is a frictional force generated due to the relative speed difference between the moving molecular chain in the flow field and the solvent molecules, in addition to the surface tension produced by the free surface at the flow's leading edge.This embodies the Stokes drag force exerted on the rigid bead.

( ( ))
(3) here  is resistance coefficient of the fluid, ̇ is the speed of the rigid bead, and () i vr represents the speed of the polymer solvent at the position of the i-th bead.Within the shear flow field, a shear stress oriented towards the direction of the mold wall is generated, inducing radial movement towards the direction of the mold wall.The shear stress is expressed as: The Brownian force instigates the rapid collisions of the polymer solvent molecules with the polymer chain over a short period of time, fulfilling: In the above formula, B k denotes Boltzmann's constant, T stands for absolute temperature,  is the delta function, and signifies ensemble average, which can be represented by the last term in equation ( 6) within the time span dt.By substituting equations ( 1), ( 3), ( 4), and (6) into equation ( 2) and subsequent rearrangement, the ensuing equation can be deduced to govern the evolution of the bead's position: In the aforementioned equation, n denotes a random two-dimensional vector associated with the Wiener process.Each component of n adopts a Gaussian random number with a mean of zero and a variance of one.During computations, a random number within the range of [-1,1] can be utilized as a replacement.This approach not only conserves computation time but also does not substantially affect the convergence of the solution.
As well as the equation for the center of mass of a single dumbbell molecule: 2 + = 12 c rr r (9) By subtracting the position evolution equation for the rigid bead, we derive the configurational control equation for the dumbbell: Wherein the orientation angle of a single dumbbell molecule relative to the flow direction can be depicted as: arctan( ) In the aforementioned equation, , xy respectively denote the flow direction and the direction of the velocity gradient.The shear stress tensor of a polymer solution is constituted by contributions from both solvent molecules and polymer molecules.The polymer stress tensor can be expressed as: In the above equation, 0 n represents the number of polymer molecules per unit volume, and δ stands for the unit tensor.

Shape Calculation
For the Finitely Extensible Nonlinear Elastic (FENE) dumbbell molecular model, we can introduce a time constant  4   , a distance constant √/, and a dimensionless parameter  =    0 2 /   .
For the equation controlling the position evolution of the dumbbell, divide both sides of Equations ( 8) and ( 10) by the distance constant √/, and then multiply the right side of the equation by the time constant   =  4   , to complete the non-dimensionalization process.After non-dimensionalization, we get , ̅ 1 , ̅ 2 .For convenience, we still denote  as  and  ̅ as , as shown in Equations ( 13), (14) The analysis leads to the equation for the horizontal conformation of the dumbbell molecule as: The equation for the vertical conformation is as follows: The non-dimensional representation of material viscosity is as follows [14]: , 00 () 1 The initial dumbbell configuration of the dumbbell is randomly obtained within the range of [0, √/2].
The initial position of 1 r is randomly obtained at x = 0 in the direction of the mold wall along the center plane, then the position vector 2 r of the bead is obtained through calculation of 21 r = r + Q .In the first calculation, the shear stress is taken as F  = , after which the iterative operation is conducted according to formulas ( 17) and (18), and Q is substituted to calculate the configuration equation of the next step Q , thereby obtaining the configuration vectors and shape vectors at all times.

Material parameters and Calculation Scheme
In this study, the fountain effect at the fluid front in the simple shear flow field during the fluid filling process is considered.It is assumed that the dilute solution is incompressible within the area affected by the fountain effect, and the gravity and inertial forces of the solute molecules are neglected, as are thermal convection and conduction during the flow process.The Euler explicit difference method is used to simulate the flow field and flow front.The simulation parameters for the material calculation in the dumbbell molecular model are as shown in Table 1.
The velocity components of the shear flow field in the polymer dilute solution are set as   = (1 − /ℎ) ⋅ ̇ and 0 y v = , where h is the length from the center plane to the wall of the mold.Considering the influence of surface tension, the velocity components of the front part of the solution are set as   = (1/2 − /ℎ) ⋅ ̇ and 0 y v = , which serve as the boundary conditions at the flow front.Taking into account the symmetry of the flow field, only the upper half is selected for computation.Given that the equations cannot be directly solved to obtain an analytical solution, the Euler numerical method described earlier is employed.The initial dumbbell conformation information is first obtained, then equations ( 14) and ( 15) are solved using the initial dumbbell conformation to get the initial state information.Subsequently, equations ( 19) and ( 20) are solved.During each iteration, the conformation information of the dumbbell molecules is sorted.The dumbbell molecules at the forefront are calculated iteratively using the forefront flow field, while the remaining dumbbell molecules are calculated according to the original flow field velocity components.This process enables the acquisition of all connection vectors as well as bead conformations.
( ) The total number of dumbbell molecules in this study is n = 10000, and the dimensionless time step is set to 1e-3.When the computational time step is not particularly small, non-convergence situations may occur due to computational errors and the action of random forces, where the spring length exceeds the maximum extension length in the simulation process.This study adopts two treatments: firstly, increasing the maximum extension length of the dumbbell spring; secondly, using a method of sequentially reducing the time step [17] when overstretching occurs.The time step is reduced by a factor of N, and the previous state is recalculated in N steps.If overlength situations persist, N is increased, and computation continues from the previous state.Afterward, the original time step is used for calculation.This approach can both avoid persistently using shorter time steps, thus saving computational time, and control the solution of vectors within the convergence interval using shorter time steps.The computational flow is illustrated in Figure 3.

Figure 3. Flow chart of simulation calculation
At the inlet cross-section, information about all dumbbell molecules at the centroid of the crosssection is first obtained at each moment.Then, all nearby dumbbell molecules are iteratively acquired for ensemble averaging, which is used to calculate the polymer stress.At the same time, the shape changes and molecular chain fluctuations of individual dumbbells during motion are tracked.The stress change at the flow forefront is obtained by ensemble averaging of the dumbbells at the forefront at any time.

Results and Discussion
By tracking the shape equation of the dumbbell molecules and randomly selecting nine dumbbell molecules for tracking under the condition of   =1, we obtain the motion trajectory as shown in Figure 4, aligns macroscopically with the fountain streamline described by E.Mitsoulis [2].It's identical to the expected streamline of the fountain effect at the flow front in theory.As the dumbbell molecules move, the velocity in the main flow direction gradually decreases while gaining a radial velocity, slowly surging towards the wall.The motion trajectory of dumbbell molecules near the center plane more closely conforms to the idealized fountain effect, manifesting a more pronounced fountain surge.6 records the conformations of certain dumbbell molecules at the forefront at four distinct moments.With the motion of the dumbbell molecules, the forefront distribution gradually transitions from a slope similar to the flow field to a frontal curve closer to the free surface, progressively becoming stable.In the course of the flow, changes such as rotation and stretching occur in the dumbbell molecular chains.As time progresses, the orientation of the dumbbell molecular chains along the direction of motion becomes increasingly evident.This is in accordance with the simulation and experimental outcomes of Jong [11], as well as Borzenko's [8] simulations in a circular tube, corroborating the effectiveness of the FENE dumbbell model for simulating the flow front in polymers.
Nonetheless, there remains a discrepancy with the idealized semi-circular shape, necessitating further in-depth study.8 depicts the curve illustrating the variation of the polymer's viscosity coefficient.As can be seen from the figure, the viscosity of the polymer exhibits a declining trend with the increasing shear rate, indicating a shear-thinning behavior.Additionally, as the temperature escalates, the intermolecular forces decrease and Brownian motion becomes more pronounced, leading to a further decline in polymer viscosity.The calculated results for the polymer stress at the flow forefront under different shear rates are presented in Figure 9.When the shear rate is relatively high, the stress at the flow forefront of the polymer swiftly grows to its peak value before gradually reaching a steady state.An overshoot of stress under high shear rates is also observed, aligning with Min Zhiyu's [14] simulation results on the rheological properties of polymer molecules.

Conclusion
In this paper, we implemented a simulation methodology for the fountain flow at the forefront in dilute polymer solutions based on the Finitely Extensible Nonlinear Elastic (FENE) dumbbell model.Our approach successfully simulated the fountain effect at the forefront, and we recorded the IOP Publishing doi:10.1088/1742-6596/2706/1/01208911 configurational distribution of the dumbbells, achieving satisfactory fountain streamline results.These findings confirm the appropriateness of the FENE dumbbell model for simulations of forefront flow in dilute polymer solutions.Furthermore, we carried out a simulation analysis of the cross-sectional and forefront stresses during the flow process.The stress exhibited complex variations in line with the fluctuations in the region's molecular weight.Stress at the cross-section, influenced by the fountain effect at the forefront, showed initial fluctuations before reaching a state of stability.As the shear rate increased, a stress overshoot occurred at the forefront, with larger peaks and smaller peak intervals at higher shear rates.Dumbbells situated closer to the center plane were more significantly influenced by the fountain effect.Initially, forefront molecules oriented along the direction of main flow.With the emergence of the fountain effect, they gradually oriented towards the direction of the fountain streamlines, eventually reaching a stable configuration.The molecular chain extension gradually reached stability, and the influence of surface tension on the forefront flow was comparatively weak.
r denotes the connecting vector of the spring, s H represents the spring constant, and 0Q signifies the maximum extendable length of each spring.

Figure 1 .
Figure 1.Schematic diagram of the dumbbell molecular model The rheological behavior of polymer molecules is predominantly affected by various factors including viscous drag, spring constraints, Brownian motion, intermolecular interactions, excluded volume, and molecular chain entanglement.Ignoring the impacts of intermolecular forces, excluded volume, and molecular chain entanglement during the flow process, and while accounting for the 'fountain' effect of leading-edge flow, the force analysis of the dumbbell molecule necessitates the calculation of the influence of shear stress on fountain flow lines.The leading-edge flow is represented in Figure 2.

Figure 2 .
Figure 2. Schematic diagram of the frontier areaDisregarding the inertial effects and mass forces of the molecular chains, the governing equation for the dynamic chain can be derived from the equation stipulating that the sum of external forces exerted on the rigid beads equals zero[16]. 0

Figure 4 .
Figure 4. Dumbbell molecular tracer streamlines Figure 5 shows the simulated stress at x=0.05 and x=0.06.The stress initially increases over time, reaching a peak, then progressively shifts to a steady-state value.Some fluctuations occur before stabilization, due to changes in the number of molecules resulting from the motion of the dumbbell molecules.The maximum fluctuation of stress at x=0.05 surpasses that at x=0.06, but the steady-state value is smaller at x=0.05.Moreover, the degree of stress fluctuation is larger at x=0.05, and it also reaches a steady-state more quickly.

Figure 5 .
Figure 5. Cross-sectional stress fluctuation with time Figure6records the conformations of certain dumbbell molecules at the forefront at four distinct moments.With the motion of the dumbbell molecules, the forefront distribution gradually transitions from a slope similar to the flow field to a frontal curve closer to the free surface, progressively becoming stable.In the course of the flow, changes such as rotation and stretching occur in the dumbbell molecular chains.As time progresses, the orientation of the dumbbell molecular chains along the direction of motion becomes increasingly evident.This is in accordance with the simulation and experimental outcomes of Jong[11], as well as Borzenko's[8] simulations in a circular tube, corroborating the effectiveness of the FENE dumbbell model for simulating the flow front in polymers.

Figure 6 .
Frontside dumbbell molecular bit shape distribution Figure7shows the analysis of the time-dependent fluctuating extension of polymer molecular chains at the forefront under varying shear rates.The extension of molecules progressively lengthens and approaches a steady state as the dumbbell molecules traverse the flow field.When H  equals 4, the extension length of the dumbbell molecules tends to 0.65, while for H  equaling 1, it converges to approximately 0.4.The stable value of the dumbbell extension is larger under high shear rates, with more pronounced fluctuations observed during the transformation process.Additionally, a longer period is needed for the dumbbell to reach a steady state under lower shear rates.

Figure 7 .
Figure 7. Dumbbell molecular chain elongation fluctuates with time

Figure 8 .
Figure 8. Variation of shear viscosity with shear rate and temperature

Figure 9 .
Figure 9. Variation of shear stress with time

Table 1 .
Material parameter values Material parametersValues