Interface instabilities and chaotic rheological responses in binary polymer mixtures under shear flow

Abstract

A unified model by combining the Rolie–Poly constitutive model and the Flory–Huggins mixing free energy functional through a two fluid approach is presented for studying flow-induced phase separation in polymer mixtures. It is numerically solved in two-dimensional flow with monotonic and non-monotonic constitutive behaviour. The results are analyzed and show that this model can capture the essential dynamic features of viscoelastic phase separation reported in literature. The steady-state shear-banding and interface instabilities are reproduced. In the case with a non-monotonic constitutive behaviour, It is observed that the band structures are strongly unstable both in time and in space. The correlations between the microstructure evolution and chaotic rheological responses have been identified. A vortex structure emerges within the central band. Numerical results obtained from this study suggest that the dynamic features of rheochaos can be captured by the proposed model without introducing extra parameters.

---

1 Introduction

Shear-induced phase separation in complex fluids has attracted great attention over the last decade. Nearly 60 years ago, it was observed¹ that phase transitions of polymer solutions are influenced by flow. Since then there have been many reports on shear-induced phase separation of polymer mixtures and other complex fluids.^2–7 As a closely related phenomenon, shear-banding was also observed in numerous publications, recent experiments include wormlike micelles,^8,9 telechelic polymers,^10,11 dispersions¹² and entangled DNA solutions.^13–15

A typical experimental system targeted by this work is the Couette cell described in ref. 8. Fluids are contained in the gap between the stationary inner cylinder and the rotating outer cylinder. Various shear rates are applied by superimposing different angular velocities to the rotating cylinder. Focused on the bulk fluid behaviour, we simulate the shear flow in a rectangular region. A commonly accepted idea is that the flow field can cause stress imbalance between different components of polymer mixtures and relative motion can thus take place. By an adequate formulation of polymer blends, and careful control of temperature and flow fields, it is possible to manipulate the morphology of materials in micromorphological engineering.

Phase separations in classical binary fluids have been extensively studied both experimentally and theoretically over several decades.¹⁶ Theoretical frameworks for modelling the effect of flow on phase transitions in polymer mixtures must incorporate both thermodynamic and viscoelastic forces into dynamic equations. Most theoretical research are based on the two fluid model, which combines thermodynamics, hydrodynamics, and viscoelasticity in a natural, physical way.

Few numerical results were presented for the two fluid model, especially regarding the effects of flow on phase diagram. Most of the early solutions of the two fluid model were obtained from approximations to the model, such as linearization of the diffusive equation and the adiabatic approximation,^17–19 in which time dependence of the viscoelastic stress tensor was not considered to seek a possible steady state solution. Some recent research^20–22 aim to numerically study the phase separation in a binary fluid, and large scale 2D and 3D simulations are carried out, but in absence of a realistic constitutive model, all of them can not accurately capture the viscoelastic effects. Modeling the phase separation in viscoelastic complex fluids with a realistic model remains elusive.

Jupp and Yuan^23–25 presented a two fluid model (FH–JS model) to investigate the flow-induced phase mixing and demixing of polymer mixtures, which coupled the Flory–Huggins mixing free energy function and the Johnson–Segalman constitutive model in a unified way. The numerical results were obtained by solving the time-dependent two fluid model without any simplification used widely in early studies. Fielding and Olmsted^26,27 proposed a different version of the two fluid model for the coupling between shear banding and shear-induced demixing, which still retained a diffusive stress term in the J–S constitutive equation. It is questionable whether or not such an ad hoc treatment for mechanical interface is justifiable as the origin of the diffusive stress is molecular diffusion across the streamlines and has already been considered in a thermodynamically consistent way in the original two fluid model.

Moreover, the J–S model is a phenomenological constitutive model. The Rolie–Poly model²⁸ is a simplified version of the full molecular model based on the tube dynamics,²⁹ yet still includes all dynamic modes considered in the full theory, including the reptation, chain stretch, contour length fluctuation, and thermal/convective constraint releases (CCR). Many numerical results^30–33 were carried out using a diffusive Rolie–Poly model, which verified that the Rolie–Poly constitutive model has the capability to accurately capture the rheological responses in a complex fluid. Recently, Cromer³⁴ predicted the steady-state shear-banding with Rolie–Poly constitutive model for a monodisperse polymer solution, which coupled the concentration and the stress through a two-fluid approach, but only one-dimensional calculations had been presented. To the best of our knowledge, this is the first attempt to couple the Rolie–Poly constitutive model with a two fluid model without adding a diffusive stress term. However, compared to Cromer's model, the framework of FH–JS model^23–25 is more general and can model not only the polymer solutions with a reduced formulation but also the binary polymer mixtures in which both of components are viscoelastic.

In this study, the FH–JS model is modified by replacing the J–S constitutive model with the Rolie–Poly model. It is not only inherent all the advantages of the original model, but also can accurately capture polymer dynamics in flow-induced phase separation. Simulation work focuses on the shear-induced phase transition of polymer mixtures in which the volume fractions of the two components are equal. Both of transient and steady-state phase transition behaviour and corresponding rheological responses will be presented. In some cases, rather than reaching a steady state, the numerical results show erratic fluctuations. This phenomenon, often named as rheochaos, has been observed in many complex fluids, including wormlike micelles,^35–37 associative polymer networks¹¹ and lamellar phase systems.^38,39 In order to capture the key features of chaotic flows, several models^40,41 have been studied by introducing some coupling parameters between flow and microstructure. Numerical results presented in this study demonstrate that the dynamic features of rheochaos can be captured by our model without introducing any extra parameters. The correlations between microstructure evolution and chaotic rheological responses are also identified.

The remainder of the paper is organized as follows. In the next section a detailed description of the two fluid model integrated with the Rolie–Poly constitutive model is presented. Detailed numerical algorithms used in this paper are described in Section 3; Simulation results and discussions are provided in Section 4, in which the dynamic phase separation and rheological behaviour for both monotonic and non-monotonic constitutive relationships are analyzed through numerical simulations implemented using finite volume method. Section 5 contains our conclusions.

2 A two-fluid model integrated with Rolie–Poly constitutive model

Here we focused on the strong flow regime, as in such regime the viscoelastic and hydrodynamic forces completely dominate over the forces originated from the thermal fluctuations, it is therefore justifiable to neglect the effects of the thermal fluctuations as we did in the previous work.^23–25

To establish a theoretical framework for describing the phase separations of polymer blends, we coupled the Flory–Huggins free energy function with the Rolie–Poly constitutive model. The two-fluid model provides a physical approach to capture the inhomogeneous polymer dynamics. In the standard framework of the two-fluid model, the system consist of two components: the polymer and the solvent. A constitutive equation is required to describe the viscoelastic stress of the polymer component, and the solvent is pure Newtonian.

Compared to the original two-fluid approach,^19,42,43 there are two main distinctions in our model. Firstly, both components are viscoelastic, and some mixing rules are introduced for calculation of the viscoelastic stress. This general formulation for binary polymer mixtures can be easily reduced to the standard two-fluid framework. Secondly, we replace the phenomenological constitutive model with a molecular-based model: Rolie–Poly model, which can provide more realistic predictions of the viscoelastic response under flow.

For polymer blends with two components, let ϕ_A(r⃑, t) and ϕ_B(r⃑, t) be the volume fractions of the component A and component B at some point r⃑ and time t, and let v⃑_A(r⃑, t) and v⃑_B(r⃑, t) be their velocities, respectively. The volume average velocity v⃑ is given by v⃑(r⃑, t) = ϕ_Av⃑_A(r⃑, t) + ϕ_Bv⃑_B(r⃑, t). The evolution equation for the volume fraction follows is given by the conservation law of mass, and may be written in terms of the volume average velocity:

The essential term of eqn (1) can be expressed by the velocity difference between the two components, which depends on the osmotic force, μ and the force due to inhomogeneity in polymer viscoelastic stress, ·σ_ve:

where ζ = ζ_Aζ_B/(ζ_A + ζ_B) is a frictional coefficient. Following the Doi and Onuki 's argument⁴² in their two fluid model for polymer mixtures, the friction coefficient may be expressed as ζ_i = ϕ_iζ_0iM_i/M_ei (i ∈ {A, B}, ζ_0i and M_ei are the microscopic friction constant and the entanglement molecular weight for component i, respectively). For simplicity, we assume ζ_0i = ζ₀.

For the polymer inhomogeneity, we introduce a dimensionless coefficient α, and different from the previous work, the ratio of relaxation times is replaced by the ratio of the entanglements number,

where the coefficient is related to the ratio of the modulus (G′ = G_B/G_A) and the ratio of the entanglements number (Z′ = Z_B/Z_A). Z_i = M_i/Mⁱ_e is the number of entanglements of each component polymer, and can be calculated by Rouse relaxation time τⁱ_R and disentanglement relaxation time τⁱ_d as Z_i = τⁱ_d/(3τⁱ_R). This coefficient vanishes for G′Z′ = 1, which includes the special case of a rheologically homogeneous system, where G′ = 1 and Z′ = 1. The viscoelastic force term, with the dimensionless coefficient α, can work either to oppose or to reinforce the thermodynamic force in the diffusive dynamics.

The osmotic force can be computed by the chemical potential difference μ, which is defined as the functional derivative of the mixing free energy with respect to local volume fraction, i.e.,

Here we take the first order approximation of Flory–Huggins free energy function form^44,45 as

[scr F, script letter F]_mix[ϕ_A(r⃑)]/k_BT = ∫dr⃑{f_mixk_BT + (Γ/2)[∇ϕ_A]²} (5)

where Γ is the interfacial tension coefficient and any dependence of Γ on ϕ_A is ignored, and the equilibrium phase behavior can be defined by the Flory–Huggins mixing free energy, f_mix/k_BT = (1/M_A)ϕ_A ln ϕ_A + (1/M_B)ϕ_B ln ϕ_B + χϕ_Aϕ_B, where M_i is the molecular weight of each component polymer and χ is the Flory–Huggins interaction parameter.

Apart from the osmotic force, another force the polymer experiences is the viscoelastic stress σ_ve(r⃑, t), which can be calculated by constitutive equations. In our model, we replace the original JS model with the Rolie–Poly model. As the stress in entangled polymer systems is supported by chain molecules, a stress gradient would create a net force on the chains resulting in their relative migration. Thus v⃑_T should be used in the viscoelastic constitutive equation rather than v⃑.

The tube velocity given by , may be expressed in terms of the volume average velocity by

where v⃑_A − v⃑_B is defined by eqn (2).

The relaxation modulus of a binary viscoelastic fluid must depend on its composition. For polymer linear blend, the relaxation modulus may be described as G(t) = G_R(t) + G_d(t), where G_R(t) describes the high-frequency polymer relaxation (Rouse process) and G_d(t) describes the disentanglement process of the chains. For the Rouse process, the linear mixing rule⁴⁶ is adopted as,

The ‘double reptation’ rule offers much better agreement with experimental data than the ‘linear mixing rule’ for disentanglement process, hence a quadratic mixing rule is suggested in the form^47,48 of

where gⁱ_R(t) and gⁱ_d(t) are the linear-viscoelastic relaxation moduli of different components, which may be expressed by the exponential functions of the elastic modulus (Gⁱ_R, Gⁱ_d) and the relaxation times (τⁱ_R, τⁱ_d).⁴⁸ After this treatment, the viscoelastic stress σ_ve(r⃑, t) is calculated from the sum of three parts, the contribution of the component A and B, and a mixing term representing the interaction of both components:

In this representation the quantity W_i is the viscoelastic strain and stress G_iW_i is parameterized by the elastic modulus G_i. According to the mixing rules, the elastic modulus can be expressed as G_AA = ϕ_AG^A_R + ϕ_A²G^A_d, G_BB = ϕ_BG^B_R + ϕ_B²G^B_d and . Accordingly the relaxation times can be written as τ^AA_d = τ^A_d, τ^BB_d = τ^B_d and τ^AB_d = 2τ^A_dτ^B_d/(τ^A_d + τ^B_d). For the Rouse time, there is no interactive part, so τ^AA_R = τ^A_R, τ^BB_R = τ^B_R and τ^AB_R = 0.

After applying the mixing rules to the Rolie–Poly constitutive model, the evolution equation of viscoelastic strain W_i may be written as,³³

where and the material derivative, , is defined in terms of the tube velocity as ; the coefficients involving the trace of the strain tensor are defined byfor W_AB, as τ^AB_R → 0 leading to a non-stretching limit,²⁸ the expressions are given by

f^AB₂ = β (14)

We assume, for simplicity, that the components are of equal density, ρ. The motion of an incompressible, isothermal viscoelastic binary fluid is governed by the continuity equation

and the osmotic pressure and the viscoelastic stress are brought into the momentum balance equation,where p is the pressure field, and the Newtonian stress with viscosity η₀ can be a contribution from the solvent or from energy dissipation processes much faster than the slowest viscoelastic dynamics governed by the time-dependent viscoelastic stress tensor σ_ve(r⃑, t).

The equations above complete the Flory–Huggins–Rolie–Poly (FH–RP) fluid model in this study. The parameter β is the CCR magnitude coefficient and δ is an exponent to quantify the CCR. For simplicity, but without much loss of generality, the same value for β and δ for the three parts of σ_ve is assumed. By exploring the parameter space of the FH–RP model in terms of η₀, τ^A_R, τ^A_d, τ^B_R, τ^B_d, G^A_R, G^B_R, G^A_d, G^B_d (G_i = Gⁱ_R + Gⁱ_d) and shear rate, the details of nonlinear dynamic behaviour in shear-induced phase transition of polymer mixtures could be systematically studied.

3 Numerical methods

The governing equations of the FH–RP model are discretized through finite volume method (FVM) using an open source OpenFOAM CFD toolbox released by OpenCFD Ltd. The finite volume method ensures that the physical conservation laws are locally satisfied, through computing each term of the governing equations by volume integral over a control volume. For those terms involving divergence operators, volume integral can be converted to surface integral over the surfaces of the control volume by applying the Gauss theorem. The resulting expressions are then discretised using appropriate differencing schemes. Gauss MINMOD and Gauss linear are applied in the discretisation scheme for the spatial discretization terms, and the temporal term is discretised using the simple Euler scheme. The above discretisation produces a system of linear equations to relate the cell face values of the unknown field variables at current time to their values at previous time and at the neighbouring cell centres. The system can be solved by one of the iterative matrix solvers, including the conjugate (PCG) and bi-conjugate gradient (PBiCG) methods, available in the OpenFOAM toolbox. The details on how to use it can be found in the OpenFOAM manual.⁴⁹

The iterative solution algorithm used in this paper is a PISO-based algorithm, which have been well tested in a numerical study for the dynamics of polymer solutions in contraction flow.⁵⁰ The numerical method involves a modification to the standard PISO algorithm by introducing polymer stress and the chemical potential unknowns into the momentum equation, as a part of the source term. In the absence flow, our simulation could reproduce the equilibrium phase behavior defined by the Flory–Huggins mixing free energy functional. Also the simulation results of homogeneous fluids under shear flow converge to the solution of the constitutive model with satisfactory accuracy. Moreover a similar OpenFOAM based viscoelastic flow solver has been critically validated under benchmark flow problems.⁵⁰

We summarised the key procedure of the iterative algorithm as follows

Step 1: assuming νⁿ⁺¹ = νⁿ, integrate the constitutive equations implicitly to get the polymer viscoelastic stress components, σ_ve(r⃑, t)ⁿ⁺¹.

Step 2: compute the force field μ and ·σ_ve, then solve eqn (1) to get ϕⁿ⁺¹_A.

Step 3: solve the momentum equation without the contribution of the pressure gradient term to obtain the estimated velocity components Uⁿ⁺¹.

Step 4: use the estimated velocity Uⁿ⁺¹ to solve the pressure-correction equation of PISO algorithm, get the pressure field pⁿ⁺¹.

Step 5: use the estimated pressure field pⁿ⁺¹ to calculate the corrected velocity field νⁿ⁺¹.

Step 6: use the corrected velocity field νⁿ⁺¹ to solve the constitutive equation and the composition field evolving equation (eqn (1)) again to get the corrected stress field and the corrected composition field.

Step 7: repeat the steps 3–6, using the corrected variable νⁿ⁺¹, pⁿ⁺¹, σ_ve(r⃑, t)ⁿ⁺¹ and ϕⁿ⁺¹_A as improved estimates until all corrections are negligibly small for the solutions at the present time.

Step 8: repeat from the step 1 and advance to the next time step.

Here the superscript n and n + 1 represent the past and the present times, respectively. The numerical scheme solves the non-linear governing equations in time and it is capable of modeling the unsteady flow.

4 Results and discussions

The maximum lattice size of two-dimensional simulation box in square shape is set to 256 × 256. The top and bottom boundaries are physical walls and periodic boundary condition is applied to the other two sides of the box. Care has been taken in all simulations to ensure that there are sufficient lattice sites (seven at least) across the sharpest interfaces and also that the simulation results are independent of any further refinement of lattice density. Note also that these simulations cannot reproduce any physical structure with length scale larger than the simulation box.

For simplicity but without much loss of generality, the parameters are set as: M_A = M_B = 1, k_BT = 1.3, Γ = 1, ζ = 0.1, η₀ = 0.1. Thus the equilibrium phase diagram is always symmetric with respect to ϕ_A = 0.5. The critical point is at χ = 2. Therefore, χ = 1.0 for the FH–RP fluid means far away from the phase-coexistent region at equilibrium and the fluid is homogeneous for the entire composition range at this temperature.

A Gaussian random noise with an intensity of 10⁻³ was superimposed on a specified initial uniform composition field ϕ⁰_A = 0.5 in all simulations reported here. Each component of the viscoelastic stress tensor σ_ve, including σ_i, = 1, 2, 3, are initially set to zero. A wall velocity of U₀ = [small gamma, Greek, dot above](L/2) is applied to the top and bottom walls in equal speed and opposite direction to achieve a required shear rate [small gamma, Greek, dot above], where L is characteristic length (distance between the walls).

By choosing appropriate values for the parameters β, δ, G_i, τⁱ_d and τⁱ_R, the FH–RP model exhibits different constitutive behaviour. Both of the monotonic and non-monotonic constitutive relationships will be analyzed in the following sections.

4.1 Interface instabilities with a monotonic constitutive equation

The CCR magnitude coefficient β is a control parameter on the monotonicity of the constitutive equation and can be fitted from experimental data. In comparison with the prediction of the full molecular theory for polymer melts, Likhtman and Graham²⁸ used β = 0.5 for the transient data and β = 1 for the steady state data. However in order to obtain a good fit of the experimental data, they suggest that a small value for the CCR parameter β should be used. Adams and Olmsted tune the parameter between 0.5 and 1 to generate either non-monotonic or monotonic constitutive behaviour.^30,33 They usually set , which is much higher than the values needed to observe phase separation experimentally, and β > 0.5 to obtain constitutive curve with a broad plateau.

Kabanemi and Chung^31,32 set a much lower value of β, and in some cases β = 0 to produce non-monotonic constitutive curve. In general, the value of β should be determined from experimental data, hence allow the FH–RP model to capture various physical features of complex fluids.

The simulation results presented in this section are obtained by setting the model parameters β = 0.2 and δ = −0.5, G_A = 1.0 (G^A_R = 0.1, G^A_d = 0.9), G_B = 4.0 (G^B_R = 0.2, G^B_d = 3.8), τ^A_R = 1.0, τ^B_R = 0.01, , and (i.e., Z_A = 6.67, Z_B = 66.67). Those parameters give a strictly monotonic constitutive equation for all values of ϕ_A, and define a binary polymeric system which is miscible under equilibrium condition and shows shear-induce phase separation in a range of shear rates. The steady state solutions of the model are calculated by solving eqn (9) to eqn (14) under a fixed composition field ϕ_A using the pre-defined ODE solvers available in MATLIB and are shown in Fig. 1. The absolute differences between the numerical results obtained from the OpenFOAM solver and the MATLAB calculation for a homogenous Rolie–Poly fluid are less than 0.003.

Fig. 1 The scaled steady-state shear stress T_xy versus [scr W, script letter W]_e = [small gamma, Greek, dot above]τ^A_d calculated by the homogeneous FH–RP model of various ϕ_A, as indicated, with parameters β = 0.2, G_A = 1.0, G_B = 4.0, and .

With this set of the parameters, all the constitutive curves of homogenous fluid under various ϕ_A are monotonic as indicated in Fig. 1. There is no intrinsic constitutive instability. Under low shear rate ([scr W, script letter W]_e ≲ 6.0), the calculated shear stress and first normal stress difference exactly match the corresponding analytical results of homogeneous fluid (ϕ_A = 0.5). The phase separation only occurs in a range of shear rates: 7.0 ≲ [scr W, script letter W]_e ≲ 38.0, in which the homogeneous phase with initial small composition fluctuation separates into A-component-rich and B-component-rich phases in a form of shear bands. In the higher shear rate region of [scr W, script letter W]_e ≳ 40.0, the simulation results again coincide with the analytical calculations for a homogeneous fluid and no phase separation is observed over the duration of our simulations. Fig. 1 shows that no stress plateau exists. It was also reported in the previous study of the FH–JS model^24,25 and a liquid crystalline model.⁵¹

The transient behaviour corresponding to several points in Fig. 1 are presented in Fig. 2(a) for shear stress and Fig. 2(b) for first normal stress difference. Within the phase separation region, the numerical results of the shear stress and the first normal difference coincide with the theoretical solution of a homogeneous FH–RP fluid but then turn into a two-phase solution as a significant drop (the second drop in shear stress) of T_xy and N₁ appears. It takes as long as 70τ_A to reach a steady state.

Fig. 2 Transient shear stress and first normal stress difference for various values of the Weissenberg number [scr W, script letter W]_e with parameters: β = 0.2, G_A = 1.0, G_B = 4.0, and .

Snapshots of the composition field at several instant during phase separation for the case of [scr W, script letter W]_e = 20.0 have been shown in Fig. 4. The evolution of the shear rate field and the stress field follow the same pattern as that of the composition field shown in the figure. The letters a to j mark the points along the transient stress curves in Fig. 2 which refer to the individual images in Fig. 4, respectively. From the evolution of composition field, there are two dynamical mechanisms by which the shear bands reach to a steady state: diffusive mechanism and convective through wave instability. They have also been observed in a 30/70 polymer mixture,²⁴ but only a purely diffusive mechanism was seen in 50/50 polymer mixtures.²⁵ At t = 0, the composition is homogeneous. The significantly stress dropping between points b and d in Fig. 2, corresponding to Fig. 4(b)–(d), shows the direct formation of fine composition bands by nucleation and unstable wave structure interacting with each other. At point d, the stress has reached a low plateau but the composition field is still unstable. As approaching to the steady state, there is an interesting dynamical process observed between point e and g referring to the snapshots shown in Fig. 4(e)–(g): a wave structure breaks up forming droplet, then the droplet interacts with the remaining bands and subsequently forms into a coarser band. This phenomenon explains the stress fluctuations in Fig. 2; It is not only seen in this specific case for [scr W, script letter W]_e = 20.0, but also for all the cases shown in Fig. 2, fluctuations always appear along with a process breaking wave into droplet then to much coarser band. From t = 30τ_A, the band structures gradually grow via diffusive mechanism, and shortly after the small oscillation around t = 85τ_A corresponding to another wave to droplet process, the shear bands get much coarser and far away from each other. From t = 100τ_A, the remaining bands change very little and reach the steady state.

Fig. 3 The composition profiles across the simulation box for a binary fluid of [scr W, script letter W]_e = 20.0 with β = 0.2, G_A = 1.0, G_B = 4.0, and . (Top) snapshots of composition along the shear gradient (y) direction; (bottom) snapshots of composition along the flow direction (x) near the band interface (y = 64) in the steady state.

Fig. 4 Snapshots of the composition field ϕ_A at various instant t for an initially homogeneous viscoelastic binary fluid of [scr W, script letter W]_e = 20.0, β = 0.2, G_A = 1.0, G_B = 4.0, and .

It is worth to note that our numerical results show a dynamic steady state, unlike the previous studies presenting a strictly smooth shear bands with J–S model,^24,25 there are wavy band structures moving with the flow without further coarsening. Fig. 3 shows the composition profiles across the simulation box. As seen in the top figure, the fluid formed into band structures and the maximum value of ϕ_A is less than 0.8 in steady state. In order to reveal structural evolution of the unstable interfaces, we plot the composition profile along the flow direction near the band interfaces, as shown in the bottom figure of Fig. 3, which presents a typical moving pattern: a wavy structure moves along the flow direction.

Recently, a growing body of data reveals that the shear bands can fluctuate. Most of the observations suggest a wavy structure oscillating between two bands. Many experimental results identify such fluctuations in complex fluids^9,52–57 and the interface instabilities have also been observed in numerical studies with a diffusive J–S model^27,55,58 and a two-species model.⁵⁹ The simulation results demonstrate that the proposed FH–RP model can capture the essential instabilities in complex fluids through a thermodynamically consistent way.

4.2 Chaotic rheological responses with a non-monotonic constitutive equation

In this section, we numerically study the dynamic phase separations in polymer mixtures with a typical non-monotonic constitutive equation. To generate a non-monotonic constitutive curve with a wide range of ϕ_A, the model parameters are set as: β = 0 and δ = −0.5, G_A = G_B = 2.5, τ^A_R = 0.2, τ^B_R = 0.01, and . This defines a binary polymeric system with a non-monotonic constitutive behaviour. It is miscible under equilibrium condition.

The intrinsic constitutive curves of various values of ϕ_A have been presented in Fig. 5, all of which are non-monotonic. The numerical simulation results of the shear stress have also been plotted in this figure. In the homogeneous range for low shear rate ([scr W, script letter W]_e ≲ 7.0) and high shear rate ([scr W, script letter W]_e ≳ 80), the simulation results exactly coincide with the FH–RP constitutive curve in a homogenous state, which are indicated by the circles in Fig. 5. Similar to the monotonic case, in the region of 8 ≲ [scr W, script letter W]_e ≲ 77, flow-induced phase separation takes place and the stress deviation is observed. Different from previous cases in which the shear stress always can reach a steady state, in this system the stress does not reach a steady value, but continuously fluctuates. The fluctuation range of the shear stress is indicated with the vertical error bars in Fig. 5. Another significant difference is that the flow curve of shear stress versus shear rate is not linear as the monotonic case, but exhibits some wavy structure within the phase separation region. This may attribute to the large fluctuations in the simulations. In order to check whether or not the final state achieved here, it is necessary to use a much larger simulation box and this will form part of our future research. However, in an experimental study of associative polymer networks, Sprakel¹¹ presented similar fluctuations of the flow curves and their work revealed rheochaos of the transient stress response, which is similar to our simulation results.

Fig. 5 (Vertical bars) fluctuation range of T_xy versus scaled shear rate [scr W, script letter W]_e = [small gamma, Greek, dot above]τ^A_d for shear-induced phase separation with parameters β = 0, G_A = G_B = 2.5, and , along with the corresponding constitutive curves for homogeneous FH–RP fluids of various values of ϕ_A, as indicated.

The transient behaviour of the stress components at various shear rates corresponding to several points in Fig. 5 are shown in Fig. 6(a) for shear stress T_xy and Fig. 6(a) for first normal stress difference N₁. Under low shear rate ([scr W, script letter W]_e = 10), the magnitude of the fluctuation is small. In the middle of the phase separation region, such as [scr W, script letter W]_e = 20 and [scr W, script letter W]_e = 30, the fluctuations become much larger. In even higher shear rate region, the stress oscillation reduces as shown in the case of [scr W, script letter W]_e = 40. If looking into the detailed evolutions of the transient behaviour of the stresses, it reveals some special patterns. For [scr W, script letter W]_e = 40, the stress continuously oscillates with the same amplitude and periodically repeats similar patterns. Surprisingly, the case of [scr W, script letter W]_e = 20 shows a quite difference behaviour. Except for the larger fluctuations, there is a quiet region in which the stress oscillation vanishes.

Fig. 6 Transient shear stress and first normal stress difference revealing rheochaos for various values of the Weissenberg number [scr W, script letter W]_e with parameters: β = 0, G_A = G_B = 2.5, and .

Snapshots of the composition field during phase separation for the case of [scr W, script letter W]_e = 20 with a non-monotonic constitutive curve are shown in Fig. 7. The letters a–j marked in Fig. 6 refer to the individual images in Fig. 7 and indicate the points along the transient stress curves corresponding to each image. At t = 0, the composition is homogeneous with a small random noise. After the initial stress overshoot, at t = 0.95τ^A_d, the composition field forms smooth band structures. Shortly it is found that the flat bands have been broken as shown in Fig. 7(c). From here, the chaotic stress response emerges, and typical phase transitions corresponding to the rheochaos can be seen in Fig. 7(d)–(g). At the beginning, the bands stay with unstable interfaces, then a vortex structure emerges within the central band, and this turbulence gradually subsides at t = 20.7τ^A_d. We plot the snapshots of streamlines in Fig. 8 between t = 17.16τ^A_d and t = 20.7τ^A_d, from which we can find a symmetric vortex clearly appearing within the central band. This process repeat for a long period and lead to large fluctuations in stress curves. Although the results in an experimental study of telechelic polymers¹¹ suggest that the fluctuating stress is caused by interface instabilities, our simulations reveal another possible cause: the vortex instability inside the band structure. We hope more visualization experiments could be conducted to confirm our observations in the future.

Fig. 7 Snapshots of the composition field ϕ_A at various times t for an initially homogeneous viscoelastic binary fluid of [scr W, script letter W]_e = 20.0, with parameters: β = 0, G_A = G_B = 2.5, and

Fig. 8 Snapshots of streamlines corresponding to typical rheochaos, which clearly show a vortex structure emerges for an initially homogeneous viscoelastic binary fluid of [scr W, script letter W]_e = 20.0

At t = 61.36τ^A_d, corresponding to a smooth region of the stress curves showed in Fig. 6, the band structures also form a more stable state in Fig. 7(h): the central band becomes uniform with a smooth interface, and the fine bands near the wall combine into darker ones after interacting with each other. This smooth region terminates at around t = 70.95τ^A_d, thus the stress curves come back to the fluctuations. The composition field has been shown in Fig. 7(i) and (j). It is seen that instabilities emerge within the central band, and it is still unstable by t = 100τ^A_d. This particular intermittent rheochaos have not been observed in the numerical studies with a FH–JS model.

The visual impression of this vortex structure shown in Fig. 7 and 8, are observed in all of the cases with large rheological fluctuations. A more careful analysis is presented in Fig. 9. It shows the average Fourier spectra of the intensity reflected by composition field along the velocity direction in the central of the simulation box, which exhibits power law decay over one decade of frequency. Similar to the case of elastic turbulence,^60,61 the power spectrum in wave vector seems to approximately follow a power law up to k_d ∼ 33(Δx)⁻¹, where is the lattice size in this study. The case of [scr W, script letter W]_e = 10 shows significant difference from other cases in Fig. 9, as indicated in the figure, the critical value k_d′ < k_d. Compared to other cases, with [scr W, script letter W]_e = 10, we find that the stress fluctuation is relatively small and the shear bands are more stable, and this also leads to a much smaller vortex structure in the streamline plot. Since k⁻¹_d is related to the stress correlation length, we may conclude that larger vortex and rheological fluctuations are relevant to a smaller stress correlation length.

Fig. 9 Average Fourier spectra for various [scr W, script letter W]_e obtained on the wave vector domain by computing Fourier transforms of the intensity reflected by the composition field along the velocity direction at the middle point, and then by averaging for different times. The two solid lines are k^−6.0 and k^−1.2 fit respectively.

5 Conclusions

The two fluid model directly couples the viscoelastic stress and the diffusive fluid composition, thus there is no need to introduce any diffusive stress term into the constitutive equations. In this study, we extend the FH–JS model by integrating with Rolie–Poly constitutive model. Two-dimensional simulations have been carried out and the numerical results for a typical monotonic and non-monotonic constitutive equation are analyzed respectively, which show that our model can capture the essential features of the shear-induced phase transitions in binary polymer mixtures.

The simulations reproduce many dynamic phenomena reported in literature, including the steady-state shear-banding, the interface instabilities between the bands, and the dynamic mechanisms to reach stable band structures. With a non-monotonic constitutive equation, we have also observed that the band structures are strongly unstable both in time and in space. The stress continuously fluctuates in a chaotic manner. The correlations between microstructure evolution and chaotic rheological responses have been identified. Other than the interface instability, our results reveal another possible cause of the rheochaos: a vortex structure emerges within the central band.

Coupled with a realistic constitutive model, the computational model presented here could provide valuable guidance to new experimental work and could also be used to quantitatively study some outstanding problems by appropriately selecting the parameters to mimic the real physical systems. Motivated by many experimental observations, we hope to address an in-depth quantitative study of the rheochaos in complex fluids in future simulations, including large-scale three-dimensional simulations.

Acknowledgements

This work is partially supported by BBSRC grant (BB/K011146/1) in UK, the National Natural Science Foundation of China (Grant no. 61221491, no. 61303071 and no. 61120106005), and Open fund from State Key Laboratory of High Performance Computing (no. 201303-01).

References

1. A. Silberberg and W. Kuhn, Nature, 1952, 170, 450–451.
2. D. Jou, J. Casasvazquez and M. Criadosancho, Phys. Prop. Polym., 1995, 120, 207–266.
3. A. Onuki, J. Phys.: Condens. Matter, 1997, 9, 6119.
4. M. Laurati, G. Petekidis, N. Koumakis, F. Cardinaux, A. B. Schofield, J. M. Brader, M. Fuchs and S. U. Egelhaaf, J. Chem. Phys., 2009, 130, 134907.
5. S. Lerouge and J. F. Berret, Polymer Characterization: Rheology, Laser Interferometry, Electrooptics, 2010, vol. 230, pp. 1–71.
6. B. Lonetti, M. Camargo, J. Stellbrink, C. N. Likos, E. Zaccarelli, L. Willner, P. Lindner and D. Richter, Phys. Rev. Lett., 2011, 106, 228301.
7. S. Ravindranath, S. Q. Wang, M. Olechnowicz, V. S. Chavan and R. P. Quirk, Rheol. Acta, 2011, 50, 97–105.
8. J. Y. Lee, G. G. Fuller, N. E. Hudson and X.-F. Yuan, J. Rheol., 2005, 49, 537–550.
9. R. Angelico, C. O. Rossi, L. Ambrosone, G. Palazzo, K. Mortensen and U. Olsson, Phys. Chem. Chem. Phys., 2010, 12, 8856–8862.
10. S. Manneville, A. Colin, G. Waton and F. Schosseler, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 061502.
11. J. Sprakel, E. Spruijt, M. A. C. Stuart, N. A. M. Besseling, M. P. Lettinga and J. van der Gucht, Soft Matter, 2008, 4, 1696–1705.
12. T. Divoux, D. Tamarii, C. Barentin and S. Manneville, Phys. Rev. Lett., 2010, 104, 208301.
13. P. E. Boukany and S. Q. Wang, J. Rheol., 2009, 53, 73–83.
14. P. E. Boukany, Y. T. Hu and S. Q. Wang, Macromolecules, 2008, 41, 2644–2650.
15. P. E. Boukany and S. Q. Wang, Macromolecules, 2010, 43, 6950–6952.
16. C. Domb, Phase transitions and critical phenomena, Access Online via Elsevier, 2000, vol. 19.
17. H. Ji and E. Helfand, Macromolecules, 1995, 28, 3869–3880.
18. N. Clarke and T. C. B. McLeish, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1998, 57, R3731–R3734.
19. S. T. Milner, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 3674–3691.
20. P. Stansell, K. Stratford, J. C. Desplat, R. Adhikari and M. E. Cates, Phys. Rev. Lett., 2006, 96, 085701.
21. K. Stratford, J.-C. Desplat, P. Stansell and M. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 030501.
22. S. M. Fielding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 021504.
23. X. F. Yuan and L. Jupp, Europhys. Lett., 2002, 60, 691–697.
24. L. Jupp and X. F. Yuan, J. Non-Newtonian Fluid Mech., 2004, 124, 93–101.
25. L. Jupp, T. Kawakatsu and X. F. Yuan, J. Chem. Phys., 2003, 119, 6361–6372.
26. S. M. Fielding and P. D. Olmsted, Phys. Rev. Lett., 2003, 90, 224501.
27. S. M. Fielding and P. D. Olmsted, Phys. Rev. Lett., 2006, 96, 104502.
28. A. E. Likhtman and R. S. Graham, J. Non-Newtonian Fluid Mech., 2003, 114, 1–12.
29. R. S. Graham, A. E. Likhtman and T. C. B. McLeish, J. Rheol., 2003, 47, 1171–1200.
30. J. M. Adams and P. D. Olmsted, Phys. Rev. Lett., 2009, 102, 067801.
31. K. K. Kabanemi and J. F. Hetu, Rheol. Acta, 2009, 48, 801–813.
32. C. Chung, T. Uneyama, Y. Masubuchi and H. Watanabe, Rheol. Acta, 2011, 50, 753–766.
33. J. M. Adams, S. M. Fielding and P. D. Olmsted, J. Rheol., 2011, 55, 1007–1032.
34. M. Cromer, M. C. Villet, G. H. Fredrickson and L. G. Leal, Phys. Fluids, 2013, 25, 051703.
35. R. Bandyopadhyay, G. Basappa and A. K. Sood, Phys. Rev. Lett., 2000, 84, 2022–2025.
36. L. Becu, S. Manneville and A. Colin, Phys. Rev. Lett., 2004, 93, 018301.
37. R. Ganapathy and A. K. Sood, Phys. Rev. Lett., 2006, 96, 108301.
38. S. Manneville, J. B. Salmon and A. Colin, Eur. Phys. J. E: Soft Matter Biol. Phys., 2004, 13, 197–212.
39. L. Gentile, B. F. B. Silva, S. Lages, K. Mortensen, J. Kohlbrecher and U. Olsson, Soft Matter, 2013, 9, 1133–1140.
40. B. Chakrabarti, M. Das, C. Dasgupta, S. Ramaswamy and A. K. Sood, Phys. Rev. Lett., 2004, 92, 055501.
41. S. M. Fielding and P. D. Olmsted, Phys. Rev. Lett., 2004, 92, 084502.
42. M. Doi and A. Onuki, J. Phys. II, 1992, 2, 1631–1656.
43. H. Ji and E. Helfand, Macromolecules, 1995, 28, 3869–3880.
44. H. Tanaka, Phys. Rev. Lett., 1993, 71, 3158–3161.
45. B. M. Forrest and R. Toral, J. Stat. Phys., 1994, 77, 473–489.
46. E. van Ruymbeke, J. Nielsen and O. Hassager, J. Rheol., 2010, 54, 1155–1172.
47. J. Descloizeaux, Europhys. Lett., 1988, 5, 437–442.
48. E. van Ruymbeke, C.-Y. Liu and C. Bailly, Rheol. Rev., 2007, 53–134.
49. OpenFOAM programmer's guide, OpenFOAM programmer's guide (Version 2.1.0), 2013.
50. S. C. Omowunmi and X. F. Yuan, Rheol. Acta, 2013, 52, 337–354.
51. P. D. Olmsted and C. Y. D. Lu, Faraday Discuss., 1999, 112, 183–194.
52. M. R. Lopez-Gonzalez, W. M. Holmes, P. T. Callaghan and P. J. Photinos, Phys. Rev. Lett., 2004, 93, 268302.
53. S. Lerouge, M. Argentina and J. P. Decruppe, Phys. Rev. Lett., 2006, 96, 088301.
54. J. P. Decruppe, L. Becu, O. Greffier and N. Fazel, Phys. Rev. Lett., 2010, 105, 258301.
55. P. Nghe, S. M. Fielding, P. Tabeling and A. Ajdari, Phys. Rev. Lett., 2010, 104, 248303.
56. M. A. Fardin and S. Lerouge, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 178303.
57. M. A. Fardin, T. J. Ober, V. Grenard, T. Divoux, S. Manneville, G. H. McKinley and S. Lerouge, Soft Matter, 2012, 8, 10072–10089.
58. S. M. Fielding and H. J. Wilson, J. Non-Newtonian Fluid Mech., 2010, 165, 196–202.
59. M. Cromer, L. P. Cook and G. H. McKinley, J. Non-Newtonian Fluid Mech., 2011, 166, 566–577.
60. A. Groisman and V. Steinberg, Nature, 2000, 405, 53–55.
61. M. A. Fardin, D. Lopez, J. Croso, G. Gregoire, O. Cardoso, G. H. McKinley and S. Lerouge, Phys. Rev. Lett., 2010, 104, 178303.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra08448a

---