The purpose of this paper is twofold. The 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 on measuring interfacial tension. To verify, the interfacial tension between molten PP and POE is evaluated using the proposed method and obviously closer result to the true value is shown.

1. 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. The 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]. Therefore, 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 [712], it still has been a challenging problem to obtain the interfacial tension experimentally due to the intrinsic viscoelasticity of polymers [1315]. These measuring technologies can be divided into two categories: the equilibrium methods and the dynamic methods [14]. The equilibrium methods, including the pendant drop method [1618], the sessile drop method [5, 6, 19], and the spinning drop method [2022], 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 [2325], imbedded fiber method [2629], and deformed drop retraction method [3034], 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. The 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].

Therefore, 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 [36]. Mukherjee et al. simulated the shape relaxation of a sheared drop in the case that either or both phases are Oldroyd-B fluids by the finite difference method and front tracking algorithm [37]. Yue et al. used the 2D diffuse-interface method to simulate retraction of a Oldroyd-B drop in quiescent matrix [38]. Yue et al. incorporated the Oldroyd-B viscoelasticity constitutive equation into the phase-field framework of Cahn–Hilliard equation with the adaptive meshing scheme and studied the interfacial dynamics of the drop dynamics [39]. Yu et al. extended Batchelor’s approach to viscoelastic fluids and calculated the creeping flow around the drop [40]. The 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]. The 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. The 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. The 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. The 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 [5052], flow through porous media [5355], capillary flow [5658], non-Newtonian flow [5961], and more importantly, the multiphase flow [6264]. The 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 [6568]. 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]. This paper combines the LB method and polymer kinetic theory to simulate the retraction of the deformed drop and analyze the effect of viscoelasticity.

The 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. The effect of viscoelasticity on drop dynamics is analyzed and a novel evaluation method of the interfacial tension is devised. Then, 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. Through comparison the better accuracy of the proposed method is verified. Finally, a conclusion is drawn in Section 4.

2. 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. The outline of the three-dimensional 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. The distribution functions of the two components, and , which interact through pseudopotential force, describe the multiphase flow and result in the drop morphology evolution. The governing equation of the distribution functions is elaborated in Appendix, as equation (A.2). The term denotes the interaction between different components and gives rise to spontaneous phase separation (shown in Figure 1).

This 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. The 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.

The conformation of the dumbbell can be defined by the tensor that contains the length and orientation of polymer chain. Similar to the idea of the LB method, the function is defined to indicate the probability distribution of given polymer dumbbell conformation at position and time . By solving the governing equation of the evolution of , equation (A.15), namely, the Fokker–Planck equation [75], the conformations of the polymer dumbbells are updated all over the space with time. Then, the viscoelastic stress is calculated according to the conformation according to equation (A.20).

The evolution of the drop morphology and the polymer chain conformation is simulated simultaneously and separately and are coupled through rewriting the equilibrium function ((A.24)). That is, the elastic stress originating from the dumbbell conformation act as a forcing term in Guo’s scheme [76].

2.2. 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. This 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. This 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.

The 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. The viscosities of both the drop and matrix components, and , are set to be 1000 Pa s, the density, and , 1000 kg/m3, and the interfacial tension between the drop and the matrix, , 0.116 mN/m.

The pseudopotential parameter of LB is related to according to the Laplace law and decided through trials and error. The spring constant is set to be 1.

These 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.

The strength of the polymer viscoelasticity is characterized by the Deborah number, as expressed below:which defines the ratio of relaxation time and the characteristic time . It is obvious that the larger the De, the stronger the viscoelasticity, and a zero De indicates purely viscous polymer.

The shape of the ellipsoid drop is defined by two parameters, equivalent radius and the deformability :

Of all the simulation cases, the initial values of and 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 , is defined:

In terms of the boundary conditions, periodic boundaries are applied on all the outer boundaries of the simulation domain.

Three 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.

The 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.

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 ), the retraction rate of the drop gradually slows down. This 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 to the time are utilized. The duration between and is the first stage of the drop retraction, and empirically, the value is more appropriate to set −1.5.

It should be noted that this treatment also agrees with the dumbbell model. The 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. Therefore, 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.

3. Experimental Study

3.1. 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.

The 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.

3.1.1. 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.

Then, the POE sheet is molded out on a compression molding press. The 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.

Through 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.

3.1.2. 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.

3.1.3. Experimental Observation and Data Processing

The 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, and , at each shot time and the radius R when the spherical drops are restored.

4. Results and Discussion

The 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, and , at different times is measured by the software Scion image and then the variation in the deformation recovery degree over time is calculated, as shown in Figure 10. A common first-fast-then-slow 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.

There 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. 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.

The 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. Their measurement results are shown in Table 3. The 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 is used to fit equation (4).

At the bottom of Table 3, the results fluctuation characterizes the impact of drop size on the measurement results:where A is the set of all the measured data, is the minimum value in A, and 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%. The 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. This measurement value is approximately the one obtained by the proposed method in all the three cases of different drop sizes. The 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. There 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 between the same polymers of different molecular weights is not the same and needs measuring separately.

5. 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. Then, a novel evaluation method of the interfacial tension is put forward that only takes advantage of the first stage of the experimental data. The 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. Therefore, 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 is the lattice speed, is the lattice spacing, and is the time step.

Using the multirelaxation time collision operator MRT, the governing equation of can be written aswhere is the collision operator and defined by

is the density distribution function associated with velocity at position and time , and their moments satisfywhere represents the interaction force between components. is the discrete forcing term accounting for the interaction force :where is the unit matrix, while is defined by

In the collision operator, is the constant transformation matrix and is a diagonal nonnegative relaxation time matrix that is decided by fluid properties.

For multicomponent flow, is the corresponding equilibrium distribution of component and is determined bywhere is the weight associated with the velocity :

is the equilibrium macroscopic velocity in consideration of the component interaction:where satisfies the momentum equation without the external force:in which and are the macroscopic density and velocity of component.

is the interaction between the particles of different components, when only considering the nearest neighbor, it can be written aswhere denotes the interaction strength between components and . When is positive, the repulsive force is generated between the particles, and when the 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:

B. FENE Dumbbell Model

Likewise, the conformation distribution function follows the discrete Fokker–Planck equation on D3Q19 lattice:

As the polymer melt is thought isotropic at equilibrium, the equilibrium conformation distribution function is equal to :

is the single relaxation time of towards accounts for the elongation and rotation of the dumbbell:

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. Then, the macroscopic elastic stress is linked with the dumbbell conformation through the Kramers formula:where is the spring force acted on the dumbbell beads:

When applied the Chapman–Enskog expansion, the equation (A.15) results in the macroscopic Oldroyd-B constitutive model:

The model parameters and are linked to that of the lattice Fokker–Planck equation by the following relations:

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 :

It is noticeable that the modified expression reflects the effect of polymer viscoelasticity through the introduction of into . After the modification, the second-order moment of the distribution function contains both viscous and elastic stress:where is the viscous stress, , so that viscoelastic Navier–Stokes equation can be recovered from equation (A.2).

The numerical procedure for each time step at each lattice site reads as follows:(1)Stream the component distribution functions between the adjacent lattice sites (left-hand side of equation (A.2)) and calculate the velocity and density (equation (A.3))(2)Calculate the convection of the dumbbell conformation (equation (A.20)) and the elongation and rotation of the dumbbells due to the polymer flow (equation (A.18))(3)Update the new conformation functions (equation (A.15)(4)Sum the viscous and elastic stress to calculate the total local stress (equation (A.25))(5)Collide and update the new component distribution functions (right-hand side of equation (A.2))(6)Repeat the procedures (1)∼(5) until the simulation ends

Data Availability

The experimental data used to support the findings of this study are included within the article. The 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

The authors declare that they have no conflicts of interest.


This work was financially supported by the National Natural Science Foundation of China (no. 51805379).