ADAPTIVE FINITE ELEMENT METHOD FOR A PHASE FIELD
SIAM J. SCI. COMPUT.Vol. 30, No. 3, pp. 1634–1657c 2008 Society for Industrial and Applied Mathematics ADAPTIVE FINITE ELEMENT METHOD FOR A PHASE FIELDBENDING ELASTICITY MODEL OF VESICLE MEMBRANEDEFORMATIONS QIANG DU† AND JIAN ZHANG†Abstract. In this paper, a three-dimensional adaptive finite element method is developed fora variational phase field bending elasticity model of vesicle membrane deformations. Using a mixedfinite element formulation, residual type a posteriori error estimates are derived for the associatednonlinear system of equations and, they are used to introduce the mesh refinement and coarsening.The resulting mesh adaptivity significantly improves the efficiency of the phase field simulation ofvesicle membranes and enhances its capability in handling complex shape and topological changes.The effectiveness of the adaptive method is further demonstrated through numerical examples.Key words. vesicle membrane, phase field, elastic bending energy, a posteriori error estimator,adaptive finite element, mixed finite elementAMS subject classifications. 65N30, 70G75, 92C05DOI. 10.1137/0606564491. Introduction. This paper presents an adaptive finite element method for thenumerical simulation of vesicle membrane deformation based on a phase field bending elasticity model. The vesicle membranes, formed by bilayers of lipid molecules,are simple forms of biological membranes which exist everywhere in life and compartmentalize living matter into cells and subcellular structures, and are essential formany biological functions [38]. The equilibrium shapes of bilayer vesicle membraneshave been successfully modeled via the minimization of certain shape energy; see, forinstance, [14, 36, 41, 48], and the references cited therein.In the isotropic case, the most relevant energetic contribution to the equilibriummembrane geometry is usually the elastic bending energy of the form [11, 12, 41]: k 2Eelastic H ds,(1.1)Γ 2where H is the mean curvature of the membrane surface. The parameter k is thebending rigidity, which can depend on the local heterogeneous concentration of thespecies (such as protein and cholesterol molecules), but it is mostly assumed to be aconstant in this manuscript. Taking a simplified description of the effect of densitychange and osmotic pressure, it is assumed that the variation of the bending energyis subject to the constraints of specified volume and surface area. More general formsof the bending elastic energy, attributed to Canham and Helfrich, also incorporateeffects of surface tension, the Gaussian curvature, and the spontaneous curvature[36, 41]. For the sake of simplicity in our presentation, we only focus on the energy(1.1), though much of our studies can be naturally extended to more general casesincluding the effect of the spontaneous curvature [20], the Gaussian curvature [21, 24],and the vesicle fluid interactions [17, 18]. Receivedby the editors April 5, 2006; accepted for publication (in revised form) October 31,2007; published electronically April 18, 2008. This research is supported in part by NSF-DMR0205232 and NSF-DMS 44.html† Department of Mathematics, Pennsylvania State University, University Park, PA 16802 (qdu@math.psu.edu, zhang j@math.psu.edu).1634Copyright by SIAM. Unauthorized reproduction of this article is prohibited.
ADAPTIVE FEM FOR PHASE FIELD MODEL OF MEMBRANES1635Computationally, there are many simulation methods developed for studying various deforming interface problems, such as the boundary integral and boundary element methods, immersed boundary and interface methods, front-tracking methods,and level-set methods (see, for instance, [34, 40, 43, 46, 49], and the references givenin [22, 23]). For bending elasticity models, applications of these types of methods canalso be found in [6, 28, 32, 53]. The phase field model can be viewed as a physicallymotivated level-set method; this is by virtue of its energy-based variational formalism.One of the main attractions to use the phase field method is its capability of easilyincorporating the complex morphological changes of the interface, in particular, thechanges in both topological and geometrical structures. For more detailed discussions,we refer to [22, 23], and the references given there.In [22], a finite difference method was used to study the energy-minimizing vesiclemembrane in the three-dimensional axis-symmetric case. In [23], a Fourier spectralmethod was used to study the full three-dimensional case. Parallel implementation ofsuch a spectral approach was also carried out to improve the computational efficiency.The various simulation examples given in the earlier studies demonstrated the effectiveness of the phase field approach. However, in the three-dimensional case, the highcomputational cost remains a formidable challenge in making the phase field simulation efficient. Indeed, the phase field function is defined on the whole physical domain,and it changes rapidly only near the transition layer around the membrane surface(the zero level set of the phase field function). Hence, uniform computational grids aregenerally not optimal, and it is natural to consider the application of adaptive finiteelement methods based on a posteriori error estimators [2, 3, 50]. We anticipate thatadaptivity based on effective error estimations could make the phase field simulationmuch more efficient computationally and yet retain the advantage of being able toavoid the explicit tracking of the interfaces. This is indeed confirmed by the presentwork.In the adaptive method presented in this paper, a mixed finite element method(FEM) is used to discretize the phase field bending elasticity model. A residualtype a posteriori error estimator is derived for the development of the adaptive FEMalgorithm. Effectively, the nodes of the adaptive mesh are concentrated near theinterface (the membrane surface) so that the number of nodes is significantly reducedcompared with the number of nodes in the uniform mesh cases, while the resolutionof the numerical resolution of the adaptive FEM remains at the same level. Thesenumerical results reveal the great potential of using the adaptive FEM to significantlyreduce the computational cost of phase field approaches.Detailed descriptions, analysis, and numerical examples of our adaptive FEMapproach are presented in the rest of the paper as follows: in section 2, a brief introduction is given to the phase field bending elasticity model for the vesicle membraneproblem. In section 3, we set up a finite element discretization for the model basedon a mixed formulation. In section 4, we derive some a posteriori error estimators.In section 5, an adaptive algorithm is outlined along with discussions on other implementation issues involved. In section 6, numerical examples are presented, and finally,in section 7, some concluding remarks are given.2. The phase field bending elasticity model. As in [22, 23], we introduce aphase function φ φ(x) defined on the physical (computational) domain Ω, which isused to label the inside and the outside of the vesicle Γ. The level set {x : φ(x) 0}gives the membrane surface Γ, while {x : φ(x) 0} represents the outside of themembrane and {x : φ(x) 0} the inside. The original elastic bending energy modelCopyright by SIAM. Unauthorized reproduction of this article is prohibited.
1636QIANG DU AND JIAN ZHANGconsists of minimizing (1.1) among all surfaces with specified surface area and enclosedvolume.Define the following modified elastic energy: 2 1k φ 2 (φ2 1)φ dx,E(φ) (2.1) Ω 2where is a transition parameter which is taken to be very small compared to thesize of the vesicle. In our numerical examples, varies in the same range of that in[23], that is, taking values generally about 1/40 to 1/200 of the size of Ω to ensure asufficient level of resolution of the interface and the phase field solution profile.For a minimizer φ of E(φ), it has been shown that 1H(φ) φ 2 (φ2 1)φ is the phase field approximation of the mean curvature H of the interface, whichbehaves like a measure concentrated in the diffuse interfacial layer of the zero levelset of φ as goes to zero [23, 25]. Notice that the definition of E(φ) is 1k 2H (φ)dx,E(φ) Ω2with the normalizing factor 1/ accounting for the contribution due to the diffuseinterfacial layer. The phase field bending elasticity model is then given by minimizingthe above energy E(φ) subject to prescribed values of 1 222 φ (φ 1) dx .(2.2)φ(x)dx and B(φ) A(φ) 4 ΩΩ 2Intuitively, it is insightful to consider a special phase field function of the form φ(x) tanh( d(x,Γ)), where d(x, Γ) is the signed distance from a point x Ω to the2 surface Γ; the geometric meanings of the energy E and the constraints (2.2) wouldthen become clear. More rigorously, the existence of the minimizer to E(φ) subject toprescribed A(φ) and B(φ) has been established in [25]. Moreover, it has been shownin [19, 51] that under some general ansatz assumptions, as 0, the minimum of thephase field energy E(φ) with the specified constraints approaches to the minimum ofthe original energy (1.1), with A(φ) approaching the difference of the outside volume and the inside volume of the membrane surface and B(φ) approaching to 2 2/3times the surface area of Γ. The previously developed finite difference, finite element,and spectral methods in [22, 23, 25, 26] are based on the above phase field bendingelasticity model. Some convergence analysis and a priori error estimates have beengiven in [25, 26].3. The finite element discretization. The variational phase field bendingelasticity model described in the previous section is conveniently stated as: 2 1 2k arg min E(φ) φ 2 (φ 1)φ dx, φ Ω 2 A(φ) φdx α,Ω(3.1) 1 2 22 φ (φ 1) dx β,B(φ) 4 Ω 2 φ Ω 1,Copyright by SIAM. Unauthorized reproduction of this article is prohibited.
ADAPTIVE FEM FOR PHASE FIELD MODEL OF MEMBRANES1637where α and β are given constants, and Ω [ 1, 1]3 . The well-posedness of the aboveproblem can be found in [19, 25]. The consistency to the original bending elasticitymodel when 0 has also been examined in [19, 51].To deal with the nonlinear constraints, a penalty formulation can be used. LetG(φ) E(φ) M1 [A(φ) α]2 M2 [B(φ) β]2 ,where M1 and M2 are penalty constants, and let us introduceλM (φ) 2M1 [A(φ) α] ,μM (φ) 2M2 [B(φ) β] .As the penalty constants M1 and M2 go to infinity, the minimizer of G goes tothe solution of the constrained problem, and λM and μM converge to the Lagrangemultipliers [22, 25]. For the sake of simplicity, we denote them by λ and μ in theanalysis. In the simulations, M1 and M2 are taken to be fixed large numbers thatassure the convergence of the Lagrange multipliers to within the given numericalaccuracy.To derive the mixed weak formulation [26], we introduce a function f asf 1 2k φ 2 (φ 1)φ . We note that f is a scaled phase field approximation of the mean curvature, andfor small other boundary conditions such as the homogeneous Neumann boundarycondition on φ may also be used.Multiplying the equation for f by a test function w H01 (Ω) and integrating overΩ, after integration by parts, we get 1 2 φ · w 2 (φ 1)φw dx Ω f wdx k Ωfor any w H01 (Ω). Note that the boundary condition φ 1 is imposed.Taking the variational derivative of the energy functional, we get 12 f · v 2 f (3φ 1)v dx λ μ ΩΩ k fk vdx 0for any v H 1 (Ω). Given any spatial region D Ω, let (u, v)D 1/2 u D (u, u)Duvdx,Ddenote the standard L2 inner product and the L2 norm on D, respectively, and u, vD u · vdx . DCopyright by SIAM. Unauthorized reproduction of this article is prohibited.
1638QIANG DU AND JIAN ZHANGLet us defineH H 1 (Ω) H 1 (Ω)andH0 H 1 (Ω) H01 (Ω).The weak form of our problem (3.1) is to find (f, φ) with (f, φ 1) H0 suchthat, for all (v, w) H0 , we have(3.2) 1 2 k f, v Ω 2 f (3φ 1), v Ω λ(1, v)Ω (f, v)Ω 0, μ k 1 2 k φ, w Ω 2 (φ 1)φ, w Ω (f, w)Ω 0 . The above weak form leads naturally to a mixed FEM for its numerical solution [26].Define the operator F : H L(H0 , R2 ) byF (f, φ)(v, w) k f, v12Ω 2 f (3φ 1), v k φ, wΩ (f, v)Ω ,k λ(1, v)Ω μW Ω 1(φ2 1)φ, w 2 Ω T (f, w)Ω,where W is a weight constant to be determined in simulations. Now the problembecomes that of solving(3.3)F (f, φ) 0in an abstract form, together with the boundary conditions.To construct the finite element approximation to the mixed formulation, we takethe discrete function space asVh Wh {v C 0 (Ω) H 1 (Ω) v K P1 (K) K Jh },where Jh is a triangulation of Ω consisting of tetrahedra K, whose diameters hK arebounded above by h maxK Jh hK , and P1 (K) denotes the linear function space onelement K. The mesh is assumed to be regular so that the standard minimum anglecondition is satisfied and the number of adjacent elements to any given element isbounded independently of h. We also assume that the family Jh is uniformly regularso that for any K Jh the ratio of the diameter of K and the diameter of the largestball enclosed in K is uniformly bounded from above by a constant independent of h.The detailed adaptive construction of Jh will be described later in the paper. Now,letX h Vh W hand Xh0 Vh0 Wh ,Copyright by SIAM. Unauthorized reproduction of this article is prohibited.
ADAPTIVE FEM FOR PHASE FIELD MODEL OF MEMBRANES1639whereVh0 Wh0 {v v Vh , v Ω 0}.For the trial solution (fh , φh 1) Xh0 satisfying the boundary conditions and thetest function (vh , wh ) Xh0 , we define the discrete version of F :Fh (fh , φh )(vh , wh ) F (fh , φh )(vh , wh ).Hence the discrete problem is to solve(3.4)Fh (fh , φh ) 0.The solutions of (3.3) and (3.4) depend on the parameters α and β. By normalization,we may always take β 1(some constant surface area); then the solution of (3.3) maybe viewed as a solution branch ofF (α; x) 0.The basic mathematical analysis of the discretized system (such as the existence ofsolutions and an a priori error estimate, specialized to a piecewise linear element) hasbeen made in [26]. In particular, if there exists a C 1 branch {α, x(α)} of nonsingularsolutions of the nonlinear variational problem (3.3) for α Λ which is a compactinterval in R, then for small h there is a unique branch {α, xh (α)} of solutions of(3.4) converging to {α, x(α)}. Moreover, if α x(α) is a C 1 function from Λ intoH (H 2 (Ω) H 2 (Ω)), we have the optimal order a priori error estimate x(α) xh (α) H Ch,where C is a constant independent of h.While the a priori error analysis can offer theoretical assurance on the convergenceof the finite element methods as the mesh size gets smaller and smaller, to implementan adaptive strategy for the finite element approximations, we need a posteriori errorestimators. We choose to work with residual-type estimators which are derived in thenext section.4. A posteriori error estimate. Adaptive methods often lead to efficient discretization to problems with solutions that are singular or have large variations insmall scales. In phase field models, the sharp interface of physical quantities is replaced by regularized phase field functions. However, for a small interfacial widthconstant , the phase field solutions may display large gradients within the diffusiveinterfacial region. Thus, adaptivity in the form of mesh refinement and coarseningas well as mesh transformation can greatly improve the efficiency of the numericalapproximations of phase field models [10, 27, 37, 44]. A posteriori error estimatorsare key ingredients in the design of adaptive methods [3]. There have been manyexisting studies on deriving such estimators for the finite element approximation oflinear and nonlinear variational problems and for standard Galerkin and mixed finiteelement formulations; see, for example, [1, 8, 45, 50], and the references cited therein.Let x0 (f0 , φ0 ) be a solution of a nonlinear operator equation F (x0 ) 0. Wecall x0 a regular solution if the Frechet derivative DF (x0 ) is well defined and is alinear homeomorphism; that is, DF (x0 ) and its inverse are bijective and continuouslinear operators. First, using the abstract approximation results from [50] (see alsoCopyright by SIAM. Unauthorized reproduction of this article is prohibited.
1640QIANG DU AND JIAN ZHANG[7], and particular applications to Ginzburg–Landau-type models [16], and the phasefield bending elasticity model of vesicle membranes presented in [26]), we immediatelyget the following.Proposition 4.1 (see [50]). Let H0 L(H0 , R2 ) and F C 1 (H, H0 ). Let x0 bea regular solution of F (x) 0 with Z DF (x0 ) L(H0 ,H 0 ) and Ẑ DF (x0 ) 1 L(H 0 ,H0 ) .Assume, in addition, that DF is Lipschitz continuous at x0 with a constant γ 0,i.e., there is a M0 such thatγ : sup x x0 H M0 DF (x) DF (x0 ) L(H0 ,H 0 ) . x x0 HSetM : min{M0 , γ 1 Ẑ 1 , 2γ 1 Z}.Then, for all x H with x x0 H M .(4.1)1 F (x) H 0 x x0 H 2Ẑ F (x) H 0 .2ZTo apply the above framework to our problem, one relevant issue is how thevarious constants γ, Z, Ẑ, and M depend on . We note that obtaining a very precisedependence in general is very difficult, and it remains largely as an open problemdue to the nonlinear nature of the problem (see [29, 39, 42] for related discussion).Some crude estimates can be obtained, for example, by noticing that for given f0 , φ0the components of the vector-valued operator DF (f0 , φ0 ) contain typical terms like k (Δ 12 (3φ20 1)I). It is then reasonable to take a typical tanh profile for φ0based on the earlier sharp interface limit analysis in [19] and to get estimates on thespectra of such operators. Similarly, one can estimate the Lipschitz constant γ. Theseestimates, however, are not sharp in general, and future studies are obviously neededto examine such dependence more carefully.We now continue to derive a posteriori error estimators. Let xh (fh , φh ) be asolution of Fh (xh ) 0 and Rh be a restriction or projection operator from H to thefinite element space. For y H0 , by Galerkin orthogonalityFh (xh )Rh y 0.HenceF (xh )y Fh (xh )y F (xh )(y Rh y).Specifically, let Rh [Ih , Ih ], where Ih denotes the Clement interpolation [13]. Wefirst note that, for a given element K and a given edge e defined by the triangulationJh , the operator Ih has the following properties: u Ih u K ChK u N (K)and1 u Ih u L2 (e) Che2 u N (e) ,where N (K) denotes the union of K and its neighbor elements and N (e) denotes theunion of elements that have e as a face.Copyright by SIAM. Unauthorized reproduction of this article is prohibited.
1641ADAPTIVE FEM FOR PHASE FIELD MODEL OF MEMBRANESNext, we have F (fh , φh )(v, w) F (fh , φh )(v Ih v, w Ih w) k fh , (v Ih v) KK k f (3φ2h 1) λ μ3/2 h K W 2k φh , (w Ih w) K 2 fh , v Ih vkKK f K k 2(φ 1)φh , w Ih w 3/2 h 2 12 K fh v N (e) Ck he e n e e Ω 2 k 2fh v N (K) hK λ 3/2 fh (3φh 1) μ k KK 1 1 φh w N (e) W 2k he2 e 2 n e e Ω 2 12 k 2 hK fh 3/2 (φh 1)φh w N (K). 1212KηKLet 11 C1 k he2 e 2 e Ω KK fh k hK λ 3/2 fh (3φ2h 1) μ n e 2 fh k K 1 2 2 1 1 φhk 2 C3 hK f(4.2) W 2 C2 k he2 e 2 , (φ 1)φ hhh n e 3/2 Ke Ω K where C1 , C2 , and C3 are weights depending on constants in the interpolation inequalities. Without loss of generality, we simply take them to be 1, and the choice ofthe weight W used in the simulations is described later. We have12 2 F (fh , φh ) L(H0 ,R2 ) CηK.KApplying Proposition 4.1, we get(4.3) (fh , φh ) (f0 , φ0 ) H C 122ηK.KWe would like to point out that the constant C depends on through Ẑ
adaptive finite element, mixed finite element AMS subject classifications. 65N30, 70G75, 92C05 DOI. 10.1137/060656449 1. Introduction. This paper presents an adaptive finite element method for the numerical simulation of vesicle membrane deformation based on a phase field bend-ing elasticity model.
Finite Element Method Partial Differential Equations arise in the mathematical modelling of many engineering problems Analytical solution or exact solution is very complicated Alternative: Numerical Solution – Finite element method, finite difference method, finite volume method, boundary element method, discrete element method, etc. 9
An adaptive mixed least-squares finite element method for . Least-squares Raviart–Thomas Finite element Adaptive mesh refinement Corner singularities 4:1 contraction abstract We present a new least-squares finite element method for the steady Oldroyd type viscoelastic fluids.
Finite element analysis DNV GL AS 1.7 Finite element types All calculation methods described in this class guideline are based on linear finite element analysis of three dimensional structural models. The general types of finite elements to be used in the finite element analysis are given in Table 2. Table 2 Types of finite element Type of .
1 Overview of Finite Element Method 3 1.1 Basic Concept 3 1.2 Historical Background 3 1.3 General Applicability of the Method 7 1.4 Engineering Applications of the Finite Element Method 10 1.5 General Description of the Finite Element Method 10 1.6 Comparison of Finite Element Method with Other Methods of Analysis
Bruksanvisning för bilstereo . Bruksanvisning for bilstereo . Instrukcja obsługi samochodowego odtwarzacza stereo . Operating Instructions for Car Stereo . 610-104 . SV . Bruksanvisning i original
The Finite Element Method: Linear Static and Dynamic Finite Element Analysis by T. J. R. Hughes, Dover Publications, 2000 The Finite Element Method Vol. 2 Solid Mechanics by O.C. Zienkiewicz and R.L. Taylor, Oxford : Butterworth Heinemann, 2000 Institute of Structural Engineering Method of Finite Elements II 2
The Generalized Finite Element Method (GFEM) presented in this paper combines and extends the best features of the finite element method with the help of meshless formulations based on the Partition of Unity Method. Although an input finite element mesh is used by the pro- . Probl
The Adventures of Tom Sawyer 4 of 353 She went to the open door and stood in it and looked out among the tomato vines and ‘jimpson’ weeds that constituted the garden. No Tom. So she lifted up her voice at an angle calculated for distance and shouted: ‘Y-o-u-u TOM!’ There was a slight noise behind her and she turned just
