Analysis of a viscoelastic phase separation model

A new model for viscoelastic phase separation is proposed, based on a systematically derived conservative two-fluid model. Dissipative effects are included by phenomenological viscoelastic terms. By construction, the model is consistent with the second law of thermodynamics, and we study well-posedness of the model, i.e., existence of weak solutions, a weak-strong uniqueness principle, and stability with respect to perturbations, which are proven by means of relative energy estimates. A good qualitative agreement with mesoscopic simulations is observed in numerical tests.


Introduction
In this work we propose and analyze a new viscoelastic phase separation model.
The mixing properties of polymer solutions are highly temperature-dependent, leading to the terms "good" and "poor" solvent.
While good solvent corresponds to temperatures at which contacts between solvent particles and macromolecules are more favorable than self-contacts, for a poor solvent the opposite is the case, i.e.
self-contacts are more favorable. Of principle interest is to understand the dynamics of demixing after a quench from good to poor solvent. The relevant processes and dynamical equations are well-known for Newtonian binary fluids, and typically described in terms of the so-called "model H" [1][2][3]. However, in polymer solutions several new effects appear which are related to the different length and time scales of the solvent and the polymer molecules, leading to so-called dynamical asymmetry, see [4].
The origin of our investigations goes back to the work of Tanaka [5], who studied dynamic asymmetry theoretically and experimentally, see also [6,7] for three-dimensional simulations. The apparent inconsistency of Tanaka's model with the second law of thermodynamics was addressed in [8] and cured by modifications based upon applying basic principles of non-equilibrium thermodynamics [9][10][11]. In this work, a closure of the system was obtained by relating the relative velocity between solvent and polymer through an ad-hoc constitutive relation and additional equations for the viscoelastic stresses. For further work on viscoelastic phase separation we refer the reader to Refs. [12][13][14][15][16], while more recent studies including applications can be found in [17][18][19][20][21][22][23].
In Section 2, we discuss an alternative derivation starting from a conservative model for a binary fluid which can be obtained by systematic coarse graining from a microscopic description [24]. Viscous effects are introduced via additional stresses of phenomenological type, and some standard simplifications lead to a simplified version of the model in [8].
The systematic derivation allows us to describe in detail the assumptions underlying the model and to justify the closure relations used in [8]. Motivated by the considerations in [25], we include additional dissipative terms in the momentum equations and utilize a change of viscoelastic variables to obtain our final model. These deviations from [8] turn out to be essential for the verification of mathematical well-posedness.
In Section 3 we discuss the consistency with the second law of thermodynamics, and the existence of weak solutions. Based on relative energy estimates, we study conditional uniqueness and stability of solutions with respect to perturbations in the problem data and the model parameters.
In Section 4 we present a comparison of structure factors obtained by simulations of the new model with those of a mesoscopic description. A good qualitative agreement is observed in these numerical tests, which may serve as a first validation. We also refer to [26,27] where results of numerical simulations in two and three space dimensions are compared with physical experiments from [5] and indicate very good agreement.

Model derivation
We start from the following special case of a binary fluid model proposed in [24], namelẏ Here c, v, u, p denote, respectively, the mixture density difference (−1 ≤ c ≤ 1), the mass-averaged velocity, the relative velocity, and the pressure, which here is the Lagrange multiplier associated to incompressibility of v. Furthermore, f is the Helmholtz free energy density of mixing and ∂f ∂c the corresponding chemical potential. As usual, we denote byġ = ∂ t g+(v·∇)g the convective derivative with respect to the mass-averaged velocity v.
We assume in the following that Ω ⊂ R d is a bounded domain and require, for simplicity of presentation, appropriate periodic boundary conditions for all field variables. Then with E kin = 1 2 |v| 2 + 1−c 2 8 |u| 2 denoting the kinetic energy density. From a thermodynamic point of view, model (1) is thus fully conservative, as it is appropriate for isothermal conditions. Note also that we consider the overall density to be unity.
In a phenomenological manner, we now add viscoelastic effects by introducing corresponding stress tensors in the two momentum equations. Thus the modified momentum equations reaḋ Here T v and T u denote the tensors related to the trace-free parts of the (symmetrized) velocity gradients of v and u, respectively, while T trace u = qI is the tensor related to the trace of the velocity gradient of u. The tensor related to the trace part of the velocity gradient of v can be neglected, since we assumed v to be incompressible. Note that the momentum balance requires a divergence form in the equation for v. The divergence form of the viscoelastic variables in the equation for u is postulated.
Constitutive relations for the viscoelastic variables are required to close the system.
As a next step, we introduce some simplifying assumptions: Since the relative velocity u is a fast non-hydrodynamic variable, we assume that u is quasistatic, and replace ∂ t u + (v · ∇)u + (u · ∇)v by a relaxation term 1 τ u. Furthermore, we assume u to be small and neglect terms of higher order in u. This leads to the simplified momentum equationṡ As a consequence of the simplifying assumptions, the kinetic energy will now only depend on v.
As a next step, we neglect the stress tensor T u and, similar to [8,25], we choose the dissipative Oldroyd-B-type model to describe the time evolution of q and T v . In that case, the evolution of a stress tensor T is described by For the special case of a tensor T = qI, one simply obtains, due to incompressibility,q By combination of the previous considerations and adding additional contributions corresponding to the linear stresses Adiv (u) and A 2 D S v induced by the associated velocities, see [8], where D S v denotes the symmetrized velocity gradient tensor, we arrive at the following intermediate model Following [8] and common practice in the rheological literature, we assume that T v is positive semidefinite, i.e. has non-negative trace. Similarly, again following [8], we postulate the free energy whose time evolution is found, by straightforward calculation, to be By setting A = 1 and n(c) = 4(1 − c 2 ) −1 , the last term becomes a square, such that we obtain which shows that, under the above assumptions, the resulting model is consistent with the second law of thermodynamics. After eliminating the relative velocity u, we arrive at the following systeṁ Let us note that the first and third equation form a cross-diffusion system for c and q, and that such systems are well-suited to describe pattern formation. We further observe that the structure of the reduced model (7) is exactly the same as that in [8].
In fact, the model of this reference can be obtained by the following minor modifications: (i) change of variable c to volume fraction φ; (ii) Newtonian part in the viscous stress in (7) 2 ; (iii) inclusion of the Korteweg interface stress for φ in (7) 2 ; (iv) extra stress q replaced by A(φ)q; (v) special choice of τ dependent on φ; (vi) ε i = 0 for i = 1, 2.
The first two points are straightforward and we omit a thorough discussion. Point three can be incorporated by extending the potential f to include a gradient term modelling the interface stress, leading to the typical Ginzburg-Landau form. Point four is based on the idea that the bulk stress qI should change rapidly when φ transits between the pure concentrations. In fact, [5] assumed that for the solvent A should vanish. In fact, the whole derivation above can be done by replacing q by A(φ)q in the relation of the relative velocity, i.e. (4). Point five can be motivated on a similar reasoning. Let us finally comment on the last modification: In most of the rheological literature, non-diffusive equations with ε i = 0 are considered. A justification for ε i > 0 has been given in [25], where it is argued that such terms are related to the centerof-mass diffusion of the polymer chains. While the parameters ε i may be rather small in general, their positive value leads to much stronger results concerning analysis and numerical treatment.
Including the modifications (i)-(v) will lead to the model of [8], and our derivations so far may be considered as an alternative derivation of their model. On the basis of the same arguments, we will further replace the Oldroyd-B model for the elastic stress tensor T v by the Peterlin model for the conformation tensor C, which can be interpreted as a nonlinear, diffusive Oldroyd-B model, where we assume that the stiffness of the spring is not constant but depends on the elongation of the conformation tensor, i.e., on tr(C); see [28] for details and the references therein. Let us mention that, in analogy to the binary fluid model, the dissipative viscoelastic phase separation model can be derived in the context of the GENERIC formalism; see [26] for details. Hence our final model for viscoelastic phase separation readṡ Here φ denotes the volume fraction of the polymers while the other state variables have the same meaning as above. Furthermore, n(φ) denotes the mobility function and h 1 is the relaxation rate for q. Similarly, h 2 (φ)B(tr(C)) denotes a generalized relaxation rate, η(φ) is a volume fraction dependent viscosity, and A(φ) is related to the dynamic asymmetry, sometimes called bulk relaxation modulus; see [5,8]. Let us further mention that the elastic stress tensor T is now defined implicitly via an evolutionary equation for the conformation tensor C, which is assumed to be positive-definite. The term µ∇φ is related to the Korteweg interface stress and can be obtained in this form by a suitable redefinition of the pressure p.

Mathematical properties
In the following, we summarize the most important mathematical properties of our new model (8) for viscoelastic phase separation. For ease of notation, we again assume that the equations are complemented by appropriate periodic boundary conditions.

Thermodynamic consistency
Let us start by investigating the thermodynamic consistency of the proposed model. The free energy of the system (8) is given by which indicates the decomposition of the free energy via into several distinct contributions. By a lengthy but straightforward calculation [26], one finds the following energy-dissipation identity which establishes the consistency with the second law of thermodynamics.

Existence of solutions
In the following we review results of [26,27], where we rigorously proved the existence of global weak solutions for (8) in two space dimension. The proofs in these works are based on energy functionals similar to (9). We will distinguish the following two particular cases: The first, so-called regular case, assumes a smooth and strictly positive mobility function n(s) paired with a general polynomial-like double well potential, including the usual Ginzburg-Landau potentials. The second case is concerned with degenerate mobility functions, which are allowed to vanish at the points of single phases, paired with logarithmic potentials, like the physically relevant Flory-Huggins potential.
Lemma 3.1 (Regular case, [26]). Under suitable assumptions on the parameter functions, there exists at least one global weak solution of (8) for the Ginzburg-Landau type potential.
Lemma 3.2 (Degenerate case, [27]). Under suitable assumptions on the parameter functions, there exists at least one global weak solution of (8) for the Flory-Huggins type potential. Furthermore, the solution satisfies the physical bounds φ ∈ [0, 1].
Recall that weak solutions are, in particular, solutions in the sense of distributions. Moreover, the spatial and temporal regularity is sufficient such that the energy-dissipation structure (9) is preserved; we refer to [29] for details. Let us finally remark that similar existence result are not proven and actually not expected to hold for the model in [5].

Relative energy estimates, weak-strong uniqueness, and stability of solutions
In addition to pure existence of weak solutions, we now report on conditional uniqueness and stability results which encode the mathematical well-posedness of the model (8). The basic tool for the investigations of this section are relative energy estimates.
We start by defining a modified free energy functional for the problem under consideration, namely These modifications for the compositional and elastic degrees of freedom are introduced in order to obtain a functional that is strictly convex. This is the case if α > 0 is sufficiently large. Let us recall [30,31] that for any strictly convex functional E(z), we may define a relative energy This amounts to the quadratic Taylor remainder which, by strict convexity of E, is proportional to the quadratic distance between z andẑ. Note that the brackets denote a suitable inner product or dual pairing. From the physical point of view this can also be interpreted as statistical distance on the phase space with respect to the energy landscape or, more generally, as some measure for information.
Abbreviating the states by z = (φ, q, v, C) and z = (φ,q,v,Ĉ), the relative energy of our phase separation model reads to abbreviate the Taylor remainder. For α sufficiently large, the energy E(z) can be seen to be strictly convex, and the relative energy E(z|ẑ) satisfies E(z|ẑ) ≥ 0 and E(z|ẑ) = 0 ⇐⇒ z =ẑ.
Under additional assumptions, one can show that the relative energy allows to bound the squared distance, i.e., there exists a positive constant λ such that Using these properties, we can show the following decay estimated for the relative energy in a similar fashion as in [29]. Theorem 3.3. Let z = (φ, q, v, C) be a weak solution of (8) andẑ = (φ,q,v,Ĉ) be another sufficiently smooth solution of (8). Then for some constant C > 0 and with relative dissipation functional defined by As a direct consequence of Theorem 3.3, we obtain the following weak-strong uniqueness principle. This result can be obtained applying Theorem 3.3 since E(z|ẑ)(0) = 0 by construction, see (11). This implies E(z|ẑ)(t) = 0 for t ∈ (0, T ). Weak solutions, provided by Lemma 3.1, are thus unique, if at least one of them is sufficiently regular. Theorem 3.3 further allows to immediately deduce stability of solutions with respect to perturbations in the initial data. The following slight modification of the argument allows to study stability also with respect to further perturbations.
Before we state our results, let us note that (8) involved an additional variable µ = −c 0 ∆φ + ∂f (c) ∂c , i.e., the chemical potential. Now let (z, µ), with z = (φ, q, v, C) as before, denote a weak solution of (8), and let (ẑ,μ) be a corresponding given set of sufficiently smooth functions. By inserting (ẑ,μ) into (8), we can define residuals r i , i ∈ {φ, µ, q, v, C}, according to The brackets here denote the inner product in the space of square-integrable functions L 2 (Ω) or corresponding duality products, and the variational characterizations of the residuals are assumed to hold for all sufficiently smooth test functions (ψ, ξ, ζ, w, D).
We now obtain the following generalization of Theorem 3.3.
In contrast to Theorem 3.3, which allowed us to "only" study the effect of varying the initial conditions, Theorem 3.5 permits us to investigate the full spectrum of stability, in particular with respect to varying parameters, but also with respect to aspects like model variability and asymptotic limits. In the context of numerical analysis, a discrete version of the above relative energy estimate allows to derive convergence error estimates for appropriate numerical schemes. Closely related to the above results are relative entropy estimates in the context of coarse-graining, that are used to measure the information loss when traversing through a hierarchy of models. Appropriate versions of the above relative energy estimates can also be applied for such an analysis by choosing the functions (ẑ,μ) appropriately.

Comparison with a mesoscopic model
In this section we will compare results from our macroscopic model with those of a mesoscopic, particle based model.
In the mesoscopic simulation approach, implemented in the simulation package Espresso++ [32], the polymer component is modeled via a Kremer-Grest type [33] bead-spring model. The quality of the solvent can be varied by adjusting the attraction strength of the pair interaction based upon a modified Lennard-Jones potential [34]. This potential defines the scales for length σ LJ , energy ε LJ , and thereby time t LJ = (mσ 2 LJ )/ε LJ , where m is the particle mass. The polymers are dissipatively coupled to a Lattice-Boltzmann (LB) solvent background in order to include hydrodynamic interactions [35,36].
Note that since we consider φ as the volume fraction the relevant range is φ ∈ [0, 1] and we chose φ * as the mean value of the initial data for φ. Let us discuss the dimensionality. Since φ is dimensionless, Note that the parametric functions are chosen such that ηh −1 1 , n 2 h −1 1 , c 0 ∼ 1, meaning that all typical length scales match, and are of the order of one lattice spacing. Similarly we obtain matching timescales via this motivates us to define the time unit via t F E = c 0 η −1 . For typical simulation results we refer to [26]. We can clearly observe that the numerical simulations agree well with the real experiment presented in [5], see Figure 8 of Ref. [26]. The experiment represents the separation of a polystyrene-poly vinyl methyl ether (PS-PVME) mixture right after the quench. It is inspected via video phase-contrast microscopy. For a more detailed discussion of the PS-PVME experiments we refer to [37]. According to this observation Tanaka introduced six different regimes of the viscoelastic phase separation. We note that we can observe the same regimes also in our numerical simulations in [26], see Figure 1 of Ref. [26] and the associated discussion. Initial configurations were produced running the mesoscopic model at temperature T = 1 enforced by a Langevin thermostat with friction constant 1. For the initial equilibration, the polymer system was not coupled to the LB fluid. The system consists of 1024 polymer chains, each of which comprises 128 beads of mass m = 1, in a box of size 512 × 512 × 4. At equilibration time the nonbonded interaction is purely repulsive, corresponding to an attraction strength of λ = 0. The FENE bond potential has a strength of K = 30 and a maximum extension of 1.5. In order to generate configurations that are effectively two-dimensional, the particles are confined in the zdirection by an external potential while making sure that fluctuations in the z-direction are small compared to σ LJ . The LB fluid has shear and bulk viscosity of 3 and its density is 1; it is coupled to the particles with a Stokes friction parameter of 20. For every 10 Molecular Dynamics integrations, one LB step is performed. We extracted mean and variance of the mesoscopic density and used this for a starting configuration for the macroscopic simulations. In the mesoscopic model, phase separation is induced by suddenly introducing an attractive well with a depth of λ = 2 in the nonbonded potential.
Then the structure factor as defined by was calculated for both approaches at successive times t. More precisely, we average the structure factor over spherical shells, i.e.
In case of the mesoscopic model, configurations were interpolated onto a discrete two-dimensional lattice in order to evaluate the structure factor via Fast Fourier Transform. Four independent simulations were performed and the average of the structure factor curves taken. 25 independent runs were used in case of the macroscopic model. Note that at q = 0 the structure factor is time-independent, because of the conservation law Ω φ(x, t) dx = const..
In Fig. 1 we can see the structure factors of both methods. Larger systems allow for a higher resolution, as the smallest resolvable wavenumber depends on system size L via q 0 = 2π/L, resulting in a better resolution for the mesoscopic simulation data. One should note, however, that the mesoscopic simulations are orders of magnitude more costly in terms of CPU time compared to the macroscopic simulations.
The value of lim q→0 S(q) is very different from S(q = 0); in thermal equilibrium it is proportional to the compressibility. The figures show that the compressibilities of the two models differ substantially, because the equation of state was not yet adjusted to match. After the quench at t = 0 a peak starts to develop, which grows in height and moves towards smaller q values as time progresses.
In the second row of Fig. 1, the scattering intensity is normalized by its maximum value S(q max ) and the wavenumber scaled by q max . This fixes the peak at (1,1) in the spirit of the analysis of Tanaka and Araki [38]. The fact that the curves do not collapse onto a master curve indicates that the dynamic scaling hypothesis is violated and the demixing behavior is non-standard. In order to confirm that finding, we have conducted simulations of simple fluids and created an analogous plot. In case of the mesoscopic model, the same system was simulated without connectivity, which is accommodated for by a larger quench depth of λ = 3. For the macroscopic model elastic stresses where neglected, i.e. the equation and the coupling terms for q and C are removed. Note that to neglect q setting A(φ) = 0 is sufficient. In row three of Fig. 1 we see that the dynamic structure factor indeed shows a tendency to collapse for simple fluids. Disagreement for small wave numbers is expected due to finite-size effects.
Since the wavenumber with highest scattering intensity q max is inversely proportional to a characteristic length scale in the system, q max (t) reflects the coarsening dynamics in phase separation. In Fig. 2 we show q max (t) as obtained from mesoscopic and macroscopic simulations in a double logarithmic plot. In the initial phase the agreement is quite good and both approaches show a growth-law with an exponent close to the Lifshitz-Slyozov exponent of −1/3, hinting at diffusive dynamics. At about t ≈ 512∆t however, the macroscopic data starts to deviate and becomes flatter, indicating a transition to a different dynamical regime. This is confirmed by Fig. 3 which shows a sharp rise in scattering intensity at this point, whereas the behavior of the mesoscopic simulation is more regular.
This deviation is not surprising, since the parameters of the macroscopic model are not really calibrated to the mesoscopic description. The peak wave number is highly determined by the interplay of ). The first row shows the normalized scattering intensity S(q, t)/S(q = 0) vs. the wavenumber q. Note that the normalization implies a value of unity at q = 0. Sinced it would be completely off scale, this value is omitted. In the second row, the scattering intensity has been normalized by its peak value S(qmax) and the wavenumber by qmax. For comparison we show another normalized plot for simple fluids without any elastic effects in the third row.
interface width c 0 and the choice of potential f (φ). However c 0 is directly related to the length scale. In the future we want to calibrate the interface width and the potential in order to obtain more accurate results. Furthermore, we want to investigate the effects of the length scales by conducting a variety of simulations based on different coarse-graining levels of the mesoscopic data, i.e. using more data than only mean and variance.

Discussion and outlook
In this work we have first considered a suitable reduction of a modified binary fluid model. The system is closed by deriving a relation between the relative velocity and suitable state variables that arise.
In the course of this model reduction we obtain a relation for the relative velocity that is quite similar to the one given in [8]. To arrive at this model, we make essentially the same (probably partly debatable)   Normalized peak-scattering intensity S(qmax, t)/S(0, t) for mesoscopic (∆t = t LJ ) and macroscopic simulations (∆t = t FE ) respectively. The difference in scale is due to discrepancies in the equation of state of both approaches.
phenomenological assumptions as [8], but with a derivation that makes the elimination of the relative velocity more explicit. Furthermore, the construction gives insight what the concrete physical interpretation of q may be, since the interpretation of T v or C is fairly clear. Inclusion of dissipative terms and suitable evolution equations for the viscoelastic effects lead to our new macroscopic model (8).
In the second part we shortly discuss the thermodynamic consistency of the proposed model and furthermore present our results on existence of weak solution in two relevant cases. Afterwards we introduced the notion of relative energy, as a suitable distance in terms of energy. Finally, we discussed the stability and weak-strong uniqueness by means of relative energy methods. All together the section implies well-posedness of the problem.
In the last part we compared the results of the microscopic description with our macroscopic model. We can see that the structure factor is a valuable tool in comparing key features of both approaches.
In order to perform a more thorough comparison however, the parameters of both models need to be adjusted to match, and to do so the interpretability of the macroscopic equations in terms of microscopic physics needs to be improved.
This was the main motivation for parallel investigations by us, which aimed at the development of a viscoelastic phase separation model from scratch, i.e. starting from a simple microscopic model and then applying coarse-graining. This alternative model, which has the advantage of being well-rooted in microscopic physics, is presented in our companion paper [24], with analytical results that look quite encouraging. While the equations share some similarities with those presented here, we note that the bulk stress is not present and the energetic structure is different. Whether the closure of the present paper can be justified from microscopic physics remains an open question, into which we will look further in the future. We also plan to study the differences between the models in order to deepen our understanding of the problem.