Numerical approximations for a smectic-A liquid crystal flow model: First-order, linear, decoupled and energy stable schemes

In this paper, we consider numerical approximations for a model of smectic-A liquid crystal flows in its weak flow limit. The model, derived from the variational approach of the de Gennes free energy, is consisted of a highly nonlinear system that couples the incompressible Navier-Stokes equations with two nonlinear order parameter equations. Based on some subtle explicit-implicit treatments for nonlinear terms, we develop an unconditionally energy stable, linear and decoupled time marching numerical scheme for the reduced model in the weak flow limit. We also rigorously prove that the numerical scheme obeys the energy dissipation law at the discrete level. Various numerical simulations are presented to demonstrate the accuracy and the stability of the scheme.

1. Introduction. Liquid crystal (LC) is often viewed as the fourth state of the matter besides the gas, liquid and solid. It may flow like a liquid, but its molecules may be oriented in a crystal-like way. There are many different types of liquidcrystal phases, which can be distinguished by their different optical properties. Thermotropic LCs can be distinguished into two main different phases: Nematic and Smectic. In Nematic phases, the rod-shaped molecules have no positional order, but molecules self-align to have a long-range directional order with their long axes roughly parallel. Thus, the molecules are free to flow and their center of mass positions are randomly distributed as in a liquid, although they still maintain their long-range directional order. In smectic phases, which are found at lower temperatures, the well-defined layers form, that can slide over one another in a manner similar to that of soap. The smectics are thus positionally ordered along one direction inside the layer. There are many different smectic phases, all characterized by different types and degrees of positional and orientational order (cf. [3, 8, 9, 14, 47, 49-51, 60, 62, 63] and the references therein). In particular, in Smectic-A phase, molecules are oriented along the normal vector of the layers, while in Smectic-C phases they are tilted away from the normal vector of the layer.
The mathematical model of liquid crystals can often be derived from an energybased variational formalism (energetic variational approach), leading to well-posed nonlinearly coupled systems that satisfy thermodynamics-consistent energy dissipation laws. This makes it possible to carry out mathematical analysis and design numerical schemes which satisfy a corresponding energy dissipation law at the discrete level. For smectic-A phase liquid crystals, in de Gennes' pioneering work [3], the phenomenonlogical free energy of smectic-A phase is presented by coupling two order parameters which characterize the average direction of molecular alignment, as well as the layer structure, respectively. In [2], the authors modified the de Gennes' model by adding a second order gradient term for the smectic order parameter to investigate the nematic to smectic-A or smectic-C phase transition, and to predict the twist grain boundary phase in chiral smectic liquid crystals. In [22], the authors used the de Gennes energy to study smectic-A liquid crystals to simulate the chevron (zigzag) pattern formed in the presence of an applied magnetic field. In [5], the authors derived the hydrodynamics coupled model for smectic-A phase by assuming that the director field is strictly equal to the gradient of the layer, thus the free energy is reduced to one order parameter.
From the numerical point of view, it is specifically desired to design numerical schemes that could preserve the thermo-dynamically consistent dissipation law (energy-stable) at the discrete level, since the preservation of such laws is critical for numerical methods to capture the correct long time dynamics. The noncompliance of energy dissipation laws may lead to spurious numerical solutions if the grid and time step sizes are not carefully controlled. To the best of the author's knowledge, although a variety of the smectic liquid crystal models had been developed, we notice that the successful attempts in designing efficient energy stable schemes are very scarce due to the complex nonlinearities. For instances, in [16], the authors present a temporal second order numerical scheme to solve the model of [5]. The scheme is energy stable, however, it is nonlinear thus the implementation is complicated and the computational cost might be high. In [22], the authors developed a temporal first order scheme. However, it does not follow the energy dissipation law even though the schemes are linear and decoupled. Therefore, the main purpose of this paper is to construct the efficient schemes to solve the reduced de Gennes type smectic-A liquid crystal model in the weak flow limit (cf. [3,22]). To solve the model, the main difficulties roughly include (i) the coupling between the velocity and director field/layer function through the convection terms and nonlinear stresses; (ii) the coupling of the velocity and pressure through the incompressibility constraint; (iii) the nonlinear coupling between the director field and layer function. In this paper, we develop a time discretization scheme which (a) is unconditionally stable; (b) satisfies a discrete energy law; and (c) leads to linear, decoupled equations to solve at each time step.
The rest of the paper is organized as follows. In Section 2, we present the whole model and present the PDE energy law. In Section 3, we develop the numerical scheme and prove the unconditional stability. In Section 4, for the reduced model, we present some numerical experiments to validate the proposed scheme. Finally, some concluding remarks are presented in Section 5.
2. The smectic-A liquid crystal fluid flow model and its energy law. The de Gennes free energy of smectic A liquid crystal is described by a unit vector (director field) d and a complex order parameter ψ, to represent the average direction of molecular alignment and the layer structure, respectively. The smectic order parameter is written as where ω(x) is the order parameter to describe the layer structure so that ∇ω is perpendicular to the layer. The smectic layer density ρ(x) is the mass density of the layers. The de Gennes free energy reads as follows, where the order parameters C, K, g, r are all fixed positive constants, h is the unit vector representing the direction of magnetic field and H 2 is the strength of the applied field. Ω = (−L, L) 2 × (−d, d). Now we consider the simple case by assuming the density ρ(x) = r/g [22], then the energy becomes Let φ(x) = ω(x) d , thus the normalized energy becomes (2.5) The dimensionless parameter η is in fact the ratio of the layer thickness to the sample thickness and thus η 1. To release the unit vector constraint of |d| = 1, a nonlinear potential G(d) = 1 4 2 (|d| 2 − 1) 2 which is a Ginzburg-Landau type penalty term, is added into the free energy to approximate the unit length constraint of d [23] , where 1 is a penalization parameter. Thus the modified total free energy with the hydrodynamics reads as follows where u is the fluid velocity field. Assuming a generalized Fick's law that the mass flux be proportional to the gradient of the chemical potential [7,24], we can derive the following system: where p is the pressure, σ d is the Caughy stress tensor, ν is the viscosity, is vorticity tensor, a is a geometry parameter of the liquid crystal molecule, M 1 , M 2 are the relaxation order parameters. The variational derivatives µ φ and µ d are . Following the work in [5], the Cauchy stress tensor σ d reads as follows, where µ i > 0 and D(u) = 1 2 ((∇u) + (∇u) T ) is the strain tensor. In this paper, we assume µ 1 , µ 5 are negligible comparing to µ 4 , thus the stress tensor is simplifed to (2.14) For simplicty, we assume the boundary conditions as follows.
To obtain the dissipation law of the system (2.7)-(2.9), we take the L 2 inner product of (2.7) with δE δφ , (2.8) with δE δd , and (2.9) with u, perform the integration by parts, and add all equalities together, we obtain 3. Numerical scheme. For the algorithms development, we only consider the simplified or reduced model where the two tensors W (u), D(u) and the elastic stress are neglected since that is extremely hard to develop linear and unconditionally energy stable schemes with these highly nonlinear terms. The interested readers can refer to a nonlinear type schemes that was constructed in [60] with these tensors. The reduced model applies to smectic-A liquid crystal flows in in the weak flow limit.
The emphasis of our algorithm development is placed on designing numerical schemes that are not only easy-to-implement, but also satisfy a discrete energy dissipation law. We will design schemes that in particular can overcome the following difficulties, namely, • the coupling of the velocity and pressure through the incompressibility condition; • the stiffness in the director equation associated with the penalty parameter ; • the nonlinear couplings among the fluid equation, the layer equation and the director equation.
We construct an energy stable scheme based on a stabilization approach [31,32]. To this end, we shall assume that G(d) satisfies the following conditions, i.e., One immediately notes that this condition is not satisfied by this double-well potential G(x). However, it is a common practice that one truncates this fourth order polynomial G to quadratic growth outside of an interval [−M, M ] without affecting the solution if the maximum norm of the initial condition |d 0 | is bounded by M . Therefore, one can (cf. [31]) consider the truncated double-well potentialG(d) by modifying this function outside a ball in {x : |d(x)| ≤ 1} ∈ R 3 of radius 1 as follows.
Hence, there exists a positive constant L G such that For convenience, we consider the problem formulated with the substituteG, but omit the˜in the notation. When deriving the energy law (2.16), we notice that the nonlinear terms in δEtot δφ and δEtot δd involve second order derivatives, and it is not convenient to use them as test functions in numerical approximations, making it difficult to prove the discrete energy dissipation law. To overcome this difficulty, we first reformulate the system (2.7)-(2.10) in an alternative form which is convenient for numerical approximations. The system reads as follows, whereφ = φ t + u · ∇φ andḋ = d t + u · ∇d. To obtain the dissipation law of the system (3.4)-(3.7), we take the L 2 inner product of (3.4) with φ t , (3.5) with d t , and (3.6) with u, perform the integration by parts, and add all equalities together. We obtain We now fix some notations. For scalar function u, v and vector function u = (u 1 , u 2 , u 3 ) and v = (v 1 , v 2 , v 3 ), we denote the L 2 inner product as follows.
In the above, S is a stabilizing parameter to be determined. We have the following remarks in order:

Remark 1.
• A pressure-correction scheme is used to decouple the computation of the pressure from that of the velocity. For the Navier-Stokes equations, there exist a large quantity of literatures concerning the subjects about theoretical/numerical analysis, scheme developments, as well as their applications in various research fields, see [1, 4, 10-13, 15, 19, 21, 25-30, 35-45, 58]. • The nonlinear term g(d) mainly takes the form like 1 2 d(|d| 2 − 1), so the explicit treatment of this term usually leads to a severe restriction on the time step δt when 1. Thus we introduce in (3.10) a linear "stabilizing" term to improve the stability while preserving the simplicity. It allows us to treat the nonlinear term explicitly without suffering from any time step constraint [31,34,48,57,59,60]. Note that this stabilizing term introduces an extra consistent error of order O(δt) in a small region near the interface, but this error is of the same order as the error introduced by treating it explicitly, so the overall truncation error is essentially of the same order with or without the stabilizing term. It is remarkable that the truncation error of the stabilizing approach is exactly the same as the convex splitting method [6].
• Inspired by [32], which deal with a phase-field model of three-phase viscous fluids or complex fluids, we introduce two new, explicit, convective velocities u n and u n in the phase equations. u n and u n can be computed directly from (3.11) and (3.13), i.e., It is easy to get det(I + c(∇φ) T ∇φ) = 1 + c∇φ · ∇φ, thus the above matrix is invertible. • The scheme (3.10)-(3.15) is a totally decoupled, linear scheme. Indeed, (3.10), (3.12) and (3.14) are respectively (decoupled) linear elliptic equations for φ n+1 , d n+1 andũ n+1 , and (3.15) can be reformulated as a Poisson equation for p n+1 − p n . Therefore, at each time step, one only needs to solve a sequence of decoupled elliptic equations that can be solved very efficiently. • The scheme (3.10)-(3.15) is first order accurate in time. To develop a linear, second order and unconditionally energy stable scheme, one can consider to use the newly developed "Invariant Energy Quadratization" approach (cf. the authors recent work in [17,33,46,[52][53][54][55][56]61]). However, IEQ type schemes will couple the computations of d, φ with the velocity u together, even though the whole system is still linear. It is quite an open problem on how to develop a linear, second order scheme that can decouple the computations of all variables efficiently. • As we shall show below, the above scheme is unconditionally energy stable.
For Term A, we apply the Taylor expansions to obtain For Term B, we have (3.34) For Term C, we have Finally, we obtain the desired result after dropping some positive terms.   The boundary conditions are Neumann type along y-axis (cf. (2.15)). We use 128 × 128 grid points to discretize the space, and perform the mesh refinement test for time. We choose the numerical solution with the time step size δt = 1 × 10 −4 as the benchmark solution (approximate exact solution) for computing errors. Figure 4.1 plots the L 2 errors for various time step sizes. We observe that the scheme is asymptotically first-order accurate in time for all variables as expected. In Figure 4.2, we compare the time evolution for three different time steps of δt = 0.0001, 0.001, 0.01. We observe that for smaller time steps of δt = 0.0001 and 0.001, the two energy curves coincide very well. For larger time steps 0.01, the energy curve deviates viewable away from others. This means the adopted time step size should not be larger than 0.01, in order to get reasonably good accuracy. Thus we choose δt = 0.001 for all simulations.

4.2.
Chevron pattern induced by the magnetic force. We now consider the effects from the magnetic force for the no flow case (u = 0). Initially, a smectic A liquid crystal is confined between two flat parallel plates and uniformly aligned in a way that the smectic layers are parallel to the bounding plates and the directors are aligned homeotropically, that is, perpendicular to the smectic layers. A magnetic field is applied in the direction parallel to the smectic layers, which induce the layer  undulation (chevron pattern) phenomena. The initial conditions read as follows.
d(t = 0) = (0, 1) + 0.001(rand(x, y), rand(x, y)), where the rand(x, y) is the small perturbation that is the random number in [−1, 1] and has zero mean. We set the Dirichlet type boundary condition for φ and d as follows, d| y=±1 = (0, 1), φ| y=1 = 1, φ| y=−1 = −1. We take δt = 0.001 to obtain better accuracy. Fig. 4.3 shows the snapshots of the layer function φ at t = 0, 0.2, 0.4 and 0.8. Initially at t = 0, the layer function take the linear profile along the y− axis. When time evolves, we observe some undulations appear at t = 0.2. The layer function quickly reaches the steady solution  at t = 0.8 with the saw tooth shape. This undulation phenomenon is called the Helfrich-Hurault effect (cf. [18,20]). Fig. 4.4 shows the snapshots of the director field d. The numerical solution presents similar features to those obtained in [22]. We also plot the energy dissipative curve in Fig. 4.5, which confirms that our algorithm is energy stable.
5. Concluding remarks. In this paper, we consider the numerical approximations for a hydrodynamically coupled smectic-A model. For the simplified model of smectic-A liquid crystal flows, we have presented an efficient, semi-discrete in time, numerical scheme with provably unconditionally stability. Based on the linear stabilized approach, the proposed scheme conquer the inconvenience from nonlinearities by linearizing all nonlinear terms explicitly. We show that the linear scheme developed is unconditionally energy stable. We verify numerically that our scheme is of first order accurate in time and present ample numerical results from some standard numerical tests in the presence of the magnetic field.