**Open Journal of Fluid Dynamics**

Vol.08 No.04(2018), Article ID:88315,18 pages

10.4236/ojfd.2018.84022

Numerical Simulation of Liquid Crystal Flow Induced by Annihilation of a Pair of Defects

Shigeomi Chono, Tomohiro Tsuji^{ }

Department of Mechanical Engineering, Kochi University of Technology, Kochi, Japan

Copyright © 2018 by authors and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: September 21, 2018; Accepted: November 3, 2018; Published: November 6, 2018

ABSTRACT

We investigate the flow induced by annihilation of a pair of defects in liquid crystals using the Doi theory with the Marrucci-Greco potential, in which the orientation state is described with the orientational distribution function. We have numerically studied both the transient behaviors of two defects with different structures and their velocity field, and estimated the magnitude of the induced velocity. A defect with positive strength moves faster than one with negative strength. The long-range order of the molecular orientation field has a large effect on the annihilation time, and the annihilation time is reduced by increasing the long-range order. We find that flows are induced during the annihilation of a pair of defects and that several vortices are generated in the vicinity of the defects. The maximum velocity is predicted to develop spatially between the two defects just after their annihilation in time. In our simulation, the maximum induced velocity reaches an order of 10 μm/s. The induced velocity increases with increasing long range-order and nematic potential strength.

**Keywords:**

Complex Fluid, Liquid Crystal, Defect, Annihilation, Generation of Flow

1. Introduction

At present, liquid crystals are widely used as displays (low-molar-mass liquid crystals) and high-strength engineering plastics (liquid crystalline polymers) by applying their anisotropy of material constants. On the other hand, because a liquid crystal has both solid and liquid properties, it has a possibly much wider range of applications, as do solid and liquid. Therefore, studies on new applications of liquid crystals are needed.

When a nematic liquid crystal is observed under a polarizing microscope, a region in which molecules orient parallel or perpendicular to the direction of polarization appears as a dark field, a region in which molecules orient ±45˚ to the direction of polarization appears as a bright field, and a region in which molecules orient in the other direction appears as a gray field. Figure 1(a) is a photograph of a nematic state right after the phase transition from an isotropic state of N-(p-methoxybenzylidene)-p’-butyl-aniline (MBBA), a typical low-molar-mass nematic liquid crystal. We observe the Schlieren texture peculiar to liquid crystals. There are some points in Figure 1(a) where black lines intersect with each other. These points are called defects [1] [2] , around which spatial distortions of molecular orientation fields are generated. Figure 1(b) shows the representative defect orientation states. At the defect cores, the direction of orientation is discontinuous, and the orientation state is singular.

The defects often generated in manufacturing liquid crystal products cause degradation of the productivity and performance of the products. However, it is experimentally known that a pair of defects such as those in Figure 1(b), which have different molecular orientations, attract and finally annihilate each other [2] . This defect annihilation changes the orientation direction, so that a flow is expected to be generated [3] - [8] . Applying thermal energy induces defect generation because the orientation order at defects is lower. Therefore, generating defects artificially at arbitrary locations by applying thermal energy would inevitably generate flows there, leading to the development of a new type of microactuator that converts thermal energy into kinetic energy.

A theory that is able to describe the molecular orientation state satisfactorily should include the effects of long-range and short-range molecular order as well as the effect of flow. Two theories, the Leslie-Ericksen theory [9] [10] and the Doi theory [11] [12] , are well-known and have been widely used for predicting the dynamics of liquid crystals. From the aforementioned viewpoint, however, the Leslie-Ericksen theory, which does not include the short-range order effect, and the Doi theory, which does not account for the long-range order effect, are out as candidates. Marrucci and Greco [13] have expanded the Maier-Saupe potential [14] used in the Doi theory to include both the short-range and long-range

Figure 1. Liquid crystal defect structures. (a) A photograph of a nematic state of MBBA. (b) Molecular orientation configurations around defect cores: Ellipses represent the molecules, and small circles are the defect cores.

order effects and applied it to the original Doi theory [15] . Feng et al. [16] have proposed an expression for the stress tensor and to complete the theory.

Our objectives are to study the flow induced by the annihilation of a pair of defects and to estimate its magnitude using the aforementioned theory as a constitutive equation. A number of studies have simulated the liquid crystal flow during the annihilation of paired defects [17] [18] [19] [20] . Experiments for the moving speed of defects under electric fields have also been performed [21] [22] . However, they do not aim to develop actuators but focus on the movement of defects rather than the liquid crystal flow.

To calculate the orientational distribution functions (ODFs) at many points in a physical space, we have to account for both computational accuracy and time. In this paper, we: 1) approximate the ODF with a series of spherical harmonic functions, 2) study the minimum number of expanded terms required to simulate the orientation state properly, and 3) finally estimate the magnitude of the induced flow during the annihilation of a pair of defects to obtain useful data that can contribute to developing new actuators.

2. Computations

2.1. Governing Equations

The evolution equation for the ODF f is written as

$\frac{\partial f}{\partial t}=\stackrel{\xaf}{D}{\nabla}_{u}\cdot \left({\nabla}_{u}f+\frac{f{\nabla}_{u}V}{kT}\right)-{\nabla}_{u}\cdot \left\{f\left(\kappa \cdot u-\kappa :uuu\right)\right\}$ . (1)

Here, t is the time, k the Boltzmann constant, T the absolute temperature, u the unit vector parallel to a liquid crystal molecule, ${\nabla}_{u}$ the differential operator on a unit sphere, and κ the velocity gradient tensor. $\stackrel{\xaf}{D}$ and V are the rotary diffusivity and Marrucci-Greco nematic potential [13] expressed as

$\stackrel{\xaf}{D}=D{\left(1-\frac{3}{2}S:S\right)}^{-2}$ , (2)

and

$V=-\frac{3}{2}UkT\left(S+\frac{{l}_{i}^{2}}{24}{\nabla}^{2}S\right):uu$ (3)

where D is the rotary diffusivity in an isotropic state, $\nabla $ the differential operator in physical space, and S the order parameter tensor defined as

$S={\displaystyle {\int}_{\left|u\right|=1}\left(uu-\frac{\delta}{3}\right)f\text{d}\Omega}\equiv \langle uu-\frac{\delta}{3}\rangle $ . (4)

δ is the unit tensor.

The conservation equations for the isothermal slow flow of liquid crystals are the continuity and linear momentum equations

$\nabla \cdot v=0$ , (5)

and

$\rho \frac{\partial v}{\partial t}=-\nabla p+\nabla \cdot \tau $ , (6)

where v is the velocity vector, ρ is the fluid density, and p is the pressure. τ is the extra stress tensor derived by Feng et al. [16] expressed as

$\begin{array}{l}\tau =3ckT[A-U\left(A\cdot A-A:Q\right)-\frac{U{l}_{i}^{2}}{24}\{A\cdot {\nabla}^{2}A-Q:{\nabla}^{2}A\underset{}{\overset{\stackrel{}{}}{{\displaystyle}}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\left(\nabla A\right):{\left(\nabla A\right)}^{{\rm T}}}{4}-\frac{\nabla \nabla A:A}{4}\}]+\frac{c}{2}\zeta \text{\hspace{0.05em}}\kappa :Q\end{array}$ (7)

Here, c is the number density of liquid crystal molecules, U the dimensionless nematic potential intensity, l_{i} a parameter indicating the long-range order effect of molecules, ζ the drag coefficient of rotary molecules. A and Q are the second and fourth moments of the ODF f, respectively, expressed as

$A={\displaystyle {\int}_{\left|u\right|=1}uu}f\text{d}\Omega \equiv \langle uu\rangle $ , (8)

and

$Q={\displaystyle {\int}_{\left|u\right|=1}uuuu}f\text{d}\Omega \equiv \langle uuuu\rangle $ . (9)

2.2. Computational Procedure and Boundary Conditions

Let us consider a two-dimensional square area with a side length of H, shown in Figure 2, and put a pair of defects with different orientational states on two points whose coordinates are P (0.3H, 0.5H) and Q (0.7H, 0.5H). For the coordinate system in Figure 2, the components of the velocity gradient tensor κ and order parameter tensor S are

$\kappa =\left|\begin{array}{ccc}\frac{\partial u}{\partial x}& 0& \frac{\partial u}{\partial z}\\ 0& 0& 0\\ \frac{\partial w}{\partial x}& 0& \frac{\partial w}{\partial z}\end{array}\right|$ , (10)

and

$S=\left|\begin{array}{ccc}{S}_{xx}& 0& {S}_{xz}\\ 0& {S}_{yy}& 0\\ {S}_{zx}& 0& {S}_{zz}\end{array}\right|$ , (11)

where u and w are the x and z components of the velocity vector v.

For computation of the orientation field, we substitute Equation (10) and Equation (11) into Equation (1) to obtain the ODF f. In this study, we approximate f with a finite series of spherical harmonic functions Y_{lm}(u) [23] [24] [25] .

$f\left(t,u,x,z\right)={\displaystyle \underset{l=0}{\overset{{l}_{\mathrm{max}}}{\sum}}{\displaystyle \underset{m=-l}{\overset{l}{\sum}}{C}_{lm}\left(t,x,z\right){Y}_{lm}\left(u\right)}}$ . (12)

Here, C_{lm}(t, x, z) are coefficients, and l_{max} is the maximum of the azimuthal quantum number on which the number of terms of the series solution depends.

Figure 2. Flow geometry and coordinate system. The coordinates of a pair of defects are P (0.3H, 0.5H) and Q (0.7H, 0.5H). A defect with positive strength exists at point P and a defect with negative strength exists at point Q. Line segments mean the initial distribution of the director (Equation (18)).

Since the head and tail of a liquid crystal molecule are indistinguishable, we have $f\left(t,u\right)=f\left(t,-u\right)$ . From the parity of the spherical harmonic functions, that is,

${Y}_{lm}\left(-u\right)={\left(-1\right)}^{l}{Y}_{lm}\left(u\right)$ , (13)

the expression (?1)^{l} =1 is obtained, which restricts l to even values. We have non-dimensionalized Equation (1) with 1/D, multiplied the resulting equation by the complex conjugate of spherical harmonic functions,
${Y}_{lm}^{*}\left(={\left(\u20131\right)}^{m}{Y}_{lm}\right)$ , and integrated over the unit sphere to get the ordinary differential equations with respect to C_{lm} (see Appendix 1). The orthonormality of the spherical harmonic functions [26]

$\int}_{\left|u\right|=1}{Y}_{lm}{Y}_{{l}^{\prime}{m}^{\prime}}^{*}\text{d}\Omega}={\delta}_{l{l}^{\prime}}{\delta}_{m{m}^{\prime$ (14)

has been used.

To compute the velocity field, we eliminate p from Equation (6) to obtain the vorticity transport equation

$\rho \frac{\partial \omega}{\partial t}=\frac{{\partial}^{2}{N}_{1}}{\partial x\partial z}-\frac{{\partial}^{2}{\tau}_{xz}}{\partial {x}^{2}}+\frac{{\partial}^{2}{\tau}_{zx}}{\partial {z}^{2}}$ , (15)

where
$\omega \left(=\partial u/\partial z-\partial w/\partial x\right)$ is the vorticity, and N_{1} (=τ_{xx} − τ_{zz}) is the first normal stress difference. Stresses N_{1}, τ_{xz}, and τ_{zx} are explained in the Appendix 2.

We have used the finite-difference scheme to discretize the equations and the implicit Euler method for time integration. Equation (15) is solved using an iterative procedure with the convergence criterion that the average relative error at each node is less than 10^{−5}. We use periodic boundary conditions.

Values of the physical quantities are the fluid density ρ = 10^{3} kg/m^{3}, the absolute temperature T = 320 K, the rotary diffusivity in an isotropic state D = 5.2 × 10^{3} s^{−1}, the number density of molecules c = 2.25 × 10^{24} m^{−3}, the drag coefficient of rotary molecules ζ = 8.89 × 10^{−24} kgm^{2}/s, and the side length of the computational domain H = 1 μm. The computational parameters we select are the nematic potential intensity U = 5.0, 5.5, and 6.0, and the long-range order effect
${l}_{i}^{\ast}$ (=l_{i}/H) = 0.02 − 0.1 with a step of 0.01. The choice of l_{max} is important, and we have determined it by accounting for both the computational accuracy and load. We report the details in the following chapter.

The mesh size is set to
$\Delta {x}^{*}\left(=\Delta x/H\right)=\Delta {z}^{*}\left(=\Delta z/H\right)={10}^{-2}$ , and the time step Δt^{*} (=ΔtD) is varied depending on
${l}_{i}^{\ast}$ ; for example, Δt^{*} = 5 × 10^{−4} when
${l}_{i}^{\ast}=0.1$ .

2.3. Initial Values

The initial velocity vector v is 0. The initial values of C_{lm} are obtained as follows: We multiply Equation (12) by
${Y}_{lm}^{*}$ , integrate over the unit sphere, and use Equation (14) to get

${C}_{lm}\left(0,x,z\right)={\displaystyle {\int}_{\left|u\right|=1}f\left(0,u,x,z\right){Y}_{lm}^{*}\left(u\right)\text{d}\Omega}$ . (16)

Assuming that f is in an equilibrium state (no flow) at t^{*} = 0, f is expressed as [12]

$f\left(0,u,x,z\right)=\mathrm{exp}\left(-\frac{V}{kT}\right)/{\displaystyle {\int}_{\left|u\right|=1}\mathrm{exp}\left(-\frac{V}{kT}\right)\text{d}\Omega}$ . (17)

We use the denominator to normalize f in Equation (17). Since f has uniaxial symmetry in the absence of both flow and external field, the order parameter tensor, Equation (4), is rewritten as

$S=S\left(nn-\frac{\delta}{3}\right)$ , (18)

where n is a unit vector describing the average local molecular orientation, called the director, and S is the scalar order parameter ranging from 0 in a random orientation state to 1 in a perfect alignment, defined as

$S=\sqrt{\frac{3}{2}S:S}$ . (19)

The symbol “:” means the double dot product of two tensors. Substituting Equation (3) and Equation (18) into Equation (17), expanding exponential terms into a power series, and expressing the power by the spherical harmonic functions give

${C}_{lm}\left(0,x,z\right)={\displaystyle \underset{p=0}{\overset{\infty}{\sum}}\frac{{\left(\frac{3}{2}US\right)}^{p}\left(2p\right)!}{p!\left(2p-l\right)!!\left(2p+l+1\right)!!}{Y}_{lm}^{*}\left(n\right)}/{\displaystyle \underset{q=0}{\overset{\infty}{\sum}}\frac{{\left(\frac{3}{2}US\right)}^{q}}{\left(2q+1\right)q!}}$ . (20)

where $\left(2n\right)!!=\left(2n\right)\left(2n-2\right)\cdots 4\cdot 2$ , $\left(2n+1\right)!!=\left(2n+1\right)\left(2n-1\right)\cdots 3\cdot 1$ , 0!! = 1 (n is 0 or even number). In Equation (20), n and S are functions of x and z. We have determined both values (see Appendix 3).

3. Results and Discussion

3.1. Determination of l_{max}

Before computing the velocity and orientation fields, we must determine the value of l_{max}. Let us consider a liquid crystal to which simple shear flow is applied. Setting l_{i} = 0, Equation (3) reduces to

$V=-\frac{3}{2}UkTS:uu$ , (21)

which is the Maier-Saupe nematic potential [14] . The coordinate system is shown in Figure 3. The flow is in the x direction, and the velocity gradient is in the z direction. The polar and azimuthal angles of the unit vector u are θ and ϕ, respectively. For simple shear flow, the velocity gradient tensor is expressed as

$\kappa =\left|\begin{array}{ccc}0& 0& \stackrel{\dot{}}{\gamma}\\ 0& 0& 0\\ 0& 0& 0\end{array}\right|$ , (22)

where
$\stackrel{\dot{}}{\gamma}$ is the shear rate. Equation (11), Equation (21), Equation (22) are substituted into Equation (1) to obtain the ODF f. Equation (20) can be used as an initial value of C_{lm}, but it is a function only of time t. At t^{*} = 0 the director n is assumed to orient in the x direction, so that we set θ_{m} = π/2 and ϕ_{m} = 0, where θ_{m} and ϕ_{m} are the polar and azimuthal angles of the director, respectively.

The azimuthal angle of the director ϕ_{m} is always 0 because the director is in the x-z plane in Figure 3. The polar angle θ_{m} is obtained from S as follows:

$\mathrm{tan}2{\theta}_{m}=\frac{2{S}_{xz}}{{S}_{zz}-{S}_{xx}}$ . (23)

Figure 4 shows the transient behaviors of θ_{m} and S at l_{max} = 4 − 10 (no convergent solution was obtained at l_{max} = 2) when the potential intensity U = 5 and the dimensionless shear rate
${\stackrel{\dot{}}{\gamma}}^{*}\left(=\stackrel{\dot{}}{\gamma}/D\right)=1$ and 40. In Figure 4(a), the director rotates and the order parameter oscillates about the equilibrium value S = 0.615. The behaviors at l_{max} = 8 are almost the same as those at l_{max} = 10, and the behaviors at l_{max} = 6 have smaller periods compared with those at l_{max} = 8 and 10. At l_{max} = 4, steady-state instead of periodic behaviors are obtained. In Figure 4(b), the director and the order parameter reach their steady state values in a short time at l_{max} = 6, 8, and 10, whereas the predictions at l_{max} = 4 are even qualitatively different from the results at other l_{max}. Thus, l_{max} = 4 is not available, l_{max} = 6 is qualitatively acceptable, and l_{max} ≥ 8 is quantitatively

Figure 3. Coordinate system.

Figure 4. Transient behavior of preferred angle θ_{m} and order parameter S for the potential intensity U = 5. No convergent solution was obtained at l_{max} = 2. The lines at l_{max} = 10 overlap with those at l_{max} = 8. (a)
${\stackrel{\dot{}}{\gamma}}^{*}=1$ . The director rotates and the order parameter oscillates about the equilibrium value S = 0.615; (b)
${\stackrel{\dot{}}{\gamma}}^{*}=40$ . Both the director and the order parameter reach their steady state values.

satisfactory. Expansion terms required to approximate the ODF depend on the potential intensity and shear rate. When U and ${\stackrel{\dot{}}{\gamma}}^{*}$ are large, more terms are necessary because the ODF becomes steep and non-uniaxial.

3.2. Orientation Fields

Since the flows induced by annihilation of a pair of defects are not large, and the selected potential intensity U in the present study is not high (5, 5.5, 6), we have used l_{max} = 6 for the following computations.

Figure 5 shows a sequence of predictions of transmitted light intensity observed under crossed nicols for ${l}_{i}^{\ast}=0.\text{1}$ . The transmitted light intensity I is evaluated using

$I={I}_{0}{\mathrm{sin}}^{2}\left(2{\theta}_{m}\right)$ . (24)

Here, I_{0} is the incident light intensity, and θ_{m} is the polar angle of the director obtained by Equation (23). A pair of defect cores attracts each other with time, and brushes connecting the defects become short. Finally the defects are annihilated at t^{*} = 3.85. Since the initial orientation profile is symmetric with respect to z^{*} = 0.5, the defect cores move along this line. After finishing the annihilation, the orientation state does not become homogeneous immediately, but slightly gray areas are discernible at t^{*} = 5. A completely homogeneous orientation state is achieved at t^{*} = 6.05.

Figure 6 shows the profiles of S for the same parameters and times as

Figure 5. A time series of the molecular orientation fields at
${l}_{i}^{\ast}=0.\text{1}$ . The transmitted light intensity is numerically predicted. The horizontal axis is x^{*} (=x/H) and the vertical axis is z^{*} (=z/H). The polarizing and analyzing directions are ±45˚ with respect to the x^{*} axis. A field where the director orients to the polarizing or analyzing directions is dark, and a field where it orients parallel to the x^{*} or z^{*} axis is bright, while the other field is gray.

Figure 5. Since the order parameters near the defect cores are low, the profile of S looks like an inverted cone. The two inverted cones gradually approach and coalesce laterally with each other as time advances. The two sharp ends merge into a single end at t^{*} = 3.85, which coincides with the time the defects

Figure 6. Time series of the order parameter fields at ${l}_{i}^{\ast}=0.\text{1}$ .

annihilate in Figure 5. Just after the annihilation, we can see slight inhomogeneity in S.

Figure 7(a) plots the locations of a pair of defect cores before annihilation for several
${l}_{i}^{\ast}$ . The defect cores approach at a constant speed, but they slow down just before the annihilation. The defect core with the strength s = +1/2 (at point P in Figure 2) moves slightly faster than that with the strength s = −1/2 (at point Q in Figure 2). Thus, the annihilation takes place at a position slightly over x^{*} = 0.5.

Figure 7(b) shows the relationship between the time ${t}_{a}^{\ast}$ required for the

Figure 7. (a) Time series of locations of a pair of defects. Symbols “×” are the positions at which the two defect cores coalesce with each other; (b) Relationship between the time ${t}_{a}^{\ast}$ required for the annihilation process and ${l}_{i}^{\ast}$ .

annihilation process and the long-range order parameter ${l}_{i}^{\ast}$ . A large ${l}_{i}^{\ast}$ gives a short annihilation time because the attractive force between defects becomes strong for large ${l}_{i}^{\ast}$ . For example, ${t}_{a}^{\ast}=190$ at ${l}_{i}^{\ast}=0.02$ , and ${t}_{a}^{\ast}=3.85$ at ${l}_{i}^{\ast}=0.\text{1}$ (refer to Figure 5 and Figure 6); that is, increasing ${l}_{i}^{\ast}$ by a factor of 5 decreases ${t}_{a}^{\ast}$ by a factor of 50. Therefore, ${l}_{i}^{\ast}$ has a large effect on annihilation.

3.3. Velocity Fields

Figure 8 shows the velocity distributions at t^{*} = 1, 3.85, and 5 for
${l}_{i}^{\ast}=0.\text{1}$ as well as the locations of defect cores. An arrow at the bottom of each figure is the reference velocity vector. It is confirmed that flows are induced during annihilation of a pair of defects. Since the molecular orientation state is symmetric with respect to z^{*} = 0.5, the velocity distribution is also symmetric. Flows to the right are observed in the vicinity of z^{*} = 0.5. As a result, a counterclockwise vortex in the upper half region and a clockwise one in the lower half region are generated. Four vortices (two in each half region) are generated at t^{*} = 1 and 3.85. By comparing the velocity vectors near the left and right defect cores at t^{*} = 1, we find that the velocity vector and vortex on the left side are larger, and that the velocity distributions depend on the defect structure. Both right and left vortices rotate in the same direction at t^{*} = 1, so that they coalesce into a larger vortex at t^{*} = 3.85. Since the induced velocity vector is maximum spatially on the line connecting the defect cores (z^{*} = 0.5), it is efficient to use the flow between defects when using the flow induced during their annihilation. We have checked the velocity profiles at times different from those of Figure 8 to verify that the velocity vector has its maximum value when annihilation finishes (t^{*} = 3.85).

We define v_{max} as the maximum velocity value in the computational region at the time when annihilation finishes. Figure 9(a) plots the relationship between v_{max} and
${l}_{i}^{\ast}$ . When
${l}_{i}^{\ast}$ increases, the induced velocity also increases. However, the effect of
${l}_{i}^{\ast}$ on the velocity is not large compared with that on the annihilation

Figure 8. Velocity profiles for ${l}_{i}^{\ast}=0.\text{1}$ . Red circles are the locations of defect cores.

Figure 9. Relationship between maximum velocity v_{max} and
${l}_{i}^{\ast}$ and U. (a) Effect of
${l}_{i}^{\ast}$ ; (b) Effect of U for
${l}_{i}^{\ast}=0.\text{1}\text{.}$

time
${t}_{a}^{\ast}$ ; the relationship between v_{max} and
${l}_{i}^{\ast}$ is almost linear. Figure 9(b) plots v_{max} against the nematic potential intensity U for
${l}_{i}^{\ast}=0.\text{1}$ . The induced velocity increases with increasing U. We explain this result as follows: Since annihilation stems mainly from the spatial gradient of S, an increase in U corresponds to an increase in S, resulting in a steep spatial gradient for S. Thus, a liquid crystal with higher liquid crystallinity (a low-temperature liquid crystal) may generate a faster flow.

4. Conclusion

In this study we have predicted flows induced by the annihilation of a pair of liquid crystal defects using the Doi theory with the Marrucci-Greco potential and the constitutive equation of Feng et al. The long-range order effect on the time required for the annihilation process of a pair of defects is remarkable; when the long-range order is large, the annihilation time becomes short. We have shown that a flow is induced by the annihilation and that several vortices are generated in the vicinity of the defects. The maximum flow is obtained on the line connecting the two defect cores in space, and at the time, the annihilation is just finished. The maximum value of the induced velocity is on the order of 10 μm/s in our study. The induced velocity becomes large when the long-range order and nematic potential strength are high.

Acknowledgements

This work was partially supported by the Japan Society for the Promotion of Science KAKENHI (Grant No. 25289035).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Chono, S. and Tsuji, T. (2018) Numerical Simulation of Liquid Crystal Flow Induced by Annihilation of a Pair of Defects. Open Journal of Fluid Dynamics, 8, 343-360. https://doi.org/10.4236/ojfd.2018.84022

References

- 1. de Gennes, P.G. and Prost, J. (1993) The Physics of Liquid Crystals. 2nd Edition, Clarendon Press, Clarendon.
- 2. Chandrasekhar, S. (1992) Liquid Crystals. 2nd Edition, Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511622496
- 3. Brochard, F. (1973) Backflow Effects in Nematic Liquid Crystals. Molecular Crystals and Liquid Crystals, 23, 51-58. https://doi.org/10.1080/15421407308083360
- 4. Pieranski, P., Brochard, F. and Guyon, E. (1973) Static and Dynamic Behavior of a Nematic Liquid Crystal in a Magnetic Field, Part II: Dynamics. Journal de Physique, 34, 35-48. https://doi.org/10.1051/jphys:0197300340103500
- 5. van Doorn, C.Z. (1975) Dynamic Behavior of Twisted Nematic Liquid-Crystal Layers in Switched Fields. Journal of Applied Physics, 46, 3738-3745. https://doi.org/10.1063/1.322177
- 6. Berreman, D.W. (1975) Liquid-Crystal Twist Cell Dynamics with Backflow. Journal of Applied Physics, 46, 3746-3751. https://doi.org/10.1063/1.322159
- 7. Chono, S. and Tsuji, T. (2008) Proposal of Mechanics of Liquid Crystals and Development of Liquid Crystalline Microactuators. Applied Physics Letters, 92, 051905. https://doi.org/10.1063/1.2840673
- 8. Zhou, Y., Tsuji, T. and Chono, S. (2016) Fundamental Study on the Application of Liquid Crystals to Actuator Devices. Applied Physics Letters, 109, 011902.https://doi.org/10.1063/1.4955267
- 9. Leslie, F.M. (1968) Some Constitutive Equations for Liquid Crystals. Archive for Rational Mechanics and Analysis, 28, 265-283. https://doi.org/10.1007/BF00251810
- 10. Leslie, F.M. (1979) Theory of Flow Phenomena in Liquid Crystals. Advances in Liquid Crystals, 4, 1-81. https://doi.org/10.1016/B978-0-12-025004-2.50008-9
- 11. Doi, M. (1981) Molecular Dynamics and Rheological Properties of Concentrated Solutions of Rodlike Polymers in Isotropic and Liquid Crystalline Phases. Journal of Polymer Sciences: Polymer Physics Edition, 19, 229-243. https://doi.org/10.1002/pol.1981.180190205
- 12. Doi, M. and Edwards, S.F. (1986) The Theory of Polymer Dynamics. Oxford University Press, Oxford.
- 13. Marrucci, G. and Greco, F. (1991) The Elastic Constants of Maier-Saupe Rodlike Molecule Nematics. Molecular Crystals and Liquid Crystals, 206, 17-30. https://doi.org/10.1080/00268949108037714
- 14. Maier, W. and Saupe, A. (1958) Eine einfache molekulare theorie des nematischen kristallinflüssigen zustandes. Zeitschrift für Naturforschung, 13a, 564-566.
- 15. Shimada, T., Doi, M. and Okano, K. (1988) Concentration Fluctuation of Stiff Polymers. III. Spinodal Decomposition. The Journal of Chemical Physics, 88, 7181-7186. https://doi.org/10.1063/1.454370
- 16. Feng, J.J., Sgalari, G. and Leal, L.G. (2000) A Theory for Flowing Nematic Polymers with Orientational Distortion. Journal of Rheology, 44, 1085-1101. https://doi.org/10.1122/1.1289278
- 17. Denniston, C., Tóth, G. and Yeomans, J.M. (2002) Domain Motion in Confined Liquid Crystals. Journal of Statistical Physics, 107, 187-202. https://doi.org/10.1023/A:1014562721540
- 18. Tóth, G., Denniston, C. and Yeomans, J.M. (2002) Hydrodynamics of Topological Defects in Nematic Liquid Crystals. Physical Review Letters, 88, Article ID: 105504. https://doi.org/10.1103/PhysRevLett.88.105504
- 19. Tóth, G., Denniston, C. and Yeomans, J.M. (2003) Hydrodynamics of Domain Growth in Nematic Liquid Crystals. Physical Review E, 67, Article ID: 051705. https://doi.org/10.1103/PhysRevE.67.051705
- 20. Svenšek, D. and Žumer, S. (2002) Hydrodynamics of Pair-Annihilating Disclination Lines in Nematic Liquid Crystals. Physical Review E, 66, Article ID: 021712. https://doi.org/10.1103/PhysRevE.66.021712
- 21. Cladis, P.E., van Saarloos, W., Finn, P.L. and Kortan, A.R. (1987) Dynamics of Line Defects in Nematic Liquid Crystals. Physical Review Letters, 58, 222-225. https://doi.org/10.1103/PhysRevLett.58.222
- 22. Blanc, C., Svenšek, D., Žumer, S. and Nobili, M. (2005) Dynamics of Nematic Liquid Crystal Disclinations: The Role of the Backflow. Physical Review Letters, 95, Article ID: 097802. https://doi.org/10.1103/PhysRevLett.95.097802
- 23. Doi, M. and Edwards, S.F. (1978) Dynamics of Rod-Like Macromolecules in Concentrated Solutions. Part 2. Journal of the Chemical Society, Faraday Transactions 2, 74, 918-932. https://doi.org/10.1039/f29787400918
- 24. Larson, R.G. (1990) Arrested Tumbling in Shearing Flows of Liquid Crystal Polymers. Macromolecules, 23, 3983-3992. https://doi.org/10.1021/ma00219a020
- 25. Larson, R.G. and Öttinger, H.C. (1991) Effect of Molecular Elasticity on Out-of-Plane Orientations in Shearing Flows of Liquid-Crystalline Polymers. Macromolecules, 24, 6270-6282. https://doi.org/10.1021/ma00023a033
- 26. Sasaki, R. (1996) Physical Mathematics. Baifukan. (In Japanese)

Appendix 1

The time evolution equation for coefficients C_{lm} is expressed as

$\begin{array}{c}\frac{\partial {C}_{lm}}{\partial t}=-u\frac{\partial {C}_{lm}}{\partial x}-w\frac{\partial {C}_{lm}}{\partial z}-l\left(l+1\right){\left(1-{S}^{2}\right)}^{-2}{C}_{lm}-\frac{3}{4}U{\left(1-{S}^{2}\right)}^{-2}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\times [\left\{\sqrt{\frac{\text{2\pi}}{15}}\left({D}_{22}+{D}_{2-2}\right)-\sqrt{\frac{\text{4\pi}}{45}}{D}_{20}\right\}{\left(-1\right)}^{m}\sqrt{\frac{\text{8\pi}}{15}}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}\{A\left(\Vert l,-m+1;2,1;{l}^{\prime},{m}^{\prime}\Vert -\Vert l,-m+1;2,-1;{l}^{\prime},{m}^{\prime}\Vert \right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-B\left(\Vert l,-m-1;2,1;{l}^{\prime},{m}^{\prime}\Vert -\Vert l,-m-1;2,-1;{l}^{\prime},{m}^{\prime}\Vert \right)-2m\left(\Vert l,-m;2,2;{l}^{\prime},{m}^{\prime}\Vert -\Vert l,-m;2,-2;{l}^{\prime},{m}^{\prime}\Vert \right)\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left\{\sqrt{\frac{\text{2\pi}}{15}}\left({D}_{22}+{D}_{2-2}\right)+\sqrt{\frac{\text{4\pi}}{45}}{D}_{20}\right\}{\left(-1\right)}^{m}\sqrt{\frac{\text{8\pi}}{15}}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}\{-A\left(\Vert l,-m+1;2,1;{l}^{\prime},{m}^{\prime}\Vert +\Vert l,-m+1;2,-1;{l}^{\prime},{m}^{\prime}\Vert \right)}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-B\left(\Vert l,-m-1;2,1;{l}^{\prime},{m}^{\prime}\Vert +\Vert l,-m-1;2,-1;{l}^{\prime},{m}^{\prime}\Vert \right)+2m\left(\Vert l,-m;2,2;{l}^{\prime},{m}^{\prime}\Vert -\Vert l,-m;2,-2;{l}^{\prime},{m}^{\prime}\Vert \right)\}\end{array}$

$\begin{array}{l}\text{\hspace{0.05em}}+\sqrt{\frac{\text{16\pi}}{45}}{D}_{20}{\left(-1\right)}^{m}\sqrt{\frac{\text{32\pi}}{15}}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}\left(A\Vert l,-m+1;2,-1;{l}^{\prime},{m}^{\prime}\Vert +B\Vert l,-m-1;2,1;{l}^{\prime},{m}^{\prime}\Vert \right)}\\ \text{\hspace{0.05em}}-2\sqrt{\frac{\text{2\pi}}{15}}\left({D}_{21}-{D}_{2-1}\right){\left(-1\right)}^{m}\sqrt{\frac{\text{8\pi}}{15}}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}\{A\Vert l,-m+1;2,-2;{l}^{\prime},{m}^{\prime}\Vert -B\Vert l,-m-1;2,2;{l}^{\prime},{m}^{\prime}\Vert}\\ \text{\hspace{0.05em}}+m\left(\Vert l,-m;2,1;{l}^{\prime},{m}^{\prime}\Vert +\Vert l,-m;2,-1;{l}^{\prime},{m}^{\prime}\Vert \right)-\sqrt{\frac{3}{2}}\left(A\Vert l,-m+1;2,0;{l}^{\prime},{m}^{\prime}\Vert -B\Vert l,-m-1;2,0;{l}^{\prime},{m}^{\prime}\Vert \right)\}]\\ \text{\hspace{0.05em}}-\frac{\partial u}{\partial x}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}(\sqrt{\frac{\text{\pi}}{30}}{A}^{\prime}{\Vert l,m;2,1;{l}^{\prime},{m}^{\prime}-1\Vert}^{*}-\sqrt{\frac{\text{\pi}}{30}}{B}^{\prime}{\Vert l,m;2,1;{l}^{\prime},{m}^{\prime}+1\Vert}^{*}-\sqrt{\frac{\text{\pi}}{30}}{A}^{\prime}{\Vert l,m;2,-1;{l}^{\prime},{m}^{\prime}-1\Vert}^{*}}\\ \text{\hspace{0.05em}}+\sqrt{\frac{\text{\pi}}{30}}{B}^{\prime}{\Vert l,m;2,-1;{l}^{\prime},{m}^{\prime}+1\Vert}^{*}-\sqrt{\frac{\text{2\pi}}{15}}{m}^{\prime}{\Vert l,m;2,2;{l}^{\prime},{m}^{\prime}\Vert}^{*}+\sqrt{\frac{\text{2\pi}}{15}}{m}^{\prime}{\Vert l,m;2,-2;{l}^{\prime},{m}^{\prime}\Vert}^{*}\end{array}$

$\begin{array}{l}\text{\hspace{0.05em}}+\sqrt{\frac{\text{4\pi}}{5}}{\Vert l,m;2,0;{l}^{\prime},{m}^{\prime}\Vert}^{*}-\sqrt{\frac{\text{6\pi}}{5}}{\Vert l,m;2,2;{l}^{\prime},{m}^{\prime}\Vert}^{*}-\sqrt{\frac{\text{6\pi}}{5}}{\Vert l,m;2,-2;{l}^{\prime},{m}^{\prime}\Vert}^{*})\\ \text{\hspace{0.05em}}-\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right){\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}\{\sqrt{\frac{\text{4\pi}}{45}}{B}^{\prime}{\Vert l,m;2,0;{l}^{\prime},{m}^{\prime}+1\Vert}^{*}-\sqrt{\frac{\text{4\pi}}{45}}{A}^{\prime}{\Vert l,m;2,0;{l}^{\prime},{m}^{\prime}-1\Vert}^{*}}\\ \text{\hspace{0.05em}}+\sqrt{\frac{\text{2\pi}}{15}}{m}^{\prime}\left({\Vert l,m;2,1;{l}^{\prime},{m}^{\prime}\Vert}^{*}+{\Vert l,m;2,-1;{l}^{\prime},{m}^{\prime}\Vert}^{*}\right)+\sqrt{\frac{\text{6\pi}}{5}}\left({\Vert l,m;2,1;{l}^{\prime},{m}^{\prime}\Vert}^{*}-{\Vert l,m;2,-1;{l}^{\prime},{m}^{\prime}\Vert}^{*}\right)\}\\ \text{\hspace{0.05em}}-\left(\frac{1}{6}\frac{\partial u}{\partial z}-\frac{1}{3}\frac{\partial w}{\partial x}\right)\left(A{C}_{lm-1}-B{C}_{lm+1}\right)\text{\hspace{0.05em}}+\frac{\partial w}{\partial z}{\displaystyle \underset{{l}^{\prime}{m}^{\prime}}{\sum}{C}_{{l}^{\prime}{m}^{\prime}}(\sqrt{\frac{\text{2\pi}}{15}}{A}^{\prime}{\Vert l,m;2,1;{l}^{\prime},{m}^{\prime}-1\Vert}^{*}}\\ +\sqrt{\frac{\text{2\pi}}{15}}{B}^{\prime}{\Vert l,m;2,-1;{l}^{\prime},{m}^{\prime}+1\Vert}^{*}+\sqrt{\frac{\text{16\pi}}{5}}{\Vert l,m;2,0;{l}^{\prime},{m}^{\prime}\Vert}^{*})\end{array}$

where

${D}_{lm}={C}_{lm}+\frac{{l}_{i}^{2}}{24}\left(\frac{{\partial}^{2}{C}_{lm}}{\partial {x}^{2}}+\frac{{\partial}^{2}{C}_{lm}}{\partial {z}^{2}}\right)={C}_{lm}+\frac{{l}_{i}^{2}}{24}{\nabla}^{2}{C}_{lm}$

$\Vert \text{\hspace{0.17em}}{l}_{1},{m}_{1};{l}_{2},{m}_{2};{l}_{3},{m}_{3}\Vert ={\displaystyle {\int}_{\left|u\right|=1}{Y}_{{l}_{1}{m}_{1}}{Y}_{{l}_{2}{m}_{2}}{Y}_{{l}_{3}{m}_{3}}\text{d}\Omega}$

${\Vert \text{\hspace{0.17em}}{l}_{1},{m}_{1};{l}_{2},{m}_{2};{l}_{3},{m}_{3}\Vert}^{*}={\displaystyle {\int}_{\left|u\right|=1}{Y}_{{l}_{1}{m}_{1}}^{*}{Y}_{{l}_{2}{m}_{2}}{Y}_{{l}_{3}{m}_{3}}\text{d}\Omega}$

${S}^{2}=\frac{3}{2}S:S$

$A=\sqrt{\left(l+m\right)\left(l-m+1\right)}$

$B=\sqrt{\left(l-m\right)\left(l+m+1\right)}$

${A}^{\prime}=\sqrt{\left({l}^{\prime}+{m}^{\prime}\right)\left({l}^{\prime}-{m}^{\prime}+1\right)}$

and

${B}^{\prime}=\sqrt{\left({l}^{\prime}-{m}^{\prime}\right)\left({l}^{\prime}+{m}^{\prime}+1\right)}$

The non-zero components of the order parameter tensor S and fourth-order tensor Q are expressed in terms of C_{lm} as

${S}_{xx}=\sqrt{\frac{\text{2\pi}}{15}}\left({C}_{22}+{C}_{2-2}\right)-\sqrt{\frac{\text{4\pi}}{45}}{C}_{20}$

${S}_{yy}=-\sqrt{\frac{\text{2\pi}}{15}}\left({C}_{22}+{C}_{2-2}\right)-\sqrt{\frac{\text{4\pi}}{45}}{C}_{20}$

${S}_{zz}=\sqrt{\frac{\text{16\pi}}{45}}{C}_{20}$

${S}_{xz}=-\sqrt{\frac{\text{2\pi}}{15}}\left({C}_{21}-{C}_{2-1}\right)$

$\begin{array}{c}{Q}_{xxxx}=\sqrt{\frac{\text{2\pi}}{315}}\left({C}_{44}+{C}_{4-4}\right)-\sqrt{\frac{\text{8\pi}}{2205}}\left({C}_{42}+{C}_{4-2}\right)+\sqrt{\frac{\text{4\pi}}{1225}}{C}_{40}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\sqrt{\frac{\text{72\pi}}{735}}\left({C}_{22}+{C}_{2-2}\right)-\sqrt{\frac{\text{16\pi}}{245}}{C}_{20}+\sqrt{\frac{\text{4\pi}}{25}}{C}_{00}\end{array}$

${Q}_{xxxz}=-\sqrt{\frac{\text{\pi}}{315}}\left({C}_{43}-{C}_{4-3}\right)+\sqrt{\frac{\text{\pi}}{245}}\left({C}_{41}-{C}_{4-1}\right)-\sqrt{\frac{\text{6\pi}}{245}}\left({C}_{21}-{C}_{2-1}\right)$

${Q}_{xxyy}=-\sqrt{\frac{\text{2\pi}}{315}}\left({C}_{44}+{C}_{4-4}\right)+\sqrt{\frac{\text{4\pi}}{11025}}{C}_{40}-\sqrt{\frac{\text{16\pi}}{2205}}{C}_{20}+\sqrt{\frac{\text{4\pi}}{225}}{C}_{00}$

$\begin{array}{c}{Q}_{xxzz}=\sqrt{\frac{\text{8\pi}}{2205}}\left({C}_{42}+{C}_{4-2}\right)-\sqrt{\frac{\text{64\pi}}{11025}}{C}_{40}+\sqrt{\frac{\text{2\pi}}{735}}\left({C}_{22}+{C}_{2-2}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\sqrt{\frac{\text{4\pi}}{2205}}{C}_{20}+\sqrt{\frac{\text{4\pi}}{225}}{C}_{00}\end{array}$

${Q}_{xyyz}=\sqrt{\frac{\text{\pi}}{315}}\left({C}_{43}-{C}_{4-3}\right)+\sqrt{\frac{\text{\pi}}{2205}}\left({C}_{41}-{C}_{4-1}\right)-\sqrt{\frac{\text{2\pi}}{735}}\left({C}_{21}-{C}_{2-1}\right)$

${Q}_{xzzz}=-\sqrt{\frac{\text{16\pi}}{2205}}\left({C}_{41}-{C}_{4-1}\right)-\sqrt{\frac{\text{6\pi}}{245}}\left({C}_{21}-{C}_{2-1}\right)$

$\begin{array}{c}{Q}_{yyzz}=-\sqrt{\frac{\text{8\pi}}{2205}}\left({C}_{42}+{C}_{4-2}\right)-\sqrt{\frac{\text{64\pi}}{11025}}{C}_{40}-\sqrt{\frac{\text{2\pi}}{735}}\left({C}_{22}+{C}_{2-2}\right)\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\sqrt{\frac{\text{4\pi}}{2205}}{C}_{20}+\sqrt{\frac{\text{4\pi}}{225}}{C}_{00}\end{array}$

and

${Q}_{zzzz}=\sqrt{\frac{\text{256\pi}}{11025}}{C}_{40}+\sqrt{\frac{\text{64\pi}}{245}}{C}_{20}+\sqrt{\frac{\text{4\pi}}{25}}{C}_{00}$

Appendix 2

The first normal stress difference N_{1} and the shear stresses τ_{xz} and τ_{zx} are expressed in terms of S and Q as

$\begin{array}{c}\frac{{N}_{1}}{3ckT}={S}_{xx}-{S}_{zz}-U({S}_{xx}^{2}-{S}_{zz}^{2}+\frac{1}{3}{S}_{xx}-\frac{1}{3}{S}_{zz}-{S}_{xx}{Q}_{xxxx}+{S}_{xx}{Q}_{xxzz}-{S}_{yy}{Q}_{xxyy}\\ \underset{}{\overset{}{}}+{S}_{yy}{Q}_{yyzz}-{S}_{zz}{Q}_{xxzz}+{S}_{zz}{Q}_{zzzz}-2{S}_{xz}{Q}_{xxxz}+2{S}_{xz}{Q}_{xzzz})\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{U{l}_{i}}{24}[{S}_{xx}{\nabla}^{2}{S}_{xx}-{S}_{zz}{\nabla}^{2}{S}_{zz}+\frac{1}{3}{\nabla}^{2}\left({S}_{xx}-{S}_{zz}\right)-\left({Q}_{xxxx}-{Q}_{xxzz}\right){\nabla}^{2}{S}_{xx}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\left({Q}_{xxyy}-{Q}_{yyzz}\right){\nabla}^{2}{S}_{yy}-\left({Q}_{xxzz}-{Q}_{zzzz}\right){\nabla}^{2}{S}_{zz}-2\left({Q}_{xxxz}-{Q}_{xzzz}\right){\nabla}^{2}{S}_{xz}\end{array}$

$\begin{array}{l}\text{\hspace{0.05em}}+\frac{1}{4}\{{\left(\frac{\partial {S}_{xx}}{\partial x}\right)}^{2}-{\left(\frac{\partial {S}_{xx}}{\partial z}\right)}^{2}+{\left(\frac{\partial {S}_{yy}}{\partial x}\right)}^{2}-{\left(\frac{\partial {S}_{yy}}{\partial z}\right)}^{2}+{\left(\frac{\partial {S}_{zz}}{\partial x}\right)}^{2}-{\left(\frac{\partial {S}_{zz}}{\partial z}\right)}^{2}\\ \text{\hspace{0.05em}}+2{\left(\frac{\partial {S}_{xz}}{\partial x}\right)}^{2}-2{\left(\frac{\partial {S}_{xz}}{\partial z}\right)}^{2}-\left(\frac{{\partial}^{2}{S}_{xx}}{\partial {x}^{2}}-\frac{{\partial}^{2}{S}_{xx}}{\partial {z}^{2}}\right){S}_{xx}-\left(\frac{{\partial}^{2}{S}_{yy}}{\partial {x}^{2}}-\frac{{\partial}^{2}{S}_{yy}}{\partial {z}^{2}}\right){S}_{yy}\\ \text{\hspace{0.05em}}-\left(\frac{{\partial}^{2}{S}_{zz}}{\partial {x}^{2}}-\frac{{\partial}^{2}{S}_{zz}}{\partial {z}^{2}}\right){S}_{zz}-2\left(\frac{{\partial}^{2}{S}_{xz}}{\partial {x}^{2}}-\frac{{\partial}^{2}{S}_{xz}}{\partial {z}^{2}}\right){S}_{xz}\}]\\ \text{\hspace{0.05em}}+\frac{\zeta}{6kT}\left\{\frac{\partial u}{\partial x}\left({Q}_{xxxx}-{Q}_{xxzz}\right)+\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)\left({Q}_{xxxz}-{Q}_{xzzz}\right)+\frac{\partial w}{\partial z}\left({Q}_{xxzz}-{Q}_{zzzz}\right)\right\}\end{array}$

$\begin{array}{c}\frac{{\tau}_{xz}}{3ckT}={S}_{xz}-U\left\{\left({S}_{xx}+{S}_{zz}\right){S}_{xz}+\frac{1}{3}{S}_{xz}-{S}_{xx}{Q}_{xxxz}-{S}_{yy}{Q}_{xyyz}-{S}_{zz}{Q}_{xzzz}-2{S}_{xz}{Q}_{xxzz}\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{U{l}_{i}^{2}}{24}\{\underset{\_}{\underset{\_}{{S}_{xx}{\nabla}^{2}{S}_{xz}+{S}_{xz}{\nabla}^{2}{S}_{zz}}}+\frac{1}{3}{\nabla}^{2}{S}_{xz}-{Q}_{xxxz}{\nabla}^{2}{S}_{xx}-{Q}_{xyyz}{\nabla}^{2}{S}_{yy}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}-{Q}_{xzzz}{\nabla}^{2}{S}_{zz}-2{Q}_{xxzz}{\nabla}^{2}{S}_{xz}+\frac{1}{4}(\frac{\partial {S}_{xx}}{\partial x}\frac{\partial {S}_{xx}}{\partial z}+\frac{\partial {S}_{yy}}{\partial x}\frac{\partial {S}_{yy}}{\partial z}+\frac{\partial {S}_{zz}}{\partial x}\frac{\partial {S}_{zz}}{\partial z}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+2\frac{\partial {S}_{xz}}{\partial x}\frac{\partial {S}_{xz}}{\partial z}-\frac{{\partial}^{2}{S}_{xx}}{\partial x\partial z}{S}_{xx}-\frac{{\partial}^{2}{S}_{yy}}{\partial x\partial z}{S}_{yy}-\frac{{\partial}^{2}{S}_{zz}}{\partial x\partial z}{S}_{zz}-2\frac{{\partial}^{2}{S}_{xz}}{\partial x\partial z}{S}_{xz})\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{\zeta}{6kT}\left\{\frac{\partial u}{\partial x}{Q}_{xxxz}+\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right){Q}_{xxzz}+\frac{\partial w}{\partial z}{Q}_{xzzz}\right\}\end{array}$

and

$\frac{{\tau}_{zx}}{3ckT}=\cdots -\frac{U{l}_{i}^{2}}{24}\left\{\underset{\_}{\underset{\_}{{S}_{xz}{\nabla}^{2}{S}_{xx}+{S}_{zz}{\nabla}^{2}{S}_{xz}}}+\cdots \right\}+\cdots $

The difference between τ_{xz} and τ_{zx} is only the underlined term.

Appendix 3

With respect to the director n, let us define the angle between n and the x axis as θ_{c}, so that the angular momentum equation of the Leslie-Ericksen theory with the one-constant approximation of the elastic constants in the molecular field reduces simply to

${\nabla}^{2}{\theta}_{c}=0$ (A1)

in an equilibrium state [2] . When defects are included, a fundamental solution of Equation (A1) is

${\theta}_{c}=s{\mathrm{tan}}^{-1}\left(\frac{z}{x}\right)$ , (A2)

where s is the defect strength. Since Equation (A1) is linear, the superposition of solutions is effective. Thus, when a defect with s = +1/2 exists at P (x_{1}, z_{1}) and a defect with s = −1/2 exists at Q (x_{2}, z_{2}) in Figure 2, the director distribution around the defect cores is

${\theta}_{c}=\frac{1}{2}{\mathrm{tan}}^{-1}\frac{z-{z}_{1}}{x-{x}_{1}}-\frac{1}{2}{\mathrm{tan}}^{-1}\frac{z-{z}_{2}}{x-{x}_{2}}$ . (A3)

Finally, we have modified the values of Equation (A3) so that they fit the periodic boundary condition. We denote the initial distribution of the director by line segments in Figure 2.

For the scalar order parameter, we set S = 0 at defect cores and S = S^{eq} at the other region. S^{eq} depends only on the nematic potential intensity U, and is obtained from Equation (1) without flow terms. For example, S^{eq} = 0.615 at U = 5.