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

1.   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 $\Omega\in R^{d}$ ($d = 2, 3$) be a bounded and open domain with smooth boundary $\partial\Omega$. $\Omega_{T} = \Omega\times [0, T]$ is used for all the following equations, where $x\in \Omega, \; t\in [0, T]$. Formulated by a unit director $\mathbf{d}$, the Ericksen-Leslie model is given as follows:

where $\mathbf{u}$ describes the velocity field, $\mathbf{d}$ represents the molecular orientation, $\mathbf{h}$ describes the molecular field, $\mathbf{W} = \dfrac{1}{2}\left(\partial_{\beta}\mathbf{u}_{\alpha}-\partial_{\alpha}\mathbf{u}_{\beta}\right)$ presents the vorticity tensor, $\mathbf{D} = \dfrac{1}{2}\left(\partial_{\beta}\mathbf{u}_{\alpha}+\partial_{\alpha}\mathbf{u}_{\beta}\right)$ denotes the rate of strain sensor, $a$ is a geometry parameter of liquid crystal molecules, $-\dfrac{a}{2}\left(\mathbf{d}\mathbf{h}+\mathbf{h}\mathbf{d}\right)+\dfrac{1}{2}\left(\mathbf{d}\mathbf{h}-\mathbf{h}\mathbf{d}\right)$ means the elastic stress tensor associated with liquid crystal dynamics. Some operators are listed as $\nabla\mathbf{d} = \partial_{j}d_{i}$, $\triangle\mathbf{d} = \sum\limits_{i = 1}^{M}\partial_{ii}\mathbf{d}$ and $\mathbf{u}\cdot\nabla\mathbf{d} = \sum\limits_{i = 1}^{M}\mathbf{u}_{i}\partial_{i}\mathbf{d}$.

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].

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

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 $-\mathbf{W}\cdot\mathbf{d}-a\mathbf{D}\cdot\mathbf{d}$ and $-\dfrac{a}{2}\left(\mathbf{d}\mathbf{h}+\mathbf{h}\mathbf{d}\right)+\dfrac{1}{2}\left(\mathbf{d}\mathbf{h}-\mathbf{h}\mathbf{d}\right)$ are neglected due to the complexity of structures. It leads to a much simpler version. The model (1.2) is reduced to

Here, $p$ is the hydrostatic pressure, the relaxation time $\gamma$ and the fluid viscosity $\eta$ 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 $\mathbf{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(\mathbf{u}, \mathbf{d})$ is the sum of the kinetic energy $E_{k}$, the elastic energy $E_{e}$ and the penalty energy $\mathbf{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 semi-implicit 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 well-known 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.

2.   The discrete scheme

2.1.   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 $\left(\cdot, \cdot\right)$ and $\|\cdot\|$ in $L^{2}$, respectively. $\left(u, v\right) = \int_{\Omega }u\left(x\right)v\left(x\right)dx$, $\|u\| = \|u\|_{L^{2}}$, $\|u\|_{H^{m}} = {\left(\sum\limits_{\mid\alpha\mid\leq m}\|D^{\alpha}u\|^{2}\right)}^{\frac{1}{2}}$.

Then, we introduce the following spaces: $\mathbf{Y} = H^{1}\left(\Omega\right)^{d},$ $\mathbf{Z} = L^{2}\left(\Omega\right)^{d}, \; \mathbf{X}$ $= H^{1}_{0}\left(\Omega\right)^{d} =$ $\{v\in H^{1}\left(\Omega\right)^{d}: v_{\partial{\Omega}} = 0\},$ $R = L^{2}_{0}\left(\Omega\right)$ $= \{r\in L^{2}\left(\Omega\right):$ $\int_{\Omega}rdx = 0\}.$

Let $N$ be a positive integer and $0 = t_{0} < t_{1} < \cdot\cdot\cdot < t_{N} = T$ be a uniform partition of $[0, T]$, where $t_{i} = i\tau$, $\tau = t_{i}-t_{i-1}$. Let $\mathbf{u}^{n} = \mathbf{u}\left(\cdot, t_{n}\right)$, $\mathbf{d}^{n} = \mathbf{d}\left(\cdot, t_{n}\right)$, $q^{n} = q\left(\cdot, t_{n}\right)$, $p^{n} = p\left(\cdot, t_{n}\right)$. We use a second-order BDF scheme for $n\geq1$, 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. For $n\geq1$, $(\overline{\mathbf{u}}, \overline{\mathbf{d}}, \mathbf{v}, e, \mathbf{s}, b)\in \mathbf{X}\times \mathbf{Z}\times \mathbf{X}\times Z\times \mathbf{Z}\times R$, given $\left(\mathbf{d}^{n-1}, \mathbf{u}^{n-1}, q^{n-1}, \mathbf{d}^{n}, \mathbf{u}^{n}, q^{n}\right)$, find $\left(\tilde{\mathbf{u}}^{n+1}, \mathbf{d}^{n+1}, q^{n+1}, \mathbf{h}^{n+1}\right)\in \left(\mathbf{X}\times\mathbf{Y}\times Z\times \mathbf{Z}\right)$, such that

Step 1:

Find $(\mathbf{u}^{n+1}, p^{n+1})$ as the solution of

Step 2:

For $n = 1$, given $\left(\mathbf{u}^{0}, \mathbf{d}^{0}, q^{0}\right)$, find $\left(\mathbf{u}^{1}, \mathbf{d}^{1}, q^{1}, \mathbf{h}^{1}\right)\in \left(\mathbf{X}\times\mathbf{Y}\times Z\times\mathbf{Z}\right)$ as the solution of

2.2.   Fully discrete scheme

Let $\Pi_{h}$ be a set of triangulations of $\Omega$ with $\overline\Omega = \cup_{k\in\Pi_{h}}k$, where $h\rightarrow0$ is assumed to be uniformly regular. Here $h = \sup\limits_{k\in\Pi_{h}}diam\left(k\right)$. We denote the finite element spaces $\mathbf{X}_{h}\in \mathbf{X}, \mathbf{Y}_{h}\in \mathbf{Y}, \mathbf{Z}_{h}\in \mathbf{Z}, R_{h}\in R$.

Pressure stability relies on the inf-sup condition. There exists $\beta > 0$ with no restriction of the mesh grid size $h$ such that $\inf\limits_{r_{h}\in R_{h}}\sup\limits_{v_{h}\in \mathbf{X}_{h}}\dfrac{\left(r_{h}, \nabla\cdot v_{h}\right)}{\|r_{h}\|_{0}\|v_{h}\|_{1}}\geq\beta$. 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 $\mathbf{X}_{h}$. For $n\geq1$, given $\left(\mathbf{u}_{h}^{n-1}, \mathbf{d}_{h}^{n-1}, q_{h}^{n-1}, \mathbf{u}_{h}^{n}, \mathbf{d}_{h}^{n}, q_{h}^{n}\right)$, find $\left(\mathbf{d}_{h}^{n+1}, \tilde{\mathbf{u}}_{h}^{n+1}, q_{h}^{n+1}, \mathbf{h}_{h}^{n+1}\right)\in \left(\mathbf{Y}_{h}\times \mathbf{X}_{h}\times Z_{h}\times \mathbf{Z}_{h}\right)$ as the solution of

Step 1:

Find $\mathbf{u}_{h}^{n+1}, p_{h}^{n+1}$ as the solution of

Step 2:

For n = 1, given $(\mathbf{u}_{h}^{0}, \mathbf{d}_{h}^{0}, q_{h}^{0})$, find $(\mathbf{u}_{h}^{1}, \mathbf{d}_{h}^{1}, q_{h}^{1}, \mathbf{h}_{h}^{1})\in \left(\mathbf{X}_{h}\times\mathbf{Y}_{h}\times Z_{h}\times\mathbf{Z}_{h}\right)$ as the solution of

3.   The analysis of stability

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

Theorem 3.1. For all $\tau > 0$ and $1\leq n\leq N-1$, the numerical scheme (2.4)-(2.5) is unconditionally energy stable and satisfies the following discrete energy law

where

Proof. By denoting auxiliary variable

the following identity holds

Choosing $\overline{\mathbf{d}}_{h} = \dfrac{3\mathbf{d}_{h}^{n+1}-4\mathbf{d}_{h}^{n}+\mathbf{d}_{h}^{n-1}}{2\tau}$, we have

Form (3.2)-(3.3), and taking $\overline{\mathbf{d}}_{h} = \tilde{\mathbf{u}}_{h}^{n+1}\cdot\nabla\left(2\mathbf{d}_{h}^{n}-\mathbf{d}_{h}^{n-1}\right)$, we arrive at

Using $\overline{\mathbf{u}}_{h} = \tilde{\mathbf{u}}_{h}^{n+1}$, we gain

Taking $\mathbf{v}_{h} = \dfrac{3\mathbf{d}_{h}^{n+1}-4\mathbf{d}_{h}^{n}+\mathbf{d}_{h}^{n-1}}{2\tau}$, we obtain

Setting $e_{h} = \dfrac{1}{\epsilon^{2}}q_{h}^{n+1}$, we easily get

According to (2.5), we easily obtain $\left(\tilde{\mathbf{u}}_{h}^{n+1}, \mathbf{v}\right) = \left(\mathbf{u}_{h}^{n+1}, \mathbf{v}\right)$. Then applying the fact $\left(a-b, a\right) = \dfrac{1}{2}\left(\|a\|^{2}-\|b\|^{2}+\|a-b\|^{2}\right)$, we have

To approach the pressure term, we rewrite (2.5) as

By squaring both sides of the above equation, we can derive

Then, we have

Noticing the fact

we deduce

Similarly, we get

Combining (3.5)–(3.7) and (3.11), then multiplying both sides by $2\tau$, we conclude

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 $\tau$

Proof. Summing the discrete energy inequality (3.1) from $n = 1$ to $N-1$, we obtain

The proof is complete.

The above second order scheme consists of three temporal levels $n+1, n, n-1$. Thus, $n\geq1$. Then, we have to start with the initial values.

Theorem 3.2. For all $\tau > 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

Choosing $\overline{\mathbf{d}}_{h} = \dfrac{\mathbf{d}_{h}^{1}-\mathbf{d}_{h}^{0}}{\tau}$, we have

Setting $\overline{\mathbf{d}}_{h} = \mathbf{u}_{h}^{1}\cdot\nabla\mathbf{d}_{h}^{0}$, we arrive at

By replacing $\overline{\mathbf{u}}_{h} = \mathbf{u}_{h}^{1}$, we get

Taking $\mathbf{v}_{h} = \dfrac{\mathbf{d}_{h}^{1}-\mathbf{d}_{h}^{0}}{\tau}$, we obtain

Setting $e_{h} = \dfrac{1}{\epsilon^{2}}q_{h}^{1}$, we get

Noticing the fact

we deduce

Similarly, we easily have

Combining (3.22)-(3.24), then multiplying $\tau$ at two sides, we obtain

Finally, we can deduce the discrete energy law. The proof is complete.

4.   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].

4.1.   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 $\xi_{\mathbf{u}}^{h} = \mathbf{u}^{h}(x, T)-\mathbf{u}^{\frac{h}{2}}(x, T)$, $\xi_{\mathbf{d}}^{h} = \mathbf{d}^{h}(x, T)-\mathbf{d}^{\frac{h}{2}}(x, T)$, and obtain the convergence rate $\log_{2}(\xi_{\mathbf{u}}^{h}/\xi_{\mathbf{u}}^{\frac{h}{2}})$, $\log_{2}(\xi_{\mathbf{d}}^{h}/\xi_{\mathbf{d}}^{\frac{h}{2}})$. The problem of Ericksen-Leslie equation is in the domain $[-1, 1]\times[-1, 1]$, and fix the parameter $T = 1$.

We compute the reference solution. The initial conditions are

Table 1 shows the spacial convergence rates of $\|\xi_{\mathbf{u}}^{h}\|_{H^{1}}$ and $\|\xi_{\mathbf{d}}^{h}\|_{H^{1}}$. Parameters are set as $\gamma = 0.35$, $\eta = 0.7$, $\varepsilon = 0.2$, $h = 0.25, 0.125, 0.0625, 0.03125, 0.015625$. The spacial convergence orders computed from the errors $\|\xi_{\mathbf{u}}^{h}\|_{H^{1}}$ and $\|\xi_{\mathbf{d}}^{h}\|_{H^{1}}$ are close to 2. We get the same conclusion in comparison to variable $\epsilon$ in Tables 2–3, where $\epsilon = 0.23, 0.25$.

Table 4 shows the temporal convergence rates of $\|\xi_{\mathbf{u}}^{h}\|_{H^{1}}$ and $\|\xi_{\mathbf{d}}^{h}\|_{H^{1}}$. Parameters are set as $\gamma = 0.11$, $\eta = 0.82$, $\epsilon = 0.197$, $\tau = 0.25, 0.125, 0.0625, 0.03125$. The temporal convergence orders computed from the errors $\|\xi_{\mathbf{u}}^{h}\|_{H^{1}}$ and $\|\xi_{\mathbf{d}}^{h}\|_{H^{1}}$ are close to 2. We get the same conclusion in comparison to variable $\epsilon$ in Tables 5 and 6, where $\epsilon = 0.198, 0.2$. These numerical results show that our numerical algorithm is correct.

4.2.   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]\times[-1, 1]$, and sets the parameters $T = 0.5$, $\gamma = \eta = 1$. The initial conditions are

Figure 1 describes the evolution of elasticity and kinetic energy under different values of variable $\varepsilon$, where $\varepsilon = 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 $\varepsilon$, The velocity is dissipated, and the elastic energy shows a linear and slow decay.

We consider problem with $\gamma = 1, 0.5$ and $n = 16, 3 2$, which are shown in Figure 2. The elastic energy and kinetic energy are sensitive to $\gamma$. With regard to $\eta$ and $\tau$, we also plot the results for different values of $\varepsilon_{d}(t)$ and $\varepsilon_{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 $\eta$ and $\gamma$ has a great influence on the results.

4.3.   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 $\mathbf{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 $\mathbf{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 set the domain as $[-1, 1]\times[-1, 1]$, and the parameters are set as $\tau = 0.0625$, $h = \dfrac{1}{10}$, $\gamma = \eta = 1$. The initial conditions are

The evolution of singularities is depicted in Figure 4. We present snapshots of the director field $\mathbf{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 the initial directors, we give the annihilation 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 $\gamma = \eta = 1$. The initial conditions are

These 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.

5.   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.

Acknowledgments

This work was supported by the Research Project Supported by Shanxi Scholarship Council of China (No.2021-029) and International Cooperation Base and platform project of Shanxi Province (No.202104041101019).

Conflict of interest

All authors declare no conflicts of interest in this paper.