A projectionbased timesplitting algorithm for approximating nematic liquid crystal flows with stretching
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 timesplitting algorithm to uncouple the calculation of the variables: the velocity and pressure are computed by using a projectionbased 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 equalorder interpolation for the velocity and the pressure. Numerical experiments concerning annihilation of singularities are presented to show the stability and efficiency of the scheme.
Mathematics Subject Classification: Nematic liquid crystal; Finite elements; Projection method; Timesplitting 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 wellknown 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 diskshaped molecules may exhibit strikingly different properties from those composed of rodshape 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 socalled nematic phase. In it, the molecules have no positional order, but they selforganize to have longrange directional order. Thus, the molecules flow freely and their center of mass positions are randomly distributed as in a liquid, while maintain their longrange directional order. The breakdown of the selfalignment 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 EricksenLeslietype 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):
(1) 
We complete this system with homogeneous Dirichlet conditions for the velocity field (nonslip) and homogeneous Neumann boundary conditions for the director field:
(2) 
and the initial conditions
(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, rodlike and disklike 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
(4) 
It should be noted that , for all , for the following scalar potential function
The following energy law for system (1) holds under some regularity assumptions for and [17, 18]:
(5) 
where
(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 NavierStokes 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 timeconsuming 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, projectionbased timesplitting 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, infsup 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 CrankNicolson 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 semiimplicit BDF2 rotational pressurecorrection method for an axisymmetric domain for studying the annihilation of a hedgehogantihedgehog 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 abovementioned techniques for the NavierStokes 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 loworder finiteelements, 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 energystability 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 thpower integrable realvalued 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
and its norm is simply denoted by . For a nonnegative integer, we denote the classical Sobolev spaces as
associated to the norm
where is a multiindex 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:
We denote by and , the closures of , in the  and norm, respectively, which are characterized (for being Lipschitzcontinuous) by (see [22])
where is the outward normal to on . Finally, we consider
2.2. Hypotheses
Herein we introduce the hypotheses that will be required along this work.

Let be a bounded domain of with a polygonal or polyhedral Lipschitzcontinuous boundary.

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

Assume three sequences of finitedimensional spaces , and associated with such that , and . Also, consider an extra finiteelement space with .

Suppose that with a.e. in .
In particular, hypothesis allows us to consider equalorder finiteelement 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
and the set of piecewise constant functions as
We choose the following continuous finiteelement spaces
for approximating the director, the velocity and the pressure, respectively. Additionally, we select the extra discontinuous finiteelement to be the space for an auxiliary variable related to the vector director.
Observe that our choice of the finiteelement spaces for velocity and pressure does not satisfy the discrete infsup condition
(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 ScottZhang interpolation operator, see [21].
Proposition 1.
Assuming hypotheses , there exists an interpolation operator satisfying
(8) 
(9)  
(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 timesplitting be very appealing for high performance computing.
The starting point to design our timesplitting method is the nonincremental velocitycorrection method for the NavierStokes 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 timestep size. To start up the sequence of approximation solutions, we consider
(12) 
and such that
(13a)  
(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:

Find satisfying
(15) for all , with
and
where is an algorithmic constant, and is the orthogonal projection operator onto .

Find satisfying
(16) for all
To enforce the skewsymmetry of the trilinear convective term in (16), we have defined
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].
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
We thus have
(17) 
where for some . Since is essentially quadratic, each component of the associated Hessian is uniformly bounded as
Hence, the Frobenius norm is bounded as
(18) 
In particular, using the consistence of the Frobenius norm gives
(19) 
Consequently, by adding to a large enough firstorder linear dissipation term, , we have, from (17) and (19),
(20) 
This inequality will play an essential role for the energystability of scheme (14). To prove this, we denote the discrete energy as
We are now in a position to prove the following result concerning a localintime discrete energy estimate.
Lemma 2.
Proof.
We take in and in and use (20) to obtain
(22) 
Next we take , as a test function, in (16) and introduce the auxiliary velocity to obtain
(23) 
Now we take in (15) and use the auxiliary velocity introduced previously to obtain
(24) 
From the definition of in (15), we write
Hence,
(25) 
Moreover, from the definitions of , and in , we deduce the following equalities:
(26)  
(27)  
(28) 
Adding equalities (23)(28) leads to
(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 localintime discrete energy estimate to a globalintime one.
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
(31) 
for some constant , then