A finite element method for nematic liquid crystals with variable degree of orientation
Abstract
We consider the simplest oneconstant model, put forward by J. Ericksen, for nematic liquid crystals with variable degree of orientation. The equilibrium state is described by a director field and its degree of orientation , where the pair minimizes a sum of Franklike energies and a double well potential. In particular, the EulerLagrange equations for the minimizer contain a degenerate elliptic equation for , which allows for line and plane defects to have finite energy.
We present a structure preserving discretization of the liquid crystal energy with piecewise linear finite elements that can handle the degenerate elliptic part without regularization, and show that it is consistent and stable. We prove convergence of discrete global minimizers to continuous ones as the mesh size goes to zero. We develop a quasigradient flow scheme for computing discrete equilibrium solutions and prove it has a strictly monotone energy decreasing property. We present simulations in two and three dimensions to illustrate the method’s ability to handle nontrivial defects.
A music video summary of the paper is available on YouTube: “Mathematical Modeling and Simulation of Nematic Liquid Crystals (A Montage),” http://www.youtube.com/watch?v=pWWw7_6cQU.
Key words. liquid crystals, finite element method, gammaconvergence, gradient flow, line defect, plane defect
AMS subject classifications. 65N30, 49M25, 35J70
1 Introduction
Complex fluids are ubiquitous in nature and industrial processes and are critical for modern engineering systems [32, 41, 15]. An important difficulty in modeling and simulating complex fluids is their inherent microstructure. Manipulating the microstructure via external forces can enable control of the mechanical, chemical, optical, or thermal properties of the material. Liquid crystals [47, 25, 21, 4, 3, 13, 7, 33, 34, 5, 46] are a relatively simple example of a material with microstructure that may be immersed in a fluid with a free interface [53, 52].
Several numerical methods for liquid crystals have been proposed in [10, 29, 23, 35, 2] for harmonic mappings and liquid crystals with fixed degree of orientation, i.e. a unit vector field (called the director field) is used to represent the orientation of liquid crystal molecules. See [28, 36, 49] for methods that couple liquid crystals to Stokes flow. We also refer to the survey paper [6] for more numerical methods.
In this paper, we consider the oneconstant model for liquid crystals with variable degree of orientation [26, 25, 47]. The state of the liquid crystal is described by a director field and a scalar function , , that represents the degree of alignment that molecules have with respect to . The equilibrium state is given by which minimizes the socalled oneconstant Ericksen’s energy (LABEL:energy).
Despite the simple form of the oneconstant Ericksen’s model, its minimizer may have nontrivial defects. If is a nonvanishing constant, then the energy reduces to the OseenFrank energy whose minimizers are harmonic maps that may exhibit point defects (depending on boundary conditions) [14, 16, 20, 34, 33, 42]. If is part of the minimization of (LABEL:energy), then may vanish to allow for line (and plane) defects in dimension [5, 46], and the resulting EulerLagrange equation for is degenerate. However, in [34], it was shown that both and have strong limits, which enabled the study of regularity properties of minimizers and the size of defects. This inspired the study of dynamics [21] and corresponding numerics [8], which are most relevant to our paper. However, in both cases they regularize the model to avoid the degeneracy introduced by the parameter.
We design a finite element method (FEM) without any regularization. We prove stability and convergence properties and explore equilibrium configurations of liquid crystals via quasigradient flows. Our method builds on [12, 9, 11] and consists of a structure preserving discretization of (LABEL:energy). Given a weakly acute mesh with mesh size (see Section LABEL:subsec:FEM), we use the subscript to denote continuous piecewise linear functions defined over , e.g. is a discrete approximation of .
Our discretization of the energy is defined in (LABEL:discrete_energy) and requires that be weakly acute. This discretization preserves the underlying structure and converges to the continuous energy in the sense of convergence [17] as goes to zero. Next, we develop a quasigradient flow scheme for computing discrete equilibrium solutions. We prove that this scheme has a strictly monotone energy decreasing property. Finally, we carry out numerical experiments and show that our finite element method, and gradient flow, allows for computing minimizers that exhibit line and plane defects.
The paper is organized as follows. In Section LABEL:sec:Discretization, we describe the Ericksen model for liquid crystals with variable degree of orientation, as well as the details of our discretization. Section LABEL:sec:consistency shows the convergence of our numerical method. A quasigradient flow scheme is given in Section LABEL:sec:gradient_flow, where we also prove a strictly monotone energy decreasing property. Section LABEL:sec:numerics presents simulations in two and three dimensions that exhibit nontrivial defects in order to illustrate the method’s capabilities.
2 Discretization of Ericksen’s model
We review the model [26] and relevant analysis results from the literature. We then develop our discretization strategy and show it is stable. The space dimension can be arbitrary.
2.1 Ericksen’s one constant model
Let the director field be a vectorvalued function with unit length, and the degree of orientation be a real valued function. The case represents the state of perfect alignment in which all molecules are parallel to . Likewise, represents the state of microscopic order in which all molecules are orthogonal to the orientation . When , the molecules do not lie along any preferred direction which represents the state of an isotropic distribution of molecules.
The equilibrium state of the liquid crystals is described by the pair minimizing a bulkenergy functional which in the simplest oneconstant model reduces to
\hb@xt@.01(2.1) 
with and double well potential , which is a function defined on that satisfies

,

for some ,

;
see [26]. Note that when the degree of orientation equals a nonzero constant, the energy (LABEL:energy) effectively reduces to the OseenFrank energy . The degree of orientation relaxes the energy of defects (i.e. discontinuities in ), which may still have finite energy if the singular set
\hb@xt@.01(2.2) 
is nonempty; in this case, .
By introducing an auxiliary variable [34, 3], we rewrite the energy as
\hb@xt@.01(2.3) 
which follows from the orthogonal splitting due to the constraint . Accordingly, we define the admissible class
\hb@xt@.01(2.4) 
We say that the pair satisfies the structural condition for the Ericksen energy if
\hb@xt@.01(2.5) 
Moreover, we may 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
\hb@xt@.01(2.6) 
for some given functions that satisfy the structural condition (LABEL:structure) on . We assume the existence of sufficiently small such that
\hb@xt@.01(2.7) 
and the potential satisfies
\hb@xt@.01(2.8) 
This is consistent with property (1) of . If we further assume that
\hb@xt@.01(2.9) 
then the function is in a neighborhood of and satisfies on .
The existence of a minimizer is shown in [34, 3], but this is also a consequence of our convergence theory. It is worth mentioning that the constant in (LABEL:energy) plays a significant role in the occurrence of defects. Roughly speaking, if is large, then dominates the energy and is close to a constant. In this case, defects with finite energy are less likely to occur. But if is small, then dominates the energy, and may become zero. In this case, defects are more likely to occur. (This heuristic argument is later confirmed in the numerical experiments.) Since the investigation of defects is of primary interest in this paper, we consider the most significant case to be .
We now describe our finite element discretization of the energy (LABEL:energy) and its minimizer .
2.2 Discretization of the energy
Let be a conforming simplicial triangulation of the domain . We denote by the set of nodes (vertices) of and the cardinality of by (with some abuse of notation). We demand that be weakly acute, namely
\hb@xt@.01(2.10) 
where is the standard “hat” function associated with node . We indicate with the patch of a node (i.e. the “star” of elements in that contain the vertex ). Condition (LABEL:weaklyacute) imposes a severe geometric restriction on [22, 45]. We recall the following characterization of (LABEL:weaklyacute) for .
Proposition 2.1 (weak acuteness in two dimensions)
For any pair of triangles , in that share a common edge , let be the angle in opposite to (for ). If for every edge , then (LABEL:weaklyacute) holds.
Generalizations of Proposition LABEL:weak_acuteness_2D to three dimensions, involving interior dihedral angles of tetrahedra, can be found in [30, 19].
We construct continuous piecewise affine spaces associated with the mesh, i.e.
\hb@xt@.01(2.11) 
Let denote the piecewise linear Lagrange interpolation operator on mesh with values in either or . We say that a pair satisfies the discrete structural condition for the Ericksen energy if there exists such that
\hb@xt@.01(2.12) 
We then let and be the discrete Dirichlet data, and introduce the discrete spaces that include (Dirichlet) boundary conditions
as well as the discrete admissible class
\hb@xt@.01(2.13) 
In view of (LABEL:gne0), we can also impose the Dirichlet condition on .
In order to motivate our discrete version of , note that for all
because in the domain ; the set of hat functions is a partition of unity. Therefore, for piecewise linear , we have
whence, exploiting and the symmetry , we get
\hb@xt@.01(2.14)  
where we define
\hb@xt@.01(2.15) 
With this in mind, we define the discrete energies to be
\hb@xt@.01(2.16) 
and
\hb@xt@.01(2.17) 
for . The second summation in (LABEL:discrete_energy_E1) does not come from applying the standard discretization of by piecewise linear elements. It turns out that this special form of the discrete energy preserves the key energy inequality (Lemma LABEL:lemma:energydecreasing) which allows us to establish our convergence analysis for the degenerate coefficient without regularization. Eventually, we seek an approximation of the pair such that the discrete pair minimizes the discrete version of the bulk energy (LABEL:energy) given by
\hb@xt@.01(2.18) 
The following result shows that definition (LABEL:discrete_energy_E1) preserves the key structure (LABEL:auxiliary_energy_identity) of [3, 34] at the discrete level, which turns out to be crucial for our analysis as well.
We first introduce and two discrete versions of the vector field
\hb@xt@.01(2.19) 
Note that both pairs satisfy (LABEL:discretestructure).
Lemma 2.2 (energy inequality)
Let the mesh satisfy (LABEL:weaklyacute). If , then, for any , the discrete energy (LABEL:discrete_energy_E1) satisfies
\hb@xt@.01(2.20) 
as well as
\hb@xt@.01(2.21) 
Proof. Since
using the orthogonality relation and (LABEL:eqn:dirichlet_integral_identity) yields
Exploiting the relations and , we obtain
\hb@xt@.01(2.22)  
whence, we infer that
\hb@xt@.01(2.23) 
The inequality (LABEL:energy_inequality) follows directly from for .
To prove (LABEL:abs_inequality), we note that (LABEL:lem:identity) still holds if we replace with :
\hb@xt@.01(2.24)  
We finally find that
where we have dropped the last term and used the triangle inequality along with to obtain
\hb@xt@.01(2.25) 
This concludes the proof.
Remark 2.3 (relation between (LABEL:energy_inequality) and (LABEL:abs_inequality))
Both (LABEL:energy_inequality) and (LABEL:abs_inequality) account for the variational crime committed when enforcing and only at the vertices, and mimics (LABEL:auxiliary_energy_identity). Since we have a precise control of the consistency error for (LABEL:energy_inequality), this inequality will be used later for the consistency (or limsup) step of convergence of our discrete energy (LABEL:discrete_energy_E1) to the original continuous energy in (LABEL:energy). On the other hand, (LABEL:abs_inequality) has a suitable structure to prove the weak lower semicontinuity (or liminf) step of convergence. This property is not obvious when , the most significant case for the formation of defects.
3 convergence of the discrete energy
In this section, we show that our discrete energy (LABEL:discrete_energy_E1) converges to the continuous energy (LABEL:energy) in the sense of convergence. To this end, we first let the continuous and discrete spaces be
We next define as in (LABEL:energy) for and for . Likewise, we define as in (LABEL:discrete_energy) for and for all
We split the proof of convergence into four subsections. In subsection LABEL:S:limsup, we use the energy to show the consistency property (recall Remark LABEL:rem:why_grad_I_h_abs_S), whereas we employ the energy in subsection LABEL:S:liminf to derive the weak lower semicontinuity property. Furthermore, our functionals exhibit the usual equicoercivity property for both pairs and , but not for the director field , which is only welldefined whenever the order parameter . We discuss these issues in subsection LABEL:S:equicoercivity and characterize the limits , and . We eventually prove convergence in subsection LABEL:S:Gammaconvergence by combining these results.
3.1 Consistency or limsup property
We prove the following: if , then there exists a sequence converging to in and a discrete director field converging to in such that
\hb@xt@.01(3.1) 
We observe that if , then and (LABEL:limsup) is valid for any sequences and in light of (LABEL:energy_inequality).
We first show that we can always assume for .
Lemma 3.1 (truncation)
Given , let be the truncations
Then and
The same assertion is true for any except that the truncations are defined nodewise, i.e. .
Proof. The fact that satisfy the Dirichlet boundary conditions is a consequence of (LABEL:Dirichletrestrict). Moreover, and the structural property (LABEL:structure) holds by construction, whence . We next observe that
[27, Ch. 5, Exercise 17]. Consequently, we obtain
as well as
because of (LABEL:potential). This concludes the proof.
To construct a recovery sequence we need point values of and thus a regularization procedure of functions in the admissible class . We must enforce both the structural property (LABEL:structure) and the Dirichlet boundary conditions and ; neither one is guaranteed by convolution. We are able to do this provided and the Dirichlet datum satisfies (LABEL:gne0).
Proposition 3.2 (regularization of functions in )
Let , and let satisfy (LABEL:gne0). Given there exists a pair such that
\hb@xt@.01(3.2) 
\hb@xt@.01(3.3) 
Proof. We construct a twoscale approximation with scales , which satisfies the boundary conditions exactly. We split the argument into several steps.
Step 1: Regularization with Dirichlet condition. Extend by zero to . Let be a smooth and nonnegative mollifier with support contained in the ball centered at with radius . Define by
which is Lipschitz in , and observe that is supported in the boundary layer
and . We consider the Lipschitz approximations of given by
Since vanishes on we readily see that on . Moreover, the following properties are valid
\hb@xt@.01(3.4) 
The last property is a consequence of the middle one via triangle inequality, and the first two are similar. It thus suffices to show the first property for . We simply write
Since , and on , we apply Poincaré’s inequality to deduce
whence
Likewise, a similar argument gives for the fourth term
On the other hand, the estimate yields
which implies, for the second term above,
Finally, for the third term we recall that equals outside and exploit the convergence in to obtain
Step 2: Structural condition. The pair does not satisfy the structural condition (LABEL:structure) unfortunately. We now construct a closely related pair that satisfies (LABEL:structure). Recall that satisfy the bounds (LABEL:Dirichletrestrict) in , whence so do the extensions of because , outside . Thus, we can show that
we only argue with because dealing with is similar. We have because and the convolution preserves constants, whence
We next introduce the second parameter and the Lipschitz approximation of the sign function
along with the twoscale approximation of
We note that by construction, whence (LABEL:structure) holds, and