Nematic order on a deformable vesicle with anchoring effects

Results in Applied Mathematics 8 (2020) 100102 Contents lists available at ScienceDirect Results in Applied Mathematics journal homepage: www.elsevier.com/locate/results-in-applied-mathematics Nematic order on a deformable vesicle with anchoring effects Francisco Guillén-Gonzáleza,1, María Ángeles Rodríguez-Bellidoa,1, Giordano Tierrab,∗,1 aDepartamento de Ecuaciones Diferenciales y Análisis Numérico and IMUS, Universidad de Sevilla, Facultad de Matemáticas, Campus de Reina Mercedes, C/ Tarfia, s/n- 41012 Seville, Spain bDepartment of Mathematics, University of North Texas, Denton TX 76203, USA article info Article history: Received 29 December 2019 Received in revised form 8 March 2020 Accepted 23 March 2020 Available online 16 April 2020 Keywords: Diffuse interface Phase field Liquid crystals Cahn–Hilliard Bending energy Vesicle membrane abstract Models to represent the interactions of liquid crystalline phases with vesicles membranes have been widely studied in recent times due to its connection with biological materials. In this work we propose a new model to represent a vesicle membrane with internal nematic order whose equilibrium states depend on the competition between the bending, elastic and anchoring energies. Moreover we present a new unconditionally energy stable numerical scheme to approximate the model. Additionally, we present several numerical simulations in order to illustrate the good performance of the proposed scheme and the influence of the nematic order in the dynamics of the system. ©2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Biological vesicle membranes can be understood as closed structures made of phospholipid bilayers (two lipid monolayers) with anchored proteins that separates an aqueous compartment from a surrounding fluid, whose study has direct applications to biology, biophysics and bioengineering [1–3]. In studying structural changes, dynamics and deformation of vesicles, one must additionally understand the coupling of these materials with the environmental conditions, including interactions with different type of flows, thermotropic effects and electric or magnetic fields [3,4]. In particular, there is a growing interest on understanding the coupling of vesicle membranes with anisotropic flows of the liquid crystal type. It has been observed that liquid crystals play an important role in the dynamics of biological components. One classical example is discussed in [5], where the similarities between cholesteric liquid crystals and the extracellular matrices of fibrous tissues in plants and animals are exposed, concluding that although these biological materials can be studied using the main concepts introduced by physicists in the field of liquid crystals, it is necessary to keep in mind that living cells are present in such systems, and this fact opens considerable perspectives for new research. Moreover, liquid crystalline phases found in many biological materials, such as actin, DNA, cellulose, and collagen can be responsible for the deformation of cell membranes [6]. In [7] it is reported that the spontaneous assembly of phospholipids at planar interfaces (like the ones in the vesicle membranes) between thermotropic liquid crystals and aqueous phases produce patterned orientations of the liquid crystals that reflect the spatial and temporal organization of the phospholipids. Moreover, it is known that in membranes with nematic liquid-crystalline order there is a geometric coupling between the nematic ∗Corresponding author. E-mail addresses: [email protected] (F. Guillén-González), [email protected] (M.Á. Rodríguez-Bellido), [email protected] (G. Tierra). 1Contributed equally to this study. https://doi.org/10.1016/j.rinam.2020.100102 2590-0374/©2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons. org/licenses/by-nc-nd/4.0/). 2F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 director and the shape of the membrane and this fact is potentially useful for designing colloidal particles for photonic applications [8]. In [9], the authors report on the dynamics of continuous anchoring transitions at interfaces formed between nematic liquid crystals (LCs, 4´ -pentyl-4-cyanobiphenyl (5CB)) and immiscible aqueous phases that are induced by either non-specific or specific interactions between phospholipid vesicles and proteins adsorbed at the liquid crystalline interfaces. There are also evidences that liquid crystalline materials, such as actin or tubulin networks, are known to be capable of deforming the shape of cells [10]. The results presented in [10] indicate that, depending on its elastic properties, the liquid crystal is indeed able to deform the vesicle until it reaches an equilibrium, anisotropic shape. The magnitude of the deformation is determined by a balance of elastic and surface forces and these predictions are confirmed by their experimental observations of spindle-like shapes in experiments with giant unilamellar vesicles with planar anchoring. Another interesting relation between liquid crystals and vesicle membranes has been presented in [11], where it has been shown that when a thin film of active, nematic microtubules and kinesin motor clusters is confined on the surface of a vesicle, four 1 2-topological defects oscillate in a periodic manner between tetrahedral and planar arrangements. In fact, the authors present a theoretical description of nematics, coupled to the relevant hydrodynamic equations to explain the dynamics of active nematic shells. It is known that the equilibrium shapes assumed by vesicles membranes (without interaction with external fields) correspond with minimizers of different surface energies, such as the bending elastic energy, and this energy can be expressed using the diffuse interface approach, introducing a phase field unknown to localize the interface [12–15]. There are already several works in the literature to study analytically and numerically the dynamics of vesicle membranes using the phase field approach with and without coupling with external fields (check [16–20] and the references therein). To point out some interesting ideas, we refer the reader to [21], where the authors explore a wide variety of patterns of closed surfaces that minimize the elastic bending energy with fixed surface area and volume by constructing phase-field functionals of bending energy with penalty terms for the constraints, although no energy stability for the schemes is provided. Another interesting work is [22], where the authors present an unconditionally energy stable numerical scheme for a vesicle membrane model that satisfies exactly the conservation of volume constraint and penalizes the surface area constraint. There are also several works studying the coupling of vesicles with external fields. In [23] a phase field method is developed to investigate the morphological evolution of a vesicle in an electric field, taking into account coupled mechanical and electric effects such as bending, osmotic pressure, surface tension, flexoelectricity, and dielectricity of the membrane, studying in detail the morphological evolution of an axisymmetric vesicle under an electric field. Another approach is considered in [24], where the authors investigate using a lattice-Boltzmann algorithm the behavior of 2D deformable vesicles immersed in a nematic liquid crystal assuming homeotropic anchoring and they present the resulting equilibrium shapes for a range of surface elasticities, and investigate the interactions between pairs of vesicles. The motivation of this work comes from studies like [3], where it has been shown that vesicles with nematic internal order are realistic situations. To represent this situation we consider an extension of the vesicle membrane model studied in [22] (that comes from the interesting work presented in [12]). The difference with the previous approach is that we are interested in studying the deformation of membrane-bound vesicles considering that inside the membrane there is component with a preferred orientation of its molecules, as studied in [6]. In fact, we are going to consider that the component inside the membrane is a nematic liquid crystal whose orientation inside the membrane and its interaction with the membrane (anchoring effects) will determine the dynamics and equilibrium configurations of the system. Additionally, the coupling between the phase field part and the liquid crystal part is derived following similar arguments to the ones presented in [25], where the authors presented and studied a model for complex fluids composed by the mixture between isotropic (Newtonian fluid) and nematic (liquid crystal) flows taking into account anchoring effects of the liquid crystal molecules on the interface between both fluids. In fact, in recent years there have been many related and relevant works studying the interactions between isotropic and nematic flows, see [26–28] and the references therein. The paper is organized as follows: In Section 2we present the model that we are considering and the main ideas to derive it. Section 3is devoted to design new numerical splitting schemes and show its unconditionally energy stability. In Section 4we present several numerical experiments in order to show the validity of our approach and to show the dynamics of the system. Finally, we state the conclusions of our work in Section 5. 2. The model We denote by φthe phase field variable which is used to localize the interior (φ=1) and the exterior (φ= −1) of the membrane. The dynamics of the membrane is derived through the energetic variational approach with respect to the total energy (Etot(φ, d)) of the system, that relates the bending, the nematic and the anchoring energies. The bending energy is associated to the phase field function used to localize the membrane, with ε > 0 being a parameter related with the interfacial width of the membrane. In particular, the bending energy is defined as Eben(φ):= ε 2∫Ω(∆φ−1 ε2G(φ))2 dx=ε 2∫Ω ω2dx,(2.1) where ω:= −∆φ+1 ε2G(φ),G(φ):= F′(φ)−εk(x)H′(φ) (2.2) F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 3 with F(φ):= 1 4(φ2−1)2,H(φ):= 1 3φ3−φ, and k(x) is a given function representing the spontaneous curvature. The volume and surface area of the vesicle are defined as A(φ):= ∫Ω φdxand B(φ):= ∫Ω(ε 2|∇φ|2+1 εF(φ))dx, respectively. We want to develop a model such that the volume (A(φ)) and surface area (B(φ)) of the vesicles remain constant in time, that is, we need to enforce somehow these constraints in our model. We are going to follow the approach considered in [22], imposing volume conservation exactly by considering a H−1-gradient flow (Cahn–Hilliard-type model) while the conservation of the surface area will be approximated via introducing a penalization term in the energy. In fact, the total energy of the system that we are going to consider is: Etot(φ, d)=Ebp(φ)+λnemEnem(d, φ)+λanch Eanch(d, φ),(2.3) where Ebp(φ) denotes the penalized bending energy, Enem(d, φ) denotes the elastic energy due to the nematic liquid crystal (that also contains a penalization part related with the unitary constraint of the director vector) and Eanch(d, φ) denotes the anchoring energy that represents the influence of the interfacial effects on the orientation of the nematic liquid crystal molecules in the membrane. Moreover, parameters λnem >0 and λanch >0 are introduced to balance the effect of each energy in the system. In particular, the energy terms reads: Ebp(φ):= Eben(φ)+1 2η(B(φ)−β)2,(2.4) with η > 0 being the penalization parameter and β > 0 the desired surface area of the system. Moreover, Enem(d, φ)=∫Ω I(φ)(1 2|∇d|2+P(d))dx, and finally, depending on the anchoring effect considered: Eanch(d, φ)=1 2∫Ω(δ1|d|2|∇φ|2+δ2|d·∇φ|2)dx, with (δ1, δ2)={(0,0) no anchoring , (0,1) parallel anch. , (1,−1) homeotropic anch. . (2.5) The functional P(d) is considered as the following double-well potential whose minimums (and consequently their equilibrium states) are located at ±1: P(d)=1 4η2 d (|d|2−1)2,(2.6) with ηd>0 and we represent its derivative as p(d):= P′(d). We have introduced the volume fraction of the liquid crystal (I(φ)∈ [0,1]) as in [25] with its derivative represented by i(φ):= I′(φ). In fact, we consider the following interpolation function: I(φ):= ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 16 (φ+1)3(3φ2−9φ+8) if φ∈(−1,1) , 1 if φ≥1, 0 if φ≤ −1, (2.7) and its derivative is defined as i(φ):= I′(φ)=⎧ ⎪ ⎨ ⎪ ⎩ 15 16 (φ+1)2(φ−1)2if φ∈(−1,1) , 0 in other case . 4F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 Now, we are going to derive the coupled system. This can be done combining ideas from the Least Action Principle (LAP) and the Maximum Dissipation Principle (MDP) [29,30], arriving at the following PDE system: ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ dt+γnem (δEtot δd)=0, φt−∇ ·(γben∇δEtot δφ )=0, (2.8) with γnem, γben >0. The expressions for each variational derivative in (2.8) will be introduced as two new variables: z:= δEtot δd=λnem δEnem δd+λanch δEanch δd =λnem(−∇ ·(I(φ)∇d)+I(φ)p(d))+λanch δEanch δd,(2.9) and µ:= δEtot δφ =δEbp δφ +λnem δEnem δφ +λanch δEanch δφ = −ε∆ω+1 εG′(φ)ω+1 η(B(φ)−β)(−ε∆φ+1 εF′(φ)) +λnemI′(φ)(1 2|∇d|2+P(d))+λanch δEanch δφ , where the anchoring terms will depend on the case considered ((δ1, δ2) as in (2.5)): δEanch δd=δ1|∇φ|2d+δ2(d·∇φ)∇φ , (2.10) and δEanch δφ = ∇ ·(δ1|d|2∇φ+δ2(d·∇φ)d).(2.11) By using the variational derivative variables zand µ, system (2.8) can be reformulated as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ dt+γnemz=0, λnem[−∇ ·(I(φ)∇d)+I(φ)p(d)]+λanch δEanch δd−z=0, φt−∇ ·(γben∇µ)=0, −ε∆ω+1 εG′(φ)ω+1 η(B(φ)−β)(−ε∆φ+1 εF′(φ)) +λnemI′(φ)(1 2|∇d|2+P(d))+λanch δEanch δφ −µ=0, εω +ε∆φ−1 εG(φ)=0. (2.12) where (δEanch/δd) and (δEanch/δφ) were previously defined in (2.10) and (2.11), respectively. The PDE system (2.12) is supplemented with the following initial and boundary conditions: d|t=0=d0, φ|t=0=φ0in Ω, (I(φ)∇d)·n⏐⏐∂Ω=0,∇µ·n|∂Ω=0 in (0,T), φ|∂Ω= −1,∇φ·n|∂Ω=0,in (0,T), (2.13) where ndenotes the outwards normal vector to the boundary ∂Ω. Lemma 2.1. System (2.12) complemented with the initial and boundary conditions proposed in (2.13) satisfies the following dissipative energy law, d dt Etot(φ(t),d(t)) +γben ∥∇µ(t)∥2 L2+γnem ∥z(t)∥2 L2=0.(2.14) Proof. Testing (2.12)1by z,(2.12)2by dt,(2.12)3by µ,(2.12)4by φt, ((2.12)4)tby ω, and adding these relations we easily derive the energy law (2.14).■ 3. Numerical scheme The aim of this section is to design unconditionally energy-stable schemes for approximating the system (2.12). Hereafter (·,·)denotes the L2(Ω)-scalar product. F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 5 3.1. A generic FE space-discrete scheme Let Dh×Zh×Φh×Mh×Wh⊂H1(Ω)×L2(Ω)×H1(Ω)×H1(Ω)×H1(Ω) be conformed finite element spaces associated to a regular and quasi-uniform triangulation Thof the domain Ωwhose polyhedric boundary is denoted by ∂Ω. For the sake of simplicity we skip the use of the subscript hto denote functions that are discrete in space. Then the problem reads: Find (d(t),z(t), φ(t), µ(t), ω(t)) ∈Dh×Zh×Φh×Mh×Wh such that d⏐⏐t=0=PDhd0, φ⏐⏐t=0=PΦhφ0,in Ω,(3.15) (with PXdenoting the L2-projection into the space X) and ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (dt,¯ z)+γnem(z,¯ z)=0, λnem(I(φ)∇d,∇¯ d)+λnem(I(φ)p(d),¯ d)+λanch (δEanch δd,¯ d)−(z,¯ d)=0, (φt,¯µ)+γben (∇µ, ∇ ¯µ)=0, ε(∇ω, ∇¯ φ)+1 ε(G′(φ)ω, ¯ φ)+1 η(B(φ)−β)(ε(∇φ, ∇¯ φ)+1 ε(F′(φ),¯ φ)) +λnem (i(φ)(1 2|∇d|2+P(d)),¯ φ)+λanch (δEanch δφ ,¯ φ)−(µ, ¯ φ)=0, ε(ω, ¯ω)−ε(∇φ, ∇ ¯ω)−1 ε(G(φ),¯ω)=0, (3.16) for any (¯ d,¯ z,¯ φ, ¯µ, ¯ω)∈Dh×Zh×Φh×Mh×Wh. Lemma 3.1. Any solution (d(t),z(t), φ(t), µ(t), ω(t)) of the space-discrete scheme (3.16) satisfies the following space-discrete version of the energy law (2.14): d dt˜ Etot(φ(t),d(t), ω(t)) +γben ∥∇µ(t)∥2 L2+γnem ∥z(t)∥2 L2=0,(3.17) where the following modified energy appears ˜ Etot(φ(t),d(t), ω(t)) =˜ Ebp(φ(t), ω(t)) +λnemEnem(d(t), φ(t)) +λanch Eanch(d(t), φ(t)) , with ˜ Ebp(φ(t), ω(t)) =∫Ω ω2(t)dx+1 2η(B(φ(t)) −β)2. Proof. Taking as test functions (¯ d,¯ z,¯ φ, ¯µ, ¯ω)=(d(t),z(t), φ(t), µ(t), ω(t)) and making the time derivative of ((3.16))5, we arrive at (3.17).■ 3.2. Fully-discrete scheme for the vesicle problem For simplicity, we assume a uniform partition of the time interval [0,T]:tn=n∆t, with ∆t=T/Ndenoting the time step, δtdenoting the discrete time derivative and an+1 2denoting the midpoint approximation like δtan+1:= an+1−an ∆tand an+1 2:= an+1+an 2.(3.18) We consider a first order semi-implicit finite difference time scheme to approximate (3.16), where we split the computation of the system into two sub-steps: Initialization: Let φ0,d0given in (3.15) and ω0∈Whsuch that ε(ω0,¯ω)−ε(∇φ0,∇ ¯ω)−1 ε(G(φ0),¯ω)=0∀ ¯ω∈Wh.(3.19) Step n+1: Given (φn,dn, wn)∈Φh×Dh×Wh. 6F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 Substep 1: Find (dn+1,zn+1)∈Dh×Zhsuch that, for each (¯ d,¯ z)∈Dh×Zh ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ (dn+1−dn ∆t,¯ z)+γnem(zn+1,¯ z)=0, λnem(I(φn)∇dn+1,∇¯ d) +λnem(I(φn)p∆t(dn+1,dn),¯ d)+λanch (Λd(dn+1, φn),¯ d)−(zn+1,¯ d)=0, (3.20) where Λd(dn+1, φn) and p∆t(dn+1,dn) denote first order approximations of δEanch δd(d(tn+1), φ(tn+1)) and p(d(tn+1)), respectively. In fact, Λd(d, φ):= δEanch δd(d, φ)=δ1|∇φ|2d+δ2(d·∇φ)∇φ, with (δ1, δ2) chosen as in (2.5). Substep 2: Find (φn+1, µn+1, ωn+1)∈Φh×Mh×Whsuch that for any (¯ φ, ¯µ, ¯ω)∈Φh×Mh×Wh: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (φn+1−φn ∆t,¯µ)+γben(∇µn+1,∇ ¯µ)=0, ε(∇ωn+1,∇¯ φ)+1 ε(G∆t sec(φn+1, φn)ωn+1,¯ φ)+λanch (Λφ(dn+1, φn+1),∇¯ φ) +1 η(B(φn+1)+B(φn) 2−β)(ε(∇φn+1 2,∇¯ φ)+1 ε(F∆t sec (φn+1, φn),¯ φ)) +λnem (i∆t(φn+1, φn)[1 2|∇dn+1|2+P(dn+1)],¯ φ)−(µn+1,¯ φ)=0, ε(ωn+1,¯ω)−ε(∇φn+1,∇ ¯ω)−1 ε(G(φn+1),¯ω)=0, (3.21) where i∆t(φn+1, φn) and −∇ · Λφ(dn+1, φn+1) represent first order approximations of i(φ(tn+1)) and δEanch δφ (d(tn+1), φ(tn+1)), respectively. In fact, Λφ(d, φ):= δ1|d|2∇φ+δ2(d·∇φ)d,(3.22) with the values of (δ1, δ2) defined in (2.5). Moreover, F∆t sec (φn+1, φn) and G∆t sec(φn+1, φn) denote secant type approximations (see [31–33] and the references therein for different ways of handling these type of potentials): F∆t sec (φn+1, φn):= F(φn+1)−F(φn) φn+1−φn=1 4(φn+1+φn)((φn+1)2+(φn)2−2) (3.23) and G∆t sec(φn+1, φn):= G(φn+1)−G(φn) φn+1−φn=(φn+1−φn)2+3(φn+1φn−1−εk(x)(φn+1+φn)).(3.24) Lemma 3.2. Scheme (3.21) satisfies the conservation of mass, that is, ∫Ω φn+1=∫Ω φn=. . . =∫Ω φ0.(3.25) Proof. Testing (3.21) by ¯µ=1. ■ 3.3. Energy stability Lemma 3.3. Scheme (3.20)–(3.21) satisfies the following discrete dissipative energy law (which is a discrete version of the energy law (2.14)), δt˜ Etot(ωn+1, φn+1,dn+1)+γben∥∇µn+1∥2 L2+γnem∥zn+1∥2 L2+NDn+1 ω+NDn+1 elast +NDn+1 penal +NDn+1 interp +NDn+1 anch =0, (3.26) F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 7 where the numerical dissipation terms are: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ NDn+1 ω=ε∆t 2∥δtωn+1∥2 L2 NDn+1 elast =λnem ∆t 2∫Ω I(φn)⏐⏐δt∇dn+1⏐⏐2dx, NDn+1 penal =λnem ∫Ω I(φn)(p∆t(dn+1,dn)·δtdn+1−δtP(dn+1))dx, NDn+1 interp =λnem ∫Ω(|∇dn+1|2 2+G(dn+1))(i∆t(φn+1, φn)δtφn+1−δtI(φn+1))dx, (3.27) and NDn+1 anch =λanch ∆t 2∫Ω[δ1(|δtdn+1|2|∇φn|2+|dn+1|2|δt∇φn+1|2)+δ2(|δtdn+1·∇φn|2+|dn+1·∇δtφn+1|2)]dx, with the values of (δ1, δ2)depending on the type of anchoring defined in (2.5). Proof. For the sake of simplicity we will only show the case of homeotropic anchoring (the parallel anchoring case can be studied using the same arguments while the no anchoring case is just a trivial generalization). Taking (¯ z,¯ d)=(zn+1, δtdn+1) in (3.20), we obtain: γnem∥zn+1∥2 L2+NDn+1 elast +NDn+1 penal +λnem ∫Ω I(φn)δt(|∇dn+1|2 2+P(dn+1))dx +λanch 2∆tδ1∫Ω[|∇φn|2|dn+1|2−|∇φn|2|dn|2]dx +λanch 2∆tδ2∫Ω[|dn+1·∇φn|2−|dn·∇φn|2]dx +λanch 2∆t∫Ω[δ1|δtdn+1|2|∇φn|2+δ2|δtdn+1·∇φn|2]dx=0. (3.28) On the other hand, taking ( ¯µ, ¯ φ)=(µn+1, δtφn+1) in (3.21)1,2, adding the resulting relations and taking into account the definition of F∆t sec (φn+1, φn) and G∆t sec(φn+1, φn) given in (3.23) and (3.24), that is, 1 ε(G∆t sec(φn+1, φn)ωn+1, δtφn+1)=1 ε∆t(G(φn+1)−G(φn), ωn+1), (ε(∇φn+1 2,∇δtφn+1)+1 ε(F∆t sec (φn+1, φn), δtφn+1))=B(φn+1)−B(φn), we obtain γben∥∇µn+1∥2 L2+ε(∇ωn+1,∇δtφn+1) +1 ε∆t(G(φn+1)−G(φn), ωn+1)+1 2ηδt(B(φn+1)−β)2 +NDn+1 interp +λnem ∫Ω(|∇dn+1|2 2+G(dn+1))δtI(φn+1)dx +λanch 2∆tδ1∫Ω[|∇φn+1|2|dn+1|2−|∇φn|2|dn+1|2]dx +λanch 2∆tδ2∫Ω[|dn+1·∇φn+1|2−|dn+1·∇φn|2]dx +λanch 2∆t∫Ω[δ1|dn+1|2|∇δtφn+1|2+δ2|dn+1·∇δtφn+1|2]dx=0. Subtracting (3.21)3and (3.21)3for previous time step, and dividing by ∆t, we obtain ε(δtωn+1,¯ω)−ε(∇δtφn+1,∇ ¯ω)−1 ε∆t(G(φn+1)−G(φn),¯ω)=0.(3.29) Then taking ¯ω=ωn+1in (3.29), we obtain ε 2δt∥ωn+1∥2 L2+NDn+1 ω−ε(∇δtφn+1,∇ωn+1)−1 ε∆t(G(φn+1)−G(φn), ωn+1)=0.(3.30) 8F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 Adding relations (3.29) and (3.30), the terms 1 ε∆t(G(φn+1)−G(φn), ωn+1)and ε(∇δtφn+1,∇ωn+1)cancel, hence we derive: δtEbp(ωn+1, φn+1)+γben∥∇µn+1∥2 L2+NDn+1 ω+NDn+1 interp +λnem ∫Ω(|∇dn+1|2 2+G(dn+1))δtI(φn+1)dx +λanch 2∆tδ1∫Ω[|∇φn+1|2|dn+1|2−|∇φn|2|dn+1|2]dx +λanch 2∆tδ2∫Ω[|dn+1·∇φn+1|2−|dn+1·∇φn|2]dx +λanch 2∆t∫Ω[δ1|dn+1|2|∇δtφn+1|2+δ2|dn+1·∇δtφn+1|2]dx=0. (3.31) Adding relations (3.28) and (3.31), we arrive at the desired relation presented in (3.26).■ Remark 3.4. From (3.26) it is clear that scheme (3.20)–(3.21) is unconditional energy-stable if we consider approximations of p∆t(dn+1,dn) and i∆t(φn+1, φn) such that NDn+1 penal ≥0 and NDn+1 interp ≥0.(3.32) There are several ways of achieving this goal, but we will consider the approximations introduced in [25] (where it is shown that these approximations satisfy (3.32)) p∆t(dn+1,dn)=˜ p(dn)+1 2∥˜ p′(dn)∥∞(dn+1−dn), i∆t(cn+1,cn)=i(cn)+5√3 12 (cn+1−cn), (3.33) with ˜ p(d)=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 2 η2 d (|d|−1)d |d|if |d| ≥ 1, 1 η2 d (|d|2−1) dif |d| ≤ 1. (3.34) Lemma 3.5. Scheme (3.20)–(3.21) considering the approximations of p∆t(dn+1,dn)and i∆t(φn+1, φn)presented in (3.33) is unconditionally energy-stable. 4. Simulations In this section we present numerical results to show the effectiveness of the numerical schemes derived in the paper. In particular, we have considered the scheme presented in (3.20)–(3.21) considering the approximations of p∆t(dn+1,dn) and i∆t(φn+1, φn) presented in (3.33) and for the sake of simplicity we only consider the case with no spontaneous curvature (k(x)=0). It is important to remark that Substep 1 is a linear scheme while Substep 2 is a nonlinear one. 4.1. Iterative scheme In order to approximate (3.21), we use a variant of the Newton’s method considered in [22] to approximate the solution: Initialization (l =0): Set (φ0,n, µ0,n, ω0,n)=(φn, µn, ωn)∈Φh×Mh×Wh Algorithm: Given (φl,n, µl,n, ωl,n)∈Φh×Mh×Wh, compute (φl+1,n, µl+1,n, ωl+1,n)∈Φh×Mh×Whsuch that for all (¯ φ, ¯µ, ¯ω)∈Φh×Mh×Wh: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (φl+1,n−φn ∆t,¯µ)+γben(∇µl+1,n,∇ ¯µ)=0, ε(∇ωl+1,n,∇¯ φ)+1 ε(G∆t sec(φl,n, φn)ωl+1,n,¯ φ)+λanch (Λφ(dn+1, φl+1,n),∇¯ φ) +1 η(B(φl,n)+B(φn) 2−β)(ε 2(∇(φl+1,n+φn),∇¯ φ)+1 ε(ˆ F∆t sec (φl+1,n, φl,n, φn),¯ φ)) +λnem (i∆t(φl+1,n, φn)[1 2|∇dn+1|2+P(dn+1)],¯ φ)−(µl+1,n,¯ φ)=0, ε(ωl+1,n,¯ω)−ε(∇φl+1,n,∇ ¯ω)−1 ε(ˆ G(φl+1,n, φl,n),¯ω)=0, (4.35) F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 9 Table 1 Parameters. Ω[0,T]h dt λnem λanch γnem γben ε ηdηtol [0,1]2[0,1.0e−3]1/120 1.0e−710.1 10.1 0.5 0.01 0.01 0.075 1.0e−31.0e−5 Table 2 Values considered for βand d0. Case 1 Case 2 Case 3 Case 4 β0.75 ×B(φ0) 0.75 ×B(φ0) 1.25 ×B(φ0) 1.25 ×B(φ0) d0(0,1) (1,0) (0,1) (1,0) where ˆ F∆t sec (φl+1,n, φl,n, φn)=F∆t sec (φl,n, φn)+(F∆t sec )′(·, φn)⏐⏐⏐φl,n(φl+1,n−φl,n) and ˆ G(φl+1,n, φl,n)=G(φl,n)+G′(φl,n)(φl+1,n−φl,n). To iterate until ∥(φl+1,n, µl+1,n, ωl+1,n)−(φl,n, µl,n, ωl,n)∥L2×L2×L2 ∥(φl+1,n, µl+1,n, ωl+1,n)∥L2×L2×L2≤tol, with tol >0 being a certain given tolerance. All the simulations have been carried out in 2Dusing FreeFem++ software [34], with an unstructured mesh and all the discrete spaces have been discretized as P1. The discrete and physical parameters are presented in Table 1. We focus on two different settings in the simulations, the first one consists on imposing as the surface area a value lower than the initial surface area (β < B(φ0)) while in the second one we impose that the desired surface area is higher than the initial surface area (β > B(φ0)). The idea in both situations is to consider the same initial condition for φ(a biconcave vesicle) with different initial orientations of the nematic liquid crystal (d0) and see how these initial conditions play a role in the dynamics of the whole system depending on the anchoring effects considered. In fact, we consider four different cases that are presented in Table 2. 4.2. Homeotropic anchoring effects (δ1, δ2)=(1,−1) We can observe on Figs. 1 and 6the dynamics for Cases 1-2 and Cases 3-4, respectively. In Fig. 1 the vesicles are evolving in such a way that the system is reducing the amount of interface because the imposed (desired) value of the surface area is lower than the initial one, although it seems that the dynamics does not change significantly when different values for d0are considered. It is clear from the graphs presented in Figs. 2–3that the total energy is decreasing in both cases and in Figs. 4–5we can observe that the value of the surface area is decreasing and that ∫Ωφalways remain constant and therefore the volume of the vesicles (volume=(1/2) ∫Ω(φ+1)) also remain constant as expected. On the other hand, when we impose that the surface area of the vesicle has to increase (Cases 3 and 4), we can observe in Fig. 6 that the dynamics of the vesicle are strongly related with the initial configuration of the director vector d0. In fact, in one case the system achieve a pronounced biconcave shape (Case 3) and in the other the system evolve to an elongated shape (Case 4). In both cases the total energy of the system decreases as expected, as it can be observed in Figs. 7–8. Moreover, it can also be observed in Figs. 9–10 that the volume of the vesicles remain constant in both cases as well as the fact that the value of the surface area is increasing as expected in both cases. 4.3. Parallel anchoring effects (δ1, δ2)=(0,1) We can observe on Figs. 11 and 16 the dynamics using parallel anchoring for Cases 1-2 and Cases 3-4, respectively. We can observe in Fig. 11 that the vesicles are evolving in such a way that the system is reducing the amount of interface because the imposed (desired) value of the surface area is lower than the initial one, although it seems that the dynamics does not change significantly when different values for d0are considered (as it happened in the homeotropic case). It is clear from the graphs presented in Figs. 12–13 that the total energy is decreasing in both cases and in Figs. 14–15 we can observe that ∫Ωφremains constant in time (therefore also the volume remains constant) and the value of the surface area is decreasing as expected. On the other hand, we have imposed that the surface area of the vesicle has to increase (Cases 3 and 4) in Fig. 16 and we can observe that the dynamics of the vesicle are strongly related with the initial configuration of the director vector d0. In fact, in one case the system achieve an elongated shape (Case 3) and in the other the system evolve to a pronounced biconcave shape (Case 4). In both cases the total energy of the system decreases as expected, as it can be observed in 16 F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 Fig. 19. Evolution of ∫Ωφ(left) and Surface Area (right) of the vesicle for Case 3. Fig. 20. Evolution of ∫Ωφ(left) and Surface Area (right) of the vesicle for Case 4. Several numerical computations using this new numerical scheme to approximate the system have been reported, showing the good performance of the proposed method considering different initial conditions and studying the different dynamics that homeotropic and parallel anchoring effects produce. In fact, the resulting dynamics are different from the ones without considering internal nematic order, that were presented in [22]. In all the cases the energy-stability is achieved and we illustrate how the different choice of the initial configuration of the nematic director vector (d0) characterize the behavior of the system. Finally, the obtained results suggest that the presented numerical scheme is the right choice to do a systematically study of all the possible configurations that can be obtained in 2D, due to its efficiency and accuracy. Moreover, the numerical scheme can be trivially extended to study the dynamics of vesicles with internal nematic order in 3D. Declaration of competing interest No author associated with this paper has disclosed any potential or pertinent conflicts which may be perceived to have impending conflict with this work. For full disclosure statements refer to https://doi.org/10.1016/j.rinam.2020.100102. Acknowledgment The research of Francisco Guillén-González and María Ángeles Rodríguez-Bellido has been supported by Proyecto PGC2018-098308-B-I00,financiado por FEDER/Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación, Spain. References [1] Campelo F, Hernández-Machado A. Shape instabilities in vesicles: A phase-field model. Eur Phys J Spec Top 2007;143:101–8. [2] Campelo F, Cruz A, Pérez-Gil J, Vázquez L, Hernández-Machado A. Phase-field model for the morphology of monolayer lipid domains. Eur Phys J E 2012;35:49. F. Guillén-González, M.Á. Rodríguez-Bellido and G. Tierra / Results in Applied Mathematics 8 (2020) 100102 17 [3] Hocine S, Brulet A, Lin J, Yang J, Di Cicco A, Bouteillerce L, et al. Structural changes in liquid crystal polymer vesicles induced by temperature variation and magnetic fields. Soft Matter 2011;7:2613–23. [4] Barthes-Biesel D. Motion and deformation of elastic capsules and vesicles in flow. Annu Rev Fluid Mech 2016;48:25–52. [5] Bouligand Y. Liquid crystals and biological morphogenesis: Ancient and new questions. C R Chim 2008;11:281–96. [6] Rofouie P, Pasini D, Rey AD. Morphology of elastic nematic liquid crystal membranes. Soft Matter 2017;13:5366. [7] Brake JM, Daschner MK, Luk YY, Abbott NL. Biomolecular interactions at phospholipid-decorated surfaces of liquid crystals. Science 2003;302(5653):2094–7. [8] Nguyen TS, Geng J, Selinger RLB, Selinger JV. Nematic order on a deformable vesicle: theory and simulation. Soft Matter 2013;9:8314. [9] Tan LT, Abbott NL. Dynamic anchoring transitions at aqueous-liquid crystal interfaces induced by specific and non-specific binding of vesicles to proteins. J Colloid Interface Sci 2015;449:452–61. [10] Zhang R, Zhou Y, Martinez-Gonzalez JA, Hernandez-Ortiz JP, Abbott NL, de Pablo JJ. Controlled deformation of vesicles by flexible structured media. Sci Adv 2016;2(8). e1600978. [11] Zhang R, Zhou Y, Rahimi M, de Pablo JJ. Dynamic structure of active nematic shells. Nature Commun 2016;7:13483. [12] Du Q, Liu C, Ryham R, Wang X. Energetic variational approaches in modeling vesicle and fluid interactions. Physica D 2009;238:923–30. [13] Du Q, Liu C, Wang X. A phase field approach in the numerical study of the elastic bending energy for vesicle membranes. J Comput Phys 2004;198:450–68. [14] Du Q, Liu C, Wang X. Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions. J Comput Phys 2006;212:757–77. [15] Du Q, Wang X. Convergence of numerical approximations to a phase field bending elasticity model of membrane deformations. Int J Numer Anal Model 2007;4:441–59. [16] Climent-Ezquerra B, Guillén-González F. Convergence to equilibrium of global weak solutions for a Cahn-Hilliard-Navier–Stokes vesicle model. Z Angew Math Phys 2019;70:125. [17] Du Q, Li M, Liu C. Analysis of a phase field navier-stokes vesicle-fluid interaction model. Discrete Contin Dyn Syst 2007;8:539–56. [18] Chen R, Ji G, Yang X, Zhang H. Decoupled energy stable schemes for phase-field vesicle membrane model. J Comput Phys 2015;302:509–23. [19] Du Q, Zhang J. Adaptive finite element method for a phase field bending elasticity model of vesicle membrane deformations. SIAM J Sci Comput 2008;30:1634–57. [20] Du Q, Zhu L. Analysis of a mixed finite element method for a phase field bending elasticity model of vesicle membrane deformation. J Comput Math 2006;24:265–80. [21] Banham T, Li B, Zhao Y. Pattern formation by phase-field relaxation of bending energy with fixed surface area and volume. Phys Rev E 2014;90(3). 033308. [22] Guillén-González F, Tierra G. Unconditionally energy stable numerical schemes for phase-field vesicle membrane model. J Comput Phys 2018;354:67–85. [23] Gao LT, Feng XQ, Gao H. A phase field method for simulating morphological evolution of vesicles in electric fields. J Comput Phys 2009;228:4162–81. [24] Mackay FE, Denniston C. Deformable vesicles interacting in a nematic liquid crystal. Soft Matter 2013;9:5285. [25] Guillén-González F, Rodríguez-Bellido MA, Tierra G. Linear unconditional energy-stable splitting schemes for a phase-field model for nematic-isotropic flows with anchoring effects. Internat J Numer Methods Engrg 2016;108:535–67. [26] Yang X, Forest G, Liu C, Shen J. Shear cell rupture of nematic liquid crystal droplets in viscous fluids. J Non-Newton Fluid Mech 2011;166:487–99. [27] Yang X, Forest G, Li H, Liu C, Shen J, Wang Q, et al. Modeling and simulations of drop pinch-off from liquid crystal filaments and the leaky liquid crystal faucet immersed in viscous fluids. J Comput Phys 2013;236:1–14. [28] Zhao J, Yang X, Shen J, Wang Q. A decoupled energy stable scheme for a hydrodynamic phase-field model of mixtures of nematic liquid crystals and viscous fluids. J Comput Phys 2016;305:539–56. [29] Yue P, Feng JJ, Liu C, Shen J. A diffuse-interface method for simulating two-phase flows of complex fluids. J Fluid Mech 2004;515:293–317. [30] Feng JJ, Liu C, Shen J, Yue P. An energetic variational formulation with phase field methods for interfacial dynamics of complex fluids: advantages and challenges. In: Modeling of soft matter. New York: Springer; 2005, p. 1–26. [31] Tierra G, Guillén-González F. Numerical methods for solving the Cahn-Hilliard equation and its applicability to related energy-based models. Arch Comput Methods Eng 2015;22:269–89. [32] Guillén-González F, Tierra G. On linear schemes for a Cahn Hilliard diffuse interface model. J Comput Phys 2013;234:140–71. [33] Guillén-González F, Tierra G. Second order schemes and time-step adaptivity for Allen-Cahn and Cahn-Hilliard models. Comput Math Appl 2014;68:821–46. [34] Hecht F. New development in FreeFem++. J Numer Math 2012;20:251–65.