A second-order numerical scheme for the Ericksen-Leslie equation

: In this paper, we consider a ﬁnite element approximation for the Ericksen-Leslie model of nematic liquid crystal. Based on a saddle-point formulation of the director vector, a second-order backward di ﬀ erentiation formula (BDF) numerical scheme is proposed, where a pressure-correction strategy is used to decouple the computation of the pressure from that of the velocity. Designing this scheme leads to solving a linear system at each time step. Furthermore, via implementing rigorous theoretical analysis, we prove that the proposed scheme enjoys the energy dissipation law. Some numerical simulations are also performed to demonstrate the accuracy of the proposed scheme.


Introduction
The study of liquid crystal has aroused an increasing interest in biology, physics and engineering owing to their optical properties.There are different theories to describe the nematic liquid crystal, including Doi-Onsager theory, Landau-de Gennes theory [1,2], and Ericksen-Leslie theory [3][4][5].
Here, we are interested in the Ericksen-Leslie model that simulates the motion of point defects.Let Ω ∈ R d (d = 2, 3) be a bounded and open domain with smooth boundary ∂Ω.Ω T = Ω × [0, T ] is used for all the following equations, where x ∈ Ω, t ∈ [0, T ].Formulated by a unit director d, the Ericksen-Leslie model is given as follows: For the Ericksen-Leslie model, some analysis results are pointed out in the literature on the existence, uniqueness, regularity, and long time asymptotic behavior of the solution [6][7][8][9][10][11][12][13][14].There also exist abundant works on the numerical methods.Spectral method has been discussed in [15].A linear fully discrete method and a semi-implicit Euler method are considered for solving a penalized nematic liquid crystal model in [16].A linearized semi-implicit Euler finite element method is proposed in [17], where the temporal and spatial errors are shown by using an error splitting technique.Two fully discrete finite element methods for the limiting Ericksen-Leslie model are presented in [18].In addition, the concave-convex decomposition method of the corresponding energy function is considered in [19].An unconditional energy stable time-splitting finite element method for approximating the Ericksen-Leslie equations is proposed in [20].
In numerical analysis, we notice that the current numerical analysis are based on the assumption that the primary component for the time invariant derivative consisting of the term −W Here, p is the hydrostatic pressure, the relaxation time γ and the fluid viscosity η are the physical positive constants.The first equation of (1.3) is the Navier-Stokes equation related to the conservation of the linear momentum.The second equation of (1.3) models the conservation of the angular momentum.The last equation of (1.3) represents the incompressibility of the liquid.The boundary conditions are where n presents the outward normal vector on the boundary.The initial conditions are Derived from a variational approach coupled with Onsager energy [21,22], the total energy E(u, d) is the sum of the kinetic energy E k , the elastic energy E e and the penalty energy E p [23] The penalized version of Ericksen-Leslie model has been extensively investigated.In [24], Lin introduces an unconditionally stable, nonlinear scheme for the penalized Ericksen-Leslie model.In [25], based on continuous finite element in space and a semi-implicit integration in time, a fully discrete scheme is studied to approximate the Ericksen-Leslie model by means of a Ginzburg-Landau penalized problem.Based on a saddle-point structure of the director vector, Santiago Badia designs a linear semiimplicit algorithm in [26].It is noticed that the above papers are first-order numerical schemes in time.There are few works on the development of high order and energy stable numerical schemes for the Ericksen-Leslie model.
In recent years, the second-order accurate numerical schemes have attracted more attention, due to the great advantage over their first-order scheme with regard to numerical accuracy.It is wellknown that the discussion of second-order scheme is more difficult than the first-order one.For these reasons, second-order energy stable schemes have been highly desirable.The BDF discrete scheme approximates each term at t n+1 .A fully nonlinear scheme is obtained by applying the BDF method directly.To overcome this difficulty, we propose a semi-implicit scheme, in which the nonlinear terms are approximated by the second-order extrapolation formula.Such a numerical algorithm leads to a linear system of equations to solve at each time step.The second-order backward differentiation formula (BDF) approximation has also been used in many time-dependent problems, such as the Landau-Lifshitz equation [27,28], Allen-Cahn equation [29], Cahn-Hilliard equation [30][31][32], and Cahn-Hilliard-Hele-Shaw equation [33].But it is still an open question for the Ericksen-Leslie model.Furthermore, the BDF method can also be used to study magneto-hydrodynamic equations in the future [34].
In this paper, we propose a second-order scheme to approximate the time derivative for the Ericksen-Leslie equation.Furthermore, a pressure-correction strategy is used to decouple the computation of the pressure from that of the velocity [35,36].Via implementing rigorous theoretical analysis, we prove that the proposed scheme enjoys the unconditionally energy stability.Finally, some numerical simulations are also performed to demonstrate the accuracy of the proposed scheme.
The outline of this article is organized as follows.In Section 2, we give the discrete finite element scheme; The unconditional stability of proposed numerical scheme is proven rigorously in Section 3; In Section 4, some numerical simulations are performed to show the accuracy of the scheme; In Section 5, conclusions are drawn.

The semi-discrete scheme
In this section, some notations are given in order to prove the following conclusions.The inner product and its associating norm are denoted by (•, •) and . Then, we introduce the following spaces: We use a second-order BDF scheme for n ≥ 1, and use a first-order scheme for n = 1.The semi-discrete numerical scheme of the penalized Ericksen-Leslie equations can be written as follows.

Fully discrete scheme
Let Π h be a set of triangulations of Ω with Ω = ∪ k∈Π h k, where h → 0 is assumed to be uniformly regular.Here h = sup k∈Π h diam (k).We denote the finite element spaces Pressure stability relies on the inf-sup condition.There exists β > 0 with no restriction of the mesh grid size h such that inf The strategy of pressure-correction is time-marching method which is composed of two steps.In the first step, the pressure is treated explicitly.In the second step, the pressure is projected by the former velocity onto the space X h .For n ≥ 1, given (2.6)

The analysis of stability
In this section, we prove that the fully discrete scheme is unconditionally energy stable.

Proof. By denoting auxiliary variable
Choosing 3), and taking , we arrive at Using u h = ũn+1 h , we gain

.6)
Setting e h = 1 2 q n+1 h , we easily get According to (2.5), we easily obtain ũn+1 To approach the pressure term, we rewrite (2.5) as By squaring both sides of the above equation, we can derive Then, we have 3 4τ Noticing the fact we deduce (3.12) Similarly, we get (3.14) Combining (3.5)-(3.7)and (3.11), then multiplying both sides by 2τ, we conclude (3.15) Finally, we can deduce the discrete energy law.The proof is complete.
Corollary 3.1.The following estimates hold for some constant C > 0 independent of τ Proof.Summing the discrete energy inequality (3.1) from n = 1 to N − 1, we obtain (3.17) The proof is complete.
The above second order scheme consists of three temporal levels n + 1, n, n − 1.Thus, n ≥ 1.Then, we have to start with the initial values.Theorem 3.2.For all τ > 0, the numerical scheme (2.6) is unconditionally energy stable and satisfies the following discrete energy law where Proof.By defining auxiliary variable the following identity holds Setting By replacing u h = u 1 h , we get Noticing the fact Similarly, we easily have (3.28) Combining (3.22)-(3.24),then multiplying τ at two sides, we obtain Finally, we can deduce the discrete energy law.The proof is complete.

Numerical simulations
In this section, we perform some numerical simulations to show the accuracy and stability for the proposed scheme.All numerical simulations are carried out by using the Freefem++ package [37].

The order of convergence
In this subsection, we perform numerical simulations about the spacial and temporal convergence rates.We use the differences between the solutions on a coarse mesh and a fine mesh to calculate the error , and obtain the convergence rate log 2 (ξ h u /ξ , and fix the parameter T = 1.
We compute the reference solution.The initial conditions are Table 1 shows the spacial convergence rates of ξ h u H 1 and ξ h d H 1 .Parameters are set as γ = 0.35, η = 0.7, ε = 0.2, h = 0.25, 0.125, 0.0625, 0.03125, 0.015625.The spacial convergence orders computed from the errors ξ h u H 1 and ξ h d H 1 are close to 2. We get the same conclusion in comparison to variable in Tables 2-3, where = 0.23, 0.25.
Table 4 shows the temporal convergence rates of ξ h u H 1 and ξ h d H 1 .Parameters are set as γ = 0.11, η = 0.82, = 0.197, τ = 0.25, 0.125, 0.0625, 0.03125.The temporal convergence orders computed from the errors ξ h u H 1 and ξ h d H 1 are close to 2. We get the same conclusion in comparison to variable in Tables 5 and 6, where = 0.198, 0.2.These numerical results show that our numerical algorithm is correct.

The dissipation of energy
In this part, we perform the energy decay for the proposed scheme.The kinetic energy and elastic energy of the proposed scheme can be written as The problem Ericksen-Leslie hydrodynamic model is in domain [−1, 1] × [−1, 1], and sets the parameters T = 0.5, γ = η = 1.The initial conditions are u 0 = 0, d = (sin(2π(cos x − sin y)), cos(2π(cos x − sin y))).
Figure 1 describes the evolution of elasticity and kinetic energy under different values of variable ε, where ε = 0.09, 0.1, 0.15, 0.2.It indicates that the proposed scheme is unconditionally energy stable.Furthermore, we find that the energies are not sensitive to ε, The velocity is dissipated, and the elastic energy shows a linear and slow decay.
We consider problem with γ = 1, 0.5 and n = 16, 32, which are shown in Figure 2. The elastic energy and kinetic energy are sensitive to γ.With regard to η and τ, we also plot the results for different values of ε d (t) and ε u (t) in Figure 3.The decay curves of elastic energy is almost identical in all cases.As can be seen from the kinetic energy curve, the choice of parameters η and γ has a great influence on the results.

Annihilation and stable defects
The numerical experiments are associated with the annihilation of singularities.We will consider two numerical simulations consisting of the motions of two singularities.We perform the evolutions of director field d over time under certain conditions.
The director field has an influence on the stress tensors that can govern the velocity fluid and the disappearance of singularity.Thus we can conclude that the director field plays an important role in the annihilation of singularities.For the different initial values of d, we give the numerical examples which are related to the annihilation of the singularities.We perform the evolution of the singularities until the simulation reaches the steady state.

Example 1:
In the test, we will consider numerical experience consisting of the motion of two singularities.We  The evolution of singularities is depicted in Figure 4. We present snapshots of the director field d displayed at times T = 0.001, 0.005, 0.1, 3. It is observed that two singularities would move towards the origin, meet and annihilation.We observe that the energy starts to have a significant decay with the annihilation.The energy has no significant change after reaching 1.We find that the energy reaches the steady state.Furthermore, the annihilation happens later.
In Figure 4, we plot the director field at four different times.The singularities devote to moving closer with the flow at T = 0.005.The singularities just approach to annihilation at T = 0.1.Finally, a steady state is reached at T = 3.

Example 2:
Changing initial directors, we give the evolution of singularities.A comparison of the annihilation for two different initial directors can be found.It is computed in the unit circle (x, y) : x 2 + y 2 < 1.The parameters are chosen as γ = η = 1.The initial conditions are u 0 = 0, d = (sin(2π(cos x − sin y)), cos(2π(cos x − sin y))).singularities annihilate in a short time.Two singularities annihilate roughly after T = 3.We also find a little difference about the molecule orientation near the boundary.Figure 5 shows the simulation results.Therefore, these numerical results are performed as expected.

Conclusions
In this paper, a second order BDF numerical scheme with Lagrange multiplier for the Ericksen-Leslie equation is presented.In addition, a pressure-correction strategy is used to decouple the computation of the pressure from that of the velocity.With a plenty of powerful proofs and calculations, we show that the proposed scheme is unconditionally stable in energy.Furthermore, the convergence rates of relative error are close to 2 in time and space.The motions of singularities are simulated with different director values.The results are performed as expected.

(1. 1 )∂
where u describes the velocity field, d represents the molecular orientation, h describes the molecular field, W = 1 2 ∂ β u α − ∂ α u β presents the vorticity tensor, D = 1 2 ∂ β u α + ∂ α u β denotes the rate of strain sensor, a is a geometry parameter of liquid crystal molecules, − a 2 (dh + hd) + 1 2 (dh − hd) means the elastic stress tensor associated with liquid crystal dynamics.Some operators are listed as ∇d = ∂ j d i , ii d and u • ∇d = M i=1 u i ∂ i d.

From a numerical
point of view, the nonconvex constraint | d |= 1 is difficult to achieve at the discrete level.Then, a penalized version is usually considered, where the constraint | d |= 1 is weakly performed by adding the Ginzburg-Landau function f (d) = 1 2 | d | 2 −1 d, where > 0 is a positive constant that dictates the interface width.It is convenient to introduce the Lagrange multiplier q =| d | 2 −1, then the penalty function is rewritten by f (d) = 1 2 qd.Therefore, based on a saddle-point structure of the director vector d, the penalized version of Ericksen-Leslie model reads

Figure 1 . 2 .εFigure 3 .
Figure 1.The elastic energy ε d (t) and the kinetic energy ε u (t).The results are shown for different values of .

Table 1 .
The spacial convergence rate.

Table 2 .
The spacial convergence rate.

Table 3 .
The spacial convergence rate.

Table 4 .
The temporal convergence rate.

Table 5 .
The temporal convergence rate.

Table 6 .
The temporal convergence rate.