Existence of Weak Solutions for a Diffuse Interface Model for Two-Phase Flow with Surfactants

Two-phase flow of two Newtonian incompressible viscous fluids with a soluble surfactant and different densities of the fluids can be modeled within the diffuse interface approach. We consider a Navier-Stokes/Cahn-Hilliard type system coupled to non-linear diffusion equations that describe the diffusion of the surfactant in the bulk phases as well as along the diffuse interface. Moreover, the surfactant concentration influences the free energy and therefore the surface tension of the diffuse interface. For this system existence of weak solutions globally in time for general initial data is proved. To this end a two-step approximation is used that consists of a regularization of the time continuous system in the first and a time-discretization in the second step.


Introduction and Main Result
Two-phase flow of two macroscopically immiscible, viscous, incompressible fluids with a soluble surfactant is very relevant in many applications and can be modelled by diffuse interface models.Surface active agents (surfactants) are chemical substances which lower the surface tension at fluid interfaces.They play an important role as detergents, emulsifiers, wetting agents, dispersants and foaming agents and hence mathematical modeling, analysis and numerical computations involving surfactants have many potential applications.In order to describe surfactants one has to model the fluid flow in the bulk regions away from the interface which is typically done with the help of the Navier-Stokes equations.In addition, diffusion and advection of the surfactants has to be taken into account both at the interface and in the bulk.At the interface transport phenomena and force balances also involving surface tension effects have to be modeled.As the presence of surfactants leads to a non-constant surface tension also Marangoni forces at the interface play an important role.
Classically the effects of surfactants at fluid interfaces are modeled within the context of sharp interface models, i.e., the interface is modeled as a hypersurface.We refer to [13,14,20] for analytical results and to [9,11,12,18,21,22,23,27,34] for numerical results on sharp interface approaches for surfactant flow.
However, in the context of sharp interface approaches a change of topology, i.e., for example a splitting of a droplet or a reconnection of droplets, is difficult to describe theoretically and typically happens in an ad-hoc way numerically.Therefore, in the last ten years diffuse interface approaches (also called phase field approaches) have been introduced in the context of surfactants in fluid flow.We refer to [17,19,24,26,30,31,32] for such approaches.
In this paper we will rely on the work of Garcke, Lam and Stinner [19] where for the first time diffuse interface models were introduced which allowed for a physically sound energy inequality.This fact is crucial for the analysis in this paper as it allows to derive important a priori estimates.The models in [19] generalize the thermodynamically consistent diffuse interface model for two-phase flow with different densities introduced in [6] to situations with surfactants.Existence of solutions for the model in [6] has been shown in [3,4].However, the case with surfactants turns out to be much more involved due to the complex, nonlinear coupling of the Navier-Stokes/Cahn-Hilliard system to the equations describing the surfactants.This paper gives a first existence result for this general situation.
In the following we consider a diffuse interface model for such a two-phase flow, where a partial mixing of the macroscopically immiscible fluids is taken into account and an order parameter, here in form of the difference of volume fractions ϕ, is introduced.The interface is no longer treated explicitly in form of a lower dimensional surface.It is replaced by an interfacial region, where the order parameter ϕ is not close to +1 or −1.In a sufficiently smooth situation the interfacial region has a thickness proportional to ε > 0, where ε > 0 is a parameter in the system.
The model that we will analyze is a variant of the "model C" in [19].In the latter contribution different models for the case of soluble and insoluble surfactants were derived.Here we consider a model, where the surfactant is soluble in both fluids and accumulates at the interface.The model leads to the system ∂t(ρv) + div(v ⊗ (ρv + J)) where Q = Ω × (0, ∞), ∂ • t := ∂t + v • ∇ is the material time derivative and J = ∂ρ(ϕ) ∂ϕ (− m(ϕ)∇µ). (1.6) Throughout the paper Ω ⊆ R d , d = 2, 3, is a bounded domain with C 2 -boundary.We close the system by prescribing initial values and boundary conditions In this model, v is the mean velocity of both fluids and Dv := 1 2 (∇v + ∇v T ).Moreover, ϕ denotes the order parameter, which is taken as the volume fraction difference of both fluids, ε > 0 is a constant parameter related to the "thickness" of the interfacial region and p is the pressure.The viscosity of the mixture is denoted by η(ϕ) > 0 and W : R → R is a potential, which is part of the free energy of the fluid mixture and of double well shape with two global minima at ±1.The chemical potential for the fluids is denoted by µ, while q is the chemical potential for the surfactants.Furthermore, m(ϕ) > 0 and m(ϕ, q) > 0 are the mobilities for the phase field and the surfactant diffusion equations, respectively.The density ρ depends on ϕ and is given by In (1.1), J is a flux of the fluid density relative to ρv.The term h(q) is the square of the surface tension.Furthermore, thermodynamic relations d ′ (q) = f ′ (q)q, d(q) = h(q) + f (q)q, G ′ (q) = g ′ (q)q.
(1.9) are required.The first two equations imply h ′ = −f and hence d = h − h ′ q.This means that d is the Legendre transform of h and this fact is related to the classical Gibbs relation from thermodynamics, see [7].In [7] it is discussed how the diffuse interface model (1.1)-(1.5)can be related to classical descriptions, see [13,14,20], involving a sharp interface.The equations (1.1) and (1.2) describe the conservation of mass and linear momentum of the fluid mixture, (1.3) describes the diffusion of the surfactants in the bulk of both fluids as well as along the diffuse interface and (1.4)-(1.5)describes the dynamics of the diffuse interface, the order parameter ϕ, respectively.Moreover, is the free energy density of the fluid mixture and 1 ε f (q)qW (ϕ) + G(q) is the free energy density of the surfactants at the diffuse interface (related to f (q)q) and in the bulk (related to G(q)).Hence the total free energy density is given by The total energy of the system is given by together with the boundary conditions (1.8) satisfy the energy identity for all t ∈ (0, T ).This is proved by testing (1.1), (1.3)-(1.5)with v, q, µ, and ∂tϕ, respectively, integration by parts and using the relations (1.9).We refer to [33, Section 2.2] for the details.
To this end it is essential that ∂ρ(ϕ) ∂ϕ for all ϕ ∈ [−1, 1] is constant and therefore the continuity equation is the diffuse interface model for a two-phase flow of incompressible, viscous fluids with different densities (without surfactants) that was derived in [6].Existence of weak solutions to the model was shown in [3] for non-degenerate mobility m(ϕ) and singular free energy density W and for degenerate mobility in [4].Existence and uniqueness of strong solutions locally in time was shown in [33].We refer to [5] for an overview of the different models for two-phase flows without surfactants and analytic results.
In a situation where the order parameter obtains values outside of the interval [−1, 1], we obtain a modified continuity equation given by where we used (1.4).In order to preserve the energy identity above we modify (1.1) to Here 1 2 Rv describes the change in the kinetic energy associated to the source term R in the modified continuity equation.We note that this modified continuity equation and the additional source term R is also used in [2].We also remark that the additional source term is zero for physically meaningful values ϕ ∈ [−1 , 1].
With (1.5), equation (1.13) can equivalently be written as where a suitable scalar function is added to p, see also [19].Before we state our main result we define a weak solution for (1.2)-(1.8),(1.13): a weak solution of (1.2)-(1.8),(1.13) if the following equations are satisfied: , where J is defined by (1.6) and Moreover, the energy inequality (1.17) has to hold for all t ∈ [s, ∞) and almost all s ∈ [0, ∞) including s = 0, where Etot is defined as in (1.10).
We refer to Section 2 below for the definition of the function spaces above.Remark: We note that (1.15) contains a weak formulation for the initial value of f (q)W (ϕ) + g(q).More precisely using (1.15) and Assumption 1.2, one can show and f (q)W (ϕ) + g(q)|t=0 = f (q0)W (ϕ0)g(q0) in H −1 (Ω).By this the initial value of q is prescribed implicitly since q → f (q)W (ϕ) + g(q) is strictly monotone, cf.Assumption 1.2 below, and the initial value for ϕ is determined by the evolution equation for ϕ.Finally, we note that (1.14) prescribes the initial value of ρv as ρ(ϕ0)v0.Divergence free test functions are sufficient for this because of the same observations as in [3,Section 5.2].
To obtain the existence of weak solutions in the sense of Definition 1.1, we make the following assumptions.
From these assumptions it follows that g is strongly monotone due to Moreover, these assumptions imply f ′ (q) = 0 for all q / ∈ [qmin, qmax] as it holds f ′ (q)q = d ′ (q).In particular, we observe that f is a bounded function.Due to h(q) = d(q) − f (q)q, we can deduce that there exists a constant C > 0 such that |h(q)| ≤ C(|q| + 1).
(1. 19) In applications physically relevant values for q lie in a bounded interval and hence the assumptions for the values of functions for |q| large give no restrictions in practice.This is due to the fact that we can modify the functions for large |q| as we like.
Our main result is:

Remark: Due to the term ∂
3) one of the main difficulties in proving Theorem 1.3 is to show compactness with respect to q in suitable approximating problems.The structural assumptions on d, f, g, h, G will first allow us to show an energy inequality related to (1.11) and then enables us to show compactness of q.Moreover, the treatment of the term R is subtle.To the end we approximate the system in two steps.First we regularize the system in order to obtain ∂tϕ ∈ L 2 ((0, ∞) × Ω)).Then we solve the regularized system with the aid of an implicit time discretization.

Preliminaries
Notation: Let X be a Banach space and let X ′ be its dual space.Then the duality product is given by x ′ , x X ′ ,X = x ′ , x = x ′ (x), where x ′ ∈ X ′ and x ∈ X.The natural numbers without 0 are denoted by N and we set N0 Function spaces: In the following, let Ω be a bounded domain with C 2 -boundary.For 1 ≤ p ≤ ∞ we denote by L p (Ω) the usual Lebesgue-space equipped with its norm If Ω = (0, T ) the Banach-space valued variants are denoted by L p (0, T ; X) and W k p (0, T ; X), respectively.Furthermore, and C ∞ (0) ([0, ∞), X) denotes the space of smooth function f : [0, ∞) → X with a support that is compactly contained in [0, ∞).
We set The operator P0 is the orthogonal projection onto L 2 (0) (Ω) given by P0f For the proof of existence of weak solutions, we need the following compactness result.Theorem 2.1.Let X ⊆ B ⊆ Y be Banach spaces with compact embedding X ֒→ B. Furthermore, let 1 ≤ p ≤ ∞ and assume that The proof of this theorem can be found in [29,Theorem 5].

Existence of Weak Solutions for the Surfactant Model
We approximate the system (1.3)-(1.5),(1.13) by the following equations: ) ) and the energy inequality has to hold for all t ∈ [s, ∞) and almost all s ∈ [0, ∞) including s = 0.
Remark 3.2.Note the difference between (1.16) and (3.8).If all appearing terms are smooth enough, both definitions are equivalent due to the modified continuity equation (1.12).We use (3.8) because this representation allows us to derive an energy inequality.But as the term ∂tρ(ϕ) appears in (3.8), we also need to estimate this term.Therefore, we insert the additional term δ∂tϕ in equation (3.11) and approximate the initial system.Hence, we first need to solve the approximating system and then pass to the limit δ → 0. But due to the additional terms δ∂tϕ and δ∆v we are then able to use (1.16) instead of (3.8) in the proof of Theorem 1.3.
In order to prove our main result we use: and q0 ∈ L 2 (Ω) be given.Then there exists a weak solution ).This theorem is proven in Section 4 with the aid of an approximation by a time-discrete problem for a time-step size h > 0 and passing to the limit h → 0. The arguments are a non-trivial generalization to the case with surfactants of the proof of the corresponding result in [3].
Theorem 3.3 yields the existence of weak solutions (v δ , ϕ δ , µ δ , q δ ) of the approximating system (3.1)-(3.5) in the sense of Definition 3.1.It remains to prove the existence of weak solutions for the system (1.1)-(1.5).To this end, we have to pass to the limit δ → 0. We can assume w.l.o.g.Ω ϕdx = 0.By changing ϕ by a constant and shifting W , we can always reduce to this case.Moreover, we will assume ε = 1 for simplicity in the following.Furthermore, we split the equation ) by considering the system ) ) In order to get uniform bounds for the solutions of (3.16)-(3.20)for ϕ δ 1 and ϕ δ 2 , we use the following lemma.
We show these bounds in detail.
Ad iv) The growth conditions for h and W ′ yield the estimate By the Gagliardo-Nirenberg inequality and the Hölder inequality we have the embedding Together with the boundedness of q δ ∈ L 2 (0, T ; L 6 (Ω)) and ∇ϕ δ ∈ L ∞ (0, T ; L 2 (Ω) d ) we get the statement.

Existence for the Approximate System
It is the aim of this section to prove the existence of weak solutions for the approximating system (3.1)-(3.5).To this end, we determine an appropriate time discretization.We solve the time-discrete problem by using the Leray-Schauder principle, cf.Theorem 4.2 below.Then we prove Theorem 3.3.
For the time discretization, we set h (Ω) and q k ∈ L 2 (Ω) be given.We determine (v k+1 , ϕ k+1 , µ k+1 , q k+1 ) as a weak solution of the system ) with boundary conditions where and H : R × R → R is defined by We call a weak solution of (4.1)-(4.6)for initial data for all ψ ∈ V (Ω), where Jk+1 is defined as in (4.7), and ) for all φ ∈ H 1 (Ω), where we define for ψ ∈ V (Ω) Note that (4.9) can equivalently be written as Theorem 4.2 (Existence of weak solutions for the time-discrete problem).
The proof generalizes the proof of [3,Lemma 4.2], where the existence of weak solutions for the model without surfactants developed in [6] is shown.

Proof of Theorem 4.2.
We start with the proof that the energy estimate (4.14) holds for any weak solution of (4.1)-(4.5) in the sense of Definition 4.1.To this end, we test equation (4.9) with v k+1 (i.e., choose ψ = v k+1 ), (4.10) with q k+1 , (4.11) with µ k+1 and (4.12) with . Then one important tool to simplify the equations is the identity Note that when we test (4.12) with , we use the identity In (4.10) tested with q k+1 , we use f (q) = −h ′ (q) for all q ∈ R and the fact that h is a concave function to obtain Moreover, we can estimate and have d(q) = h(q) + f (q)q for all q ∈ R. Altogether this yields (4.14).Now we need to prove the existence of weak solutions for the time-discrete problem (4.1)-(4.5).To this end, we define two operators L k , F k : X → Y and apply the Leray-Schauder principle, where Moreover, we define for , where Jk+1 is given as in (4.7).Note that for the definition of F k : X → Y we used (4.13), which is equivalent to (4.9).Then it holds if and only if w k+1 = (v k+1 , q k+1 , µ k+1 , ϕ k+1 ) ∈ X is a weak solution of (4.1)-(4.5).For the equation with F ∈ Y , the Lax-Milgram theorem yields the existence of unique weak solutions v k+1 ∈ V (Ω) and ϕ k+1 , µ k+1 , q k+1 ∈ H 1 (Ω).By elliptic regularity we obtain that µ k+1 ∈ H 2 n (Ω).Hence L k : X → Y is invertible and the open mapping theorem yields the boundedness of To this end, we introduce the Banach space These estimates follow from Sobolev embeddings, the boundedness of f together with the growth conditions |h(q)| ≤ C(|q| + 1), |W (q)| ≤ C(|q| 3 + 1), |W ′ (q)| ≤ C(|q| 2 + 1) and |G ′ (q)| ≤ C(|q| + 1) for all q ∈ R and a constant C > 0.Moreover, we use the identities f = −h ′ , G ′ (q) = g ′ (q)q and the fact that f is constant outside an interval [qmin, qmax].
The continuity of the linear terms of F k : X → Ỹ follows from their boundedness.For the nonlinear terms, the continuity can be proven with the aid of the theorem on continuity of Nemyckii operators on L p -spaces and with the multilinear structures.Since F k : X → Ỹ is a continuous and bounded operator and since the embeddings (Ω) ֒→ L 2 (Ω) are compact, we can conclude that F k : X → Y is a compact operator.
In the following, we want to apply the Leray-Schauder principle.We already noted that w k+1 ∈ X is a weak solution of (4.1) -(4.5) if and only if We set : Y → Y and note that proving the existence of a weak solution for (4.1) -(4.5) is equivalent to proving the existence of a fixed-point of We can prove the existence of such a fixed-point with the Leray-Schauder principle, cf.[35, Theorem 1.D.].We have to show: To this end, let g k+1 ∈ Y and 0 ≤ λ < 1 be given such that g k+1 = λK k g k+1 .As Testing (4.9) with v k+1 , (4.10) with q k+1 , (4.11) with µ k+1 and (4.12) with , using the same identities as in the derivation of the energy inequality and omitting some non-negative terms, we can estimate For more details, we refer to [33,Section 3.2.3].Thus (4.15) is fulfilled and the Leray-Schauder principle yields the existence of Finally, we need to show higher regularity for ϕ k+1 .From L k (w k+1 ) = F k (w k+1 ) with w k+1 = (v k+1 , q k+1 , µ k+1 , ϕ k+1 ) it follows where the right-hand side is bounded in the L 2 -norm.Thus elliptic regularity theory yields ϕ k+1 ∈ H 2 n (Ω).Hence, there exists a weak solution for the time-discrete problem (4.1)-(4.5) in the sense of Definition 4.1, which fulfills the discrete energy estimate (4.14).
Using Theorem 4.2, we can prove the existence of weak solutions for the approximating system (3.1)-(3.5).
The proof of the strong convergence of (q N ) N∈N in L 2 (0, T ; L 2 (Ω)) for a suitable subsequence is a bit more complicated than the proof for (q δ ) δ>0 in Theorem 1.3.As before, we show the compactness of (q N ) N∈N in L 2 (QT ) with the aid of Theorem 2.1.First, we note that (q N ) N∈N is bounded in L 2 (0, T ; H 1 (Ω)) and therefore condition i) is fulfilled for every 0 < T < ∞.For the proof of condition ii) we show ≤ C(T )s for any λ > 0 and a constant C(T ) > 0 independent of λ and N .This shows condition (ii) of Theorem 2.1 and we obtain the compactness for the sequence (q N ) N∈N , Therefore let s = mh be given for m ∈ N and h = 1 N .Moreover, we define F (ϕ N , q N ) := 1 ε f (q N )W (ϕ N ) + g(q N ), f (t) := F (ϕ N (t), q N (t)).
From the strong monotonicity of g and the monotonicity of f it follows that F (ϕ N , •) is strongly monotone as before.Thus there exists a constant C > 0 such that for every t ∈ [kh, (k + 1)h).Multiplying these inequalities with |q N (t + s) − q N (t)|, integrating from (k − 1)h to kh with respect to t, integrating over the domain Ω and summing over Since f is a bounded function, we can conclude for the first term in (4.22) , a.e. in (0, ∞) for every 0 < T < ∞ and a constant C(T ) > 0 depending on T .The proof for the latter inequality is similar to the proof of (3.24), where we use that there exists k ∈ N such that t ∈ [kh, (k + 1)h) and t ) for all 0 < T < ∞, the statement follows similarly as in (3.24).
For the second term in (4.22), we set l(t) := t h and t(t) := h t h .Then it holds t(t) = t k for t ∈ [kh, (k + 1)h).Hence, we have We use this identity in (4.22) and obtain T −s 0 Ω F (ϕ N s+ , q N s+ ) − F (ϕ N , q N ) q N s+ − q N dx dt In the following, we study these three terms separately.The first term can be estimated by For the second term we get in the same way g(q N (τ + h)v N (τ + h) d τ ∇q N s+ − ∇q N dx dt ≤ C(T )s 1 4 .

Definition 4 . 1 (
Note that H(a, b)(a − b) = W (a) − W (b) for every a, b ∈ R. It remains to define a weak solution for the time-discrete problem (4.1)-(4.6).Weak solution of the time-discrete problem).

1 4 ( 4
.20) for all s = mh with m ∈ N and a constant C(T ) > 0 independent of s and N ∈ N. Then[10,  Lemma 9