Research on the Mesoscopic Characteristics of Kelvin–Helmholtz Instability in Polymer Fluids with Dissipative Particle Dynamics

Abstract

In this paper, the two-dimensional Kelvin–Helmholtz (KH) instability occurring in the shear flow of polymer fluids is modeled by the dissipative particle dynamics (DPD) method at the coarse-grained molecular level. A revised FENE model is proposed to properly describe the polymer chains. In this revised model, the elastic repulsion and tension are both considered between the adjacent beads, the bond length of which is set as one segment’s equilibrium length. The entanglements between polymer chains are described with a bead repulsive potential. The characteristics of such a KH instability in polymer fluid shear flow can be successfully captured in the simulations by the use of the modified FENE model. The numerical results show that the waves and vortexes grow more slowly in the shear flow of the polymer fluids than in the Newtonian fluid case, these vortexes become flat, and the polymer impedes the mixing of fluids and inhibits the generation of turbulence. The effects of the polymer concentration, chain length, and extensibility are also investigated regarding the evolution of KH instability. It is shown that the mixing of two polymer fluids reduces, and the KH instability becomes more suppressed as the polymer concentration increases. The vortexes become much longer with the evolution of the elongated interface as the chain length turns longer. As the extensibility increases, the vortexes become more flattened. Moreover, the roll-up process is significantly suppressed if the polymer has sufficiently high extensibility. These observations show that the polymer and its properties significantly influence the formation and evolution of the coherent structures such as the waves and vortexes in the KH instability progress.

1. Introduction

It is well known that Kelvin–Helmholtz (KH) instability occurs at the interface between two parallel streams with a relative velocity difference [1]. KH instability has been an interesting and challenging subject in that the interface undergoes a complex evolution. It should be noted that the KH instability can be encountered in the flow of both Newtonian [2] and non-Newtonian fluids [3,4,5]. This phenomenon is of fundamental importance in science research and engineering applications, such as marine science [6], turbulent mixing [7,8], air-assisted atomization [9], material friction [10], nuclear engineering [11], and so on.

So far, the majority of the KH instability studies are limited to Newtonian fluids [12,13,14,15,16,17,18,19,20,21]. That occurring in the flow of non-Newtonian fluids has rarely been researched. Cadot and Kumar [22] experimentally studied the viscoelastic effect on shear instabilities. Their results showed that the viscoelastic solutions drastically modified the instabilities. Azaiez and Hornsy [23] numerically investigated the effect of viscoelasticity on the free shear flow of viscoelastic liquids. They also simulated the roll-up of non-Newtonian free shear flow at high Reynolds numbers [24]. Kumar and Homsy’s simulations of the KH instability [25] showed how the addition of a polymer affected the transition to turbulence in free shear layers. In Yu and Phan-Thien’s simulations [26], the functions of the added polymer were discussed in detail. Richter et al. [27] simulated the three-dimensional viscoelastic flow pasting a circular cylinder and investigated the roll-up instability formed in the wake. Tiwari et al. [28] simulated the evolution of a shear flow structure in viscoelastic fluids with a generalized hydrodynamic model. They observed that the growth rate of the KH instability diminished with the increasing phase velocity of the transverse wave. In these numerical simulations, the fluid flows are just roughly described with a macroscopic continuum model.

As a great deal of interest is focused on the physical phenomena and industrial applications that involve the complex flow at a small scale, it is of great importance to understand the KH instability occurring in the non-Newtonian fluids at a mesoscale. As the system scale decreases, thermal fluctuations have a greater influence on the fluid dynamics. What is worse, continuum formulations based on partial differential equations, i.e., the Navier–Stokes equations, are even unable to describe the KH instability problem at such scale without modifications. As a mesoscopic simulation method, the lattice Boltzmann method (LBM) may be employed to model the vortex formation at a small scale for the KH instability problem [29,30,31,32]. However, it does not work when modeling the transition from laminar flow to turbulence due to the restriction of the lattice. These characteristics of complex fluids at a small scale may also be modeled by the use of molecular dynamics (MD) [33,34]. However, MD simulation needs too high of a calculation load to model routine problems [35]. Hence, dissipative particle dynamics (DPD) are considered to be the first choice for investigating the KH instability that occurs in the flow of polymer fluids in this work.

DPD is a coarse-grained method that connects microscopic and macroscopic worlds. This method was originally conceived as a mesoscale technique for simulating fluid systems by Hoogerbrugge and Koelman [36]. Subsequently, Español and Warren [37] and Marsh [38] established a wholesome theoretical basis for DPD based on statistical mechanics. In recent years, as the small-scale technologies developed fast, DPD has attracted increasing attention, and its applications have been rapidly extended to various research fields. These applications include colloidal suspension [39], liquid–liquid two-phase flow [40], gas–liquid two-phase flow [41], electroosmotic flow [42], liquid nanojets [43], solid–liquid interactions [44], particle sedimentation [45], blood flow [46], colloid transport in porous media [47], oil-in-water emulsion [48], magnetorheological fluids [49], etc. In contrast to these aspects, the applications of DPD in polymer systems are more extensive, such as polymer solutions [50], macromolecular suspension [51,52], topological entanglements in polymer melts [53,54], the phase behavior of polymer melts [55], the scaling behavior of polymers [56], and so on. For more details about the DPD method, see Moeendarbary et al. [57] and Liu et al. [58]. It should be noted that the DPD method has also been applied to study the phenomenon of fluid instability at the mesoscopic particle level [59,60,61,62,63,64], but most of them just considered the instability occurring in the flow of simple fluids.

Polymer fluids are non-Newtonian fluids with polymers. Polymers usually have long polymer chains and a large molecular weight. According to the concentration of polymers, polymer fluids can be divided into polymer melts, polymer concentrated solutions, and polymer dilute solutions. Polymer chains in polymer fluids can produce interactions such as entanglement and interpenetration, which affect the rheological properties of polymer fluids. For the high-concentration polymer fluids, the interaction forces between molecules are strengthened. The entanglement and interpenetration between polymer chains are also enhanced, making the rheological properties of polymer fluids more complex.

Studying the KH instability that occurs in the flow of non-Newtonian fluids such as polymer fluids will gain good cognition of non-Newtonian fluids’ effect on the instability in complicated flows, especially at a small scale. For investigating the dynamic behaviors of polymer fluids, such as KH instability, it is important to construct an effective model for describing the polymers. In simulations of polymers at micro- or meso-scales, the extendable non-linear elastic (FENE) chain is widely used owning to its vivid description of the polymer chain motion. The drawback of the traditional FENE chain model is that it only considers the elastic tension between adjacent beads of a chain; the elastic repulsion, however, is ignored [65,66].

This article presents a modified FENE chain model in which the interaction between adjacent beads is determined with consideration of both elastic tension and elastic repulsion. The entanglement model is given and employed in the simulations.

2. Problem Setup

The two-dimensional rectangle computational domain is $\Omega =\left[-L/2,L/2\right] \times$ $[-H/2,H/2]$, where the length of the computational domain $L$ is chosen to be quadruple the domain height $H$. The computational domain for the KH instability problem is represented by a set of particles with an equidistant particle distribution. At the beginning of the simulation, the computational domain is halved by a midline $y=0$, where each half represents the different fluid regions with the same number of particles. Initially, the particles of two-layer fluids are set into thermal motion and superposed shear velocities in opposite directions with the same magnitude, that is:

$$U=\left\{\begin{array}{l}{u}_0 y\ge 0\\ -u_0 y<0\end{array}\right.$$

(1)

where $u_0$ is the magnitude of shear velocities. Periodic boundary conditions are implemented at $x=-L/2$ and $x=L/2$, and wall boundary conditions are implemented at $y=-H/2$ and $y=H/2$, where each solid wall is represented by two rows of wall boundary particles, setting their velocities to $u_0$ and $-u_0$, correspondingly, throughout the simulation. The no-slip boundary conditions on boundary walls are implemented by Maxwellian reflection. In addition, it should be pointed out that in the DPD simulation, thermal fluctuations can spontaneously initiate the interface instability in shear flow, unlike the macroscopic continuum model in which initial sinusoidal perturbations are artificially implemented to the fluid–fluid interface.

3. DPD Method and Polymer Model

3.1. The DPD Method

In DPD simulation, the fluids are modeled by a system of particles, and the motions of these DPD particles are governed by Newton’s second law of motion. The motion of particle $i$ is determined as [67]

$$\frac{d{r}_i}{dt}=v_i$$

(2)

$$m_i\frac{d{v}_i}{dt}=F_i$$

(3)

where $r_i$ is the particle position, $v_i$ is the particle velocity, $m_i$ is the particle mass, and $F_i$ is the resultant force exerted on the particle $i$ by all the other particles. In the standard DPD model, the resultant force $F_i$ contains three parts: the conservative force $F_i^C$, the dissipative force $F_i^D$ and the random force $F_i^R$, which are pairwise and center-to-center. The resultant force $F_i$ is written as [67]

$$F_i={\displaystyle \sum _{j\left(\ne i\right)}{F}_{ij}^C}+{\displaystyle \sum _{j\left(\ne i\right)}{F}_{ij}^D}+{\displaystyle \sum _{j\left(\ne i\right)}{F}_{ij}^R}$$

(4)

where the conservative force $F_{ij}^C$ is a soft repulsive force acting along the line of centers. $F_{ij}^C$ is calculated as [67]

$$F_{ij}^C=a_{ij}{w}_C\left(r_{ij}\right)e_{ij}$$

(5)

$$w_C\left(r_{ij}\right)=\left\{\begin{array}{l}1-\frac{r_{ij}}{r_c} r_{ij}\le r_c\\ 0 r_{ij}>r_c\end{array}\right.$$

(6)

where $a_{ij}$ is a parameter representing the strength of a repulsive force, $r_{ij}=r_i-r_j$ is the relative position vector between particle $i$ and particle $j$, $r_{ij}=\left|r_{ij}\right|$, $e_{ij}$ is a unit vector given by $e_{ij}=r_{ij}/r_{ij}$, $r_c$ is the cutoff radius, and $w_C\left(r_{ij}\right)$ is the weight function for conservative force.

In order to simulate liquid–liquid interface with the DPD method, Groot and Warren [67] have made a link between the repulsion parameter and Flory–Huggins theory. For the number densities $n=3$ and $n=5$, they have also obtained the relationships between the Flory–Huggins parameter $\chi$ and the repulsion parameter $a_{ij}$ for two-type particles. Subsequently, AlSunaidi et al. [68] have obtained the relationship between $\chi$ and $a_{ij}$ for the number density $n=4$ by interpolating the equations provided for number densities $n=3$ and $n=5$, that is [62]:

$$a_{ij}=\left\{\begin{array}{l}\alpha \mathrm{if}\mathrm{particles}i\mathrm{and}j\mathrm{are}\mathrm{the}\mathrm{same}\mathrm{type}\mathrm{particle}\\ \alpha +2.05\chi \mathrm{if}\mathrm{particles}i\mathrm{and}j\mathrm{are}\mathrm{different}\mathrm{type}\mathrm{particles}\end{array}\right.$$

(7)

So, as far as two-dimensional simulation is considered, initially, two particles are used for each square cell, which corresponds to the face-centered cubic with four particles per cubic lattice in a three-dimensional system. The size of the square cells and cubic lattices keeps unit, which means the side length of all the cells and lattices is equal to 1, and the corresponding number density is 4. Therefore, in our simulation, the number density can be deemed as $n=4$, and Equation (7) is taken to model the binary fluids. When the Flory–Huggins parameter $\chi$ is positive, the binary fluids are immiscible, while they are miscible when it is zero or negative.

In Equation (4), $F_{ij}^D$ is calculated as [67]

$$F_{ij}^D=-\gamma w_D\left(r_{ij}\right)\left(e_{ij}\cdot v_{ij}\right)e_{ij}$$

(8)

where, $\gamma$ is the strength of the dissipative force, $w_D\left(r_{ij}\right)$ is the weight function for $F_{ij}^D$, and $v_{ij}=v_i-v_j$ is the relative velocity. $F_{ij}^R$ is calculated as [67]

$$F_{ij}^R=\sigma w_R\left(r_{ij}\right){\zeta }_{ij}{e}_{ij}$$

(9)

where $\sigma$ is the strength of random force, $w_R\left(r_{ij}\right)$ is the weight function for $F_{ij}^R$, $v_{ij}=v_i-v_j$ is the relative velocity, and ${\zeta }_{ij}$ is the stochastic variable with the properties [67]

$$〈{\zeta }_{ij}〉=0$$

(10)

$$〈{\zeta }_{ij}\left(t\right){\zeta }_{i^{\prime }{j}^{\prime }}\left(t^{\prime }\right)〉=\left({\delta }_{i{i}^{\prime }}{\delta }_{j{j}^{\prime }}+{\delta }_{i{j}^{\prime }}{\delta }_{j{i}^{\prime }}\right)\delta \left(t-t^{\prime }\right)$$

(11)

where $〈〉$ denotes an ensemble average and ${\delta }_{ij}$ is the Kronecker delta.

This is according to the fluctuation dissipation theorem [67]

$$w_D\left(r_{ij}\right)=w_R^2\left(r_{ij}\right)$$

(12)

and

$${\sigma }^2=2\gamma k_{B}T$$

(13)

where $k_{B}T$ is the Boltzmann temperature of the system. The conventional dissipative weight function is adopted as [67]

$$w_D\left(r_{ij}\right)=w_R^2\left(r_{ij}\right)=\left\{\begin{array}{l}{\left(1-\frac{r_{ij}}{r_c}\right)}^2 r_{ij}\le r_c\\ 0 r_{ij}>r_c\end{array}\right.$$

(14)

As can be seen from (12) and (14), the cutoff radii for the dissipative force should be equal to that for the random.

3.2. The Polymer Model

The adjacent particles on a polymer backbone can be linked together by the FENE spring force as follows [66].

$$F_{ij}^S=-\frac{H_S{r}_{ij}}{1-{\left(r_{ij}/r_{\max }\right)}^2}$$

(15)

where $H_S$ is the spring coefficient, $r_{ij}=r_i-r_j$, $r_{ij}=\left|r_{ij}\right|$, $r_i$, and $r_j$ denote the position vectors of bead $i$ and bead $j$, and $r_{\max }$ refers to the maximum extensibility of one chain segment. However, the spring force (15) describes only the elastic tension, and the elastic repulsion is neglected.

To model the polymer chain more veritably, a modified FENE spring force [69] is adopted as follows.

$${\overline{F}}_{ij}^S=-\frac{H_S\left(r_{ij}-r_{bond}\right)r_{ij}}{\left[1-{\left(r_{ij}/r_{\max }\right)}^2\right]\left(r_{\max }-r_{bond}\right)}$$

(16)

where $r_{bond}$ represents the bond length, and other physical quantities have the same meaning as the corresponding parameters in the original model (15). Figure 1 shows the schematic diagram of the elastic tension and repulsion in a polymer chain.

In polymer fluids, the uncrossability constraint creates entanglements when the chains try to cross each other, as shown in Figure 2. At present, there are some effective models that have been developed to deal with the bond crossing problem. Goujon et al. [70] and Sirk et al. [71] used a segmental repulsive potential to model polymer entanglements in DPD simulation. In order to avoid bond crossing in a two-dimensional system, Liu et al. [72] added a rigid core to each chain bead to prevent the beads from penetrating each other.

In this work, the above method is also employed, and a repulsive force with the same functional form as the DPD conservative potential is given as follows [63]:

$$F_{ij}^E=\left\{\begin{array}{l}{a}_{ij}^E\left(1-\frac{r_{ij}}{r_c^E}\right)e_{ij} r_{ij}\le r_c^E\\ 0 r_{ij}>r_c^E\end{array}\right.$$

(17)

where $F_{ij}^E$ is a force acting between beads $i$ and $j$ which belong to different chains with the distance $r_{ij}$, and it is called the entanglement force, $a_{ij}^E$ is the controlling parameter for the entanglement force, $r_c^E$ represents the cutoff value for the entanglement interactions, and the entanglement cutoff distance should be less than the cutoff distance for traditional DPD forces.

3.3. Integration Algorithms

Newton’s equations of motion (Equations (2) and (3)) are discretized with a modified version of the velocity–Verlet algorithm, including Equations (18)–(21):

$$r_i\left(t+\Delta t\right)=r_i\left(t\right)+v_i\left(t\right)\Delta t+\frac{1}{2}\Delta t^2{\dot{v}}_i\left(t\right)$$

(18)

$${\tilde{v}}_i\left(t+\Delta t\right)=v_i\left(t\right)+\lambda \Delta t{\dot{v}}_i\left(t\right)$$

(19)

$${\dot{v}}_i\left(t+\Delta t\right)={\dot{v}}_i\left(r_i\left(t+\Delta t\right),{\tilde{v}}_i\left(t+\Delta t\right)\right)$$

(20)

$$v_i\left(t+\Delta t\right)=v_i\left(t\right)+\frac{1}{2}\Delta t\left[{\dot{v}}_i\left(t\right)+{\dot{v}}_i\left(t+\Delta t\right)\right]$$

(21)

where $\Delta t$ is the DPD time step, and ${\dot{v}}_i$ is the particle acceleration caused by the forces adding to particle $i$ from all other particles. The intermediate velocity ${\tilde{v}}_i\left(t+\Delta t\right)$ is predicted before the forces are updated, and then this value is corrected at the end of each integration step; $\lambda$ is an empirical parameter, and $\lambda =0.65$ is taken as the optimum value [67]. Here, the modified velocity–Verlet algorithm ((18)–(21)) is based on the Euler method to discretize Newton’s second law equations, so it is stable. To analyze the convergence rate of the modified velocity–Verlet algorithms (18)–(21), Taylor expansions were performed on displacement and velocity, respectively

$$r\left(t+\Delta t\right)=r\left(t\right)+v\left(t\right)\Delta t+\frac{\Delta t^2}{2}\dot{v}\left(t\right)+O\left(\Delta t^3\right)$$

(22)

$$\mathbf{v}\left(t+\Delta t\right)=v\left(t\right)+\dot{v}\left(t\right)\Delta t+O\left(\Delta t^2\right)$$

(23)

It can be seen from (22) and (23) that the displacement is third-order convergence, and the velocity is at least first-order convergence. In addition, the mass of each particle is fixed, and the number of particles is constant in the simulation. Thus, the calculation system is mass conservation. The momentum of the system is related to the dissipated force, which consumes some momentum, and the boundary velocity, which supplements some momentum. Although the system does not maintain absolute conservation of momentum, the momentum fluctuations are so small that it can be considered basically conserved. In addition, for all calculations conducted in this study, the control equations and the results are represented by dimensionless variables using the reference mass, length, and energy, i.e., the particle mass $\left[m\right]$, the side length $\left[l\right]$ of the unit mesh, and the Boltzmann temperature $\left[k_{B}T\right]$.

4. Simulation Results and Discussion

4.1. KH Instability in Newtonian Fluids

Consider the KH instability occurring in the flow of Newtonian fluids with twofold purposes. One is to verify the DPD code for the KH instability problem, and the other is to provide a counterpart as a contrast for the case of polymer fluids that is mainly discussed in this work. The KH instability in two-component Newtonian fluids is simulated with the initial stream velocity:

$$U=\left\{\begin{array}{l}1 y\ge 0\\ -1 y<0\end{array}\right.$$

(24)

It is on the domain Ω = [–500, 500] × [−126, 126]. In this two-dimensional system, the triangular lattice with two particles per unit cell is adopted in the initial coordinates. The right and left boundaries are the periodic boundary; the top and bottom boundaries are the wall boundary, and each wall boundary is composed of two layers of particles. The fluid particle number is 500,000, and the wall particle number is 4000. In this simulation, the two-component fluids have the same density, and the remaining parameter values are $\alpha =18.75$, $\gamma =4.5$, $r_c=1.125$, $k_{B}T=1$, $\chi =5$, and $\Delta t=0.01$. As discussed in the literature [73], the gravitation is neglected in this test.

Figure 3 gives the temporal evolution of the fluid particles to illustrate the evolution of the KH instability occurring in the flow of Newtonian fluids. The upper layer of the fluids moves in the positive direction, while the lower layer of the fluids moves in the negative direction. There exists distinct penetration between both fluids. Early meso-level inter-diffusion causes a broadening of the interface, where meso-level thermal fluctuations perturb the interface and trigger the growth of small waves. As time goes on, these small waves grow in amplitude. Then, both fluids tend to gain horizontal velocity opposite to their initial bulk velocities and eventually crest to form small vortices. As the KH instability grows in the y-axis direction, the wavelength of the waves and vortices grows in the x-axis direction. Initially, small vortices grow with short wavelength modes, but ultimately, the larger vortex structures grow at the expense of the small ones. This process results in the formation of vortices, and the characteristic form of KH instability becomes much more pronounced at a later time. The roll-up interface dynamics are qualitatively consistent with those in the previous experiments [12] and simulations [33].

4.2. KH Instability in Polymer Fluids

The KH instability occurring in the flow of polymer fluids is a more complex physical process than that in the Newtonian fluids, which will be discussed in the following content. As mentioned earlier, Equations (16) and (17) are used for the polymer solution. Here, the spring coefficient $H_S=4.5$, the bond length $r_{bond}=\sqrt{2}/2$, the maximum extensibility of one chain segment $r_{\max }=2.5r_{bond}$, the chain length ${\rm X}n=4$, the polymer concentration $c=40\%$, the entanglement force controlling parameter $a_{ij}^E=8$, the entanglement cutoff distance $r_c^E=0.8r_c$, and other parameter values are similar to those used in the Newtonian case.

Figure 4 shows the instantaneous snapshots of the fluid particles to illustrate the evolution of the KH instability occurring in the flow of polymer fluids. In the initial time $t=0$, the shear flow starts with the upper-layer fluid moving in a positive horizontal direction with bulk velocity $U=1$ and the lower-layer fluid moving in the opposite direction with bulk velocity $U=-1$. Subsequently, the interface disturbance generates spontaneously due to the thermal fluctuations, and the flow system becomes unstable. Similar to the Newtonian case, because of the shear action and the thermal disturbance, the lower-layer fluid starts moving in a positive vertical direction, whereas the upper-layer fluid starts moving in the opposite direction. As a result, both polymer fluids begin to penetrate each other. During this initial period, unlike the Newtonian case that generates many small vortexes near the original interface, it is difficult to form small vortexes in the shear flow of polymer fluids on account of the inhibition of the polymer. Instead, the interface instability occurring in the flow of polymer fluids manifests itself as some finger-like structures and needle-like structures developing at the interface. As time progresses, these finger-like structures and needle-like structures break and merge. Gradually, some waves generate at the interface with a long wavelength, which subsequently roll up and form vortexes. In contrast with the Newtonian case, the two polymer fluids have a lower mixing level, and the curly interface is significantly elongated and does not break up. This is because the polymer chains make the liquid micelles connect together. In the Newtonian case, under shear flow, the interface becomes more incompact and manifests plume-like structures, namely, the stream gradually changes from laminar flow to turbulence. However, the vortex structures in the polymer fluids are compact cluster structures instead of plume-like patterns. Thus, it can be seen that the polymer suppresses the growth of turbulence. Due to the inhibiting effect of the polymer on the interface instability, it needs more time to generate the waves and vortexes, and these vortex structures are flattened compared to the Newtonian ones. Our simulation results are consistent with previous experimental and numerical research revealing that the polymer inhibits the roll-up process [22,25,27].

Figure 5 shows the flow velocity in the vertical direction of Newtonian fluids and polymer fluids. In particular, the suppression influence of the polymer on the formation of small vortexes is clearly reflected in vertical velocity (see the top frames of Figure 5). In the initial, for the Newtonian case, many small vortexes appear near the original interface, where a small red zone (left, up-flow velocity) and a small blue zone (right, down-flow velocity) correspond to a vortex (see the top frame of Figure 5a). However, the vertical velocity is very small in the polymer case since there is almost no small vortex appearing at the interface (see the top frame of Figure 5b). In the late, the interface forms some large vortex structures, and the high vertical velocity is generated in the polymer fluids, which is just slightly lower than the Newtonian case (see the bottom frames of Figure 5). In conclusion, due to the inhibition of the polymer, the vortex of Newtonian fluid is formed earlier, while the vortex of polymer fluid is formed later. The vortex velocity of polymer fluid is smaller than that of Newtonian fluid at the same time point. Here, in order to reduce the random noise, the values are averaged over multiple steps to obtain the results at the midpoint of the time interval for the vertical velocity, and so are the horizontal velocity and velocity vector-field in subsequent sections.

Figure 6 presents the vortex growth occurring in Newtonian fluid and polymer fluid, which are mean values that are obtained by averaging the vortex growth for several instances. From this figure, it can be seen that the vortexes develop more quickly in Newtonian fluid than they do in polymer fluid.

The above phenomenon can be explained by the following. When the molecular chain of a polymer is in a flow field with a velocity gradient caused by shear rate changes, each part of the long chain may be in different velocity regions: one end may be in the faster velocity region, while the other end may be in the slower velocity region, and even the flow velocity direction may be opposite in different regions where different parts of the molecular chain are located. At this point, relative movement occurs at both ends of the molecular chain, which may cause stretching and orientation of the molecular chain. The molecular chain interacts with the fluid during stretching and orientation, resulting in the HK vortex being stretched and flattened.

4.3. Effect of Polymer Concentration

To investigate the effect of polymer concentration, the same initial velocity conditions and other parameter values were used to create Figure 4, except for the concentration. Consider the fluid particles, velocity vector-field, horizontal velocity, and vertical velocity for the concentrations of $c=10\%$, $c=30\%$, and $c=70\%$, respectively. In our simulations, the vortexes that reach saturation can be captured before $t=700$. However, it is much later for $c=70\%$ to form a vortex.

Figure 7 shows the fluid particles, velocity vector-field, horizontal velocity, and vertical velocity. In the case of $c=10\%$, as the upper shear layer moves in a positive level direction and the lower shear layer moves in the opposite direction, the waves generate, grow, coalesce, get larger, and form vortexes, as expected. In this case, the wave evolution and the vortex roll-up process are similar to what are seen in the Newtonian case, but the vortexes are slightly flattened relative to the Newtonian ones. As the polymer concentration increases, $c=30\%$, the polymer stabilizing effect on the interface evolution becomes obvious. The formation of the wave and vortex patterns is slightly postponed in contrast to the foregoing case, and the vortexes become more flattened. For a higher polymer concentration, $c=70\%$, at the beginning, the waves develop in the vertical direction with a small increasing amplitude. As time progresses, the vertical growth of the waves becomes slower and slower. Instead, finger-like structures begin to develop and are almost parallel to the horizontal direction. Gradually, these finger-like structures are elongated. At last, these finger-like structures break up and form some drops. The vortex patterns are not seen at any time. One can easily notice that, with the polymer concentration increasing, the mixing of two polymer fluids reduces, the KH instability becomes more suppressed, and the vortex formation is delayed to a later simulation time. Peculiarly, for the case of a high polymer concentration, the vortex structures are completely inhibited at the mesoscopic level. The effect of the polymer concentration on the KH instability is also reflected in the fluid velocity. As the polymer concentration increases, eddy velocities weaken. In particular, when the polymer concentration rises to $c=70\%$, the vertical velocity is very small, and almost no eddy current can be observed in the velocity vector-field. Here, the simulation results are consistent with the notion that the polymer can inhibit the formation of concentrated regions of the vortex [25].

This arises due to the complex interaction between the molecular chains. As the concentration of the polymer increases, the number of polymer molecules increases, and the interactions between molecular chains become frequent. It is the process of entanglement and unwinding in the molecular chains that hinders the formation and development of KH vortices.

4.4. Effect of Chain Length

Consider the effect of chain length on the evolution of KH instability. Here, the polymer concentration is 30%. The same initial velocity conditions and other parameter values were utilized to produce Figure 4, except for the chain length. Consider in the following the fluid particles, velocity vector-field, horizontal velocity, and vertical velocity at the moment that the vortexes reach saturation for the chain lengths of ${\rm X}n=2$, ${\rm X}n=8$, and ${\rm X}n=30$, respectively.

Figure 8 shows snapshots of the fluid particles, velocity vector-field, horizontal velocity, and vertical velocity at the moment that the vortexes reach saturation. In the case of ${\rm X}n=2$, i.e., the FENE dumbbell model, at the early time, some small waves form at the interface, which grow in the vertical direction with increasing amplitude. Gradually, these waves merge and evolve into vortexes with a larger amplitude. This process is very similar to the Newtonian case. In the case of longer chains, ${\rm X}n=8$, it can be observed that the waves grow and roll up with a slower rate relative to the case of ${\rm X}n=2$. The vortexes grow up with a small altitude, which are elongated by following the evolution of the shear flow. When the chain length becomes much longer, ${\rm X}n=30$, the interface roll-up is seriously retarded, the formation of vortexes needs more time, and the vortexes become longer with the evolution of the elongated interface. From the velocity vector-field, horizontal velocity, and vertical velocity, it can be seen that with the chain length increasing, the eddy currents weaken markedly.

The flow characteristics of polymer fluids are closely related to the structure of molecular chains. Long-chain molecular structures involve a wide range of flow regions, and long-chain molecules are more prone to entanglement. Therefore, large vortices are easily generated under the interaction of these long-chain molecules and fluid molecules. In addition, these long-chain molecules undergo stretching and orientation under shear, which simultaneously react on the fluid, causing the vortices to flatten.

4.5. Effect of Polymer Extensibility

When the polymer concentration is 30%, the same initial conditions and other parameter values are adopted to create Figure 4, except for the parameter $r_{\max }$. Figure 9 shows snapshots of the fluid particles, velocity vector-field, horizontal velocity, and vertical velocity for the extensibilities of $r_{\max }=2r_{bond}$, $r_{\max }=7r_{bond}$, and $r_{\max }=20r_{bond}$, respectively. For $r_{\max }=2r_{bond}$ and $r_{\max }=7r_{bond}$, these snapshots are captured at the moment that the vortexes reach saturation. For $r_{\max }=20r_{bond}$, the snapshots are captured at a certain moment in a later stage. For the case of $r_{\max }=2r_{bond}$, the roll-up vortexes and the mixing of two fluids are similar to what are seen in the lower-concentration and Newtonian cases. However, for the higher extensibility, $r_{\max }=7r_{bond}$, it is obviously postponed to form the vortex structures, and the roll-up vortexes appear to be strongly flattened relative to the vortexes at $r_{\max }=2r_{bond}$. This flattening is somewhat similar to what is seen for an increasing chain length. The simulations are consistent with the results of Kumar and Homsy [25], which show that the roll-up process is significantly suppressed if the polymer has sufficiently long extensibility. When the parameter $r_{\max }$ further increases, $r_{\max }=20r_{bond}$, the roll-up process of the interface is significantly inhibited. For this case, at the early time, the waves grow in the vertical direction with a small amplitude and, at the same time, merge in the horizontal direction. As time goes by, the vertical growth of the waves gradually slows down. Their horizontal growth goes on as the finger-like structures are elongated longer and longer. It is also observed that a high extensibility reduces eddy velocities.

For molecular chains with high extensibility, the stretching deformation scale of the molecular chain is larger under shear, and the deformation relaxation time is longer, which means that the molecular chains with orientation during flow do not easily restore their original shape. Therefore, under the interaction of these highly extensible molecular chains with fluids, vortices become flatter and even develop into finger-like structures.

5. Summary and Conclusions

The KH instability occurring in the shear flow of polymer fluids is numerically researched using DPD at the mesoscopic molecular level. The polymer chains are modeled with the modified FENE model that describes both the elastic tension and the elastic repulsion between the adjacent beads with the bond length as the equilibrium length of one segment. Meanwhile, the entanglements between polymer chains are modeled with a simplified bead repulsive potential. The following results are obtained through numerical simulation.

(1)

The DPD algorithm is validated by simulating the KH instability that occurs in the shear flow of Newtonian fluids, and the mesoscopic results of the present simulation for Newtonian fluids are qualitatively consistent with previous experimental and numerical results.

(2)

The characteristics of the KH instability that occurs in the flow of polymer fluids are explored with the validated DPD code. It is shown that, in contrast with their Newtonian counterparts, the waves and vortexes forming in the shear flow of polymer fluids grow more slowly, and these vortexes are flattened. In particular, the roll-up interface is conspicuously elongated and does not break up due to the reason that the polymer chains make the liquid micelles join together. Thus, the polymer hinders the mixing of fluids and suppresses the generation of turbulence. It is also observed that there is almost no small vortex appearing in the development process of the KH instability due to the inhibition of the polymer.

(3)

The effects of the polymer concentration on the KH instability are researched. The results show that with the polymer concentration increasing, the inhibitory effect becomes more conspicuous, the mixing of two polymer fluids reduces, and the transition from laminar to turbulence is even more restrained.

(4)

The effects of the chain length and extensibility of the polymer on the KH instability are also investigated. For the case of a long chain or high extensibility, the effect of the polymer on the flow structure is more profound. When the chain length gets longer, even the concentration remains constant, the formation of the waves and vortexes is delayed to a later simulation time, and the vortexes become longer and more flattened. This kind of situation also occurs when the polymer extensibility increases. The most striking effect of the high extensibility is the complete suppression of all vortex structures.