The Ericksen Model of Liquid Crystals with Colloidal and Electric Effects
We present a robust discretization of the Ericksen model of liquid crystals with variable degree of orientation coupled with colloidal effects and electric fields. The total energy consists of the Ericksen energy, a weak anchoring (or penalized Dirichlet) energy to model colloids, and an electrical energy for a given electric field. We describe our special discretization of the total energy along with a method to compute minimizers via a discrete quasi-gradient flow algorithm which has a strictly monotone energy decreasing property. Numerical experiments are given in two and three dimensions to illustrate that the method is able to capture non-trivial defect patterns, such as the Saturn ring defect. We conclude with a rigorous proof of the -convergence of our discrete energy to the continuous energy.
keywords:liquid crystals, finite element method, gamma-convergence, gradient flow, line defect, plane defect, Saturn ring defect
Msc:65N30, 49M25, 35J70
This paper presents a method for solving the Ericksen model of liquid crystals (1); (2), with additional effects due to colloidal domains and electric fields. Liquid crystals are a work-horse technology enabling electronic displays (3); (4); (5), for instance. Moreover, they have a host of potential applications in material science (6); (7); (8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20); (21). One avenue is to use external fields (e.g. electric fields) and colloidal dispersions to build new materials through directed self-assembly (22); (7); (12); (23); (24); (25); (26); (27); (28); (29); (16); (19); (30); (31).
A significant amount of mathematical analysis has been done on liquid crystals (32); (33); (34); (35); (36); (37); (38); (39); (40); (41); (42); (43); (44). Moreover, a host of numerical methods have been developed for statics and dynamics (45); (46); (47); (48); (49); (50); (51). In particular, the methods in (52); (43); (53); (54); (55) are for harmonic mappings and liquid crystals with fixed degree of orientation, i.e. a unit vector field (called the director field) represents the orientation of liquid crystal molecules. See (56); (57); (58); (59); (60) for methods that couple liquid crystals to Stokes flow. We also refer to the survey paper (45) for more numerical methods.
The method we present (61); (62) is for the one-constant model of liquid crystals with variable degree of orientation (1); (2); (32) (Ericksen’s model). The state of the liquid crystal is described by a director field and a scalar function , the so-called degree-of-orientation, which minimize the energy
Hereafter, is a material constant, is a bounded Lipschitz domain in with , and is a double well potential (defined below).
Minimizers of the Ericksen model may exhibit non-trivial defects (depending on boundary conditions) (63); (10); (64); (39); (38); (65). The presence of in (1) leads to an Euler-Lagrange equation for that is degenerate. This allows for line and plane defects (singularities of ) in dimension when vanishes; these types of defects are important for applications, especially defects that lie on three dimensional space curves (7); (66). Regularity properties of minimizers, and the size of defects, were studied in (39). This leads to the study of dynamics (33) and corresponding numerics (46), which are relevant to our paper. But in both cases they regularize the model to avoid the degeneracy associated with the order parameter vanishing.
Our finite element method (FEM) does not require any regularization. We discussed the mathematical foundation of our method in (61); (62): we proved stability and convergence properties via -convergence (67) (as the mesh size goes to zero) and developed a quasi-gradient flow method to solve the discrete problem. Our discretization of the energy (1), defined in (16), requires that the mesh be weakly acute (or the stronger condition of having non-obtuse angles). This discretization preserves the underlying structure and robustly handles the unit length constraint on and the degeneracy present when vanishes. Our previous paper (61) showed a variety of simulations of minimizers with interesting defect structures.
The present paper demonstrates the ability of the Ericksen model, and of our method, to capture defect structures induced by colloidal inclusions (i.e. holes in the domain) and effects due to electric fields. We are able to recreate the famous Saturn ring defect (68), which occurs around colloidal particles in different situations, by using both a conforming (non-obtuse mesh) and a non-conforming cube mesh with an immersed boundary approach to model the colloid. In addition, we include electric field effects by incorporating an electric energy term into the total energy (1), and demonstrate the classic Freedericksz transition (32); (2); (41); (69); (70). We also investigate the coupling of colloidal and electric effects.
The paper is organized as follows. In Section 2, we recall the Ericksen model for liquid crystals with variable degree of orientation, and describe our discretization of the continuous energy. Section 3 recalls properties of the discretization, and our initial minimization scheme. Section 4 describes the details for properly implementing our method. Section 5.1 illustrates our method in the presence of a colloidal inclusion with a conforming non-obtuse mesh. Section 5.2 shows an alternative way to model colloids by an immersed boundary approach (along with supporting simulations). In Section 6, we show how to include electric field effects in the model and describe a modified minimization procedure to compute minimizers. Section 7 presents the monotone energy decreasing property of the quasi-gradient flow algorithm to compute discrete minimizers. Section 8 provides a summary of the -convergence for our discrete energy. We close in Section 9 with some discussion.
2 Ericksen’s model
Let the director field be a vector-valued function with unit length (see Figure 1 for a description of the meaning of ). The degree-of-orientation is a real valued function (see Figure 2 for a description of the meaning of ). The variable , by itself, cannot properly describe a “loss of order” in the liquid crystal material because it has unit length. The variable provides a way to characterize the local order (see Figure 2).
2.1 Ericksen’s one constant energy
The equilibrium state of the liquid crystal material is described by the pair minimizing a bulk-energy functional (1) which we split as
for some ,
We may also enforce boundary conditions on , possibly on different parts of the boundary. Let be open subsets of where we set Dirichlet boundary conditions for . Then we have the following restricted admissible class
for some given functions that satisfy the following in a neighborhood of : and , for some . Note that if we further assume
then the function is in a neighborhood of and satisfies on .
When the degree of orientation is a non-zero constant, the energy in (2) effectively reduces to the Oseen-Frank energy . The purpose of the degree of orientation is to relax the energy of defects. In fact, discontinuities in (i.e. defects) may still occur in the singular set
with finite energy: . The existence of minimizers in the admissible class, subject to Dirichlet boundary conditions, was shown in (35); (39). Minimizers with defects are constructed explicitly in (32) or discovered numerically in (61).
The parameter in (2) plays a major role in the occurrence of defects. Assuming the boundary condition for is a positive constant well away from zero, if is large, then dominates the energy and stays close to a positive constant within the domain . Thus, defects are less likely to occur. If is small (say ), then dominates the energy, and may vanish in regions of and induce a defect. This was confirmed by our numerical experiments in (61); (62). The physically relevant case is the more difficult case with regard to proving -convergence (see (61)) because the energy is no longer convex.
2.2 Discretization of the energy
Let be a conforming simplicial triangulation of . The set of nodes (vertices) of is denoted and has cardinality . We demand that be weakly acute, namely
where is the standard “hat” basis function associated with node . We indicate with the patch of a node (i.e. the “star” of elements in that contain the vertex ). Of course, (8) imposes a severe geometric restriction on (71); (72) (especially in three dimensions). We recall the following characterization of (8) for .
Lemma 1 (weak acuteness in two dimensions).
For any pair of triangles , in in two space dimensions that share a common edge , let be the angle in opposite to (for ). Then (8) holds if and only if for every edge .
Generalizations of Lemma 1 to three dimensions, involving interior dihedral angles of tetrahedra, can be found in (73); (74). Note that a non-obtuse mesh (one where all interior angles are bounded by ) is automatically weakly-acute.
The method uses the following finite element spaces:
where imposes the unit length constraint at the vertices of the mesh.
Let denote the piecewise linear Lagrange interpolation operator on mesh with values in either or . We have the following discrete version of the admissible class:
Next, we let and be the discrete Dirichlet data, and introduce the discrete spaces that include (Dirichlet) boundary conditions
as well as the discrete admissible class with boundary conditions:
In view of (6), we can also impose the Dirichlet condition on .
Our discrete version of is “derived” by invoking basic properties of the stiffness matrix entries . First note for all , and for piecewise linear we have
Thus, using and the symmetry , we get
where we define
Therefore, we define the discrete energy to be
where the second term is a first order approximation of , which is novel in the finite element literature (61). The double well energy is discretized in the usual way:
The specific form of (14) lies on the fact that it makes the nodal values of and readily accessible for analysis. The identity in (3) is obtained (at the continuous level) by taking advantage of the unit length constraint . However, at the discrete level, we only impose the unit length constraint at the nodes of the mesh, and we cannot hope for much more because is a piecewise polynomial. Hence, in order to obtain a similar identity to (3) (see Lemma 2 below), we need access to nodal values. In (61), we show that the discrete energy (14) allows us to handle the degenerate coefficient without regularization.
The discrete formulation is as follows. Find such that the following energy is minimized:
3 Review of the method
3.1 Energy inequality
Note that both , are in . We now state a discrete version of (3).
Lemma 2 (discrete energy inequality).
3.2 Minimization scheme
We summarize the discrete quasi-gradient flow scheme in (61) which we use to compute discrete minimizers.
In the continuous setting, Dirichlet boundary conditions are enforced in the space. Let
where is either , and , . We assume and (6) to be valid on . The trace is thus well defined on .
The superscript will stand for an iteration counter. Therefore, and indicate iterates satisfying Dirichlet boundary conditions where and . We will further simplify the notation in some places upon writing:
First order variation
We start with the energy . Due to the unit length constraint at the nodes in (see (9)), we introduce the space of discrete tangential variations:
Next, the first order variation of in the direction reads
The first order variation of in the direction consists of two terms
We next consider the energy . In order to guarantee a monotonically energy decreasing scheme, we employ the convex splitting technique in (75); (76); (77), i.e. we split the double well potential into a convex and concave part. Let and be both convex for all so that , and set
which yields the inequality
for any and in (61). Note that
Discrete quasi-gradient flow algorithm
Our scheme for minimizing , defined in (16), is an alternating direction method, which minimizes with respect to and evolves separately in the steepest descent direction during each iteration. Therefore, this algorithm is not a standard gradient flow but rather a quasi-gradient flow.
Algorithm: Given in , iterate Steps (a)-(c) for .
Step (a): Minimization. Find such that minimizes the energy for all in , i.e. satisfies
Step (b): Projection. Normalize at all nodes .
Step (c): Gradient flow. Using , find in such that
In the numerical experiments, we impose Dirichlet boundary conditions for both and on subsets of the boundary. Note that the scheme has no restriction on the time step thanks to the implicit Euler method in Step (c).
Theorem 3 (monotonicity (61)).
We implemented our method using the MATLAB/C++ finite element toolbox FELICITY (78). In this section, we give details on forming the ensuing discrete systems, and how to solve part (a) of the quasi-gradient flow algorithm in 3-D using the tangent space. For all 3-D simulations, we used the algebraic multi-grid solver (AGMG) (79); (80); (81); (82) to solve the linear systems in parts (a) and (c) of the quasi-gradient flow algorithm. In 2-D, we simply used the “backslash” command in MATLAB.
4.1 Finite element matrices
Implementing the algorithm requires construction of the discrete energy, as well as its variational derivative. This requires the symmetric mass and stiffness finite element matrices: , , where
and is the set of basis functions of the space .
4.2 Finite element functions and coefficient vectors
The function is represented by a linear combination of . If the dimension of is , then the vector field has components, where each component is written as a linear combination of . The nodal values of and , at node , are denoted by and .
The corresponding coefficient vectors (arrays) are denoted with non-italicized capital letters, i.e. , , such that
for , where are the canonical basis vectors of . In other words, we store the coefficients of so that the components are first, followed by the components, and so on. Therefore,
4.3 Discrete variations
Let us write from (25) in a different form:
where is the coefficient vector corresponding to , and is the symmetric matrix defined by
Let be an arbitrary matrix, and , be arbitrary column vectors. Then
where is a column vector and is a diagonal matrix formed from .
Next, we write out from (24) in a different form:
where we defined
with the coefficient vector corresponding to .
Let us now focus on :
where is the symmetric matrix defined by
Using symmetry gives
where denotes the first components of , etc., and
4.4 Discrete quasi-gradient flow
Given , we have the corresponding coefficient vectors . We now rewrite the Algorithm in Section 3.2.3 in terms of the matrices and vectors introduced earlier.
where in is the coefficient vector corresponding to in . Note that (43) must be modified to enforce Dirichlet boundary conditions (if necessary).
Remark 5 (solving a degenerate system).
The system matrix in (43) is symmetric positive semi-definite, which is easily verified from the properties of . Moreover, it is positive definite if everywhere. Hence, the system can be solved by any method for symmetric positive definite matrices.
When at a sufficient number of nodes, the matrix will be singular. In this case, one could use a conjugate gradient method (83); note that the right-hand-side of (43) is guaranteed to be in the column space of the system matrix.