On the effects of Ericksen and Deborah numbers on the flow in microfluidic capillaries

ABSTRACT The use of nematics in microfluidic devices offers methods of flow control that are inaccessible to isotropic fluids. While nematics have been used in channel flows found in lab on chip devices, their flow and microstructure formation still needs to be fully understood in geometries with a cylindrical cross-section. This paper investigates liquid crystal flows through capillaries via 1) the Beris-Edwards (BE) model and 2) the Leslie-Ericksen (LE) model, which emerges as a limiting case (uniaxial, constant order parameter) of the BE model. We show that the microstructure distribution is controlled by the Ericksen ( ) and Deborah ( ) numbers in the Beris-Edwards model. Defects at low Ericksen numbers arise when the elastic energy dominates the thermotropic energy. Defect destruction provides a new shear-thinning mechanism, which cannot be observed in the Leslie-Ericksen theory. GRAPHICAL ABSTRACT


Microfluidics
Microfluidics is a rapidly emerging science focused on manipulating fluid motion in micro-scale domains [1,2].Microfluidic devices can integrate and automate several laboratory processes into a single low-footprint instrument; they are often considered as portable replacements for traditional, significantly larger chemical analysis tools [3].Micro-scale devices are advantageous in several respects [2,4]: 1) analysis only requires low volumes of process fluids (i.e.droplets) compared to traditional large-scale devices -this is particularly beneficial when the cost of processing material is high; 2) measurements can be performed almost anywhere and 3) analysis takes a much shorter time while maintaining high resolution and sensitivity.The development of miniaturised sensors is particularly important for medical applications, as low-cost and massmanufacturability are essential to detect infection and to control the spread of disease in large populations [5,6].Other potential applications of microfluidic devices include drug delivery [7], optofluidic modulators [8], colour filters [9], velocimetry devices [10] and small-scale valves [11].
Lab-on-chip devices consist of microscopic interconnected channels: T-junctions, crossed channels and manifolds [12].Their exploitation enables highly accurate flow tuning that can be used for mixing, separation and transportation purposes [12,13].Microfluidic devices often employ non-Newtonian fluids such as emulsions, foams or colloids to provide the CONTACT Kamil Fedorowicz kamil.fedorowicz@manchester.ac.uk functionality for advanced sensors [14,15].The exploitation of these devices requires a thorough understanding of the physics governing micro-scale fluid flow, where the flow dynamics are often driven by different processes to those encountered in macro-scale applications [3].Although microfluidic flows are typically laminar, complexities arise due to constitutive behaviours driven by viscoelastic, electrostatic or thermodynamic effects; these complicate the design and development of new instruments [16,17].
The functionality of many lab-on-chip devices can be further enhanced by using materials with anisotropic microstructure, as found in liquid crystals (LC) [12,18].The material properties of liquid crystals (viscosity, electrical permeability and thermal conductivity) can be controlled through surface treatment of the flow geometry or via an application of external magnetic or electric fields [19].The use of liquid crystals provides an additional degree of freedom in tuning the flow compared to devices based on isotropic fluids [18,20].For example, during mixing, the flow rate in each branch of a crossed microchannel can be controlled through the electric field, thus enabling a more accurate modulation of the chemical composition of the resultant mixture [12,21].

Liquid crystals
Liquid crystals consist of anisotropic molecules such as rods or discs, which form structures with long-range orientational order.In the case of nematic liquid crystals, the rod-like molecules lack positional order, but their long axes tend to align in the same direction (called the director n) [22].Many properties of nematic liquid crystals are strongly anisotropic; electrical permeability, thermal conductivity and optical birefringence depend on the alignment of the rod-like molecules [23].
In the absence of external fields, the director orientation in the bulk of a liquid crystal fluid is governed by the competition between flow and elastic effects, whose relative importance is quantified through the Ericksen number ðErÞ [19].Er < < 1 corresponds to an elasticitydominated limit, where the director field is insensitive to the fluid motion.The bulk director distribution is almost entirely governed by the conditions imposed on the flow boundaries, which can be controlled through surface treatment [19].Microstructure ordering becomes more intricate in complex geometries (i.e.junctions [24], contractions [25,26] or channels with obstacles [27]), where each wall may impose a different director alignment.The geometry thus drives distortions in the microstructure [16,22], and non-zero static stresses may result.By increasing Er, the region of elastically dominated microstructure shrinks.As Er > > 1, n aligns in the flow direction over most of the geometry; the director aligns in the wall-imposed direction only in a narrow near-wall layer whose size scales as Er À 1 2 [19].An opportunity for enhanced flow control exists via management of the interactions between microstructure, boundary conditions, flow and external fields [12].In [21], Na et al. exploited the anisotropic viscosity of liquid crystals in the design of microfluidic separation channels to control the flow rate in each channel.A similar approach has been adopted by Mo et al. [28], who used the same strong viscosity variations to manage the active ingredient flow rate in the design of drug delivery devices.
The interaction between flow and microstructure gives rise to elastic stresses.These have little impact on fully developed flow dynamics in straight channels and pipes.Flow and microstructure interaction do, however, affect flow development [29], secondary motion in curved pipes [30] and flow separation in more complex configurations [25].Thus, it is essential to account for effects on microstructure in order to control the fluid motion.
The motivation for early works on liquid crystals was their application in liquid crystal displays and the control of director patterns for the purpose of light modulation [23,31].Initial investigations aimed to understand the flow and microstructure evolution in a channel flow through analytical [32,33], and numerical techniques [34][35][36].Studies of liquid crystals in channels reveal many interesting phenomena (i.e.backflow [37,38] or tunable flow shaping [20,39]) that can be relevant for developing microfluidic devices.In configurations with homeotropic anchoring, Denniston et al. [36] have found that the viscosity need not be a unique function of the shear rate, because the steady-state director distribution and velocity profile depend on whether the liquid crystal is initially in a nematic or isotropic state.A similar phenomenon has been found for tumbling nematics by Manneville [40].However, little work has been conducted in capillary flows -an experimentally relevant geometry where defects can be induced through boundary conditions [41,42].Studies of discotic liquid crystals in capillaries were initiated by Lima and Rey [43], who found that the LE theory admits multiple predictions of director distribution and velocity profile for a fixed Ericksen number; shear-thinning and shear-thickening are also possible.Their work was extended for rod-like liquid crystals by Steffen et al. [44], who also found that the Leslie-Ericksen theory allows for multiple director and velocity distributions at the same Ericksen number.In this paper, we aim to extend the previous works using a higher-order tensorial theory to 1) investigate the interplay of elastic and bulk energies on the flow and 2) examine the effect of deformation on the flow and director profiles.

Paper outline
This paper extends the work conducted by Lima and Rey [43] and Steffen et al. [44] to analyse the capillary flow of liquid crystals using the Beris-Edwards model.We investigate the effect of Ericksen and Deborah numbers on the existence of defects and their influence on the order parameter.Section 2 introduces the governing equations of the Beris-Edwards and Leslie-Ericksen models and provides relevant dimensionless scales governing liquid crystal flow.Results are discussed in section 3, where we demonstrate the importance of the defect size on the director arrangement, the flow and the effective viscosity in the capillary.A summary and perspectives for future work are provided in section 4.

Methodology
The fluid motion is described by the conservation equations of mass and linear momentum.The density of liquid crystals is assumed constant, and fluid incompressibility is ensured by the continuity equation ũ is the velocity vector, whose evolution is described through the linear momentum equation where ρ is the density, D D t denotes the material derivative, p is the pressure and τ represents the total (viscoelastic) stress tensor, whose exact form is determined by the choice of the model.Tildes are used to denote dimensional quantities.
In the case of liquid crystals, the viscoelastic stress tensor arises in response to microstructure distortions induced by boundary conditions and by the flow-director interaction.The exact form of τ depends on the model.The following sections present an overview of the Beris-Edwards and Leslie-Ericksen theories.

Beris-Edwards model
The stress tensor in the Beris-Edwards model is given by [35,45]: μ is the Newtonian viscosity, � is a dimensionless parameter that controls the flow-director alignment at high Ericksen numbers, and D ¼ Ñũþð ÑũÞ T 2 denotes the symmetric part of the velocity gradient tensor.H is the molecular field, which arises due to microstructure distortions.The microstructure in Equation ( 3) is encapsulated in the order parameter tensor: where a denotes the orientation of a single nematic molecule and I is the identity tensor.ψðaÞ is a probability distribution function that quantifies the chance of finding a molecule aligned in a direction parallel to a.The system description through ψðaÞ is computationally expensive and thus is rarely used to model liquid crystal rheology [46].Instead, the mean orientation of the system can be expressed through the director field n, which represents the eigenvector corresponding to the largest eigenvalue of Q. Uniaxiality or biaxiality is determined by the relative magnitudes of the two smaller eigenvalues.Since QðaÞ ¼ QðÀ aÞ, the tensorial approach represents the director as a headless rod, which removes the directionality issue encountered in vectorial models [47].The degree of alignment between rod-like molecules is quantified through the order parameter Q ¼ 3 2 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Q : Q p .Q ¼ 1 corresponds to a perfectly aligned liquid crystal, and Q ¼ 0 describes an isotropic system [48].The variable degree of order is important when the microstructure contains defects [19].In the limit of a uniaxial, constant order parameter, the order parameter tensor can be expressed in terms of the director field The elastic stress component of Equation (3) arises due to microstructure distortions and is given by The viscoelastic stress in the Beris-Edwards model depends on flow-microstructure interaction, with the Q À tensor evolution specified through the angular momentum equation: where ðQ þ I 3 Þ □ denotes the Gordon-Schowalter derivative [16,49]: Alternatively, the Beris-Edwards transport equation can be expressed in terms of the Gordon-Schowalter derivative of Q, for which an additional term 2 D 3 appears on the RHS of eq. ( 7).The molecular field H in the Q À tensor evolution Equation ( 7) is defined as where fLdG denotes the Helmholtz free energy, and H measures the deviation of the system from the lowest free energy state.fLdG in the BE model, it has a distortional (elastic) component and a bulk-free energy term [45]: where KQ is an elastic constant and ã; b; c are phenomenological parameters of the bulk-free energy [50].In the absence of flow and distortions of Q, the system maintains a uniform order parameter that minimises the bulk-free energy [23]

Leslie-Ericksen model
The viscoelastic stress tensor in the Leslie-Ericksen model is given by [16,22] where Ñ ¼ Dn D t À n � ω is the corotational time derivative of the director and ω ¼ ÑũÀ ð ÑũÞ T 2 is the vorticity tensor.Leslie viscosities in the LE model can be expressed as a function of the Beris-Edwards material parameters [35] with the mapping equations provided in appendix 5. Apart from the Newtonian viscosity α4 , other Leslie viscosities lack a clear physical interpretation; however, their combinations can be used to describe the effective viscosity depending on the flow/director alignment [22,51].The director orientation is computed from the balance equation for angular momentum: where γ1 ¼ α3 À α2 , γ2 ¼ α3 þ α2 .h is the molecular field which quantifies the system's departure from the minimum Helmholtz free energy state.The order parameter is fixed in the LE model, so the free energy has only the distortional component, which in the one-constant approximation reads [22,52] KLE is the Frank constant in the LE model; this can be expressed in terms of the equilibrium order parameter and the Frank constant of the Beris-Edwards model.A mapping equation is provided in appendix A.

LE Model as a special case of the BE model
Material properties of liquid crystals are traditionally measured in the Leslie-Ericksen framework [53].By combining Equation ( 5) with the stress and angular momentum (Equations ( 3) and (7), respectively), the Leslie-Ericksen theory (Equations ( 12) and ( 13)) emerges as the uniaxial limiting case of the Beris-Edwards model with a constant (scalar) order parameter [35].The six Leslie viscosities of the LE theory can then be expressed in terms of four constants of the BE model (μ, �, Γ, Q eq ), as given in appendix 5. We note that, in general, it may not always be possible to match experimentally measured constants of the LE model with BE parameters.

Geometry
We consider a fully developed, axisymmetric pipe flow of liquid crystals with a pressure gradient in the axial direction; the system is depicted in Figure 1.The fully developed, axisymmetric assumption enforces @τ @z ¼ @τ @ϕ ¼ 0, and since the flow takes place in the axial direction ũr ¼ ũϕ ¼ 0. The continuity equation is thus automatically satisfied, and only the axial component of the linear momentum balance (Equation ( 2)) is nontrivial, which in cylindrical coordinates takes the form @p @z is the axial pressure gradient and τrz is the shear stress in the deformation r À z plane.We assume that the director lies in the (r À z) plane, and the axisymmetric assumption enables n to be expressed in terms of the polar angle θ visualised in Figure 1:

Boundary and initial conditions
A no-slip boundary condition is imposed on the pipe wall.The effect of director wall anchoring on the flow is investigated by implementing both homeotropic and homogeneous planar (aligned in the axial direction) boundary conditions.The director angle at the wall can be expressed in the most general case through the weak anchoring boundary condition [22,54] θ B denotes the wall-imposed director angle, and σ0 measures the strength of effects aligning n in the θ B direction.A ¼ σ0 r0 KLE quantifies the relative strength of wall anchoring with respect to elastic effects; A ¼ 0 corresponds to the zero-gradient boundary condition, and A ! 1 represents infinite anchoring type [22].For typical microfluidic flows r0 � 10 À 5 À 5 � 10 À 4 m [12,55], KLE � 4 � 10 À 11 N [20], σ0 � 10 À 7 À 10 À 3 Jm 2 [50], so using consistent non-dimensionalisation, A � 10 À 2 À 10 4 .We concentrate on the upper limit of wall anchoring strength, at which strong and infinite types of anchoring are practically indistinguishable; thus, we impose θðr 0 Þ ¼ θ B .The use of infinite anchoring in these simulations is consistent with other research groups (e.g.Yeomans and co-workers [20,34,36], Rey and co-workers [43,46,56], and Steffen et al. [44]) who have investigated the behaviour of liquid crystals in straight pipes and channels.
Anchoring in the Leslie-Ericksen theory is prescribed by setting θ B : θ B ¼ ðk þ 1  2 Þπ with integer k corresponds to homeotropic anchoring; θ B ¼ kπ corresponds to the homogeneous planar boundary condition.Due to the polar nature of the director field, we consider two configurations with homeotropic anchoring: θ B ¼ � π 2 , corresponding to a radially outward (inward) facing director for θ B positive (negative).Despite the same notional interpretation of these boundary conditions, Anderson et al. [54] have shown that the corresponding flow fields need not be identical.
Infinite anchoring is also imposed for the equivalent Beris-Edwards model, calculated based on the boundary director orientation n B : where Q B ¼ 1 is the boundary order parameter; experimental measurements actually suggest that 0:3 < Q B < 0:7, depending on the surface treatment method [57,58].For our simulations, we have imposed perfect ordering through Equation ( 19) and, while not realised in practice, this choice of ordering condition most clearly shows the effect of the wall on the rheology.
A negative pressure gradient is applied at t > 0, which induces a flow in the positive z À direction.The solution is initialised with a static fluid (ũ ¼ 0) and with the director field oriented in the axial direction (θðt ¼ 0Þ ¼ 0) in the bulk.Imposing θðt ¼ 0Þ ¼ � π 2 as an initial condition would result in arbitrarily large distortion energies at the axis in the Leslie-Ericksen theory [22].Therefore, an axially aligned initial condition in the bulk ensures consistency between the vectorbased and tensor-based approaches.Experimentally, the initial state can be achieved by applying a strong external field to a static fluid before the pressure gradient is imposed [22].

Scaling
We introduce the following scaling where r0 is the capillary radius.There is no explicit velocity scale associated with the flow, so we have chosen À @p @z r2 0 4μ , which is the mean velocity of a Newtonian fluid pipe flow with viscosity μ.The director deformation length-scale depends on the geometry size and the distortion component of the molecular field ( Hðf Q d Þ) and thus scales as is non-dimensionalised through b, as this constant is typically the largest in magnitude [23,24].
Using the scaling given by Equation ( 20), the dimensionless momentum balance in the axial direction (Equation ( 16)) becomes where the Reynolds number is defined as

Beris-Edwards model (dimensionless)
Combining the Beris-Edwards Equations ( 3) and (7) with the scaling parameters (eq.( 20)), we obtain the dimensionless form of the BE model.The stress is given by and the Q À tensor transport equation reads where the dimensionless Gordon-Schowalter derivative is given by Since no explicit velocity scale is associated with the flow, the Ericksen number in Equation ( 24) is defined in terms of the pressure gradient: A similar definition was used by Lima and Rey [43].The microstructure distribution is also controlled by the Deborah number, defined as De quantifies the relative strength of bulk and flow effects acting on the order parameter [59].Typically, De < < Er [60], so the microstructure is in the flowaligned state when De > 1, and shearing leads to the increase in the order parameter.
The contribution of flow effects to the Q À tensor transport equation is negligible when Er < < 1; De < < 1; this enables the following simplification of Equation ( 24) in the Beris-Edwards framework to

Since the Hðf
, the strength of elastic and nematic bulk effects is comparable when significant changes in the order parameter occur over a distance which can be understood as a defect size.

Leslie-Ericksen model (dimensionless)
The dimensionless stress tensor in the Leslie-Ericksen theory is given by The director orientation is found from the dimensionless angular momentum equation: The Ericksen number is defined as where the second equality was obtained with the use of equations (47) and (48).The definition of Ericksen number given by Equation (32) ensures that the same Er in the LE and BE models represents an identical partitioning of viscous and distortional effects.

Solution procedure
The transport equations of the Leslie-Ericksen model can be simplified by assuming a fully developed and axisymmetric flow/director distribution.Expressing n in terms of the polar angle and combining Equations ( 17) and ( 31) yields a single partial differential equation for the evolution of the director angle The system is described by the solution of linear and angular momentum balances (Equation ( 21) and (33), respectively); variables of the problem are the axial velocity u z and the director angle θ.Dimensionless equations are solved with the pdepe MATLAB solver designed to solve one-dimensional parabolic and elliptic partial differential equations [61].The angular momentum balance in the Beris-Edwards model cannot be simplified in a similar manner to the LE theory.Therefore, the solution of the Beris-Edwards model is obtained by solving the linear momentum Equation ( 21) for the axial velocity; the solution of the angular momentum balance (Equation ( 24)) provides rr, θθ, zz, and rz components of the Q À tensor.
All simulations are run until the flow and microstructure reach a steady state; the total system's energy remains constant because the work injected into the system through the pressure-work term is balanced by the viscous dissipation.

Material parameters
In numerical investigations, we have chosen parameters similar to Sengupta et al. [20] for 5CB.The following dimensionless combinations are used: The order parameter that minimises the bulk-free energy can be calculated from Equation (11) to give Q eq ¼ 0:6208.Material parameters of the Leslie-Ericksen theory are mapped through the equations provided in Equation (A1).Low and high Er regimes are explored.Additionally, the effect of defect size on the BE predictions is examined: � N ¼ 10 À 2 corresponds to defects much smaller than the flow geometry; � À 1 N corresponds to relatively large defect sizes.

Low Ericksen number (Er ¼ 0:1) flows with homeotropic anchoring
Pipe flows at low Er computed via the Leslie-Ericksen theory have previously been discussed by Lima and Rey [43] and Steffen et al. [44], who show how the free energy is minimised by n escaping into the axial direction (Figure 2).Although θðr ¼ 1Þ ¼ � π 2 represents the same homeotropic anchoring conditions in principal, Figure 2 shows that opposing escape directions can arise in practice on the pipe axis for the two boundary conditions.The direction of director escape depends on the relative orientations of n at the wall and at the pipe axis, such that the director distortion is minimised.Since flow effects are negligible, the Helmholtz free energy of both solutions is identical, and there is no clear indication as to which solution is the correct one.However, Figure 3 shows that the differing director solutions produce first and second normal stress differences (N 1 ¼ τ zz À τ rr , N 2 ¼ τ rr À τ θθ [16]) of opposing magnitudes.N 1 and N 2 are irrelevant in straight pipe flows, but their sign and magnitude become important when microfluidic domains contain bends and curved pipes [29,30].
In the Beris-Edwards model, the distribution of Q across the pipe strongly depends on the dimensionless defect size � N .For � N ¼ 10 À 2 , the order parameter changes significantly in the elasticity-dominated nearwall boundary layer (r > 0:95, Figure 4(a)).Outside the boundary layer, the order parameter is governed by the free energy, which drives the system towards the fixed order parameter Q eq .The Q À tensor has a pair of repeating eigenvalues (Figure 4(b)), indicating that the system is in the uniaxial state.These microstructure characteristics (constant Q and uniaxiality) constrain the predictions of the Beris-Edwards model to be very close to those of the Leslie-Ericksen theory outside of the boundary layer.Further calculations (not reported here) demonstrate linear growth of the boundary layer with � N .Therefore, in configurations where the defect size is insignificant with respect to the geometry dimension (� N !0), the liquid crystal flow can be modelled with the lower-order Leslie-Ericksen theory.
The similarity of director fields predicted by the Leslie-Ericksen theory and the Beris-Edwards model with � N ¼ 10 À 2 is reflected in the distribution of the effective viscosity η eff ð¼ j τ rz @uz @r jÞ, shown in Figure 5. Changes in η eff arise due to spatial variations of director orientation and -in the case of the Beris-Edwards model -via changes in the order parameter.For small defect sizes ð� N ¼ 10 À 2 Þ, the Beris-Edwards model and the Leslie-Ericksen theory exhibit good agreement in their predictions of effective viscosity outside of the elastic boundary layer.Although both models predict that η eff increases as the wall is approached (and where n aligns in the velocity gradient direction), the effective viscosities predicted by the BE model at the wall are ,50% higher than those of the LE model.The effect is caused by the elevated order parameter adjacent to the wall in the BE model, which increases the magnitude of implied Leslie viscosities.Here, the imposition of Q ¼ 1 at the walls represents stronger anchoring than typically encountered   (0:3 < Q < 0:7); in the latter case, less severe viscosity perturbations are expected than those reported here, which -in effect -represent a worst case estimate.As the distance from the axis decreases, the director increasingly aligns in the flow direction, reducing the effective viscosity.The similarity of viscosity profiles is reflected in the velocity distributions shown in Figure 6.The peak velocity predicted by the BE model with � N ¼ 10 À 2 is marginally lower than that predicted by the LE theory.That happens because the near-wall order parameter is larger than Q eq in the BE model, which results in a higher effective viscosity than in the LE theory.
When the defect size is comparable with the geometry size (� N � 10 À 1 ), the escaped director configuration becomes energetically unfavourable.The Helmholtz free energy is minimised in the planar-radial configuration via a reduction of the order parameter [42], as shown in Figure 7(a).Omitting the pipe axis and the capillary wall, the order parameter tensor has three distinct eigenvalues (Figure 7(b)), so the liquid crystal transforms from the uniaxial state near the wall to a biaxial state at r � 0:3, and back to the uniaxial state at the axis (Figure 8) [42].At r ¼ 0, the eigenvector corresponding to the non-repeating eigenvalue points in the axial direction, indicating a planar-uniaxial field in the capillary core.Increasing the defect size drives the microstructure towards the isotropic state, which may be an undesired effect in applications where the product functionality depends on anisotropic properties of liquid crystals, e.g.nematic valves [37].For typical values KQ ¼ 4 � 10 À 11 N, b ¼ 2:12 MJ m 3 [20], � N ¼ 10 À 1 corresponds to a capillary size of 40nm; the dimension can thus be interpreted as a limiting size for the design of microfluidic devices exploiting anisotropic effects.
An additional shear-thinning mechanism is provided by the Beris-Edwards model for the large defect case (� N ¼ 10 À 1 ), where the microstructure has a planar radial arrangement.Even though n aligns in the velocity gradient direction throughout the capillary, a reduction of the order parameter decreases the magnitude of the implied Leslie viscosities which, in turn, leads to a smaller effective viscosity (Figure 5).At all points in the capillary, the effective viscosity of the planar radial configuration is larger than that in the escaped radial configuration; the effect is reflected in the velocity profile (Figure 6), where the peak velocity for the BE model with � N ¼ 10 À 1 is ,20% smaller than the maximum velocity in the BE model with � N ¼ 10 À 2 .

Microstructure and flow arrangement
The biaxial structures observed for flows with large defects at low Ericksen numbers do not appear at intermediate Ericksen numbers.This is because the strength of shearing effects becomes comparable to distortional effects.In this latter case, the Beris-Edwards model predicts that the director reorients from homeotropic conditions at r ¼ 1 to flow aligned at the axis.The defect size has only a weak effect on the director distribution both at the wall and the pipe axis.In the near-wall region, larger values of � N increase the width of the elastic boundary layer (defined as the region of increased order parameter -Figure 9).The magnitude of the implied Frank constant increases (Equation( 48)), resulting in larger resistance to material distortions at the wall.Close to the pipe axis, flow shearing leads to the director escaping into the flow direction irrespective of the defect size (Figure 10). Figure 10 demonstrates the sensitivity of the Leslie Ericksen model to the choice of boundary and initial conditions.Contrary to the Beris-Edwards model, significant director gradients may appear near the pipe axis for particular combinations of initial and boundary conditions.While the BE results may be replicated by the choice θðr ¼ 1Þ ¼ À π 2 ; θðt ¼ 0Þ ¼ 0, Figure 10 shows that at least two other solutions are possible.One of these solutions provides the same qualitative behaviour as the Beris-Edwards solution, shifted by π; in the other solution, strong near-axis gradients in the director are apparent.Since the fluid stresses arise via director gradients, this solution leads to significantly different stress and velocity distributions.
The effective viscosities arising from the differing predictions are shown in Figure 11.η eff tends to increase with distance from the pipe axis as the director rotates towards the velocity gradient direction.Significant nearwall differences between LE and BE viscosity predictions arise from an increased order parameter there, which increases the magnitude of implied Leslie viscosities.This effect is not present in the LE model.The largest predicted viscosity disagreement occurs for the Leslie-Ericksen theory with θðr ¼ 1; tÞ ¼ π 2 , θðr ¼ 1; tÞ ¼ 0. Strong local rotation of the director near the pipe axis leads to shear-thickening via the same mechanism as in the near-wall behaviour: n locally aligns in the velocity gradient direction.
The predicted velocity distributions are given in Figure 12.For most cases, the velocity profiles are similar -changes in the near-wall effective viscosity filter through to the velocity, altering the peak predicted values and the boundary layer profiles.It is only the solution with θðr ¼ 1Þ ¼ π 2 , where a significant director rotation near the axis causes shear-thickening, and a plug-like region emerges there.
Figure 13 shows that the LE solution with significant director gradients near the axis results in a larger Helmholtz free energy, and the difference is most noticeable in the intermediate Er regime.However, there is some ambiguity as to whether this represents actual flow physics or is a consequence of the axisymmetric constraint imposed on the director at the axis coupled to the choice of boundary/initial conditions.This point is further explored in the next section.

History effects in the Beris-Edwards model
Previous results show that the Beris-Edwards model provides robust predictions with respect to the boundary and     initial conditions.Here, we use the BE model as a basis of comparison to explore whether LE solutions (with seemingly artificial director rotations) are plausible but otherwise inaccurate manifestations of the flow history, or whether they arise solely as a function of model deficiencies.To that end, we explore memory effects in the Beris-Edwards model, which may arise from the interplay between elasticity and the choice of initial conditions.
Simulations are initialised with a director arrangement arising at low Ericksen numbers, which exhibits a clockwise tilt (Figure 14a).Examining the change in the director distribution initially as a function of increasing Er, we distinguish the following structures: (1) When Er < < 1, viscous effects are weak; the director arrangement shown in Figure 14(a) is maintained.
(2) At Er � 5, viscous contributions are sufficiently strong to rotate the director only in the near-wall region (Figure 14(b)).( 3) For Er > 20, effects arising from the initial conditions are erased as the viscous contribution forces n to escape into the flow direction (Figure 14(c)).If instead, the flow is relaxed from high Er, then a new director distribution may be revealed: (a) Starting with the distribution shown in Figure 14(c) and reducing Er maintains the escape direction and weakly affects the director angle (Figure 14(d)).Thus, previous director arrangements (i.e.memory effects) are carried out by the Beris-Edwards model [24].A qualitatively similar phenomenon was observed by Denniston et al. [36], who analysed the channel flow of the Beris-Edwards model with homeotropic anchoring; their results indicate that the director orientation at the axis depends on the initial condition.A wall-normal initial condition aligns the director in the velocity gradient direction at the channel centreline; conversely, an initial isotropic state in the bulk of the channel aligns the director with the flow.The solution dependence on the initial state disappears at higher Ericksen numbers, and the configuration with the flow-aligned director appears.The effect of director orientability on the steady-state flow and microstructure distribution in a channel flow was also investigated by Anderson et al. [54].Their results show that the director orientation at the channel centreline in the Leslie-Ericksen theory is controlled solely by the relative orientation of n at the upper and lower wall: (1) the director aligns in the flow direction when walls impose opposing orientations (n upper ¼ À n lower ) and (2) the director aligns in the velocity gradient direction when both walls impose the same orientation (n upper ¼ n lower ).
In the latter case, even high shear rates cannot align the director in the flow direction at the centreline.Thus, similarly to capillary flows, the interplay between the infinite anchoring boundary condition and director orientability produces infinite memory effects in a channel flow.That need not be the case in flows with weak anchoring, where shearing effects at high Ericksen numbers can alter the director orientation at the boundary [31].

High Ericksen numbers (Er ¼ 200) with homeotropic anchoring
For high Ericksen numbers, the microstructure is governed mostly by flow effects, and the director rotates towards the Leslie angle.The polarity of the director in the LE model leads to two potential (nominally identical) Leslie angles.The angle actually adopted in the solution depends on the choice of boundary condition.Suppose the molecular field is momentarily omitted in equation (31).In that case, it is straightforward to show that the local shear/vorticity terms act in the opposite sense on a director of the form n ¼ ðcos θ; 0; sin θÞ T than they do on a director of the form n ¼ ðcosðθ þ πÞ; 0; sinðθ þ πÞÞ T .This leads to the Leslie angles θ L � À 23 o for θðr ¼ 1Þ ¼ À π 2 and θ L � 157 o for θðr ¼ 1Þ ¼ π 2 as shown in Figure 15.Despite differing director angles, director distributions in the near-wall region share similar director gradients (and, hence, similar stresses).Shear gradients reduce near the axis (r < Er À 1 3 � 0:17), and the residual elasticity forces the director to adopt the axial orientation specified through the initial conditions [19].For the case where θ L � À 23 o , the director distortion is consistent with that of the BE model.For the θ L � 157 o case, the director is forced to the axial boundary condition, resulting in a much greater director distortion at the axis.The director gradient is thus seen to be a result of the competition between infinite homeotropic anchoring and the axisymmetry introduced through the cylindrical coordinate system.Infinite anchoring need not be satisfied in real applications, and weak anchoring actually militates against the extreme distributions shown in Figure 15 [22].A derivation presented in appendix 49 shows that in configurations with weak anchoring, flow effects rotate n away from the boundary-imposed orientation and allow the solution to approach the 'principal' Leslie angle found in the other solutions.The phenomenon has been demonstrated by Cousins et al. [31], who found that weak anchoring reduces the director distortion in the near-wall region by aligning the director closer to the Leslie angle at the boundary.
When Er > > 1, the director aligns at the Leslie angle to the flow across the majority of the pipe.The resulting viscosity is essentially constant, and a Newtonian-like velocity profile is predicted by both Leslie-Ericksen and Beris-Edwards models (Figure 16).The Leslie angle is fixed in the LE theory, and depends on the relative magnitude of α 2 and α 3 [16].The Leslie angle in the Beris-Edwards model is the same as in the Leslie-Ericksen theory, provided that the order parameter takes the equilibrium value Q eq .However, shearing increases the order parameter above Q eq (Figure 17), which reduces the Leslie angle Smaller θ L leads to a decrease in the effective viscosity, which explains why the peak velocity obtained with the BE model with � N ¼ 10 À 1 takes the largest value.

Flows with homogeneous planar anchoring
The no-flow, homogeneous planar configuration does not introduce director distortions, as n aligns in the axial direction throughout the whole domain.The observed behaviour is qualitatively the same for all Ericksen numbers, so only results for Er ¼ 10 are presented below.The Leslie-Ericksen theory predicts two distinct solutions whose existence depends on the initial and boundary conditions (Figure 18).The boundary condition θðr ¼ 1Þ ¼ π is incompatible with the initial condition θðt ¼ 0; r < 1Þ ¼ 0, and produces a high director gradient near the axis (dashed line in Figure 18); the effect is not predicted in the BE model.Conversely, the solution with θðr ¼ 1Þ ¼ 0 weakly differs from the BE model prediction; the difference is largest in the near-wall region.This is because of the elevated order parameter at the wall, which increases resistance to material distortions via the mechanism discussed in section 3.2.1.In addition, the steady-state distribution of n depends only on shearing in the BE model; the history effects found in homeotropic configurations (and arising from hysteresis in the director distribution at low Ericksen numbers) do not play a role here.
The director angle is small in configurations with homogeneous planar anchoring ( À θ L < θ < θ L ).The effective viscosity is nearly constant across the domain, and the velocity profile has a parabolic shape (Figure 19).In this case, there is little difference between infinite and weak anchoring, as the director nearly aligns with the flow.The peak velocity is the largest in the Beris-Edwards model due to a reduced director angle, which decreases the effective viscosity in the same way as was discussed in section 3.3.The smallest peak velocity is predicted by the LE theory with θðr ¼ 1Þ ¼ π because the opposing boundary conditions at the axis produce a director alignment in the velocity gradient direction at r � 0:07; the effect is manifested through a flatter velocity profile near the axis (dashed line in Figure 19).

Effective viscosity
The variation of the mean effective viscosity (defined as ηðrÞrdr) as a function of Ericksen number and the defect size is illustrated in Figure 20.Elevated � η eff in the BE model with homeotropic anchoring is caused by the higher-order parameter in the near-wall region (discussed in section 3.1).Shear-thinning occurs at intermediate Ericksen numbers, and the driving mechanism depends on the director distribution when Er < < 1.For an initial configuration with a large defect at the axis, shear-thinning is initially caused by the destruction of the defect, followed by director reorientation towards the flow direction.At low Er with � N !0, only the reorientation mechanism is present in geometries where n escapes to the axial direction.In homeotropic configurations, shear-thickening occurs at intermediate Ericksen numbers for the case with boundary conditions θðr ¼ 1Þ ¼ π 2 .The effect is caused by an increased director alignment in the velocity gradient direction near the pipe axis (discussed in section 3.2.1).A shear-thickening mechanism is present in configurations with homogeneous planar anchoring.� η eff is the smallest at Er < < 1, where n aligns parallel to the flow direction; viscosity increases with Ericksen number because strong shearing imposes a small flow/ director misalignment at the Leslie angle.

Infinite Ericksen number, finite Deborah number flow
We now consider infinite Ericksen number flows, where the shearing effects are strong enough to erase the effect of wall anchoring on the director orientation.In this regime, the order parameter and director orientation are controlled only by the Deborah number.Since the elastic contribution vanishes, the flow aligns the local Leslie angle θ L [35].
Neglecting the contribution of long-range elasticity, the Q À tensor evolution equation ( 24) simplifies to Assuming that the liquid crystal is in a uniaxial state and the director is confined to the r À z plane, the Q À tensor can be expressed in terms of the director angle and scalar order parameter Combining and manipulating Equations ( 35) and ( 36), we can decompose the reduced Q À tensor into separate evolution equations for the director angle and the order parameter: A complete derivation of the above equations is presented in appendix 7. Equations ( 37) and ( 38) demonstrate explicitly the coupling between θ and Q.The steady-state order parameter is affected by the competition between bulk and viscous effects; the relative importance of these effects is measured by the term @u z @r De.The flow contribution can be neglected at small Deborah numbers, which recovers the static limit of the equilibrium order parameter [23]: Conversely, the bulk contribution can be neglected in the limit @u @r De > > 1, and eq.( 38) reaches equilibrium when Q ¼ 1; the corresponding Leslie angle is given by The analysis above shows that, even when the effect of wall anchoring (which imposes a perfect director alignment Q ¼ 1) is not important at high Ericksen numbers, the order parameter (and thus the director angle, Equation ( 34)) change spatially due to flow effects.Q is maximum at the wall (Figure 21), where the shear rate is the largest; conversely, @u z @r ¼ 0 at the axis, so the order parameter assumes the static solution is given by Equation (39).

Summary
We have compared predictions of the Beris-Edwards and Leslie-Ericksen models in a capillary flow for microfluidic applications.The Beris-Edwards model reveals new features of microstructure distribution, which cannot be captured by the LE theory: small geometries (large � N ) at low Ericksen numbers are prone to the formation of a defect at the axis.The structure is destroyed at higher Ericksen numbers ( � 3), where viscous effects rotate the director towards the axial direction.Defects are absent when � N < < 1, and the Beris-Edwards model predicts a uniaxial state.The order parameter is constant except for the elasticitydominated near-wall boundary layer, whose size scales with � N .The boundary layer disappears as � N !0, which recovers the Leslie-Ericksen limit.
The deformation history controls the microstructure distribution in the Beris-Edwards model at intermediate Ericksen numbers; n can escape along or against the flow direction, the latter of which leads to high director gradients at the axis.Memory effects in the BE model disappear at higher Ericksen numbers, where shearing contribution dominates the director orientation, and n aligns at the Leslie angle to the flow direction.In contrast, the initial condition in the Beris-Edwards model controls the steady-state solution due to the assumption of a fully-developed, axisymmetric flow.Distortional effects can be neglected when Er > > 1, and the microstructure characteristics in the BE model depend on the Deborah number, which affects both the order parameter and the Leslie angle.In the Leslie-Ericksen theory, the combination of a vector description, infinite anchoring, and initial conditions results in a locked director orientation at the axis.The feature is particularly unrealistic in the high Ericksen number regime, where despite the dominance of shear, there is a narrow reorientation layer near the axis, which increases the Helmholtz free energy of the system.A comparison with the existing literature shows that the problem can be mitigated by imposing a more realistic weak type.
Investigations of the Beris-Edwards model reveal a new mechanism of viscosity control at low Ericksen numbers: in configurations with large � N , shearthinning is caused by defect destruction.A further viscosity reduction occurs via a director reorientation, similar to the Leslie-Ericksen model.Geometry size and boundary conditions become irrelevant at higher Ericksen numbers, where n aligns at the Leslie angle in the whole domain.As a result of the uniform director orientation, the velocity profile can be approximated by the Newtonian solution.
Future work aims to investigate the flow of liquid crystals in more complex microfluidic geometries.In particular, we aim to explore the effect of deformation history on the steady-state flow and microstructure distributions in curved channels and T-junctions.

Open access statement
For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.
Let us consider the limit of strong anchoring with σ 0 > 10, so that the underlined terms in the equation above can be neglected.In the case of homeotropic anchoring (θ B ¼ � π 2 ) we can distinguish two limiting scenarios, which recover zero gradient and Dirichlet boundary conditions: • Er σ 2 0 !0, sin½4ðθ À θ B Þ� ¼ 0, which is satisfied by θ ¼ �θ B .

Figure 1 .
Figure 1.Pipe geometry considered in this paper.

Figure 2 .
Figure 2. (Colour online) Distribution of the director angle predicted by the Leslie-Ericksen and Beris-Edwards models at Er ¼ 0:1.The vector plot compares the director arrangement in the LE model when n on the boundary points towards/away from the pipe axis.

Figure 3 .
Figure 3. Distribution of the viscous component of normal stress differences predicted by the Leslie-Ericksen theory at Er ¼ 0:1.

Figure 4 .
Figure 4. (a) Order parameter distribution and; (b) eigenvalues of the Q À tensor predicted for the Beris-Edwards model with � N ¼ 10 À 2 .The microstructure is in the uniaxial state since Q has a pair of repeating eigenvalues.

Figure 6 .
Figure 6.(Colour online) Velocity profiles predicted by the Leslie Ericksen model with different boundary conditions and the Beris-Edwards model with different defect sizes, Er ¼ 0:1.

Figure 7 .
Figure 7. (a) Order parameter distribution and; (b) eigenvalues of the Q À tensor predicted for the Beris-Edwards model with � N ¼ 10 À 1 .The microstructure is in the biaxial state (except at the capillary wall and axis) since all eigenvalues of Q are different.

Figure 11 .
Figure 11.(Colour online) Viscosity profiles predicted by the Leslie Ericksen model with different boundary conditions and the Beris-Edwards model with different defect sizes, Er ¼ 10.

Figure 10 .
Figure 10.(Colour online) Distribution of the director angle for Er ¼ 10. Results of the Beris-Edwards model nearly collapse on each other, indicating a minor role of � N .Quiver plots correspond to black continuous and dashed lines, respectively.

Figure 12 .
Figure 12. (Colour online) Velocity profiles predicted by the Leslie Ericksen model with different boundary conditions and the Beris-Edwards model with different defect sizes, Er ¼ 10.

Figure 13 .
Figure 13.Helmholtz free energy per unit length of the developed flow predicted by different solutions of the Leslie-Ericksen model.

Figure 15 .
Figure 15.(Colour online) Distribution of the director angle for Er ¼ 200.

Figure 16 .
Figure 16.(Colour online) Velocity profiles predicted by the Leslie Ericksen model with different boundary conditions and the Beris-Edwards model with different defect sizes, Er ¼ 200.

Figure 18 .
Figure 18.(Colour online) The director angle predicted by the Leslie-Ericksen and Beris-Edwards models with homogeneous planar anchoring at Er ¼ 10.

Figure 19 .
Figure 19.(Colour online) Velocity profile predicted by the Leslie-Ericksen and Beris-Edwards models with homogeneous planar anchoring at Er ¼ 10.

Figure 20 .
Figure 20.(Colour online) Viscosity as a function of Ericksen number for the LE and BE models.(h) and (p) denote the homeotropic and homogeneous (planar) anchoring conditions.