Computational fluid dynamics for nematic liquid crystals

Due to recent advances in fast iterative solvers in the field of computational fluid dynamics, more complex problems which were previously beyond the scope of standard techniques can be tackled. In this paper, we describe one such situation, namely, modelling the interaction of flow and molecular orientation in a complex fluid such as a liquid crystal. Specifically, we consider a nematic liquid crystal in a spatially inhomogeneous flow situation where the orientational order is described by a second rank alignment tensor. The evolution is determined by two coupled equations: a generalised Navier–Stokes equation for flow in which the divergence of the stress tensor also depends on the alignment tensor and its time derivative, and a convection-diffusion type equation with non-linear terms that stem from a Landau-Ginzburg-DeGennes potential for the alignment. In this paper, we use a specific model with three viscosity coefficients that allows the contribution of the orientation to the viscous stress to be cast in the form of an orientation-dependent force. This effectively decouples the flow and orientation, with each appearing only on the right-hand side of the other equation. In this way, difficulties associated with solving the fully coupled problem are circumvented and a stand-alone fast solver, such as the state-of-the-art preconditioned iterative solver implemented here, can be used for the flow equation. A time-discretised strategy for solving the flow-orientation problem is illustrated using the example of Stokes flow in a lid-driven cavity.


Introduction
In recent years, significant advances have been made in the development of effective preconditioned iterative solvers for finite element models in incompressible fluid dynamics, such as solution of Stokes and Navier-Stokes equations [5]. The ready availability of these methods in public domain codes such as IFISS [3] and TRILI-NOS [10] has extended the range of possible applications by making it easier for the practitioner to apply these fast solvers to specific situations. In particular, the MAT-LAB package IFISS, as well as being a useful source of benchmark problems, also provides a convenient starting point for developing solvers dedicated to a particular application [4, §3.4]. In this paper, we describe one such situation where an IFISSbased fast solver is used as part of an algorithm designed to compute flow in a liquid crystal cell.
Computing flow in complex fluids such as liquid crystals and polymers is very challenging because of the altered structure of the flow equation: the underlying Navier-Stokes problem contains additional terms representing the interaction between the flow and the orientation of the molecules within the fluid. In liquid crystal applications, the usual form of the stress tensor is very complicated. Our method relies on reformulating the time derivative in the stress tensor in a way which effectively decouples the flow and orientation, with each appearing only on the right-hand side of the other equation. In this way, difficulties associated with solving the fully coupled problem are circumvented and a stand-alone solver can be used for the flow equation.
The flow of a nematic liquid crystal can be described in various ways. While the most common approach uses the Ericksen-Leslie theory for the nematic director, a more general description using the second rank alignment tensor is needed for problems that involve defects. Different constitutive theories for the alignment tensor have been derived [11,[16][17][18]25,26] , and numerical solutions for some special cases have been produced [6,7,28]. The creation of backflow and its influence on the annihilation of defects in two space dimensions has been examined in [27,29]. In this example, the reorientation of the alignment is the driving force. Also the impact of flow on the orientation has received much attention. Possibly the earliest application was given by Leslie: the flow alignment of the director in a simple shear [13]. The behaviour of the alignment tensor under shear in a monodomain was also extensively studied in [8,19], and other studies have considered lid-driven cavity flow [9,14,33].
Even in a homogeneous, simple shear flow, many different types of behaviour can be found, such as flow aligning, tumbling, and chaos. Furthermore, to obtain a complete picture, spatially inhomogeneous situations have to be considered. The evolution is determined by two equations: the flow is governed by a generalised Navier-Stokes equation, in which the divergence of the stress tensor also depends on the alignment tensor and its time derivative, and the evolution of the orientation is governed by a convection-diffusion type equation that contains non-linear terms that stem from a Landau-DeGennes potential [22].
In this paper we consider a specific model with three viscosity coefficients that allows us to write the contribution of the orientation to the viscous stress in the form of an orientation dependent force. As an example application, we consider the standard fluid flow test problem of Stokes flow in a lid-driven cavity. We propose a time-discretised strategy for solving the flow-orientation problem that involves two alternating steps. First, for a given flow field, one time step for the orientation equation is carried out according to the methods described in [20]. Then, the flow field of the Stokes flow is computed for the given orientation field. This is done using state-of-the-art Krylov subspace and multigrid iteration techniques implemented in IFISS [3].

General underlying equations
We consider a nematic liquid crystal whose orientational order is described by the second rank alignment tensor Q. If u denotes a unit vector parallel to the symmetry axis of an effectively uniaxial molecule, Q can be defined as the local average where I is the identity tensor and . . . denotes the symmetric traceless part of a tensor. Equations of motion for incompressible flow and alignment can conveniently be formulated in terms of a frame-independent invariant rate of the alignment tensor [23].
Here we use the co-rotational time derivative is the skew part of the velocity gradient, with v satisfying the incompressibility constraint If the free energy connected with the alignment is given as a function W = W (Q, ∇Q), the dissipation is specified as a function R = R( • Q, Q, D) that is a quadratic form in • Q, and the symmetric part of the velocity gradient is D = 1 2 (∇v + (∇v) T ), then the equations for flow and alignment take the general form [22,25] where the stress tensor T is given by This tensor contains an isotropic contribution from the hydrostatic pressure p, a viscous stress with symmetric part ∂ R/∂D and skew part Q ∂ R/∂

Q Q, and an elastic stress
which is analogous to the Ericksen elastic stress in a director based description.

Specific model
To obtain a specific model, we choose the free energy to be of the form is a Landau-deGennes potential, and a curvature elastic energy with one elastic constant L 1 is used. Although other models involving additional elastic constants exist (see, for example, [2]), we choose this commonly-used one-constant approximation for simplicity in the equations below. For an alignment tensor theory to be consistent with Ericksen-Leslie theory (in the case of uniaxial alignment with constant scalar order parameter), the dissipation function R needs to contain at least five terms. The choice with five phenomenological viscosity coefficients ζ 1 , ζ 2 , ζ 3 , ζ 31 , and ζ 32 leads to the stress tensor proposed in [18]. Although more general forms for R are available (see [24], [25, eq. (4.23)]), omitting terms other than those in (3.2) simply amounts to neglecting higher-order corrections to the Ericksen-Leslie viscosity coefficients, so we retain the simpler form here. Using (3.1) and (3.2) in (2.4) yields the equation for the alignment where is the derivative ∂φ/∂Q of the Landau-deGennes potential φ. The different contributions to the stress tensor (2.5) then take the following explicit forms: the skew-symmetric part is and the symmetric traceless part of the viscous stress is In the one elastic constant approximation used here in (3.1), the elastic contribution to the stress is symmetric and given by

Solution strategy
We begin by writing the stress tensor in a more convenient form, namely, we remove its explicit dependence on the time derivative of the alignment tensor. To this end, we observe that on a solution Using this in expression (3.3) for the skew part of the viscous stress, we find that where the last equality holds because is simply a polynomial in Q and hence commutes with Q. Applying the same procedure to the symmetric part of the viscous stress yields where we have introduced a renormalised isotropic viscosity ζ 4 according to ζ 4 := ζ 3 − ζ 2 2 /ζ 1 . From now on, we will neglect the last two terms in (4.3), that is, we will assume that ζ 31 = ζ 32 = 0. In terms of the Leslie viscosities, this amounts to making the assumptions α 1 = 0 and α 5 = −α 6 , see [22]. We note that while these assumptions are reasonable for small molecule liquid crystals, they will have to be modified for polymeric liquid crystals (see Sect. 7). The advantage of making these assumptions is that with ζ 31 = ζ 32 = 0, the divergence of the stress tensor takes a very convenient form, namely, We non-dimensionalise by expressing all lengths in terms of the nematic coherence length ξ = 9C L 1 /(2B 2 ) and all times in terms of the relaxation time τ 1 = 9Cζ 1 /(2B 2 ). In addition, we rescale the alignment tensor according toQ = 3C/(2B) Q. This leads to the dimensionless Landau-deGennes potential is a dimensionless temperature parameter. In these units, the clearing point T c and the pseudo critical temperature T * correspond to ϑ = 1 and ϑ = 0, respectively [12]. Note that, for convenience, the tildes are dropped in all subsequent formulae. The final dimensionless equations then are for the orientation and for the flow, together with the incompressibility constraint (2.3). Note that here we have introduced two dimensionless parameters: the backflow parameter, Bf = 4Bζ 2 /(3Cζ 4 ), measures the impact of the orientation on the flow, and the tumbling parameter, Tu = 3Cζ 2 /(2Bζ 1 ), measures the relative strength of the viscosities ζ 2 and ζ 1 . In a simple shear one can expect flow alignment for Tu > 1, where the liquid crystal aligns at an angle of cos 2φ a = −1/Tu to the direction of the flow gradient [13]. For values of Tu < 1 some dynamic state such as tumbling should prevail.

Implementation and solver details
In what follows, we will assume flow at a low Reynolds number, that is, we assume that flow inertia can be neglected sov = 0. Equation (4.8) is then simply a Stokes equation with a force equal to the divergence of the tensor F in (4.9) that depends only on Q and its spatial derivatives. This suggests the following iterative solution strategy: Coupled flow-orientation algorithm Note that, within this framework, any two stand-alone solvers (one for the orientation equations and one for the Stokes equation) can be used. That is, the algorithm structure is independent of the discretisation methods used within each solver, and specific details of the underlying flow and orientation problems (such as shape of the domain and boundary conditions). Furthermore, changing the form of the elastic energy in (3.1) would change only the right-hand side of the flow problem (through F) and would not affect the Stokes iteration matrix (see below).
For the orientation solver, we note that a symmetric, traceless second-rank tensor has five independent degrees of freedom. Once Q is expressed in terms of a suitable basis [20,21], the orientation equation (4.7) takes the form of five coupled non-linear partial differential equations. In the numerical experiments which follow, these were discretised in space using finite differences on a uniform grid and in time using an explicit Euler method. Although stability considerations mean that the size of timestep which can be used is limited with an explicit method, that is not a concern here as small time-steps are already needed for accuracy in terms of modelling flow evolution. Also, the complexity and computational expense of implementing a matrix-based nonlinear iterative solver for the system of five coupled equations in (4.7) makes a fully implicit approach impractical.
As highlighted in the introduction, the Stokes equation (4.8) was solved at each time-step using a finite element based iterative solver adapted from the public domain MATLAB package IFISS [3,4]. The particular finite element discretisation used was a Q 2 − Q 1 Taylor-Hood mixed approximation (that is, quadratic elements for velocity and linear elements for pressure). The resulting linear equations take the form of a saddle point system  [15]. This is clearly not practical for realistic problems, because it involves explicitly inverting A and the Schur complement several times. However, it suggests that choosing P and S to be approximations to A and the Schur complement which are cheap to invert will result in an effective preconditioner M.
For P, any good preconditioner for the Laplacian A can be used. For the Schur complement approximation, we use S = M P , where M P is the finite element mass matrix corresponding to the pressure, which is spectrally equivalent to the Schur complement. Note that, although we invert M P explicitly in the experiments below, if necessary the action of M −1 P can be effectively approximated using a small number of steps of Chebyshev iteration (see [5,Remark 4.5], [31]).
In the numerical experiments of §6, we show results obtained using three different preconditioners of the form (5.2): -Diagonal preconditioning (DP): P = diag(A), S = diag(M P ). This basic preconditioner should offer a modest reduction in iteration counts but is clearly very cheap to invert. -Ideal preconditioning (IP): P = A, S = M P . This represents the best possible preconditioner of this form, as A is inverted exactly. It can be shown that the eigenvalues of M −1 A lie in small intervals that are uniformly bounded away from ±∞ and the origin, meaning that MINRES will converge rapidly and in a number of iterations which is independent of the discrete problem size [30]. Having access to good flow solvers is a key ingredient of our approach, as efficient solution of system (5.1) is critical to the overall practicality of the coupled floworientation algorithm.

Numerical experiments
To illustrate how the coupled flow-orientation algorithm in §5 can be applied in practice, in this section we present the results of some numerical experiments on a lid-driven cavity problem [1]. The lid-driven cavity is a classic test problem in fluid dynamics where flow in a square cavity is driven by the lid moving from left to right, see Fig. 1.
The flow boundary conditions are of Dirichlet type everywhere, with the velocity fixed at some positive rate in the x-direction on the lid and zero along all other cavity walls.
Here we use a 'watertight' cavity, that is, the velocity is fixed to be zero at the top corner points on both left and right boundaries. The resulting discontinuous horizontal velocity generates a strong singularity in the pressure solution, but away from these corners the pressure is essentially constant. Dirichlet boundary conditions are also used for the alignment tensor. The same uniaxial alignment with equilibrium order parameter is prescribed at all boundaries and also as an initial condition in the bulk. In this way, without driving flow a homogeneous uniaxial orientation, as given by the initial condition (see the left of Fig. 2), would result. As mentioned above, alternative problem domains and boundary conditions could be implemented directly in the flow and orientation solvers.

In-plane orientation
For a pure in-plane evolution, we used lid velocity v = 10 and cavity length L = 8. This corresponds to a Reynolds number of Re = V Lρ/ζ 4 ≈ 10 −5 for a typical small-molecule liquid crystal. The Ericksen number is then Er = ζ 1 V L/L 1 ≈ 80, and we chose Bf = 2/3 and Tu = 1/5. The temperature was chosen as ϑ = 0, corresponding to the pseudo-critical temperature T * . The time-step used in the explicit Euler method (for the orientation equations) was t = 0.0001. This ensures stability of the method for the range of spatial discretisation parameters used (from h = 1/16 to h = 1/256 for the experiments reported on below). For orientation boundary conditions, we used homeotropic anchoring on the top and bottom of the cavity and planar anchoring on the lateral sides. The initial orientation is shown on the left of Fig. 2. The boxes shown lie parallel to the eigensystem of the alignment tensor, and the lengths of the edges correspond to the respective eigenvalues, see [20]. The shading of the box shows its degree of biaxiality: a white box corresponds to uniaxial alignment, where two eigenvalues are equal (such as in the initial configuration), and a black box corresponds to perfectly biaxial alignment, where one eigenvalue is zero. In general, the level of darkness of a particular box is proportional to the biaxiality measure β 2 = 1 − 6(tr Q 3 ) 2 /(tr Q 2 ) 3 used in [12]. Note that the number of boxes plotted has been chosen for clear representation of the solutions, and does not correspond to the number of degrees of freedom used in the calculations.  Fig. 3 Flow field during the evolution. The left picture shows the streamlines of the flow field for the initial homogeneous configuration: they are symmetric with respect to a vertical axis through the centre of the cavity. The picture in the middle shows the streamlines at a later time, which are no longer symmetric but, due to the changes in the orientation field, are shifted slightly to the right. The right picture shows a contour plot of the difference between the two flow fields The right picture in Fig. 2 shows a snapshot of the alignment field after the flow has developed. The evolution displayed shows two distinct types of orientation. On the one hand, in the lower part of the cavity, the orientation is dominated by the elastic forces and a stationary state of aligned flow is found. On the other hand, close to the lid where the velocity gradient is large, a periodic solution of in-plane tumbling orientation is found. Furthermore, because of the fixed boundary conditions, the orientation necessarily shows defects. In the alignment tensor description, these defects are characterised by a planar uniaxial orientation. They are generated close to the upper right corner of the lid and travel towards the centre of the cavity and from there to the upper left corner.
For the given choice of the parameters, the flow field is only slightly affected by the orientation (see Fig. 3). Initially, with a homogeneous orientation, the stream lines are symmetric about a vertical axis through the centre of the cavity. This reflects the time-reversal symmetry of the Stokes equation. When the orientation is no longer homogeneous, however, this symmetry is broken and the streamlines shift to the right. This is an effect similar to that found in isotropic fluids at high Reynolds numbers. It is found here in a linear flow equation because of the influence of the orientation on the flow.
To illustrate the efficiency of the flow solver, in Table 1 we present a summary of the performance of the three preconditioners discussed in §5 (as compared to results with unpreconditioned MINRES, which are in the column labelled 'none').
For each method, two quantities are tabulated: k is the average number of MIN-RES iterations required at each time-step to compute the flow field (with convergence tolerance 0.0001), and s is the amount of time associated with this computation (in seconds, as calculated using the MATLAB commands tic and toc). In both cases, the results have been averaged over the first 200 time-steps, as this initial phase poses the greatest challenge for the flow solver. As expected, the number of MINRES iterations required with no preconditioning grows with the problem size. Although using diagonal preconditioning (DP) reduces the iteration count slightly, it can be seen by comparing the values of s that the expense involved in constructing the preconditioner outweighs its benefits as h decreases. For ideal preconditioning, we see that k is, as predicted by the theory in [30], essentially independent of the discretisation parameter h, although the cost of explicitly inverting A grows rapidly. When this inversion is avoided by replacing the action of A −1 by using one multigrid V-cycle based on A, as in MGP, there is a slight growth in the number of iterations needed but, crucially, the method is still essentially grid-independent, and at a much reduced cost. In this framework, MGP is clearly extremely efficient, as is necessary for the overall coupled flow-orientation algorithm to be practical.

Out-of-plane orientation
To obtain an out-of-plane evolution, both the boundary and initial conditions were tilted by an angle of 15 • out of the shear plane. The initial orientation is again homogeneous; with respect to the in-plane orientation on the left of Fig. 2, the top of the alignment tensor is simply tilted by 15 • out of the plane towards the observer. Here we used v = 15 and L = 16, which corresponds to a Reynolds number of 3 × 10 −5 and an Ericksen number of 240. As before, Bf = 2/3 and ϑ = 0, but this time we chose Tu = 4/5 to facilitate the occurrence of out-of-plane periodic solutions (see, for example, the phase diagrams for monodomains in [8]). The time-step used for the orientation solver was t = 0.0001 as before.
The resulting evolution (as illustrated by the snapshot on the left of Fig. 4) again shows in the lower part of the cavity a region of flow alignment that here is out of the plane. Close to the lid, periodic kayaking is found. A close-up view as seen from the top right of the cavity is shown on the right of Fig. 4, where this periodic out-of-plane behaviour is clearly visible. This is again accompanied by the creation of defects in the upper right corner and their annihilation in the upper left corner. A notable difference from the in-plane evolution, however, is that the reduction of the scalar order parameter is far less pronounced. Here, it takes place mostly around the defects: the orientation can go out of the plane to avoid the frustration induced by the flow gradient.
When the orientation has components that lie out of the plane, the force density div F that is due to the orientation-related contributions to the stress can have a component that lies out of the plane even when both flow field and orientation field are assumed to be homogeneous in that direction. This out-of-plane component of div F is shown in Fig. 5. It is most noticeable at the corners where the pressure is divergent, but it is present throughout the cavity. This shows that, for out-of-plane evolutions, truly three-dimensional flow fields will arise and that two-dimensional computations are therefore only of limited value in this context.
The relative performance of the various preconditioned Stokes solvers discussed in §5 is very similar to that for the in plane orientation example of §6.1, so iteration counts and timings have not been displayed here. The MGP preconditioner again significantly outperformed the other methods.

Conclusions
In this paper we have described a highly efficient algorithm for the computation of flow and orientation in nematic liquid crystals. Writing the influence of the orientation field on the flow in the form of a force density allows us to solve the flow equation by using well established fast solvers. This aspect of the modelling dominates the computational time required, so that the overhead added to the computational fluid dynamics by the anisotropic liquid is rather small. We note that, although here we have focussed on the use of a fast solver from the IFISS package, other existing software could also be used.
One disadvantage of the method that we have presented is that only three viscosity coefficients enter the viscous stress. However, when the co-rotational time derivative that we have used here is replaced by a general co-deformational time derivative Q =Q − 2 WQ − 2σ DQ , the same numerical procedure as before can be employed. As long as only the terms proportional to ζ 1 , ζ 2 , and ζ 3 are considered in the dissipation function, the influence of the orientation on the flow still takes the form of a force density This makes the type of algorithm presented in this paper suitable for a more general class of materials, such as polymeric liquid crystals.
Generalisation to high Reynolds numbers is also straightforward: it requires the retention of the inertial term ρv on the left hand side of (4.8) and solution of the resulting Navier-Stokes equation with a specified force term. At each time-step, the latter leads to a saddle-point system of a form similar to (5.1) but with the diffusion component replaced by a discrete convection-diffusion operator. Such systems could be solved efficiently using advanced preconditioned iterative techniques for finite element Navier-Stokes approximations, such as the pressure convection-diffusion and least-squares commutator preconditioners described in [5] and implemented in IFISS.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.