Liquid crystals (LCs) display many of the flow characteristics of liquids but exhibit long range orientational order. In the nematic phase, the coupling of structure and flow leads to complex hydrodynamic effects that remain to be fully elucidated. Here, we consider the hydrodynamics of a nematic LC in a hybrid cell, where opposite walls have conflicting anchoring boundary conditions, and we employ a 3D lattice Boltzmann method to simulate the time-dependent flow patterns that can arise. Due to the symmetry breaking of the director field within the hybrid cell, we observe that at low to moderate shear rates, the volumetric flow rate under Couette and Poiseuille flows is different for opposite flow directions. At high shear rates, the director field may undergo a topological transition which leads to symmetric flows. By applying an oscillatory pressure gradient to the channel, a net volumetric flow rate is found to depend on the magnitude and frequency of the oscillation, as well as the anchoring strength. Taken together, our findings suggest several intriguing new applications for LCs in microfluidic devices.

In recent years, liquid crystals (LCs) have found an increasing range of applications that go beyond traditional display technologies to colloidal self-assembly,1,2 biosensing,3 active biomaterials,4 and state-imprinted memory effects.5 Their static properties in the bulk, such as topological defects and optical characteristics, have been studied extensively. Less is known, however, about their behavior under confinement and their hydrodynamic response to applied fields.

Advances in microfabrication have led to sophisticated technologies that include ink-jet printing and lab-on-a-chip devices. In these technologies, microfluidic devices can collect, separate, and transport a variety of heterogeneous materials, such as bacteria, colloids, or proteins.4,6–8 The liquid medium is typically provided by a simple Newtonian fluid. In this work, we examine the usage of a nematic LC as the medium, where the flow and transport properties may be considerably different from those of the traditional, simple liquids that have been considered to date.

Experimental studies of the flow of nematic LCs have primarily focused on issues related to alignment and tumbling.9–16 Jewell et al. found two steady state director fields for Poiseuille flow in a homeotropic channel.17 Earlier experiments had been conducted by Fishers and Fredrickson,18 and the two-steady state phenomenon had been proposed by Ericksen19 and simulated by Marenduzzo et al.20 and Batista et al.21 Horn and Winter and Holmes et al. discovered that a small surface pretilt can induce distinct director profiles for shear and Poiseuille flow.22–24 Sengupta et al. found that a temperature gradient across the channel can tune the flow shape of Poiseuille flow.25 A review of recent experiments on LC microfluidics can be found in Ref. 26.

Theories of nematic LC hydrodynamics date back to the 1960s, when Ericksen, Leslie, and Parodi introduced a nematodynamics model in the director representation, the so-called ELP theory.27–30 The average orientation of the local LC molecules is represented by a unit vector n, which is the central quantity of any director model. The ELP theory has been widely used to interpret the flow behavior of nematic LCs. It does, however, suffer from limitations, most notably the fact that it cannot describe topological defects. Subsequently, Beris and Edwards proposed a tensorial model based on the Poisson bracket formalism31–34 in which a symmetric traceless tensor, Q, was introduced to include information about the local average molecular orientation n and the degree of order. The advantage of the so-called Q approach is that it can describe order parameter variations in the LC. For completeness, we note that other tensorial treatments of nematodynamics are available,35,36 some of which were also derived from the Poisson bracket formalism.37,38

Within the structure of a tensor-based description, one can resort to various numerical approaches to describe LC hydrodynamics. A radial basis function method and the Stark-Lubensky formalism have been used to investigate cavity flow.39 A finite element method has been applied to solve the Qian-Sheng model.40 The lattice Boltzmann method has also been used to solve the Beris-Edwards, Qian-Sheng, or even ELP equations in the context of various LC fluid mechanics problems.20,21,41–46 A review of the lattice Boltzmann method for liquid crystal hydrodynamics can be found in Ref. 47.

Although simple (Couette and Poiseuille) and pulsatile flows of nematic LCs have been studied experimentally, numerically, and analytically,16,48–51 several key issues remain to be addressed. First, systematic studies of flow in hybrid (splay-bend) cells have not been carried out. Second, systematic numerical studies of unsteady flows are not available. In this manuscript, we develop a surface evolution equation and solve the Beris-Edwards equation, coupled with the Navier-Stokes equation, within the framework of a hybrid lattice Boltzmann/finite difference method.52,53 We find that by applying opposite flows to a hybrid cell, the resulting field is asymmetric at low Ericksen numbers. For high Ericksen numbers, the steady-state flow is symmetric, but topological transitions may occur for certain initial director orientations. The anchoring strength of the surface determines the flow-rate-difference of the opposing flows, as well as the transition Ericksen number. We further simulate an oscillatory pressure-driven flow, namely, pulsatile flow, in a hybrid cell and find that a net volumetric flow rate develops, even for very small Reynolds number. The net flow is a function of the pressure gradient, oscillation frequency, and surface anchoring.

As mentioned above, our calculations rely on solution of the Beris-Edwards and Navier Stokes equations. The hydrodynamics of the system are described by the density ρ, the velocity field u, and the pressure P. The LC phase is described in a tensorial form. An order parameter is introduced, Q=nn13I, where the unit vector n represents the molecular orientation, I is the identity tensor, and 〈〉 denotes an ensemble average. By introducing the velocity gradient Wij = ∂jui, the Beris-Edwards equation reads34 

(t+u)QR(W,Q)=ΓH,
(1)

which governs the behavior of the LC phase under flow. The part of advection term for rod-like molecules R is defined as

R(W,Q)=(ξD+Ω)(Q+I/3)+(Q+I/3)(ξDΩ)2ξ(Q+I/3)Tr(QW),

where D = (W + WT)/2 and Ω = (WWT)/2 are the symmetric and antisymmetric parts of W, respectively. Here, ξ is a material constant related to the molecular aspect ratio, and Γ is related to the rotational viscosity γ1 of the LC by Γ = 2S2/γ1,41 where S is the scalar order parameter of the nematic LC, obtained from the eigenvalues of Q. The molecular field H, i.e., the free energy driven term is given by

H=δFδQst,

where st refers to the symmetric and traceless projection, and F is the total free energy of the system, given by

F=bulkdVf+surfacedSfsurf,

where f and fsurf are the bulk and surface energy densities, respectively. Here, f takes the Landau-de Gennes form (Einstein notation assumed)54,55

f=A02(1U3)QαβQαβA0U3QαβQβγQγα+A0U4(QαβQαβ)2+12κγQαβγQαβ,
(2)

where κ is the elastic constant of the LC (here, we adopt a single-elastic-constant approximation) and is related to the Frank elastic constant K by K = 2S2κ. Parameter U controls the magnitude of S of a homogeneous static system via

S=14+34183U.

The nematic coherence length is ξN=K/A0 and determines the size of a defect in the nematic LC. It is also the natural length scale of the system, as will be seen later.

The Navier-Stokes (momentum) equation for LCs reads

ρ(t+uββ)uα=βΠαβ+ηβ[αuβ+βuα+(13ρP0)γuγδαβ].
(3)

The stress tensor Π is decomposed into a viscous stress Πv and an elastic stress Πe,56 

Παβ=Παβv+Παβe,Παβv=ρTδαβξHαγ(Qγβ+13δγβ)ξ(Qαγ+13δγβ)Hγβ+2ξ(Qαβ+13δαβ)QγϵHγϵ+QαγHγβHαγQγβ,
(4)
Παβe=(ff0)δαββQγϵδFδαQγϵ,

where f0 is the bulk free-energy density of an undistorted nematic and the temperature T is related to the speed of the sound cs by T=cs2. Isotropic pressure is P0 = ρTf + f0.

Two types of anchoring can be introduced, namely, non-degenerate or degenerate planar. Non-degenerate surfaces exhibit a preferred Q field, Qsurf = S(nsnsI/3), with a preferred surface molecular orientation ns. In the so-called Rapini-Papoular form,57 the energy penalty associated with deviations of Q from Qsurf is given by

fsurf=12W(QQsurf)2.
(5)

Degenerate planar surfaces penalize any off-plane distortion of the Q tensor; in the Fournier-Galatola form,58 they are quantified through

fsurf=W(Q̃Q̃)2,
(6)

where Q̃=Q+(S/3)I and Q̃=ZQ̃Z. Here, Z is the projection operator associated with the surface normal ν as Z = Iνν. The anchoring strength W has values that range from 10−7 to 10−3 N/m,59 depending on surface treatments. Similar to the derivation of Beris-Edwards equation (1), one can show that the governing equation of the evolution of the surface Q-field is given by

Qt=ΓsκνQ+fsurfQst,
(7)

where Γs = Γ/ξN. The above evolution equation is similar to Eq. (9) in Ref. 43. For steady flow, it reduces to the mixed boundary condition given in Ref. 21.

We use a finite-difference method to solve evolution equations (1) and (7). The lattice Boltzmann method is employed to solve Navier-Stokes equation (3) over a D3Q15 grid60 (15 discrete velocities in a 3D regular lattice). Stress implementation follows the approach proposed by Guo et al.61 A no-slip boundary condition is assumed at all surfaces. Model validation was performed by comparing our simulation results to predictions using the ELP theory.27–30 

The characteristic time scale is set by the relaxation time of the liquid crystal, τ = γ1/A0. Typically, τ ≈ 1 μs for 4-Cyano-4’-pentylbiphenyl (5CB). The characteristic length scale is set by the nematic coherence length ξN, which corresponds to the size of a LC defect. For 5CB, ξN = 6.63 nm.59 The two equations are non-dimensionalized according to these time and length scales.

As illustrated in Fig. 1, we consider a channel bounded by two walls in the z direction. To form a hybrid cell, we impose homeotropic and planar anchoring at the top and bottom walls, respectively. In the absence of flow with single-elastic-constant approximation, minimization of the elastic energy yields30 

K2θz2=0.

Thus, θ is linear in terms of the channel depth z. There are two degenerate states associated with the above solution for any ϕ0, namely, a counterclockwise state (para-field), with θ decreasing from 90° at the bottom wall to 0° at top wall, and a clockwise state (dia-field), with θ increasing from 90° at the bottom wall to 180° at the top wall. The flow is introduced by either moving the top wall in the x-direction at velocity v0 (Couette flow) or by imposing a pressure gradient ∂xP in the x-direction (Poiseuille flow). Applying the same flow type to the two degenerate states is equivalent to applying two opposite flow directions to a single state. Therefore, we focus on the counterclockwise director field (para-field) here and apply flow in the two opposite directions. An Ericksen number Er is introduced to describe the ratio of the viscous force and the elastic force

Er=ηvLzK,

where η is the viscosity and v is the characteristic velocity. When Er ≪ 1, the director field is close to that of a static system. When Er ≫ 1, viscous effects dominate and the directors align with the flow. The Reynolds number Re = ηvLz/ρ in our simulations ranges from 10−6 to 0.2. We are particularly interested in measuring the orientation of the directors described by θ and ϕ, as this information can be directly translated into optical images under cross polarized light that could then be compared to experiments. An additional, experimentally relevant parameter is the volumetric flow rate Φ, which is defined as

Φ=Avxdydzt.

Here, A is the cross section of the channel (a plane whose normal is in the x direction), and 〈〉t represents a time average (for unsteady flow). In this work, we assume that the LC is nearly incompressible (Mach number Ma ≪ 1; thus, density ρ is a constant). Therefore, the volumetric flow rate is equivalent to the mass flow rate. In the following, we occasionally refer to them simply by “flow rate.” We say that a flow is symmetric if the opposite flow directions yield identical flow fields (flow rates); otherwise, we say that the flow is asymmetric.

FIG. 1.

Channel geometry. The LC is bounded by two walls separated in z direction by Lz in an xyz Cartesian coordinate system. The two walls are parallel to the xy plane. Flow is imposed by either moving the top wall at velocity v0 or by applying a pressure gradient ∂xP along the x axis. The director field of the hybrid cell is within a plane perpendicular to the walls. The director n is projected onto the xy plane as n′. n is described by two angles θ and ϕ. θ is the polar angle between n and the z-axis. ϕ is the azimuthal angle between n′ and x axis. Before flow is turned on, the azimuthal angle of the director on the bottom wall is ϕ0.

FIG. 1.

Channel geometry. The LC is bounded by two walls separated in z direction by Lz in an xyz Cartesian coordinate system. The two walls are parallel to the xy plane. Flow is imposed by either moving the top wall at velocity v0 or by applying a pressure gradient ∂xP along the x axis. The director field of the hybrid cell is within a plane perpendicular to the walls. The director n is projected onto the xy plane as n′. n is described by two angles θ and ϕ. θ is the polar angle between n and the z-axis. ϕ is the azimuthal angle between n′ and x axis. Before flow is turned on, the azimuthal angle of the director on the bottom wall is ϕ0.

Close modal

Simulation parameters are chosen so that the system is in the flow-aligning regime. The channel height is Lz = 51ξN, the molecular constants are ξ = 0.8 and U = 3.5, implying that the equilibrium order parameter is Sbulk ≃ 0.62. The LBGK (lattice Boltzmann Bhatnagar–Gross–Krook) single relaxation time is chosen to be τf = 0.84. Surface anchoring sets an additional length scale, namely, the Kleman-de Gennes length ξKDK/W.30 To avoid uniform director fields in our hybrid cell, ξKD < Lz requires that we work with anchoring strengths of magnitude W ≥ 10−4 N/m.

We first examine Couette flow in a hybrid cell. The bottom wall is stationary and the top wall moves with velocity v0 or −v0 in the x-direction (v0 > 0 is assumed). We only consider rigid (infinite) anchoring for Couette flow. We measure the flow rates Φ+ for para-field and Φ for dia-field and plot the difference ΔΦ = |Φ| − |Φ+| against Ericksen number Er. Our results are shown in Fig. 2(a). All flow rates are scaled by Φ0, which is the flow rate of an isotropic liquid with the same viscosity in the same channel. At moderate and low Ericksen numbers, shearing the LC along opposite directions with identical shear rate yields different flow rates, i.e., ΔΦ≠0, and it is linear with respect to Er. ΔΦ is largest when the azimuthal angle is ϕ0 = 0°. Before |Φ| and |Φ+| become identical, at high Er number, |Φ+| surpasses |Φ| in the transient regime. Fig. 2(b) shows the velocity fields with different azimuthal angles ϕ0 at Er = 5. The velocity profiles of the two states are different: the dia-field has higher velocities in the lower region of the channel. The difference stems from the different responses of the director field to the flow. In Fig. 2(c), we plot the angle profile θ in terms of channel depth z at Er = 5. One observes that θ is no longer linear with respect to channel depth. When the flow is weak compared to the elastic response (Er < 1), the directors of the para-field tend to align with the flow by increasing θ, whereas for the dia-field, the directors rotate against the flow and θ decreases. For infinite anchoring conditions, we see super-linear behavior of θ for the para-field and sub-linear behavior for the dia-field. When ϕ0 = 0, the effective shear viscosity can be determined from θ via

ηeff=α1sin2θcos2θ+ηbsin2θ+ηccos2θ,

where α1 is a Leslie viscosity coefficient and ηb and ηc are the Miesowicz viscosities with ηb < ηc. The apparent shear viscosity, ηeff, is plotted in the inset of Fig. 2(b). Due to the viscous anisotropy of the LC, the shear viscosity is inhomogeneous across the channel. In a hybrid cell filled with 5CB at room temperature (25 °C), the shear viscosity can vary by a factor of 5.62 A greater θ leads to a lower viscosity, which in turn supports a larger velocity gradient. Thus, the velocity drops from v0 at the top wall to 0 at the bottom wall earlier in the para-field than in the dia-field, leading to a smaller flow rate in para-field than that in dia-field. In Figs. 2(b) and 2(c), we also calculate the velocity and angle profiles from ELP theory (solid and dashed lines). Our model is in agreement with the ELP nematodynamics theory.

FIG. 2.

Steady-state Couette flow in LC hybrid cell with infinite nondegenerate anchoring at both walls. (a) Flow rate difference versus Ericksen number between para-field (Φ+) and dia-field (Φ), normalized by the flow rate Φ0 of a Newtonian liquid with the same shear rate. Inset is the same plot at Er < 5. (b) Velocity profile at Er = 5 from simulations (symbols) and ELP theory (lines). The inset is the apparent shear viscosity against the polar angle θ. (c) Profiles of polar angle θ at Er = 5 from simulations (symbols) and ELP theory (thick lines). The thin dashed line is when flow is absent: θ = (π/2)(z/Lz); symbols and lines have the same meaning as in (b). (d) Time evolution of the directors with ϕ0 = 0° at Er = 50. The directors are colored according to the local order parameter.

FIG. 2.

Steady-state Couette flow in LC hybrid cell with infinite nondegenerate anchoring at both walls. (a) Flow rate difference versus Ericksen number between para-field (Φ+) and dia-field (Φ), normalized by the flow rate Φ0 of a Newtonian liquid with the same shear rate. Inset is the same plot at Er < 5. (b) Velocity profile at Er = 5 from simulations (symbols) and ELP theory (lines). The inset is the apparent shear viscosity against the polar angle θ. (c) Profiles of polar angle θ at Er = 5 from simulations (symbols) and ELP theory (thick lines). The thin dashed line is when flow is absent: θ = (π/2)(z/Lz); symbols and lines have the same meaning as in (b). (d) Time evolution of the directors with ϕ0 = 0° at Er = 50. The directors are colored according to the local order parameter.

Close modal

When viscous effects overcome elastic effects (Er ≥ 10), the dia-field becomes unstable and transits to para-field, which leads to identical flow rates. The time evolution of the directors is shown in Fig. 2(d). It is known that at high Er, the directors adopt the Leslie angle30 

θLeslie=±arctanα2α3,

where α2 and α3 are Leslie viscosity coefficients. For our system, θLeslie ≈ ± 76° (compared to θLeslie ≈ ± 78° for 5CB at room temperature62). At the early stage of evolution, directors near the top wall adopt the Leslie angle. They conflict with the directors near the bottom. Thus, a high elastic distortion appears and gradually approaches the bottom wall. Later, the high-elastic region becomes unstable. The order parameter drops and the directors adopt Leslie angles thereafter. At steady state, the order parameter recovers the bulk value and the channel becomes defect free. We should notice that dia-field and para-field states are topologically different. The defect has to form during the transition from a dia-field to a para-field. The critical Ericksen number Erc ≈ 50 is consistent with a similar experimental result, which reports that a uniform homeotropic liquid crystal (5CB) subjected to Poiseuille flow shows a topological transition at Erc ≈ 25.17 

We next consider Poiseuille flow in a hybrid cell. Here, we limit our analysis to ϕ0 = 0° with finite anchoring. Again we plot the flow-rate-difference ΔΦ = (|Φ| − |Φ+|)/Φ0 between para-field and dia-field states against Ericksen number Er in Fig. 3(a). ΔΦ increases to a maximum at low Er and decreases to 0 at high Er. The threshold Ericksen number, Erc, at which ΔΦ vanishes is dependent on the anchoring strength W: higher W leads to higher Erc. As shown in Fig. 3(b), velocity fields at ∂P/∂x = 105 Pa/m (Er ≈ 5) are plotted for both para- and dia-fields. Anchoring strength W has the opposite effect on the two fields. For para-field, the weakest anchoring leads to the lowest flow rate, whereas for dia-field, lower W significantly increases the velocities. The underlying reason can be understood by relying on Fig. 3(c). For the para-field, θ decreases throughout the channel due to flow alignment, and the effective shear viscosity also increases. The resulting flow rates have decreased. Dia-field director fields respond to the flow in an opposite way: the flow increases θ and the effective shear viscosity has therefore decreased, leading to increased flow rates. When W = 10−4 N/m, the relative weak anchoring exacerbates the above flow alignment effect and, consequently, the flow-rate-difference becomes the greatest. Note that the shape of the velocity fields tilts to the bottom wall in most cases, due to the fact that the shear viscosity near the bottom wall is lower than that near the top wall. The position of the velocity peak varies with anchoring strength. This implies that one can shape the velocity field by tuning surface chemistry.25 

FIG. 3.

Steady state of Poiseuille flow in LC hybrid cell. (a) Flow rate difference versus Ericksen number for various anchoring strengths. Data points are connected by lines to guide the eyes. Velocity profiles (b) and polar angle profiles (c) with ∂P/∂x = 105 Pa/m (Er ≈ 5). Solid lines are predicted by ELP theory. (b) and (c) share the same legend. (d) Stability diagram: “the blue circle”: both dia and para-fields are stable; “the green square”: para-field is stable; and “the red diamond”: dia-field is stable.

FIG. 3.

Steady state of Poiseuille flow in LC hybrid cell. (a) Flow rate difference versus Ericksen number for various anchoring strengths. Data points are connected by lines to guide the eyes. Velocity profiles (b) and polar angle profiles (c) with ∂P/∂x = 105 Pa/m (Er ≈ 5). Solid lines are predicted by ELP theory. (b) and (c) share the same legend. (d) Stability diagram: “the blue circle”: both dia and para-fields are stable; “the green square”: para-field is stable; and “the red diamond”: dia-field is stable.

Close modal

We next address the origin of the peak in the flow-rate-difference plot of Fig. 3(a). When Er ≤ 1, the flow rate increases as the pressure gradient increases and, therefore, the flow-rate-difference between opposite flowing directions depends monotonically on the Ericksen number. When Er approaches the transition value Erc, the relation between Er and Φ± is no longer linear. As the directors approach the Leslie angle, they have less space to rotate, and the difference between the two fields decreases. Recall that, for Couette flow, the dia-field is stable only at low Er and, for para-field, it is stable for the entire range of Er’s. In Fig. 3(d), we show the stability diagram for Poiseuille flow. The dia-field is unstable in a narrow window of Er numbers. The transition between both fields being stable and para-field being stable depends on the anchoring strength W. This transition is what we observe from Fig. 3(a) when the flow-rate-difference vanishes. However, the transition between the para-field being stable and the dia-field being stable at a higher Er is independent of W. In Fig. 4(a), we show how the directors of a para-field transit to a dia-field at Er ≈ 500, where transient defects emerge and disappear in the channel center. The upper and lower regions of the channel adopt Leslie angles of opposite sign. The conflict happens at the lowest velocity gradient region (near the channel center), which becomes unstable and eventually rotates the directors from vertical to horizontal. In Fig. 4(b), we show the time evolution of the directors of a dia-field at Er ≈ 250. For this smaller velocity gradient, the directors in the lower region of the channel also adopt the Leslie angle. Frustration happens near the channel center, where the order parameter drops from its bulk value of 0.62 to 0.5 and stabilizes at this specific value. For the dia-field state at even lower Er, the lower region of the channel does not adopt a Leslie angle and remains at an opposite direction with respect to the velocity gradient. The anchoring of the walls controls the angle at steady state and therefore controls at what shear rate the directors conflict with the surface preferred angle. Therefore, the first transition, between two stable fields and one stable para-field, is related to the anchoring strength. However, for the second transition, the stability of the frustration region is only dependent on the local velocity shear and thus independent of the anchoring strengths of the walls.

FIG. 4.

Time evolution of the director fields for a Poiseuille flow with W = 10−3 N/m at Er ≈ 500 (a) and Er ≈ 250 (b). Directors are colored by order parameter.

FIG. 4.

Time evolution of the director fields for a Poiseuille flow with W = 10−3 N/m at Er ≈ 500 (a) and Er ≈ 250 (b). Directors are colored by order parameter.

Close modal

We have examined the temporal velocity field and the pressure profile before the flow becomes steady for both Couette and Poiseuille flows. When the system is free of defects, the 2D velocity field gradually approaches steady-state configuration. Whereas when the system undergoes a topological transition, the velocity near the transient defect shows a sudden jump, followed by a decaying oscillation about the steady solution. Provided that our calculations are focused on a creeping flow condition, we do not see any evidence of cross flow or other complex flow patterns in our results. The associated pressure profile behaves correspondingly. The low hydrostatic pressure region arises near the domain where the LC exhibits a large distortion, which leads to high local elastic stress.

Building on the idea that flow in a hybrid cell is asymmetric, we impose an oscillatory pressure gradient ∂xP = Gsin(ωt) to the channel, where G is the amplitude and ω = 2π/T is the frequency and T is the period. It is known that a Newtonian fluid at low Reynolds number (Stokes flow) shows no net flow subjected to slow oscillatory pressure drop. For a viscoelastic liquid, this is not true even at low Reynolds number. The Womersley number, α, is defined as α=Lz/2ν/ω.63 The kinematic viscosity ν = γ1/ρ. Here, α ranges from 0.12–6.0. In Fig. 5, we plot the net flow rate |Φ| against the oscillation period T. For small pressure gradient, |Φ| is inversely dependent on α. Rapid oscillations lead to a low net flow, and slow oscillations yield a high net flow. This observation can be understood from the velocity and director profiles in Fig. 6. When T is comparable to the LC relaxation time τ, the director has little time to respond to the pressure gradient and thus remains linear; the associated velocities are small as well. As the frequency decreases, or when T > 10τ, one starts to see the asymmetric oscillations of the director, as well as the velocity fields, which contributes to a greater net flow. In the extreme case, when Tτ, the outlines of the velocity and director angle profiles are essentially the steady-state solutions of the dia and para-fields. The magnitude of the extreme velocity reaches a maximum when T ≈ 10τ. This is because when the oscillation period equals the viscous time (≈5τ), the LC directors resonate with the pressure oscillations. At the limit of high frequency (T is small or α is large), weaker anchoring leads to higher net flow. This is consistent with what we have found in steady Poiseuille flow, where weaker surface anchoring allows more elastic deformation, resulting in higher asymmetry of the flow fields. For a moderate or high pressure gradient, the net flow only exists at high frequencies. As shown in Figs. 5(b) and 5(c), the transition period at which |Φ| vanishes is dependent on the pressure gradient, but the surface anchoring has little influence. Fig. 7 shows the details of the underlying mechanism. For rapid pressure oscillations, the director has no time to perform a topological transition (as in Fig. 4) and |Φ| remains nonzero. However, when the pressure oscillation slows down, a topological transition occurs and the director of the bulk becomes horizontal. The transition of the directors being vertical to horizontal drastically lowers the effective shear viscosity, and the extreme velocity increases accordingly. As a consequence of the topological transition, only the director near the top wall oscillates according to the pressure, the responding flow becomes symmetric, and the net flow vanishes. As found in the steady flow, the higher the pressure gradient, the shorter time the LC needs to undergo a topological transition. Therefore, for smaller oscillation magnitudes, the transition period is longer than that observed at higher oscillation magnitudes.

FIG. 5.

Net flow rate |Φ| versus oscillation period T with G = − 3 × 106 Pa/m (a), G = − 3 × 107 Pa/m (b), and G = − 3 × 108 Pa/m (c). They share the same legend as shown in (c).

FIG. 5.

Net flow rate |Φ| versus oscillation period T with G = − 3 × 106 Pa/m (a), G = − 3 × 107 Pa/m (b), and G = − 3 × 108 Pa/m (c). They share the same legend as shown in (c).

Close modal
FIG. 6.

Velocity fields (top rows) and director fields (bottom rows) for a Pulsatile flow at W = 10−3 N/m and G = − 3 × 105 Pa/m. Oscillation period is T = 2τ (a), 10τ (b), 250τ (c), and 1000τ (d). Different lines are from different times. Dashed lines in (d) correspond to vx = ± 0.5v0.

FIG. 6.

Velocity fields (top rows) and director fields (bottom rows) for a Pulsatile flow at W = 10−3 N/m and G = − 3 × 105 Pa/m. Oscillation period is T = 2τ (a), 10τ (b), 250τ (c), and 1000τ (d). Different lines are from different times. Dashed lines in (d) correspond to vx = ± 0.5v0.

Close modal
FIG. 7.

Velocity fields (top rows) and director fields (bottom rows) for a Pulsatile flow at W = 10−3 N/m and G = − 3 × 107 Pa/m. Oscillation period is T = 0.4τ (a), 2τ (b), 10τ (c), and 250τ (d). Different lines are from different times.

FIG. 7.

Velocity fields (top rows) and director fields (bottom rows) for a Pulsatile flow at W = 10−3 N/m and G = − 3 × 107 Pa/m. Oscillation period is T = 0.4τ (a), 2τ (b), 10τ (c), and 250τ (d). Different lines are from different times.

Close modal

In this work, we have developed a surface evolution equation and coupled it with a lattice Boltzmann solver to examine the flow of liquid crystals described by the Beris-Edwards equation. Finite anchoring boundary conditions and time dependent dynamics can be handled by our method. We observe that at low and moderate Ericksen numbers (Er ≤ 1), the Couette and Poiseuille flows in a hybrid cell are asymmetric and dependent on the direction of the shear or pressure gradient. The difference in flow rates is linear in terms of the Ericksen number. At high Ericksen numbers (Er ≥ 10), the director fields may transit from one state to the other through a topological transition during which a transient defect emerges. The resulting flow becomes symmetric at steady state. The transition Ericksen number is dependent on the surface anchoring strength. Our low Ericksen calculations are in agreement with predictions of the ELP theory, which cannot describe topological transitions. Our high Ericksen number observations are consistent with an experimental result in which a critical Ericksen number of the same order of magnitude was reported. We then extend the method to simulate pulsatile flow in a LC channel. A sinusoidal pressure gradient is applied to the hybrid cell. We show that the net flow depends on the oscillation frequency and the anchoring strength, as well as the pressure magnitude. At high pressures and low frequency, the LC undergoes a topological transition, the resulting flow field becomes symmetric, and the net flow rate vanishes.

Our surface evolution equation enables us to systematically study the temporal behavior of LC hydrodynamics. The asymmetry of the hybrid cell due to flow stems from a spontaneous symmetry breaking of the splay-bend structure. For more complicated director configurations, one should expect similar behaviors. Unlike Newtonian fluids, which have no such internal structure, viscoelastic liquids could find useful applications in future microfluidic devices.

The studies of transient, non-equilibrium behavior of flowing liquid crystals presented here were supported by the National Science Foundation, Grant No. DMR-1410674. The authors are grateful for the support from the University of Chicago and Argonne National Laboratory for development of innovative collaborative efforts. The development of robust algorithms and computational codes for simulation of nanostructured complex fluids was supported by the Department of Energy, Basic Energy Sciences, Division of Materials Research, through the Midwest Integrated Center for Computational Materials (MICCoM). The research of I.S.A. was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering.

1.
P.
Poulin
,
H.
Stark
,
T. C.
Lubensky
, and
D. A.
Weitz
,
Science
275
,
1770
(
1997
).
2.
I.
Muševič
,
M.
Škarabot
,
U.
Tkalec
,
M.
Ravnik
, and
S.
Žumer
,
Science
313
,
954
(
2006
).
3.
I.-H.
Lin
,
D. S.
Miller
,
P. J.
Bertics
,
C. J.
Murphy
,
J. J.
de Pablo
, and
N. L.
Abbott
,
Science
332
,
1297
(
2011
).
4.
S.
Zhou
,
A.
Sokolov
,
O. D.
Lavrentovich
, and
I. S.
Aranson
,
Proc. Natl. Acad. Sci. U. S. A.
111
,
1265
(
2014
).
5.
I.
Muševič
and
S.
Žumer
,
Nat. Mater.
10
,
266
(
2011
).
6.
L. R.
Huang
,
E. C.
Cox
,
R. H.
Austin
, and
J. C.
Sturm
,
Science
304
,
987
(
2004
).
7.
D. L.
Huber
,
R. P.
Manginell
,
M. A.
Samara
,
B.-I.
Kim
, and
B. C.
Bunker
,
Science
301
,
352
(
2003
).
8.
T.
Bowman
,
J.
Frechette
, and
G.
Drazer
,
Lab Chip
12
,
2903
(
2012
).
9.
P. T.
Mather
,
D. S.
Pearson
, and
R. G.
Larson
,
Liq. Cryst.
20
,
527
(
1996
).
10.
P. T.
Mather
,
D. S.
Pearson
, and
R. G.
Larson
,
Liq. Cryst.
20
,
539
(
1996
).
11.
A. E.
White
,
P. E.
Cladis
, and
S.
Torza
,
Mol. Cryst. Liq. Cryst.
43
,
13
(
1977
).
12.
K.
Hongladarom
and
W. R.
Burghardt
,
Macromolecules
26
,
785
(
1993
).
13.
I.
Soga
,
A.
Dhinojwala
, and
S.
Granick
,
Jpn. J. Appl. Phys.
38
,
6118
(
1999
).
14.
A. V.
Zakharov
and
R. Y.
Dong
,
Phys. Rev. E
65
,
052701
(
2002
).
15.
A. V.
Zakharov
and
R. Y.
Dong
,
J. Chem. Phys.
116
,
6348
(
2002
).
16.
S.
Pasechnik
,
I.
Nasibullayev
,
D.
Shmeliova
,
V.
Tsvetkov
,
L.
Zhijian
, and
V.
Chigrinov
,
Liq. Cryst.
33
,
1153
(
2006
).
17.
S. A.
Jewell
,
S. L.
Cornford
,
F.
Yang
,
P. S.
Cann
, and
J. R.
Sambles
,
Phys. Rev. E
80
,
041706
(
2009
).
18.
J.
Fishers
and
A. G.
Fredrickson
,
Mol. Cryst.
8
,
267
(
1969
).
19.
J. L.
Ericksen
,
Trans. Soc. Rheol.
13
,
9
(
1969
).
20.
D.
Marenduzzo
,
E.
Orlandini
, and
J. M.
Yeomans
,
J. Chem. Phys.
121
,
582
(
2004
).
21.
V. M. O.
Batista
,
M. L.
Blow
, and
M. M.
Telo da Gama
,
Soft Matter
11
,
4674
(
2015
).
22.
B. L. V.
Horn
and
H. H.
Winter
,
Rheol. Acta
39
,
294
(
2000
).
23.
B. L. V.
Horn
,
D. M.
Boudreau
, and
H. H.
Winter
,
Rheol. Acta
42
,
585
(
2003
).
24.
C. J.
Holmes
,
S. L.
Cornford
, and
J. R.
Sambles
,
Phys. Rev. Lett.
104
,
248301
(
2010
).
25.
A.
Sengupta
,
U.
Tkalec
,
M.
Ravnik
,
J. M.
Yeomans
,
C.
Bahr
, and
S.
Herminghaus
,
Phys. Rev. Lett.
110
,
048303
(
2013
).
26.
A.
Sengupta
,
S.
Herminghaus
, and
C.
Bahr
,
Liq. Cryst. Rev.
2
,
73
(
2014
).
27.
J. L.
Ericksen
,
Mol. Cryst. Liq. Cryst.
7
,
153
(
1969
).
28.
F. M.
Leslie
,
Q. J. Mech. Appl. Math.
19
,
357
(
1966
).
30.
M.
Kleman
and
O. D.
Lavrentovich
,
Soft Matter Physics
(
Springer
,
2001
).
31.
B. J.
Edwards
and
A. N.
Beris
,
J. Non-Newtonian Fluid Mech.
35
,
51
(
1990
).
32.
B. J.
Edwards
,
A. N.
Beris
,
M.
Grmela
, and
R. G.
Larson
,
J. Non-Newtonian Fluid Mech.
36
,
243
(
1990
).
33.
A. N.
Beris
and
B. J.
Edwards
,
J. Rheol.
34
,
55
(
1990
).
34.
A. N.
Beris
and
B. J.
Edwards
,
Thermodynamics of Flowing Systems with Internal Microstructure
(
Oxford University Press, Inc.
,
New York, Oxford
,
1994
).
35.
N.
Kuzuu
and
M.
Doi
,
J. Phys. Soc. Jpn.
52
,
3486
(
1983
).
36.
P. D.
Olmsted
and
P.
Goldbart
,
Phys. Rev. A
41
,
4578
(
1990
).
37.
T.
Qian
and
P.
Sheng
,
Phys. Rev. E
58
,
7475
(
1998
).
38.
H.
Stark
and
T. C.
Lubensky
,
Phys. Rev. E
67
,
061709
(
2003
).
39.
J. P.
Hernández-Ortiz
,
B. T.
Gettelfinger
,
J.
Moreno-Razo
, and
J. J.
de Pablo
,
J. Chem. Phys.
134
,
134905
(
2011
).
40.
R.
James
,
E.
Willman
, and
F. A.
Fernández
,
IEEE Trans. Electron Devices
53
,
1575
(
2006
).
41.
C.
Denniston
,
E.
Orlandini
, and
J. M.
Yeomans
,
Phys. Rev. E
63
,
056702
(
2001
).
42.
C.
Denniston
,
D.
Marenduzzo
,
E.
Orlandini
, and
J. M.
Yeomans
,
Philos. Trans. R. Soc., A
362
,
1745
(
2004
).
43.
T. J.
Spencer
and
C. M.
Care
,
Phys. Rev. E
74
,
061708
(
2006
).
44.
D.
Marenduzzo
,
E.
Orlandini
, and
J. M.
Yeomans
,
Europhys. Lett.
64
,
406
(
2003
).
45.
C.
Denniston
,
E.
Orlandini
, and
J. M.
Yeomans
,
Comput. Theor. Polym. Sci.
11
,
389
(
2001
).
46.
N.
Sulaiman
,
D.
Marenduzzo
, and
J. M.
Yeomans
,
Phys. Rev. E
74
,
041708
(
2006
).
47.
C. M.
Care
and
D. J.
Cleaver
,
Rep. Prog. Phys.
68
,
2665
(
2005
).
48.
M. G.
Clark
,
F. C.
Saunders
,
I. A.
Shanks
, and
F. M.
Leslie
,
Mol. Cryst. Liq. Cryst.
70
,
195
(
1981
).
49.
P. T.
Mather
,
D. S.
Pearson
, and
W. R.
Burghardt
,
J. Rheol.
39
,
627
(
1995
).
50.
P.
Tóth
,
A. P.
Krekhov
,
L.
Kramer
, and
J.
Peinke
,
Europhys. Lett.
51
,
48
(
2000
).
51.
L. R. P.
de Andrade Lima
and
A. D.
Rey
,
Liq. Cryst.
33
,
711
(
2006
).
52.
J. S.
Lintuvuori
,
D.
Marenduzzo
,
K.
Stratford
, and
M. E.
Cates
,
J. Mater. Chem.
20
,
10547
(
2010
).
53.
M. E.
Cates
,
O.
Henrich
,
D.
Marenduzzo
, and
K.
Stratford
,
Soft Matter
5
,
3791
(
2009
).
54.
P.
de Gennes
and
J.
Prost
,
The Physics of Liquid Crystals
(
Oxford University Press, Inc.
,
Oxford
,
1995
).
55.
L.
Landau
and
E.
Lifshitz
,
Statistical Physics
, 3rd ed. (
Pergamon-Press
,
Oxford
,
1980
).
56.
J. i.
Fukuda
,
H.
Yokoyama
,
M.
Yoneya
, and
H.
Stark
,
Mol. Cryst. Liq. Cryst.
435
,
63/[723]
(
2005
).
57.
M.
Papoular
and
A.
Rapini
,
Solid State Commun.
7
,
1639
(
1969
).
58.
J.
Fournier
and
P.
Galatola
,
Europhys. Lett.
72
,
403
(
2005
).
59.
M.
Ravnik
and
S.
Žumer
,
Liq. Cryst.
36
,
1201
(
2009
).
60.
Z.
Guo
and
C.
Shu
,
Lattice Boltzmann Method and Its Applications in Engineering
, 1st ed. (
World Scientific Publishing Company
,
Singapore
,
2013
).
61.
Z.
Guo
,
C.
Zheng
, and
B.
Shi
,
Phys. Rev. E
65
,
046308
(
2002
).
62.
I. W.
Stewart
,
The Static and Dynamic Continuum Theory of Liquid Crystals
(
Taylor and Francis
,
London, New York
,
2004
).