A finite element method for nematic liquid crystals with variable degree of orientation

A finite element method for nematic liquid crystals with variable degree of orientation

Ricardo H. Nochetto,   Shawn W. Walker,   Wujun Zhang
Abstract

We consider the simplest one-constant 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 Frank-like energies and a double well potential. In particular, the Euler-Lagrange 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 quasi-gradient 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 non-trivial 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_6cQ-U.

11footnotetext: rhn@math.umd.edu,  walker@math.lsu.edu,  wujun@umd.edu.

Key words. liquid crystals, finite element method, gamma-convergence, 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 one-constant 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 so-called one-constant Ericksen’s energy (LABEL:energy).

Despite the simple form of the one-constant Ericksen’s model, its minimizer may have non-trivial defects. If is a non-vanishing constant, then the energy reduces to the Oseen-Frank 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 Euler-Lagrange 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 quasi-gradient 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 quasi-gradient 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 quasi-gradient 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 non-trivial 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 vector-valued 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 bulk-energy functional which in the simplest one-constant model reduces to

 E[s,n]:=∫Ω(κ|∇s|2+s2|∇n|2)dx=:E1[s,n]+∫Ωψ(s)dx=:E2[s], \hb@xt@.01(2.1)

with and double well potential , which is a function defined on that satisfies

1. ,

2. for some ,

3. ;

see [26]. Note that when the degree of orientation equals a non-zero constant, the energy (LABEL:energy) effectively reduces to the Oseen-Frank energy . The degree of orientation relaxes the energy of defects (i.e. discontinuities in ), which may still have finite energy if the singular set

 S:={x∈Ω,s(x)=0} \hb@xt@.01(2.2)

is non-empty; in this case, .

By introducing an auxiliary variable [34, 3], we rewrite the energy as

 E1[s,n]=˜E1[s,u]:=∫Ω((κ−1)|∇s|2+|∇u|2)dx, \hb@xt@.01(2.3)

which follows from the orthogonal splitting due to the constraint . Accordingly, we define the admissible class

 A:={(s,u):Ω→(−1/2,1)×Rd: (s,u)∈[H1(Ω)]d+1,u=sn,n∈Sd−1}. \hb@xt@.01(2.4)

We say that the pair satisfies the structural condition for the Ericksen energy if

 u=sn,  −1/2

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

 A(g,r):={(s,u)∈A: s|Γs=g,u|Γu=r}, \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

 −12+δ0≤g(x),r(x)⋅ξ≤1−δ0∀x∈Rd,ξ∈Rd,|ξ|=1, \hb@xt@.01(2.7)

and the potential satisfies

 ψ(s)≥ψ(1−δ0)for s≥1−δ0,ψ(s)≥ψ(−12+δ0)for s≤−12+δ0. \hb@xt@.01(2.8)

This is consistent with property (1) of . If we further assume that

 g≥δ0 on ∂Ω, \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

 kij:=−∫Ω∇ϕi⋅∇ϕjdx≥0for % all i≠j, \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:weakly-acute) imposes a severe geometric restriction on [22, 45]. We recall the following characterization of (LABEL:weakly-acute) 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:weakly-acute) 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

 uh=Ih[shnh],−12

We then let and be the discrete Dirichlet data, and introduce the discrete spaces that include (Dirichlet) boundary conditions

 Sh(Γs,gh):={sh∈Sh:sh|Γs=gh},Uh(Γu,rh):={uh∈Uh:uh|Γu=rh},

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

 N∑j=1kij=−N∑j=1∫Ω∇ϕi⋅∇ϕjdx=0

because in the domain ; the set of hat functions is a partition of unity. Therefore, for piecewise linear , we have

 ∫Ω|∇sh|2dx=−N∑i=1kiish(xi)2−N∑i,j=1,i≠jkijsh(xi)sh(xj),

whence, exploiting and the symmetry , we get

 ∫Ω|∇sh|2dx =N∑i,j=1kijsh(xi)(sh(xi)−sh(xj)) \hb@xt@.01(2.14) =12N∑i,j=1kij(sh(xi)−sh(xj))2=12N∑i,j=1kij(δijsh)2,

where we define

 δijsh:=sh(xi)−sh(xj),δijnh:=nh(xi)−nh(xj). \hb@xt@.01(2.15)

With this in mind, we define the discrete energies to be

 Eh1[sh,nh]:=κ2N∑i,j=1kij(δijsh)2+12N∑i,j=1kij(sh(xi)2+sh(xj)22)|δijnh|2, \hb@xt@.01(2.16)

and

 Eh2[sh]:=∫Ωψ(sh(x))dx, \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

 Eh[sh,nh]:=Eh1[sh,nh]+Eh2[sh]. \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

 uh:=Ih[shnh]∈Uh,˜uh:=Ih[˜shnh]∈Uh. \hb@xt@.01(2.19)

Note that both pairs satisfy (LABEL:discrete-structure).

Lemma 2.2 (energy inequality)

Let the mesh satisfy (LABEL:weakly-acute). If , then, for any , the discrete energy (LABEL:discrete_energy_E1) satisfies

 Eh1[sh,nh]≥(κ−1)∫Ω|∇sh|2dx+∫Ω|∇uh|2dx=:˜Eh1[sh,uh], \hb@xt@.01(2.20)

as well as

 Eh1[sh,nh]≥(κ−1)∫Ω|∇˜sh|2dx+∫Ω|∇˜uh|2dx=:˜Eh1[˜sh,˜uh]. \hb@xt@.01(2.21)

Proof. Since

 sh(xi)nh(xi)−sh(xj)nh(xj) =sh(xi)+sh(xj)2(nh(xi)−nh(xj)) +(sh(xi)−sh(xj))nh(xi)+nh(xj)2,

using the orthogonality relation and (LABEL:eqn:dirichlet_integral_identity) yields

 ∫Ω|∇uh|2dx=12N∑i,j=1kij|sh(xi)nh(xi)−sh(xj)nh(xj)|2 = 12N∑i,j=1kij(sh(xi)+sh(xj)2)2|δijnh|2+12N∑i,j=1kij(δijsh)2∣∣∣nh(xi)+nh(xj)2∣∣∣2.

Exploiting the relations and , we obtain

 ∫Ω|∇uh|2dx =12N∑i,j=1kijsh(xi)2+sh(xj)22|δijnh|2 \hb@xt@.01(2.22) +12N∑i,j=1kij(δijsh)2−N∑i,j=1kij(δijsh)2∣∣∣nh(xi)−nh(xj)2∣∣∣2,

whence, we infer that

 Eh1[sh,nh]=∫Ω((κ−1)|∇sh|2+|∇uh|2)dx+N∑i,j=1kij(δijsh)2∣∣∣δijnh2∣∣∣2. \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 :

 ∫Ω|∇˜uh|2dx =12N∑i,j=1kij˜sh(xi)2+˜sh(xj)22|δijnh|2 \hb@xt@.01(2.24) +12N∑i,j=1kij(δij˜sh)2−N∑i,j=1kij(δij˜sh)2∣∣∣nh(xi)−nh(xj)2∣∣∣2.

We finally find that

 ˜Eh1[˜sh,˜uh] =∫Ω(|∇˜uh|2+(κ−1)|∇˜sh|2)dx=12N∑i,j=1kij˜sh(xi)2+˜sh(xj)22|δijnh|2 +κ2N∑i,j=1kij(δij˜sh)2−N∑i,j=1kij(δij˜sh)2∣∣∣nh(xi)−nh(xj)2∣∣∣2≤Eh1[sh,nh],

where we have dropped the last term and used the triangle inequality along with to obtain

 ∥∇˜sh∥2L2(Ω)=12N∑i,j=1kij(δij˜sh)2≤12N∑i,j=1kij(δijsh)2=∥∇sh∥2L2(Ω). \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 lim-sup) 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 semi-continuity (or lim-inf) 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

 X:=L2(Ω)×[L2(Ω)]d,Xh:=Sh×Uh.

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:lim-sup, 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:lim-inf to derive the weak lower semi-continuity property. Furthermore, our functionals exhibit the usual equi-coercivity property for both pairs and , but not for the director field , which is only well-defined whenever the order parameter . We discuss these issues in subsection LABEL:S:equi-coercivity and characterize the limits , and . We eventually prove -convergence in subsection LABEL:S:Gamma-convergence by combining these results.

3.1 Consistency or lim-sup property

We prove the following: if , then there exists a sequence converging to in and a discrete director field converging to in such that

 E1[s,n]≥limsuph→0Eh1[sh,nh]≥limsuph→0˜Eh1[sh,uh]. \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

 ^s(x)=min{1−δ0,max(−12+δ0,s(x))},^u(x)=^s(x)n(x)a.e.x∈Ω.

Then and

 E1[^s,n]≤E1[s,n],E2[^s]≤E2[s].

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:Dirichlet-restrict). Moreover, and the structural property (LABEL:structure) holds by construction, whence . We next observe that

 ∇^s=χΩ0∇s,Ω0:={x∈Ω:−12+δ0≤s(x)≤1−δ0};

[27, Ch. 5, Exercise 17]. Consequently, we obtain

 E1[^s,n]=∫Ωκ|∇^s|2+|^s|2|∇n|2≤∫Ωκ|∇s|2+|s|2|∇n|2=E1[s,n],

as well as

 E2[^s]=∫Ωψ(^s)≤∫Ωψ(s)=E2[s],

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 A(g,r))

Let , and let satisfy (LABEL:gne0). Given there exists a pair such that

 ∥(s,u)−(sϵ,uϵ)∥H1(Ω)≤ϵ, \hb@xt@.01(3.2)
 −12+δ0≤sϵ(x),uϵ(x)⋅ξ≤1−δ0∀x∈Ω,ξ∈Rd,|ξ|=1. \hb@xt@.01(3.3)

Proof. We construct a two-scale 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 non-negative mollifier with support contained in the ball centered at with radius . Define by

 dδ(x):=χΩ(x)min{δ−1dist(x,∂Ω),1},

which is Lipschitz in , and observe that is supported in the boundary layer

 ωδ:={x∈Ω:dist(x,∂Ω)≤δ},

and . We consider the Lipschitz approximations of given by

 sδ:=dδ(s∗ηδ)+(1−dδ)g,uδ:=dδ(u∗ηδ)+(1−dδ)r.

Since vanishes on we readily see that on . Moreover, the following properties are valid

 sδ→s,uδ→u,|uδ|→|u|a.e. and in H1(Ω). \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

 ∇(sδ−s)=∇dδ(s−g)∗ηδ+∇dδ(g∗ηδ−g)+dδ∇(s∗ηδ−s)+(dδ−1)∇(s−g).

Since , and on , we apply Poincaré’s inequality to deduce

 ∥s−g∥L2(ωδ)≤Cδ∥∇(s−g)∥L2(ωδ),

whence

 ∥∇dδ(s−g)∗ηδ∥L2(Ω)≤Cδ−1∥s−g∥L2(ωδ)≤C∥∇(s−g)∥L2(ωδ)→0  as δ→0.

Likewise, a similar argument gives for the fourth term

 ∥∥(dδ−1)∇(s−g)∥∥L2(Ω)≤C∥∇(s−g)∥L2(ωδ)→0 as δ→0.

On the other hand, the estimate yields

 ∥g∗ηδ−g∥L2(ωδ)≤|ωδ|1/2δ∥∇g∥L∞(Rd)≤Cδ32∥∇g∥L∞(Rd),

which implies, for the second term above,

 ∥∥∇dδ(g∗ηδ−g)∥∥L2(Ω)≤Cδ12∥∇g∥L∞(Rd).

Finally, for the third term we recall that equals outside and exploit the convergence in to obtain

 ∥dδ∇(s∗ηδ−s)∥L2(Ω)→0 as δ→0.

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:Dirichlet-restrict) in , whence so do the extensions of because , outside . Thus, we can show that

 −12+δ0≤sδ(x),uδ(x)⋅ξ≤1−δ0∀x∈Ω,ξ∈Rd,|ξ|=1;

we only argue with because dealing with is similar. We have because and the convolution preserves constants, whence

 sδ≤dδb+(1−dδ)b=b,sδ≥dδa+(1−dδ)a=a in Ω.

We next introduce the second parameter and the Lipschitz approximation of the sign function

 ρσ(t)=min{1,max{−1,t/σ}},

along with the two-scale approximation of

 sσ,δ:=ρσ(sδ)|uδ|,uσ,δ:=|ρσ(sδ)|uδ.

We note that by construction, whence (LABEL:structure) holds, and