Investigation of Shear-Driven and Pressure-Driven Liquid Crystal Flow at Microscale: A Quantitative Approach for the Flow Measurement

The liquid crystal-based method is a new technology developed for flow visualizations and measurements at microscale with great potentials. It is the priority to study the flow characteristics before implementation of such a technology. A numerical analysis has been applied to solve the simplified dimensionless two-dimensional Leslie–Ericksen liquid crystal dynamic equation. This allows us to analyze the coupling effect of the LC’s director orientation and flow field. We will be discussing two classic shear flow cases at microscale, namely Couette and Poiseuille flow. In both cases, the plate drag speed in the state of Couette flow are varied as well as the pressure gradients in Poiseuille flow state are changed to study their effects on the flow field distributions. In Poiseuille flow, with the increase of applied pressure gradient, the influence of backflow significantly affects the flow field. Results show that the proposed method has great advantages on measurement near the wall boundaries which could complement to the current adopted flow measurement technique. The mathematical model proposed in this article could be of great potentials in the development of the quantitatively flow measurement technology.


Introduction
Thermotropic liquid crystals (LCs) are a mesophase exhibiting the appearance of an anisotropic liquid whose molecules represent in rodlike or discotic shape [1]. Its unique features are the temperature dependent nature, whose phase transitions are controlled by the temperature. For thermotropic LCs, the liquid crystals in nematic phase attracts tremendous attention [2]. Among them, the rod like molecules with large length-to-diameter ratio, and no long range order in the molecular centroid are the most common ones. Due to the presence of the rod-shaped molecules, the fluid translational motion is coupled with the internal orientational motion of molecules. Hence, their flow properties can be much richer than simple Newtonian fluids. Conversely, the change of molecular orientation order affects the flow in return. This coupling has important practical consequences upon the application of LCs [3]. Therefore, investigations on mechanism of hydrodynamic shear induced reorientation of the liquid crystal molecules are necessary. The research of the T.G. Anderson team focuses on the physical change of liquid crystal fluid under microfluid channel with pressure gradient [4]. The team studies the transformation process of the two states of liquid crystal fluid by numerical calculation and compares the results of the transition of the flow state with the experiment of the Sengupta team [5]. The dimensionless equation is obtained through normalization, which makes the results more intuitive and more universal. In addition, the elastic free energy of liquid crystals in different states is presented. The reasons for the change of liquid crystal flow state are well explained by the comparison of the elastic free energy [5]. Compared to the fluid flow at marcoscale, flow at microscale (mostly in microchannels here) is fundamentally different due to the limited inertial effects and dominating viscous stress and interfacial tension.
In this paper, we choose 4-Cyano-4 -pentylbiphenyl (5CB) as which stays in nematic phase under room temperature. 5CB is a type of nematic liquid crystal (NLC) called "calamitics" [6] which represents the rod-like molecular structure. The structure can be observed experimentally using the polarized microscopy configuration. Recent literatures have been focusing more on the characteristic response to the electric field [7][8][9][10]. The well-known orientation Fredericks transition has been extensively studied and applied in practical applications [11]. Furthermore, the combination of the delicate microscale fluidic control with LCs, especially the nematic LCs has allowed the possibility of topological studies [12], phase transitions [13][14][15][16]and the unique stripe structure [17,18]. More importantly, it also leads to the diverse application in fields of microfabrications, 3D printing and bioengineering [12,13,19,20]. Nonetheless, very few works have been reported on the reorientation influenced by the hydrodynamic pressure at microscale [21][22][23]. Currently, the most commonly used model to study directional vector deflection and flow of nematic liquid crystals is the Leslie-Ericksen equations (L-E equation). We apply the Leslie-Ericksen formalism for nematodynamics to investigate the shear flow model. First, we describe the model of the liquid crystal by means of L-E equations [24]. The analysis is restricted to the steady state deformations of the nematic phase where α 3 /α 2 > 0. The flow behavior of nematic liquid crystals can be generally categorized into two major types through the signs of the two Leslie viscosities α 2 and α 3 . If α 3 /α 2 > 0, then the liquid crystal is flow aligning, otherwise it is flow tumbling [25]. Next, we addressed two typical steady 2D cases: the shear-driven (Couette) and the pressure-driven (Poiseuille) flow of nematic liquid crystal in a parallel microchannel respectively. Then we applied justified assumptions for simplification purpose before solving them numerically. The director profile is calculated for various pressure gradients and shear stresses. For pressure-driven flow cases, we find that under the limitation of strong anchoring and weak flow effects, flow alignment is not presented. In fact, the director field is majorly determined by the boundary conditions. For other cases, the results clearly show the influence of flow condition on reorientation of director field, which provide guideline for flow measurements at microscale [26][27][28].

Mathematical Models
There are two major differences considering the hydrodynamics of simple Newtonian fluids and that of the LCs. First, LCs' molecules can be rotated by the pressure gradient due to its unique geometries. Second, the equilibrium free energy of LCs is more complex. This coupling between the elastic energy and the flow, which is called backflow, leads to rich hydrodynamic behaviors. In order to seize the main physical properties of LCs, the details of the molecules are neglected and the ideal approximation is processed. For example, for nematic LCs, LC molecules are often idealized as a long rod with a symmetric head and tail, and hence uniaxial symmetry. As mentioned earlier, LCs molecules are anisotropic, therefore we need to introduce variables called order parameter to describe the alignment and orientation of LCs. We adopt the well accepted Leslie-Ericksen model to describe the director orientation and LC flow. In the continuum theory, we define the orientation of the LCs at a point by a unit vector n called the director.
The full equation of nematodynamics consists of Representing mass, energy, and momentum conservations respectively. We consider a steady 2D plane Poiseuille-Couette flow of nematic liquid crystal in a microfluidic channel with two parallel boundaries and a thin film geometry. The upper plate travels along the x direction at a constant speed V (shear-driven flow) or a constant pressure gradient is applied within the LCs (pressure-driven flow) as shown in Figures 1 and 2 respectively.

Couette Flow
It will be assumed that the director and the velocity take the forms We now investigate the coquette flow at a fixed distance, d = 2h where the lower plate is at rest and the upper plate is moving at a constant velocity V as shown in Figure 1. Three assumptions are listed as follows (Macsithigh, G. P, and P. K. Currie, 1977): The one-constant approximation applied to the Oseen-Zocher-Frank elastic free energy equation. In details K 11 = K 22 = K 33 ≡ K, where K 11 is the splay coefficient, K 22 is the twist coefficient, and K 33 is the bend coefficient. (Currie, P. K., and G. P. Macsithigh, 1979).
Parodi's relation [29] The zero viscosity coefficient α 1 = 0. From the (2), (3) and the assumptions, it follows that For simplicity we nondimensionalize the equations by using the scaling lengths with half-width of the channel h, the velocity v with v = V/2, and all viscosities γ i α i with α 4 . The dimensionless governing equations become: where γ 1 = α 3 − α 2 , γ 2 = α 6 − α 5 , the α i are constant viscosities for the NLC. The director is strongly anchored parallel to the plates at the boundaries (Leslie, F.M. 1979), and that the solution for Φ is symmetric in z. A more detailed derivation process can be refereed to Appendices A and B for Couette flow.
We applied the boundary conditions as:

Poiseuille Flow
Here we consider the steady 2D Poiseuille flow in a microfluidic channel −h ≤ z ≤ h as schemed in Figure 2. To investigate a typical, two-dimensional shear flow behavior, the director and the velocity take the forms In order to simplify the boundary conditions, it is assumed that Φ is the angle between z axis and the molecule axis ( Figure 2). We use the first two assumption α 6 = α 5 + α 3 + α 2 and K1 = K3 = K. Therefore, the following formulas can be obtained 2K where To nondimensionalize the equations, we adopted the following strategies: scale the physical length with the half width of the channel h, the velocity v with v = 2Gh 2 /3α 4 (Anderson, T. G. et al., 2015), and all viscosities γ i , α i with α 4 . The dimensionless form of the governing equations become, where g = Gh 3 /K is a dimensionless pressure gradient. It is assumed that the director is strongly anchored vertically and uniform to the plates at the boundaries before the flow turned on and the solution for Φ is symmetric in z.
The director profile across the dimensionless channel width −1≤ z ≤ 1 under strong anchoring and g = 25 is calculated in Figure 3a

Couette Flow
In this section, we conducted a systematic parametric investigation on velocity and director fields of a nematic liquid crystal under the shear controlled by the velocity V of the upper plate. The height of the channel was setup as 10 µm. The computations were carried out for different velocities of the upper plane. Figure 4a shows director profiles when the velocity of the upper plate is 10 µm/s, 100 µm/s, 250 µm/s, and 350 µm/s respectively. It can be seen in Figure 4a that the director profile is almost an axisymmetric parabola at the center. With increasing of the velocity, the angle of the director becomes smaller which is more flow orientated. In Figure 4b, the mid-plane angles ∅ m are plotted as functions of velocity V. We found that the mid-plane angles ∅ m tend to ∅ c =arctan[(α 3 /α 2 ) 0.5 ] (∅ c = 0.22, in this article) with the increasing velocity, which agrees with the continuum theory proposed by Leslie and Ericksen [30]. Leslie defined the Leslie angle as: "for which there is no hydrodynamic torque on the director in simple shear flow of an infinitely thick sample". It is found in this article that the Leslie angle still exists in microchannel. Figure 5 shows the solutions at Poiseuille flow with the change of dimensionless pressure gradient g. It can be found that for both weak flow and strong flow anchoring boundary conditions the solutions exist for all flow rates ( Figure 5). From T.G. Anderson's founding on the elastic energies of the two solutions, there is a critical dimensional pressure gradient g*. The energy is lower when the value of g is lower. Therefore, it explains that for cases of g < g*, the weak flow would occur, while for the cases of g > g* the strong flow solution would be expected. In this part, g* is around 40 for strong anchoring boundary condition.   Figure 6 shows weak flow solution at g = 10 and g = 25. We find out that the solution to the anchoring-dominated case has two regions symmetrically placed around the channel centerline. Near the walls and at the centerline, the director follows the strong homeotropic anchoring condition. The velocity profiles are approximately parabolic. The greater the dimensionless pressure gradient is, the greater the velocity and the perturbation of the director are. Strong flow solutions at g = 50 are plotted in Figure 7, where we see more complex phenomena. The velocity front is a sharper parabolic shape and the director profile is also symmetrical. This differs from the weak flow, where the direction filed is mostly aligned with the flow, regulated by the homeotropic anchoring only near the walls. Therefore, when the pressure gradient further increases, the nematic profile evolves into a flow-aligned state where the director field is majorly flow-aligned. For weak flow, the director orientates normally to the flow orientation at the center w while for strong flow, the director aligns along the flow direction at the center region. Above, we discussed two kinds of boundary conditions. However, during numerical investigations, a third occurrence may happen, arising from the boundary condition Φ(−1) = π instead of the condition Φ(−1) = −π. The boundary conditions drive the director to rotate through an angle + π across the channel, but the flow near the wall pushes the director in the opposite direction.

Discussion
We have calculated the derivative of the directional profile with respect to the channel positions, which represents the sensitivity of the proposed method for flow measurement, as shown in Figure 8. The results indicate that for both the Couette flow and Poiseuille flow, the maximum sensitivity always occurs near the wall boundary (that is the two plates), where the traditionally adopted tracing particle-based measurement technique (namely particle image velocimetry) shows vulnerability. In the cases of Couette flow, we can always find the highest sensitivities near the two plates (see Figure 8a,b) and the lowest sensitivity in the middle region. Similarly, we obtained the maximum sensitivities near the two boundaries (Figure 8c,d). However, when it is under the relatively weak pressure gradient for Poiseuille flow (g = 10~30), Figure 8c indicates a good sensitivity can be obtained in the middle region while for a strong flow case (Figure 8d) g = 50), we get a similar conclusion with scenario in (a) and (b) that the method fails to sense the flow in the center region. Overall, in all of the cases, the region near the boundary favors the measurement sensitivity hence offering the privilege over conventional particle based measurement method. We do believe the no slip boundary condition and the homeotropic anchoring condition near the boundary contribute to the large direct field gradient and hence the relatively high sensitivity.

Conclusions
In this paper, we have numerically studied director field and velocity field of twodimensional Poiseuille flow and Couette flow. The director profile was calculated for various pressure gradient and shear stress. By solving simplified nondimensional L-E equations we obtained numerical results of pressure-driven flow and shear-driven flow.
Firstly, for pressure-driven flow, we find that (i) Limited by weak flow effects, flow alignment does not occur and the director orientation is influenced predominantly by the boundary conditions. In such cases, the velocity profile is like a parabola and symmetric about z = 0. (ii) When the pressure gradient arrives at a value, a qualitative change in the director profile is given explicitly.
Flow alignment occurs and orientational boundary layers exist near the planes. For shear-driven flow, (i) the results show that if the velocity is less or equal to the threshold, the director field is similar to a parabola and the velocity field is presented in a linear manner and the mid-plane angles ∅ m tend to ∅ c = arctan[(α 3 /α 2 ) 0.5 ] with the increasing velocity.
(ii) When the velocity exceeds the threshold, the profiles will lose its stability. Vthreshold is about 10 −3 -10 −4 m/s (depend on different parameters). If the velocity exceeds the threshold, the director deviates from the plane of shear. The deformation takes a form of director rotation about the axis perpendicular to the layer plane. As a result, transverse flow arises. The method of nondimensionalizing and the numerical approach proposed in this article will be a potential technique in liquid crystal research. The coupling influence what we analyze is might be able to give clues to the flow dynamics.
Flow measurement by means of the LC reorientation has been proposed recently and in this article, we have theoretically analyzed the viability of such method in flow sensing. The results show that the proposed method has great advantages for sensing near the boundaries which could complement to the technique currently adopted. In this article, we also showed that the orientation angle of the LC follows a non-linear relationship to the flow field distribution. The mathematical model proposed here could be of great potentials when the quantitatively measurement strategy is to be developed for practical applications.

Data Availability Statement:
The data published is available upon reasonable request with the corresponding authors.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
It will be assumed that the director and the velocity take the forms Clearly, the constraints (1) are satisfied and because of the dependence of z on n and v, the governing LE dynamic equations in the absence of any forces F and G become where W F is the usual nematic free energy density. The non-zero contribution to the rate of strain tensor A and vorticity tensor W are where a prime denotes differentiation with respect to z. t ij = α 1 n k A kp n p n j + α 2 N i n j + α 3 n i N j + α 4 A ij + α 5 n j A ik n k + α 6 n i A jk n p (A6) Equation (7) is It is clear from Equations (12) and (14), and the dependence of t 12 upon z, that Inserting (15) into (13) shows that For some constants a and c because the right-hand side of (10) is a function of z only. (taking the derivative of (12) with respect to z and the derivative of (10) with respect x will also lead to the same result.) Therefore It is now seen from (10) and (14) that For some constant p 0 , and hence, by (12), the pressure takes the form Using the form for t 12 given by (7), the result (13) may be formulated conveniently as Φ is the angle between z axis. Equation (3) now has been reduced to (17), the pressure being available via (16) to allow this reduction.
Our attention is now turned to the remaining Equation (4). Using The equations are simplified to K 2 n 1,22 + (K 3 − K 2 )(n 2 2 n 1,2 ) ,2 + g 1 = λn 1 (A19) K 1 n 2,22 + (K 3 − K 2 )[(n 2 2 n 2,2 ) ,2 − n 2 n p,2 n p,2 ] + g 2 = λn 2 (A20) The Lagrange multiplier λ can be eliminated from these two equations by multiplying (20) by n 2 , (21) by n 1 and then subtracting the resulting equations to produce, after much tedious algebra, the single equation Further substitutions for g i and n i using (8) and (1) finally reduce this equation to This last inequality being valid because K 1 and K 3 are necessarily non-negative. The flow v(z) in the z-direction may be induced in a sample of nematic liquid crystal confined between two parallel plates by the horizontal motion of the upper plate. In this case, the pressure p given by Equation (16) may be independent of z so that the constant a can be set to zero. Equations (3) and (4) then where γ 1 = α 3 − α 2, γ 2 = α 6 − α 5 , the α i are constant viscosities for the NLC. These two equations form the key starting point for the two shear problems in the following problems subsections.