A Finite Element Approximation for Nematic Liquid Crystal Flow with Stretching Effect Based on Nonincremental Pressure-Correction Method

In this paper, a new decoupling method is proposed to solve a nematic liquid crystal flow with stretching effect. In the finite element discrete framework, the director vector is calculated by introducing a new auxiliary variable w, and the velocity vector and scalar pressure are decoupled by a nonincremental pressure-correction projection method. Then, the energy dissipation law and unconditional energy stability of the resulting system are given. Finally, some numerical examples are given to verify the effects of various parameters on the singularity annihilation, stability and accuracy in space and time.


Introduction
Liquid crystal is a kind of material with excellent properties, which has been widely used in many new advanced technical fields. For example, display devices of the electronic industry and skin cancer examination in medicine. The defects, phase transition phenomenon, molecular distribution regularity, and dynamic behavior observed in liquid crystal are related to the quality of liquid crystal equipment, which are very important issues in liquid crystal technology and have attracted extensive attention of a large number of engineers and scientists.
Generally, according to the formation conditions of liquid crystals, they can be divided into thermotropic and lyotropic. The thermotropic liquid crystal is subdivided into smectic, nematic, and cholesteric based on symmetry. Among them, the nematic liquid crystal is widely used at present. Its molecules are rod-shaped. The molecular long axes are parallel to each other but not arranged in layers. They can slide up and down, to the left and right sides, front and back, or only be parallel or nearly parallel to each other in the molecular long-axis direction; the short-range interaction between the molecules is weak. The arrangement and movement of nematic liquid crystal molecules are relatively free, and they are quite sensitive to external forces. Nematic liquid crystal is currently the main material for making liquid crystal display devices.
From a statistical average point of view, the rods locally tend to be ordered, which can be described by a unit vector that represents the average direction, that is, the director vector d d d (x). The local directivity of liquid crystal material is easily changed by the external influences, such as the molecular orientation at the interface of different materials, electric field and magnetic field, etc., which leads to a change in material properties. With the change in molecular local direction, due to the discontinuity of the molecular local direction, defects (or singularities) often appear in the material.
Physicists and mathematicians have established various mathematical models for the study of liquid crystal. In the 1960s, Ericksen [1] and Leslie [2] proposed the hydrodynamics theory of liquid crystal. The Ericksen-Leslie nematic liquid crystal model is derived from a macroscopic point of view and involves many coupling terms between two vector domains. Because the whole system is too complicated, most of the research works are based on the simplification and approximation of the model.
When the fluid is incompressible, a simplified Ericksen-Leslie model is obtained [3]. It is composed of Navier-Stokes equations coupled with anisotropic elastic stress tensor and a convective harmonic mapping heat flow equation mapped onto the sphere. The model contains linear differential constraints and nonlinear algebraic constraints. The nonlinear terms in the model bring great difficulties to the theoretical analysis and the design of numerical examples. Especially, the restriction on director vector | d d d |= 1 must be satisfied everywhere. It is very difficult to satisfy the constraint condition almost everywhere in numerical simulation. Therefore, in order to relax the constraint condition | d d d |= 1, the penalty method and the saddle point method are proposed.
Lin-Lin proposed the penalized Ericksen-Leslie model. In the literature [4], the Galerkin method was used by them to prove the local existence of the classical solution and the global existence of the weak solution for the Dirichlet problem of the two-dimensional and three-dimensional penalized Ericksen-Leslie model, and the energy estimate was also given. For the numerical solution of the penalized Ericksen-Leslie equation, ref. [5] proposed a semi-implicit, first-order linear scheme, in which the function f (d d d) was fully explicitly processed, and the scheme was conditionally stable. In [6], a nonincremental pressure-correction projection scheme [7,8] was used, a pressure stabilization term [9] was added at the same time, and finally a linear and decoupled scheme was obtained. For the penalized Ericksen-Leslie equation with stretching effect, ref. [10] proposed a corrected modified midpoint scheme by using the finite difference method to study the formation of defects at the interface of liquid crystal. Ref. [11] used finite element method to obtain spatial discretization and the time-splitting method based on a nonincremental velocity correction projection scheme [8] to decouple variables. The literature in [12] used a spectral method to study this kind of problem. At the same time, there are many papers introducing different auxiliary variables ω ω ω to analyze the model [13][14][15]. Ref. [13] presents a conditionally stable fully discrete scheme when the time step satisfies certain constraints. On this basis, the paper by [14] introduces an auxiliary variable different from [13], giving a semi-implicit Euler time discrete scheme explicitly dealing with the Ginzburg-Landau penalized function, obtaining a coupled, linear, and unconditionally stable scheme.
The saddle point method [16], through introducing the Lagrange multiplier q to the equation of director vector d, can exert spherical constraint | d d d |= 1. Based on this idea, the penalty Ericksen-Leslie model can be unified. The authors give two unconditionally stable, conserved implicit schemes and a linear semi-implicit unconditionally stable scheme, semi-implicit with respect to the nonlinear term.
Compared with some existing works, a lot of works are devoted to study the nematic liquid crystal without stretching effect, which cannot be widely used in deformable devices. The research on the liquid crystal materials with stretching effect is of great practical significance, which can be widely used in deformable devices. However, there are few studies on this model in practice.
In this paper, we turned our attention to the Ericksen-Leslie model with stretching effect. First, we use the projection scheme of nonincremental pressure correction [8] to decouple the variables. We also give a proof of the energy stability of our decoupling system and discuss the annihilation of singularities. More precisely, we give details to take care of the relation among the stabilization constant H F , the viscosity parameter ν, the geometrical parameter β, and the penalization parameter ε. Moreover, we give the evolution of the director field and velocity field of a rotating flow. Clearly, the annihilation time is much smaller than that obtained in reference [11]. Finally, the convergence of the numerical solution in time and space is analyzed. This article is arranged as follows. Some necessary notations are given in Section 2. Some necessary hypotheses are given, and the fully discrete decoupling scheme is proposed in Section 3. In Section 4, we give a proof of a priori energy estimate for the algorithm, which provides an unconditional energy stability property. Finally, a number of numerical examples are provided to demonstrate the effects of the parameters and the performance of the scheme; the numerical accuracy of the proposed system in time and space also are given in Section 5. The conclusion is reported in Section 6.

Notations and Preliminaries
In this section, we present the essential notations and preliminaries which are necessary for further consideration.

Notations
As usual, L p (Ω)(p ≥ 1) denotes the space of pth-power integrable functions defined on Ω, and the norm is · L p = ( Ω | · | p dx x x) 1/p or · L ∞ =ess sup x∈Ω | · |. If p = 2, the L 2 norm is · , and we denote the inner product in L 2 by (·, ·). For example, if For m as a non-negative integer, we denote the classical Sobolev spaces as Let C ∞ 0 be the space of infinitely times differentiable function with compact support on Ω. Then, H m 0 (Ω) is introduced as the closure of C ∞ 0 in H m (Ω). We now introduce the following function spaces in the context.

A Penalized Ericksen-Leslie Model with Stretching Effect
A general penalty version of the Ericksen-Leslie model with stretching effect to enforce the sphere constraint reads as where : is the penalty function related to the constraint | d d d |= 1; we define it by where ε > 0 is the penalty parameter.
The function f f f ε (d d d) is the gradient of the following scalar potential function F F F ε (d d d), It is easy to verify that . Furthermore, the parameters γ, ν, and λ > 0 represent relaxation time, viscosity, and elasticity, respectively. Additionally, parameter β ∈ [−1, 0] is a constant that determines the geometry of the molecule. For instance, when β = −1, −1/2, 0, the molecules are rod-shaped, spherical, and disk-shaped [10,11,18], respectively.
To system (1)-(3), we add homogeneous Dirichlet conditions for the velocity field and homogeneous Neumann boundary conditions for the director field, and the initial conditions In order to better understand our proposed decoupling scheme, we hereby give an energy law for system (1)-(3), which is true under some regularity assumptions for d d d and u u u (see [10,13,19] for details). First, we give equations Then, taking the inner product of (1) and (2) with here, 1 2 u u u 2 represents the kinetic energy, λ 2 ∇d d d 2 represents the elastic energy, represents the penalty energy, and obviously, the total energy is E(u u u, d d d).

Hypotheses
We introduce the hypotheses that are required in the following.
(H2) Let K represent any subregion after dividing Ω into finite subregions, also called element domain. Additionally, its a bounded closed set with nonempty interior and piecewise smooth boundary. We use {T h } h>0 for all of K, so Ω = K∈T h K. In general, K is a triangle or quadrilateral in a two-dimensional space and a tetrahedron or hexahedron in a three-dimensional space.
(H5) Let P 1 (K) denote the set of linear polynomials on K. Under Hypotheses 4, the space of continuous, piecewise polynomial functions associated to T h are denoted as follows: and denote the space of piecewise constant function as The triangulation of Ω and the discrete spaces satisfy ( [5,6,11]): (a.) The approximation properties: (b.) The stability properties: where I I I h is an interpolation operator into D D D h , and C 1 , C 2 , and C 3 are constant independents of h.
In this paper, we, respectively, choose D h = X h , V h = X h H 1 0 (Ω), and P h = X h L 2 0 (Ω) as the approximate spaces of direction, velocity, and pressure. We choose discontinuous finite element space W h = Y h as the approximate space of an auxiliary variable. The finite element spaces that we choose for velocity and pressure do not satisfy the discrete inf-sup condition for α > 0 independent of h.

The Fully Discrete Scheme
Here, we use the nonincremental pressure-correction method to obtain the fully discrete scheme of system (1)- (3).
The stabilization term j(p, q) can be defined as where S is an algorithmic constant and Π h is a standard L 2 -projection operator onto Y h . This stabilization term is also used in [20].
Step(n + 1) We discretize system (1) be given. For n + 1, n ≥ 0, perform the following steps: herein, and with M being the space dimension ( [6]). (2) By using the nonincremental pressure-correction method for all We set the trilinear convective term has skew-symmetry, so c(u u u n h , u u u n+1 h , u u u n+1 h ) = 0.
Since scheme (14)- (19) is linear, we can easily prove its existence and the uniqueness of the solution by using the same techniques as in [6].

Energy Estimate
In this section, we give the proof about energy estimates for schemes (14)- (19). First of all, we give an inequality in the following, the proof of which already exists in [6].
Theorem 1. Assume that the hypotheses in Section 3.1 are satisfied. In addition, let Then, the numerical solutions (d d d n+1 h , ω ω ω n+1 h , u u u n+1 h , and p n+1 h ) of the schemes (14)- (19) have the property (14), and take into account (20), obtaining Moreover, selecting u u u h = ∆tu u u n+1 h in (17) and defining u u u h = u u u h1 +u u u h2 +u u u h3

3
, it follows that Next, choosing q h = p n+1 h in (19) and by using the equality (21), we have therefore, taking the inner product of both sides of (21) with itself, we obtain Hence, we obtain, by adding Equations (24) and (26), In addition, taking the inner product of both sides of u u u h1 , u u u h2 , and u u u h3 in (15) with u u u h1 , u u u h2 , and u u u h3 , respectively, we have From the definition of u u u h , we have therefore, Adding (23), (27)-(30), and (32) implies that which implies the assertion (22).
Theorem 1 is a result concerning a local discrete energy estimate. Then, we give a global energy estimate about time for schemes (12)- (19) in the following theorem.
where (u u u 0 h andd d d 0 h ) are defined in (12). Additionally, if for some constant K > 0, (h, ε) satisfy h/ε ≤ K, then there exists a constant C 0 > 0 such that Proof. The proof of (34) follows easily from Theorem 1 and by summing over n. Using the same method as Lemma 4.4 in [6], we can prove that E(u u u 0 h , d d d 0 h ) ≤ C 0 . Then, we can easily obtain (35).

Numerical Experiments
In this section, we present a numerical example to carry out a sensitivity study of schemes (12)- (19). More precisely, we give details to take care of the relation among the stabilization constant H F defined in (16), the viscosity parameter ν, the geometrical parameter β, and the penalization parameter ε. Then, we consider a rotating flow and give the evolution of the director field and velocity field for the annihilation of two and four singularities. At the end of this paper, we investigate the numerical accuracy of the proposed system in space and time.
The numerical solutions are implemented by FreeFem++ [21] and Matlab.

Annihilation of Singularities
We consider the initial conditions of the nematic liquid crystal with stretching effect (1)-(3) are u u u 0 = 0 0 0, where d d d = (x 2 + y 2 − 0.25, y). This example was also used to research nematic liquid crystal in [5,6,10,11,13,22]. In this experiment, we choose computational domain Ω = (−1, 1) × (−1, 1), time step size ∆t = 0.001, and use a 32 × 32 grid in the computation. First of all, let β = −1 and ε = 0.05. We research the stability , annihilation time, and energies of two singularities at M = 0, 1, 2, 3 under different viscosity coefficient ν; see Table 1. We are concerned with the dependence of stability on the parameters M and ν. In particular, we present snapshots of the director and velocity fields at M = 2 and the time T = 1.0 in Figure 1. As can be seen from Figure 1, when the numerical solution is in a stable state, different ν values have no great influence on the director field,where the stability refers to the ability of fluid motion of a certain form to recover its original form after initial disturbance. However, it is obvious that the velocity field is somewhat chaotic when ν = 0.001, while the trend of the velocity field becomes regular when ν gradually increases. Especially at ν = 1.0, the velocity field has largely quieted down. A, we also show the behavior of the energies in Figure 2. We can find that with the annihilation of the singularity, the energy decreases rapidly, which is consistent with the results obtained in [10].  Secondly, we investigate the dependence of stability on parameters H F and β. We take ν = 1.0, and ε = 0.05 and vary β = 0, −0.5, −1, and M = 0, 1, 2, 3. The results are shown in Table 2. In order to more intuitively observe the influence of different β values on the annihilation, we specially give snapshots of the director field in the initial state and the stable state in Figure 3. We do not find any significant difference between these figures with these different β values. Then, we investigate the dependence of stability on parameters H F and ε. We take ν = 1.0 and β = −1 and vary ε = 0.1, 0.05, 0.01, and 0.001 and M = 0, 1, 2, 3. The results are shown in Table 3. The results are similar to those presented in Table 4 of [11]. It is found that when ε = 0.1, and 0.05, the annihilation time becomes longer and longer with the increase in M, and the maximum value of the corresponding kinetic energy becomes smaller and smaller. Even when ε = 0.01 and 0.001, where there is no longer annihilation, this rule can still be followed.  Figure 4 shows that evolution in time of the kinetic energy for different ε values when M = 2. By observation, we find that the kinetic energy behavior is obviously different under the conditions that singularities can annihilate (a) or not annihilate (b). As noted in references [6,11], a possible explanation for this behavior is that the velocity field generated by the elastic tensor is insufficient to move the singularity though the convective term in the director equation. Additionally, especially at ε = 0.01, and 0.001, the kinetic energy goes down to zero at the beginning. One might think that if the kinetic energy associated with the velocity field is large enough to move the singularities, then they will move towards each other until annihilation.  Furthermore, to visualize the difference in energy in the stable state and unstable state, we present Figure 5. In addition, in Figure 6, we show the evolution in time of energy for the unstable state. Obviously, they are completely different from the evolution of the energy for the stable state. The snapshots of the unstable state (ε = 0.001, and M = 0) and the no annihilation (ε = 0.001, and M = 1) at t = 0.001 and t = 1.0 are given in Figures 7 and 8, respectively.  Lastly, we give the snapshots of the director and the velocity field displayed at times t = 0.1, 0.2, 0.3, and 1.0 in Figure 9. Here, we choose the parameters ε = 0.1, β = −1, ν = 1.0, M = 2, λ = 1.0, γ = 1.0, and ∆t = 0.001.

The Behavior under Rotating Flow
In this example, we consider a rotating flow in a square domain Ω = (−1, 1) × (−1, 1). The initial director field of two singularities is the same as Example 5.1. The initial director field of four singularities is the same as Section 4.2.2 in [11]. The initial velocity field is u u u 0 = (−20y, 20x). We present the snapshots of the director field for the annihilation of two singularities at times t = 0.001, 0.4, 0.6, and 1.0 and four singularities at times t = 0.001, 0.1, 0.2, and 1.0 (see Figures 10 and 11). Here, we choose the parameters ε = 0.05, β = −1, ν = 1.0, M = 2, λ = 1.0, and γ = 1.0. Additionally, ∆t = 0.001. Clearly, the annihilation times are around t = 0.555 and t = 0.160, respectively. They are all smaller than those obtained in [11]. Figure 10. The snapshots of the director field of two singularities for the rotating flow at times t = 0.001, 0.4, 0.6, and 1.0. Figure 11. The snapshots of the director field of four singularities for the rotating flow at times t = 0.001, 0.1, 0.2, and 1.0.
Here, v v v can be u u u, d d d, and p. In Figure 12 and Table 4, we present the time errors and convergence rates for the director, velocity, and pressure measured in the L 2 − norm and H 1 − norm at the final time T = 0.1, respectively. The time error for the director, velocity, and pressure in the L 2 − norm and H 1 − norm are of O(∆t), respectively. Here, we run the code with the spatial mesh size: 64 × 64 and time steps ∆t = 10 −3 , 5 × 10 −4 , 2.5 × 10 −4 , 1.25 × 10 −4 , and 6.25 × 10 −5 , respectively. Table 4. Time convergence rates for the director, velocity, and pressure.  In Figure 13 and Table 5, we present the space errors and space convergence rates for the director, velocity, and pressure measured in the L 2 − norm and H 1 − norm at the final time T = 0.1. The error for the director, velocity, and pressure in the L 2 − norm and H 1 − norm are of O(h 2 ) and O(h), respectively. Compared with the result obtained in [6], ours is obviously better. Here, we run the code with the spatial mesh sizes: 16 × 16, 32 × 32, 64 × 64, 128 × 128, and 256 × 256 and time step ∆t = 0.001, respectively. Remark 1. Compared with [6], which used the same method as this paper to research the Ericksen-Leslie model without stretching effect, the time and space convergence rates obtained in this paper are better.

Conclusions
In this paper, a new decoupling scheme based on the nonincremental pressurecorrection projection method was proposed to approximate the penalized Ericksen-Leslie model with stretching effect. The scheme is a linear and unconditionally stable system. One bright spot is that equal low-order finite element spaces were used for our scheme. That is, P1 − P1 − P1 finite element spaces were used for the director, velocity, and pressure.
In our numerical experiment, the sensitivity of the viscosity parameter ν, the geometrical parameter β, the penalization parameter ε, and the stabilization constant H F in the proposed scheme were studied. It was found that the sensitivity of singularity annihilation to parameters β, ε, and H F is very consistent with the result in [11]. Notably, we studied the effect of different viscosity coefficients on the annihilation and found that our method works for different ν (related to the Reynolds number). However, when ν = 0.001, the velocity field looks messy near the singularity region. We will consider how to overcome this problem in future studies. In addition, experiments on the annihilation of two singularities and four singularities in a rotating flow were performed. What is gratifying is that by comparing with the existing methods, the results we obtain are better than those in [11]. Furthermore, we verified the numerical accuracy in time and space.