A Nitsche-type variational formulation for the shape deformation of a single component vesicle

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. xxx (xxxx) xxx www.elsevier.com/locate/cma

A Nitsche-type variational formulation for the shape deformation of a single component vesicle Tae-Yeon Kima ,∗, Wen Jiangb , Sungmun Leec , Jeong-Hoon Songd , Chan Yeob Yeune , Eun-Jae Parkf a

Civil Infrastructure and Environmental Engineering, Khalifa University of Science and Technology, Abu Dhabi, 127788, United Arab Emirates b Computational Mechanics and Materials, Idaho National Laboratory, Idaho Falls ID 83415, USA c Department of Biomedical Engineering, Khalifa University of Science and Technology, Abu Dhabi 127788, United Arab Emirates d Department of Civil, Environmental and Architectural Engineering, University of Colorado at Boulder, Boulder, CO 80309, USA e Center for Cyber-Physical Systems, Department of Electrical Engineering and Computer Science, Khalifa University of Science and Technology, Abu Dhabi 127788, United Arab Emirates f Department of Computational Science and Engineering, Yonsei University, Seoul 03722, Korea Received 28 November 2018; received in revised form 17 September 2019; accepted 29 September 2019 Available online xxxx

Abstract This paper concerns the development of a finite-element formulation using Nitsche’ method for the phase-field model to capture an equilibrium shape of a single component vesicle. The phase-field model derived from the minimization of the curvature energy results in a nonlinear fourth-order partial differential equation. A standard conforming Galerkin formulation thus requires C 1 -elements. We derive a nonconforming finite-element formulation that can be applied to C 0 -elements and prove its consistency. Continuity of the first derivatives across interelement boundaries is weakly imposed and stabilization of the method is achieved via Nitsche’s method. The capability of the proposed finite-element formulation is demonstrated through numerical study of the equilibrium shapes of axisymmetric single component vesicles along with budding and fission phenomena. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Continuous–discontinuous Galerkin; Giant unilamellar vesicle; Discontinuous Galerkin; Nonconforming

1. Introduction Vesicles are small sacs or cavities containing liquid. In cell biology, vesicles consist of liquid enclosed by a lipid bilayer and they are involved in many biological functions such as transportation of material throughout the cell, metabolism, regulating protein activity, and replication of viruses [1]. The lipid bilayer with the thickness of a few nanometers can form vesicles ranging from 50 nm to tens of micrometers in diameter [2]. Owing to high flexible nature of vesicles, they show a variety of shapes of different symmetry and topology. Also, morphological transitions such as budding and fission can be induced by temperature changes due to membrane fluidity. The ∗ Corresponding author.

E-mail address: [email protected] (T.-Y. Kim). https://doi.org/10.1016/j.cma.2019.112661 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

2

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

equilibrium vesicle shapes arise from the competition between curvature energy, which derives from the bending elasticity of the membrane, geometrical constraints such as fixed surface area and fixed enclosed volume [3]. The objective of this study is to develop a nonconforming finite-element formulation of a phase-field bending elastic model to capture the equilibrium shape of a single component vesicle. Capturing the stationary shape of vesicles can be achieved from minimization of Canham–Helfrich–Evans curvature energy [4–6]. In this study, we employ the phase-field based curvature energy that has been commonly used to derive the phase-field model for modeling the shape dynamics of vesicles [7–10]. One main advantage of the phase-field model is that the complex topological and geometrical changes can be easily captured without explicitly tracking the interfaces. The phase-field model arising from the phase-field based curvature energy gives rise to a nonlinear fourth-order partial differential equation. As a consequence, with a standard conforming Galerkin method, it requires C 1 -basis functions such that the phase-field function and its first derivatives are continuous across element boundaries. A traditional example of C 1 -basis functions is Hermite elements that were developed for two- and threedimensional problems [11–13]. While achieving C 1 -continuity is feasible for one-dimensional problems, it becomes complicated for two- or three-dimensional problems. They can also be relatively quite expensive due to additional derivative unknowns per element, for example, 16 unknowns for the 2D Hermite cubic rectangular element [11]. Mixed-finite element methods can be alternative. However, their formulations involve separate approximations of both primary and secondary fields, resulting in the computational inefficiency and stability issue arising from the combination of interpolation functions for different fields [14]. Several methods based on a phase-field approach have been recently developed for modeling the shape dynamics of vesicles. Du et al. [7] employed the finite-difference method to study the shape deformation of axisymmetric cases of single component vesicle shapes. A three-dimensional modeling of a vesicle was performed using the Fourier spectral method [8]. Du and Zhang [10] developed an adaptive mixed finite-element method for threedimensional vesicle membrane deformations. Most recently, Embar et al. [15] used B-splines to capture the formation and evolution of microdomains on the deforming surface of giant unilamellar vesicles. Two-dimensional and axisymmetric numerical examples of domain evolution coupled to vesicle shape deformation were presented. In this study, we propose a nonconforming finite-element formulation such that C 0 -basis functions can be used to discretize the nonlinear fourth-order phase-field model for the shape deformation of vesicles. Nitsche’s [16] method is employed to weakly impose continuity of the first derivatives of the phase-field variable across interelements along with an additional term for stabilization. Baker [17] first employed Nitsche’s method to develop a non-conforming Galerkin formulation for a fourth-order elliptic problem. The approach in this study is relevant to a consistent C 0 -interior penalty method that was presented for beam and plate theories [18] and a second-gradient theory [19–24]. Recently, application of Nitsche’s method has been further extended for weakly imposing Dirichlet boundary conditions in standard Galerkin methods for second- and fourth-order partial differential equations [25–27] and B-splines [28–35]. The remainder of the paper is organized as follows. In Section 2, we briefly present the governing equation based on the phased-field model for the equilibrium shape of vesicles. In Section 3, we introduce a variational formulation based on Nitsche’s method for the phase-field model along with the proof of its consistency. Notice that this formulation can be applied to any types of C 0 -elements. In Section 4, the capability of the proposed variational formulation is tested for vesicles with several equilibrium shapes. Budding followed by fission is also studied with this formulation. Finally, in Section 5, we provide conclusions and a summary of directions for future work. 2. Phase-field model for the equilibrium shape of vesicles In this section, we provide a governing equation based on a phase-field model for the equilibrium shapes of vesicles. Canham [4], Helfrich [5], and Evans [6] first proposed a purely mechanical model to study the equilibrium shapes of single component vesicles such as red blood cells. The equilibrium shape of a vesicle can be determined by minimizing its shape energy, usually taken to be its bending energy. The Canham–Helfrich–Evans elastic bending energy including the spontaneous curvature is given by ∫ 1 E= κ(H − Hsp )2 ds (2.1) 2 S where κ is the bending rigidity which can depend on the local heterogeneous concentration of the species (such as protein molecules on the blood cells) and H is the mean curvature of the membrane surface (or vesicle) S defined Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

3

by H = (1/R1 + 1/R2 )/2, with R1 and R2 being the principal radii of curvature at any spatial location. Here, Hsp is the spontaneous curvature to represent the natural curvature of the membrane surface, reflecting a possible asymmetry in the bilayer at any spatial location. Some of the possible origins of the spontaneous curvature could be explained by different chemical environments on either sides of the membrane, different chemical compositions of the bilayer, protein insertions, locally varying structure of the lipid molecules, and so on. The original bending energy model consists of minimizing (2.1) with respect to the membrane surface S along with the specified surface area and enclosed volume. In the present work, we use the phase-field model of (2.1) for easy tracking of the shape deformation of the membrane surface. In doing so, the phase-field function φ(x) is defined on the computational domain Ω as described in Du and his colleagues [7,8,10]. In other words, the zero level set, i.e., {x : φ(x) = 0}, gives the membrane surface S, while {x : φ(x) < 0} represents the inside of the membrane and {x : φ(x) > 0} the outside. To derive the phase-field model, we introduce the phase-field bending energy functional [7,8] )2 ( ∫ 1 1 2 κϵ ∆φ − 2 (φ − 1)(φ + Cϵ) dv (2.2) E(φ) = ϵ Ω 2 √ which is equivalent to (2.1). Here, ϵ is the transition parameter related to the interface thickness and C = 2Hsp is the parameter involving the effect of the spontaneous curvature. Notice that (2.2) is defined on the whole domain Ω , not on the membrane surface S, indicating the increase of the dimension from two to three. This gives rise to the increase of the computational cost. Then, the phase-field model for the equilibrium configuration of a vesicle can be derived by minimizing the energy functional (2.2) subject to the prescribed surface area and the prescribed bulk volume. From the work of Du and his colleagues [7,8,10], we define the enclosed volume as (∫ ) 1 V (φ) = φ(x) dv (2.3) 2 Ω for the prescribed volume constraint and the surface area as ] ∫ [ ϵ 1 2 3 2 2 |∇φ| + (φ − 1) dv A(φ) = √ (2.4) 4ϵ 2 2 Ω 2 for the prescribed surface area. To capture the stationary shape of a vesicle, we need to define the initial phase-field function. In this study, we use the following initial phase-field function ) ( d(x, S) (2.5) φ(x) = tanh √ 2ϵ that was commonly used by Du and his colleagues [7,8,10]. Here, d(x, S) is the signed distance function from a point x ∈ Ω to the membrane surface S. Using (2.5), an initial configuration for a bi-concave shape vesicle can be defined as displayed in Fig. 1. The phase-field function φ(x) = −1 for the inside of the membrane (blue) and φ(x) = 1 for its outside (red). The plot clearly shows the transition layer with green color. As a consequence of above statements, as described in [10], the variational phase-field elastic bending model can be stated as ( )2 ∫ 1 2 1 κϵ ∆φ − 2 (φ − 1)(φ + Cϵ) dv (2.6) min E(φ) = φ ϵ Ω 2 subject to constraints A(φ) = A0 ,

V (φ) = V0

(2.7)

and a boundary condition φ = φ0

on ∂Ω ,

(2.8)

where V0 and A0 represent the prescribed volume and the prescribed surface area, respectively. Here, φ0 can be the boundary value. Notice that φ0 can be one with the definition of (2.5), as shown in Fig. 1. Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

4

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 1. Example of an initial phase-field for a bi-concave shape vesicle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

To incorporate surface area and volume constraints into (2.2), the following modified total energy is introduced [10]: αv αa E T (φ) = E(φ) + (V (φ) − V0 )2 + (A(φ) − A0 )2 (2.9) 2 2 in which αv and αa are the penalty parameters for volume and surface area constraints, respectively. The problem to find the equilibrium configuration of a vesicle membrane can be formulated as finding the phase-field function φ(x) on the whole domain Ω that minimizes the energy functional (2.9). Notice that the existence of the minimizer to E P (φ) has been established in [9]. To find a minimizer of (2.9), we use a gradient flow approach that has been successfully used for solving the phase-field model of single component vesicles [7,8,10]. The variational problem can then be solved via the gradient flow δ E T (φ) ∂φ =− (2.10) ∂t δφ where δ/δφ denotes the first variation of the energy functional (2.9). The meaning of (2.10) is that, given an initial guess, as t → ∞, the dynamic solution φ(x, t) converges to a steady state which is a critical point of the energy E T (φ). Taking the first variation of (2.9) and substituting it into (2.10) give rise to the nonlinear fourth order time-dependent partial differential equation ∂φ 1 = ∆Q − 2 Q(3φ 2 + 2Cϵφ − 1) + αv (V (φ) − V0 ) + ∂t ϵ [ ( )] 1 3 αa (A(φ) − A0 ) − √ ϵ ∆φ − 2 (φ 2 − 1)φ ϵ 2 2

(2.11)

where Q = ϵ∆φ − 1ϵ (φ 2 − 1)(φ + Cϵ). Notice that the phase-field model (2.11) with initialization (2.5) is used to find the equilibrium shape of a vesicle. Solving (2.11) requires a nonlinear solver and a time-stepping algorithm. To construct the weak form of (2.11), we introduce the space of an admissible solution field as V ⊂ H 2 (Ω ), where H 2 (Ω ) denotes the Sobolev space of order 2. For convenience, we define the space of a test field w as W ⊂ H 2 (Ω ) with W = {w ∈ H 2 (Ω ) | w = ∇w · n = 0

on ∂Ω }.

(2.12)

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

5

By multiplying a test function w ∈ W, taking integration over Ω , and applying integration by parts, the weak form of (2.11) can be stated as: Find φ ∈ V such that, for all w ∈ W, ∫ ∫ ∫ ∂φ 2 w (∇w · ∇φ)(3φ 2 + 2Cϵφ − 1) dv dv = ϵ ∆w∆φ dv + ∂t ϵ Ω Ω Ω ∫ 1 + w(6φ + 2Cϵ)(∇φ · ∇φ) dv ϵ Ω ∫ 1 w(3φ 2 + 2Cϵφ − 1)(φ 2 − 1)(φ + Cϵ) dv + 3 (2.13) ϵ Ω [ ( )] ∫ 3 1 + αa w(A(φ) − A0 ) − √ ϵ ∆φ − 2 (φ 2 − 1)φ dv ϵ 2 2 Ω ∫ + αv w(V (φ) − V0 ) dv. Ω

Notice that (2.13) is valid for any C 1 -basis functions such as Hermite elements and B-splines. 3. Discretization In this section, we introduce a Nitsche-type variational formulation for the gradient flow (2.11). Since (2.11) is the fourth-order nonlinear partial differential equation, a standard finite-element approximation requires the solution space to be globally H 2 -conforming as described in (2.13). In other words, basis functions are required to be C 1 -continuous, thereby ensuring the continuity of the phase-field function and its first derivatives. In the present study, we derive the non-conforming formulation based on the idea of Engel et al. [18]. This formulation uses C 0 -basis functions so that their first derivatives are discontinuous. Continuity of these first derivatives on element boundaries is weakly imposed based on Nitsche’s method [16]. Notice that the number of unknowns per element arising for this method is considerably fewer than for alternatives based on traditional strategies such as C 1 -basis functions [18]. N To construct the bases, we consider a regular finite-element partition Ω h = ∪e=1 Ωe , with Ω h ≈ Ω and N the total number of elements in the mesh. We choose an approximation function that is continuous on the entire domain, but discontinuous in its first and higher-order derivatives across element boundaries, i.e., we define the finite-dimensional solution and test function spaces as V h = {φh ∈ H 1 (Ω ) | φh |Ωe ∈ Pk (Ωe ), 1 ≤ e ≤ N },

(3.1)

W h = {wh ∈ H 1 (Ω ) | wh |Ωe ∈ Pk (Ωe ), 1 ≤ e ≤ N , wh = ∇wh · n = 0 on ∂Ω },

(3.2)

in which Pk (Ωe ) is the space of complete polynomials with the order less than or equal to k defined over an element Ωe . The approximation φh (x) to the phase-field function φ(x) is given by ∑ φh (x) = N I (ξ (x))φ I (3.3) I

where φ I is the nodal value at node I and ξ (x) is the position in a reference element. Here, N I is the basis function that belongs to the set of Lagrangian iso-parametric functions, i.e., {N I } = {N I ∈ C 0 (Ω h ) : N I |Ωe ∈ Pk (Ωe )},

I = 1, . . . , M

(3.4)

in which M denotes the number of nodes in the mesh. For the sake of brevity, we denote the set of element interiors ˜ and the set of all interior edges by Γ ˜. In two dimensions, Γ ˜ refers only to those element edges that are shared by Ω by two spatially adjacent elements, and does not include edges along the physical boundary ∂Ω . Then, we define ˜ as the integration over Ω ∫ ∑ ∫ φ dv = φ dv ˜ Ω

˜ Ωe ∈ Ω

Ωe

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

6

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

˜ as and the integration over Γ ∫ ∑∫ φ ds φ ds = ˜ Γ

˜ Γi ∈Γ

Γi

with Γi being the element edge. For a scalar function f , the jump operator of f across the interior boundary is defined as JfK = f+ − f−

where f ± = limδ→0 f (x ∓ δn) for δ > 0 with n being any of the two unit normals to the interior boundary. The average operator of f is defined as 1 + ( f + f − ). 2 From the definitions of the jump and average operators, we have the useful identity ⟨⟨ f ⟩⟩ =

(3.5)

J f gK = J f K⟨⟨g⟩⟩ + ⟨⟨ f ⟩⟩JgK.

The Nitsche-type finite-element formulation to approximate the solution of (2.11) can be stated as: Find φh ∈ V h such that, for all wh ∈ W h , ∫ ∫ ∫ ∂φh 2 wh dv = ϵ ∆wh ∆φh dv + (∇wh · ∇φh )(3φh2 + 2Cϵφh − 1) dv ∂t ϵ ˜ ˜ Ω Ω Ω ∫ 1 + wh (6φh + 2Cϵ)(∇φh · ∇φh ) dv ϵ Ω ˜ ∫ 1 + 3 wh (3φh2 + 2Cϵφh − 1)(φh2 − 1)(φh + Cϵ) dv ϵ Ω ˜ ∫ ∫ − ϵ J∇wh · nK⟨⟨∆φh ⟩⟩ ds − ϵ J∇φh · nK⟨⟨∆wh ⟩⟩ ds (3.6) ˜ Γ

˜ Γ

+τ

∫ ˜ Γ

J∇wh · nKJ∇φh · nK ds [

)] ( 3 1 2 wh (A(φh ) − A0 ) − √ ϵ ∆φh − 2 (φh − 1)φh dv ϵ ˜ 2 2 Ω

+ αa

∫

+ αv

∫

wh (V (φh ) − V0 ) dv

˜ Ω

in which τ is the stabilization parameter which is inversely proportional to the mesh size h, i.e., τ ∼ 1/ h. The particular difference of (3.6) with (2.13) is additional terms in the fourth and fifth lines included based on the idea of Nitsche’s method. They aim to weakly enforce continuity of the first derivatives to the normal direction at element interior edges, allowing for the use of C 0 -basis functions. In the following sections, we use φ and w instead of φh and wh in (3.6) for convenience. 3.1. Proof of consistency To prove consistency of (3.6), taking integration by parts of the first term yields ∫ ∫ ∫ ∫ ϵ ∆w∆φ dv = − ϵ ∇w∇(∆φ) dv + ϵ (∇w · n)∆φ ds + ϵ J(∇w · n)∆φK ds ˜ ˜ ˜ Ω ∂Ω Γ ∫ Ω ∫ ∫ 2 =ϵ w∆ φ dv − ϵ w∇(∆φ) · n ds + ϵ (∇w · n)∆φ ds ˜ Ω ∂Ω ∂Ω ∫ ∫ + ϵ J(∇w · n)∆φK ds − ϵ Jw∇(∆φ) · nK ds. ˜ Γ

(3.7)

˜ Γ

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Upon using the identities ∫ ∫ J(∇w · n)∆φK ds = (J∇w · nK⟨⟨∆φ⟩⟩ + ⟨⟨∇w · n⟩⟩J∆φK) ds, ˜ Γ

(3.8)

˜ Γ

∫ ˜ Γ

7

∫ Jw∇(∆φ) · nK ds =

(JwK⟨⟨∇(∆φ) · n⟩⟩ + ⟨⟨w⟩⟩J∇(∆φ) · nK) ds,

(3.9)

˜ Γ

˜, (3.7) reduces to and w = ∇w · n = 0 on ∂Ω and JwK = 0 on Γ ∫ ∫ ∫ ϵ ∆w∆φ dv = ϵ w∆2 φ dv + ϵ (J∇w · nK⟨⟨∆φ⟩⟩ + ⟨⟨∇w · n⟩⟩J∆φK) ds ˜ Ω

˜ Ω

−ϵ

˜ Γ

∫

(3.10)

⟨⟨w⟩⟩J∇(∆φ) · nK ds. ˜ Γ

For convenience, let M(φ) = 3φ 2 + 2Cϵφ − 1. Similarly, applying integration by parts to the second term in (3.6) results in ∫ ∫ ∫ 1 1 1 (∇w · ∇φ)M(φ) dv = − w∇ · [(∇φ)M(φ)] dv + w(∇φ · n)M(φ) ds ϵ Ω ϵ Ω ϵ ∂Ω ˜ ˜ ∫ 1 Jw(∇φ · n)M(φ)K ds + ϵ Γ˜ ∫ ∫ 1 1 =− w∆φ M(φ) dv − w(∇φ) · ∇ M(φ) dv ϵ Ω ϵ Ω ˜ ˜ ∫ ∫ 1 1 + Jw(∇φ · n)M(φ)K ds w(∇φ · n)M(φ) ds + (3.11) ϵ ∂Ω ϵ Γ˜ ˜ ∫ ∫ 1 1 w∆φ M(φ) dv − w(6φ + 2Cϵ)(∇φ · ∇φ) dv =− ϵ Ω ϵ Ω ˜ ˜ ∫ ∫ 1 1 + Jw(∇φ · n)M(φ)K ds w(∇φ · n)M(φ) ds + ϵ ∂Ω ϵ Γ˜ ˜ ∫ ∫ 1 1 3 2 w(∆φ − ∆φ + Cϵ∆φ ) dv + ⟨⟨wM(φ)⟩⟩J(∇φ · n)K ds. =− ϵ Ω ϵ Γ˜ ˜ In the last equality, we used identities ∆φ 3 = 6φ(∇φ · ∇φ) + 3φ 2 ∆φ,

∆φ 2 = 2(∇φ · ∇φ + φ∆φ)

˜. Using the sum of the second and third terms, we can obtain and the continuity of w and φ on Γ ∫ ∫ 1 1 (∇w · ∇φ)M(φ) dv + w(6φ + 2Cϵ)(∇φ · ∇φ) dv ϵ Ω ϵ Ω ˜ ˜ ∫ ∫ 1 1 ∂ M(φ) = (∇w · ∇φ)M(φ) dv + w ∇φ · ∇φ dv ϵ Ω ϵ Ω ∂φ ˜ ˜ ∫ 1 = ∇(wM(φ)) · ∇φ dv ϵ Ω ˜ ∫ ∫ 1 1 =− wM(φ)∆φ dv + wM(φ)(∇φ · n) ds ϵ Ω ϵ ∂Ω ˜ ˜ ∫ 1 + Jw(∇φ · n)M(φ)K ds ϵ Γ˜ ∫ ∫ 1 1 =− 2 wM(φ)(ϵ∆φ) dv + ⟨⟨wM(φ)⟩⟩J(∇φ · n)K ds. ϵ Ω ϵ Γ˜ ˜

(3.12)

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

8

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

On substituting (3.10), (3.11), and (3.12) into (3.6), we obtain ∫ ∫ ∫ ∂φ 1 2 w w(∆φ 3 − ∆φ + Cϵ∆φ 2 ) dv dv = ϵ w∆ φ dv − ∂t ϵ Ω ˜ ˜ ˜ Ω Ω ∫ ∫ 1 1 wM(φ)(ϵ∆φ) dv + 3 wM(φ)(φ 2 − 1)(φ + Cϵ) dv − 2 ϵ Ω ϵ Ω ˜ ˜ ∫ ∫ 2 + ⟨⟨w⟩⟩⟨⟨M(φ)⟩⟩J(∇φ · n)K ds + ϵ (J∇w · nK⟨⟨∆φ⟩⟩ + ⟨⟨∇w · n⟩⟩J∆φK) ds ϵ ∫ Γ˜ ˜ Γ ∫ ∫ −ϵ

J∇w · nK⟨⟨∆φ⟩⟩ ds − ϵ J∇φ · nK⟨⟨∆w⟩⟩ ds ˜ ˜ Γ Γ ∫Γ˜ ∫ + τ J∇w · nKJ∇φ · nK ds + αv w(Vi (φ) − V0 ) dv ˜ ˜ [ (Ω )] ∫Γ 3 1 2 + αa w(A(φ) − A0 ) − √ ϵ ∆φ − 2 (φ − 1)φ dv. ϵ ˜ 2 2 Ω Then, (3.13) can be rewritten as ∫ ∫ ∫ ∂φ 1 w 0= − w∆Q dv − 2 dv + ϵ wM(φ)Q dv ∂t ϵ Ω ˜ ˜ ˜ Ω Ω [ ( )] ∫ ∫ 3 1 + αv w(Vi (φ) − V0 ) dv + αa w(A(φ) − A0 ) − √ ϵ ∆φ − 2 (φ 2 − 1)φ dv ϵ ˜ ˜ 2 2 Ω Ω ) ( ∫ 2 ⟨⟨w⟩⟩⟨⟨M(φ)⟩⟩ − ϵ⟨⟨∆w⟩⟩ + τ J∇w · nK J∇φ · nK ds + ϵ ˜ ∫Γ ∫ + ϵ (⟨⟨∇w · n⟩⟩J∆φK) ds − ϵ ⟨⟨w⟩⟩J∇(∆φ) · nK ds. ˜ Γ

(3.13)

⟨⟨w⟩⟩J∇(∆φ) · nK ds − ϵ

(3.14)

˜ Γ

where 1 ∆Q = ϵ∆2 φ − (∆φ 3 − ∆φ + Cϵ∆φ 2 ). ϵ ˜ and the remaining terms on the third From the first and second lines, we can obtain the gradient flow (2.11) on Ω and fourth lines are jump terms to impose weakly continuities of derivatives on Γ˜ . This proves the consistency of the method. 3.2. Algorithm linearization In this section, we describe the time-stepping algorithm and the nonlinear solver used for the present study. Unless otherwise specified, we use two indices n and i to represent time-step and iteration, respectively. The variational problem (3.6) is nonlinear. A Newton–Raphson iteration scheme is employed to resolve the nonlinear system of equations. The volume and area constraint terms are non-local. As a result, high computational cost is required to compute these terms. For more efficient simulation, we consider these terms as explicit and other terms as implicit, i.e., ∫ φn+1 − φn w dv = Fn+1 + Cn (3.15) ∆t ˜ Ω where ∫ ∫ 2 2 Fn+1 =ϵ ∆w∆φn+1 dv + (∇w · ∇φn+1 )(3φn+1 + 2Cϵφn+1 − 1) dv ϵ Ω ˜ ˜ Ω ∫ 1 + w(6φn+1 + 2Cϵ)(∇φn+1 · ∇φn+1 ) dv ϵ Ω ∫˜ (3.16) 1 2 2 + 3 w(3φn+1 + 2Cϵφn+1 − 1)(φn+1 − 1)(φn+1 + Cϵ) dv ϵ∫ Ω ˜ ∫ −ϵ

˜ Γ

(J∇w · nK⟨⟨∆φn+1 ⟩⟩ + J∇φn+1 · nK⟨⟨∆w⟩⟩) ds + τ

˜ Γ

J∇w · nKJ∇φn+1 · nK ds

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

9

and Cn = αv

∫

w(Vi (φn ) − V0 ) dv [ ( )] ∫ 3 1 + αa w(A(φn ) − A0 ) − √ ϵ ∆φn − 2 (φn2 − 1)φn dv. (3.17) ϵ ˜ 2 2 Ω The variational problem (3.15) can be restated as a nonlinear finite-element problem in terms of the residual R as ∫ φn+1 − φn R(φn+1 , w) = w dv − Cn − Fn+1 = 0, ∀w ∈ W h . (3.18) ∆t ˜ Ω For a given solution φn at the ith iteration, the residual (3.18) is linearized as ⏐ ∂R ⏐⏐ i+1 i δφ, ∀w ∈ W h (3.19) 0 = R(φn+1 , w) ≃ R(φn+1 , w) + ∂φ ⏐ ˜ Ω

n+1 i

for increment δφ. The increment found in (3.19) is used to update solution φn+1 at time-step n + 1: i+1 i φn+1 = φn+1 + δφ.

(3.20)

By taking the variational derivative of R, we can obtain the second term in (3.19) as ∫ ∂R δφn+1 ∂Fn+1 δφ = w dv − ∂φn+1 ∆t ∂φn+1 ˜ Ω where ∫ ∫ 2 ∂Fn+1 2 ∆w∆δφn+1 dv + (∇w · ∇δφn+1 )(3φn+1 + 2Cϵφn+1 − 1) dv =ϵ ∂φn+1 ϵ Ω ˜ ˜ Ω ∫ 2 + (∇w · ∇φn+1 )(6φn+1 + 2Cϵ)δφn+1 dv ϵ Ω ∫˜ ∫ 1 1 + 6wδφn+1 (∇φn+1 · ∇φn+1 ) dv + 2w(6φn+1 + 2Cϵ)(∇φn+1 · ∇δφn+1 ) dv ϵ Ω ϵ Ω ˜ ∫˜ 1 2 + 3 wδφn+1 (6φn+1 + 2Cϵ)(φn+1 − 1)(φn+1 + Cϵ) dv ϵ Ω ˜ ∫ 1 2 wδφn+1 (3φn+1 + 2Cϵφn+1 − 1)(2φn+1 )(φn+1 + Cϵ) dv + 3 ϵ Ω ˜ ∫ 1 2 2 + 3 wδφn+1 (3φn+1 + 2Cϵφn+1 − 1)(φn+1 − 1)dv ϵ∫ Ω ˜ ∫ −ϵ

˜ Γ

(J∇w · nK⟨⟨∆δφn+1 ⟩⟩ + J∇δφn+1 · nK⟨⟨∆w⟩⟩) ds + τ

˜ Γ

(3.21)

(3.22)

J∇w · nKJ∇δφn+1 · nK ds.

Upon considering the finite-element approximation in (3.3) and substituting (3.21) into (3.19), a linear system of equations Kd = −R

(3.23)

is obtained, where the tangent stiffness matrix K originated from (3.21) and the residual vector R from (3.18). Here, d represents the increment of the phase-field function φ. We obtain the increment δφ using (3.23) and update the phase-field solution (3.20) until it is convergent. 4. Numerical study The finite-element formulation (3.6) is implemented in Multi-Physics Object Oriented Simulation Environment (MOOSE), which is an open source finite-element tool to solve partial differential equations [36]. For the implementation of (3.6), we use DG kernel with the standard Lagrangian shape functions and a Newton–Raphson iteration as a nonlinear solver as described in the previous section. The present study is restricted to the axisymmetric equilibrium shapes of the single component vesicles and thus the three-dimensional problem is simplified to the two-dimensional one. Unless otherwise specified, all simulations are performed using the thickness of the transition layer ϵ = 0.2 and the time step ∆t = 10−5 . In MOOSE, QUAD8 and QUAD9 elements can be used as C 0 -basis functions. Thus, the capability of our algorithm is tested for these two elements. Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

10

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 2. Influence of the stabilization parameter τ : zero iso-contours of equilibrium shapes with (a) QUAD8 for Cτ = 100 (red) and Cτ = 0.1, 1, 10, 50 (blue) and (b) QUAD9, Cτ = 1, 10, 50 (blue). All iso-contours are almost indistinguishable. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.1. Influence of parameters In this section, we examine the sensitivity of the numerical solutions on the values of various parameters as shown in Du et al. [7]. Of particular interest are the stabilization parameter τ and the penalty parameters αv and αa for volume and area constraints involved in the variational form (3.6). Moreover, the mesh size h must be chosen based on ϵ to ensure enough resolution in the transition layer. We study the effect of these parameters on the performance of the proposed method using a bi-concave shape vesicle. The equilibrium shape of the bi-concave shape is obtained from the initial shape shown in Fig. 1. We first begin by examining the influence of the stabilization parameter τ . Notice that τ is inversely proportional to h, i.e., τ = Cτ / h with a constant Cτ . Hence, our study is performed by varying Cτ . We choose h = ϵ/4 and αv = αa = 104 based on the following study. Various values of Cτ are investigated for both QUAD8 and QUAD9 elements. The equilibrium shapes of the vesicle obtained from our simulations are displayed in Fig. 2. These shapes are almost indistinguishable with the choices of Cτ displayed in the figure. Our numerical study shows that the nonlinear solver is not convergent if Cτ is too small or large. To get convergent results, our numerical tests show that the lower bound of Cτ for QUAD8 and QUAD9 is 0.1 and 1, respectively, and the upper bound of Cτ for QUAD8 and QUAD9 is 100. Insufficient Cτ would cause stability issue. In Fig. 3, we display examples of the unstable solutions with Cτ = 0.01 for QUAD8 and Cτ = 0.1 for QUAD9. Too large Cτ would worse the conditioning of the global stiffness matrix and make a solution converge difficult. The influence of τ on the evolution of the total energy is examined for both QUAD8 and QUAD9 in Fig. 4. Interestingly, the total energy evolution for QUAD9 is not sensitive to the stabilization parameter. On the other hand, QUAD8 with Cτ = 0.1 is not consistent with other cases. These results indicate that the consistent equilibrium shape of a vesicle can be obtained by choosing 1 ≤ Cτ ≤ 100, i.e., 1/ h ≤ τ ≤ 100/ h for both QUAD8 and QUAD9 elements. In Fig. 5, we examine the influence of the mesh size h within the transition layer ϵ on the equilibrium shape of the bio-concave vesicle structure. We choose Cτ = 1 which is proven to be stable in the previous sensitivity study on the stabilization parameter. Zero iso-contours (i.e., equilibrium shapes) are displayed for various cases of the mesh size (h = ϵ/2, ϵ/4, ϵ/6, ϵ/8) with ϵ = 0.2. For both QUAD8 and QUAD9, while contours obtained using h = ϵ/4, ϵ/6, ϵ/8 are almost indistinguishable, significant deviation of the contour with h = ϵ/2 is observed. This result indicates that, with the choice of h ≤ ϵ/4, we can avoid high sensitivity of the stabilization parameter and ensure enough resolution to resolve the transition layer. In Table 1, the effect of the penalty parameters αv and αa on the change of the relative volume and area is examined. Based on the previous study, we fix h = ϵ/4 and τ = 1/ h. As expected, the relative errors for volume Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

11

Fig. 3. Examples of unstable solutions with insufficient choices of Cτ .

Fig. 4. The evolution of the total energy with various choices of τ .

Table 1 Relative volume and area errors with various values of αv and αa . αv , αa

1000

2500

5000

7500

10,000

QUAD8

(V − V0 )/V0 (A − A0 )/A0

6.1661% 0.4359%

2.3996% 0.1645%

1.1914% 0.0808%

0.7926% 0.0536%

0.5935% 0.0402%

QUAD9

(V − V0 )/V0 (A − A0 )/A0

6.1697% 0.4359%

2.4000% 0.1644%

1.1915% 0.0808%

0.7926% 0.0535%

0.5937% 0.0400%

and area decrease for both QUAD8 and QUAD9 with increasing αv and αa . With αv = αa = 104 , the relative errors for volume and area are approximately 0.6% and 0.04% for both QUAD8 and QUAD9, respectively. 4.2. Shape equilibrium of vesicle structures In this section, the capability of the variational form (3.6) is tested for single component vesicles. From the previous parameter studies, unless otherwise specified, all the simulations are performed using τ = 1/ h, h = ϵ/4, Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

12

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 5. Influence of the mesh size h within the transition layer ϵ for fixed Cτ = 1.

Fig. 6. A three-dimensional view of the shape deformation of a bi-concave vesicle structure obtained from an axisymmetric simulation in two dimensions: (a) initial shape and (b) equilibrium shape.

αv = 104 , and αa = 104 . Notice that other choices of these parameters are possible based on the previous studies. All the simulations are performed using QUAD8 element. 4.2.1. Bi-concave and peanut shapes Our method is first applied to the vesicle with a bi-concave equilibrium shape. The initial phase field used for this study is displayed in Fig. 1, which is defined using (2.5) with ϵ = 0.2. For a better visual effect, we show a three-dimensional view in Fig. 6. The initial configuration of an oblate spheroid is displayed in panel (a) of the figure. The spontaneous curvature is not allowed by setting C = 0 in (3.6). The bi-concave shape vesicle obtained from the simulation is presented in Fig. 6(b). In biological term, this shape may be referred to as a discocyte (or erythrocyte) which is the usual and favored human red blood cell shape. We display the evolution of the total curvature energy in Fig. 7. As shown in the previous section, the curvature energy decreases until it reaches to the steady state, corresponding to the equilibrium shape of the bi-concave disc. Next, our algorithm is applied to obtain a vesicle structure with a peanut or longan shape. Same with the previous example, the zero spontaneous curvature is assumed by setting C = 0. In Fig. 8, a three-dimensional view is presented for better visualization. Fig. 8(b) depicts the ‘peanut’ stationary shape that evolves from an initial prolate spheroid in Fig. 8(a). This peanut shape can be obtained for relatively higher volume-to-surface ratio than the Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

13

Fig. 7. The evolution of total energy (2.9) for a bi-concave shape.

Fig. 8. A three-dimensional view of the shape deformation of a peanut vesicle structure obtained from an axisymmetric simulation in two dimensions: (a) initial shape and (b) equilibrium shape.

bi-concave discocyte (see Figure 14 in Seifert [3] for more details). Although we do not display the evolution of the total energy curvature energy, we notice that it is similar to the previous one in Fig. 7. 4.2.2. Budding and fission In this section, the proposed algorithm is applied to budding and fission transitions of a vesicle. The budding transition is the most prominent example of a vesicle shape transformation with increasing temperature as shown in Figures 1 and 2 of Seifert [3]. Unlike the previous example, by setting positive nonzero values of C in (3.6), the spontaneous curvature effect is included. In Fig. 9, a single case of vesicle budding accompanied by fission is modeled by varying a non-dimensional spontaneous curvature C = 5, 10, 15, 20, 25. We choose a prolate shape as an initial configuration based on Seifert [3]. Budding transition is observed by the choice of C = 5 as shown in Fig. 9(b). With other choices of C, multiple fission transitions are observed in Fig. 9(c)–(f). These results show that our algorithm can model budding and fission phenomena with a non-zero spontaneous curvature. The evolutions of the total energy for these various spontaneous curvatures are displayed in Fig. 10. As expected from (2.2) and (2.9), the total energy increases with increasing C and reaches to the steady state for all the values of C chosen. Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

14

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

Fig. 9. Equilibrium shapes for different spontaneous curvatures.

Fig. 10. The evolution of curvature energy for different spontaneous curvatures.

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

15

Table 2 Comparison of the computational time with the Hermite cubic element. Element type

# of dofs

Total time (s)

Ki (s)

Ke (s)

Solver (s)

HERMITE QUAD8 QUAD9

161,604 120,801 160,801

74.895 58.075 91.086

39.226 5.555 7.171

0 2.467 3.048

25.29 39.092 67.165

4.3. Computational efficiency In this section, the computational cost of the present algorithm is compared with the Hermite cubic element which is available on MOOSE. Notice that the Hermite cubic element has 16 degrees of freedom (dofs) per element due to additional derivative unknowns at each node. On the other hand, QUAD8 and QUAD9 have only 8 and 9 dofs per element, respectively. All calculations for this comparison are carried out in parallel with 4 MPI processes on a workstation Xeon E5-2637 v3 CPU. A directLU solver is employed to obtain the solutions of the linear algebraic system of equations for each nonlinear step. For accurate numerical integration, Gaussian quadrature points we used are 3 × 3 for QUAD8 and QUAD9 and 4 × 4 for the Hermite cubic element. For all three element types, we choose the same number of elements with the relation h = ϵ/4 between the mesh size and the transition layer thickness. In Table 2, we display the computational time for only one time step which takes five nonlinear iterations. For convenience, we divide the tangent stiffness matrix K into two parts: Ki and Ke . In (3.21), Ki and Ke represent ˜ and edges Γ ˜, respectively. Notice that the standard the numerical integration of the terms over element interiors Ω Galerkin formulation using C 1 -elements does not need Ke . This is the additional computational cost of the proposed algorithm in this study. For better understanding, we also display the computational time of assembling Ki and Ke and the solver. The total computational time of QUAD8 is relatively smaller than others, indicating the computational efficiency of our algorithm. Assembling Ki for Hermite takes more time than others due to the large number of dofs and Gaussian quadrature points. Interestingly, while the computational time to assemble Ke for QUAD8 and QUAD9 is not significant, solving the system of equations using the directLU solver takes more time than Hermite. We attribute this to the increase of the bandwidth of the stiffness matrix K in QUAD8 and QUAD9 due to the presence of Ke . This result shows that the proposed C 0 -element based algorithm can improve the computational efficiency. We expect that using a more efficient solver can reduce the computational cost of our algorithm. Moreover, further numerical studies for the error behaviors and convergence rate of the proposed method remain as our future work. 5. Summary and conclusion In this paper, we presented a nonconforming finite-element method for modeling the shape deformation of vesicles using the phase-field model. The governing equation arising from the phase-field curvature energy is nonlinear and fourth-order with respect to the phase-field variable. Rather than employing C 1 -basis functions or a fully mixed approach, we base our method on a continuous–discontinuous Galerkin formulation using C 0 -basis functions. Continuity of higher-order derivatives of the phase-field variable is then enforced between elements using an idea of Nitsche’s method, involving jump quantities across interelement boundaries. Using our finite-element formulation, we captured the equilibrium shapes of a bi-concave and peanut type vesicles. Moreover, budding and fission phenomenon of a vesicle was modeled by considering the non-zero spontaneous curvature. Numerical solutions are sensitive to the parameters involving in the formulation. As a result, we studied the influence of these parameters including the mesh size to ensure enough resolution on the transition layer. The results show that the stabilization parameter must be chosen 1/ h ≤ τ ≤ 100/ h for both QUAD8 and QUAD9 to get convergence solutions for the Newton–Raphson iteration. Moreover, choosing the mesh size h ≤ ϵ/4 is less sensitive to the stabilization parameter τ . As expected, the relative volume and area errors decrease with increasing the penalty parameters for volume and area constraint. The phase-field model is defined on the whole three-dimensional physical domain, not on the surface of the vesicle in two dimensions. Owing to this increase of the dimension, one of main challenges using the phase-field model is high computational cost. In addition, the phase-field function rapidly changes near the transition layer Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

16

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

around the vesicle membrane surface. Thus, highly refined mesh is needed for accurate modeling of the transition layer. Therefore, adaptive refinement near the transition layer is more ideal than uniform meshes on the whole domain. Developing an algorithm with easy adaptive refinement is important to model the shape deformation of the vesicle. We expect that application of adaptivity to our proposed C 0 -elements based formulation can be easily achieved. Future work will focus on developing an adaptive technique of the proposed formulation. Moreover, the current study was restricted to the axisymmetric cases of the single component vesicles. In the future, our study will be extended to model the equilibrium shapes of the multi-component vesicles and eventually applied to a full three-dimensional analysis. We will also develop a more efficient algorithm based on the meshfree point collocation method [37–39]. Acknowledgments T.-Y. Kim and S.M. Lee gratefully appreciate the support from Khalifa University, United Arab Emirates Internal Research Funding (KURIF Level 2 No. 8474000002) and ADEK Award for Research Excellence (AARE) 2017, United Arab Emirates (AARE17-069). E.-J. Park was supported by NRF, Korea-2015R1A5A1009350 and NRF, Korea-2019R1A2C2090021. References [1] O. Mouritsen, Life - As a Matter of Fat: The Emerging Science of Lipidomics, Springer, Berlin, 2005. [2] L. Deseri, M.D. Piccioni, G. Zurlo, Derivation of a new free energy for biological membranes, Contin. Mech. Thermodyn. 20 (2008) 255–273. [3] U. Seifert, Configurations of fluid membranes and vesicles, Adv. Phys. 46 (1997) 13–137. [4] P. Canham, The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell, J. Theoret. Biol. 26 (1970) 61–81. [5] W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Z. Naturforsch. C 28 (1973) 693–703. [6] E.A. Evans, Bending resistance and chemically induced moments in membrane bilayers, Biophys. J. 14 (1974) 923–931. [7] Q. Du, C. Liu, X. Wang, A phase field approach in the numerical study of the elastic bending energy for vesicle membranes, J. Comput. Phys. 198 (2004) 450–468. [8] Q. Du, C. Liu, X. Wang, Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions, J. Comput. Phys. 212 (2006) 757–777. [9] Q. Du, X. Wang, Convergence of numerical approximations to a phase field bending elasticity model of membrane deformations, Int. J. Numer. Anal. Model. 4 (2007) 441–459. [10] Q. Du, J. Zhang, Adaptive finite-element method for a phase field bending elasticity model of vesicle membrane deformations, SIAM J. Sci. Comput. 30 (2008) 1634–1657. [11] F.K. Bogner, R.L. Fox, L.A. Schmit, The generation of interelement-compatible stiffness and mass matrices by the use of interpolation formula, in: Proceedings of 1st Conference on Matrix Methods in Structural Mechanics AFFDL-TR-66-80, 1966, pp. 397–443. [12] J. Petera, J.F.T. Pittman, Isoparametric hermite elements, Internat. J. Numer. Methods Engrg. 37 (1994) 3489–3519. [13] S.A. Papanicolopulos, A. Zervos, I. Vardoulakis, A three-dimensional C1 finite element for gradient elasticity, Internat. J. Numer. Methods Engrg. 77 (2009) 1396–1415. [14] M. Fortin, F. Brezzi, Mixed and Hybrid Finite Element Methods, Springer, New York, 1991. [15] A. Embar, J. Dolbow, E. Fried, Microdomain evolution on giant unilamellar vesicles, Biomech. Model. Mechanobiol. 12 (2012) 597–615. [16] J.A. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Semin. Univ. Hambg. 36 (1970/71) 9–15. [17] G.A. Baker, Finite element methods for elliptic equations using nonconforming elements, Math. Comp. 31 (1977) 45–59. [18] G. Engel, K. Garikipati, T.J.R. Hughes, M.G. Larson, L. Mazzei, R.L. Taylor, Continuous/discontinuous finite element approximations of fourth-order elliptic problems in structural and continuum mechanics with applications to thin beams and plates, and strain gradient elasticity, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3669–3750. [19] T.-Y. Kim, J.E. Dolbow, E. Fried, A numerical method for a second-gradient theory of incompressible fluid flow, J. Comput. Phys. 223 (2007) 551–570. [20] T.-Y. Kim, J.E. Dolbow, An edge-bubble stabilized finite element method for fourth-order parabolic problems, Finite Elem. Anal. Des. 45 (2009) 485–494. [21] T.-Y. Kim, J.E. Dolbow, E. Fried, Numerical study of the grain-size dependent Young’s modulus and Poisson’s ratio of bulk nanocrystalline materials, Int. J. Solids Struct. 49 (2012) 3942–3952. [22] T. Singhal, E. Kim, T.-Y. Kim, J. Yang, Weak bond detection in composites using highly nonlinear solitary waves, Smart Mater. Struct. 26 (5) (2017) 055011. [23] A. Schiffer, A.I. Alkhaja, J. Yang, E.N. Esfahani, T.-Y. Kim, Interaction of highly nonlinear solitary waves with elastic solids containing a spherical void, Int. J. Solids Struct. 118–119 (2017) 204–212.

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.

T.-Y. Kim, W. Jiang, S. Lee et al. / Computer Methods in Applied Mechanics and Engineering xxx (xxxx) xxx

17

[24] T.-Y. Kim, L.G. Rebholz, E. Fried, A deconvolution enhancement of the Navier–Stokes-αβ model, J. Comput. Phys. 231 (11) (2012) 4015–4027. [25] T.-Y. Kim, E.-J. Park, D.-W. Shin, A C 0 -discontinuous Galerkin method for the stationary quasi-geostrophic equations of the ocean, Comput. Methods Appl. Mech. Engrg. 300 (2016) 225–244. [26] D. Shin, Y. Kang, E.J. Park, C 0 -Discontinuous Galerkin methods for a wind-driven ocean circulation model: Two-grid algorithm, Comput. Methods Appl. Mech. Engrg. 328 (2018) 321–339. [27] A. Embar, J. Dolbow, I. Harari, Imposing Dirichlet boundary conditions with Nitsche’s method and spline-based finite elements, Internat. J. Numer. Methods Engrg. 83 (2010) 877–898. [28] T.-Y. Kim, E. Puntel, E. Fried, Numerical study of the wrinkling of a stretched thin sheet, Int. J. Solids Struct. 49 (2012) 771–782. [29] T.-Y. Kim, T. Iliescu, E. Fried, B-spline based finite-element method for the stationary quasi-geostrophic equations of the ocean, Comput. Methods Appl. Mech. Engrg. 286 (2015) 168–191. [30] W. Jiang, T.-Y. Kim, Spline-based finite-element method for the stationary quasi-geostrophic equations on arbitrary shaped coastal boundaries, Comput. Methods Appl. Mech. Engrg. 299 (2016) 144–160. [31] W. Jiang, C. Annavarapu, J.E. Dolbow, I. Harari, A robust Nitsche’s formulation for interface problems with spline-based finite elements, Internat. J. Numer. Methods Engrg. 104 (7) (2015) 676–696. [32] W. Jiang, J.E. Dolbow, Adaptive refinement of hierarchical B-spline finite elements with an efficient data transfer algorithm, Internat. J. Numer. Methods Engrg. 102 (2015) 233–256. [33] I. Al Balushi, W. Jiang, G. Tsogtgerel, T.-Y. Kim, Adaptivity of a B-spline based finite-element method for modeling wind-driven ocean circulation, Comput. Methods Appl. Mech. Engrg. 332 (2018) 1–24. [34] D. Kim, T.-Y. Kim, E.J. Park, D. Shin, Error estimates of B-spline based finite-element methods for the stationary quasi-geostrophic equations of the ocean, Comput. Methods Appl. Mech. Engrg. 335 (2018) 255–272. [35] N. Rotundo, T.-Y. Kim, W. Jiang, L. Heltai, E. Fried, Error analysis of a B-spline based finite-element method for modeling wind-driven ocean circulation, J. Sci. Comput. 69 (1) (2016) 430–459. [36] D. Gaston, C. Newman, G. Hansen, D. Lebrun-Grandié, MOOSE: A parallel computational framework for coupled systems of nonlinear equations, Nucl. Eng. Des. 239 (2009) 1768–1778. [37] L.C. Berselli, T.-Y. Kim, L.G. Rebholz, Analysis of a reduced-order approximate deconvolution model and its interpretation as a Navier–Stokes–Voigt regularization, Discrete Contin. Dyn. Syst. Ser. B 21 (4) (2016) 1027–1050. [38] J.-H. Song, Y. Fu, T.-Y. Kim, Y.-C. Yoon, J.G. Michopoulos, T. Rabczuk, Phase field simulations of coupled microstructure solidification problems via the strong form particle difference method, Int. J. Mech. Mater. Des. 14 (4) (2018) 491–509. [39] A. Beel, T.-Y. Kim, W. Jiang, J.-H. Song, Strong form-based meshfree collocation method for wind-driven ocean circulation, Comput. Methods Appl. Mech. Engrg. 351 (2019) 404–421.

Please cite this article as: T.-Y. Kim, W. Jiang, S. Lee et al., A Nitsche-type variational formulation for the shape deformation of a single component vesicle, Computer Methods in Applied Mechanics and Engineering (2019) 112661, https://doi.org/10.1016/j.cma.2019.112661.