Flow and nematic director profiles in a microfluidic channel: the interplay of nematic material constants and backflow

Flow and nematic director profiles in a microfluidic channel: the interplay of nematic material constants and backflow Sourav Mondala, Ian M. Griffithsa, Florian Charletb and Apala Majumdarc,∗ a Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK b ENSTA ParisTech, 828, Boulevard des Maréchaux, 91120 Palaiseau, France c Department of Mathematical Sciences, University of Bath, Bath BA2 7AY, UK * Correspondence: a.majumdar@bath.ac.uk; Tel.: +44(0)1225 385321, +44(0)1865 615137 Academic Editor: name Version June 1, 2018 submitted to Fluids; Typeset by LATEX using class file mdpi.cls Abstract: We numerically and analytically study the flow and nematic order parameter profiles in 1 a microfluidic channel, within the Beris–Edwards theory for nematodynamics, with two different 2 types of boundary conditions – strong anchoring/Dirichlet conditions and mixed boundary 3 conditions for the nematic order parameter. We primarily study the effects of the pressure gradient, 4 the effects of the material constants and viscosities modelled by a parameter L2 and the nematic 5 elastic constant L∗, along with the effects of the choice of the boundary condition. We study 6 continuous and discontinuous solution profiles for the nematic director and these discontinuous 7 solutions have a domain wall structure, with a layered structure that offers new possibilities. 8 Our main results concern the onset of flow reversal as a function of L∗ and L2, including the 9 identification of certain parameter regimes with zero net flow rate. These results are of value in 10 tuning microfluidic geometries, boundary conditions and choosing liquid crystalline materials for 11 desired flow properties. 12


Introduction
Nematic liquid crystals are classical examples of partially ordered complex liquids for which the constituent molecules have translational freedom but exhibit a degree of long-range orientational ordering or certain preferred directions of averaged molecular alignment, that vary in space and time [1].The nematic hydrodynamics is particularly rich because of the intrinsic coupling between fluid motion and nematic molecular orientations i.e., the fluid motion influences the nematic orientational ordering and equally, the inhomogeneities in the orientational ordering have a kick-back effect on the fluid flow, a phenomenon known as "backflow" [2].Backflow has no counterpart in conventional isotropic Newtonian fluids.Consequently, nematics can offer unusual mechanical and rheological properties compared to their Newtonian counterparts, such as complex wetting transitions, surface effects and stable topological defects.Backflow is of fundamental scientific interest and equally, has practical consequences on switching rates of liquid crystal display devices and their refresh times [3,4].
There are two popular hydrodynamic theories for nematic liquid crystals in the literature-the Ericksen-Leslie theory [5,6] and the Beris-Edwards model [7].In the Ericksen-Leslie framework, we have two variables-the flow field and the nematic director, which is interpreted as the direction of preferred averaged molecular alignment in space.The typical mathematical framework comprises the incompressibility constraint and evolution equations for the flow field and the nematic director.The evolution equations for the flow field and the nematic director are intrinsically coupled with new anisotropic stresses, compared to the isotropic Newtonian counterpart, and the solution landscapes depend on flow parameters (such as the pressure gradient) and nematic parameters (nematic material constants, temperature, boundary conditions for the nematic director and nematic viscosities) [8].The Ericksen-Leslie theory for nematodynamics is based on the premise that the nematic is purely uniaxial, with one single distinguished direction of orientational ordering, referred to as "director", with a constant degree of orientational order.As said before, the director is interpreted as the single preferred direction of molecular alignment in the sense that all directions perpendicular to the director are physically equivalent.Hence, the Ericksen-Leslie theory is hugely successful for modelling situations that are expected to have a uniform degree of nematic ordering; this usually holds for defect-free situations or for certain choices of material constants that promote perfect nematic ordering such as the vanishing elastic constant limit of the Landau-de Gennes theory studied in [9].However, the Ericksen-Leslie theory cannot capture sharp variations in the degree of nematic ordering, complicated topological defects and biaxiality, for which there is a primary and secondary direction of preferred molecular alignment, since the Ericksen-Leslie theory only has two dependent variables.The Beris-Edwards theory is more general than the Ericksen-Leslie since it employs the Landau-de Gennes Q-tensor order parameter to describe the nematic orientational ordering.The Landau-de Gennes Q-tensor order parameter is a symmetric, traceless 3 × 3 matrix that contains information about the preferred directions of nematic molecular alignment and the degree of ordering about these directions within its eigenvectors and eigenvalues respectively [9][10][11].The Landau-de Gennes Q-tensor can capture both uniaxial and biaxial states, along with variable orientational order and is hence better suited to capture finer structural details and topological defects.The evolution equations for the flow field and the Q-tensor are again coupled through "coupling stresses".A detailed discussion of the Beris-Edwards model and its connections to closely related models can be found in the literature [12][13][14].
We work in a reduced Beris-Edwards framework to model a microfluidic channel, with an applied pressure gradient to induce fluid flow, and different types of boundary conditions for a reduced Q-tensor on the channel walls with the usual no-slip boundary conditions for the flow field.We use a reduced Q-tensor, which only has two degrees of freedom, in contrast to the usual five degrees of freedom in a three-dimensional approach.This reduced approach has been successfully used for severely confined systems elsewhere [15,16] and can be related to the usual Landau-de Gennes Q-tensor explicitly [17].In particular, we model the microfluidic channel as a two-dimensional domain and this reduced approach is successful in capturing the in-plane system characteristics.The two degrees of freedom of the reduced Q-tensor are an angle θ that describes the preferred in-plane alignment of the nematic molecules or the direction of the nematic director n in the plane, and a scalar order parameter s, that is a measure of the degree of orientational order about the director n.We note that the Ericksen-Leslie framework does not include the order parameter s.The Beris-Edwards system can be recast as a coupled system for s, θ and the flow field parameterised by u (since we assume unidirectional channel flow).We study the effects of certain key variables on the long-time or equilibrium profiles for s, θ and u.Namely, we look at the effects of the pressure gradient p x , the nematic elastic constant, the nematic material constants and viscosities (modelled by L 2 ) and the anchoring conditions for θ (modelled by either the winding number ω or the surface anchoring coefficient B).For small L 2 , the evolution equation for the flow field effectively reduces to the Navier-Stokes equation and we recover the usual Poiseuille flow.However, the flow field does influence the θ profiles in this regime and we carry out some explicit analysis to compute the s and θ profiles in this limit, for both small and large L * .The analysis captures both continuous and discontinuous solution profiles for θ; the discontinuous profiles are featured by domain walls with layered structures such that θ jumps abruptly across an interface.Again, the discontinuous solutions cannot be captured by the Ericksen-Leslie approach.The discontinuous solutions are the analogue of the well studied "order-reconstruction" solutions [18], with the novel feature of flow effects.For small p x , the flow is negligible as expected.The most interesting regime is when p x and L 2 are of comparable magnitude and there is two-way coupling between the flow and nematic order, where the effects of backflow are most pronounced.We expect the asymptotic analysis in this paper to be useful for subsequent detailed analysis of the Beris-Edwards system in this interesting regime.
There is a large body of existing literature on the Ericksen-Leslie theory and the Beris-Edwards theory.For example, an extensive review of existence and regularity results in the Ericksen-Leslie framework can be found in [19].In [20], the authors use perturbation methods in the Ericksen-Leslie framework to study the effects of backflow on defect dynamics.In parallel, there are several papers that focus on the role of backflow in the hydrodynamics of defects, in the Beris-Edwards framework, see for example [12,21,22].In recent years, there are rigorous existence and regularity results for the Beris-Edwards framework too [11,14,23] and numerical simulations for microfluidic set-ups in [24,25].The various dynamical theories of nematic liquid crystals and the key results are surveyed in [26] and in [27]; the authors rigorously derive the Ericksen-Leslie equations from the Beris-Edwards model.In [28], the authors use a lattice-Boltzmann algorithm to study nematodynamics in a microfluidic channel, in the Beris-Edwards framework, with both Dirichlet and mixed boundary conditions on the channel walls.The emphasis is on the flow rate as a function of the applied pressure gradient and the qualitative effects of the boundary conditions on the director profiles.Our setting is similar but not the same as in [28].For example, our Dirichlet conditions are inhomogeneous i.e., different fixed boundary conditions on the two bounding surfaces, whereas the authors employ the same Dirichlet condition on both surfaces in [28].A large part of the elegant asymptotics in [28] is carried out in the L * → 0 limit, for which we expect a uniform degree of nematic order or s ≈ 1 almost everywhere.This limit cannot capture the discontinuous solutions for θ described above.Importantly, our emphasis is on "flow reversal" as a function of the pressure gradient, the material and temperature-dependent parameter L 2 and the re-scaled elastic constant L * and this is not addressed in [28].In fact, flow reversal or flow in the direction of increasing pressure gradient is a distinct manifestation of backflow, only observable for L 2 large enough; this warrants further study in the future.
Our main findings can be summarised as follows.
(a) We compute a phase plane in terms of L * and L 2 , for a fixed p x , which demarcates regions of fluid flow in the direction of decreasing pressure from regions of fluid flow in the direction of increasing pressure and this flow reversal is a clear manifestation of backflow.(b) We compute the total flow rates in different parameter regimes.In particular, we show that backflow can be attained for a window of values of L * i.e., L * crit,1 < L * < L * crit,2 and these critical values depend on p x , L 2 , the boundary conditions and other material parameters.(c) We study two different kinds of boundary conditions for θ-Dirichlet and mixed boundary conditions.The mixed boundary conditions are phrased in terms of an anchoring coefficient B on the bottom surface and accompanied by a Dirichlet condition on the top surface.The mixed boundary conditions offer greater scope for tuning the solution landscape.(d) We perform some investigations on how we can choose a suitable initial condition to attain the discontinuous solution for θ at long times, and this may be useful for studying multistability in such model settings.
The paper is organised as follows.In Section 2, we present the governing equations, the boundary conditions and the initial conditions.In Section 3, we present our results on the effects of p x , L 2 , ω and L * .In Section 4, we perform the explicit analysis in the small L 2 limit and in Section 5, we give some conclusions.

Theory
We study spatio-temporal pattern formation in a long microfluidic channel where L/H 1, filled with nematic liquid crystals under the action of a pressure gradient applied at the end x = −H in the x direction.This pressure gradient naturally induces a fluid flow and we assume a unidirectional channel flow in the x-direction.There are two main macroscopic variables of interest: the flow field u = (u (x, y, t) , 0, 0), where t is time, and the nematic order parameter, which is a measure of the nematic ordering.We work in a reduced two-dimensional Landau-de Gennes framework, similar to the setting in [15,16] for which the nematic order parameter q is a symmetric, traceless 2 × 2 matrix, with two degrees of freedom.Equivalently, we can write q = s (n n − I/2) where I is the 2D identity matrix, with s = √ 2 |q| being the scalar order parameter and n = (cos θ, sin θ) being the two-dimensional director.The general Landau-de Gennes nematic order parameter, Q, is a symmetric, traceless 3 × 3 matrix with five degrees of freedom but for severely confined systems, where the vertical z-dimension is much smaller than the lateral dimensions; it is reasonable to assume that structural details are independent of the z-coordinate and we can relate the reduced tensor, q in this paper to the full Landau-de Gennes Q tensor as has been done in [17].A reduced approach, such as the one employed in this paper and others, is analytically and computationally more efficient and is a physically relevant approach for severely confined systems.
We work within the standard and powerful Beris-Edwards theoretical framework for nematodynamics [24].There are three governing equations: the incompressibility constraint, an evolution equation for the flow field, which is essentially the Navier-Stokes equation with an additional stress (σ) due to the nematic ordering and an evolution equation for q, which has an additional stress induced by the fluid vorticity.The governing partial differential equations are given below.
where D Dt = ∂ ∂t + u • ∇ is the material derivative, ρ and µ are the density and viscosity of the fluid medium, p is the hydrodynamic pressure, u is the fluid velocity, and prime denotes the transpose.The nematic stress (σ) is [9,24,29] where s = √ 2|q| is the scalar order parameter, h is the molecular field, κ is the nematic elastic constant and a and c are parameters related to the temperature and material properties.We work with temperatures below the nematic-isotropic transition temperature and hence we take a > 0. The evolution equation for the q tensor is given by [11,29] where γ is the rotational diffusion constant [10,21] and ξ is the anti-symmetric part of the velocity gradient tensor.We can also identify q with a two-dimensional vector: q = (q 11 , q 22 ) where q 11 = s 2 cos 2θ; We assume that all quantities of interest only depend on y i.e., we work in a reduced one-dimensional setting, which is a physically relevant setting for very long channels with small height.Equations ( 2) and ( 4) can be recast in terms of the order parameter (s) and the director angle (θ) in one spatial dimension as Using the following scalings, Equations ( 5)-( 7) can be non-dimensionalised as where and L is the half-height of the channel.Physically, L * is the scaled elastic constant.The parameter L 1 = Re/Er * where Re is the Reynolds number and Er * = u 0 Lγ/κ is analogous to the Ericksen number (Er = u 0 Lµ/κ) in terms of the rotational diffusion constant, γ, rather than the viscosity µ.It can also be interpreted as the ratio of the inertial to rotational forces.The parameter L 2 = (2a/c) (Er * /Er) is the product of ratio of the temperature and material constants and the ratio of the rotational to momentum diffusion.The boundary conditions for s and ũ are This simply means that we assume the nematic molecules are perfectly ordered at ỹ = ±1 and we impose the typical no-slip boundary conditions on ỹ = ±1.For the nematic director, we look at two different cases: (i) symmetric Dirichlet conditions for θ on ỹ = ±1 consistent with strong anchoring and in the spirit of [8], (ii) a Neumann-type boundary condition modelling weak anchoring on ỹ = −1 accompanied by a Dirichlet condition on y = 1 as shown below. Asymmetric: where is the winding number and B is a rescaled anchoring strength (see Figure 1 for a sketch of these configurations).It is worth pointing out that positive B models tangential or planar boundary conditions on the bottom substrate i.e., it originates from a surface energy of the form A sin 2 θ that is minimised by either θ = 0 or θ = π, for a positive anchoring coefficient A that measures the strength of the anchoring and we integrate over the surface ỹ = −1.The initial conditions for the system above Equations ( 9)-( 11) are given by s ( ỹ, 0) = 1, ( 18) where θ (−1, 0) is the root of the equation Bθ (−1, 0) − sin [2θ (−1, 0)] = 0 for the asymmetric case.
Figure 1.Schematic of the director orientation in equilibrium when applying (a) symmetric anchoring conditions Equation ( 15) and (b) asymmetric anchoring conditions Equation ( 16).
We will often make comparisons between situations with no flow to situations with fluid flow.In the no-flow case, we simply set ũ( ỹ, t) = 0 in Equations ( 9) and (10) and analyse the resulting system with the same boundary (Equations ( 13), ( 15) and ( 16)) and initial conditions (Equations ( 18) and ( 19)).

Results
The numerical computations are carried out using the finite-element-based commercial package COMSOL v5.2 [30].Recent numerical work shows that the system ( 9)-( 11) subject to the boundary conditions ( 13), (14) and either (15) or ( 16) along with initial conditions ( 18)- (20) has multiple solutions, which can be interpreted as multiple steady-state solutions or locally stable states for the long-time dynamics.Hence, our numerical results are not exhaustive but serve to illustrate interesting solutions that can be attained under certain model situations, particularly with reference to flow-reversal solutions for the velocity profiles.

Comparison of the Flow and No-Flow Situation
We neglect time dependence or transient dynamics in this section and focus on the long-time equilibrium profiles of s, ũ and θ in this section.We fix the parameters L * = 10 −3 and L 1 = 10 −6 , which are physically relevant values from the typical values of material constants reported in the literature and investigate the effects of the parameters, px and L 2 on the solution profiles for Equations ( 9)- (11) [9,24].The results are presented in Figures 2 and 3 where we plot the no-flow profiles for s and θ for reference and then compare these profiles to the distorted profiles with a non-zero pressure gradient px .In Figure 2, we study the effect of the ratio px /L 2 on the spatial profiles of s, ũ and θ with the symmetric Dirichlet boundary conditions, (Equation ( 15)).Here, the solution profiles are symmetric around ỹ = 0 due to the imposed symmetry of the boundary conditions.For L 2 /| px | 1, it is relatively straightforward to see that the flow profile is simply the parabolic Poiseuille form.We confirm this observation in Section 4, where we determine the precise form of ũ, along with s and θ, via a systematic asymptotic analysis of the set-up in the limit L 2 /| px | 1.For very small values of px , the flow is weak as expected.In Figure 3, we do the same for asymmetric boundary conditions (Equation ( 16)).Naturally, the profiles of θ are not symmetric around ỹ = 0 in this case.The asymmetric behaviour in θ is weak, but more pronounced for larger values of L 2 .The profiles of s and ũ are largely unaffected by the asymmetric boundary conditions for θ, at least for the parameter values employed in this section.Further, we can also compute the total fluid flow rate as well as the wall shear stress, which is related to the skin friction coefficient, C s , by The skin friction coefficient represents the friction drag exerted by the wall, which resists the fluid movement.
In Figure 4, we plot the total volumetric flow rate as a function of L 2 for two different values of px = ±1, Dirichlet conditions on ỹ = ±1 and two different values of B. The results suggest that, for L 2 ≥ 1, the net flow rate is greater for px = −1 compared to px = 1.This is the richest regime where both the pressure gradient and the liquid crystal parameter L 2 influence fluid flow and it would be interesting to investigate how p x and L 2 couple together in fluid flow profiles.Moreover, the effect of the wall-alignment constant B on the flow properties can be inferred from Figure 4.This suggests that by altering the wall-anchoring properties (manifested through B) one can manipulate the flow rate and the skin-friction losses.For example, positive B corresponds to preferred tangential anchoring on ỹ = −1 and negative B indicates preferred normal/homeotropic boundary conditions on ỹ = −1.
Since θ = π 2 on ỹ = 1 (ω = 1 2 ), we have homeotropic boundary conditions on ỹ = 1.These results suggest that the net flow rate and the wall shear stress are enhanced by Dirichlet conditions or mixed tangential conditions on ỹ = −1 along with normal boundary conditions on ỹ = 1.We emphasise that the wall shear stress τw is computed on ỹ = −1 where the mixed boundary condition is imposed and the other boundary wall will have a different magnitude (or even direction) of τw associated with it.A phase-space plot of the parameters that correspond to net zero flow rate is shown in Figure 5, also demonstrating the effect of the wall anchoring conditions.For illustrative purposes, we take px = 1.The combination of parameters in the region below the curves in Figure 5 corresponds to fluid flow in the direction of decreasing hydrodynamic pressure, while the region above the curves is in the direction of increasing hydrodynamic pressure, which is a manifestation of backflow.A similar situation, in which a net zero fluid flow rate can be observed, is in electro-osmotic flows, for a critical electrical field strength that exactly balances the hydrodynamic driving pressure [31].The results in Figure 5 provide quantitative estimates for the onset of flow reversal for a specific choice of parameters.
A more exhaustive study on these lines can predict the onset of flow reversal for experimentally relevant or applications-oriented modelling scenarios and flow reversal or tunable flow directions offer new possibilities for topological defects and transport phenomena in microfluidic channels.We do not explore this further in this manuscript.The curve (solid) corresponding to B = 1/3 is for the asymmetric boundary condition (Equation ( 16)).The dotted curve is for the case of the symmetric (Dirichlet) boundary condition (Equation ( 15)).The values of the parameters used are px = 1, ω = 1/2 and L 1 = 10 −6 .

Effect of the Winding Number ω
The impact of the winding number ω on the long-time profiles for the symmetric (Equation ( 15)) and asymmetric case (using Equation ( 16)) respectively, are shown in Figures 6 and 7 for L * = 10 −3 .As ω increases, the energetic penalties for distortions in θ increases (this can be seen by re-scaling θ = ω θ in Equation ( 9)), so the long-time θ-profiles become more linear as ω increases.The gradient θ ỹ is usually maximum in magnitude at ỹ = 0 and consequently, s is a minimum at ỹ = 0 since the energetic penalty is proportional to s2 .Further, the minimum value of s decreases as ω increases, again for similar reasons in the sense that the order decreases to compensate for more distortion in the θ profiles.It is interesting that the total flow rate decreases as ω increases since the flow meets greater resistance from the increasingly distorted θ profiles.It seems difficult to extract this behaviour from a simple analysis of the governing partial differential equations.

Effect of the Parameter L *
The regime of small L * is well understood in the no-flow case.Here, s is approximately unity everywhere and θ is linear with Dirichlet boundary conditions without flow.When a pressure gradient is imposed, one can perform heuristic calculations to predict that the flow profile is approximately parabolic for small L 2 , as in Figures 8c and 9c.(See Section 4 for a more detailed analysis).In the case of symmetric Dirichlet boundary conditions, as L * becomes larger, θ becomes approximately constant everywhere except for a jump at ỹ = 0 to enable the boundary conditions to be satisfied.This can be seen from Equation ( 9) that as L * increases, θ ỹ tends to zero almost everywhere and there is reduced energetic penalty associated with deviations from s = 1.In fact, s = 0 when θ has a jump discontinuity, to regularise the solution.This is referred to as order reconstruction in the liquid-crystal literature when the system interpolates between two fixed boundary conditions by not rotating the eigenframe but by switching the leading eigenvalue at the centre of the cell.However, whilst the order-reconstruction phenomenon is relatively well understood without flow effects, it is far less studied with flow effects.We make certain observations here.For the parameter choices in Figures 8, the θ profile switches from a continuous solution to a discontinuous solution at ỹ = 0 at L * ≈ 0.1.This has a very interesting effect on the flow profile (see Figure 8c) in the sense that there is a distinct region in the channel interior where ũ < 0 for L * = 0.1 and the flow field has a cusp-like minimum at ỹ = 0.For larger values of L * , when the system has settled into the order-reconstruction regime, the flow profiles are less surprising and have the usual parabolic-like form.We point out that θ is not a constant on either side of the jump discontinuity for order reconstruction with flow, in contrast to order reconstruction without flow.The qualitative features are unchanged with asymmetric boundary conditions, see the results in Figure 9c.
For a select range of values of L * , both the continuous and discontinuous solutions for θ are attainable, with the state achieved dependent on the initial condition.We analyse this further in Section 3.5.
There is evidence that the local fluid flow switches direction (at least locally) for certain choices of L * and we have investigated the impact of L * on the net fluid flow rate, as shown in Figure 10.As we have seen in Figure 5, for a given L 2 large enough, there exists a critical L * (L * crit,1 ) for which there is zero net flow.However, here we find that there is a second critical L * (L * crit,2 ), beyond which the flow switches back to the direction of the decreasing pressure.The critical scaled elastic constants L * crit,1 and L * crit,2 have almost the same values for both the symmetric and asymmetric boundary conditions, for the parameter choices in Figure 10.The critical values are relatively large however, and hence unlikely to be attained in most applications.
Asymmetric Symmetric

Dynamic Evolution of the Spatial Profiles
We briefly examine the dynamic evolution of the director profile, order parameter and the fluid flow profiles in Figure 11 (symmetric case) and Figure 12 (asymmetric case).We note that, even though L 1 1 in our simulations, the velocity is time dependent because s and θ are time dependent.
The dynamics are not particularly interesting for this choice of parameters but illustrate how s assumes a U-shaped profile with a shallow minimum as θ evolves from the perfectly linear initial condition, under the effect of flow.The initial flow profile is Poiseuille and the nematic effects suppress the flow profile and distort the parabolic shape.16)).The legends of all the sub-figures are the same as in (a).

Effect of the Initial Condition
Next, we make some preliminary comments on the effect of the initial condition on the equilibrium solution.As noted in Section 3.3, as we increase the parameter L * , the θ solution transitions from a continuous to a discontinuous profile, at a critical value of L * referred to as L * switch .This motivates the question of whether, by an appropriate choice of initial condition, there are parameter regimes that admit multiple steady-state solutions with a basin of attraction.
In Figure 13, we consider a specific case of no fluid flow and symmetric Dirichlet boundary conditions in θ (Equation (15)) and show that the system can indeed exist in multiple steady states.We find that there is a window of values of L * < L * switch ≈ 0.0335 for ω = 1/2 (and L * switch decreases as we increase ω) for which a continuous and discontinuous steady-state solution can be achieved, depending on the initial conditions, indicating that multiple steady states may only be possible in some parameter regimes.The continuous solution is stable in this parameter regime and the discontinuous solution is unstable with respect to perturbations near the centre of the cell.This is consistent with theoretical work in the field.The order-reconstruction or discontinuous solution exists for all values of L * , for our choice of symmetric Dirichlet conditions, with no flow.It is the unique solution for suitably large L * and unstable for suitably small L * [18].However, the instability only manifests in certain directions, so that, for an appropriate choice of initial condition, we can recover the discontinuous solution for smaller values of L * .As L * increases, we recover the order-reconstruction or discontinuous solution for all initial conditions.These results are promising in the context of bistable devices, particularly if the order-reconstruction or discontinuous solution can be "stabilised" by an appropriate control and we have two stable solutions-the continuous and the discontinuous solution for small values of L * .

Steady-State Analysis
In this section, we analytically study the system Equations ( 9)-( 11) in steady state.We assume that L 2 /| px | 1 so that the flow affects the nematic orientational ordering but not vice versa, and that a uniform pressure gradient p x = −G is applied.In particular, this regime does not capture backflow where the nematic order affects fluid flow.
We integrate Equation (11) with respect to ỹ twice and apply boundary conditions (14) to give the following leading-order Poiseuille-type solution for ũ, Substituting Equation (26) into Equation ((10)) and rearranging we obtain We are thus left to solve Equations ( 9) and ( 27), subject to Equation ( 13), and either Equation (15) or Equation ( 16).The numerics have uncovered the possibility for two types of steady solution: continuous or discontinuous solutions in θ.We will study each of these in turn in the following subsections.

Continuous Solutions in θ
We first study the symmetric strong-anchoring regime or Dirichlet boundary conditions for θ.Integrating Equation ( 27) with respect to ỹ twice and applying the boundary conditions (15) gives an explicit expression for θ in terms of s: where Substituting for θ ỹ in Equation ( 9) using Equation ( 28), we obtain the integro-differential equation for s, This must be solved subject to the boundary conditions ( 13).An analogous procedure follows for the asymmetric anchoring conditions (Equation ( 16)), but we do not present the details here.

Small-L * Limit
It is observed that continuous solutions are stable for L * 1.We thus explore the system in this reduced regime.In this case, the leading-order solution in L * to Equation (30) can immediately be seen to be s = 1.As a result, Equation (28) yields the corresponding leading-order solution for θ, A similar method in the asymmetric case, Equation ( 16), yields the leading-order solution with c 2 satisfying the transcendental equation We note also that in the small-L * limit, we may relax the assumption that L 2 is small.In this case, the flow profile is still parabolic to leading order, but is given by We recall that L 2 > 0 is positive since we are working with low temperatures, so a > 0. Negative values of L 2 describe higher temperatures for which s ≈ 1 does not hold.We also note that the flow profiles in the preceding section do not agree with the perfectly parabolic profile described above.This is largely because L * is not sufficiently small in the simulations for the sake of computational efficiency.

Discontinuous Solutions in θ
We now study the case to allow for discontinuities in θ.On physical grounds, s = 0 vanishes at such discontinuities to "regularise" the discontinuities.While such point discontinuities may appear anywhere within the domain, for illustrative purposes we consider the case where a single point discontinuity in θ is present, at ỹ = 0. We focus on the symmetric strong-anchoring regime, but again note that similar methods apply to the asymmetric boundary conditions.We solve in the domain 0 < ỹ ≤ 1 and replace the boundary conditions (Equations ( 13) and ( 15 Since L 2 1, the velocity profile is not influenced by the discontinuities and is still given by Equation (26).
Integrating Equation ( 27) and applying the modified boundary conditions, we find that θ is now given by θ = ωπ + G where c 3 is given by for a given L * .Note that, equivalently, Equation (40) could be viewed as providing explicitly the value of L * corresponding to a particular chosen value of c 3 .The solution for −1 ≤ ỹ < 0 is found by an odd reflection of the solution in 0 < ỹ ≤ 1.
When the pressure gradient G 1, the system possesses a distinguished limit when L * = O(1/G 2 ).(Note we assume that L 2 G so that the second term on the right-hand side of Equation (11) can still be ignored).In this relatively simple case, the equations are amenable to asymptotic analysis, and we are able to write the solution for s and θ explicitly as

Conclusions
In this paper, we investigate the nematic order parameter (captured by the director orientation θ and the order parameter s) and flow profiles in a one-dimensional microfluidic channel, with Dirichlet boundary conditions and mixed boundary conditions for θ, as a function of the pressure gradient, the boundary conditions themselves (in terms of winding number ω and a measure of the anchoring strength B), the nematic elastic constant (L * ) and the scaled viscosities (L 2 ) in a reduced Beris-Edwards setting.For small L 2 , we can analyse the system and obtain at least semi-explicit solutions for the nematic order parameter and the flow profile, both with and without an applied pressure gradient.We consider continuous and discontinuous profiles for θ separately, again including the effect of the pressure gradient.In the discontinuous case, θ is effectively piecewise constant (without flow) for such solutions and discontinuities in θ are regularised by isotropic points with s = 0. We can analytically construct solutions with multiple discontinuities although we suspect that these solutions lose stability with respect to higher-dimensional perturbations.The analytical results set the scene for some interesting control problems on how to stabilise discontinuous solutions for small L * and these discontinuous solutions could offer interesting examples of domain walls with s = 0 in three dimensions.
Our most interesting observations include the onset of flow reversal in these model microfluidic systems.We compute specific criteria for flow reversal (flow in the direction of the pressure gradient) as a function of L * and L 2 and in particular, based on the results in Figure 10, we expect the curve in Figure 5 to fold back on itself, so that for a given L 2 large enough, flow reversal only occurs for a certain range of values L * and not in the entire region above the dotted and solid curves in Figure 5.The observed flow reversal is a distinct manifestation of backflow and only occurs for large enough L 2 .We plan to investigate discontinuous order-reconstruction solutions in the presence of flow and backflow in microfluidic channels as a function of temperature (treating a as a parameter or accounting for cases when L 2 changes sign), geometrical dimensions and the anchoring coefficient B in subsequent work.

Figure 2 .Figure 3 .
Figure 2. The effect of the fluid flow on (a) the director orientation (θ), (b) the order parameter (s), and (c) the velocity ( ũ), at equilibrium, for the case of symmetric boundary condition (Equation (15)).The values of the parameters used are ω = 1/2, L * = 10 −3 and L 1 = 10 −6 .(Here, and elsewhere, we plot the profiles at t = 10, after which time we find the solutions have relaxed to a steady state from the initial configuration, Equations (18)-(20).)Analytic solutions are given in Sections 4 and 4.1.

Figure 4 . 1 − 1
Figure 4. Plot of (a) the total volumetric flow rate and (b) the wall shear stress (which relates to the skin friction coefficient through Equation (25)) as a function of L 2 for different values of the constant B in the symmetric (Dirichlet) case, Equation (15), and the asymmetric case, Equation (16), for θ.The total flow rate is scaled with the equivalent Poiseuille flow rate for a Newtonian fluid, 1 −1 ũ| L 2 =0 d ỹ = −2 px /3.The solid and the dotted lines correspond to the negative and positive values of px respectively.The values of the parameters used are | px | = 1, ω = 1/2, L * = 10 −3 and L 1 = 10 −61.

Figure 5 .
Figure 5. Phase space plot of the parameters (L * and L 2 ) for no overall mass flow rate.Here px > 0.The curve (solid) corresponding to B = 1/3 is for the asymmetric boundary condition (Equation (16)).The dotted curve is for the case of the symmetric (Dirichlet) boundary condition (Equation (15)).The values of the parameters used are px = 1, ω = 1/2 and L 1 = 10 −6 .

Figure 6 .
Figure 6.The effect of the winding number ω on (a) the director orientation (θ), (b) the order parameter (s) and (c) the velocity ( ũ), at equilibrium ( t = 10).In this case, we have considered the symmetric boundary condition in θ (Equation (15)).The values of the parameters used are L * = 10 −3 , px = −10, L 2 = 1, and L 1 = 10 −6 .The legends of all the sub-figures are the same as in (a).

Figure 8 .
Figure 8.The effect of the parameter L * on (a) the director orientation (θ), (b) the order parameter (s) and (c) the velocity ( ũ), at equilibrium ( t = 10), in the case of the symmetric boundary conditions for θ (Equation (15)).The values of the parameters used are ω = 1/2, L 2 = 1, px = −10 and L 1 = 10 −6 .The legends of all the sub-figures are the same as in (a).

Figure 13 .
Figure13.Equilibrium profiles of (a) the director orientation, θ and (b) the order parameter s for two different initial conditions for s and a linear initial profile for θ, without any fluid flow.The blue curves are the equilibrium profiles for θ and s for s( ỹ, 0) = ỹ2 and the green curves are the equilibrium profiles for s( ỹ, 0) = 1.We have considered the symmetric condition for θ given by Equation(15).The values of the parameters used are ω = 1/2 and L * = 0.03.
(35).When there is no external flow, px = G = 0 and Equation (37) gives simply θ = ωπ.The solution for s is then given implicitly from Equation(30