Numerical Investigation of the Effect of Viscoelasticity on Drop Retraction and the Evaluation of Interfacial Tension between Polymer Melts

*e purpose of this paper is twofold. *e first is to numerically investigate and reveal the effect of polymer viscoelasticity on the retraction of a deformed drop using the lattice Boltzmann (LB) method and polymer kinetic theory. More importantly, the second is to propose a novel method to evaluate the interfacial tension between polymer melts based on the numerical study. Compared with the conventional deformed drop retraction method (DDRM), the present method is designed to greatly reduce the impact of polymer viscoelasticity onmeasuring interfacial tension. To verify, the interfacial tension betweenmolten PP and POE is evaluated using the proposed method and obviously closer result to the true value is shown.


Introduction
By blending two or more polymers, it is an inexpensive and convenient way of obtaining new excellent materials that are complementary in performance. In recent years, polymer blends have been increasingly ubiquitous in the plastics industry, i.e., semiconductors [1], nanofiltration [2], fibers [3], and printing [4], and thus have become a very active area of research in materials science. e final properties of the polymer blends depend on not only that of each polymer component but also its internal microstructure, which is the result of the processing conditions as well as the interfacial tension [5,6]. erefore, it is of nontrivial science and engineering significance to measure the interfacial tension between polymers accurately for predicting and controlling the structure and properties of polymer blends.
Although the measurement technology for polymer has developed at increasingly faster pace in recent years [7][8][9][10][11][12], it still has been a challenging problem to obtain the interfacial tension experimentally due to the intrinsic viscoelasticity of polymers [13][14][15]. ese measuring technologies can be divided into two categories: the equilibrium methods and the dynamic methods [14]. e equilibrium methods, including the pendant drop method [16][17][18], the sessile drop method [5,6,19], and the spinning drop method [20][21][22], establish the equation of the drop shape in a mechanically balanced state and the extrapolate the interfacial tension between the polymers. Although they are applicable for both purely viscous and viscoelastic fluids with very good accuracy, it requires the matrix transparency, comparatively low viscosity of the polymer, and long time to reach equilibrium, with increased risk of thermal degeneration of the polymers. On the contrary, the dynamic methods, including the breaking thread method [23][24][25], imbedded fiber method [26][27][28][29], and deformed drop retraction method [30][31][32][33][34], obtain the interfacial through the morphology evolution equation. For example, the deformed drop retraction method establishes an equation that describes the shape variation of a deformed drop during retraction. e dynamic methods are comparatively more straightforward and fit for very viscous polymers. Nonetheless, the underlying theory of the existing dynamic methods is only valid for purely viscous fluids. If applied to polymers, the dynamic method will inevitably produce an error that is difficult to estimate. Some research studies on characterization technology are reported to be an inspiration for this problem, but the effect is not clear yet [35]. erefore, prior to proposing an improved dynamic method for polymers, the effect of viscoelasticity on the drop dynamics, i.e., deformed drop retraction, should be understood first, and then, attempt is made to reduce the effect.
To the author's knowledge, many simulation research studies have focused on the dynamics of viscoelastic drop. Yu and Zhou integrated the Boussinesq-Scriven viscoelastic interfacial constitutive equation in the perturbation analysis on the flow field inside and outside the drop [ [40]. e common idea of these studies is the combination of the Navier-Stokes equations, interface tracking algorithm, and the viscoelasticity constitutive model of polymer, i.e., the Maxwell model [41], the Voigt-Kelvin model [42,43], the transient network model [44,45], the Oldroyd-B model [46], the upper convected Maxwell model [46], and the FENE-P model [46,47]. e above review highlights the complex dynamics of viscoelastic drops and the difficulties in modeling. However, most studies are restricted to low viscosity cases and fail to simulate high Debora number (viscoelasticity), due to the limitation of their numerical model and the constitutive model.
In this paper, the effect of viscoelasticity on drop retraction is studied by the LB method and polymer kinetic theory, and then, a novel method for measuring the interfacial tension between polymer melts is proposed. e LB method is a special version of the Boltzmann equation that describes evolution of particles interacting on fixed lattices [48]. It has a kinetic, microscopic origin with a characteristic scale between nanometer and millimeter. On the contrary, the LB method can be Chapman-Enskog expanded to recover the full time-dependent incompressible and compressible Navier-Stokes equations. e LB method models the fluid by an ensemble of particles, so the macroscopic properties such as density and velocity can be easily constructed once an LB solution is obtained. e advantages of the LB method include the ease in dealing with arbitrary geometries and intricate multiphase flows, while its intrinsic parallel algorithm can be solved efficiently and affordably in massively parallel computers. During the past few decades, the LB method has been proved to be a promising tool for the simulation of complex fluid flow [49], i.e., microfluidics [50][51][52], flow through porous media [53][54][55], capillary flow [56][57][58], non-Newtonian flow [59][60][61], and more importantly, the multiphase flow [62][63][64]. e automatic-phase separation mechanism and straightforward interface tracking technology are both attractive characteristics of the LB method for simulating drop morphology evolution and has been successfully applied in research [65][66][67][68]. Different from the constitutive model of polymer viscoelasticity, the polymer kinetic theory explains the origin of viscoelasticity in a microscopic way, as the net effect of the dynamics of a large collection of constituent molecules with internal degrees of freedom. Interestingly, the macroscopic constitutive models, such as the upper-convected Maxwell (UCM) model and Oldroyd-B, can be derived based on the kinetic theory [46]. For the microscopic nature of both, Onishi et al. for the first time integrated the LB method and kinetic theory and succeeded in modeling the polymeric fluids [69,70]. Osmanlic and Körner applied the LB-kinetic theory framework to the simulation of Oldroyd-B fluids [71]. is paper combines the LB method and polymer kinetic theory to simulate the retraction of the deformed drop and analyze the effect of viscoelasticity. e rest of the paper is organized as follows. In Section 2, the retraction process of a deformed drop is simulated using a coupled model of the pseudo-potential multiphase LB scheme and a dumb-bell model of polymer molecule. e effect of viscoelasticity on drop dynamics is analyzed and a novel evaluation method of the interfacial tension is devised. en, in Section 3, the interfacial tension between molten polypropylene (PP) and polyolefin elastomer (POE) is measured by the proposed method and the original deformed drop retraction method, respectively. rough comparison the better accuracy of the proposed method is verified. Finally, a conclusion is drawn in Section 4.

Numerical Investigation
In essence, the LB method is a pseudoparticle method, which streams and collides on the lattice sites over discretized time and space. By tracking the evolution of the distribution function of the pseudoparticles, the fluid flow and drop morphology evolution is captured. On the contrary, the kinetic theory of polymer chains is rewritten in the lattice scheme. So the LB method-based computational framework consists of three ingredients: the pseudopotential multiphase model [72] for drop morphology evolution, the FENE dumbbell model for the polymer viscoelasticity [70,71], and an appropriate coupling strategy. e outline of the threedimensional numerical method is briefed in this section, while the details are presented in Appendix. Using the established model, several simulation cases of the deformed drop retraction is implemented and the effect of the viscoelasticity on the retraction process is analyzed. Finally, a novel evaluation method of interfacial tension is proposed.

2.1.
Modeling. Among all the multiphase models of the LB method, the pseudopotential approach first proposed by Shan and Chen [72] and improved by Shan and Doolen [73] is the most widely used by virtue of its simplicity and versatility and thus employed in this paper. e distribution functions of the two components, f d,i (x, t) and f m,i (x, t), which interact through pseudopotential force, describe the multiphase flow and result in the drop morphology evolution. e governing equation of the distribution functions is elaborated in Appendix, as equation (A.2). e term F i ′ denotes the interaction between different components and gives rise to spontaneous phase separation (shown in Figure 1).
is automatic mechanism of phase separation is an attractive feature of the pseudopotential model because the two-phase interface is no longer a mathematic boundary, but a postprocessing variable that can be characterized by monitoring the density change of the fluid components so that any specific interface tracking or interface capturing technology is not required as in conventional CFD methods.
From the microscopic point of view, the multiphase phenomenon and polymer viscoelasticity are a result of the intermolecular forces and the stretch and relaxation of polymer chain [74], and thus, it is feasible to model the macroscopic dynamics of viscoelastic polymer blends through the mesoscale polymer kinetics. e viscoelasticity of polymer is modeled based on the polymer kinetics theory; herein, the finitely extensible nonlinear elastic (FENE) dumbbell model is chosen, where two beads are connected by a string, as illustrated in Figure 2.
e conformation of the dumbbell can be defined by the tensor q that contains the length and orientation of polymer chain. Similar to the idea of the LB method, the function ψ(q, x, t) is defined to indicate the probability distribution of given polymer dumbbell conformation q at position x and time t. By solving the governing equation of the evolution of ψ(q, x, t), equation (A.15), namely, the Fokker-Planck equation [75], the conformations of the polymer dumbbells are updated all over the space with time. en, the viscoelastic stress is calculated according to the conformation according to equation (A.20). e evolution of the drop morphology and the polymer chain conformation is simulated simultaneously and separately and are coupled through rewriting the equilibrium function f eq ((A.24)). at is, the elastic stress originating from the dumbbell conformation act as a forcing term in Guo's scheme [76].

Simulation.
Before implementing the simulation, it is noteworthy that the simulation of the drop retraction process is based on two fundamental assumptions about the initial condition and the drop shape. Firstly, in terms of initial conditions, there is no residual stress in the drop/ matrix sample before the drop is set to retract. is is because in common practice of experiments, the drop is always heated after preparation in order to eliminate the residual stress. Secondly, the drop maintains an axisymmetric ellipsoid during the retraction process with two minor axes B and W of equal length, as illustrated in Figure 3. is seemingly ideal situation can be achieved as long as the sample is well fabricated and the density difference between the drop and matrix is negligible. e primary aim of the simulation is to study the effect of polymer viscoelasticity on retraction of the deformed drop. For this purpose, all the physical parameters but viscoelasticity strength of the polymers are equal in each simulation case so as to exclude the influence of other factors. e viscosities of both the drop and matrix components, μ d and μ m , are set to be 1000 Pa s, the density, ρ d and ρ m , 1000 kg/m 3 , and the interfacial tension between the drop and the matrix, σ, 0.116 mN/m. e pseudopotential parameter g σσ of LB is related to σ according to the Laplace law and decided through trials and error. e spring constant H is set to be 1.
ese parameters are kept constant during the drop retraction process. Since the drop retracts very slowly during the process, and the ambient temperature is maintained unchanged somehow, the influence of temperature could be neglected in the simulation. e strength of the polymer viscoelasticity is characterized by the Deborah number, as expressed below: which defines the ratio of relaxation time λ H and the characteristic time t c . It is obvious that the larger the De, the stronger the viscoelasticity, and a zero De indicates purely viscous polymer. e shape of the ellipsoid drop is defined by two parameters, equivalent radius R and the deformability D: Of all the simulation cases, the initial values of R and D are set to be 0.12 mm and 0.14, respectively. In addition, for the sake of characterization of the extent to which the ellipsoid drop retracts to sphere, a new parameter, the deformation recovery degree φ(t), is defined: In terms of the boundary conditions, periodic boundaries are applied on all the outer boundaries of the simulation domain.
ree numerical cases that only differ in viscoelasticity strength, De � 0, De � 1, and De � 5, are carried out. As expected, in all cases, the drop shape gradually changes from ellipsoid to sphere under the force of interfacial tension. Take the purely viscous case (De � 0) as an instance, and the shape evolution of the drop is shown in Figure 4. e major and minor axis length of the ellipsoid drop is calculated according to the component distribution of the drop by the eigenvalue method [77]. In all cases, the drop gradually retracts from an ellipsoid to a sphere, while φ continuously decreases from 0 to − ∞ overtime, as plotted in Figure 5.

Advances in Polymer Technology
From the above simulation results, it can be seen that when the polymer is purely viscous (De � 0), the deformation recovery degree φ is in good linearity with the time, consistent with the small deformation theory (equation (4)). It means that, in this condition, the original DDRM is valid: When the drop is viscoelastic (De � 1 and De � 5), the larger the De, the more curved the φ-t results. With the purely viscous case as a reference, the retraction process of the viscoelastic case can be divided into two distinct stages. In the first stage, t < 3000 s, the drop retracts at a comparatively faster rate, and φ decreases rapidly, but when the time elapses after about 3000 s∼3500 s (1.856∼2.165 t c ), the retraction rate of the drop gradually slows down. is is because the viscoelasticity of the polymer has a certain inhibitory effect on the retraction process. At the initial moment, the polymer chains do not elastically deform, so this resistance is relatively small. As the drop continues to retract, the elastic stress keeps growing due to the extension and compression of the dumbbells. In the second stage, the stress has grown to a considerable level so that its inhibitory effect obviously reduces the retraction rate of the drop. Using the data of different stages to implement parameter fitting with equation (4), the values of the obtained interfacial tension are as shown in Figure 6.
In Figure 6, label "t < 3500" indicates that the interface tension is obtained by fitting with data from before 3500 s "t > 3500," after 3500 s, and "whole stage" all the data. It should be noted that when De � 0, the DDRM is correct and the obtained value of the interfacial tension represents the true two-phase interfacial tension set in the simulation. When De ≠ 0, the viscoelasticity of the matrix affects the shape of the curve and the interfacial tension obtained using different stages of the simulation data has apparent discrepancy with the true value. Specifically, using the data of t < 3500 produces the fitted interfacial tension much closer to the true one, and the gap is nearly insignificant.
From the above simulations and analysis, it is obvious that, in case of viscoelastic polymers, the original DDRM faces inherent error caused by polymer viscoelasticity which will be augmented with the strengthening of viscoelasticity, whereas at the initial retraction stage, its effect is very weak in comparison to that of the second stage. If only picking the data of the first stage of the drop retraction for fitting with equation (4), then the deviation between the measurement and the true value can be greatly reduced so that a credible result can be obtained.   Herein, a novel method to evaluate the interfacial tension is proposed according to the deformed drop retraction. Other than fitting the data of the whole process with equation (4) during which the drop retracts from an ellipsoid to a sphere, only the data from the onset of retraction t 0 to the time φ(t u ) � φ u are utilized. e duration between t 0 and t u is the first stage of the drop retraction, and empirically, the value φ u is more appropriate to set − 1.5.

Advances in Polymer Technology
It should be noted that this treatment also agrees with the dumbbell model. e dumbbell model of polymer molecules suggests that only the terminal polymer relaxation dynamics is relevant. It is just the elasticity that makes the original DDRM invalid. erefore, the neglect of polymer dynamics at intermediate time scales and thereafter (which actually gives rise to the elastic component of viscoelasticity) will alleviate the effect of elasticity on drop retraction.

Materials and Set-Up.
Polypropylene (PP) is a widely produced, versatile, commodity polymer with a series of desirable properties, and it is often blended with polyolefin elastomer (POE) (see Figure 7) to improve its toughness and impact strength, alleviating notch sensitivity and low-temperature brittleness. In this section, the interfacial tension between PP and POE melt is measured using the proposed method in an optical microscope equipped with a heating stage. PP, with the trademark T36 F, supplied by Chevron (Zhangjiagang) Chemical Co., Ltd, composes the drop. POE, with the trademark 8150, supplied by Dow DuPont (USA) Co., Ltd, composes the matrix. POE is a typical kind of viscoelastic polymer, while the viscoelasticity of PP is negligible compared to POE, so PP is considered purely viscous.
Before the experiment, necessary major properties of the two polymers are measured and shown in Table 1.
e experiment consists of three steps, that is, the preparation of experimental samples, getting PP and POE in full contact, and observing the drop retraction process sequentially.

Preparation of Experimental Samples. First, the PP fiber is drawn in molten state on a melt indexer instrument.
Since the melting temperature of the PP is about 170°C, the heating temperature of the apparatus is set to 230°C. By adding different weights to the melt indexer to apply different drawing forces, PP fibers of different diameters, 0.2 mm, 0.3 mm, and 0.4 mm, are obtained and then cut them into short fibers with the aspect ratio of roughly 3∼4. en, the POE sheet is molded out on a compression molding press. e temperature of the upper and lower boards of the machine is set to be 180°C to make two POE sheets with the thickness of about 0.5∼0.6 mm.
Finally, one piece of PP short fiber is embedded between two POE sheets to make a sample by stepwise heating and packing, and the parameters are shown in Table 2. rough the above treatment, the air bubbles trapped between the POE sheets are pressed out, and the sample is compact enough, the structure of which is similar to the hot dog bread, as illustrated in Figure 8.
On the contrary, in order to relieve the stress generated in the preparation of the sample, the sample is subsequently placed in a vacuum drier, annealed at 100°C for 3 hours, and then cooled in air.

Getting PP Short Fiber and POE Sheet in Full Contact.
After the above preparations, there may be still small amount of air remaining in the vicinity of the contact surface of the fiber and sheets. If not eliminated, bubbles will generate and aggregate when the sample melts, causing interference to the observation. In order to remove the air totally, the sample is put into the vacuum dryer again and heated at 200°C under vacuum condition for 10 min. In this process, the air largely escapes out of the sample when the POE sheets soften and warps the PP short fiber.

Experimental Observation and Data Processing.
e test sample is placed on a slide of the hot table, and the temperature is raised to 120°C at the rate of 40°C/min by program setting and kept for 10 min to further eliminate the residual stress in the sample. Next, the heating temperature of the sample is increased to 220°C, at which the PP and POE are melted, and the shape evolution of PP short fibers with diameters of 0.2 mm, 0.3 mm, and 0.4 mm, is observed, respectively, while a camera takes consecutive shots of the process.
When all the cases have finished, the software Scion image is used to measure the length of the major and minor axes of the ellipsoid drop, L and B, at each shot time and the radius R when the spherical drops are restored.

Results and Discussion
e retraction process of the 0.2 mm diameter PP fiber in POE sheets is shown in Figure 9. From the figure, it is clear that the retraction process of the drop can be divided into two phases: at the first phase, shape of the fiber changes from a rod to an ellipsoid and at the second phase, from the ellipsoid to a sphere.
Likewise, the PP short fiber with the diameter of 0.3 mm and 0.4 mm has the similar morphology evolution process.
From the recorded images, the length of the major and minor axes of the ellipsoid drop, L and B, at different times is measured by the software Scion image and then the variation in the deformation recovery degree φ(t) over time is calculated, as shown in Figure 10. A common first-fast-thenslow retraction process is observed, and from the simulation results, it could be inferred that even if the experiments were performed much longer, a second stage with slower retraction will probably not occur.
ere are twofold reasons for using a set of three PP short fiber of different diameters. Easy to see, measuring several times and averaging the results are good for improving accuracy. More important, it is reported that the influence of polymer viscoelasticity on drop dynamics decreases with the increase in the drop size [78,79], so trying experiments with different drop sizes from small to large will gradually approach the true value of the interfacial tension. 6 Advances in Polymer Technology In the following, it is demonstrated that, using the proposed evaluation method of interfacial tension, even with the comparative small-sized drop, a desirable accuracy is attained. e principal difference between the proposed method and the original DDRM lies in the way of deriving the interfacial tension from the shape evolution of drop retraction. eir measurement results are shown in Table 3. e label "DDRM" indicates that the interfacial tension is obtained by fitting with all the experimental data, while the label "I-DDRM" indicates that only the data spanning from the beginning of the drop retraction to the time when φ reduces to φ u is used to fit equation (4).
At the bottom of Table 3, the results fluctuation fluc(A) characterizes the impact of drop size on the measurement results: where A is the set of all the measured data, max(A) is the minimum value in A, and ave(A) is the average of A. From Table 3, it is noteworthy that using the conventional DDRM, the measurement result is apparently dependent on the drop size and the fluctuation is as high as 28.66%. e reason behind this phenomenon is the effect of polymer viscoelasticity. By contrast, when using the proposed evaluation method, the measurement results are very stable with whatever size drop and the fluctuation of the measurement results is strikingly reduced to 3.98%. On the contrary, the interfacial tension obtained by the conventional DDRM with the smallest drop (0.4 mm diameter) is thought to be the closest to the true value. is measurement value is approximately the one obtained by the proposed method in all the three cases of different drop sizes. e experimental comparison and analysis agree well with the revelation in Section 2 and validate the interfacial tension evaluation method proposed by this paper.
To supplement, the proposed method has another implication. ere is no need to try a beforehand unknown number of experiments with different drop sizes to dampen the effect of polymer viscoelasticity and approximate the true value step by step, as is done now in the conventional DDRM. By employing the novel evaluation method proposed in this paper, the value of the interfacial tension between viscoelastic polymers can be evaluated close to the true value reliably and once for all.
Finally, the effects of temperature and polymer molecular weights on the drop retraction process and interfacial tension should be mentioned. Kamal et al. [80] and Biresaw et al. [81] showed that higher molecular weight systems showed a weaker dependence of interfacial tension on temperature than lower molecular weight systems. So, the interfacial tension    Advances in Polymer Technology 7 between the same polymers of different molecular weights is not the same and needs measuring separately.

Conclusion
In this paper, the deformed drop retraction process is simulated taking account of the fluid viscoelasticity and its effect on the shape evolution of the drop is revealed. Compared with the case of purely viscous blend, the drop retraction process follows the "first fast then slow" law. en, a novel evaluation method of the interfacial tension is put forward that only takes advantage of the first stage of the experimental data. e evaluation method is successfully validated by measuring the interfacial tension between PP and POE. Compared to the original DDRM, the obtained interfacial tension is close to the true value, but the experimental times are reduced to only once. Due to the complexity of polymer viscoelasticity, the dumbbell chain model is a simplified model and has its own limitations. It cannot describe more complicated spatial configurations of the polymer molecule, let alone the entanglement phenomenon of polymer chain in the melt. erefore, to advance the measuring for interfacial tension between polymers, making more thorough insights into polymer dynamics is desired.

A. Pseudopotential Lattice Boltzmann Model
In this paper, three-dimensional LB method is implemented on the D3Q19 lattice, which is one of the most commonly used lattice types nowadays ( Figure 11).
As illustrated in Figure 11, in D3Q19, the velocity space is discretized into a set of 19 discrete velocities: where c � δx/δt is the lattice speed, δx is the lattice spacing, and δt is the time step. Using the multirelaxation time collision operator MRT, the governing equation of f i can be written as where Ω i is the collision operator and defined by f t is the density distribution function associated with velocity c i at position x and time t, and their moments satisfy where F represents the interaction force between components. F i ′ is the discrete forcing term accounting for the interaction force F: where I is the unit matrix, while F is defined by In the collision operator, M is the constant transformation matrix and S is a diagonal nonnegative relaxation time matrix that is decided by fluid properties.
For multicomponent flow, f eq σ,i is the corresponding equilibrium distribution of component σ and is determined by where w i is the weight associated with the velocity e i : u eq σ is the equilibrium macroscopic velocity in consideration of the component interaction: where u ′ satisfies the momentum equation without the external force: in which ρ σ and μ σ are the macroscopic density and velocity of component. F is the interaction between the particles of different components, when only considering the nearest neighbor, it can be written as

Advances in Polymer Technology
where g σσ denotes the interaction strength between components σ and σ. When g σσ is positive, the repulsive force is generated between the particles, and when the g σσ is negative, the attraction between the particles is generated. With the modified equilibrium distribution and pseudopotential forcing term, the Navier-Stokes equations can be recovered from the lattice Boltzmann equation in the limit of small Mach number by the Chapman-Enskog analysis and Taylor expansion: as well as the momentum equation of the blend flow: where υ is the kinematic viscosity of the blend: (A.14)

B. FENE Dumbbell Model
Likewise, the conformation distribution function ψ j (x, t) follows the discrete Fokker-Planck equation on D3Q19 lattice: As the polymer melt is thought isotropic at equilibrium, the equilibrium conformation distribution function ψ eq j is equal to w j : τ ψ is the single relaxation time of ψ j (x, t) towards ψ eq j · M j accounts for the elongation and rotation of the dumbbell: (A.17) Δψ j (x, t) is the convection of dumbbells between lattice sites due to the flow and can be written as Solving the equation (A.2) and equation (A.15) simultaneously, the conformations of the dumbbells can be updated during the flow simulation. en, the macroscopic elastic stress is linked with the dumbbell conformation through the Kramers formula: where F c j is the spring force acted on the dumbbell beads: As mentioned above, the flow field changes the conformation of the dumbbells; in turn, the conformation of the dumbbells contributes to the local stress and thus alters the flow field. In the modeling, they are coupled through the rewriting of f eq : It is noticeable that the modified expression reflects the effect of polymer viscoelasticity through the introduction of 9/2ω i (e i e i /c 2 − (1/3)I) : τ p /k B T m into f eq . After the modification, the second-order moment of the distribution function contains both viscous and elastic stress: Data Availability e experimental data used to support the findings of this study are included within the article. e software code data used to support the findings of this study are available on request to the Lin Deng via hust.smart@gmail.com.

Conflicts of Interest
e authors declare that they have no conflicts of interest.