Time-splitting approximation of nematic liquid crystal flows

# A projection-based time-splitting algorithm for approximating nematic liquid crystal flows with stretching

R.C. Cabrales F. Guillén-González  and  J. V. Gutiérrez-Santacreu
###### Abstract.

A numerical method is developed for solving a system of partial differential equations modeling the flow of a nematic liquid crystal fluid with stretching effect, which takes into account the geometrical shape of its molecules. This system couples the velocity vector, the scalar pressure and the director vector representing the direction along which the molecules are oriented. The scheme is designed by using finite elements in space and a time-splitting algorithm to uncouple the calculation of the variables: the velocity and pressure are computed by using a projection-based algorithm and the director is computed jointly to an auxiliary variable. Moreover, the computation of this auxiliary variable can be avoided at the discrete level by using piecewise constant finite elements in its approximation. Finally, we use a pressure stabilization technique allowing a stable equal-order interpolation for the velocity and the pressure. Numerical experiments concerning annihilation of singularities are presented to show the stability and efficiency of the scheme.

Departamento de Ciencias Básicas, Universidad del Bío-Bío, Casilla 447, Chillán, Chile. E-mail: roberto.cabrales@gmail.com. Partially supported under Chilean grant fondecyt 1140074.
Dpto. E.D.A.N. and IMUS, Universidad de Sevilla, Aptdo. 1160, 41080 Sevilla, Spain. E-mail: guillen@us.es. Partially supported by Ministerio de Economía y Competitividad under Spanish grant MTM2015-69875-P with the participation of FEDER
Dpto. de Matemática Aplicada I, Universidad de Sevilla, E. T. S. I. Informática. Avda. Reina Mercedes, s/n. 41012 Sevilla, Spain. juanvi@us.es. Partially supported by Ministerio de Economía y Competitividad under Spanish grant MTM2015-69875-P with the participation of FEDER

Mathematics Subject Classification: Nematic liquid crystal; Finite elements; Projection method; Time-splitting method.

Keywords: 35Q35, 65M60, 76A15

## 1. Introduction

For a long time we have all believed that matter only existed in three states: solid, liquid, and gas. However, liquid crystals are substances that combine features of both isotropic liquids and crystalline solids, exhibiting intermediate transitions between solid and liquid phases, called mesophases. This behavior is due, in part, to the fact that liquid crystals are made up of macromolecules of similar size. Moreover it is well-known that the shape of the molecules plays an important role in the flow regimes of liquid crystal fluids. Thus, in general, different behaviors are expected for the dynamics of liquid crystal fluids being constituted by molecules of different shapes. For instance, liquid crystal fluids built from disk-shaped molecules may exhibit strikingly different properties from those composed of rod-shape materials.

The mathematical theory for the hydrodynamics of liquid crystal fluids was initiated by Ericksen [8, 9], and Leslie [15, 16]. Such a theory describes liquid crystals according to the different degrees of positional and orientational ordering of their molecules. The former refers to the relative position, on average, of the molecules or groups of molecules, while the latter alludes to the fact that the molecules tend to be locally aligned toward a specific and preferred direction, described by a unit vector, called director, defined according to the form of the molecules.

Within the liquid crystal phases, we find the so-called nematic phase. In it, the molecules have no positional order, but they self-organize to have long-range directional order. Thus, the molecules flow freely and their center of mass positions are randomly distributed as in a liquid, while maintain their long-range directional order. The breakdown of the self-alignment causes the appearance of defects or singularities which are able to significantly influence the flow behavior.

In this paper we are interested in the numerical solution of the flow of a nematic liquid crystal governed by an Ericksen-Leslie-type system which incorporates the stretching effect. This stretching effect comes from the kinematic transport of the director field, which depends on the shape of the molecules.

Let be a fixed time and let or , be a bounded open set with boundary . Set and . Then the equations are written as follows (see [17], [18] for more physical background, derivation, and discussion):

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∂td+(u⋅∇)d+β(∇u)d+(1+β)(∇u)Td+γ(f(ε,d)−δd)=0 in Q,∂tu+(u⋅∇)u−νδu+∇p+λ∇⋅((∇d)T∇d)+λ∇⋅(β(f(ε,d)−δd)dT+(1+β)d(f(ε,d)−δd)T)=0 in Q,∇⋅u=0 in Q. (1)

We complete this system with homogeneous Dirichlet conditions for the velocity field (non-slip) and homogeneous Neumann boundary conditions for the director field:

 u(x,t)=0,∂nd(x,t)=0% for (x,t)∈Σ, (2)

and the initial conditions

 d({x},0)=d0({x}),u({x},0)=u0({x})% for x∈Ω. (3)

In system (1), is the velocity of the liquid crystal flow, is the pressure, and is the orientation of the molecules. The physical parameters , and stand for positive constants which represent viscosity, elasticity, and relaxation time, respectively. The geometrical parameter is a constant associated with the aspect ratio of the ellipsoid particles. For instance, and , corresponds to spherical, rod-like and disk-like liquid crystal molecules [14], respectively. Moreover, is the penalty parameter and is the penalty function related to the constraint . This penalty term can also be physically meaningful and represents a possible extensibility of molecules. It is defined by

 f(ε,d)=1ε2(|d|2−1)d. (4)

It should be noted that , for all , for the following scalar potential function

 F(ε,d)=14ε2(|d|2−1)2.

The following energy law for system (1) holds under some regularity assumptions for and [17, 18]:

 ddtE(u,d)+ν∫Ω|∇u|2+λγ∫Ω|−δd+f(ε,d)|2=0, (5)

where

 E(u,d)=12∫Ω|u|2+λ2∫Ω|∇d|2+λ∫ΩF(ε,d). (6)

The energy (6) expresses the competition between the kinetic and elastic energies. It should be noted that (5) is independent of and . Moreover, it provides a uniform bound for which leads to the unit sphere constraint in the limit as .

The mathematical structure of system (1) consists of the Navier-Stokes equations with some extra stress tensors taking into account the geometric of the molecules coupled with a gradient flow equation similar to that of harmonic maps into the sphere. The numerical resolution of system (1) is a nontrivial task due to the nonlinear nature of the system, the coupling terms between orientation and flow and the presence of the incompressibility constraint.

It is well known that time discretizations for the Navier–Stokes equations based on implicit time strategies give rise to highly time-consuming linear solver steps due to the coupled computation of both velocity and pressure. Furthermore, choices of stable finite element pairings are restricted by the inf–sup constraint. As a result, projection-based time-splitting techniques [11] turn to be a very favorable alternative to cutting down the computational cost of current iteration by successively updating velocity and pressure. It is obvious that such a strategy is desirable to solve system (1) in order to decouple the pressure computation from the velocity field as well as the computation of the director field. Instead, inf-sup conditions can be avoided by adding a pressure stabilizing term at the projection step, such that equal interpolation spaces for velocity and pressure are allowed. The stabilizing term is devised in such a way that the stability and convergence rate are not compromised [5].

There is an extensive literature on the mathematical analysis of finite element methods approximating system (1) without stretching effect, which corresponds to neglecting all the terms involving the paramter . The interested reader is referred to the works [20, 19, 4, 10, 2, 3, 12, 6]. Very little has appeared for system (1) in the context of finite elements. In [17], Lin, Liu and Zhang presented a modified Crank-Nicolson scheme for which a discrete energy law is derived. As a solver, the author devised a fixed iterative method which gave rise to a matrix being symmetric and independent of time and the number of the fixed point iterations at each iteration. In [18], Liu, Lin and Zhang proposed a numerical algorithm based on a semi-implicit BDF2 rotational pressure-correction method for an axi-symmetric domain for studying the annihilation of a hedgehog-antihedgehog pair of defects. An extra term related to was added so that the proposed scheme was somewhat unconditional, but no energy estimate was provided.

The goal of this paper is to use the above-mentioned techniques for the Navier-Stokes equations so as to be able to potentially enhance the performance of previous algorithms for system (1). In particular, we look for a numerical scheme using low-order finite-elements, which is linear at each time step and unconditionally stable and decouples the computation of the all primary variables. Moreover, we investigate the interplay of the flow of a nematic fluid and the geometric shape of its molecules in several numerical experiments.

The outline of the rest of this paper is the following. In section 2 we establish the function spaces, the notation and the hypotheses used in the work. Then the new numerical method are introduced. In section 3, we prove a priori estimates for the algorithm which provides the unconditional energy-stability property. The paper finishes with section 4, where some details about the implementation and numerical simulations that illustrate the performance of the scheme concerning the evolution of singularities are presented.

## 2. The numerical algorithm

### 2.1. Notation

For , denotes the space of th-power integrable real-valued functions defined on for the Lebesgue measure. This space is a Banach space endowed with the norm for or for . In particular, is a Hilbert space with the inner product

 (u,v)=∫Ωu(x)v(x)dx,

and its norm is simply denoted by . For a non-negative integer, we denote the classical Sobolev spaces as

 Hm(Ω)={v∈L2(Ω);∂kv∈L2(Ω) ∀ |k|≤m},

associated to the norm

 ∥v∥Hm(Ω)=⎛⎝∑0≤|k|≤m∥∂kv∥2⎞⎠1/2,

where is a multi-index and , which is a Hilbert space with the obvious inner product. We will use boldfaced letters for spaces of vector functions and their elements, e.g. in place of .

Let be the space of infinitely times differentiable functions with compact support on . The closure of in is denoted by . We will also make use of the following space of vector fields:

 V={vv∈D(Ω):∇⋅vv=0 in Ω}.

We denote by and , the closures of , in the - and -norm, respectively, which are characterized (for being Lipschitz-continuous) by (see [22])

 H = {u∈L2(Ω):∇⋅u=0 in Ω,u⋅n=0% on ∂Ω}, V = {u∈H1(Ω):∇⋅u=0 in Ω,u=0 on ∂Ω},

where is the outward normal to on . Finally, we consider

 L20(Ω)={p∈L2(Ω): ∫Ωp(x)dx=0}.

### 2.2. Hypotheses

Herein we introduce the hypotheses that will be required along this work.

1. Let be a bounded domain of with a polygonal or polyhedral Lipschitz-continuous boundary.

2. Let be a family of regular, quasi-uniform subdivisions of made up of triangles or quadrilaterals in two dimensions and tetrahedra or hexahedra in three dimensions, so that .

3. Assume three sequences of finite-dimensional spaces , and associated with such that , and . Also, consider an extra finite-element space with .

4. Suppose that with a.e. in .

In particular, hypothesis allows us to consider equal-order finite-element spaces for velocity and pressure. For instance, let be the set of linear polynomials on a triangle or tetrahedron . Thus the space of continuous, piecewise polynomial functions associated to is denoted as

 Xh={vh∈C0(¯¯¯¯Ω):vh|K∈P1(K), ∀K∈Th},

and the set of piecewise constant functions as

 Yh={wh∈L∞(Ω):wh|K∈R, ∀K∈Th}.

We choose the following continuous finite-element spaces

 dh=Xh,Vh=Xh∩H10(Ω)andPh=Xh∩L20(Ω),

for approximating the director, the velocity and the pressure, respectively. Additionally, we select the extra discontinuous finite-element to be the space for an auxiliary variable related to the vector director.

Observe that our choice of the finite-element spaces for velocity and pressure does not satisfy the discrete inf-sup condition

 ∥ph∥L20(Ω)≤αsupvvh∈Vh∖{0}(qh,∇⋅vvh)∥vvh∥H1(Ω)∀ph∈Ph, (7)

for independent of .

The following proposition is concerned with an interpolation operator associated with the space . In fact, we can think of as the Scott-Zhang interpolation operator, see [21].

###### Proposition 1.

Assuming hypotheses -, there exists an interpolation operator satisfying

 ∥d−Ihd∥≤Capph∥∇d∥∀d∈H1(Ω), (8)
 ∥Ihd∥L∞(Ω)≤Csta∥d∥L∞(Ω)∀d∈L∞(Ω), (9) ∥Ihd∥H1(Ω)≤Csta∥d∥H1(Ω)∀d∈H1(Ω), (10)

where and are constants independent of .

### 2.3. Description of the scheme

As explained in the introduction, we aim to construct a numerical solution to system (1) that, at each time step, one only needs to solve a sequence of decoupled elliptic equations for director, velocity and pressure. In particular, the linear systems associated to director and pressure are symmetric; therefore, scalable parallel solvers can be defined. Instead, the linear system associated to velocity is block diagonal, which means that each component of velocity can be computed in parallel. This makes our time-splitting be very appealing for high performance computing.

The starting point to design our time-splitting method is the non-incremental velocity-correction method for the Navier-Stokes equations. Moreover, it is also rather standard to take an essentially quadratic truncated potential instead of the “quartic” potential . To be more precise, one considers

 (11)

for which

Let and let denote the time-step size. To start up the sequence of approximation solutions, we consider

 d0h=Ihd0 (12)

and such that

 (u0h,¯uh)+(∇p0h,¯uh) =(u0,¯uh), (13a) (∇⋅u0h,¯ph)+j(p0h,¯ph) =0, (13b)

for all . It should be noted that holds.

Then the algorithm reads as follows. Let be given. For the time step, do the following steps:

1. Find satisfying

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(dn+1h−dnhk,¯wh)+(u⋆h,(∇dnh)T¯wh)−β(u⋆⋆h,∇⋅(¯wh(dnh)T))−(1+β)(u⋆⋆⋆h,∇⋅(dnh¯wTh))+γ(wn+1h,¯wh)=0,(∇dn+1h,∇¯dh)+(˜fε(dnh)+HF2ε2(dn+1h−dnh),¯dh)−(wn+1h,¯dh)=0, (14)

for all , where (see (18) below)

 HF:=(M32+(M2−M)22)1/2

and

 u⋆h = unh+3λk(∇dnh)Twn+1h, u⋆⋆h = unh−3λβk∇⋅(wn+1h(dnh)T), u⋆⋆⋆h = unh−3λ(1+β)k∇⋅(dnh(wn+1h)T).
2. Find satisfying

 k(∇pn+1h,∇¯ph)+j(pn+1h,¯ph)−(˜un+1h,∇¯ph)=0, (15)

for all , with

 ˜un+1h = u⋆h+u⋆⋆h+u⋆⋆⋆h3 = unh+λk((∇dnh)Twn+1h−β∇⋅(wn+1h(dnh)T)−(1+β)∇⋅(dnh(wn+1h)T))

and

 j(pn+1h,¯ph)=S1ν(pn+1h−π0(pn+1h),¯ph−π0(¯ph)),

where is an algorithmic constant, and is the -orthogonal projection operator onto .

3. Find satisfying

 (un+1h−unhk,¯uh)+c(unh,un+1h,¯uh)+ν(∇un+1h,∇¯uh)+(∇pn+1h,¯uh)+(−λ(∇dnh)Twn+1h+λβ∇⋅(wn+1h(dnh)T)+λ(1+β)∇⋅(dnh(wn+1h)T),¯uh)=0, (16)

for all

To enforce the skew-symmetry of the trilinear convective term in (16), we have defined

 c(uh,vvh,wh)=((uh⋅∇)vvh,wh)+12(∇⋅uh,vvh⋅wh)

for all . Thus, for all .

Since scheme (14)-(16) is linear, it suffices to prove its uniqueness, which follows easily by comparing two solutions [6].

The main characteristic of scheme (14)-(16) is that the approximations , and are performed successively. The use of the auxiliary variable is just to be able to derive a priori energy estimates, although the computation of can be avoided as will be seen in section 4.

Concerning the stabilizing term , other choices are feasible if one wants to improve the spatial convergence rate:

 j(ph,¯ph)=Sh2ν(∇ph−π(∇ph),∇¯ph−π(∇¯ph)),

where could be the -orthogonal protection operator onto [7] or the Scott-Zhang operator into [1]. The latter is more appealing since no auxiliary variable is required to compute it.

## 3. A priori energy estimates

Let us begin by noting that the -component of the Hessian matrix of the truncated potential with respect to is given by

 Hd˜F(ε,d)ij=1ε2⎧⎪ ⎪⎨⎪ ⎪⎩2didj+(|d|2−1)δij, if |d|≤1,2didj|d|3+2|d|−1|d|δij, if |d|>1.

We thus have

 ˜F(ε,dn+1)−˜F(ε,dn)=∇d˜F(ε,dn)⋅(dn+1−dn)+12(dn+1−dn)THd˜F(ε,dn+θ)(dn+1−dn), (17)

where for some . Since is essentially quadratic, each component of the associated Hessian is uniformly bounded as

 ∥Hd˜F(ε,⋅)ij∥2L∞(IRM)≤1ε2(2+δij).

Hence, the Frobenius norm is bounded as

 (M∑i,j∥Hd˜F(ε,⋅)ij∥2L∞(IRM))1/2≤1ε2HF, where HF=(M32+(M2−M)22)1/2. (18)

In particular, using the consistence of the Frobenius norm gives

 12(dn+1−dn)THd˜F(ε,dn+θ)(dn+1−dn)≤HF2ε2|dn+1−dn|2. (19)

Consequently, by adding to a large enough first-order linear dissipation term, , we have, from (17) and (19),

 (˜f(ε,dn)+HF2ε2(dn+1−dn))⋅(dn+1−dn)≥˜F(ε,dn+1)−˜F(ε,dn). (20)

This inequality will play an essential role for the energy-stability of scheme (14). To prove this, we denote the discrete energy as

 E(uh,dh)=12∥uh∥2+λ2∥∇dh∥2+λ∫Ω˜F(ε,dh).

We are now in a position to prove the following result concerning a local-in-time discrete energy estimate.

###### Lemma 2.

Under hypotheses , it follows that, for any , and , the corresponding solution of scheme (14)-(16) satisfies the following inequality:

 E(un+1h,dn+1h)−E(unh,dnh)+k(ν∥∇un+1h∥2+λγ∥wn+1h∥2)+λ2∥∇(dn+1h−dnh)∥2+12(∥un+1h−ˆun+1h∥2+∥ˆun+1h−˜un+1h∥2)+kj(pn+1h,pn+1h)+12∥˜un+1h−u⋆h∥2+∥˜un+1h−u⋆⋆h∥2+∥˜un+1h−u⋆⋆⋆h∥23+12∥u⋆h−unh∥2+∥u⋆⋆h−unh∥2+∥u⋆⋆⋆h−unh∥23≤0. (21)
###### Proof.

We take in and in and use (20) to obtain

 λ2(∥∇dn+1h∥2−∥∇dnh∥2+˜F(ε,dn+1h)−˜F(ε,dnh))+λγk∥wn+1h∥2+λ2∥∇(dn+1h−dnh)∥2+λk((u⋆h,(∇dnh)Twn+1h)−λkβ(u⋆⋆h,∇⋅(wn+1h(dnh)T))−λk(1+β)(u⋆⋆⋆h,∇⋅(dnh(wn+1h)T)≤0. (22)

Next we take , as a test function, in (16) and introduce the auxiliary velocity to obtain

 12∥un+1h∥2−12∥ˆun+1h∥2+12∥un+1h−ˆun+1h∥2+νk∥∇un+1h∥2=0. (23)

Now we take in (15) and use the auxiliary velocity introduced previously to obtain

 12∥ˆun+1h∥2−12∥˜un+1h∥2+12∥ˆun+1h−˜un+1h∥2+kj(pn+1h,pn+1h)=0. (24)

From the definition of in (15), we write

 ˜un+1h−u⋆h3+˜un+1h−u⋆⋆h3+˜un+1h−u⋆⋆⋆h3=0.

Hence,

 12∥˜un+1h∥2−12∥u⋆h∥2+∥u⋆⋆h∥2+∥u⋆⋆⋆h∥23+12∥˜un+1h−u⋆h∥2+∥˜un+1h−u⋆⋆h∥2+∥˜un+1h−u⋆⋆⋆h∥23=0. (25)

Moreover, from the definitions of , and in , we deduce the following equalities:

 16∥u⋆h∥2−16∥unh∥2+16∥u⋆h−unh∥2−λk((∇dnh)Twn+1h,u⋆h) = 0, (26) 16∥u⋆⋆h∥2−16∥unh∥2+16∥u⋆⋆h−unh∥2+λkβ(u⋆⋆h,∇⋅(wn+1h(dnh)T)) = 0, (27) 16∥u⋆⋆⋆h∥2−16∥unh∥2+16∥u⋆⋆⋆h−unh∥2+λk(1+β)(u⋆⋆⋆h,∇⋅(dnh(wn+1h)T) = 0, (28)

 12∥un+1h∥2−12∥unh∥2+νk∥∇un+1h∥2+λk((∇dnh)Twn+1h,u⋆h)+12(∥un+1h−ˆun+1h∥2+∥ˆun+1h−˜un+1h∥2)+kj(pn+1h,pn+1h)+12∥˜un+1h−u⋆h∥2+∥˜un+1h−u⋆⋆h∥2+∥˜un+1h−u⋆⋆⋆h∥23+12∥u⋆h−unh∥2+∥u⋆⋆h−unh∥2+∥u⋆⋆⋆h−unh∥23−λk((∇dnh)Twn+1h,u⋆h)+λkβ(u⋆⋆h,∇⋅(wn+1h(dnh)T))+λk(1+β)(u⋆⋆⋆h,∇⋅(dnh(wn+1h)T)=0. (29)

Finally, we add (22) and (29); thus the terms , and cancel out, and hence (21) holds. This finishes the proof. ∎

Now, it is not difficult to extend the previous local-in-time discrete energy estimate to a global-in-time one.

###### Theorem 3.

Assume that - are satisfied. The discrete solution of scheme (14)-(16) satisfies

 maxr∈{0,⋯,N−1}{E(ur+1h,dr+1h)+kr∑n=0(ν∥∇un+1h∥2+λγ∥wn+1h∥2)}≤E(u0h,d0h). (30)
###### Proof.

The proof follows easily from Lemma 2 and by summing over . ∎

It remains to prove that the initial energy is bounded independent of .

###### Lemma 4.

Assume that hypotheses - hold. If are chosen satisfying

 hε≤K, (31)

for some constant , then