Lattice-Boltzmann simulation of free nematic-isotropic interfaces

TL;DR

This paper employs a hybrid lattice Boltzmann and finite difference approach to simulate and analyze the behavior of nematic-isotropic interfaces in liquid crystals, focusing on interface stability, velocity, and domain stabilization.

Contribution

It introduces a hybrid simulation method to study nematic-isotropic interfaces, providing new insights into interface stability, velocity, and domain stabilization near coexistence.

Findings

Interface is static at coexistence temperature

Profile width matches theoretical predictions

Circular domains stabilized by temperature shifts

Abstract

We use a hybrid method of lattice Boltzmann and finite differences to simulate flat and curved interfaces between the nematic and isotropic phases of a liquid crystal described by the Landau-de Gennes theory. For the flat interface, we measure the interfacial velocity at different temperatures around the coexistence. We show that the interface is completely static at the coexistence temperature and that the profile width is in line with the theoretical predictions. The interface is stable in a range of temperatures around coexistence and disappears when one of the two phases becomes mechanically unstable. We stabilize circular nematic domains by a shift in temperature, related to the Laplace pressure, and estimate the spurious velocities of these lattice Boltzmann simulations.

Equations

$${\partial }_t{Q}_{\alpha \beta }+u_{\gamma }{\partial }_{\gamma }{Q}_{\alpha \beta }-S_{\alpha \beta }(W_{\alpha \beta },Q_{\alpha \beta })=\Gamma H_{\alpha \beta }.$$

$$S_{\alpha \beta }=$$

$$\displaystyle\mathcal{F}=\int_{V}d^{3}r\Big{[}\frac{A_{0}}{2}\left(1-\frac{\gamma}{3}\right)Q_{\alpha\beta}Q_{\alpha\beta}-\frac{A_{0}\gamma}{3}Q_{\alpha\beta}Q_{\beta\gamma}Q_{\gamma\alpha}+\frac{A_{0}\gamma}{4}(Q_{\alpha\beta}Q_{\alpha\beta})^{2}+\frac{K}{2}(\partial_{\gamma}Q_{\alpha\beta})(\partial_{\gamma}Q_{\alpha\beta})\Big{]},$$

$$H_{\alpha \beta }=-A_0\left(1-\frac{\gamma }{3}\right)Q_{\alpha \beta }-A_0\gamma \left(Q_{\nu \gamma }^2\frac{{\delta }_{\alpha \beta }}{3}-Q_{\alpha \gamma }{Q}_{\beta \gamma }\right)-A_0\gamma Q_{\alpha \beta }{Q}_{\gamma \lambda }^2+K{\partial }_{\gamma }{\partial }_{\gamma }{Q}_{\alpha \beta }.$$

$${\partial }_t\rho +{\partial }_{\alpha }(\rho u_{\alpha })=0,\rho {\partial }_t{u}_{\alpha }+\rho u_{\beta }{\partial }_{\beta }{u}_{\alpha }={\partial }_{\beta }{\Pi }_{\alpha \beta }+\eta {\partial }_{\beta }\left({\partial }_{\alpha }{u}_{\beta }+{\partial }_{\beta }{u}_{\alpha }\right).$$

$${\Pi }_{\alpha \beta }=$$

$$-{\partial }_{\alpha }{Q}_{\gamma \nu }\frac{\delta \mathcal{F}}{\delta ({\partial }_{\beta }{Q}_{\gamma \nu })}+Q_{\alpha \gamma }{H}_{\gamma \beta }-H_{\alpha \gamma }{Q}_{\gamma \beta }.$$

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)-f_i(\mathbf{x},t)={\left(\partial f_i/\partial t\right)}_{\text{coll}}+S_i,$$

$$f_i^{eq}=\rho w_i\left[1+\frac{{\mathbf{c}}_i\cdot \mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_i\cdot \mathbf{u}{)}^2}{2c_s^4}-\frac{{\mathbf{u}}^2}{2c_s^2}\right],$$

$$\rho =\sum _i{f}_i,\rho \mathbf{u}=\sum _i{\mathbf{c}}_i{f}_i+\frac{{\mathbf{F}}_i\Delta t}{2}.$$

$$S=\frac{S_n}{2}\left[1-\tanh \left(\frac{x-x_i}{2\lambda }\right)\right],$$

Full text

1 Introduction

Liquid crystals are systems that flow as liquids but where the particles are partially organized in crystal-like ways p1995physics ; beris1994thermodynamics . They are formed by elongated particles the orientations of which may point in a preferential direction. At low temperatures, these liquid crystals are in the nematic phase (orientationally ordered) while at high temperatures, they are in the isotropic phase, where the orientation of the particles is random (orientationally disordered). At the nematic-isotropic temperature, the two phases coexist and a stable interface may be observed. The standard theory that describes this static interface is based on a tensor order parameter expansion, with gradient terms that account for the order parameter variations, proposed by de Gennes p1995physics .

To simulate the dynamics of liquid crystals, it is common to solve numerically the continuum equations of Beris-Edwards and Navier-Stokes marenduzzo07steady ; C9SM00859D ; Doostmohammadi2016 . A hybrid method of lattice Boltzmann kruger2016lattice and finite differences vesely2001computational is used. Special care has to be taken in order to simulate interfaces due to spurious currents that arise in these numerical methods, which can lead to interfacial motion or changes in the director field kruger2016lattice ; PhysRevE.82.046708 .

We analyze the applicability of the method of Ref. C9SM00859D to simulate free interfaces, which are very sensitive to these spurious numerical effects. We show that it is possible to simulate static flat interfaces at coexistence, with the profiles predicted by the Landau-de Gennes theory. We estimate the magnitude of the spurious velocities at nearly static circular interfaces stabilized by temperature shifts, related to the Laplace pressure arising from the curvature.

1.1 Equations of motion

To simulate the liquid crystal dynamics, we solve the Beris-Edwards equation coupled to the Navier-Stokes equation. The first describes the evolution of the order parameter $Q_{\alpha\beta}=S(n_{\alpha}n_{\beta}-\delta_{\alpha\beta}/3)$

[TABLE]

Here, $\xi$ stands for the alignment parameter, $\Gamma$ is the rotational diffusive constant and the co-rotational term is given by:

[TABLE]

where $W_{\alpha\beta}=(\partial_{\beta}u_{\alpha}-\partial_{\alpha}u_{\beta})/2$, $D_{\alpha\beta}=(\partial_{\beta}u_{\alpha}+\partial_{\alpha}u_{\beta})/2$. The Landau-de-Gennes free energy reads

[TABLE]

where $A_{0}$ is a constant, $\gamma$ is a temperature related parameter and $K$ is the elastic constant. Thus, the molecular field, $H_{\alpha\beta}=-\delta\mathcal{F}/\delta Q_{\alpha\beta}+(\delta_{\alpha\beta}/3)\operatorname{Tr}(\delta\mathcal{F}/\delta Q_{\gamma\epsilon})$, is

[TABLE]

The continuity and Navier-Stokes equations describe the evolution of the density and velocity of the fluid

[TABLE]

where, the stress-tensor is

[TABLE]

$P_{0}=\rho c_{s}^{2}$ is the hydrostatic pressure. Using Eq. (3), we obtain $\delta\mathcal{F}/\delta(\partial_{\beta}Q_{\gamma\nu})=K\partial_{\beta}Q_{\gamma\nu}$, which can be replaced in Eq 6.

1.2 Hybrid method

Our numerical scheme solves Eq. (1) using finite differences and Eq. (5) using lattice Boltzmann. Both methods are solved using the same grid (spatial discretization) and the same time step.

1.2.1 Lattice Boltzmann

The lattice Boltzmann method is a numerical technique which solves the Boltzmann equation and recovers the Navier-Stokes equation in the macroscopic limit. The space is discretized in a regular grid and the velocity space is discretized according to the lattice (Gaussian quadrature). Here, we use the D3Q19 lattice, which has 19 velocity vectors $\mathbf{c}_{i}$ isotropically distributed in three dimensions. The discrete version of the Boltzmann equation is implemented in the method

[TABLE]

where $f_{i}$ is the distribution function corresponding to the $\mathbf{c}_{i}$ vector. The simplest collision operator is the BhatnagarâGrossâKrook (BGK) one, which assumes that the full-equilibrium distribution $f_{i}$ relaxes to the equilibrium $f_{i}^{eq}$ with a characteristic relaxation time $\tau$: $\left(\partial f_{i}/\partial t\right)_{\text{coll}}=-(\Delta t/\tau)(f-f^{eq})$. The multi-relaxation time operator, used in this work, is a generalization of the BGK: $\left(\partial f_{i}/\partial t\right)_{\text{coll}}=\mathbf{M}^{-1}\mathbf{R}\mathbf{M}[f_{i}(\mathbf{x},t)-f_{i}^{eq}(\mathbf{x},t)]\Delta t$, where the matrix $\mathbf{M}$ transforms from the distribution function space to the hydrodynamic moments space and $\mathbf{R}$ is the relaxation matrix. It assumes that the hydrodynamic moments may relax with different relaxation times. This approach is known to improve the accuracy and stability in many problems kruger2016lattice . The equilibrium distribution is the second order expansion in Hermite polynomials of the Maxwell-Boltzmann distribution COELHO2018144 :

[TABLE]

where $c_{s}$ is the speed of sound in the given lattice. For the D3Q19, it is $c_{s}=1/\sqrt{3}$. The density and velocity fields are calculated as follows:

[TABLE]

The source term, for the BGK operator, reads: $S_{i}=w_{i}\left(1-\frac{1}{2\tau}\right)\left[\frac{\mathbf{c}_{i}}{c_{s}^{2}}\left(1+\frac{\mathbf{c}_{i}\cdot\mathbf{u}}{c_{s}^{2}}\right)-\frac{\mathbf{u}}{c_{s}^{2}}\right]\cdot\mathbf{F}$. For the MRT operator, this source term must be relaxed in the moments space as the distribution function (see Refs. C9SM00859D ; kruger2016lattice for more details of the MRT method). The force is calculated using the stress tensor of Eq. (6): $F_{\alpha}=\partial_{\beta}(\Pi_{\alpha\beta}+P_{0}\delta_{\alpha\beta})$.

1.2.2 Finite differences

The Eq. (1) is solved explicitly using a predictor-corrector finite difference method. All the differences are second order accurate. For instance, the first derivative of a vector $V$ in the $x$-direction is: $\partial_{x}V(x)=[V(x+\Delta x)-V(x-\Delta x)]/(2\Delta x)+\mathcal{O}(\Delta x^{2})$. To calculate the time evolution of $Q_{\alpha\beta}$, we follow the following steps: 1) Calculate the time derivative, $(\partial Q_{\alpha\beta}/\partial t)^{\prime}$, using Eq. (1) with the current $Q_{\alpha\beta}$; 2) Calculate the predictor: $Q^{P}_{\alpha\beta}=Q_{\alpha\beta}+\Delta t(\partial Q_{\alpha\beta}/\partial t)^{\prime}$; 3) Calculate the time derivative using the predictor, $(\partial Q_{\alpha\beta}^{P}/\partial t)$, with Eq. (1); 4) Calculate the corrected derivative: $\frac{\partial Q_{\alpha\beta}}{\partial t}=\frac{1}{2}\left[\left(\frac{\partial Q_{\alpha\beta}}{\partial t}\right)^{\prime}+\left(\frac{\partial Q_{\alpha\beta}^{P}}{\partial t}\right)\right]$; 5) Calculate the corrector: $Q_{\alpha\beta}(t+\Delta t)=Q_{\alpha\beta}(t)+\Delta t\frac{\partial Q_{\alpha\beta}}{\partial t}$; 6) Return to step 1. Notice that the velocity used in Eq. (1) is the actual one (and not the lattice Boltzmann one: $\rho\mathbf{u}_{LB}=\sum_{i}\mathbf{c}_{i}f_{i}$), given by Eq. (9), which is second order accurate.

2 Free interface

In this section, we show the results from numerical simulations of unconfined nematic-isotropic interfaces. First, we consider a flat interface in two dimensions and then we simulate a circular interface. In the simulations, we apply periodic boundary conditions in both directions. The initial velocity is set to zero and the density to $\rho=1$ everywhere. Other parameters are: $\tau=1.5$, $K=0.04$, $\xi=0.7$, $A_{0}=0.1$ and $\Gamma=0.34$. Our results are given in lattice units: the distance between nodes is $\Delta x=1$ and the time step is $\Delta t=1$.

2.1 Flat interface

In order to apply periodic conditions, we initialize our system as follows: the nematic phase, with directors in the vertical direction, is in the center of the domain with isotropic phase elsewhere. Thus, we simulate two interfaces, but we focus our analysis on the right interface. At the coexistence temperature ($\gamma_{NI}=2.7$) the interface is completely static and there are no spurious velocities. Eq. (3) allows an interface solution of the form MESQUITA1998

[TABLE]

where $S_{n}$ is the nematic order parameter and $\lambda=\sqrt{27K/A_{0}\gamma}$ is the correlation length. At $\gamma=2.7$, $\lambda=2$ and $S_{n}=1/3$. Fig. 1 D shows that the simulation reproduces accurately the interfacial profile.

Next, we change the temperature slightly away from coexistence. In this case, the interface moves with constant velocity, which can be to the left (negative) or to the right (positive). Figs. 1 A and C, show the interfaces at two different temperatures after an interval of time. For $\gamma<\gamma_{NI}$ the domain of the nematic phase shrinks and for $\gamma>\gamma_{NI}$ the nematic domain expands. In the inset of Fig. 1D, we plot the interface velocity at different temperatures. It is approximately linear close to $\gamma_{NI}$ but is deviates from the linear behavior as the shift from coexistence increases. Above a certain shift from $\gamma_{NI}$ one of the two phases becomes unstable and the entire domain becomes nematic or isotropic in the absence of interfacial propagation. Below $\gamma_{<}=2.65$ the system becomes isotropic and above $\gamma_{>}=3.1$ it becomes nematic.

2.2 Circular interface and spurious velocities

Here we estimate the spurious velocities generated at a circular interface. The spurious velocities arise at curved interfaces as a result of the discretization of the velocity space. However, it is difficult to distinguish the spurious velocities from the physical ones in this multiphase model, because curved interfaces are usually unstable. Fig. 2E shows the interface position at three different temperatures. In order to identify the spurious velocities, we iteratively search for the temperature where the interface velociy is of the order of the spurious ones. We found that for $\gamma=2.707681$ the interface velocity at $t=35000$ is approximately $u_{i}\approx 1.89\times 10^{-7}$. Fig. 2A and B show the oder parameter and velocity field of this almost static interface. From the vorticity field, we notice that the there are vortices between the directions of the velocity vectors of the D3Q19 lattice, which is an indication that these are spurious currents. In Fig. 2 D we plot the velocity field for the same system rotated by 20∘. One notes that the intensity of the vortices changed due to the physical component of these velocities (the small interface velocity), but the positions of the vortices are the same confirming that they are mostly due to the velocity discretization. From this analysis, we estimate that the spurious velocities are around $10^{-6}$ in lattice units for this setup. If necessary, one could reduce the spurious velocities by increasing the isotropy of the lattice PhysRevE.73.047701 or using more advanced collision operators 2019arXiv190412948C .

3 Conclusions

We examined the applicability of the hybrid method of lattice Boltzmann and finite differences for the simulation of nematic-isotropic interfaces. The shape and width of the interface is correctly reproduced and the free interface is static at the coexistence temperature. At temperatures slightly different from the coexistence, the interface moves with constant velocity that depends on the temperature. The interface exists only for a certain range of temperatures: for temperatures much lower or higher than the coexistence, one of the two phases becomes unstable and the interface disappears without propagation. The circular interfaces are not static at the coexistence temperature. Thus, we stabilize them by a temperature shift in a way that the interface velocity is comparable to the spurious velocities. This allows us to visualize the profile of theses spurious currents and to estimate their magnitude.

Acknowledgements

We acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under the contracts: PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-028146) and UID/FIS/00618/2019. Margarida Telo da Gama (MTG) would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the program âThe mathematical design of new materialsâ where most of this work was carried out. This program was supported by EP-SRC Grant Number: EP/R014604/1. MTG participation in the program was supported in part by a Simons Foundation Fellowship.