Numerical simulation of the 2D lid-driven cavity flow of chiral liquid crystals

ABSTRACT In this study, the two-dimensional (2D), lid-driven cavity flow of chiral liquid crystals (LCs) was modelled using the Landau de-Gennes (LdG) theory. Parametric studies investigating the effect of Er (Ericksen number) and Θ (chiral strength) on the microstructure of chiral LCs were performed. In this study, we observed that an increase in Θ caused the chiral texture to have more striations and shorter pitch lengths as causes for an increased number of defects. When Θ was held constant, an increase in Er disrupted the chiral structure and even broke it at very high Er. Interestingly, a transition from low Er (10) to moderate Er (1,000) increased the number of defects; however, further increases in Er reduced the number of defects since much of the chiral structure was destroyed by the high viscous flow effects. In particular, even at very high Er, the chiral structure and defects were still present as a vortex was always present with lower velocities, where the viscous flow effect was smaller. We also found that a hexagonal structure with penta-hepta defects formed at high chiral strengths. GRAPHICAL ABSTRACT


Introduction
Liquid crystals (LCs) are a state of matter that is intermediate between crystalline solid and amorphous liquid [1].For years after their discovery, LC materials were thought to have no practical use, but after the 1960s, LCs were found to have broad applications in engineering and sciences, which generated an interest to investigate their behaviour and modelling [1].Traditionally, LCs are known for their successful usage in flat panel displays [2].Nevertheless, in recent years, the research focus has shifted to optics, biotechnology, nanotechnology, and novel composites [2].For example, LCs have been used in energy conversion devices and organic electronics [2].In addition, due to their ability to form ordered layers in the vicinity of boundaries, LCs are used as industrial lubricants to reduce wear and the coefficient of friction, to improve energy efficiencies [3].
Various theories have been used to model the flow of LCs, such as the continuum-based vectorial Ericksen-Leslie (EL) theory and the mesoscopic tensorial Landaude Gennes (LdG) theory.In the LdG theory, a secondorder, symmetric, and traceless tensor order parameter Q is used to describe the microstructure of the LCs [4].The LdG theory is a complete theory since it considers both short-range order elasticity due to intermolecular forces and long-range order elasticity due to surface effects.Therefore, LdG equations coupled with Navier-Stokes equations with a modified stress tensor to account for the anisotropic and viscoelastic contributions can model complex patterns and phenomena in CONTACT Dana Grecov degrecov@mech.ubc.caSupplemental data for this article can be accessed online at https://doi.org/10.1080/02678292.2023.2179121.
LCs, including defect nucleation, which is not possible using the continuum-based EL theory [5][6][7].Other tensor-based models have also been used to describe the flow of LCs [8,9].In the Beris-Edwards model [8], for example, the Poisson bracket formulation is used to derive the governing equations.The Qian-Sheng model [9] can also be used to describe the hydrodynamics of LC flows.
Cholesteric LCs, also known as chiral nematic LCs, are a type of LC that display a helical structure with a twist axis perpendicular to the local molecular director [4,10].The cholesteric LCs possess two structural parameters: the helical pitch (p o ) and the twitch sense.The pitch is defined as the distance along the axis corresponding to a 360° rotation in the orientation of the molecules; the twist sense controls whether the helix is left-or right-handed [10].Much like other types of LCs, cholesteric LCs can have defects.A defect occurs when the local symmetry of the director field is broken [11].Cholesteric LCs have a wide range of applications in engineering and sciences [12], such as thermochromic systems and flexible cholesteric displays.One example of a cholesteric LC is the aqueous suspension of cellulose nanocrystals (CNCs) at certain concentrations [13].Cholesteric LCs exist in most organisms since the helical structure of cholesterics is a fundamental design in convergent evolution [14].For example, DNA, chitin, silk, and cholesterol are examples of cholesteric structures [14].Numerous studies have been conducted on the formation of the chiral structure with or without flow [13,[15][16][17][18][19][20][21][22].Noroozi et al. [13] estimated the viscosity coefficients and rheological functions of CNC aqueous suspensions by modelling the flow of cholesteric LCs with the LdG theory.In [15], Rey used the Ericksen-Leslie theory to model the shear flow of cholesterics and found that the helix of non-aligning cholesteric LCs would not uncoil, and that the helix of flow-aligning cholesterics would only uncoil when the shear rate reached a threshold value.In the case of no flow, Luca and Rey [16] simulated the formation process of biological fibrous composites, which possessed a planar twist structure similar to the chiral structure, with the LdG model.In [17], Luca and Rey developed a mathematical model based on the LdG theory to study the formation of a composite material self-assembled in a lyotropic chiral nematic liquid-crystalline mesophase.Marenduzzo et al. [18] solved the Beris-Edwards equations of motion (a variant of the LdG equations) to simulate the permeative flows of cholesteric LCs.
Venhaus et al. [19] modelled the shear flow of cholesteric LCs using the LdG theory with a chiral term added to the evolution equation of the tensor order parameter Q.In [20], Wiese et al. investigated the rheology and flow-induced morphological changes of cholesteric LC structures under the Poiseuille flow at various anchoring conditions at the wall by solving the Beris-Edwards model.In [21], Pospisil et al. investigated the orientation relaxation dynamics of CNC dispersions and found that parameters such as the chiral strength and the gap confinement played an important role in the formation of the chiral configuration.In [22], Pospisil et al. studied the evaporation front and the development of chiral organisations in the course of casting sulphuric acidhydrolysed CNC films by modelling the system with an LdG model coupled with a mass-transfer equation.In [23] and [24], a hexagonal structure was found when the concentration of DNA molecules was high.In [25], using the LdG theory, Nikzad and Grecov studied the flow of chiral LCs between concentric and eccentric cylinders under various flow conditions and found that the hexagonal structure formed by chiral LCs between cylinders were significantly correlated to chiral strength.
As one of the most researched problems in fluid mechanics, lid-driven cavity flow is regarded as a classical benchmark for computational fluid dynamics (CFD) [26].One of its industrial applications is coating and mixing devices [26].For example, a short-dwell coater for making high-grade paper and the photographic film is described in [27].Lid-driven cavity flow also has broad applications in navigation, meteorology, mining, etc [28].First investigated by Kawaguti [29] and Burggraf [30], lid-driven cavity flow has been studied for both Newtonian and non-Newtonian fluids.In [31], the cavity flow of nematic LCs was simulated using a molecular model based on the LdG theory and it was found that the relaxation of the system to the stable defect configuration present without flows was facilitated by moderate flows.Haque et al. [26] studied cavity flow that extended infinitely in the spanwise direction of non-Newtonian shear-thinning and shear-thickening fluids with viscosity modelled by the Carreau model.They found that the critical Reynolds number initially decreases as the power-index n decreases, but it later increases for highly pseudoplastic fluids.Syrakos et al. [32] modelled the creeping lid-driven cavity flow of a Bingham plastic to evaluate the performance of the finite volume method for solving viscoplastic flows.In [33], Hoang-Trong et al. investigated the cavity flow of a kaolinite suspension, modelled using the Bingham-Papanastasiou approach with various aspect ratios of the cavity.In [34], Yang et al. simulated the 2D liddriven cavity of LCs and identified defects in the flow domain.Nevertheless, to the best of our knowledge, no studies have investigated the 2D lid-driven cavity flow of chiral LCs.Therefore, the objective of this research is to model, simulate, and evaluate the 2D lid-driven cavity flow of cholesteric LCs using the LdG model, and to perform parametric studies on the flow.

Theory and governing equations
In the LdG theory, the microstructure and orientation of LCs are described by a second-order tensor Q, known as the tensor order parameter.All Q tensors at each spatial point form a Q tensor field.The Q tensor can be described as follows [4,13]: where In Equations ( 1) and ( 2), I is the unit tensor; μ n , μ m , and μ l are the eigenvalues of Q, with μ n being the largest, μ m the second largest, and μ l the smallest; n, m, and l are the corresponding eigenvalues.n is the uniaxial director representing the preferred or average orientation of LC molecules.S is the scalar order parameter, whose value is typically between zero and one, with zero describing the absolutely random state (isotropic) and one representing perfect orientational order.P is the biaxial order parameter.Equations ( 1) and (2) show that the scalar order parameter S can be computed from the tensor order parameter Q by finding its eigenvalues and eigenvectors.
As defined originally by de Gennes and Prost [4], an alternative definition of Q is as follows: where u is a unit vector representing the orientation of LC molecules and ω is an orientation distribution function.
In general, the full three-dimensional (3D) Q has the following components: where Q 11 , Q 12 ; Q 13 ,Q 22 ; and Q 23 are independent components of Q tensor and Q is symmetric.
The chiral term is [21,38,39]: where � is the Levi-Civita symbol.The chiral term P � contributes to the formation of director twisting patterns, which are competing factors of the elasticity terms [13,19,39].
In Equation ( 6), the following non-dimensional numbers are chosen: The Deborah number is the ratio of the molecular timescale to the flow timescale: De ¼ U 6HD r ; the Ericksen number is the ratio of the viscous flow effect to longrange elasticity: Er ¼ UHcKT 6L 1 D r ; chiral strength: Θ ¼ L 2 L 2 ; the energy ratio is the ratio of short-range elasticity to longrange elasticity: R ¼ Er De ; Reynolds number is Re ¼ ρU 2 cKT , where U is a characteristic velocity scale; H is a characteristic length scale; D r is the microstructural rotational diffusivity; c is the molecular concentration; K is the Boltzmann constant; T is the absolute temperature; L 1 is the Landau coefficient; L is the isotropic longrange length scale; L is the anisotropic long-range length scale; and ρ is the characteristic density.Note: chiral strength Θ is inversely related to helical pitch p o [21].
Equations ( 8) and ( 9) below are the dimensionless incompressible continuity equation and the momentum equation of the Navier-Stokes equations, respectively.
where u � is the velocity vector and not the unit vector of the LC molecules; p � is the dimensionless pressure.
The dimensionless total stress tensor σ � is defined as the sum of various stress contributions: the symmetric component (τ s � ), anti-symmetric component (τ a � ) and the Ericksen stress term (τ Er � ).The symmetric term � has viscous and elastic contributions.The dimensionless total stress tensor has the following form [36]: In Equation ( 10), the total stress tensor is nondimensionalised, using: where ν 1 � , ν 2 � and ν 4 � are the dimensionless Landau viscosity coefficients.

Methodology
COMSOL Multiphysics ® 5.6 was used to solve the dimensionless LdG equations coupled with the dimensionless Navier-Stokes equations.The full LdG equations, as a set of coupled PDEs, were implemented in the general form PDE (g) module and the momentum equations were solved in the laminar flow (spf) module in COMSOL.Multifrontal massive parallel sparse direct solver (MUMPS) was used to solve for the solution at each time-step and the backward differentiation formula (BDF) was used for time-stepping in the timedependent solver.The governing equations and the simulation set-up for this study had been validated in a previous study [13] in our research group.
Figure 1 shows the computation domain of the 2D lid-driven cavity flow of chiral LCs, which is a square domain with sides H = 1 unit.The lid was driven by a known velocity and all four sides were stationary.Mesh independence was established by comparing velocity and the final mesh had 5,625 elements with an average mesh quality of 1.0.
The tolerance was set to 10 −6 .The initial condition of Q was set to zero for absolute randomness and the initial flow velocity was zero.As suggested in [34], to suppress corner defects in the flow domain, uniform orientational anchoring conditions were prescribed so that LC molecules were perpendicular to the top and bottom sides but parallel to the other sides.No-slip conditions were enforced on all four sides of the flow domain.
To show the evolution of the chiral structure as time progressed, simulations with the following parameters were performed: Θ = 60 with Er = 100.Simulations for Θ = 100 with Er = 100 were conducted to discuss the hexagonal structures.To  investigate the effect of the Ericksen number Er and chiral strength Θ on the flow, the following two sets of parametric studies were conducted: 1) Same chiral strength Θ = 60, different Er (Er = 10, 100, 1,000, 5,000, and 10,000).2) Same Er = 10, different values of chiral strength Θ (Θ = 0, 10, 20, 30, 40, 50, and 60).For all simulations, except for Er and Θ, the following parameters were used: β ¼ 1; U ¼ 4; L � ¼ À 0:3, where β is the molecular shape factor and L � is the dimensionless Landau coefficient [13].The Reynolds number, Re, was 1, which was chosen to ensure that the inertial effect was low so that the viscous effect was dominant over the inertial effect.The Re in this article was defined differently from the conventional Re; hence, it is not suitable to compare its value to those of the conventional Re.The energy ratio, R, was 10 5 .The dimensionless Landau viscosities were as follows [13]: ν 1 � ¼ 0:5057, ν 2 � ¼ 0:3138, and ν 4 � ¼ À 0:4383.

Results and discussion
This  section, the magnitude of the out-of-plane component of the director |n z |, twist angle θ, and tilt angle φ [40] are shown in the Supplemental Online Material.The twist angle θ, and the tilt angle φ [40] are defined in Figure 2. All results in this section including the times are dimensionless.

Time evolution of the chiral structures
Figure 3 is a schematic representation of chiral LC molecules and the scalar order parameter S at different times when Θ = 60 and Er = 100.As discussed in [13], because the eigenvalues of Q can be negative, Q is augmented by 1/3 to ensure that the eigenvalues can no longer be negative.Therefore, a new tensor M ¼ Q þ 1 3 I was created to describe the microstructure of LCs.With this technique, the chiral LC molecules are represented by triaxial ellipsoids.To draw the schematic representation of chiral LC molecules, eigenvalues and eigenvectors of M are first found and then the eigenvectors are ranked according to the value of their corresponding eigenvalues.The triaxial ellipsoids were constructed so that the major axis is along the eigenvector corresponding to the largest eigenvalue.The eigenvectors corresponding to the second largest and the smallest eigenvalues determine the directions of the first and the second minor axes of the ellipsoids, respectively [13].
As Figure 3 shows, the location of defect cores corresponds to reduced scalar order parameter S, which is discussed in [41] and [42].As identified in Figure 3(a), all defects present in the flows are of strength s = −1/2.Defects are seen to migrate over time.The presence of defects in lid-driven cavity flows of LCs has been demonstrated by Yang and Forest [34].Over time, although velocity streamlines become more and more stable, the microstructure still changes slightly.This phenomenon was observed in [3] and [34] and is described as the 'quasi-stationary condition'.

The hexagonal structure
Figure 4 shows the chiral structural visualisations and the scalar order parameter S when Θ = 100 and Er = 100 at t = 500.With a higher chiral strength, the striations in texture become shorter, indicating a smaller chiral pitch length, which is consistent with the findings in [21].In addition, many more regions with low S values are seen with the weaker chiral strength (Θ = 60), meaning that more defects occur when Θ is high.The existence of multiple hexagonal structures at Θ = 100 is also noticeable in Figure 4 as a result of the high chiral strength [25,43].As shown in Figure 4, in addition to hexagons, pentagons and heptagons also exist, which are the typical penta-hepta defects of hexagonal structures discussed in [17] and [24].

Er effect
The effect of Er on the chiral structure was also investigated.Since Er is the ratio of the viscous flow effect to long-range elasticity, increasing Er while holding R (the energy ratio) constant corresponds to an increase in the viscous flow effect.Figure 5 shows how an increase in Er affects the chiral structure and the scalar order parameter when the chiral strength Θ is constant.Figure 6 shows the number of defect cores at various Er.shows the streamlines corresponding to Figure 5(e) with colorbars indicating the velocity magnitudes.From the plots, as Er increases from 10 to 100 to 1,000, the number of defects also increases.In this range of low to medium Er, the increased viscous flow effect disrupts the chiral structure, resulting in more defects.Nevertheless, as Er increases further, the number of defects in the chiral structure begins to decrease due to defect annihilation and breakage of the chiral structure.
In particular, at Er = 10,000, most defects beyond the centre of the vortex disappear, as the complete breakage of the chiral structure occurs.The remaining chiral structure near the centre of the vortex is the result of low viscous effects associated with low velocity in that region (Figure 7).Our results agree with previous reports that when the viscous flow effect is strong, the chiral structure weakens or even disappears [13].
Notably, in our results, the chiral structure still exists at high Er, which is not observed in the simple shear flow of cholesteric LCs [13].

Chirality effect
In addition to Er, chiral strength (Θ) is another important parameter in the cavity flow of chiral LCs.In this section, we discuss how Θ affects the chiral structure and the flow when Er = 10, a relatively small value.
A small Er together with a small De, while the energy ratio R is held constant, means that both elasticity and viscous flow effects are small, which makes the effect of Θ on the flow more noticeable.Figure 8 shows the effect of increasing Θ on the chiral structure and the scalar order parameter when Er is held constant.In Figure 8, the chiral structure begins to emerge when Θ is increased from 0. As the chiral strength Θ increases, the striations shorten, indicating a smaller chiral pitch length, which agrees with the findings in [21].Figure 9 shows the number of defect cores when Θ increases from 10 to 60.The general trend is that the number of defects gradually increases as the chiral strength Θ becomes stronger.Figure 10 shows the streamlines corresponding to Figure 8(a-e) with colorbars indicating the velocity magnitudes.Figures 8(a  not only increases the velocity magnitude but also changes its distribution (Figure 10(b)).By comparing Figure 8 with Figure 10, the flow-structure coupling can be observed as the striations of the chiral structure are roughly aligned with the streamlines.The flow-structure coupling is expected as the LdG equations are coupled with the Navier-Stokes equations in the simulations.

Conclusion
In this research, the 2D lid-driven cavity flow of chiral LC was modelled, simulated, and evaluated with a variety of parametric studies to investigate the effect of chiral strength Θ and the viscous flow effect Er on the microstructure of the aforementioned flow.We found that an increase in the chiral strength caused the chiral texture to have more striations and shorter pitch lengths, leading to an increased number of defects.When the chiral strength was held constant, an increase in Er disrupted the chiral structure and even breaks it at very high Er.Interestingly, a transition from low Er (10) to moderate Er (1,000) increased the number of defects, but a further increase in the viscous flow contribution reduced the number of defects since much of the chiral structure was destroyed by the very high viscous flow effects.Notably, even at very high Er, the chiral structure and defects were still present since a vortex was always present with lower velocities, when the viscous flow effect was smaller.At high chiral strength, we also noticed the hexagonal structure with penta-hepta defects formed, and as expected, the flow-structure coupling was observed in all simulations.
To the best of our knowledge, this research is the first study to simulate the 2D lid-driven cavity flow of chiral LCs modelled by coupling LdG equations with Navier-Stokes equations with a modified stress tensor to account for the anisotropic viscoelastic properties of LCs.The results of the parametric studies of the fundamental flow are a helpful benchmark for the flow of chiral LCs.The lid-driven cavity flow also has industrial applications such as with coating and mixing devices.This study is important since the parametric studies allow for a reduction in defects in the flow of chiral LCs, which are generally undesirable in the application of LCs.To further evaluate the results, future experiments on the flow of chiral LCs should be conducted to compare the experimental results with the simulation results.In particular, the microstructure of LC flows can be visualised by polarised optical microscopy (POM) and the velocities can be measured by particle image velocimetry (PIV).In addition, parametric studies investigating the effect of Re on the microstructure of chiral LCs can be performed.
section shows the results of a range of parametric studies conducted on the 2D lid-driven cavity flow of chiral LCs.A range of snapshots of the chiral structure and the scalar order parameter S showing the evolution of the chiral structure at Θ = 60 and Er = 100 are shown first.The simulation results for Θ = 100 and Er = 100 are then presented to discuss the appearance of the hexagonal structures.The simulation results for the two set of parametric studies are also shown to investigate the effects of Er and Θ.In addition to the results presented in this