Gradient flow dynamics of two-phase biomembranes: Sharp interface variational formulation and finite element approximation
GGradient flow dynamics of two-phase biomembranes: Sharpinterface variational formulation and finite element approximation
John W. Barrett Harald Garcke Robert N¨urnberg Department of Mathematics, Imperial College London, London, SW7 2AZ, UK
E-mail address : [email protected] Fakult¨at f¨ur Mathematik, Universit¨at Regensburg, 93040 Regensburg, Germany
E-mail address : [email protected] Department of Mathematics, Imperial College London, London, SW7 2AZ, UK
E-mail address : [email protected] . Abstract.
A finite element method for the evolution of a two-phase membrane in a sharp interface for-mulation is introduced. The evolution equations are given as an L –gradient flow of an energy involvingan elastic bending energy and a line energy. In the two phases Helfrich-type evolution equations areprescribed, and on the interface, an evolving curve on an evolving surface, highly nonlinear boundaryconditions have to hold. Here we consider both C – and C –matching conditions for the surface at theinterface. A new weak formulation is introduced, allowing for a stable semidiscrete parametric finiteelement approximation of the governing equations. In addition, we show existence and uniqueness fora fully discrete version of the scheme. Numerical simulations demonstrate that the approach can dealwith a multitude of geometries. In particular, the paper shows the first computations based on a sharpinterface description, which are not restricted to the axisymmetric case. Keywords. parametric finite elements, Helfrich energy, spontaneous curvature, multi-phase mem-brane, line energy, C – and C –matching conditions. Math. classification. Introduction
Two-phase elastic membranes, consisting of coexisting fluid domains, have received a lot of attentionin the last 20 years. The interest in two-phase membranes in particular was triggered by the multitudeof different shapes observed in experiments with inhomogeneous biomembranes and vesicles. Biomem-branes are typically formed as a lipid bilayer, and often multiple lipid components are involved, whichlaterally can separate into coexisting phases with different properties. Among the complex morpholo-gies that appear are micro-domains, which resemble lipid rafts, and these are of huge interest in biologyand medicine. As the thickness of the membrane is much smaller than its lateral length scale, typicallythe membrane is modelled as a two-dimensional hypersurface in three dimensional Euclidean space.The equilibrium shape of the membrane is obtained by minimizing an energy which –besides othercontributions– contains bending energies involving the mean curvature and the Gaussian curvatureof the membrane. If different phases occur, parameters in the curvature energy are inhomogeneous,leading to an interesting free boundary problem as well as to a plethora of different shapes. We refer to[12], where multi-component giant unilamellar vesicles (GUVs) separating into different phases werestudied. These authors were able to optically resolve interactions between the different phases, itscurvature elasticity and the line tension of its interface.There have been several studies on theoretical and numerical aspects of two-phase membranes takingcurvature elasticity and line energy into account, see e.g. [28, 29, 39, 11, 40, 15, 30, 22, 23, 24, 25, 26,13, 32, 14, 9], which we discuss in the following. 1 a r X i v : . [ m a t h . NA ] A p r . W. Barrett, H. Garcke, R. N¨urnberg The by now classical model for a one-phase membrane rests on the Canham–Helfrich–Evans elasticbending energy α (cid:90) Γ ( κ − κ ) d H + α G (cid:90) Γ K d H , where Γ is a closed two-dimensional hypersurface and H denotes the two-dimensional Hausdorffmeasure. The mean curvature of Γ is denoted by κ , and K is its Gaussian curvature. The constants α and α G are bending rigidities, while κ is the spontaneous curvature reflecting asymmetry in themembrane introduced, for instance, by different environments on both sides of the membrane.In a fundamental work, J¨ulicher and Lipowsky ([28, 29]) generalized the Canham–Helfrich–Evansmodel to two-phase membranes. The geometry is now given by two smooth surfaces Γ and Γ , witha common boundary γ . In general, the constants α , α G and κ take different values in the two phasesΓ and Γ , which we will denote with an index i . On the curve γ line tension effects play an importantrole, and the total energy introduced in [28, 29] is given as E ((Γ i ) i =1 ) = (cid:88) i =1 (cid:20) α i (cid:90) Γ i ( κ i − κ i ) d H + α Gi (cid:90) Γ i K i d H (cid:21) + ς H ( γ ) , (1.1)where the constant ς ∈ R ≥ denotes a possible line tension, and where an index i ∈ { , } states thatquantities such as the curvatures and physical constants are evaluated with respect to Γ i . Of course, H denotes the one-dimensional Hausdorff measure.In [29] it is assumed that the surface Γ = Γ ∪ γ ∪ Γ is a C –surface, meaning in particular thatthe normal to Γ is continuous across the phase boundary γ . The works [25, 26, 27], on the other hand,also allow for discontinuities of the normal at γ . The first variation of the energy E in (1.1) has beenderived in [22] for the C –case and in [41] for the C – and the C –case. It is the goal of this paperto develop a numerical method for a gradient flow evolution of the energy E . To be more precise, wewill consider an evolution of the form (cid:68) (cid:126) V , (cid:126)χ (cid:69) Γ + (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ = (cid:20) δδ Γ E ((Γ i ) i =1 ) (cid:21) ( (cid:126)χ ) . (1.2)Here (cid:126) V is the velocity of the surface, δδ Γ E is the first variation of the energy, (cid:126)χ is a test vector fieldon the surface related to directions in which one perturbs the given surface Γ, and (cid:37) ≥ (cid:104)· , ·(cid:105) Γ and (cid:104)· , ·(cid:105) γ denote the L –inner products on the surface Γ and on the curve γ , respectively. The evolution of the surface is hence given as a steepest descent dynamics with respectto a weighted L –inner product that combines contributions from the surface and the curve. It willturn out that the governing equations in the case where the surface is restricted to be C are (cid:126) V = [ − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | ] (cid:126)ν i on Γ i , (1.3)together with the boundary conditions on γ = ∂ Γ i : α ( κ − κ ) + α G (cid:126) κ γ . (cid:126)ν = α ( κ − κ ) + α G (cid:126) κ γ . (cid:126)ν , (1.4a)[ α i ( ∇ s κ i )] . (cid:126)µ − [ α Gi ] τ s + ς (cid:126) κ γ . (cid:126)ν = (cid:37) (cid:126) V . (cid:126)ν , (1.4b) − [ α i ( κ i − κ i ) ] + [ α i ( κ i − κ i ) ( κ i − (cid:126) κ γ . (cid:126)ν )] + [ α Gi ] τ + ς (cid:126) κ γ . (cid:126)µ = (cid:37) (cid:126) V . (cid:126)µ . (1.4c)Equation (1.3), with ∆ s and ∇ s denoting the surface Laplacian and the surface gradient on Γ i , respec-tively, is Willmore flow taking spontaneous curvature effects into account. The boundary condition(1.4a), with (cid:126) κ γ denoting the curvature vector on γ ( t ), generalizes the equation for the mean curvaturein Navier boundary conditions, appearing for example in [18, (6)]. The equations (1.4b,c), with τ being the geodesic torsion of the curve γ ( t ) on Γ( t ) and with [ a i ] = a − a denoting the jump of a across γ ( t ), appear in the case (cid:37) = 0 in [23, (3.17), (3.18)], where additional terms to fix the surfaceareas and the enclosed volume appear. In the axisymmetric case, the equations (1.4a–c) reduce to2 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES the equations studied in [29]. Similar conditions have been derived in [39], and it has already beendiscussed in [23, Appendix B] that these authors miss one term. For positive (cid:37) the equations (1.4b,c)give rise to dynamic boundary conditions taking into account an additional dissipation mechanism atthe boundary. A similar condition for semi-free boundary conditions has been analyzed in [1, (1.3)].For evolutions where the surface areas of Γ and Γ , as well as the volume enclosed by Γ, are conserved,additional terms appear in (1.3) and (1.4c), see (2.16) and (2.20c), below. Moreover, in the case thatthe surface Γ is just continuous, the boundary conditions (1.4a–c) have to be replaced, and we referto (2.19a,b), below, for the relevant equations.Numerically mainly the C –case has been studied, with the exception of [26], where C –surfaceswith kinks in the axisymmetric case were studied numerically with the help of a phase field method. Inthe C –case already in [29] several two-phase equilibrium shapes in the axisymmetric case were com-puted by solving a governing boundary value problem for a system of ordinary differential equations.Based on research on model membranes, see [12], it has now become possible to perform a systematicanalysis of the influence of parameters also in the case of two-phase coexistence. We refer to [11],where experimental vesicle shapes were compared with shapes obtained by solving numerically theaxisymmetric shape equations derived in [29]. In this context, we also refer to [14], where, in contrastto the above works, also the effect of spontaneous curvature is taken into account in the axisymmetriccase. These authors were able to show that spontaneous curvatures already in an axisymmetric setupgive rise to a multitude of morphologies not seen in the case without spontaneous curvature.Almost all numerical results mentioned so far were for a sharp interface setup. Another successfulapproach uses a phase field to describe the two phases on the membrane. Line energy in this contextis replaced by a Ginzburg–Landau energy like in the classical Cahn–Hilliard theory. We refer to[40, 30, 22, 23, 24, 25, 26, 32, 31] for numerical results based on the phase field approach. The abovepapers use a gradient flow approach to obtain equilibrium shapes in the large time limit. An evolutionlaw using a Cahn–Hilliard equation on the membrane coupled to surface and bulk (Navier–)Stokesequations has been studied by the present authors in [9].Rigorous analytical results for two-phase elastic membranes are very limited. So far only results forthe axisymmetric case are known. We refer to the work [13], where the existence of global minimizersfor axisymmetric multi-phase membranes was shown, and the works [25, 26, 27], where the sharpinterface limit of the phase field approach in an axisymmetric situation was studied. Existence resultsfor the evolution problem are not available in the literature so far and should be addressed in thefuture.It is the goal of this paper to introduce a finite element approximation for a gradient flow dynamicsof the membrane energy E , which is based on a sharp interface approach. Instead of using a phasefield on the membrane, we will directly discretize the curve γ separating the two phases Γ and Γ . Inthree dimensions the total surface Γ will be discretized with the help of polyhedral surfaces consistingof a union of triangles. The curve γ is discretized as a polygonal curve in R fitted to the discretizationof Γ in the sense that the polygonal curve is the boundary of the open polyhedral sets Γ and Γ . Theboundary conditions (1.4a–c) are highly nonlinear and involve derivatives of an order up to three whenformulated with the help of a parameterization. It is hence highly non-trivial to discretize them in apiecewise linear setup. In this work, a splitting method is used, which basically uses the position vectorsof the nodes and an approximation of the mean curvature vector as unknowns. The approach in thispaper relies on a discretization of mean curvature leading to good mesh properties. This discretizationwas introduced by the present authors in [2, 3] and has been previously used for closed and openmembranes, see [8, 10] and for elastic curvature flow of curves with junctions, see [6].We will use the variational structure of the problem to derive a discretization which will turn outto be stable in a spatially discrete and continuous-in-time semidiscrete formulation. In order to do so,we will make use of an appropriate Lagrangian and will use ideas of PDE constrained optimization.3 . W. Barrett, H. Garcke, R. N¨urnberg (cid:126) id s (cid:126)µ (cid:126)µ (cid:126)ν (cid:126)ν Γ Γ γ Figure 1.
Sketch of Γ = Γ ∪ γ ∪ Γ with outer unit normals (cid:126)ν i , conormals (cid:126)µ i andtangent vector (cid:126) id s on γ for the case d = 3.The outline of this paper is as follows. In the subsequent section we will formulate the governingequations with all the details. In Section 3 a weak formulation is introduced using the calculus of PDEconstrained optimization. A semidiscrete discretization is formulated in Section 4. For this scheme alsoenergy decay properties and conservation properties are shown. In Section 5 a fully discrete version ofthe scheme is introduced, leading to a linear system at each time level, which is shown to be uniquelysolvable. In Section 6 we discuss ideas on how to solve the resulting linear algebra problems numerically.In Section 7 we present several numerical results showing that the new approach allows to approximatesolutions to the governing equations also in highly nontrivial geometries. In an appendix we show thatthe weak formulation derived in this work yields in fact the strong formulation for sufficiently smoothevolutions. The governing equations
In this section we precisely formulate the governing equations both for the C – and the C –case. Wealways assume that (Γ( t )) t ∈ [0 ,T ] is an evolving hypersurface without boundary in R d , d = 2 ,
3, that isparameterized by (cid:126)x ( · , t ) : Υ → R d , where Υ ⊂ R d is a given reference manifold, i.e. Γ( t ) = (cid:126)x (Υ , t ).Then (cid:126) V ( (cid:126)q, t ) := (cid:126)x t ( (cid:126)z, t ) ∀ (cid:126)q = (cid:126)x ( (cid:126)z, t ) ∈ Γ( t ) (2.1)defines the velocity of Γ( t ). In order to introduce the two-phase aspect, we consider the decompositionΓ( t ) = Γ ( t ) ∪ γ ( t ) ∪ Γ ( t ), where the interiors of Γ ( t ) and Γ ( t ) are disjoint and γ ( t ) = ∂ Γ ( t ) = ∂ Γ ( t ).We assume that each Γ i ( t ) is smooth, with outer unit normal (cid:126)ν i ( t ). See Figure 1 for a sketch of thesetup in the case d = 3. In particular, we parameterize the two parts of the surface over fixed oriented,compact, smooth reference manifolds Υ i ⊂ Υ, i.e. we let Γ i ( t ) = (cid:126)x (Υ i , t ), i = 1 ,
2. Throughout thispaper we will investigate two different types of junction conditions on γ ( t ): C –case : γ ( t ) = ∂ Γ ( t ) = ∂ Γ ( t ) , (2.2a) C –case : γ ( t ) = ∂ Γ ( t ) = ∂ Γ ( t ) and (cid:126)ν = (cid:126)ν on γ ( t ) . (2.2b)Of course, in the case (2.2b) it also holds that (cid:126)µ = − (cid:126)µ , where (cid:126)µ i denotes the outer conormal to Γ i ( t )on γ ( t ). 4 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
In order to formulate the governing problems in more detail, we denote by ∇ s = ( ∂ s , . . . , ∂ s d ) thesurface gradient on Γ i , and then define ∇ s (cid:126)χ = (cid:0) ∂ s j χ k (cid:1) dk,j =1 , as well as the Laplace–Beltrami operator∆ s = ∇ s . ∇ s = (cid:80) dj =1 ∂ s j . We then introduce the mean curvature vector as (cid:126) κ i = κ i (cid:126)ν i = ∆ s (cid:126) id on Γ i , (2.3)where (cid:126) id is the identity function on R d , and κ i is the mean curvature of Γ i , i.e. the sum of theprincipal curvatures of Γ i . In particular, the principal curvatures κ i,j , j = 1 , . . . , d −
1, togetherwith the eigenvalue zero for the eigenvector (cid:126)ν i , are the d eigenvalues of the symmetric linear map −∇ s (cid:126)ν i : R d → R d ; see e.g. [17, p. 152], where a different sign convention is used. The map −∇ s (cid:126)ν i isalso called the Weingarten map or shape operator. The mean curvature κ i and the Gaussian curvature K i of Γ i can now be stated as κ i = d − (cid:88) j =1 κ i,j = − tr( ∇ s (cid:126)ν i ) = −∇ s . (cid:126)ν i and K i = d − (cid:89) j =1 κ i,j . (2.4)Throughout the paper the main case we are interested in is d = 3, but it is often convenient to alsodiscuss the case d = 2 at the same time. To this end, we generalize the free energy (1.1) to E ((Γ i ( t )) i =1 ) = (cid:88) i =1 (cid:34) α i (cid:90) Γ i ( t ) ( κ i − κ i ) d H d − + α Gi (cid:90) Γ i ( t ) K i d H d − (cid:35) + ς H d − ( γ ( t )) , (2.5)where κ i and K i are the mean and Gaussian curvatures of Γ i ( t ), i = 1 , ς ∈ R ≥ denotes a possibleline tension, and α i ∈ R > and α Gi ∈ R denote the bending and Gaussian bending rigidities of Γ i ( t ), i = 1 ,
2, respectively. Here and throughout H k , k = 0 , ,
2, denotes the k -dimensional Hausdorffmeasure in R d .In the case d = 2, we always assume that ς = α G = α G = 0. For the case d = 3, on the other hand,we mention that the contributions (cid:88) i =1 (cid:34) α i (cid:90) Γ i ( t ) κ i d H + α Gi (cid:90) Γ i ( t ) K i d H (cid:35) (2.6)to the energy (2.5) are positive semidefinite with respect to the principal curvatures if α Gi ∈ [ − α i , , i = 1 , . (2.7)In the C –case, recall (2.2b), adding multiples of (cid:80) i =1 K i d H to the energy only changes the energyby a constant which follows from the Gauss–Bonnet theorem, see (2.12) below. Hence we obtain thatthe energy (2.5) can be bounded from below if α Gi ≥ max { α G , α G } − α i for i = 1 ,
2, which will holdwhenever min { α , α } ≥ | α G − α G | . (2.8)Variational problems for integrals including the energy (2.6) require that the energy is definite, seee.g. [33, p. 364], in order to be able to show a priori estimates. As discussed in [33], the condition ofdefiniteness leads to the constraints (2.7) and (2.8), and it is likely that these conditions also haveimplications for the existence and regularity theory of gradient flows for (2.5) in the C – and C –case,respectively.In the case d = 3, similarly to (2.3), fundamental to many approaches, which numerically approxi-mate evolving curves in a parametric way, is the identity (cid:126) id ss = (cid:126) κ γ on γ ( t ) , (2.9)where (cid:126) κ γ is the curvature vector on γ ( t ). Here we choose the arclength s of the curve γ ( t ) such that (cid:126)µ i = ( − i (cid:126)ν i × (cid:126) id s on γ ( t ) , (2.10)5 . W. Barrett, H. Garcke, R. N¨urnberg for i = 1 ,
2, denote the outer conormals to Γ i ( t ) on γ ( t ). Note that (cid:126)µ i is a vector that is perpendicularto the unit tangent (cid:126) id s on ∂ Γ i ( t ) and lies in the tangent space of Γ i ( t ). Now (2.9) can be rewritten as (cid:126) id ss = (cid:126) κ γ = ( (cid:126) κ γ . (cid:126)µ i ) (cid:126)µ i + ( (cid:126) κ γ . (cid:126)ν i ) (cid:126)ν i on γ ( t ) , (2.11)where (cid:126) κ γ . (cid:126)µ i is the geodesic curvature and (cid:126) κ γ . (cid:126)ν i is the normal curvature of γ ( t ) on Γ i ( t ), i = 1 ,
2. Itthen follows from the Gauss–Bonnet theorem, (cid:90) Γ i ( t ) K i d H = 2 π m (Γ i ( t )) + (cid:90) γ ( t ) (cid:126) κ γ . (cid:126)µ i d H , (2.12)where m (Γ i ( t )) ∈ Z denotes the Euler characteristic of Γ i ( t ), that the energy (2.5), is equivalent to E ((Γ i ( t )) i =1 ) = (cid:88) i =1 (cid:34) α i (cid:90) Γ i ( t ) ( κ i − κ i ) d H + α Gi (cid:34)(cid:90) γ ( t ) (cid:126) κ γ . (cid:126)µ i d H + 2 π m (Γ i ( t )) (cid:35)(cid:35) + ς H ( γ ( t )) . (2.13)We note that we use a sign for the conormal that is different from many authors in differential geometry,and hence we obtain a different sign in the Gauss–Bonnet formula.In some cases, in particular in applications for biomembranes, cf. [38], the surface areas of Γ ( t )and Γ ( t ) need to stay constant during the evolution, as well as the volume enclosed by Γ( t ). Here andthroughout we use the terminology “surface area” and “enclosed volume” also for the case d = 2, whenthe former is really curve length, and the latter means enclosed area. In this case one can consider E λ ((Γ i ( t )) i =1 ) = E ((Γ i ( t )) i =1 ) + λ V ( t ) L d (Ω( t )) + (cid:88) i =1 λ Ai ( t ) H d − (Γ i ( t )) , (2.14)where Ω( t ) denotes the interior of Γ( t ) and L d denotes the Lebesgue measure in R d . Here, λ Ai ( t ) areLagrange multipliers for the area constraints, which can be interpreted as a surface tension, and λ V ( t )is a Lagrange multiplier for the volume constraint which might be interpreted as a pressure difference.For the convenience of the reader, we end this section by stating the strong formulations of the L –gradient flows for (2.5) in the presence of the matching conditions (2.2a) and (2.2b), respectively.These strong formulations directly follow from the weak formulation introduced in Section 3, as weshow rigorously in the appendix.The weighted L –gradient flow, (1.2), of (2.13), for d = 2 or d = 3, then leads to the evolution law (cid:126) V . (cid:126)ν i = − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | on Γ i ( t ) , i = 1 , . (2.15)See (A.8) in the appendix for a derivation of (2.15). We remark that if the more general energy (2.14)is considered, then (2.15) is replaced by (cid:126) V . (cid:126)ν i = − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | + λ Ai κ i − λ V on Γ i ( t ) , (2.16)for i = 1 ,
2, see (A.13) in the appendix.In the case d = 3 we introduce the second fundamental form II i of Γ i ( t ), which is given asII i ( (cid:126) t ,(cid:126) t ) = − [ ∂ (cid:126) t (cid:126)ν i ] .(cid:126) t = − [( ∇ s (cid:126)ν i ) (cid:126) t ] .(cid:126) t on Γ i ( t ) , (2.17)for all tangential vectors (cid:126) t j , j = 1 ,
2. We note that II i ( · , · ) is a symmetric bilinear form, as ∇ s (cid:126)ν i issymmetric. In addition, we define τ i = II i ( (cid:126) id s , (cid:126)µ i ) on γ ( t ) , (2.18)i.e. τ i = − ( (cid:126)ν i ) s . (cid:126)µ i on γ ( t ). 6 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
Still considering the case d = 3, in the C –junction case, the boundary conditions on γ ( t ) are givenby α i ( κ i − κ i ) + α Gi (cid:126) κ γ . (cid:126)ν i = 0 on γ ( t ) , i = 1 , , (2.19a) (cid:88) i =1 (cid:2) (( α i ( ∇ s κ i ) . (cid:126)µ i − α Gi ( τ i ) s ) (cid:126)ν i − ( α i ( κ i − κ i ) + α Gi K i + λ Ai ) (cid:126)µ i (cid:3) + ς (cid:126) κ γ = (cid:37) (cid:126) V on γ ( t ) , (2.19b)see (A.12a), (A.14) in the appendix. We note that (2.19a) are two scalar conditions, while (2.19b) givesrise to two conditions as (cid:126)µ i , (cid:126)ν i and (cid:126) κ γ are all perpendicular to the tangent space to γ ( t ). ExpressingΓ and Γ locally as two graphs, we also obtain one condition for the height functions stemming fromthe C –condition. Altogether we have five conditions, as is to be expected for a free boundary probleminvolving fourth order operators on both sides of the free boundary. In this context we also refer toRemark 2.1 in [6].In the C –junction case, when (cid:126)ν = (cid:126)ν = (cid:126)ν and (cid:126)µ = (cid:126)µ = − (cid:126)µ on γ ( t ), the boundary conditions on γ ( t ) for the dissipation dynamics (1.2), with E replaced by E λ , are given by[ α i ( κ i − κ i )] + [ α Gi ] (cid:126) κ γ . (cid:126)ν = 0 on γ ( t ) , (2.20a)[ α i ( ∇ s κ i )] . (cid:126)µ + ς (cid:126) κ γ . (cid:126)ν − [ α Gi ] τ s = (cid:37) (cid:126) V . (cid:126)ν on γ ( t ) , (2.20b)[ − α i ( κ i − κ i ) + α i ( κ i − κ i ) ( κ i − (cid:126) κ γ . (cid:126)ν ) − λ Ai ] + [ α Gi ] τ + ς (cid:126) κ γ . (cid:126)µ = (cid:37) (cid:126) V . (cid:126)µ on γ ( t ) , (2.20c)where τ = τ = − τ is the geodesic torsion of the curve γ ( t ) on Γ( t ). We note that (2.20a–c), in thecase (cid:37) = 0, agree with (3.16)–(3.18) in [23], see also [24, (2.7b,a,c)]. In terms of counting the number ofequations, we see that (2.20a–c) are three conditions, together with one condition coming from (cid:126)ν = (cid:126)ν and one condition from the requirement that the two phases match up continuously, leading to fiveconditions in total. We refer to (A.15a), (A.24a,b) in the appendix for a derivation of (2.20a–c). Remark 2.1.
We note that although the conditions (2.19a,b) and (2.20a–c) were derived for the case d = 3, they are also valid in the case d = 2 on recalling that in this case we set ς = α G = α G = 0.In particular, (2.19a) then simplifies to κ i = κ i on γ ( t ), i = 1 ,
2, which is the same as the condition[6, (2.13c)] that was derived by the authors for a C –junction between two curves meeting in 2d. Inaddition, (2.19b) for d = 2 and (cid:37) = 0 collapses to [6, (2.13b)], modulo the different sign conventionemployed there.Similarly, (2.20a) for d = 2 simplifies to α ( κ − κ ) = α ( κ − κ ) on γ ( t ), which is the same asthe condition [6, (2.18e)], modulo the different sign convention employed there, that was derived bythe authors for a C –junction between two curves meeting in 2d. In addition, (2.20b,c) for d = 2 and (cid:37) = 0, collapse to [6, (2.18b,c)]. Weak formulation
In this section we derive a weak formulation of a generalized L –gradient flow of E ((Γ i ( t )) i =1 ). Theweak formulation of the standard L –gradient flow is given by (3.29), below, with (cid:37) = 0, where f Γ represents the first variation of the energy E ((Γ i ( t )) i =1 ) formulated in a suitable weak form. In whatfollows we will define a Lagrangian involving the energy and suitable constraints, which for examplerelate the curvatures to the parametrizations of the surfaces. This Lagrangian will allow us to derive(3.28a–f), below, which defines f Γ in a weak formulation involving only first order derivatives. Thisformulation will be suitable for a numerical approximation based on continuous, piecewise linear finiteelements, and such an approximation will be considered in Section 4.7 . W. Barrett, H. Garcke, R. N¨urnberg On recalling (2.1), we define the following time derivative that follows the parameterization (cid:126)x ( · , t )of Γ( t ). Let ( ∂ ◦ t φ ) | Γ i ( t ) = ( φ t + (cid:126) V . ∇ φ ) | Γ i ( t ) ∀ φ ∈ H (Γ i,T ) , (3.1)where we have defined the space-time surfacesΓ i,T := (cid:91) t ∈ [0 ,T ] Γ i ( t ) × { t } , i = 1 , , and Γ T := (cid:91) t ∈ [0 ,T ] Γ( t ) × { t } . Here we stress that (3.1) is well-defined, even though φ t and ∇ φ do not make sense separately for afunction φ ∈ H (Γ i,T ). We note thatdd t (cid:104) ψ i , φ i (cid:105) Γ i ( t ) = (cid:104) ∂ ◦ t ψ i , φ (cid:105) Γ i ( t ) + (cid:104) ψ i , ∂ ◦ t φ i (cid:105) Γ i ( t ) + (cid:68) ψ i φ i , ∇ s . (cid:126) V (cid:69) Γ i ( t ) ∀ ψ i , φ i ∈ H (Γ i,T ) , (3.2)see Lemma 5.2 in [21]. Here (cid:104)· , ·(cid:105) Γ i ( t ) denotes the L –inner product on Γ i ( t ), and (cid:104)· , ·(cid:105) Γ( t ) = (cid:80) i =1 (cid:104)· , ·(cid:105) Γ i ( t ) . It immediately follows from (3.2) thatdd t H d − (Γ i ( t )) = (cid:68) ∇ s . (cid:126) V , (cid:69) Γ i ( t ) = (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) V (cid:69) Γ i ( t ) . (3.3)Moreover, on recalling Lemma 2.1 from [17], it holds thatdd t L d (Ω( t )) = (cid:88) i =1 (cid:68) (cid:126) V , (cid:126)ν i (cid:69) Γ i ( t ) . (3.4)In this section we would like to derive a weak formulation for the L –gradient flow of E ((Γ i ( t )) i =1 ).To this end, we need to consider variations of the energy with respect to Γ( t ) = (cid:126)x (Υ , t ). Let H γ (Γ( t )) := { η ∈ L (Γ( t )) : η | Γ i ( t ) ∈ H (Γ i ( t )) , i = 1 , , ( η | Γ ( t ) ) | γ ( t ) = ( η | Γ ( t ) ) | γ ( t ) =: η | γ ( t ) ∈ H ( γ ( t )) } . In addition, for any given (cid:126)χ ∈ [ H γ (Γ( t ))] d and for any ε ∈ (0 , ε ) for some ε ∈ R > , letΓ ε ( t ) := { (cid:126) Ψ( (cid:126)z, ε ) : (cid:126)z ∈ Γ( t ) } , where (cid:126) Ψ( (cid:126)z,
0) = (cid:126)z and ∂ (cid:126) Ψ ∂ε ( (cid:126)z,
0) = (cid:126)χ ( (cid:126)z ) ∀ (cid:126)z ∈ Γ( t ) . (3.5)Of course, we have that Γ ε ( t ) = Γ ,ε ( t ) ∪ γ ε ( t ) ∪ Γ ,ε ( t ), whereΓ i,ε ( t ) := { (cid:126) Ψ( (cid:126)z, ε ) : (cid:126)z ∈ Γ i ( t ) } , i = 1 , , and γ ε ( t ) = ∂ Γ ,ε ( t ) = ∂ Γ ,ε ( t ) . Similarly to (3.3), the first variation of H d − (Γ i ( t )) with respect to Γ( t ) in the direction (cid:126)χ ∈ [ H γ (Γ( t ))] d is given by (cid:20) δδ Γ H d − (Γ i ( t )) (cid:21) ( (cid:126)χ ) = dd ε H d − (Γ i,ε ( t )) | ε =0 = lim ε → ε (cid:104) H d − (Γ i,ε ( t )) − H d − (Γ i ( t )) (cid:105) = (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) , (3.6)see e.g. the proof of Lemma 1 in [20].In order to derive a suitable weak formulation, we formally consider the first variation of (2.5)subject to the following side constraint, which is inspired by the weak formulation of (2.3), (cid:68) Q i,θ (cid:126) κ (cid:63)i , (cid:126)η (cid:69) Γ i ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ i ( t ) = (cid:104) (cid:126) m i , (cid:126)η (cid:105) γ ( t ) ∀ (cid:126)η ∈ [ H (Γ i ( t ))] d , i = 1 , , (3.7)where θ ∈ [0 ,
1] is a fixed parameter, and where Q i,θ are defined by Q i,θ = θ Id + (1 − θ ) (cid:126)ν i ⊗ (cid:126)ν i on Γ i ( t ) . (3.8)8 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
Of course, (3.7) holds trivially on the continuous level for (cid:126) κ (cid:63)i = (cid:126) κ i and for (cid:126) m i being the conormal (cid:126)µ i , independently of the choice of θ ∈ [0 , θ = 1. However, under discretization that formulation wouldlead to undesirable mesh effects. Hence, in line with the authors previous work in [10], we also allow θ ∈ [0 , θ = 0,in general. In rare cases we may need to dampen the tangential motion that occurs in the case θ = 0.To this end, we allow for the full range of values θ ∈ [0 , (cid:10) (cid:126) κ (cid:63)γ , (cid:126)η (cid:11) γ ( t ) + (cid:68) (cid:126) id s , (cid:126)η s (cid:69) γ ( t ) = 0 ∀ (cid:126)η ∈ [ H ( γ ( t ))] d . (3.9)Finally, in order to model a C – or C –contact we require C ( (cid:126) m + (cid:126) m ) = (cid:126) γ ( t ) , (3.10)where C = 0 for C and C = 1 for C .We now define the Lagrangian L ((Γ i ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) = (cid:88) i =1 (cid:104) α i (cid:104) (cid:126) κ (cid:63)i − κ i (cid:126)ν i , (cid:126) κ (cid:63)i − κ i (cid:126)ν i (cid:105) Γ i ( t ) + α Gi (cid:10) (cid:126) κ (cid:63)γ , (cid:126) m i (cid:11) γ ( t ) (cid:105) + ς H d − ( γ ( t )) − (cid:10) (cid:126) κ (cid:63)γ , (cid:126)z (cid:11) γ ( t ) − (cid:68) (cid:126) id s , (cid:126)z s (cid:69) γ ( t ) + C (cid:68) (cid:126) m + (cid:126) m , (cid:126)φ (cid:69) γ ( t ) − (cid:88) i =1 (cid:20)(cid:68) Q i,θ (cid:126) κ (cid:63)i , (cid:126)y i (cid:69) Γ i ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)y i (cid:69) Γ i ( t ) − (cid:104) (cid:126) m i , (cid:126)y i (cid:105) γ ( t ) (cid:21) , where (cid:126)y i ∈ [ H (Γ i ( t ))] d and (cid:126)z ∈ [ H ( γ ( t ))] d are Lagrange multipliers for (3.7) and (3.9), respectively.Similarly, (cid:126)φ ∈ [ L ( γ ( t ))] d is a Lagrange multiplier for (3.10). We now want to compute the firstvariation f Γ of E ((Γ i ( t )) i =1 ), subject to the side constraints (3.7), (3.9) and (3.10). This means that f Γ needs to fulfill f Γ ( (cid:126)χ ) = − (cid:20) δδ Γ E ( t ) (cid:21) ( (cid:126)χ ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d . (3.11)In particular, on using ideas from the formal calculus of PDE constrained optimization, see e.g. [37],we can formally compute f Γ by requiring that (cid:20) δδ Γ L (cid:21) ( (cid:126)χ ) = lim ε → ε (cid:104) ( L (Γ i,ε ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) − L ((Γ i ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) (cid:105) = − f Γ ( (cid:126)χ ) (3.12a) (cid:20) δδ(cid:126) κ (cid:63) L (cid:21) ( (cid:126)ξ ) = lim ε → ε (cid:104) L (Γ ( t ) , (cid:126) κ (cid:63) + ε (cid:126)ξ , (cid:126) m , (cid:126)y , Γ ( t ) , (cid:126) κ (cid:63) , (cid:126) m , (cid:126)y , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) − L ((Γ i ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) (cid:105) = 0 , (3.12b) (cid:20) δδ (cid:126) m L (cid:21) ( (cid:126)ζ ) = lim ε → ε (cid:104) L (Γ ( t ) , (cid:126) κ (cid:63) , (cid:126) m + ε (cid:126)ζ , (cid:126)y , Γ ( t ) , (cid:126) κ (cid:63) , (cid:126) m , (cid:126)y , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) − L ((Γ i ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) (cid:105) = 0 , (3.12c) (cid:20) δδ(cid:126)y L (cid:21) ( (cid:126)η ) = lim ε → ε (cid:104) L (Γ ( t ) , (cid:126) κ (cid:63) , (cid:126) m , (cid:126)y + ε (cid:126)η , Γ ( t ) , (cid:126) κ (cid:63) , (cid:126) m , (cid:126)y , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ )9 . W. Barrett, H. Garcke, R. N¨urnberg − L ((Γ i ( t ) , (cid:126) κ (cid:63)i , (cid:126) m i , (cid:126)y i ) i =1 , (cid:126) κ (cid:63)γ , (cid:126)z, (cid:126)φ ) (cid:105) = 0 , (3.12d)for variations (cid:126)χ ∈ [ H γ (Γ( t ))] d , (cid:126)ξ ∈ [ L (Γ ( t ))] d , (cid:126)ζ ∈ [ L ( γ ( t ))] d and (cid:126)η ∈ [ L (Γ ( t ))] d ; and similarlyfor the variations for (cid:126) κ (cid:63) , (cid:126) m , (cid:126)y , (cid:126) κ (cid:63)γ , (cid:126)z and (cid:126)φ .In order to calculate (3.12a–d), we note that generalized variants of (3.6) also hold. Namely, wehave that (cid:20) δδ Γ (cid:104) w i , (cid:105) Γ i ( t ) (cid:21) ( (cid:126)χ ) = dd ε (cid:104) w i,ε , (cid:105) Γ i,ε ( t ) | ε =0 = (cid:68) w i ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) ∀ w i ∈ L ∞ (Γ i ( t )) , (3.13)where w i,ε ∈ L ∞ (Γ i,ε ( t )), for any w i ∈ L ∞ (Γ i ( t )), is defined by w i,ε ( (cid:126) Ψ( (cid:126)z, ε )) = w i ( (cid:126)z ) ∀ (cid:126)z ∈ Γ i ( t ) , and similarly for (cid:126)w ∈ [ L ∞ (Γ i ( t ))] d . This definition of w i,ε yields that ∂ ε w i = 0, where ∂ ε w i ( (cid:126)z ) = dd ε w i,ε ( (cid:126) Ψ( (cid:126)z, ε )) | ε =0 ∀ (cid:126)z ∈ Γ i ( t ) . (3.14)Of course, (3.13) is the first variation analogue of (3.2) with w i = ψ i φ i and ∂ ε ψ i = ∂ ε φ i = 0. Similarly,it holds that (cid:20) δδ Γ (cid:104) (cid:126)w i , (cid:126)ν i (cid:105) Γ i ( t ) (cid:21) ( (cid:126)χ ) = dd ε (cid:104) (cid:126)w i,ε , (cid:126)ν i,ε (cid:105) Γ i,ε ( t ) | ε =0 = (cid:68) ( (cid:126)w i . (cid:126)ν i ) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) + (cid:10) (cid:126)w i , ∂ ε (cid:126)ν i (cid:11) Γ i ( t ) ∀ (cid:126)w i ∈ [ L ∞ (Γ i ( t ))] d , (3.15)where ∂ ε (cid:126)w i = (cid:126) (cid:126)ν i,ε ( t ) denotes the unit normal on Γ i,ε ( t ). Moreover, we will make use of thefollowing result concerning the variation of (cid:126)ν i , with respect to Γ( t ), in the direction (cid:126)χ ∈ [ H γ (Γ( t ))] d : ∂ ε (cid:126)ν i = − [ ∇ s (cid:126)χ ] T (cid:126)ν i on Γ i ( t ) ⇒ ∂ ◦ t (cid:126)ν i = − [ ∇ s (cid:126) V ] T (cid:126)ν i on Γ i ( t ) , (3.16)see [35, Lemma 9]. We also note that for (cid:126)η i ∈ [ H (Γ i ( t ))] d it holds that (cid:20) δδ Γ (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η i (cid:69) Γ i ( t ) (cid:21) ( (cid:126)χ ) = dd ε (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η i,ε (cid:69) Γ i,ε ( t ) | ε =0 = (cid:104)∇ s . (cid:126)η i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) + d (cid:88) l, m =1 (cid:104) (cid:104) ( (cid:126)ν i ) l ( (cid:126)ν i ) m ∇ s ( (cid:126)η i ) m , ∇ s ( (cid:126)χ ) l (cid:105) Γ i ( t ) − (cid:104) ( ∇ s ) m ( (cid:126)η i ) l , ( ∇ s ) l ( (cid:126)χ ) m (cid:105) Γ i ( t ) (cid:105) = (cid:104)∇ s (cid:126)η i , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . (cid:126)η i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)η i ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) , (3.17)where ∂ ε (cid:126)η i = (cid:126)
0, see Lemma 2 and the proof of Lemma 3 in [20]. Here D ( (cid:126)χ ) := ∇ s (cid:126)χ + ( ∇ s (cid:126)χ ) T , and we note that our notation is such that ∇ s (cid:126)χ = ( ∇ Γ (cid:126)χ ) T , with ∇ Γ (cid:126)χ = ( ∂ s l χ m ) dl,m =1 defined as in[20]. It follows from (3.17) thatdd t (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ i ( t ) = (cid:68) ∇ s (cid:126)η, ∇ s (cid:126) V (cid:69) Γ i ( t ) + (cid:68) ∇ s . (cid:126)η, ∇ s . (cid:126) V (cid:69) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)η ) T , D ( (cid:126) V ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) ∀ (cid:126)η ∈ { (cid:126)ξ ∈ H (Γ i,T ) : ∂ ◦ t (cid:126)ξ = (cid:126) } . (3.18)Similarly to (3.13) it holds that (cid:20) δδ Γ (cid:104) w, (cid:105) γ ( t ) (cid:21) ( (cid:126)χ ) = dd ε (cid:104) w ε , (cid:105) γ ε ( t ) | ε =0 = (cid:68) w (cid:126) id s , (cid:126)χ s (cid:69) γ ( t ) ∀ w ∈ L ∞ ( γ ( t )) , (cid:126)χ ∈ [ H γ (Γ( t ))] d , (3.19)10 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES where ∂ ε w = 0. Moreover, similarly to (3.17), we note that for (cid:126)η ∈ [ H γ (Γ( t ))] d it holds that (cid:20) δδ Γ (cid:68) (cid:126) id s , (cid:126)η s (cid:69) γ ( t ) (cid:21) ( (cid:126)χ ) = (cid:10) P γ (cid:126)η s , (cid:126)χ s (cid:11) γ ( t ) , (3.20)where ∂ ε (cid:126)η = (cid:126)
0, and where P γ = Id − (cid:126) id s ⊗ (cid:126) id s on γ ( t ) . (3.21)Now combining (3.12a–d), on noting (3.13)–(3.21), yields that f Γ ( (cid:126)χ ) = (cid:88) i =1 (cid:20) (cid:104)∇ s (cid:126)y i , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . (cid:126)y i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)y i ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) − (cid:68) [ α i | (cid:126) κ (cid:63)i − κ i (cid:126)ν i | − (cid:126) κ (cid:63)i . Q i,θ (cid:126)y i )] ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) + α i κ i (cid:10) (cid:126) κ (cid:63)i , ∂ ε (cid:126)ν i (cid:11) Γ i ( t ) + (cid:68) ∂ ε [ Q i,θ (cid:126) κ (cid:63)i ] , (cid:126)y i (cid:69) Γ i ( t ) (cid:21) − ς (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ ( t ) + (cid:42) (cid:126) κ (cid:63)γ . (cid:126)z − C ( (cid:126) m + (cid:126) m ) . (cid:126)φ − (cid:88) i =1 ( α Gi (cid:126) κ (cid:63)γ + (cid:126)y i ) . (cid:126) m i , (cid:126) id s . (cid:126)χ s (cid:43) γ ( t ) + (cid:10) P γ (cid:126)z s , (cid:126)χ s (cid:11) γ ( t ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d , (3.22a) α i ( (cid:126) κ (cid:63)i − κ i (cid:126)ν i ) − Q i,θ (cid:126)y i = (cid:126) i ( t ) , i = 1 , , (3.22b) α Gi (cid:126) κ (cid:63)γ + (cid:126)y i + C (cid:126)φ = (cid:126) γ ( t ) , i = 1 , , (3.22c) (cid:88) i =1 α Gi (cid:126) m i − (cid:126)z = (cid:126) γ ( t ) , i = 1 , , (3.22d)with (3.7), (3.10) and (3.9). As ∂ ε (cid:126) κ (cid:63)i = (cid:126)
0, we have that ∂ ε [ Q i,θ (cid:126) κ (cid:63)i ] = (1 − θ ) (cid:2) ( (cid:126) κ (cid:63)i . ∂ ε (cid:126)ν i ) (cid:126)ν i + ( (cid:126) κ (cid:63)i . (cid:126)ν i ) ∂ ε (cid:126)ν i (cid:3) . (3.23)We observe that (3.22b,c) imply that Q i,θ (cid:126)y i = α i (cid:126) κ (cid:63)i − α i κ i (cid:126)ν i on Γ i ( t ) and (cid:126)y i + C (cid:126)φ = − α Gi (cid:126) κ (cid:63)γ on γ ( t ) . (3.24)Let us now recover (cid:126) κ (cid:63)i and (cid:126) κ (cid:63)γ in terms of the geometry again. To this end, we first recall the identity (cid:90) Γ i ( t ) ∇ s g d H d − = − (cid:90) Γ i ( t ) g κ i (cid:126)ν i d H d − + (cid:90) γ ( t ) g (cid:126)µ i d H d − ∀ g ∈ H (Γ i ( t )) , (3.25)see e.g. Theorem 2.10 in [21] and Proposition 4.5 in [36, p. 334]. It immediately follows from (3.7),(2.3) and (3.25) that (cid:126) m i = (cid:126)µ i and Q i,θ (cid:126) κ (cid:63)i = (cid:126) κ i = κ i (cid:126)ν i , with the latter implying that (cid:126) κ (cid:63)i . (cid:126)ν i = κ i . (3.26)Hence we immediately get (cid:126) κ (cid:63)i = (cid:126) κ i for θ ∈ (0 , θ = 0, on the other hand, it follows from (3.24)and (3.26) that α i (cid:126) κ (cid:63)i = [ (cid:126)y i . (cid:126)ν i + α i κ i ] (cid:126)ν i , and so (cid:126) κ (cid:63)i = κ i (cid:126)ν i = (cid:126) κ i . Moreover, combining (3.9) and (2.9)yields that (cid:126) κ (cid:63)γ = (cid:126) κ γ . Overall, we obtain from (3.24) that Q i,θ (cid:126)y i = α i ( κ i − κ i ) (cid:126)ν i on Γ i ( t ) and (cid:126)y i + C (cid:126)φ = − α Gi (cid:126) κ γ on γ ( t ) . (3.27)However, if θ ∈ (0 , α Gi (cid:54) = 0, since thefirst condition in (3.27) yields that (cid:126)y i = α i ( κ i − κ i ) (cid:126)ν i . If C = 1 then the two conditions are in generalincompatible even if α Gi = 0. Hence for general boundaries γ ( t ) and α Gi (cid:54) = 0 we need to take θ = 0, atleast locally at the boundary. Therefore it may be desirable to consider a variable θ ∈ L ∞ (Γ( t )). The11 . W. Barrett, H. Garcke, R. N¨urnberg calculation (3.22a–d) remains valid provided that ∂ ε θ = 0. We will make this more rigorous on thediscrete level, see (4.8) below.Using (3.16), (3.23) and (3.22c,d) in (3.22a) yields the condensed version f Γ ( (cid:126)χ ) = (cid:88) i =1 (cid:20) (cid:104)∇ s (cid:126)y i , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . (cid:126)y i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)y i ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) − (cid:68) [ α i | (cid:126) κ i − κ i (cid:126)ν i | − (cid:126) κ i . Q i,θ (cid:126)y i )] ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) − α i κ i (cid:10) (cid:126) κ i , [ ∇ s (cid:126)χ ] T (cid:126)ν i (cid:11) Γ i ( t ) − (1 − θ ) (cid:10)(cid:2) ( (cid:126) κ i . [ ∇ s (cid:126)χ ] T (cid:126)ν i ) (cid:126)ν i + ( (cid:126) κ i . (cid:126)ν i ) [ ∇ s (cid:126)χ ] T (cid:126)ν i (cid:3) , (cid:126)y i (cid:11) Γ i ( t ) (cid:105) − ς (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ ( t ) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126) κ γ . (cid:126) m i , (cid:126) id s . (cid:126)χ s (cid:69) γ ( t ) + (cid:10) P γ ( (cid:126) m i ) s , (cid:126)χ s (cid:11) γ ( t ) (cid:21) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d , (3.28a) α i ( (cid:126) κ i − κ i (cid:126)ν i ) − Q i,θ (cid:126)y i = (cid:126) i ( t ) , i = 1 , , (3.28b) α Gi (cid:126) κ γ + (cid:126)y i + C (cid:126)φ = (cid:126) γ ( t ) , i = 1 , , (3.28c) (cid:68) Q i,θ (cid:126) κ i , (cid:126)η (cid:69) Γ i ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ i ( t ) = (cid:104) (cid:126) m i , (cid:126)η (cid:105) γ ( t ) ∀ (cid:126)η ∈ [ H (Γ i ( t ))] d , i = 1 , , (3.28d) C ( (cid:126) m + (cid:126) m ) = (cid:126) γ ( t ) , (3.28e) (cid:104) (cid:126) κ γ , (cid:126)η (cid:105) γ ( t ) + (cid:68) (cid:126) id s , (cid:126)η s (cid:69) γ ( t ) = 0 ∀ (cid:126)η ∈ [ H ( γ ( t ))] d . (3.28f) Remark 3.1.
We recall from (3.27) and the discussion below that in general we require θ = 0. If C = 0 then it follows from (3.28c) that (cid:126)y i = − α Gi (cid:126) κ γ on γ ( t ), for i = 1 ,
2. Combining this with (3.28b)for θ = 0 then yields that (2.19a) holds.On the other hand, in the case of a C –junction, when C = 1, then (3.28e) implies that (cid:126)µ + (cid:126)µ = (cid:126) (cid:126)ν = (cid:126)ν = (cid:126)ν on γ ( t ), and so it follows from (3.28b,c) with θ = 0 that α i ( κ i − κ i ) + α Gi (cid:126) κ γ . (cid:126)ν + (cid:126)φ . (cid:126)ν = 0 on γ ( t ) , i = 1 , , which means that (2.20a) holds.The weak formulation of a generalized L –gradient flow of E ((Γ i ( t )) i =1 ) can then be formulated asfollows. Given Γ i (0), i = 1 ,
2, for all t ∈ (0 , T ] find Γ i ( t ) = (cid:126)x i (Υ i , t ), i = 1 ,
2, with (cid:126) V ( t ) ∈ [ H (Γ( t ))] d ,and (cid:126) κ i ( t ) ∈ [ L (Γ i ( t ))] d , (cid:126)y i ( t ) ∈ [ H (Γ i ( t ))] d , (cid:126) m i ( t ) ∈ [ H ( γ ( t ))] d , i = 1 ,
2, as well as (cid:126) κ γ ∈ [ L ( γ ( t ))] d , (cid:126)z ∈ [ L ( γ ( t ))] d , (cid:126)φ ∈ [ L ( γ ( t ))] d such that (cid:68) (cid:126) V , (cid:126)χ (cid:69) Γ( t ) + (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = f Γ ( (cid:126)χ ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d (3.29)and (3.28a–f) hold. Here we note that (cid:37) = 0 recovers a weak formulation for the standard L –gradientflow. As stated in (1.2), we allow for (cid:37) ≥ γ ( t ). In numerical simulations such a damping often proves beneficial, as it suppressespossible oscillations at the contact line. On the other hand, such a dissipation mechanism at theboundary is probably also relevant in applications. Semidiscrete finite element approximation
It is the aim of this section to introduce a semidiscrete continuous-in-time finite element approxi-mation of the weak formulation (3.29), (3.28a–f) derived in the previous section. Our finite element12
RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES discretization will be given by (4.27a–f) below, and the main result of this section is the stability proofin Theorem 4.1 below.Similarly to [3], we introduce the following discrete spaces. Let Γ h ( t ) ⊂ R d be ( d − polyhedral surfaces , i.e. unions of non-degenerate ( d − d = 3), approximating the surfaces Γ( t ). In particular, let Γ h ( t ) = (cid:83) Jj =1 σ hj ( t ), where { σ hj ( t ) } Jj =1 is a family mutually disjoint open ( d − { (cid:126)q hk ( t ) } Kk =1 . In analogy tothe continuous setting, we write Γ h ( t ) = Γ h ( t ) ∪ γ h ( t ) ∪ Γ h ( t ), where γ h ( t ) = ∂ Γ h ( t ) = ∂ Γ h ( t ). Here welet Γ hi ( t ) = (cid:83) J i j =1 σ hi,j ( t ), with vertices { (cid:126)q hi,k ( t ) } K i k =1 , i = 1 ,
2. We also assume that γ h ( t ) has the vertices { (cid:126)q hγ,k ( t ) } K γ k =1 . Clearly, it holds that J = J + J and K = K + K − K γ . Then let V h (Γ hi ( t )) = { (cid:126)χ ∈ [ C (Γ hi ( t ))] d : (cid:126)χ | σ hi,j is linear ∀ j = 1 , . . . , J i } = [ W h (Γ hi ( t ))] d , i = 1 , , where W h (Γ hi ( t )) is the space of scalar continuous piecewise linear functions on Γ hi ( t ), with { χ hi,k ( · , t ) } K i k =1 denoting the standard basis of W h (Γ hi ( t )), i.e. χ hi,k ( (cid:126)q hi,l ( t ) , t ) = δ kl ∀ k, l ∈ { , . . . , K i } , t ∈ [0 , T ] . (4.1)In addition, let V h (Γ h ( t )) = { (cid:126)χ ∈ [ C (Γ h ( t ))] d : (cid:126)χ | Γ hi ( t ) ∈ V h (Γ hi ( t )) , i = 1 , } = [ W h (Γ h ( t ))] d . We denote the basis functions of W h (Γ h ( t )) by { χ hk ( · , t ) } Kk =1 . Moreover, let V h ( γ h ( t )) := { (cid:126)ψ ∈ [ C ( γ h ( t ))] d : ∃ (cid:126)χ ∈ V h (Γ h ( t )) (cid:126)χ | γ h ( t ) = (cid:126)ψ } =: [ W h ( γ h ( t ))] d , (4.2a) V h (Γ h ( t )) := { (cid:126)χ ∈ V h (Γ h ( t )) : (cid:126)χ | γ h ( t ) = (cid:126) } , (4.2b) V h (Γ hi ( t )) := { (cid:126)χ ∈ V h (Γ hi ( t )) : (cid:126)χ | γ h ( t ) = (cid:126) } . (4.2c)We denote the basis functions of W h ( γ h ( t )) by { φ hk ( · , t ) } K γ k =1 . We require that Γ hi ( t ) = (cid:126)X h (Γ hi (0) , t )with (cid:126)X h ∈ V h (Γ h (0)), and that (cid:126)q hk ∈ [ H (0 , T )] d , k = 1 , . . . , K .We denote the L –inner products on Γ h ( t ), Γ hi ( t ) and and γ h ( t ) by (cid:104)· , ·(cid:105) Γ h ( t ) , (cid:104)· , ·(cid:105) Γ hi ( t ) and (cid:104)· , ·(cid:105) γ h ( t ) ,respectively. In addition, for piecewise continuous functions, with possible jumps across the edges of { σ hi,j } J i j =1 , we also introduce the mass lumped inner product (cid:104) η, φ (cid:105) h Γ hi ( t ) = J i (cid:88) j =1 (cid:104) η, φ (cid:105) hσ hi,j ( t ) := J i (cid:88) j =1 1 d H d − ( σ hi,j ( t )) d (cid:88) k =1 ( η φ )(( (cid:126)q hi,j k ( t )) − ) , where { (cid:126)q hi,j k ( t ) } dk =1 are the vertices of σ hi,j ( t ), and where we define η (( (cid:126)q hi,j k ( t )) − ) := lim σ hj ( t ) (cid:51) (cid:126)p → (cid:126)q hi,jk ( t ) η ( (cid:126)p ).We naturally extend this definition to vector and tensor functions. We also define the mass lumpedinner products (cid:104)· , ·(cid:105) h Γ h ( t ) and (cid:104)· , ·(cid:105) hγ h ( t ) in the obvious way.Let (cid:126)ν hi denote the the outward unit normal to Γ hi ( t ), i = 1 ,
2, and similarly let (cid:126)ν h denote the theoutward unit normal to Γ h ( t ). Then we introduce the vertex normal functions (cid:126)ω hi ( · , t ) ∈ V h (Γ hi ( t ))with (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) := 1 H d − (Λ hi,k ( t )) (cid:88) j ∈ Θ hi,k H d − ( σ hi,j ( t )) (cid:126)ν hi | σ hi,j ( t ) , (4.3)where for k = 1 , . . . , K i we define Θ hi,k := { j : (cid:126)q hi,k ( t ) ∈ σ hi,j ( t ) } and set Λ hi,k ( t ) := ∪ j ∈ Θ hi,k σ hi,j ( t ). Herewe note that (cid:68) (cid:126)z, w (cid:126)ν hi (cid:69) h Γ hi ( t ) = (cid:68) (cid:126)z, w (cid:126)ω hi (cid:69) h Γ hi ( t ) ∀ (cid:126)z ∈ V h (Γ hi ( t )) , w ∈ W h (Γ hi ( t )) . (4.4)13 . W. Barrett, H. Garcke, R. N¨urnberg In the analogous fashion, we introduce the vertex normal function (cid:126)ω h ( · , t ) ∈ V h (Γ h ( t )), i.e. we set (cid:126)ω h ( (cid:126)q hk ( t ) , t ) := 1 H d − (Λ hk ( t )) (cid:88) j ∈ Θ hk H d − ( σ hj ( t )) (cid:126)ν h | σ hj ( t ) , (4.5)where for k = 1 , . . . , K we define Θ hk := { j : (cid:126)q hk ( t ) ∈ σ hj ( t ) } and set Λ hk ( t ) := ∪ j ∈ Θ hk σ hj ( t ). Of course, itholds that (cid:68) (cid:126)z, w (cid:126)ν h (cid:69) h Γ h ( t ) = (cid:68) (cid:126)z, w (cid:126)ω h (cid:69) h Γ h ( t ) ∀ (cid:126)z ∈ V h (Γ h ( t )) , w ∈ W h (Γ h ( t )) . (4.6)It clearly follows from (4.4) and (4.6) that (cid:68) (cid:126)z, (cid:126)ω h (cid:69) h Γ h ( t ) = (cid:88) i =1 (cid:68) (cid:126)z, (cid:126)ω hi (cid:69) h Γ hi ( t ) ∀ (cid:126)z ∈ V h (Γ h ( t )) . (4.7)In addition, for a given parameter θ ∈ [0 ,
1] we introduce θ h ∈ W h (Γ h ( t )) and θ h(cid:63) ∈ W h (Γ h ( t )) suchthat θ h ( (cid:126)q hk ( t ) , t ) = (cid:40) (cid:126)q hk ( t ) ∈ γ h ( t ) ,θ (cid:126)q hk ( t ) (cid:54)∈ γ h ( t ) , and θ h(cid:63) ( (cid:126)q hk ( t ) , t ) = (cid:40) (cid:126)q hk ( t ) ∈ γ h ( t ) ,θ (cid:126)q hk ( t ) (cid:54)∈ γ h ( t ) . (4.8)Then, similarly to (3.8), we introduce Q hi,θ h ∈ [ W h (Γ hi ( t ))] d × d and Q hi,θ h(cid:63) ∈ [ W h (Γ hi ( t ))] d × d by setting,for k ∈ { , . . . , K i } , Q hi,θ h ( (cid:126)q hi,k ( t ) , t ) = θ h ( (cid:126)q hi,k ( t ) , t ) Id + (1 − θ h ( (cid:126)q hi,k ( t ) , t )) (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) ⊗ (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) | (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) | , (4.9)and similarly for Q hi,θ h(cid:63) , where here and throughout we assume that (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) (cid:54) = (cid:126) k = 1 , . . . , K i and t ∈ [0 , T ]. Only in pathological cases could this assumption be violated, and in practice this neveroccurred. We note that (cid:68) Q hi,θ h (cid:126)z, (cid:126)v (cid:69) h Γ hi ( t ) = (cid:68) (cid:126)z, Q hi,θ h (cid:126)v (cid:69) h Γ hi ( t ) and (cid:68) Q hi,θ h (cid:126)z, (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:68) (cid:126)z, (cid:126)ω hi (cid:69) h Γ hi ( t ) (4.10)for all (cid:126)z, (cid:126)v ∈ V h (Γ hi ( t )), and analogously for θ h in (4.10) replaced by θ h(cid:63) . In addition, similarly to (4.7),it holds that (cid:88) i =1 (cid:68) (cid:126)z, Q hi,θ h(cid:63) (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:68) (cid:126)z, Q hi,θ h(cid:63) (cid:126)ω h (cid:69) h Γ hi ( t ) ∀ (cid:126)z ∈ V h (Γ h ( t )) . (4.11)Following the approach in the continuous setting, recall (2.13), (3.7), (3.9), we consider the firstvariation of the discrete energy E h ((Γ hi ( t )) i =1 ) := (cid:88) i =1 (cid:20) α i (cid:68) | (cid:126)κ hi − κ i (cid:126)ν hi | , (cid:69) h Γ hi ( t ) + α Gi (cid:20)(cid:68) (cid:126)κ hγ , (cid:126) m hi (cid:69) hγ h ( t ) + 2 π m (Γ hi ( t )) (cid:21)(cid:21) + ς H d − ( γ h ( t )) , (4.12)where (cid:126)κ hi ∈ V h (Γ hi ( t )), (cid:126) m hi ∈ V h ( γ h ( t )), i = 1 ,
2, and (cid:126)κ hγ ∈ V h ( γ h ( t )), subject to the side constraints (cid:68) Q hi,θ h (cid:126)κ hi , (cid:126)η (cid:69) h Γ hi ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ hi ( t ) = (cid:68) (cid:126) m hi , (cid:126)η (cid:69) hγ h ( t ) ∀ (cid:126)η ∈ V h (Γ hi ( t )) , i = 1 , , (4.13a) (cid:68) (cid:126)κ hγ , (cid:126)χ (cid:69) hγ h ( t ) + (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ h ( t ) = 0 ∀ (cid:126)χ ∈ V h ( γ h ( t )) , (4.13b) C ( (cid:126) m h + (cid:126) m h ) = (cid:126) γ h ( t ) . (4.13c)14 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
In particular, we define the Lagrangian L h ( t ) = (cid:88) i =1 (cid:20) α i (cid:68) | (cid:126)κ hi − κ i (cid:126)ν hi | , (cid:69) h Γ hi ( t ) + α Gi (cid:68) (cid:126)κ hγ , (cid:126) m hi (cid:69) hγ h ( t ) (cid:21) + ς H d − ( γ h ( t )) − (cid:68) (cid:126)κ hγ , (cid:126)Z h (cid:69) hγ h ( t ) − (cid:68) (cid:126) id s , (cid:126)Z hs (cid:69) γ h ( t ) − (cid:88) i =1 (cid:20)(cid:68) ( Q hi,θ h (cid:126)κ hi , (cid:126)Y hi (cid:69) h Γ hi ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)Y hi (cid:69) Γ hi ( t ) − (cid:68) (cid:126) m hi , (cid:126)Y hi (cid:69) hγ h ( t ) (cid:21) + C (cid:68) (cid:126) m h + (cid:126) m h , (cid:126) Φ h (cid:69) hγ h ( t ) , (4.14)where (cid:126)κ hi ∈ V h (Γ hi ( t )) and (cid:126)κ hγ ∈ V h ( γ h ( t )) satisfy (4.13a) and (4.13b), respectively, with (cid:126)Y hi ∈ V h (Γ hi ( t )) and (cid:126)Z h ∈ V h ( γ h ( t )) being the corresponding Lagrange multipliers. Similarly, (cid:126) Φ h ∈ V h ( γ h ( t ))is a Lagrange multiplier for (4.13c). It turns out that when mimicking the continuous approach fromSection 3 on the discrete level, we need several technical definitions to make the arguments rigorous.We present the majority of the necessary definitions and properties next, before proceeding with takingsuitable variations of (4.14).Following [21, (5.23)], we define the discrete material velocity for (cid:126)z ∈ Γ h ( t ) by (cid:126) V h ( (cid:126)z, t ) := K (cid:88) k =1 (cid:20) dd t (cid:126)q hk ( t ) (cid:21) χ hk ( (cid:126)z, t ) . We also introduce the finite element spaces W hT (Γ hi,T ) := { φ ∈ C (Γ hi,T ) : φ ( · , t ) ∈ W h (Γ hi ( t )) ∀ t ∈ [0 , T ] ,φ ( (cid:126)q hi,k ( · ) , · ) ∈ H (0 , T ) ∀ k ∈ { , . . . , K i }} , where Γ hi,T := (cid:83) t ∈ [0 ,T ] Γ hi ( t ) × { t } , as well as the vector valued analogue V hT (Γ hi,T ). In a similar fashion,we introduce W hT ( σ hj,T ) and V hT ( σ hj,T ) via e.g. W hT ( σ hj,T ) := { φ ∈ C ( σ hj,T ) : φ ( · , t ) is linear ∀ t ∈ [0 , T ] , φ ( (cid:126)q hj k ( · ) , · ) ∈ H (0 , T ) k = 1 , . . . , d } , where { (cid:126)q hj k ( t ) } dk =1 are the vertices of σ hj ( t ), and where σ hj,T := (cid:83) t ∈ [0 ,T ] σ hj ( t ) × { t } , for j ∈ { , . . . , J } .Moreover, we define the analogue variants W hT (Γ hT ) and V hT (Γ hT ) on Γ hT = (cid:83) t ∈ [0 ,T ] Γ h ( t ) × { t } , as wellas W hT ( γ hT ) and V hT ( γ hT ) on γ hT := (cid:83) t ∈ [0 ,T ] γ h ( t ) × { t } , with the scalar space for the latter e.g. beinggiven by W hT ( γ hT ) := { ψ ∈ C ( γ hT ) : ∃ χ ∈ W hT (Γ hT ) χ ( · , t ) | γ h ( t ) = ψ ( · , t ) ∀ t ∈ [0 , T ] } . Then, similarly to (3.1), we define the discrete material derivatives on Γ h ( t ) element-by-element viathe equations ( ∂ ◦ ,ht φ ) | σ hj ( t ) = ( φ t + (cid:126) V h . ∇ φ ) | σ hj ( t ) ∀ φ ∈ W hT ( σ hj,T ) , j ∈ { , . . . , J } . On differentiating (4.1) with respect to t , it immediately follows that ∂ ◦ ,ht χ hk = 0 ∀ k ∈ { , . . . , K } , (4.15)see also [21, Lem. 5.5]. It follows directly from (4.15) that ∂ ◦ ,ht φ ( · , t ) = K (cid:88) k =1 χ hk ( · , t ) dd t φ k ( t ) on Γ h ( t ) (4.16)for φ ( · , t ) = (cid:80) Kk =1 φ k ( t ) χ hk ( · , t ) ∈ W h (Γ h ( t )). 15 . W. Barrett, H. Garcke, R. N¨urnberg We recall from [21, Lem. 5.6] thatdd t (cid:90) σ hj ( t ) φ d H d − = (cid:90) σ hj ( t ) ∂ ◦ ,ht φ + φ ∇ s . (cid:126) V h d H d − ∀ φ ∈ W hT ( σ hj,T ) , j ∈ { , . . . , J } . (4.17)Similarly, we recall from [7, Lem. 3.1] thatdd t (cid:104) η, φ (cid:105) hσ hj ( t ) = (cid:104) ∂ ◦ ,ht η, φ (cid:105) hσ hj ( t ) + (cid:104) η, ∂ ◦ ,ht φ (cid:105) hσ hj ( t ) + (cid:104) η φ, ∇ s . (cid:126) V h (cid:105) hσ hj ( t ) ∀ η, φ ∈ W hT ( σ hj,T ) , (4.18)for all j ∈ { , . . . , J } . Moreover, it holds thatdd t (cid:104) η, φ (cid:105) hγ h ( t ) = (cid:104) ∂ ◦ ,ht η, φ (cid:105) hγ h ( t ) + (cid:104) η, ∂ ◦ ,ht φ (cid:105) hγ h ( t ) + (cid:104) η φ, (cid:126) id s . (cid:126) V hs (cid:105) hγ h ( t ) ∀ η, φ ∈ W hT ( γ hT ) . (4.19)We also note the discrete version of the time derivative variant of (3.17),dd t (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ hi ( t ) = (cid:68) ∇ s (cid:126)η, ∇ s (cid:126) V h (cid:69) Γ hi ( t ) + (cid:68) ∇ s . (cid:126)η, ∇ s . (cid:126) V h (cid:69) Γ hi ( t ) − (cid:68) ( ∇ s (cid:126)η ) T , D ( (cid:126) V h ) ( ∇ s (cid:126) id) T (cid:69) Γ hi ( t ) ∀ (cid:126)η ∈ { (cid:126)ξ ∈ V hT (Γ hi,T ) : ∂ ◦ ,ht (cid:126)ξ = (cid:126) } , (4.20)as well as the corresponding version for γ h ( t ),dd t (cid:68) (cid:126) id s , (cid:126)η s (cid:69) γ h ( t ) = (cid:68) P hγ (cid:126)η s , (cid:126) V hs (cid:69) γ h ( t ) ∀ (cid:126)η ∈ { (cid:126)ξ ∈ V hT ( γ hT ) : ∂ ◦ ,ht (cid:126)ξ = (cid:126) } , (4.21)which follows similarly to (3.20). Here, similarly to (3.21), we have defined P hγ = Id − (cid:126) id s ⊗ (cid:126) id s on γ h ( t ) . (4.22)Finally, when taking variations of (4.14), we need to compute variations of the discrete vertexnormals (cid:126)ω hi . To this end, for any given (cid:126)χ ∈ V h (Γ h ( t )) we introduce Γ hε ( t ) as in (3.5) and ∂ ,hε definedby (3.14), both with Γ( t ) replaced by Γ h ( t ). We then observe that it follows from (4.4) with w = 1and the discrete analogue of (3.15) that (cid:68) (cid:126)z, ∂ ,hε (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:68) (cid:126)z, ∂ ,hε (cid:126)ν hi (cid:69) h Γ hi ( t ) + (cid:68) ( (cid:126)z . ( (cid:126)ν hi − (cid:126)ω hi )) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ hi ( t ) ∀ (cid:126)z ∈ V h (Γ hi ( t )) , (cid:126)χ ∈ V h (Γ h ( t )) , (4.23)where ∂ ,hε (cid:126)z = (cid:126)
0. In addition, we note that for all (cid:126)ξ , (cid:126)η ∈ V h (Γ hi ( t )) with ∂ ,hε (cid:126)ξ = ∂ ,hε (cid:126)η = (cid:126) ∂ ,hε π hi (cid:20)(cid:18) (cid:126)ξ . (cid:126)ω hi | (cid:126)ω hi | (cid:19) (cid:18) (cid:126)η . (cid:126)ω hi | (cid:126)ω hi | (cid:19)(cid:21) = π hi (cid:104) (cid:126)G hi ( (cid:126)ξ, (cid:126)η ) . ∂ ,hε (cid:126)ω hi (cid:105) on Γ hi ( t ) , (4.24)where (cid:126)G hi ( (cid:126)ξ, (cid:126)η ) = 1 | (cid:126)ω hi | (cid:32) ( (cid:126)ξ . (cid:126)ω hi ) (cid:126)η + ( (cid:126)η . (cid:126)ω hi ) (cid:126)ξ − (cid:126)η . (cid:126)ω hi ) ( (cid:126)ξ . (cid:126)ω hi ) | (cid:126)ω hi | (cid:126)ω hi (cid:33) , (4.25)and where π hi ( t ) : C (Γ hi ( t )) → W h (Γ hi ( t )) is the standard interpolation operator at the nodes { (cid:126)q hi,k ( t ) } K i k =1 . It follows that (cid:126)G hi ( (cid:126)ξ, (cid:126)η ) . (cid:126)ω hi = 0 ∀ (cid:126)ξ, (cid:126)η ∈ V h (Γ hi ( t )) . (4.26)We are now in a position to formally derive the L –gradient flow of E h ( t ) subject to the sideconstraints (4.13a–c). In particular, on recalling the formal calculus of PDE constrained optimization,we set [ δδ Γ h L h ]( (cid:126)χ ) = − (cid:80) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)χ (cid:69) h Γ hi ( t ) for (cid:126)χ ∈ V h (Γ h ( t )), [ δδ(cid:126)κ hi L h ]( (cid:126)ξ ) = 0 for (cid:126)ξ ∈ V h (Γ hi ( t )),[ δδ(cid:126)Y hi L h ]( (cid:126)η ) = 0 for (cid:126)η ∈ V h (Γ hi ( t )), [ δδ (cid:126) m hi L h ]( (cid:126)ϕ ) = 0 for (cid:126)ϕ ∈ V h ( γ h ( t )), [ δδ(cid:126)κ hγ L h ]( (cid:126)φ ) = 0 for (cid:126)φ ∈ RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES V h ( γ h ( t )), leading to (cid:126)Z h = (cid:80) i =1 α Gi (cid:126) m hi , [ δδ (cid:126)Z h L h ]( (cid:126)φ ) = 0 for (cid:126)φ ∈ V h ( γ h ( t )) and [ δδ(cid:126) Φ h L h ]( (cid:126)η ) = 0 for (cid:126)η ∈ V h ( γ h ( t )). Here we recall the definition of θ h(cid:63) in (4.8). We employ this doctored version of θ h inorder to obtain existence and uniqueness for the fully discrete approximation introduced in the nextsection. See also Remark 3.1 in [6], where the analogue to our situation here corresponds to two curvesmeeting in the plane, i.e. N = d = 2 in their notation.Overall this gives rise to the following semidiscrete finite element approximation of the gradient flow(3.29), where we have noted the discrete version of (3.16), (4.23), (4.24), (4.26), variational versionsof (4.18)–(4.21) and that ∂ ,hε θ h = 0. Given Γ h (0), find (Γ h ( t )) t ∈ (0 ,T ] such that (cid:126) id | Γ h ( · ) ∈ V hT (Γ hT ). Inaddition, for all t ∈ (0 , T ] find ( (cid:126)κ hi , (cid:126)Y hi ) ∈ [ V h (Γ hi ( t ))] , i = 1 , (cid:126)κ hγ ∈ V h ( γ h ( t )), (cid:126) m hi ∈ V h ( γ h ( t )), i = 1 ,
2, and C Φ h ∈ V h ( γ h ( t )) such that (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)χ (cid:69) h Γ hi ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126)χ (cid:69) hγ h ( t ) = (cid:88) i =1 (cid:20)(cid:68) ∇ s (cid:126)Y hi , ∇ s (cid:126)χ (cid:69) Γ hi ( t ) + (cid:68) ∇ s . (cid:126)Y hi , ∇ s . (cid:126)χ (cid:69) Γ hi ( t ) − (cid:68) ( ∇ s (cid:126)Y hi ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ hi ( t ) − (cid:68) [ α i | (cid:126)κ hi − κ i (cid:126)ν hi | − (cid:126)Y hi . Q hi,θ h (cid:126)κ hi )] ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ hi ( t ) − α i κ i (cid:68) (cid:126)κ hi , [ ∇ s (cid:126)χ ] T (cid:126)ν hi (cid:69) h Γ hi ( t ) + (cid:68) (1 − θ h ) ( (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) . (cid:126)ν hi ) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , [ ∇ s (cid:126)χ ] T (cid:126)ν hi (cid:69) h Γ hi ( t ) (cid:21) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126)κ hγ . (cid:126) m hi , (cid:126) id s . (cid:126)χ s (cid:69) hγ h ( t ) + (cid:68) P hγ ( (cid:126) m hi ) s , (cid:126)χ s (cid:69) γ h ( t ) (cid:21) − ς (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ h ( t ) ∀ (cid:126)χ ∈ V h (Γ h ( t )) , (4.27a) (cid:68) Q hi,θ h (cid:126)κ hi , (cid:126)η (cid:69) h Γ hi ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ hi ( t ) = (cid:68) (cid:126) m hi , (cid:126)η (cid:69) hγ h ( t ) ∀ (cid:126)η ∈ V h (Γ hi ( t )) , i = 1 , , (4.27b) (cid:68) (cid:126)κ hγ , (cid:126)χ (cid:69) hγ h ( t ) + (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ h ( t ) = 0 ∀ (cid:126)χ ∈ V h ( γ h ( t )) , (4.27c) C ( (cid:126) m h + (cid:126) m h ) = (cid:126) γ h ( t ) , (4.27d) α Gi (cid:126)κ hγ + (cid:126)Y hi + C (cid:126) Φ h = (cid:126) γ h ( t ) , i = 1 , , (4.27e) (cid:68) α i ( (cid:126)κ hi − κ i (cid:126)ν hi ) − Q hi,θ h (cid:126)Y hi , (cid:126)ξ (cid:69) h Γ hi ( t ) = 0 ∀ (cid:126)ξ ∈ V h (Γ hi ( t )) , i = 1 , . (4.27f)We observe that choosing (cid:126)ξ = α − i (cid:126)π hi [ Q hi,θ h (cid:126)η ] in (4.27f) and combining with (4.27b), on recalling (4.10)and (4.4), yields that α − i (cid:68) Q hi,θ h (cid:126)Y hi , Q hi,θ h (cid:126)η (cid:69) h Γ hi ( t ) + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ hi ( t ) = (cid:68) (cid:126) m hi , (cid:126)η (cid:69) hγ h ( t ) − κ i (cid:68) (cid:126)ω hi , (cid:126)η (cid:69) h Γ hi ( t ) ∀ (cid:126)η ∈ V h (Γ hi ( t )) . (4.28)Here (cid:126)π hi ( t ) : [ C (Γ hi ( t ))] d → V h (Γ hi ( t )) is the standard interpolation operator at the nodes { (cid:126)q hi,k ( t ) } K i k =1 .17 . W. Barrett, H. Garcke, R. N¨urnberg In order to be able to consider area and volume preserving variants of (4.27a–f), we introduce theLagrange multipliers λ A,hi ( t ) ∈ R , i = 1 ,
2, and λ V,h ( t ) ∈ R for the constraintsdd t H d − (Γ hi ( t )) = (cid:68) ∇ s . (cid:126) V h , (cid:69) Γ hi ( t ) = (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) Γ hi ( t ) = 0 , (4.29)where we recall (4.17), anddd t L d (Ω h ( t )) = (cid:68) (cid:126) V h , (cid:126)ν h (cid:69) Γ h ( t ) = (cid:68) (cid:126) V h , (cid:126)ω h (cid:69) h Γ h ( t ) = 0 , (4.30)where we note a discrete variant of (3.4) and (4.6). Here Ω h ( t ) denotes the interior of Γ h ( t ). Onrecalling (4.7), (4.10) and (4.11), we can rewrite the constraint (4.30) as0 = (cid:68) (cid:126) V h , (cid:126)ω h (cid:69) h Γ h ( t ) = (cid:88) i =1 (cid:68) (cid:126) V h , (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)ω h (cid:69) h Γ hi ( t ) . (4.31)Hence, on writing (4.27a) as (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)χ (cid:69) h Γ hi ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126)χ (cid:69) hγ h ( t ) = (cid:68) (cid:126)r h , (cid:126)χ (cid:69) h Γ h ( t ) , we consider (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)χ (cid:69) h Γ hi ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126)χ (cid:69) hγ h ( t ) = (cid:68) (cid:126)r h , (cid:126)χ (cid:69) h Γ h ( t ) − λ V,h (cid:68) (cid:126)ω h , (cid:126)χ (cid:69) h Γ h ( t ) − (cid:88) i =1 λ A,hi (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ hi ( t ) (4.32)for all (cid:126)χ ∈ V h (Γ h ( t )), where λ V,h ( t ) ∈ R and λ A,hi ( t ) ∈ R , i = 1 ,
2, need to be determined. Of course,if we consider a volume preserving variant only, then we let λ A,h ( t ) = λ A,h ( t ) = 0 and λ V,h ( t ) = (cid:20)(cid:68) (cid:126)r h , (cid:126)ω h (cid:69) h Γ h ( t ) − (cid:37) (cid:68) (cid:126) V h , (cid:126)ω h (cid:69) hγ h ( t ) (cid:21) / (cid:68) (cid:126)ω h , (cid:126)ω h (cid:69) h Γ h ( t ) , (4.33)which we derived on choosing (cid:126)χ = (cid:126)ω h in (4.32), and noting (4.31).For the general volume and area preserving flow, we introduce the projection (cid:126) Π h : V h (Γ h ( t )) → V h (Γ h ( t )) onto V h (Γ h ( t )), recall (4.2b), and similarly (cid:126) Π hi, : V h (Γ hi ( t )) → V h (Γ hi ( t )). We introduce thesymmetric bilinear forms a hi,θ : V h (Γ hi ( t )) × V h (Γ hi ( t )) → R by setting a hi,θ ( (cid:126)ζ, (cid:126)η ) = (cid:68) Q hi,θ h (cid:126)ζ, (cid:126) Π hi, (cid:126)η (cid:69) h Γ hi ( t ) ∀ (cid:126)ζ, (cid:126)η ∈ V h (Γ hi ( t )) , i = 1 , , (4.34)where we have noted (4.10). It holds that a hi,θ ( (cid:126)ζ, (cid:126)ζ ) ≥ (cid:126)ζ ∈ V h (Γ hi ( t )), with the inequality beingstrict if (cid:126) Π hi, [ (cid:126)Q hi,θ h (cid:126)ζ ] (cid:54) = (cid:126)
0. Hence the Cauchy–Schwarz inequality holds, i.e. | a hi,θ ( (cid:126)ζ, (cid:126)η ) | ≤ [ a hi,θ ( (cid:126)ζ, (cid:126)ζ )] [ a hi,θ ( (cid:126)η, (cid:126)η )] ∀ (cid:126)ζ, (cid:126)η ∈ V h (Γ hi ( t )) , i = 1 , , (4.35)with strict inequality if (cid:126) Π hi, [ (cid:126)Q hi,θ h (cid:126)ζ ] and (cid:126) Π hi, [ (cid:126)Q hi,θ h (cid:126)η ] are linearly independent. Then we note, onrecalling (4.7), (4.27b) and (4.10), that − (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) Π h (cid:126)ω h (cid:69) Γ hi ( t ) = − (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) Π hi, (cid:126)ω hi (cid:69) Γ hi ( t ) = a hi,θ ( (cid:126)κ hi , (cid:126)ω hi ) = (cid:68) (cid:126)ω hi , (cid:126) Π hi, (cid:126)κ hi (cid:69) h Γ hi ( t ) (4.36a)and − (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) Π hi, (cid:126)κ hi (cid:69) Γ hi ( t ) = a hi,θ ( (cid:126)κ hi , (cid:126)κ hi ) . (4.36b)18 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
In addition, it follows from (4.7), (4.10) and (4.34) that (cid:68) (cid:126)ω h , (cid:126) Π h (cid:126)ω h (cid:69) h Γ h ( t ) = (cid:88) i =1 (cid:68) (cid:126)ω hi , (cid:126) Π h (cid:126)ω h (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:68) (cid:126)ω hi , (cid:126) Π hi, (cid:126)ω hi (cid:69) h Γ hi ( t ) = (cid:88) i =1 a hi,θ ( (cid:126)ω hi , (cid:126)ω hi ) . (4.37)Then (4.32), (4.37) and (4.36a,b) yield that ( λ V,h , λ
A,h , λ A,h )( t ) are such that (cid:80) i =1 a hi,θ ( (cid:126)ω hi , (cid:126)ω hi ) a h ,θ ( (cid:126)κ h , (cid:126)ω h ) a h ,θ ( (cid:126)κ h , (cid:126)ω h ) a h ,θ ( (cid:126)κ h , (cid:126)ω h ) a h ,θ ( (cid:126)κ h , (cid:126)κ h ) 0 a h ,θ ( (cid:126)κ h , (cid:126)ω h ) 0 a h ,θ ( (cid:126)κ h , (cid:126)κ h ) − λ V,h ( t ) λ A,h ( t ) λ A,h ( t ) = b ( t ) b ( t ) b ( t ) , (4.38a)where b ( t ) = (cid:88) i =1 (cid:68) (cid:126) Π h (cid:126) V h − (cid:126) V h , (cid:126)ω h (cid:69) h Γ hi ( t ) − (cid:68) (cid:126)r h , (cid:126) Π h (cid:126)ω h (cid:69) h Γ h ( t ) , (4.38b) b i ( t ) = (cid:68) (cid:126) Π hi, (cid:126) V h − (cid:126) V h , Q hi,θ h (cid:126)κ hi (cid:69) h Γ hi ( t ) + (cid:68) (cid:126) m hi , (cid:126) V h (cid:69) hγ h ( t ) − (cid:68) (cid:126)r h , (cid:126) Π hi, (cid:126)κ hi (cid:69) h Γ h ( t ) , i = 1 , . (4.38c)On recalling (4.35), we observe that the matrix in (4.38a) is symmetric and positive definite as longas (cid:126) Π hi, (cid:126)ω hi and (cid:126) Π hi, [ (cid:126)Q hi,θ h (cid:126)κ hi ] are linearly independent, for i = 1 ,
2. The right hand sides (4.38b,c) areobtained by recalling (4.32), and on noting that (4.10) and (4.31) imply that (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) Π h (cid:126)ω h (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:20)(cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) Π h (cid:126)ω h − (cid:126)ω h (cid:69) h Γ hi ( t ) + (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126)ω h (cid:69) h Γ hi ( t ) (cid:21) = (cid:88) i =1 (cid:68) (cid:126) V h , (cid:126) Π h (cid:126)ω h − (cid:126)ω h (cid:69) h Γ hi ( t ) = (cid:88) i =1 (cid:68) (cid:126) Π h (cid:126) V h − (cid:126) V h , (cid:126)ω h (cid:69) h Γ hi ( t ) , (4.39)while (4.8), (4.10), (4.27b) and (4.29) yield that (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) Π hi, (cid:126)κ hi (cid:69) h Γ hi ( t ) = (cid:68) (cid:126) Π hi, (cid:126) V h , Q hi,θ h (cid:126)κ hi (cid:69) h Γ hi ( t ) = (cid:68) (cid:126) Π hi, (cid:126) V h − (cid:126) V h , Q hi,θ h (cid:126)κ hi (cid:69) h Γ hi ( t ) + (cid:68) (cid:126) m hi , (cid:126) V h (cid:69) hγ h ( t ) − (cid:68) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) Γ hi ( t ) = (cid:68) (cid:126) Π hi, (cid:126) V h − (cid:126) V h , Q hi,θ h (cid:126)κ hi (cid:69) h Γ hi ( t ) + (cid:68) (cid:126) m hi , (cid:126) V h (cid:69) hγ h ( t ) . (4.40)We see that on removing the last two rows and columns in (4.38a), we obtain an expression similar to(4.33) for λ V,h ( t ), but here we test with (cid:126) Π h (cid:126)ω h as opposed to (cid:126)ω h . Analogously, if we want to considerphase area preservations only, then removing the first row and column in (4.38a) yields a reducedsystem for the two Lagrange multipliers λ A,hi ( t ), i = 1 , L –gradient flow of E h ( t ) subject to the side constraints (4.13a–c). We will also show that for θ = 0the scheme produces conformal polyhedral surfaces Γ ( t ) and Γ ( t ). Here we recall from [10], see also[3, § hi ( t ), i = 1 ,
2, are conformal polyhedral surfaces if (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ hi ( t ) = 0 ∀ (cid:126)η ∈ (cid:110) (cid:126)ξ ∈ V h (Γ hi ( t )) : (cid:126)ξ ( (cid:126)q hi,k ( t )) . (cid:126)ω hi ( (cid:126)q hi,k ( t ) , t ) = 0 , k = 1 , . . . , K i (cid:111) , i = 1 , . (4.41)We recall from [3, 10] that conformal polyhedral surfaces exhibit good meshes. Moreover, we recallthat in the case d = 2, conformal polyhedral surfaces are equidistributed polygonal curves, see [2, 5].19 . W. Barrett, H. Garcke, R. N¨urnberg Theorem 4.1.
Let θ ∈ [0 , , (cid:37) ≥ and let { (Γ h , (cid:126)κ h , (cid:126)κ h , (cid:126)Y h , (cid:126)Y h , (cid:126)κ hγ , (cid:126) m h , (cid:126) m h , (cid:126) Φ h )( t ) } t ∈ [0 ,T ] be a solutionto (4.27a–f) . In addition, we assume that (cid:126)κ hγ ∈ V hT ( γ hT ) , (cid:126)κ hi , (cid:126)π hi [ Q hi,θ h (cid:126)κ hi ] ∈ V hT (Γ hi,T ) , (cid:126) m hi ∈ V hT ( γ hT ) , i = 1 , . Then dd t E h ((Γ hi ( t )) i =1 ) = − (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) . (4.42) Moreover, if θ = 0 then Γ h ( t ) and Γ h ( t ) are open conformal polyhedral surfaces for all t ∈ (0 , T ] . Proof.
Taking the time derivative of (4.13a), where we choose discrete test functions (cid:126)η such that ∂ ◦ ,ht (cid:126)η = (cid:126)
0, yields for i = 1 , (cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) , (cid:126)η (cid:69) h Γ hi ( t ) + (cid:68) [( Q hi,θ h (cid:126)κ hi ) . (cid:126)η ] ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) + (cid:68) ∇ s (cid:126) V h , ∇ s (cid:126)η (cid:69) Γ hi ( t ) + (cid:68) ∇ s . (cid:126) V h , ∇ s . (cid:126)η (cid:69) Γ hi ( t ) − (cid:68) ( ∇ s (cid:126)η ) T , D ( (cid:126) V h ) ( ∇ s (cid:126) id) T (cid:69) Γ hi ( t ) = (cid:68) ∂ ◦ ,ht (cid:126) m hi , (cid:126)η (cid:69) hγ h ( t ) + (cid:68) (cid:126) m hi . (cid:126)η, (cid:126) id s . (cid:126) V hs (cid:69) hγ h ( t ) , (4.43)where we have noted (4.18), (4.19), (4.20) and that (cid:126)π hi [ Q hi,θ h (cid:126)κ hi ] ∈ V hT (Γ hi,T ), (cid:126) m hi ∈ V hT ( γ hT ), i = 1 , ∂ ◦ ,ht (cid:126)χ = (cid:126) (cid:126)κ hγ ∈ V hT ( γ hT ), that (cid:68) ∂ ◦ ,ht (cid:126)κ hγ , (cid:126)χ (cid:69) hγ h ( t ) + (cid:68) (cid:126)κ hγ . (cid:126)χ, (cid:126) id s . (cid:126) V hs (cid:69) hγ h ( t ) + (cid:68) P hγ (cid:126)χ s , (cid:126) V hs (cid:69) γ h ( t ) = 0 . (4.44)Choosing (cid:126)χ = (cid:126) V h in (4.27a), (cid:126)η = (cid:126)Y hi in (4.43), i = 1 ,
2, and combining yields, on noting the discretevariant of (3.16), that (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ h ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) + (cid:88) i =1 (cid:20) (cid:68) [ α i | (cid:126)κ hi − κ i (cid:126)ν hi | − (cid:126)Y hi . Q hi,θ h (cid:126)κ hi ] ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − α i κ i (cid:68) (cid:126)κ hi , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) + (cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) , (cid:126)Y hi (cid:69) h Γ hi ( t ) + (cid:68) ( Q hi,θ h (cid:126)κ hi . (cid:126)Y hi ) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) ( (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) . (cid:126)ν hi ) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) (cid:21) + ς (cid:68) (cid:126) id s , (cid:126) V hs (cid:69) γ h ( t ) − (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126)κ hγ . (cid:126) m hi , (cid:126) id s . (cid:126) V hs (cid:69) hγ h ( t ) + (cid:68) P hγ ( (cid:126) m hi ) s , (cid:126) V hs (cid:69) γ h ( t ) (cid:21) = (cid:88) i =1 (cid:20)(cid:68) ∂ ◦ ,ht (cid:126) m hi , (cid:126)Y hi (cid:69) γ h ( t ) + (cid:68) (cid:126) m hi . (cid:126)Y hi , (cid:126) id s , (cid:126) V hs (cid:69) γ h ( t ) (cid:21) . (4.45)Choosing (cid:126)χ = (cid:80) i =1 α Gi (cid:126) m hi in (4.44) and recalling (4.27d,e) and (4.19), it follows from (4.45) that (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ h ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) + (cid:88) i =1 (cid:20) α i (cid:68) | (cid:126)κ hi − κ i (cid:126)ν hi | ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − α i κ i (cid:68) (cid:126)κ hi , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) + (cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) , (cid:126)Y hi (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) ( (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) . (cid:126)ν hi ) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) (cid:21) RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES + ς (cid:68) (cid:126) id s , (cid:126) V hs (cid:69) ∂ Γ h ( t ) = − (cid:88) i =1 α Gi (cid:20)(cid:68) ∂ ◦ ,ht (cid:126) m hi , (cid:126)κ hγ (cid:69) hγ h ( t ) + (cid:68) (cid:126) m hi , ∂ ◦ ,ht (cid:126)κ hγ (cid:69) hγ h ( t ) + (cid:68) (cid:126)κ hγ . (cid:126) m hi , (cid:126) id s . (cid:126) V hs (cid:69) hγ h ( t ) (cid:21) = − dd t (cid:88) i =1 α Gi (cid:68) (cid:126)κ hγ , (cid:126) m hi (cid:69) hγ h ( t ) . (4.46)We have from (4.10), (4.27f) and (4.4) that (cid:88) i =1 (cid:20)(cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) , (cid:126)Y hi (cid:69) h Γ hi ( t ) − α i κ i (cid:68) (cid:126)κ hi , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) (cid:21) = (cid:88) i =1 (cid:20)(cid:68) ∂ ◦ ,ht (cid:126)κ hi , Q hi,θ h (cid:126)Y hi (cid:69) h Γ hi ( t ) − α i κ i (cid:68) (cid:126)κ hi − κ i (cid:126)ν hi , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) + (cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) − Q hi,θ h ∂ ◦ ,ht (cid:126)κ hi , (cid:126)Y hi (cid:69) h Γ hi ( t ) (cid:21) = (cid:88) i =1 (cid:20) α i (cid:68) ∂ ◦ ,ht | (cid:126)κ hi − κ i (cid:126)ν hi | , (cid:69) h Γ hi ( t ) + (cid:68) ∂ ◦ ,ht ( Q hi,θ h (cid:126)κ hi ) − Q hi,θ h ∂ ◦ ,ht (cid:126)κ hi , (cid:126)Y hi (cid:69) h Γ hi ( t ) (cid:21) . (4.47)Combining (4.46) and (4.47), on noting (4.18), (4.19), (4.12), ∂ ◦ ,ht θ h = 0 (which follows from (4.16)and (4.8)), (cid:126)κ hi ∈ V hT (Γ hi,T ), i = 1 , (cid:126)ν hi | σ hi,j ( · ) ∈ V hT ( σ hi,j,T ), j = 1 , . . . , J i , i = 1 ,
2, (which follows from thediscrete analogue of (3.16) and as (cid:126) id | Γ h ( · ) ∈ V hT (Γ hT )) and the invariance of m (Γ hi ( t )) under continuousdeformations, yields that (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ hi ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) + dd t E h ((Γ hi ( t )) i =1 ) + (cid:88) i =1 P i = 0 , where, on noting (4.9), P i := (cid:42) (1 − θ h ) (cid:126)κ hi . ∂ ◦ ,ht (cid:126)ω hi , (cid:126)Y hi . (cid:126)ω hi | (cid:126)ω hi | (cid:43) h Γ hi ( t ) + (cid:28) (1 − θ h ) (cid:126)Y hi . ∂ ◦ ,ht (cid:126)ω hi , (cid:126)κ hi . (cid:126)ω hi | (cid:126)ω hi | (cid:29) h Γ hi ( t ) − (cid:42) (1 − θ h ) ( (cid:126)κ hi . (cid:126)ω hi ) ( (cid:126)Y hi . (cid:126)ω hi ) , (cid:126)ω hi . ∂ ◦ ,ht (cid:126)ω hi | (cid:126)ω hi | (cid:43) h Γ hi ( t ) − (cid:68) (1 − θ h ) ( (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) . (cid:126)ν hi ) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) , i = 1 , . (4.48)It remains to show that P and P as defined in (4.48) vanish. To see this, we observe that it followsfrom (4.26), (4.25) and the time derivative version of (4.23) that P i = (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , ∂ ◦ ,ht (cid:126)ω hi (cid:69) h Γ hi ( t ) + (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) . ( (cid:126)ω hi − (cid:126)ν hi ) ∇ s (cid:126) id , ∇ s (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:68) (1 − θ h ) (cid:126)G hi ( (cid:126)Y hi , (cid:126)κ hi ) , ∂ ◦ ,ht (cid:126)ν hi (cid:69) h Γ hi ( t ) = 0 , i = 1 , . This proves the desired result (4.42).Finally, if θ = 0 then it immediately follows from (4.27b) that (4.41) holds. Hence Γ h ( t ) and Γ h ( t ) areopen conformal polyhedral surfaces. Theorem 4.2.
Let θ ∈ [0 , , (cid:37) ≥ and let { (Γ h , (cid:126)κ h , (cid:126)κ h , (cid:126)Y h , (cid:126)Y h , (cid:126)κ hγ , (cid:126) m h , (cid:126) m h , (cid:126) Φ h , λ V,h , λ
A,h , λ A,h )( t ) } t ∈ [0 ,T ] be a solution to (4.32 ), (4.27b–f) and (4.38a) . In addition, we assume that (cid:126)κ hγ ∈ V hT ( γ hT ) , . W. Barrett, H. Garcke, R. N¨urnberg (cid:126)κ hi , (cid:126)π hi [ Q hi,θ h (cid:126)κ hi ] ∈ V hT (Γ hi,T ) , (cid:126) m hi ∈ V hT ( γ hT ) , i = 1 , . Then it holds that dd t E h ((Γ hi ( t )) i =1 ) = − (cid:88) i =1 (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ hi ( t ) − (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) , (4.49) as well as dd t H d − (Γ hi ( t )) = 0 , i = 1 , , dd t L d (Ω h ( t )) = 0 , (4.50) where Ω h ( t ) denotes the region bounded by Γ h ( t ) . Moreover, if θ = 0 then Γ h ( t ) and Γ h ( t ) are openconformal polyhedral surfaces for all t ∈ (0 , T ] . Proof.
We recall that on choosing ( λ V,h , λ
A,h , λ A,h ) solving the system (4.38a) yields that (4.29) and(4.30) hold, and hence the desired results (4.50) hold. The stability result (4.49) directly follows fromthe proof of Theorem 4.1. In particular, choosing (cid:126)χ = (cid:126) V h in (4.32), on noting (4.29) and (4.30), yieldsthat (cid:68) Q hi,θ h(cid:63) (cid:126) V h , (cid:126) V h (cid:69) h Γ hi ( t ) + (cid:37) (cid:68) (cid:126) V h , (cid:126) V h (cid:69) hγ h ( t ) = (cid:68) (cid:126)r h , (cid:126) V h (cid:69) h Γ h ( t ) . Combining this with (4.43) yields that (4.45) holds, and the rest of the proof proceeds as that ofTheorem 4.1. Finally, as in the proof of Theorem 4.1, for θ = 0 it follows from (4.27b) that Γ h ( t ) andΓ h ( t ) are conformal polyhedral surfaces. Fully discrete finite element approximation
In this section we consider a fully discrete variant of the scheme (4.27a–f) from Section 4. To thisend, let 0 = t < t < . . . < t M − < t M = T be a partitioning of [0 , T ] into possibly variable timesteps ∆ t m := t m +1 − t m , m = 0 , . . . , M −
1. Let Γ m be a ( d − h ( t m ), m = 0 , . . . , M , with the two parts Γ mi , i = 1 , γ m . Following [19], we now parameterize the new surface Γ m +1 over Γ m . Hence, we introduce thefollowing finite element spaces. Let Γ m = (cid:83) Jj =1 σ mj , where { σ mj } Jj =1 is a family of mutually disjointopen triangles with vertices { (cid:126)q mk } Kk =1 . Then for m = 0 , . . . , M −
1, let V h (Γ m ) := { (cid:126)χ ∈ [ C (Γ m )] d : (cid:126)χ | σ mj is linear ∀ j = 1 , . . . , J } =: [ W h (Γ m )] d . (5.1)We denote the standard basis of W h (Γ m ) by { χ mk } Kk =1 . In addition, similarly to the semidiscrete settingin Section 4, we introduce the spaces W h (Γ mi ) and V h (Γ mi ), denoting the standard basis of W h (Γ mi )by { χ mi,k } K i k =1 , as well as V h ( γ m ), and the interpolation operators π mi : C (Γ mi ) → W h (Γ mi ) and similarly (cid:126)π mi : [ C (Γ mi )] d → V h (Γ mi ).We also introduce the L –inner products (cid:104)· , ·(cid:105) Γ m , (cid:104)· , ·(cid:105) Γ mi and (cid:104)· , ·(cid:105) γ m , as well as their mass lumpedinner variants (cid:104)· , ·(cid:105) h Γ m , (cid:104)· , ·(cid:105) h Γ mi and (cid:104)· , ·(cid:105) hγ m . Similarly to (4.3) and (4.5) we introduce the discrete vertexnormals (cid:126)ω mi := (cid:80) K i k =1 χ mi,k (cid:126)ω mi,k ∈ V h (Γ mi ) and (cid:126)ω m := (cid:80) Kk =1 χ mk (cid:126)ω mk ∈ V h (Γ m ).We make the following mild assumption.( A ) We assume for m = 0 , . . . , M − H d − ( σ mj ) > j = 1 , . . . , J , and that (cid:126) (cid:54)∈ { (cid:126)ω mi,k : k =1 , . . . , K i , i = 1 , } . Moreover, in the case C = 1 and θ = 0 we assume that dim span { (cid:126)ω mi,k : k = 1 , . . . , K i , i = 1 , } = d , for m = 0 , . . . , M − RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
In addition, and similarly to (4.8) and (4.9), we introduce θ m and θ m(cid:63) ∈ W h (Γ m ), and then Q mi,θ m , Q mi,θ m(cid:63) ∈ [ W h (Γ mi )] d × d , by setting Q mi,θ m ( (cid:126)q mi,k ) = θ m ( (cid:126)q mi,k ) (cid:126) Id + (1 − θ m ( (cid:126)q mi,k )) | (cid:126)ω mi,k | − (cid:126)ω mi,k ⊗ (cid:126)ω mi,k and Q mi,θ m(cid:63) ( (cid:126)q mi,k ) = θ m(cid:63) ( (cid:126)q mi,k ) (cid:126) Id + (1 − θ m(cid:63) ( (cid:126)q mi,k )) | (cid:126)ω mi,k | − (cid:126)ω mi,k ⊗ (cid:126)ω mi,k for k = 1 , . . . , K i , i = 1 ,
2. Similarly to(4.25) and (4.22), we let (cid:126)G mi ( (cid:126)ξ, (cid:126)η ) = 1 | (cid:126)ω mi | (cid:32) ( (cid:126)ξ . (cid:126)ω mi ) (cid:126)η + ( (cid:126)η . (cid:126)ω mi ) (cid:126)ξ − (cid:126)η . (cid:126)ω mi ) ( (cid:126)ξ . (cid:126)ω mi ) | (cid:126)ω mi | (cid:126)ω mi (cid:33) and P mγ = Id − (cid:126) id s ⊗ (cid:126) id s on γ m . On recalling (4.28), we consider the following fully discrete approximation of (4.27a–f). For m =0 , . . . , M −
1, find (cid:126)X m +1 ∈ V h (Γ m ), ( (cid:126)Y m +1 i , (cid:126) m m +1 i ) i =1 ∈ V h (Γ m ) × V h ( γ m ) × V h (Γ m ) × V h ( γ m ), (cid:126)κ m +1 γ ∈ V h ( γ m ) and C Φ m +1 ∈ V h ( γ m ) such that (cid:88) i =1 (cid:42) Q mi,θ m(cid:63) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 i , ∇ s (cid:126)χ (cid:69) Γ mi + α Gi (cid:10) ( (cid:126) m m +1 i ) s , (cid:126)χ s (cid:11) γ m + ς (cid:68) (cid:126)X m +1 s , (cid:126)χ s (cid:69) γ m + (cid:37) (cid:42) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) hγ m = (cid:88) i =1 (cid:20)(cid:68) ∇ s . (cid:126)Y mi , ∇ s . (cid:126)χ (cid:69) Γ mi − (cid:68) ( ∇ s (cid:126)Y mi ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ mi − (cid:68) [ α i | (cid:126)κ mi − κ i (cid:126)ν mi | − (cid:126)Y mi . Q mi,θ m (cid:126)κ mi )] ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ mi − α i κ i (cid:10) (cid:126)κ mi , [ ∇ s (cid:126)χ ] T (cid:126)ν mi (cid:11) h Γ mi + (cid:68) (1 − θ m ) ( (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) . (cid:126)ν mi ) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ mi − (cid:68) (1 − θ m ) (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) , [ ∇ s (cid:126)χ ] T (cid:126)ν mi (cid:69) h Γ mi (cid:21) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126)κ mγ . (cid:126) m mi , (cid:126) id s . (cid:126)χ s (cid:69) hγ m + (cid:10) (Id + P mγ ) ( (cid:126) m mi ) s , (cid:126)χ s (cid:11) γ m (cid:21) − λ V,m (cid:104) (cid:126)ω m , (cid:126)χ (cid:105) h Γ m − (cid:88) i =1 λ A,mi (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ mi ∀ (cid:126)χ ∈ V h (Γ m ) , (5.2a) α − i (cid:68) Q mi,θ m (cid:126)Y m +1 i , Q mi,θ m (cid:126)η (cid:69) h Γ mi + (cid:68) ∇ s (cid:126)X m +1 , ∇ s (cid:126)η (cid:69) Γ mi = (cid:10) (cid:126) m m +1 i , (cid:126)η (cid:11) hγ m − κ i (cid:104) (cid:126)ω mi , (cid:126)η (cid:105) h Γ mi ∀ (cid:126)η ∈ V h (Γ mi ) , i = 1 , , (5.2b) (cid:10) (cid:126)κ m +1 γ , (cid:126)χ (cid:11) hγ m + (cid:68) (cid:126)X m +1 s , (cid:126)χ s (cid:69) γ m = 0 ∀ (cid:126)χ ∈ V h ( γ m ) , (5.2c) C ( (cid:126) m m +11 + (cid:126) m m +12 ) = (cid:126) γ m , (5.2d) α Gi (cid:126)κ m +1 γ + (cid:126)Y m +1 i + C (cid:126) Φ m +1 = (cid:126) γ m , i = 1 , , (5.2e)and set (cid:126)κ m +1 i = α − i (cid:126)π mi [ Q mi,θ m (cid:126)Y m +1 i ] + κ i (cid:126)ω mi and Γ m +1 i = (cid:126)X m +1 (Γ mi ), i = 1 ,
2. For m ≥ (cid:126)κ mi the function (cid:126)z ∈ V h (Γ mi ),defined by (cid:126)z ( (cid:126)q mi,k ) = (cid:126)κ mi ( (cid:126)q m − i,k ), k = 1 → K i , where (cid:126)κ mi ∈ V (Γ m − i ) is given, and similarly for e.g. (cid:126)Y mi , (cid:126) m mi and (cid:126)κ mγ . 23 . W. Barrett, H. Garcke, R. N¨urnberg We note that if C = α G = α G = 0 then the weak conormals (cid:126) m m +1 i play no role in the evolution.However, for surface area conservation they do play a role also in that case, see (5.3c) below. We alsoremark that the parameter (cid:37) ≥ γ m . In practice, this wasparticularly useful for simulations involving surface area preservation, and for C experiments withGaussian curvature energy contributions.Of course, (5.2a–e) with λ V,m = λ A,m = λ A,m = 0 corresponds to a fully discrete approximationof (4.27a–f), on recalling (4.28). For a fully discrete approximation of the volume and/or surface areapreserving flow, on recalling (4.38a–c), we let ( λ V,m , λ
A,m , λ A,m ) be the solution of (cid:80) i =1 a mi,θ ( (cid:126)ω mi , (cid:126)ω mi ) a m ,θ ( (cid:126)κ m , (cid:126)ω m ) a m ,θ ( (cid:126)κ m , (cid:126)ω m ) a m ,θ ( (cid:126)κ m , (cid:126)ω m ) a m ,θ ( (cid:126)κ m , (cid:126)κ m ) 0 a m ,θ ( (cid:126)κ m , (cid:126)ω m ) 0 a m ,θ ( (cid:126)κ m , (cid:126)κ m ) − λ V,m λ A,m λ A,m = b m b m b m , (5.3a)where, on noting the fully discrete variant of (4.7), b m = (cid:88) i =1 (cid:42) ( (cid:126) Π mi, − Id) (cid:126) id − (cid:126)X m − ∆ t m − , (cid:126)ω m (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y mi , ∇ s ( (cid:126) Π m (cid:126)ω m ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π m (cid:126)ω m (cid:69) h Γ m = (cid:88) i =1 (cid:42) ( (cid:126) Π mi, − Id) (cid:126) id − (cid:126)X m − ∆ t m − , (cid:126)ω mi (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y mi , ∇ s ( (cid:126) Π mi, (cid:126)ω mi ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π mi, (cid:126)ω mi (cid:69) h Γ mi , (5.3b) b mi = (cid:42) ( (cid:126) Π mi, − Id) (cid:126) id − (cid:126)X m − ∆ t m − , Q mi,θ m (cid:126)κ mi (cid:43) h Γ mi + (cid:42) (cid:126) m mi , (cid:126) id − (cid:126)X m − ∆ t m − (cid:43) hγ m − (cid:68) ∇ s (cid:126)Y mi , ∇ s ( (cid:126) Π mi, (cid:126)κ mi ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π mi, (cid:126)κ mi (cid:69) h Γ m = (cid:42) ( (cid:126) Π mi, − Id) (cid:126) id − (cid:126)X m − ∆ t m − , Q mi,θ m (cid:126)κ mi (cid:43) h Γ mi + (cid:42) (cid:126) m mi , (cid:126) id − (cid:126)X m − ∆ t m − (cid:43) hγ m − (cid:68) ∇ s (cid:126)Y mi , ∇ s ( (cid:126) Π mi, (cid:126)κ mi ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π mi, (cid:126)κ mi (cid:69) h Γ mi , i = 1 , . (5.3c)Here, for convenience, we have re-written (5.2a) as (cid:88) i =1 (cid:42) Q mi,θ m(cid:63) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 i , ∇ s (cid:126)χ (cid:69) Γ mi + α Gi (cid:10) ( (cid:126) m m +1 i ) s , (cid:126)χ s (cid:11) γ m + ς (cid:68) (cid:126)X m +1 s , (cid:126)χ s (cid:69) γ m + (cid:37) (cid:42) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) hγ m = (cid:68) (cid:126)f m , (cid:126)χ (cid:69) h Γ m − λ V,m (cid:104) (cid:126)ω m , (cid:126)χ (cid:105) h Γ m − (cid:88) i =1 λ A,mi (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ mi ∀ (cid:126)χ ∈ V h (Γ m ) , (5.4)and, analogously to (4.34), we have defined a mi,θ : V h (Γ mi ) × V h (Γ mi ) → R by setting a mi,θ ( (cid:126)ζ, (cid:126)η ) = (cid:68) Q mi,θ m (cid:126)ζ, (cid:126) Π mi, (cid:126)η (cid:69) h Γ mi ∀ (cid:126)ζ, (cid:126)η ∈ V h (Γ mi ) , i = 1 , . RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES
As before, we note that the matrix in (5.3a) is symmetric positive definite as long as (cid:126) Π mi, (cid:126)ω mi and (cid:126) Π mi, [ (cid:126)Q mi,θ m (cid:126)κ mi ] are linearly independent, for i = 1 , Theorem 5.1.
Let θ ∈ [0 , , (cid:37) ≥ and α , α > . Let the assumptions ( A ) hold. Then there existsa unique solution (cid:126)X m +1 ∈ V h (Γ m ) , ( (cid:126)Y m +1 i , (cid:126) m m +1 i ) i =1 ∈ V h (Γ m ) × V h ( γ m ) × V h (Γ m ) × V h ( γ m ) , (cid:126)κ m +1 γ ∈ V h ( γ m ) and C Φ m +1 ∈ V h ( γ m ) to (5.2a–e) . Proof.
As (5.2a–e) is linear, existence follows from uniqueness. To investigate the latter, we considerthe system: Find (cid:126)X ∈ V h (Γ m ), ( (cid:126)Y i , (cid:126) m i ) i =1 ∈ V h (Γ m ) × V h ( γ m ) × V h (Γ m ) × V h ( γ m ), (cid:126)κ γ ∈ V h ( γ m )and C Φ ∈ V h ( γ m ) such that (cid:88) i =1 (cid:20) t m (cid:68) Q mi,θ m(cid:63) (cid:126)X, (cid:126)χ (cid:69) h Γ mi − (cid:68) ∇ s (cid:126)Y i , ∇ s (cid:126)χ (cid:69) Γ mi + α Gi (cid:104) ( (cid:126) m i ) s , (cid:126)χ s (cid:105) γ m (cid:21) + ς (cid:68) (cid:126)X s , (cid:126)χ s (cid:69) γ m + (cid:37) (cid:68) (cid:126)X, (cid:126)χ (cid:69) hγ m = 0 ∀ (cid:126)χ ∈ V h (Γ m ) , (5.5a) α − i (cid:68) Q mi,θ m (cid:126)Y i , Q mi,θ m (cid:126)η (cid:69) h Γ mi + (cid:68) ∇ s (cid:126)X, ∇ s (cid:126)η (cid:69) Γ mi = (cid:104) (cid:126) m i , (cid:126)η (cid:105) hγ m ∀ (cid:126)η ∈ V h (Γ mi ) , i = 1 , , (5.5b) (cid:104) (cid:126)κ γ , (cid:126)χ (cid:105) hγ m + (cid:68) (cid:126)X s , (cid:126)χ s (cid:69) γ m = 0 ∀ (cid:126)χ ∈ V h ( γ m ) , (5.5c) C ( (cid:126) m + (cid:126) m ) = (cid:126) γ m , (5.5d) α Gi (cid:126)κ γ + (cid:126)Y i + C (cid:126) Φ = (cid:126) γ m , i = 1 , . (5.5e)Choosing (cid:126)χ = (cid:126)X in (5.5a), (cid:126)η = (cid:126)Y i in (5.5b), (cid:126)χ = (cid:126) m i in (5.5c) and noting (5.5d,e), leads to (cid:88) i =1 t m (cid:68) Q mi,θ m(cid:63) (cid:126)X, (cid:126)X (cid:69) h Γ mi + ς (cid:68) (cid:126)X s , (cid:126)X s (cid:69) γ m + (cid:37) (cid:68) (cid:126)X, (cid:126)X (cid:69) hγ m = (cid:88) i =1 (cid:20)(cid:68) (cid:126) m i , (cid:126)Y i (cid:69) hγ m − α − i (cid:68) Q mi,θ m (cid:126)Y i , Q mi,θ m (cid:126)Y i (cid:69) h Γ mi − α Gi (cid:68) ( (cid:126) m i ) s , (cid:126)X s (cid:69) γ m (cid:21) = (cid:88) i =1 (cid:20)(cid:68) (cid:126)Y i , (cid:126) m i (cid:69) hγ m + α Gi (cid:104) (cid:126) m i , (cid:126)κ γ (cid:105) γ m − α − i (cid:68) Q mi,θ m (cid:126)Y i , Q mi,θ m (cid:126)Y i (cid:69) h Γ mi (cid:21) = − (cid:88) i =1 α − i (cid:68) Q mi,θ m (cid:126)Y i , Q mi,θ m (cid:126)Y i (cid:69) h Γ mi . (5.6)It follows from (5.6) and the definition of Q mi,θ m(cid:63) that (cid:126)X = (cid:126) γ m , and so (5.5c) implies that (cid:126)κ γ = (cid:126) (cid:126)π mi [ Q mi,θ m (cid:126)Y i ] = (cid:126) i = 1 ,
2. Hence, on adding the two equations in(5.5b) with (cid:126)η = (cid:126)X , we obtain that (cid:68) ∇ s (cid:126)X, ∇ s (cid:126)X (cid:69) Γ m = 0, and so (cid:126)X = (cid:126)
0. Then (5.5b) implies that (cid:126) m = (cid:126) m = (cid:126)
0. Next, we have from (5.5e), on recalling that (cid:126)κ γ = (cid:126)
0, that there exists a (cid:126)Y ∈ V h (Γ m )such that (cid:126)Y i = (cid:126)Y | Γ mi , i = 1 ,
2. Choosing (cid:126)η = (cid:126)Y in (5.5a) yields that (cid:68) ∇ s (cid:126)Y , ∇ s (cid:126)Y (cid:69) Γ m = 0, and hence (cid:126)Y is constant. If C = 0 we immediately obtain from (5.5e) that (cid:126)Y = (cid:126)
0. If C = 1, on the other hand,it follows that Q mi,θ m ( (cid:126)q mi,k ) (cid:126)Y = (cid:126) k = 1 , . . . K i , i = 1 ,
2, and hence (cid:126)Y . (cid:126)ω mi ( (cid:126)q mi,k ) = 0 k = 1 , . . . K i , i = 1 , . (5.7)The definition of Q mi,θ m , recall the fully discrete version of (4.9), and (5.7) then yield that θ (cid:126)Y = (cid:126) θ ∈ (0 ,
1] we immediately obtain that (cid:126)Y = (cid:126)
0, while in the case θ = 0 it follows fromassumption ( A ) and (5.7) that (cid:126)Y = (cid:126)
0. Finally, we obtain that (cid:126)
Φ = (cid:126) . W. Barrett, H. Garcke, R. N¨urnberg
Implicit treatment of volume and area conservation
In practice it can be advantageous to consider implicit Lagrange multipliers ( λ V,m +1 , λ A,m +11 , λ
A,m +12 )in order to obtain better discrete volume and surface area conservations. In particular, we replace(5.4) with (cid:88) i =1 (cid:42) Q mi,θ m(cid:63) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 i , ∇ s (cid:126)χ (cid:69) Γ mi + α Gi (cid:10) ( (cid:126) m m +1 i ) s , (cid:126)χ s (cid:11) γ m + ς (cid:68) (cid:126)X m +1 s , (cid:126)χ s (cid:69) γ m + (cid:37) (cid:42) (cid:126)X m +1 − (cid:126) id∆ t m , (cid:126)χ (cid:43) hγ m = (cid:68) (cid:126)f m , (cid:126)χ (cid:69) h Γ m − λ V,m +1 (cid:104) (cid:126)ω m , (cid:126)χ (cid:105) h Γ m − (cid:88) i =1 λ A,m +1 i (cid:68) ∇ s (cid:126)X m +1 , ∇ s (cid:126)χ (cid:69) Γ mi ∀ (cid:126)χ ∈ V h (Γ m ) , (5.8)and require the coupled solutions (cid:126)X m +1 ∈ V h (Γ m ), ( (cid:126)Y m +1 i , (cid:126) m m +1 i ) i =1 ∈ V h (Γ m ) × V h ( γ m ) × V h (Γ m ) × V h ( γ m ), (cid:126)κ m +1 γ ∈ V h ( γ m ), C Φ m +1 ∈ V h ( γ m ) and ( λ V,m +1 , λ A,m +11 , λ
A,m +12 ) ∈ R to satisfy thenonlinear system (5.8), (5.2b–e) as well as an adapted variant of (5.3a–c), where the superscript m isreplaced by m + 1 in all occurrences of (cid:126) m mi , (cid:126)κ mi , (cid:126)Y mi , λ V,m and λ A,mi . In addition, (cid:126) id − (cid:126)X m − ∆ t m − in (5.3b,c)is replaced by (cid:126)X m +1 − (cid:126) id∆ t m . In practice this nonlinear system can be solved with a fixed point iteration asfollows. Let ( λ V,m +1 , , λ A,m +1 , , λ A,m +1 , ) = ( λ V,m , λ
A,m , λ A,m ) and (cid:126)X m +1 , = (cid:126) id | Γ m . Then, for j ≥ (cid:126)X m +1 ,j +1 , (cid:126)Y m +1 ,j +1 , (cid:126)κ m +1 ,j +1 ∂ Γ , (cid:126) m m +1 ,j +1 ) to the linear system (5.8), (5.2b–e), whereany superscript m + 1 on left hand sides is replaced by m + 1 , j + 1, and by m + 1 , j on the right handside of (5.8). Then let (cid:126)κ m +1 ,j +1 i = α − i (cid:126)π mi [ Q mi,θ m (cid:126)Y m +1 ,j +1 i ] + κ i (cid:126)ω mi be defined as usual, and compute( λ V,m +1 ,j +1 , λ A,m +1 ,j +11 , λ A,m +1 ,j +12 ) as the unique solution to (cid:80) i =1 a mi,θ ( (cid:126)ω mi , (cid:126)ω mi ) a m ,θ ( (cid:126)κ m +1 ,j +11 , (cid:126)ω m ) a m ,θ ( (cid:126)κ m +1 ,j +12 , (cid:126)ω m ) a m ,θ ( (cid:126)κ m +1 ,j +11 , (cid:126)ω m ) a m ,θ ( (cid:126)κ m +1 ,j +11 , (cid:126)κ m +1 ,j +11 ) 0 a m ,θ ( (cid:126)κ m +1 ,j +12 , (cid:126)ω m ) 0 a m ,θ ( (cid:126)κ m +1 ,j +12 , (cid:126)κ m +1 ,j +12 ) − λ V,m +1 ,j +1 λ A,m +1 ,j +11 λ A,m +1 ,j +12 = b m +1 ,j +10 b m +1 ,j +11 b m +1 ,j +12 , where b m +1 ,j +10 = (cid:88) i =1 (cid:42) ( (cid:126) Π mi, − Id) (cid:126)X m +1 ,j +1 − (cid:126) id∆ t m , (cid:126)ω mi (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 ,j +1 i , ∇ s ( (cid:126) Π mi, (cid:126)ω mi ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π mi, (cid:126)ω mi (cid:69) h Γ mi (cid:21) ,b m +1 ,j +1 i = (cid:42) ( (cid:126) Π mi, − Id) (cid:126)X m +1 ,j +1 − (cid:126) id∆ t m , Q mi,θ m (cid:126)κ m +1 ,j +1 i (cid:43) h Γ mi + (cid:42) (cid:126) m m +1 i , (cid:126)X m +1 ,j +1 − (cid:126) id∆ t m (cid:43) hγ m − (cid:68) ∇ s (cid:126)Y m +1 ,j +1 i , ∇ s ( (cid:126) Π mi, (cid:126)κ m +1 ,j +1 i ) (cid:69) Γ mi − (cid:68) (cid:126)f m , (cid:126) Π mi, (cid:126)κ m +1 ,j +1 i (cid:69) h Γ mi , i = 1 , RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES and continue the iteration until | λ V,m +1 ,j +1 − λ V,m +1 ,j | + | λ A,m +1 ,j +11 − λ A,m +1 ,j | + | λ A,m +1 ,j +12 − λ A,m +1 ,j | < − . We remark that the implicit scheme is chosen such that no new system matrices need to be assembledduring the fixed point iteration. In particular, all integrals are evaluated on the old interfaces Γ mi .But all the quantities that are calculated during the linear solves are treated implicitly, i.e. (cid:126)X m +1 ,( (cid:126)Y m +1 i , (cid:126)κ m +1 i , (cid:126) m m +1 i ) i =1 , (cid:126)κ m +1 γ , C Φ m +1 , as well as the Lagrange multipliers. Solution methods
Let us briefly outline how we solve the linear system (5.2a–e) in practice. First of all, similarly toour approach in [4] for the numerical approximation of surface clusters with triple junction lines, wereformulate (5.2a) as follows.On introducing the following equivalent characterization of V h (Γ m ), recall (5.1), (cid:98) V h (Γ m ) = { ( (cid:126)η , (cid:126)η ) ∈ × i =1 V h (Γ mi ) : (cid:126)η | γ m = (cid:126)η | γ m } , we can rewrite (5.2a–e) equivalently as: Find ( (cid:126)X m +11 , (cid:126)X m +12 ) ∈ (cid:98) V h (Γ m ), ( (cid:126)Y m +1 i , (cid:126) m m +1 i ) i =1 ∈ V h (Γ m ) × V h ( γ m ) × V h (Γ m ) × V h ( γ m ), (cid:126)κ m +1 γ ∈ V h ( γ m ) and C Φ m +1 ∈ V h ( γ m ) such that (cid:88) i =1 (cid:42) Q mi,θ m(cid:63) (cid:126)X m +1 i − (cid:126) id∆ t m , (cid:126)χ i (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 i , ∇ s (cid:126)χ i (cid:69) Γ mi + α Gi (cid:10) ( (cid:126) m m +1 i ) s , [ (cid:126)χ i ] s (cid:11) γ m + ς (cid:68) [ (cid:126)X m +1 i ] s , [ (cid:126)χ i ] s (cid:69) γ m + (cid:37) (cid:42) (cid:126)X m +1 i − (cid:126) id∆ t m , (cid:126)χ i (cid:43) hγ m = (cid:88) i =1 (cid:20)(cid:68) ∇ s . (cid:126)Y mi , ∇ s . (cid:126)χ i (cid:69) Γ mi − (cid:68) ( ∇ s (cid:126)Y mi ) T , D ( (cid:126)χ i ) ( ∇ s (cid:126) id) T (cid:69) Γ mi − (cid:68) [ α i | (cid:126)κ mi − κ i (cid:126)ν mi | − (cid:126)Y mi . Q mi,θ m (cid:126)κ mi )] ∇ s (cid:126) id , ∇ s (cid:126)χ i (cid:69) h Γ mi − α i κ i (cid:10) (cid:126)κ mi , [ ∇ s (cid:126)χ i ] T (cid:126)ν mi (cid:11) h Γ mi + (cid:68) (1 − θ m ) ( (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) . (cid:126)ν mi ) ∇ s (cid:126) id , ∇ s (cid:126)χ i (cid:69) h Γ mi − (cid:68) (1 − θ m ) (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) , [ ∇ s (cid:126)χ i ] T (cid:126)ν mi (cid:69) h Γ mi (cid:21) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126)κ mγ . (cid:126) m mi , (cid:126) id s . [ (cid:126)χ i ] s (cid:69) hγ m + (cid:10) (Id + P mγ ) ( (cid:126) m mi ) s , [ (cid:126)χ i ] s (cid:11) γ m (cid:21) − λ V,m (cid:88) i =1 (cid:104) (cid:126)ω mi , (cid:126)χ i (cid:105) h Γ mi − (cid:88) i =1 λ A,mi (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ i (cid:69) Γ mi ∀ ( (cid:126)χ , (cid:126)χ ) ∈ (cid:98) V h (Γ m ) (6.1)and (5.2b–e) hold, where in (6.1) we have used the fully discrete version of (4.7).The above reformulation is crucial for the construction of fully practical solution methods, as itavoids the use of the global finite element space V h (Γ m ). With the help of (6.1), it is now possibleto work with the basis of the simple product finite element space (cid:98) V h (Γ m ), on employing suitableprojections in the formulation of the linear problem. This construction is similar to e.g. the standardtechnique used for an ODE with periodic boundary conditions.We recall from [10, (4.4a–d)] the following finite element approximation for Willmore flow of asingle open surface Γ i ( t ) with free boundary conditions for ∂ Γ i ( t ). For m = 0 , . . . , M −
1, find27 . W. Barrett, H. Garcke, R. N¨urnberg ( δ (cid:126)X m +1 i , (cid:126)Y m +1 i ) ∈ V h (Γ mi ) × V h (Γ mi ), with (cid:126)X m +1 i = (cid:126) id | Γ mi + δ (cid:126)X m +1 i , and ( (cid:126)κ m +1 ∂ Γ i , (cid:126) m m +1 i ) ∈ [ V h ( ∂ Γ mi )] such that (cid:42) Q mi,θ m(cid:63) (cid:126)X m +1 i − (cid:126) id∆ t m , (cid:126)χ (cid:43) h Γ mi − (cid:68) ∇ s (cid:126)Y m +1 i , ∇ s (cid:126)χ (cid:69) Γ mi + ς (cid:68) [ (cid:126)X m +1 i ] s , (cid:126)χ s (cid:69) ∂ Γ mi + α Gi (cid:10) [ (cid:126) m m +1 i ] s , (cid:126)χ s (cid:11) ∂ Γ mi = (cid:68) ∇ s . (cid:126)Y mi , ∇ s . (cid:126)χ (cid:69) Γ mi − (cid:68) ( ∇ s (cid:126)Y mi ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ mi − α i κ i (cid:10) (cid:126)κ mi , [ ∇ s (cid:126)χ ] T (cid:126)ν mi (cid:11) h Γ mi − (cid:68)(cid:104) α i | (cid:126)κ mi − κ i (cid:126)ν mi | − (cid:126)Y mi . Q mi,θ m (cid:126)κ mi (cid:105) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ mi + (cid:68) (1 − θ m ) ( (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) . (cid:126)ν mi ) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) h Γ mi − (cid:68) (1 − θ m ) (cid:126)G mi ( (cid:126)Y mi , (cid:126)κ mi ) , [ ∇ s (cid:126)χ ] T (cid:126)ν mi (cid:69) h Γ mi + α Gi (cid:68) (cid:126)κ m∂ Γ i . (cid:126) m mi , (cid:126) id s . (cid:126)χ s (cid:69) h∂ Γ mi + α Gi (cid:10) (Id + P m∂ Γ i ) [ (cid:126) m mi ] s , (cid:126)χ s (cid:11) ∂ Γ mi − λ A,mi (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ mi ∀ (cid:126)χ ∈ V h (Γ mi ) , (6.2a) α − i (cid:68) Q mi,θ m (cid:126)Y m +1 i , Q mi,θ m (cid:126)η (cid:69) h Γ mi + (cid:68) ∇ s (cid:126)X m +1 i , ∇ s (cid:126)η (cid:69) Γ mi = (cid:10) (cid:126) m m +1 i , (cid:126)η (cid:11) h∂ Γ mi − κ (cid:104) (cid:126)ω mi , (cid:126)η (cid:105) h Γ mi ∀ (cid:126)η ∈ V h (Γ mi ) , (6.2b) (cid:68) α Gi (cid:126)κ m +1 ∂ Γ i + (cid:126)Y m +1 i , (cid:126)ϕ (cid:69) h∂ Γ mi = 0 ∀ (cid:126)ϕ ∈ V h ( ∂ Γ mi ) , (6.2c) (cid:68) (cid:126)κ m +1 ∂ Γ i , (cid:126)η (cid:69) h∂ Γ mi + (cid:68) [ (cid:126)X m +1 i ] s , (cid:126)η s (cid:69) ∂ Γ mi = 0 ∀ (cid:126)η ∈ V h ( ∂ Γ mi ) . (6.2d)The corresponding linear system from [10, (5.1)] is then given by (cid:126)A − t m (cid:126) M Q (cid:63) − (cid:126)A ς − α Gi (cid:126)A ∂ Γ , Γ (cid:126) M Q (cid:126)A − (cid:126)M ∂ Γ , Γ ( (cid:126)M ∂ Γ , Γ ) T α Gi (cid:126)M ∂ Γ
00 ( (cid:126)A ∂ Γ , Γ ) T (cid:126)M ∂ Γ (cid:126)Y m +1 i δ (cid:126)X m +1 i (cid:126)κ m +1 ∂ Γ i (cid:126) m m +1 i = [ (cid:126) B (cid:63) − (cid:126) B + (cid:126) R ] (cid:126)Y mi + ( (cid:126)A θ + (cid:126)A ς + λ A,mi (cid:126)A ) (cid:126)X mi + (cid:126)b θ − (cid:126)b α − (cid:126)A (cid:126)X mi − κ (cid:126)M (cid:126)ω mi (cid:126) − ( (cid:126)A ∂ Γ , Γ ) T (cid:126)X mi . (6.3)On replacing (cid:126)A ς with ( (cid:126)A ς + t (cid:126)M (cid:37) ), where the definition of (cid:126)M (cid:37) is clear from (6.1), and similarlyadapting the first entry in the right hand side of (6.3) to account for the term involving λ V,m , we write(6.3) as B i Z i = g i . Hence we can rewrite the linear system for (6.1), (5.2b–e) as P B B P Z Z Z (cid:126) Φ m +1 = P B g g , (6.4a)28 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES where B = B C (cid:126)M γ B C (cid:126)M γ (0 0 0 C (cid:126)M γ ) (0 0 0 C (cid:126)M γ ) 0 , (6.4b)and where (cid:126)M γ is a mass matrix on γ m . Moveover, P B and P Z are the orthogonal projections thatencode the test and trial space (cid:98) V h (Γ m ) in (6.1), i.e. they act on the first and fifth block row in (6.4b),and on the second entries of Z and Z , respectively.The system (6.4a) can be efficiently solved in practice with a preconditioned BiCGSTAB or GMRESiterative solver, where we employ the preconditioners P Z B − B −
00 0 Id P B and P Z B − P B for the cases C = 0 and C = 1, respectively. Here we recall from [10] that B and B are invert-ible. The inverses B − and B − can be computed with the help of the sparse factorization packageUMFPACK, see [16]. Similarly, the inverse B − , which existed in all our numerical tests, can also becomputed with the help of UMFPACK.In practice we note that the preconditioned Krylov subspace solvers usually take fewer than teniterations per time step to converge. We stress that the chosen preconditioners are crucial, as withoutappropriate preconditioning the iterative solvers do not converge. This suggests that the linear systems(6.4a) are badly conditioned. Numerical results
We implemented our fully discrete finite element approximations within the finite element toolboxALBERTA, see [34]. The arising systems of linear equations were solved with the help of the sparsefactorization package UMFPACK, see [16]. For the computations involving surface area preservingWillmore flow, we always employ the implicit Lagrange multiplier formulation discussed in § (cid:126)κ i , (cid:126)Y i , (cid:126) m i , i = 1 ,
2, and (cid:126)κ γ . Given the initialtriangulation Γ i , we let (cid:126) m i ∈ V h ( γ ) be such that (cid:10) (cid:126) m i , (cid:126)η (cid:11) hγ = (cid:10) (cid:126)µ i , (cid:126)η (cid:11) γ ∀ (cid:126)η ∈ V h ( γ ) , with (cid:126)µ i denoting the conormal on ∂ Γ i , i = 1 ,
2. In addition, we let (cid:126)κ i = − R (cid:126)ω i for simulations where Γ i (0) is part of a sphere of radius R , i.e. Γ i (0) ⊂ ∂B R ( (cid:126) (cid:126)κ i ∈ V h (Γ i ) to be the solution of (cid:10) (cid:126)κ i , (cid:126)η (cid:11) h Γ i + (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)η (cid:69) Γ i = (cid:10) (cid:126) m i , (cid:126)η (cid:11) hγ ∀ (cid:126)η ∈ V h (Γ i ) . Then we define (cid:126)Y i = α i [ (cid:126)κ i − κ i (cid:126)ω i ] . . W. Barrett, H. Garcke, R. N¨urnbergFigure 2. ( C : κ = κ = 0, ς = 0 .
1) A plot of (Γ mi ) i =1 at times t = 0 , . , , E m +1 ((Γ m ) i =1 ).Moreover, we let (cid:126)κ γ ∈ V h ( γ ) be such that (cid:10) (cid:126)κ γ , (cid:126)η (cid:11) hγ + (cid:68) (cid:126) id s , (cid:126)η s (cid:69) γ = 0 ∀ (cid:126)η ∈ V h ( γ ) . Throughout this section we use uniform time steps ∆ t m = ∆ t , m = 0 , . . . , M −
1, and set ∆ t = 10 − unless stated otherwise. In addition, unless stated otherwise, we fix α i = 1 and κ i = α Gi = 0, i = 1 , ς = 0. At times we will discuss the discrete energy of the numerical solutions, which, similarlyto (4.12), is defined by E m +1 ((Γ mi ) i =1 ) := (cid:88) i =1 (cid:104) α i (cid:10) | (cid:126)κ m +1 i − κ i (cid:126)ν mi | , (cid:11) h Γ mi + α Gi (cid:104)(cid:10) (cid:126)κ m +1 γ , (cid:126) m m +1 i (cid:11) hγ m + 2 π m (Γ mi ) (cid:105)(cid:105) + ς H d − ( γ m ) . Finally, we fix θ = 0 throughout, unless otherwise stated.For the visualization of our numerical results we will use the colour red for Γ m , and the colouryellow for Γ m .7.1. The C –case In Figure 2 we show the evolution of the outer shell of a torus joined with two spherical caps. Here thetwo caps make up phase 1, with the remainder representing phase 2. The initial surface Γ satisfies( J , J ) = (2048 , K , K ) = (1090 , × ×
6, i.e. up totranslations, the smallest cuboid containing Γ is [0 , . For the parameters κ = κ = 0 and ς = 0 . κ = − κ = − .
5, which is now markedly different. The same evolution with (cid:37) = 2, which shows theslowing influence of (cid:37) >
0, is shown in Figure 4. In both experiments the effect of the two differentspontaneous curvature values for the two phases can clearly be seen. The same evolution as in Figure 4,but now for surface area preserving flow, is shown in Figure 5. Here the observed relative surface arealoss is 0 . κ i , the surface area constraints, and the C –attachment condition lead to an interesting evolution. A completely different evolution is obtained30 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANESFigure 3. ( C : κ = − κ = − . ς = 0 .
1) A plot of (Γ mi ) i =1 at times t =0 , . , ,
2. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 4. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2) A plot of (Γ mi ) i =1 at times t = 0 , . , ,
2. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).when we replace surface area conservation with volume conservation. This new simulation is visualizedin Figure 6, where the observed relative volume loss is 0 . satisfies ( J , J ) = (1816 , K , K ) = (1000 , . × . × .
1. The evolution for the parameters κ = κ = 0 and ς = 1 goes towards a fournoid.We now consider surface area preserving experiments for setups where phase 1 is represented by sixor eight disconnected components on the unit sphere. For these experiments we use the time step size∆ t = 10 − and let κ = − κ = − ς = 1 and (cid:37) = 2. The initial surface Γ in Figure 8 satisfies( J , J ) = (1032 , K , K ) = (614 , . . W. Barrett, H. Garcke, R. N¨urnbergFigure 5. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2) Surface area preserving flow. Aplot of (Γ mi ) i =1 at times t = 0 , . , ,
2. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 6. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2) Volume preserving flow. A plotof (Γ mi ) i =1 at times t = 0 , . , ,
2. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).simulation with eight disconnected components for phase 1 is shown in Figure 9. The initial surface Γ satisfies ( J , J ) = (2048 , K , K ) = (1184 , . satisfies ( J , J ) = (2274 , K , K ) = (1188 , . × . × . κ = κ = − ς = 1 and (cid:37) = 2. The relative surface area loss for thisexperiment is 0 . . is made up of two halves of an approximation of the unit sphere and satisfies ( J , J ) =32 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANESFigure 7. ( C : κ = κ = 0, ς = 1) A plot of (Γ mi ) i =1 at times t = 0 , . , . , E m +1 ((Γ m ) i =1 ). Figure 8. ( C : κ = − κ = − ς = 1, (cid:37) = 2) Surface area preserving flow. A plotof (Γ mi ) i =1 at times t = 0 , . , . , .
35. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).(2274 , K , K ) = (1188 , κ = κ = 0, ς = 1 and (cid:37) = 2 is shownin Figure 11. The evolution eventually reaches a slowly shrinking disk. Choosing the parameters α G = α G = −
1, and using the time step size ∆ t = 10 − , we obtain the simulation in Figure 12.We remark that the conditions (2.7) trivially hold. Moreover, and in contrast to the C –case, anonzero Gaussian bending energy coefficient has an influence on the evolution even if α G = α G . Inthis example we observe that for a negative α G = α G , the term (cid:80) i =1 α Gi (cid:82) Γ i K i d H for the initialsphere is negative, and hence the evolution remains convex throughout, in contrast to the evolution inFigure 11. Moreover, the evolution in Figure 12 is generally slower, since large values of the Gaussiancurvatures make (cid:80) i =1 α Gi (cid:82) Γ i K i d H more negative. Repeating the computation for α G = − . W. Barrett, H. Garcke, R. N¨urnbergFigure 9. ( C : κ = − κ = − ς = 1, (cid:37) = 2) Surface area preserving flow. A plotof (Γ mi ) i =1 at times t = 0 , . , . , .
11. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 10. ( C : κ = κ = − ς = 1, (cid:37) = 2) Volume and surface area preservingflow. A plot of (Γ mi ) i =1 at times t = 0 , . , . , .
2. Below a plot of the discreteenergy E m +1 ((Γ m ) i =1 ). α G = − . (cid:82) Γ K d H decrease the energy.7.2. The C –case We remark that in the C –case, with uniform data α = α = α , κ = κ = κ and ς = (cid:37) = α G = α G = 0, our finite element approximation collapses to the scheme from [8] for the Willmore flow of34 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANESFigure 11. ( C : κ = κ = 0, ς = 1, (cid:37) = 2) A plot of (Γ mi ) i =1 at times t =0 , . , . , . ,
1. At time t = 1 the evolution has reached a disk. Below a plot ofthe discrete energy E m +1 ((Γ m ) i =1 ). Figure 12. ( C : κ = κ = 0, ς = 1, α G = α G = − (cid:37) = 2) A plot of (Γ mi ) i =1 attimes t = 0 , . , . , .
2. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).closed surfaces. Indeed, as a numerical check we confirmed that Table 1 in [8], for the approximationof the nonlinear ODE [8, (5.1)], is reproduced exactly by our implementation of the scheme (5.2a–e).35 . W. Barrett, H. Garcke, R. N¨urnbergFigure 13. ( C : κ = κ = 0, ς = 1, α G = − α G = − . (cid:37) = 2) A plot of (Γ mi ) i =1 at times t = 0 , . , . , .
3. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 14. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2) A plot of (Γ mi ) i =1 at times t = 0 , . , . , .
5. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).A repeat of the simulation in Figure 4 in the context of a C –condition on γ is shown in Figure 14,where for this experiment we use the time step size ∆ t = 10 − . The evolution goes towards a cylinderwith two round caps, which is dramatically different to the evolution in the C –case.If we project the initial surface from Figure 14 to the unit sphere, we obtain the evolution shownin Figure 15. The evolution goes towards a cylinder with two round caps. Using the same parametersas in Figure 15 to simulate surface area preserving flow, we obtain the evolution shown in Figure 16,where here we have chosen (cid:37) = 2. The evolution goes towards a more elongated cylinder with tworound caps. Here the observed relative surface area loss is 0 . θ = 0 .
05. Here theobserved relative volume loss is − . RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANESFigure 15. ( C : κ = − κ = − . ς = 0 .
1) A plot of (Γ mi ) i =1 at times t =0 , . , . , .
03. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 16. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2) Surface area preserving flow. Aplot of (Γ mi ) i =1 at times t = 0 , . , . ,
1. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ).A repeat of the simulation in Figure 7, now in the context of a C –condition on γ , is shown inFigure 18. Once again, we observe that the C –condition has a dramatic effect on the evolution.In the next experiments we investigate the possible influence of the Gaussian curvature energy. Ifwe choose the initial surface as in Figure 12, and running with κ = κ = − . ς = α G = α G = 0,then we obtain an expanding sphere, with symmetric phases 1 and 2, which approximates the solutionto the nonlinear ODE [8, (5.1)]. On the continuous level, thanks to the Gauss–Bonnet theorem, thesame solution is obtained when choosing α G = α G = 0 .
5, and this is also replicated by our numericalapproximation. Choosing α G = 0 . α G = 1, on the other hand, leads to phase 1 growing on the37 . W. Barrett, H. Garcke, R. N¨urnbergFigure 17. ( C : κ = − κ = − . ς = 0 . (cid:37) = 2, θ = 0 .
05) Volume preservingflow. A plot of (Γ mi ) i =1 at times t = 0 , . , . ,
1. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 18. ( C : κ = κ = 0, ς = 1) A plot of (Γ mi ) i =1 at times t = 0 , . , . , . E m +1 ((Γ m ) i =1 ).expanding surface. Here we remark that (2.8) clearly holds, and that reducing the relative size of phase2 is energetically favourable. See Figure 19 for the evolution.A well known phenomenon is the moving of the phase boundary in relation to the neck of a dumbbellfor different values of the Gaussian bending rigidities, see e.g. [24, § α G ∈ {− , , } ,while α G = 0 and ς = 9. For these experiments we choose ∆ t = 10 − and let (cid:37) = 1. See Figure 20 forthe different evolution. Here we can clearly see that for α G = 2 the interface moves down relative to38 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANESFigure 19. ( C : κ = κ = − . ς = 0, α G = 0 . α G = 1) A plot of (Γ mi ) i =1 attimes t = 0 , . , . , .
25. Below a plot of the discrete energy E m +1 ((Γ m ) i =1 ). Figure 20. ( C : κ = κ = 0, ς = 9, α G = 0, (cid:37) = 4) The initial shape on the left,and (Γ mi ) i =1 at time t = 0 .
01 for α G = −
2, 0 and 2, respectively.the neck of the dumbbell, while for α G = − v r ∈ { . , . , . } , where v r = 3 L (Ω )4 π ( H (Γ )4 π ) = 6 π L (Ω )( H (Γ )) , with Ω denoting the interior of Γ . In addition, the two phases are chosen such that they have a surfacearea ratio of 0 .
1. See Figure 21 for the initial shapes, where in each case we have that the initial discretesurface Γ satisfies ( J , J ) = (2274 , K , K ) = (1188 , H (Γ ) = 4 π . For theseexperiments we set ς = 9 and (cid:37) = 4. Choosing a time step size of ∆ t = 10 − , we integrate the volumeand surface area conserving flow to a final time of t = 0 .
25 and report on the obtained shapes inFigure 22. These configurations appear to agree well with the computed shapes in [29, Fig. 8].
Appendix A. Derivation of strong formulation and boundary conditions
We recall from Section 3 that our numerical method is based on the weak formulation (3.29) and(3.28a–f) of the generalized L –gradient flow of the energy E ((Γ i ( t )) i =1 ), see (2.13). It follows from39 . W. Barrett, H. Garcke, R. N¨urnbergFigure 21. The initial shapes for v r = 0 .
95, 0 .
91 and 0 .
9, respectively.
Figure 22. ( C : κ = κ = 0, ς = 9, (cid:37) = 4) A plot of (Γ mi ) i =1 at time t = 0 .
25 forthe reduced volumes v r = 0 .
95, 0 .
91 and 0 .
9, respectively.(2.3), (3.28b) and (3.16) that (cid:126) κ i . ∂ ε (cid:126)ν i = 0 , ∂ ε ( Q i,θ (cid:126) κ i ) = − (1 − θ ) κ i [ ∇ s (cid:126)χ ] T (cid:126)ν i , (cid:104) α i | (cid:126) κ i − κ i (cid:126)ν i | − Q i,θ (cid:126)y i . (cid:126) κ i (cid:105) = − α i ( κ i − κ i ) on Γ i ( t ) , i = 1 , . We recall that on the continuous level (cid:126) m i = (cid:126)µ i and that θ ∈ [0 ,
1] is a fixed parameter. Here weneed to choose θ = 0, as otherwise the two conditions in (3.22b,c) are incompatible in general. Thenthis weak formulation can be formulated as follows. Given Γ i (0), for all t ∈ (0 , T ] find Γ i ( t ) and (cid:126)y i ( t ) ∈ [ H (Γ i ( t ))] d such that (cid:68) (cid:126) V , (cid:126)χ (cid:69) Γ( t ) + (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = (cid:88) i =1 (cid:20) (cid:104)∇ s (cid:126)y i , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . (cid:126)y i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)y i ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) + α i (cid:68) ( κ i − κ i ) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) − (1 − θ ) (cid:10) κ i (cid:126)y i , [ ∇ s (cid:126)χ ] T (cid:126)ν i (cid:11) Γ i ( t ) (cid:21) − ς (cid:68) (cid:126) id s , (cid:126)χ s (cid:69) γ ( t ) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126) κ γ . (cid:126)µ i , (cid:126) id s . (cid:126)χ s (cid:69) γ ( t ) + (cid:10) P γ ( (cid:126)µ i ) s , (cid:126)χ s (cid:11) γ ( t ) (cid:21) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d , (A.1)with γ ( t ) = ∂ Γ ( t ) = ∂ Γ ( t ), (cid:126)y i = y i (cid:126)ν i + (cid:126)u i , where y i = α i ( κ i − κ i ) and (cid:126)u i . (cid:126)ν i = 0 , on Γ i ( t ) , i = 1 , . (A.2)Of course, (3.28b) implies that (cid:126)u i = (cid:126) θ ∈ (0 , (cid:126)ν i . [ ∇ s (cid:126)χ i ] T (cid:126)ν i = ([ ∇ s (cid:126)χ i ] (cid:126)ν i ) . (cid:126)ν i = (cid:126) . (cid:126)ν i = 0,it holds that(1 − θ ) (cid:10) κ i (cid:126)y i , [ ∇ s (cid:126)χ i ] T (cid:126)ν i (cid:11) Γ i ( t ) = (1 − θ ) (cid:10) κ i y i (cid:126)ν i , [ ∇ s (cid:126)χ i ] T (cid:126)ν i (cid:11) Γ i ( t ) + (1 − θ ) (cid:10) κ i (cid:126)u i , [ ∇ s (cid:126)χ i ] T (cid:126)ν i (cid:11) Γ i ( t ) = (1 − θ ) (cid:10) κ i (cid:126)u i , [ ∇ s (cid:126)χ i ] T (cid:126)ν i (cid:11) Γ i ( t ) = (cid:10) κ i (cid:126)u i , [ ∇ s (cid:126)χ i ] T (cid:126)ν i (cid:11) Γ i ( t ) (A.3)for all θ ∈ [0 , κ i are defined by (2.3), the curve curvature vector (cid:126) κ γ is given by(2.9), and the conormals (cid:126)µ i ( t ) are defined by (2.10) and satisfy C ( (cid:126)µ + (cid:126)µ ) = (cid:126)
0. In addition, we have40
RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES from (3.28c) that (cid:126)y i = − α Gi (cid:126) κ γ − C (cid:126)φ on γ ( t ) , i = 1 , . (A.4)Starting from the weak formulation (A.1), we will now recover the corresponding strong formulationtogether with the boundary conditions that are enforced by it. It follows from (A.1) and (A.3) that (cid:68) (cid:126) V , (cid:126)χ (cid:69) Γ( t ) + (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = (cid:88) i =1 (cid:104) (cid:104)∇ s ( y i (cid:126)ν i ) , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . ( y i (cid:126)ν i ) , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) [ ∇ s ( y i (cid:126)ν i )] T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) + α i (cid:10) ( κ i − κ i ) , ∇ s . (cid:126)χ (cid:11) Γ i ( t ) + (cid:104)∇ s (cid:126)u i , ∇ s (cid:126)χ (cid:105) Γ i ( t ) + (cid:104)∇ s . (cid:126)u i , ∇ s . (cid:126)χ (cid:105) Γ i ( t ) − (cid:68) ( ∇ s (cid:126)u i ) T , D ( (cid:126)χ ) ( ∇ s (cid:126) id) T (cid:69) Γ i ( t ) − (cid:10) κ i (cid:126)u i , [ ∇ s (cid:126)χ ] T (cid:126)ν i (cid:11) Γ i ( t ) (cid:21) + ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126) κ γ . (cid:126)µ i , (cid:126) id s . (cid:126)χ s (cid:69) γ ( t ) + (cid:10) P γ [ (cid:126)µ i ] s , (cid:126)χ s (cid:11) γ ( t ) (cid:21) =: (cid:88) i =1 8 (cid:88) (cid:96) =1 T ( i ) (cid:96) + ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 α Gi (cid:20)(cid:68) (cid:126) κ γ . (cid:126)µ i , (cid:126) id s . (cid:126)χ s (cid:69) γ ( t ) + (cid:10) P γ [ (cid:126)µ i ] s , (cid:126)χ s (cid:11) γ ( t ) (cid:21) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d . (A.5)In order to identify the first term on the right hand side in (A.5), we now recall (A.21) and (A.30) in[10], where we note that in our situation β = 0, and that the results there are for d = 3, but are alsotrue for d = 2, where we always assume that ς = α G = α G = 0. Hence we have that (cid:88) (cid:96) =1 T ( i ) (cid:96) = (cid:10) − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | , (cid:126)χ . (cid:126)ν i (cid:11) Γ i ( t ) + α i (cid:104) ( ∇ s κ i ) . (cid:126)µ i , (cid:126)χ . (cid:126)ν i (cid:105) γ ( t ) − α i (cid:104) ( κ i − κ i ) ( ∇ s (cid:126)ν i ) (cid:126)µ i , (cid:126)χ (cid:105) γ ( t ) − α i (cid:10) ( κ i − κ i ) , (cid:126)χ . (cid:126)µ i (cid:11) γ ( t ) + B ( i ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d , (A.6)where B ( i ) = (cid:104)∇ s . (cid:126)u i , (cid:126)χ . (cid:126)µ i (cid:105) γ ( t ) − (cid:104)∇ s (cid:126)u i , (cid:126)µ i ⊗ (cid:126)χ (cid:105) γ ( t ) − (cid:104) κ i (cid:126)u i . (cid:126)µ i , (cid:126)χ . (cid:126)ν i (cid:105) γ ( t ) + (cid:104) ( ∇ s (cid:126)u i ) (cid:126)µ i , ( (cid:126)χ . (cid:126)ν i ) (cid:126)ν i (cid:105) γ ( t ) =: (cid:88) (cid:96) =1 D ( i ) (cid:96) . (A.7)It immediately follows from (A.6) that the strong formulation of the flow equation is (cid:126) V = (cid:2) − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | (cid:3) (cid:126)ν i on Γ i ( t ) . (A.8)Collecting the boundary terms arising in (A.5) and (A.6), similarly to [10, (A.32)], gives: (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 6 (cid:88) (cid:96) =1 B ( i ) (cid:96) = ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 (cid:104) (cid:104) α i ( ∇ s κ i ) . (cid:126)µ i , (cid:126)χ . (cid:126)ν (cid:105) γ ( t ) − (cid:10) α i ( κ i − κ i ) , (cid:126)χ . (cid:126)µ i (cid:11) γ ( t ) − (cid:104) α i ( κ i − κ i ) ( ∇ s (cid:126)ν i ) (cid:126)µ i , (cid:126)χ (cid:105) γ ( t ) + α Gi (cid:68) (cid:126) κ γ . (cid:126)µ i , (cid:126) id s . (cid:126)χ s (cid:69) γ ( t ) + α Gi (cid:10) P γ (cid:126)χ s , [ (cid:126)µ i ] s (cid:11) γ ( t ) + B ( i ) (cid:105) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d . (A.9)41 . W. Barrett, H. Garcke, R. N¨urnberg We now investigate the boundary conditions arising from (A.9) in the case C = 0. To this end, werecall from (A.2), (A.4) and (2.11) that α i ( κ i − κ i ) + α Gi (cid:126) κ γ . (cid:126)ν i = 0 and (cid:126)u i = − α Gi ( (cid:126) κ γ . (cid:126)µ i ) (cid:126)µ i on γ ( t ) , i = 1 , . (A.10)Using the simplifications in (A.34)–(A.41) in [10], as well as (A.10), the right hand side of (A.9) canbe simplified to obtain (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 6 (cid:88) (cid:96) =1 B ( i ) (cid:96) = ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 (cid:10) ( − α i ( κ i − κ i ) − α Gi K i ) (cid:126)µ i + (( α i ( ∇ s κ i ) . (cid:126)µ i − α Gi ( τ i ) s ) (cid:126)ν i , (cid:126)χ (cid:11) γ ( t ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d . (A.11)It follows from (A.10) and (A.11) that the necessary boundary conditions are α i ( κ i − κ i ) + α Gi (cid:126) κ γ . (cid:126)ν i = 0 on γ ( t ) , i = 1 , , (A.12a) ς (cid:126) κ γ + (cid:88) i =1 ( − α i ( κ i − κ i ) − α Gi K i ) (cid:126)µ i + (( α i ( ∇ s κ i ) . (cid:126)µ i − α Gi ( τ i ) s ) (cid:126)ν i = (cid:37) (cid:126) V on γ ( t ) . (A.12b)In the case of surface area preservation, there is an extra term − (cid:88) i =1 λ Ai (cid:68) ∇ s (cid:126) id , ∇ s (cid:126)χ (cid:69) Γ i ( t ) = − (cid:88) i =1 λ Ai (cid:104)∇ s . (cid:126)χ, (cid:105) Γ i ( t ) = (cid:88) i =1 λ Ai (cid:104) (cid:104) κ i (cid:126)ν i , (cid:126)χ (cid:105) Γ i ( t ) − (cid:104) , (cid:126)χ . (cid:126)µ i (cid:105) γ ( t ) (cid:105) on the right hand side of (A.1), on recalling (2.14), (3.6), (3.11) and (3.25). Similarly, in the case ofvolume conservation, there is an extra term − λ V (cid:80) i =1 (cid:104) (cid:126)ν i , (cid:126)χ (cid:105) Γ i ( t ) on the right hand side of (A.1), onrecalling (2.14), a variational variant of (3.4) and (3.25). Hence overall we obtain (cid:126) V . (cid:126)ν i = − α i ∆ s κ i + α i ( κ i − κ i ) κ i − α i ( κ i − κ i ) |∇ s (cid:126)ν i | + λ Ai κ i − λ V on Γ i ( t ) , (A.13)in place of (A.8), as well as (cid:88) i =1 (cid:2) (( α i ( ∇ s κ i ) . (cid:126)µ i − α Gi ( τ i ) s ) (cid:126)ν i − ( α i ( κ i − κ i ) + α Gi K i + λ Ai ) (cid:126)µ i (cid:3) + ς (cid:126) κ γ = (cid:37) (cid:126) V on γ ( t ) , (A.14)in place of (A.12b).We now investigate the boundary conditions arising from (A.9) in the case C = 1, where we recallthat in this situation (cid:126)ν = (cid:126)ν = (cid:126)ν and (cid:126)µ = (cid:126)µ = − (cid:126)µ on γ ( t ). To this end, we obtain from (A.2), (A.4)and (2.11) that (cid:126)y − (cid:126)y = − ( α G − α G ) (cid:126) κ γ on γ ( t ), and hence α ( κ − κ ) + α G (cid:126) κ γ . (cid:126)ν = α ( κ − κ ) + α G (cid:126) κ γ . (cid:126)ν on γ ( t ) , (A.15a) (cid:126)u − (cid:126)u = − ( α G − α G ) ( (cid:126) κ γ . (cid:126)µ ) (cid:126)µ on γ ( t ) . (A.15b)We now rewrite some of the terms in (A.9). To this end, we first note that it follows from (2.11), (2.4)and (2.17) that (cid:126) κ γ . (cid:126)ν i = (cid:126) id ss . (cid:126)ν i = − (cid:126) id s . [ (cid:126)ν i ] s = II i ( (cid:126) id s , (cid:126) id s ) and κ i = II i ( (cid:126) id s , (cid:126) id s ) + II i ( (cid:126)µ i , (cid:126)µ i ) on γ ( t ) . (A.16)Moreover, it follows from (2.18) and (2.10) that[ (cid:126)ν i ] s × (cid:126) id s = − τ i (cid:126)µ i × (cid:126) id s = ( − i τ i (cid:126)ν i on γ ( t ) , (A.17)42 RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES where we have observed that [ (cid:126)ν i ] s is perpendicular to (cid:126)ν i . We also note from (A.17) and (2.10) that[ (cid:126)µ i ] s = ( − i (cid:16) [ (cid:126)ν i ] s × (cid:126) id s + (cid:126)ν i × (cid:126) id ss (cid:17) = τ i (cid:126)ν i + ( − i ( (cid:126) κ γ . (cid:126)µ i ) (cid:126)ν i × (cid:126)µ i = τ i (cid:126)ν i − ( (cid:126) κ γ . (cid:126)µ i ) (cid:126) id s on γ ( t ) . (A.18)Now it follows from (2.17) that( ∇ s (cid:126)ν i ) (cid:126)µ i . (cid:126)χ = − II i ( (cid:126)µ i , (cid:126) id s ) (cid:126)χ . (cid:126) id s − II i ( (cid:126)µ i , (cid:126)µ i ) (cid:126)χ . (cid:126)µ i , and so B ( i )3 = α i (cid:68) ( κ i − κ i ) , II i ( (cid:126)µ i , (cid:126) id s ) (cid:126)χ . (cid:126) id s + II i ( (cid:126)µ i , (cid:126)µ i ) (cid:126)χ . (cid:126)µ i (cid:69) γ ( t ) = α i (cid:68) ( κ i − κ i ) , τ i (cid:126)χ . (cid:126) id s + ( κ i − (cid:126) κ γ . (cid:126)ν i ) (cid:126)χ . (cid:126)µ i (cid:69) γ ( t ) , (A.19)where we have noted (2.18) and (A.16). It follows from (A.36) and (A.37) in [10] that B ( i )4 + B ( i )5 = α Gi (cid:68) [ (cid:126) κ γ . (cid:126)ν i τ i − ( (cid:126) κ γ . (cid:126)µ i ) s ] (cid:126) id s + (cid:2) τ i − ( (cid:126) κ γ . (cid:126)µ i ) (cid:3) (cid:126)µ i , (cid:126)χ (cid:69) γ ( t ) − α Gi (cid:104) [( τ i ) s + ( (cid:126) κ γ . (cid:126)µ i ) (cid:126) κ γ . (cid:126)ν i ] (cid:126)ν i , (cid:126)χ (cid:105) γ ( t ) . (A.20)We have from (A.7) above and (A.39) in [10] that D ( i )1 + D ( i )2 = (cid:68) ( (cid:126)u i ) s , ( (cid:126)χ . (cid:126)µ i ) (cid:126) id s − ( (cid:126)χ . (cid:126) id s ) (cid:126)µ i (cid:69) γ ( t ) . (A.21)As (cid:126)u i . (cid:126)ν i = 0, we have from (A.7), (2.17) and (A.16) that D ( i )3 + D ( i )4 = − (cid:104) κ i (cid:126)µ i + ( ∇ s (cid:126)ν i ) (cid:126)µ i , ( (cid:126)χ . (cid:126)ν i ) (cid:126)u i (cid:105) γ ( t ) = (cid:68) II i ( (cid:126)µ i , (cid:126) id s ) , ( (cid:126)u i . (cid:126) id s ) (cid:126)χ . (cid:126)ν i (cid:69) γ ( t ) + (cid:104) II i ( (cid:126)µ i , (cid:126)µ i ) − κ i , ( (cid:126)u i . (cid:126)µ i ) (cid:126)χ . (cid:126)ν i (cid:105) γ ( t ) = (cid:68) τ i (cid:126)u i . (cid:126) id s − ( (cid:126) κ γ . (cid:126)ν i ) (cid:126)u i . (cid:126)µ i , (cid:126)χ . (cid:126)ν i (cid:69) γ ( t ) . (A.22)Therefore (A.7), (A.21) and (A.22) yield that (cid:88) i =1 B ( i )6 = (cid:68) ( (cid:126)u − (cid:126)u ) s , ( (cid:126)χ . (cid:126)µ ) (cid:126) id s − ( (cid:126)χ . (cid:126) id s ) (cid:126)µ (cid:69) γ ( t ) + (cid:68) ( (cid:126)u − (cid:126)u ) . (cid:126) id s , τ ( (cid:126)χ . (cid:126)ν ) (cid:69) γ ( t ) − (cid:104) ( (cid:126)u − (cid:126)u ) . (cid:126)µ ( (cid:126) κ γ . (cid:126)ν ) (cid:126)χ . (cid:126)ν (cid:105) γ ( t ) , where τ = τ = − τ . Hence we obtain from (A.15b) and (A.18) that (cid:88) i =1 B ( i )6 = [ α Gi ] (cid:68) ( (cid:126) κ γ . (cid:126)µ ) s (cid:126) id s + ( (cid:126) κ γ . (cid:126)µ ) (cid:126)µ + ( (cid:126) κ γ . (cid:126)µ ) ( (cid:126) κ γ . (cid:126)ν ) (cid:126)ν, (cid:126)χ (cid:69) γ ( t ) . (A.23)Combining (A.9), (A.19), (A.20) and (A.23) yields, on noting (A.15a) and (2.11), that (cid:37) (cid:68) (cid:126) V , (cid:126)χ (cid:69) γ ( t ) = ς (cid:104) (cid:126) κ γ , (cid:126)χ (cid:105) γ ( t ) + (cid:88) i =1 6 (cid:88) (cid:96) =1 B ( i ) (cid:96) = (cid:10) [ α i ( ∇ s κ i )] . (cid:126)µ − [ α Gi ] τ s + ς (cid:126) κ γ . (cid:126)ν, (cid:126)χ . (cid:126)ν (cid:11) γ ( t ) + (cid:10) − [ α i ( κ i − κ i ) ] + [ α i ( κ i − κ i ) ( κ i − (cid:126) κ γ . (cid:126)ν )] + [ α Gi ] τ + ς (cid:126) κ γ . (cid:126)µ, (cid:126)χ . (cid:126)µ (cid:11) γ ( t ) ∀ (cid:126)χ ∈ [ H γ (Γ( t ))] d . This yields the boundary conditions[ α i ( ∇ s κ i )] . (cid:126)µ − [ α Gi ] τ s + ς (cid:126) κ γ . (cid:126)ν = (cid:37) (cid:126) V . (cid:126)ν on γ ( t ) , (A.24a) − [ α i ( κ i − κ i ) ] + [ α i ( κ i − κ i ) ( κ i − (cid:126) κ γ . (cid:126)ν )] + [ α Gi ] τ + ς (cid:126) κ γ . (cid:126)µ = (cid:37) (cid:126) V . (cid:126)µ on γ ( t ) , (A.24b)43 . W. Barrett, H. Garcke, R. N¨urnberg as well as (cid:37) (cid:126) V . (cid:126) id s = 0 on γ ( t ), which has no effect on the evolution of (Γ i ( t )) i =1 . Clearly, (A.15a) and(A.24a,b) yield the conditions (2.20a–c), on accounting for the surface area constraints analogously tothe case C = 0. Acknowledgements
The authors gratefully acknowledge the support of the Regensburger Universit¨atsstiftung Hans Viel-berth.
Bibliography [1] H. Abels, H. Garcke, and L. M¨uller. Local well-posedness for volume-preserving mean curvatureand Willmore flows with line tension.
Math. Nachr. , 289:136–174, 2016.[2] J. W. Barrett, H. Garcke, and R. N¨urnberg. A parametric finite element method for fourth ordergeometric evolution equations.
J. Comput. Phys. , 222:441–462, 2007.[3] J. W. Barrett, H. Garcke, and R. N¨urnberg. On the parametric finite element approximation ofevolving hypersurfaces in R . J. Comput. Phys. , 227:4281–4307, 2008.[4] J. W. Barrett, H. Garcke, and R. N¨urnberg. Parametric approximation of surface clusters drivenby isotropic and anisotropic surface energies.
Interfaces Free Bound. , 12:187–234, 2010.[5] J. W. Barrett, H. Garcke, and R. N¨urnberg. The approximation of planar curve evolutions by sta-ble fully implicit finite element schemes that equidistribute.
Numer. Methods Partial DifferentialEquations , 27:1–30, 2011.[6] J. W. Barrett, H. Garcke, and R. N¨urnberg. Elastic flow with junctions: Variational approximationand applications to nonlinear splines.
Math. Models Methods Appl. Sci. , 22:1250037, 2012.[7] J. W. Barrett, H. Garcke, and R. N¨urnberg. On the stable numerical approximation of two-phaseflow with insoluble surfactant.
M2AN Math. Model. Numer. Anal. , 49:421–458, 2015.[8] J. W. Barrett, H. Garcke, and R. N¨urnberg. Computational parametric Willmore flow with spon-taneous curvature and area difference elasticity effects.
SIAM J. Numer. Anal. , 54:1732–1762,2016.[9] J. W. Barrett, H. Garcke, and R. N¨urnberg. Finite element approximation for the dynamics offluidic two-phase biomembranes.
M2AN Math. Model. Numer. Anal. , 51:2319–2366, 2017.[10] J. W. Barrett, H. Garcke, and R. N¨urnberg. Stable variational approximations of boundary valueproblems for Willmore flow with Gaussian curvature.
IMA J. Numer. Anal. , 37:1657–1709, 2017.[11] T. Baumgart, S. Das, W. W. Webb, and J. T. Jenkins. Membrane elasticity in giant vesicles withfluid phase coexistence.
Biophys. J. , 89:1067–1080, 2005.[12] T. Baumgart, S. T. Hess, and W. W. Webb. Imaging coexisting fluid domains in biomembranemodels coupling curvature and line tension.
Nature , 425:821–824, 2003.[13] R. Choksi, M. Morandotti, and M. Veneroni. Global minimizers for axisymmetric multiphasemembranes.
ESAIM Control Optim. Calc. Var. , 19:1014–1029, 2013.44
RADIENT FLOW DYNAMICS OF TWO-PHASE BIOMEMBRANES [14] G. Cox and J. Lowengrub. The effect of spontaneous curvature on a two-phase vesicle.
Nonlin-earity , 28:773–793, 2015.[15] S. L. Das, J. T. Jenkins, and T. Baumgart. Neck geometry and shape transitions in vesicles withco-existing fluid phases: Role of Gaussian curvature stiffness vs. spontaneous curvature.
Europhys.Lett. , 86:48003, 2009.[16] T. A. Davis. Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method.
ACM Trans. Math. Software , 30:196–199, 2004.[17] K. Deckelnick, G. Dziuk, and C. M. Elliott. Computation of geometric partial differential equa-tions and mean curvature flow.
Acta Numer. , 14:139–232, 2005.[18] K. Deckelnick, H.-C. Grunau, and M. R¨oger. Minimising a relaxed Willmore functional for graphssubject to boundary conditions.
Interfaces Free Bound. , 19:109–140, 2017.[19] G. Dziuk. An algorithm for evolutionary surfaces.
Numer. Math. , 58:603–611, 1991.[20] G. Dziuk. Computational parametric Willmore flow.
Numer. Math. , 111:55–80, 2008.[21] G. Dziuk and C. M. Elliott. Finite element methods for surface PDEs.
Acta Numer. , 22:289–396,2013.[22] C. M. Elliott and B. Stinner. Modeling and computation of two phase geometric biomembranesusing surface finite elements.
J. Comput. Phys. , 229:6585–6612, 2010.[23] C. M. Elliott and B. Stinner. A surface phase field model for two-phase biological membranes.
SIAM J. Appl. Math. , 70:2904–2928, 2010.[24] C. M. Elliott and B. Stinner. Computation of two-phase biomembranes with phase dependentmaterial parameters using surface finite elements.
Commun. Comput. Phys. , 13:325–360, 2013.[25] M. Helmers. Snapping elastic curves as a one-dimensional analogue of two-component lipid bilay-ers.
Math. Models Methods Appl. Sci. , 21:1027–1042, 2011.[26] M. Helmers. Kinks in two-phase lipid bilayer membranes.
Calc. Var. Partial Differential Equa-tions , 48:211–242, 2013.[27] M. Helmers. Convergence of an approximation for rotationally symmetric two-phase lipid bilayermembranes.
Q. J. Math. , 66:143–170, 2015.[28] F. J¨ulicher and R. Lipowsky. Domain-induced budding of vesicles.
Phys. Rev. Lett. , 70:2964–2967,1993.[29] F. J¨ulicher and R. Lipowsky. Shape transformations of vesicles with intramembrane domains.
Phys. Rev. E , 53:2670–2683, 1996.[30] J. S. Lowengrub, A. R¨atz, and A. Voigt. Phase-field modeling of the dynamics of multicomponentvesicles: Spinodal decomposition, coarsening, budding, and fission.
Phys. Rev. E , 79:0311926,2009.[31] M. Mercker and A. Marciniak-Czochra. Bud-neck scaffolding as a possible driving force in ESCRT-induced membrane budding.
Biophys. J. , 108:833–843, 2015.45 . W. Barrett, H. Garcke, R. N¨urnberg [32] M. Mercker, A. Marciniak-Czochra, T. Richter, and D. Hartmann. Modeling and computing ofdeformation dynamics of inhomogeneous biological surfaces.
SIAM J. Appl. Math. , 73:1768–1792,2013.[33] J. C. C. Nitsche. Boundary value problems for variational integrals involving surface curvatures.
Quart. Appl. Math. , 51:363–387, 1993.[34] A. Schmidt and K. G. Siebert.
Design of Adaptive Finite Element Software: The Finite Ele-ment Toolbox ALBERTA , volume 42 of
Lecture Notes in Computational Science and Engineering .Springer-Verlag, Berlin, 2005.[35] S. Schmidt and V. Schulz. Shape derivatives for general objective functions and the incompressibleNavier–Stokes equations.
Control Cybernet. , 39:677–713, 2010.[36] M. E. Taylor.
Partial differential equations I. Basic theory , volume 115 of
Applied MathematicalSciences . Springer, New York, 2011.[37] F. Tr¨oltzsch.
Optimal Control of Partial Differential Equations: Theory, Methods and Applica-tions , volume 112 of
Graduate Studies in Mathematics . American Mathematical Society, Provi-dence, RI, 2010.[38] Z.-C. Tu. Challenges in theoretical investigations of configurations of lipid membranes.
Chin.Phys. B , 22:28701, 2013.[39] Z. C. Tu and Z. C. Ou-Yang. A geometric theory on the elasticity of bio-membranes.
J. Phys. A ,37:11407–11429, 2004.[40] X. Wang and Q. Du. Modelling and simulations of multi-component lipid membranes and openmembranes via diffuse interface approaches.
J. Math. Biol. , 56:347–371, 2008.[41] C. Wutz.