Defect-influenced particle advection in highly confined liquid crystal flows

Abstract

We study the morphology of the Saturn ring defect and director structure around a colloidal particle with normal anchoring conditions and within the flow of the nematic host phase through a rectangular duct of comparable size to the particle. The changes in the defect structures and director profile influence the advection behaviour of the particle, which we compare to that in a simple Newtonian host phase. These effects lead to a non-monotonous dependence of the differential velocity of particle and fluid, also known as retardation ratio, on the Ericksen number.

---

1 Introduction

Microfluidics is concerned with the manipulation and control of fluid flow at the microscale and sits as a multidisciplinary field at the intersection of physics, engineering, and biology. During the last two decades it has seen a tremendous rise in importance as it entered the main stream with many practical applications ranging from medical diagnostics^1,2 and drug delivery^3,4 to chemical synthesis⁵ and lab-on-a-chip technologies.^6–8

The essential physics of a microfluidic system is dictated by a competition between various phenomena, which is captured by a series of dimensionless numbers expressing their relative importance.⁹ The Reynolds number, for instance, is often considered to be small in microfluidic applications. However, relatively recently a focal point has been on inertial microfluidics,^10–12 which gives rise to some interesting and counter-intuitive phenomena.

When particle-laden flows are considered, aspects of confinement become central to microfluidics. Characteristic dimensions of channels and chambers are in the range of tens to hundreds of micrometres, and can be comparable to the size of the particles being transported, and interactions between the fluid and solid boundaries become increasingly important. In the simplest case the fluid is Newtonian, without internal structure, and the geometry is a uniform duct. The first theoretical results for migrating, rigid spheres in unidirectional, two-dimensional flow were provided by Ho and Leal.¹³ These were later extended to three dimensions and refined by Ganatos¹⁴ and Staben,¹⁵ who also verified their theoretical results with an experimental study.¹⁶ Owing to the rise of computer power and the advent of sophisticated simulation methodologies, increasingly complex geometries can now be investigated.¹⁷

For similar reasons, the study of particle advection (i.e. movement in flow direction) in non-Newtonian fluids with internal order structure has evolved only slowly. Suspensions of colloidal particles in a nematic liquid crystalline host phase may serve here as prototype of systems that cannot be described with standard continuum theories, but require additional order parameters to capture the microstructure and its change under flow conditions. After first theoretical studies on the drag of colloidal particles in nematic hosts by Stark,¹⁸ a main research focus in colloid-liquid crystal suspension has been on confinement effects^19–21 and topology properties.^22,23 Confinement effects, often in combination with the behaviour in external electric fields, were also investigated in liquid crystalline emulsions²⁴ as well as in droplets and shells.²⁵ The study of flowing liquid crystals covered primarily pure and confined phases^23,26–29 as they appear frequently in microfluidic setups.

On the theoretical side, studies have been extended to multiple, explicitly resolved particles and the full nemato-hydrodynamic problem that solves for the tensor order parameter and velocity field.³⁰ These approaches are now complemented by various simulation methods.^31–33 Recently, the dynamics of anisotropic particles in nematic liquid crystals under shear flow was investigated.³⁴ However, the advection of colloidal particles in pressure driven flow and extreme confinement, which forms the focus of this work, has to our knowledge so far not been addressed.

Our paper is organised as follows: Section 2 introduces our theoretical modelling framework, the Landau–de Gennes free energy and Beris–Edwards model, while Section 3 gives details of our lattice Boltzmann simulation method. Section 4 shows simulation results for scalar order parameter and director field at various confinement ratios and flow velocities and them to those obtained for simple Newtonian fluids. Section 5 summarises our results and conclusions.

2 Theory

2.1 Landau–de Gennes free energy

The local order of the liquid crystal is described by a traceless and symmetric second-order tensor Q(r,t).^35,36 Its largest eigenvalue q < 2/3 is referred to as the scalar order parameter and provides a measure of the liquid crystalline order at a certain position and time. The eigenvector, d, associated with q is called the director and describes the corresponding average orientation of the liquid crystal molecules.

In equilibrium, the liquid crystal order is determined through minimisation of its free energy, commonly described by the Landau–de Gennes free energy functional

which includes the volume contribution f = f_b + f_g, that itself consists of a bulk contribution f_b and a gradient contribution f_g, and a surface contribution f_s. The bulk free energy density is given bywhere we use the Einstein summation convention, so that Greek indices that appear twice are summed over. In eqn (2), A₀ is a constant that sets the overall energy scale and the parameter γ controls the temperature difference from the isotropic–nematic transition, and is related to a reduced temperature τ byFor γ > 3 or τ < 0 the ordered, nematic state is the equilibrium phase, whereas for 2.7 ≤ γ ≤ 3 or 0 ≤ τ ≤ 1 the nematic state is metastable. For γ < 2.7 or τ > 1 the isotropic state is the equilibrium phase.

The gradient free energy density f_g contains the contributions of splay, bend and twist deformations of the director field as well as order-elastic effects due to gradients of the scalar order parameter,

where ∂_α = ∂/∂x_α and ε_ασν is the Levi-Civita symbol in three dimensions. In principle, the elastic constants κ₀ for splay and bend deformations, and κ₁ for twist deformations can be different. However, in our simulations we use the one-elastic-constant approximation κ₀ = κ₁.

The director is assumed to have a preferred normal orientation to the wall surfaces and to the surface of the colloidal particle, known as a homeotropic anchoring, and is described using a surface free energy term

where w is the surface anchoring strength with values w_wall and w_part at the wall and particle surfaces, respectively. The preferred orientation Q⁰_αβ is assumed to be uniaxial and is given bywhere n is the surface unit normal, δ_αβ is the Kronecker delta and S₀ is the preferred surface scalar order parameter given by

The anchoring strength at the surface of a colloidal particle is often compared to the bulk fluid elastic constant by means of the dimensionless parameter

where R and κ are the radius of the particle and the elastic constant, respectively. For small values of this parameter, the presence of a particle surface should have little impact on the local bulk liquid crystalline ordering.

2.2 Beris–Edwards model

The time evolution of Q_αβ is governed by the Beris–Edwards equation³⁷

∂_tQ_αβ + ∂_π(u_πQ_αβ) + S_αβ(W,Q) = ΓH_αβ, (9)

where ∂_t = ∂/∂t, u is the flow velocity, S(W,Q) denotes the response to shear, W is the velocity gradient tensor, H is the molecular field and Γ is a mobility parameter. The shear term is given bywhere and are the symmetric and antisymmetric contributions to the velocity gradient tensor , respectively, and ξ is the so-called flow alignment parameter, a material constant representing an effective molecular aspect ratio which determines whether the liquid crystal exhibits a flow-aligned state at the Leslie angle or tumbling state. The molecular field H is the functional derivative of the free energy functional with respect to the order parameter,The second term in eqn (11) involving the trace ensures tracelessness of the tensor order parameter as it evolves through eqn (9). This leads to the following molecular field:The governing equations of hydrodynamic motion are the equation of mass conservation, also known as the continuity equation, and the Navier–Stokes equation that describes the balance between the rate of change of linear momentum density and the gradients of the pressure and viscous stresses. In tensor notation they read

∂_tρ + ∂_α(ρu_α) = 0 (13)

and

∂_t(ρu_α) = ∂_βΠ^(LC)_αβ + ∂_βΠ^(HD)_αβ, (14)

respectively. Eqn (13) relates the local rate of change of the density ρ to the advection of mass by the fluid velocity u. Eqn (14) is Newton's second law of momentum change for the fluid and involves the thermotropic stress tensor Π^(LC)_αβ and the hydrodynamic stress tensor Π^(HD)_αβ. The thermotropic stress arises due to the liquid crystal and is given byIn eqn (15), σ_αβ and τ_αβ are the symmetric and antisymmetric stress contributions, respectively, defined aswhere is the isotropic contribution from the nematic liquid crystal to the total pressure, and

τ_αβ = Q_ασH_σβ − H_ασQ_σβ. (17)

The final term in eqn (15) may be expanded asThe hydrodynamic stress tensor is defined as

Π^(HD)_αβ = −pδ_αβ − ρu_αu_β + μ(∂_βu_α + ∂_αu_β) + ζ∂_σu_σδ_αβ, (19)

where μ and ζ are the dynamic and bulk viscosity, respectively. The hydrostatic pressure p is related to the density via an ideal gas equation of state as p = c_s²ρ with c_s as lattice speed of sound as is standard in the lattice Boltzmann method. The last term vanishes in incompressible fluids as eqn (13) becomes ∂_αu_α = 0.

No-slip and no-penetration boundary conditions are applied on the walls and particle surfaces, and the boundary conditions for Q are found from the minimisation of the free energy³⁸

where Q_αβ,γ = ∂Q_αβ/∂x_γ.

3 Simulation method

3.1 Simulation setup

Fig. 1 shows a diagram of the three-dimensional computational geometry, which consists of a duct of L_x = 24, 32 or 48 and L_y × L_z = 256 × 384 lattice sites. Solid walls are positioned at x = 0 and x = L_x, y = 0 and y = L_y. We define the measure of confinement as the ratio of the particle diameter to the height of the duct, which leads in our case to confinement ratios 2R/L_x = 0.8, 0.6, and 0.4. The value of L_y means that, with the particle at the centre of the duct, the system is effectively unconfined in y-direction since 2R/L_y = 0.075. Periodic boundary conditions are applied in the z-direction with the z-boundaries acting as inlet and outlet of the duct. A pressure gradient Ψ = Δp/L_z is applied in z-direction, leading to a body force density acting on all sites.

Fig. 1 Overview of the computational geometry: we apply no-slip and no-penetration boundary conditions and homeotropic anchoring conditions at the walls which bound the region in x-direction and y-direction and at the particle surface, with periodic boundary conditions at the z-boundaries. The top part shows the top view, looking along the x-direction, and the bottom part shows the side view, looking along the y-direction.

We use a hybrid lattice Boltzmann scheme³⁹ that applies a finite-difference method for the dynamics of the Q-tensor order parameter and solves the hydrodynamic part of the problem by means of the lattice Boltzmann method. The colloidal particle is discretised as a solid, mobile particle with a radius R = 9.6. The longitudinal and angular momenta of the colloidal particle are evolved according to Newton's second law of motion. We use a mixed explicit–implicit velocity update, which minimises the number of linear equations that must be solved, while maintaining absolute stability.⁴⁰ On both the walls and the particle surface no-slip and no-penetration boundary conditions are applied by using a bounce-back on links scheme.^41–43 Lubrication corrections are applied normal to the walls within a distance of 0.1 lattice sites.⁴⁰ The surface free energy in eqn (5) invokes a homeotropic anchoring condition with a preferred orientation of the director normal to the surfaces.

There are technical limitations to our model that should be borne in mind. While the centre of mass of the particle is integrated off-grid according to Newton's equation, the particle itself is discretised using a stair-case geometry. This can result in some inaccuracies, especially for highly confined regimes. For instance in our case 2R/L_x = 0.8 only about two lattice sites are between the particle and the walls surfaces at its narrowest point. The ideal gas equation of state that the pressure obeys, as well as the modelling of the constant pressure gradient through an additional body force density on all sites are both common treatments in the lattice Boltzmann methodology and allow for an accurate modelling of a weakly compressible fluid. But the assumption of a constant pressure gradient represents a simplification over the real situation. Thermal fluctuations have not been included since our simulations were carried out using a temperature well away from the isotropic–nematic transition line, and so elastic forces from the anchoring of the liquid crystal dominate over thermal forces by orders of magnitude. The one-elastic-constant approximation is commonly used as first approach and does not compromise our results qualitatively. However, relaxing this approximation could lead to quantitative differences, and potentially also richer phenomenology. The Beris–Edwards model uses a simplified approach to viscosities compared to the Ericksen–Leslie theory, which has six viscosity coefficients α₁, …, α₆³⁵ (only five are independent as the Parodi relation applies⁴⁴). The viscosities in the Beris–Edwards model are implicitly given through the isotropic dynamic shear viscosity μ, the rotational diffusion constant Γ, the flow alignment parameter ξ, and the scalar order parameter S₀. They can be directly related to the Ericksen–Leslie viscosities α₁, …, α₆,³⁹ but parameterise only a subset of possible values.

The simulations were run with our lattice Boltzmann code for complex fluids Ludwig version 0.15.0.⁴⁵ A typical simulation is first initialised with no applied pressure gradient for 5 × 10⁴ iteration steps for each anchoring strength. After this initial equilibration phase, the simulations are restarted with various pressure gradients that are kept constant for 4 × 10⁵ iteration steps. Typical runtimes are approximately 26 hours using a hybrid MPI/OpenMP parallelisation with 2 MPI tasks each running on 20 OpenMP threads. The overview of the simulation parameters is included in Table 1. For further information about the exact implementation used in this work, we guide the reader to the Ludwig code repository and related literature.^46,47

Table 1 Overview of simulation parameters

Bulk energy scale	A 0	0.01
Effective temperature	τ	−0.29
Elastic constants	κ 0, κ1	0.01
Wall anchoring strength	w wall	0.02
Particle anchoring strength	w part	0, 0.001, 0.01, 0.05
Anchoring parameter	ω	0, 0.96, 9.6, 48
Flow alignment parameter	ξ	0.7
Mobility parameter	Γ	0.5
Density	ρ	1.0
Shear viscosity	μ	5/6
Bulk viscosity	ζ	5/6
Particle radius	R	9.6

---

3.2 Parameter mapping

Our simulation units can be mapped to physical units by calibrating the units of pressure, time, and length. To achieve this, we assign the lattice spacing Δx, the algorithmic time step Δt and the reference pressure p* from unity in lattice Boltzmann units (LBU) to their corresponding values in SI units. The principle of this parameter mapping was also shown in our previous work³³ using a different characteristic length scale.

The calibration of the length scale is straightforward as it is simply set by considering the dimensions of the microfluidic duct. If we associate the narrowest gap size L_x = 24 in LBU corresponds to in SI units, we obtain an LBU length of in SI units.

To obtain the pressure scale, we use the measurements of the Landau–de Gennes parameters³⁶ (Appendix D therein) which give

Using A₀ = 0.01 and γ = 3.1 in our simulations results in a reference pressure of in SI units.

For the calibration of the timescale we use the following formula, which relates the rotational viscosity γ₁ of the director to the equilibrium scalar order parameter q and the order parameter mobility Γ:

We use Γ = 0.5 in LBU and bulk energy density parameters that lead to q ≈ 0.5 since it is assumed that the system is well within the nematic phase. Therefore, the rotational viscosity γ₁ = 1 in LBU. Typical values for liquid crystals in SI units are γ₁ = 0.1 Pa s.³⁵ Together with 1 Pa equating to a pressure of 10⁻⁸ in LBU, we obtain for the algorithmic time step .

The Ericksen number characterises the ratio of viscous to elastic forces and is defined as

where u is a characteristic flow velocity, in our case the velocity at the centre of the duct U_c, η is the dynamic viscosity, Λ is a characteristic length scale which is set by the narrowest gap size L_x (see Table 2 for Λ = 2R, which allows direct comparison with our previous work³³), and κ is the bulk elastic constant of the liquid crystal.

Table 2 Conversion of Ericksen numbers for different confinement ratios using different characteristic length scales, namely the size of the channel (odd columns) or as in ref. 33 the diameter of the particle (even columns)

2R/Lx = 0.4		2R/Lx = 0.6		2R/Lx = 0.8
Er (Lx)	Er (2R)	Er (Lx)	Er (2R)	Er (Lx)	Er (2R)
1.65	0.658	4.38	2.63	6.15	4.92
8.30	3.32	9.84	5.90	10.37	8.30
18.10	7.24	21.25	12.75	17.95	14.36
35.55	14.22	39.18	23.51	22.32	17.86
52.72	21.09	51.86	31.11	64.85	52.16
69.77	27.91	102.12	61.27	86.06	69.10

---

The dynamic viscosity η is calculated as an apparent viscosity, defined as the ratio, η = μΦ₀/Φ of the volumetric flux Φ₀ of a simple Newtonian fluid and the volumetric flux of the liquid crystalline system Φ, through a plane perpendicular to the flow in the z-direction, namely

with the flow being driven through the pressure gradient Ψ = Δp/L_z with Δp being the pressure difference between inlet and outlet. The volumetric flow rate Φ₀ of a Newtonian fluid with dynamic viscosity μ through a gap L_x driven by a pressure gradient Ψ in plane Poiseuille flow can be calculated as

Fig. 2 shows fluid velocity profiles for a representative confinement ratio of 2R/L_x = 0.6 that have been normalised to the peak flow velocity of a simple Newtonian fluid in Poiseuille flow at the same pressure gradient and scaled using the x-dimension of the duct. The apparent viscosity η is the ratio of the areas under the Poiseuille curve and the curves at finite Ericksen numbers.

Fig. 2 Scaled magnitude of the fluid velocity u(x) = |u(x)| normalised against the peak flow velocity u_c,Poiseuille of a simple Newtonian fluid in Poiseuille flow at the centre line of the duct at x = L_x/2. The image shows representative results for the confinement ratio 2R/L_x = 0.6. The black line is the parabolic flow profile of Poiseuille flow. Away from the walls the velocity profiles of the flowing nematic are parabolic and deviations from the parabolic profile occur only close to the walls. The inset shows the dependence of the centre line fluid velocity u_c on the Ericksen number Er. Profiles for the other confinement ratios are not shown as they look very similar once normalised.

More specifically, in Fig. 2 the flow velocities have been normalised against the maximum flow velocity of the Poiseuille flow u_c,Poiseuille(x = L_x/2) = L_x²Δp/8μL_z at the centre line of the duct with μ as dynamic viscosity and Δp as pressure difference between inlet and outlet, respectively. Away from the walls at x/L_x = 0 and x/L_x = 1 the velocity profiles of the flowing nematic are parabolic and deviations form the parabolic profile occur only close to the walls. This is a result of shear thinning as the director field flow-aligns further away from the walls, which is prevented by the normal wall anchoring in the vicinity of the walls. In all simulations that contain a colloidal particle the fluid flow velocity u_c at the centre line was taken at x = L_x/2 and a point in a distance L_z/2 upstream/downstream from the particle, which is the point farthest away from the particle in the z-direction due to the periodic boundary conditions. However, owing to this large distance the values we obtained for u_c in this manner are virtually identical to those of a pure liquid crystal without particle. Profiles for other confinement ratios are not shown as they look very similar.

4 Results and discussion

We study the advection behaviour of a single particle moving in a nematic host phase in highly confining ducts and investigate the effect that varying pressure gradient Ψ, confinement ratio 2R/L_x and homeotropic anchoring strength have. In a simple Newtonian fluid, or in a liquid crystal at temperatures above the isotropic–nematic transition point, the motion of a freely suspended spherical particle between two parallel plane walls has been studied previously theoretically,^14,15 with simulations¹⁷ and experimentally.¹⁶ The main effect is that the retardation of the particle motion to the fluid motion is primarily independent of the applied pressure gradient, but greater for particles closer to either of the walls, and therefore more so for highly confined particles due to the proximity to the walls.

In a nematic liquid crystal with homeotropic anchoring conditions at the walls the director orientation is forced to be parallel to the wall normals. The degree of alignment depends on the strength of the anchoring, but also on the velocity gradient, and therefore the pressure gradient. At low pressure gradients, the nematic order will be enforced throughout the duct. But for higher pressure gradients the director field flow-aligns at the Leslie angle. Two conformations are persistent in flowing nematics, namely the so-called bend state or H-state and the splay state or V-state. For both H- and V-state the director flow-aligns to a positive (negative) Leslie angle in the lower (upper) half of the channel. The difference between the two states is determined by the way the director rotates between the positive and negative Leslie angles at the centre: in the bend state, the director at the centre is perpendicular to the walls, whereas in the splay state the director is almost parallel to the walls at the centre. The bend state is generally adopted at low flow velocities, whereas the nematic transitions to the splay state at higher flow velocities.

Fig. 3 shows the director field, defect structure and magnitude of the fluid velocity at a medium confinement ratio 2R/L_x = 0.6 and different Ericksen numbers. The left column displays the bend state at Er = 4.38 prior to the transition to the splay state, whereas the right column shows the splay state at Er = 51.86 after transitioning from the bend state. The two top rows contain slices in the xz-plane (narrowest duct dimension and flow direction) at y = L_y/2 with walls at the x-boundaries at the top and bottom, whereas the two bottom rows show slices in the yz-plane (widest duct dimension and flow direction) at x = L_x/2 cropped to the vicinity of the colloidal particle.

Fig. 3 Director field, defect structure and fluid velocity profiles for confinement ratio 2R/L_x = 0.6 and anchoring parameter ω = 48 before and after the bend-to-splay transition. The left column shows the bend state (H-state), while theright column shows the splay state (V-state). The first and third row show the director field d with the magnitude d_z of its z-component indicated through the colour code. The second and fourth row show the magnitude of the fluid velocity u(x,z) and u(y,z) through the centre of the particle, normalised to the maximum velocity u_c at the centre line of the duct, where arrows give a sense of the vectorial dependence of the fluid velocity field. The images in the two top rows represent slices through the middle of the channel in the xz-plane (narrowest duct dimension and flow direction) and have the view along the negative y-dimension. Those in the two bottom rows show slices in the yz-plane (widest and narrowest duct dimension) and have the view in positive x-direction. The flow direction is from left to right in positive z-direction. The opacity of the defect rings (green isosurfaces) has been slightly reduced to enhance the visibility of the local director field.

The director field in the first and third row is colour-coded with red indicating an orientation parallel to the flow or z-direction and blue indicating an orientation perpendicular to the flow direction or in xy-plane. The bend state at low Ericksen number shows the Saturn ring defect oriented parallel to the walls with only very minor deformations, while the splay state has the Saturn ring defect oriented approximately perpendicular to the walls and displaced slightly downstream from the meridian of the particle in positive z-direction.

The second and fourth row show the magnitude of the fluid velocity u(x,z) = |u(x,z)| in the xz-plane and u(y,z) = |u(y,z)| in the yz-plane normalised to the maximum velocity u_c at the centre line of the duct. It is interesting to see that despite the striking differences in the director field structure and defect ring orientation at the two different Ericksen numbers both flow profiles are very similar. A minor exception is that at the lower Ericksen number the peak velocity is attained very close to the particle, whereas at the higher Ericksen number the relative fluid velocity is slightly reduced around the particle. This is a consequence of the different differential velocities between the colloidal particle and the fluid in both cases (see Fig. 8).

As a quantitative overview of our findings, we include in Fig. 4, snapshots of the particle and its defect in the steady state for varying confinement ratios and Ericksen numbers. In each cell the left images show the side view looking in the negative y-direction with walls at the top and bottom. The images on the right show the view from the top looking in the positive x-direction. The confinement increases from left to right from confinement ratios 2R/L_x = 0.4 to 2R/L_x = 0.6 to 2R/L_x = 0.8, and Ericksen numbers increase from top to bottom. The defect is shown as a green isosurface defined by a local order parameter q ≤ 0.188 and the particle anchoring strength and dimensionless anchoring parameter are w_part = 0.05 and ω = 48, respectively, as lower anchoring strengths do not result in defects that could be distinctively visualised.

Fig. 4 Snapshots of the director field and defect structure in the steady state at various Ericksen numbers and strongest particle anchoring parameter ω = 48. The bright green region corresponds to the defect where liquid crystalline order is reduced. In each cell the left images in each cell show the side view looking in the negative y-direction with walls at the top and bottom. The images on the right in each cell show the view from the top looking in the positive x-direction. The flow is in the horizontal positive z-direction from left to right. The confinement increases from left to right, and Ericksen numbers increase from top to bottom. While in most cases the particle stays at the centre of the duct throughout the simulation, there are a small number of cases where they migrate away from it. Specifically, there are two cases where the particle migrates fully to a wall (for the two highest Ericksen number and confinement ratio 2R/L_x = 0.8), and three cases where the particle migrates to a stable position between the wall and the centre (for Er = 18.10 and 2R/L_x = 0.4, for Er = 9.84 and 2R/L_x = 0.8, and for Er = 10.37 and 2R/L_x = 0.8).

At low Ericksen numbers, below the bend-to-splay transition, the particle has a Saturn ring defect whose ring plane remains parallel to the walls. This is the case for all confinement ratios and Ericksen numbers below Er = 10.37, as shown in the first and second row of Fig. 4. Two aspects are noteworthy: firstly, there is a slight increase of the defect isosurface radius downstream of the particle, for which both Ericksen number (see images for 2R/L_x = 0.4 with Er = 1.65 and Er = 8.30) and confinement ratio (see image for 2R/L_x = 0.6, Er = 4.38 and image and Movie S1 (ESI†) for 2R/L_x = 0.8, Er = 6.15) seem responsible. However, confinement appears to play a more important role in this context.

Secondly, at slightly increased Ericksen numbers (see images for 2R/L_x = 0.6, Er = 9.84 and 2R/L_x = 0.8, Er = 10.37), the Saturn ring becomes angled such that the part downstream of the particle is closer to the bottom wall, while the other part upstream of the particle remains virtually unchanged. These two particular cases reached steady state positions that are offset somewhere between the centre of the duct and the walls in the x-direction, which contributes to this asymmetric appearance. This can be explained with the migration (i.e. lateral movement perpendicular to the flow direction) to the weak attractor region that we observed in our previous work on controllable particle migration³³ in practically unconfined conditions using a much wider duct and lower confinement ratio 2R/L_x = 0.15. For direct comparison we provide in Table 2 an approximate conversion between particle Ericksen numbers, as used our previous publication, and Ericksen numbers based on the smallest duct dimension, as used in this work.

Fig. 5 shows a direct comparison of the defect rings around the particle for the lowest and highest simulated Ericksen numbers below the bend-to-splay transition, at confinement ratios (a) 2R/L_x = 0.4 (b) 0.6 and (c) 0.8. The defects at the lowest Ericksen numbers, depicted in grey, are distinctive Saturn rings that are oriented parallel to the walls in the x-direction. As previously mentioned, increasing confinement leads to a defect ring that is thicker at the downstream side of the particle, while it remains oriented parallel to the wall at the x-boundary. Increasing the Ericksen number alone does not change the orientation of the defect ring, but leads to a very slight shift in position upstream (see grey and yellow defect rings in Fig. 5(a)). However, increasing the Ericksen number and confinement ratio induces a noticeable tilt of the defect ring, shown in Fig. 5(b) and (c), as the particle migrates into the off-centre steady state position somewhere between the centre of the duct and one of the walls. Despite the difference in confinement ratio and Ericksen number (2R/L_x = 0.6 and 0.8, Er = 9.84 and 10.37, respectively) the shape of the defect rings is almost the same.

Fig. 5 Disclination lines around the particle for the lowest and highest simulated Ericksen numbers below the bend-to-splay transition at confinement ratios (a) 2R/L_x = 0.4, (b) 2R/L_x = 0.6 and (c) 2R/L_x = 0.8, respectively. The top row has the view in negative y-direction, the widest duct dimension, with the walls in the narrowest duct dimension at the x-boundaries situated closely above and below the particle. The bottom row shows the view in positive x-direction. The flow is in the horizontal positive z-direction from left to right.

Upon increasing the Ericksen number, a bend-to-splay transition takes place somewhere between 8.30 < Er < 18.10 (for 2R/L_x = 0.4), 9.84 < Er < 21.25 (for 2R/L_x = 0.6) and 10.37 < Er < 17.95 (for 2R/L_x = 0.8). The defect ring is now reoriented with its ring plane approximately perpendicular to the walls and flow direction, as shown in Fig. 4, for instance in the third row, and retains a similar shape at higher Ericksen numbers (see images for 2R/L_x = 0.4, Er = 18.10, 2R/L_x = 0.6, Er = 21.25, and 2R/L_x = 0.8, Er = 22.32). The case for 2R/L_x = 0.4, Er = 18.10 forms an exception in that the particle moves very slightly away from the centre into a stable off-centre position, while in the other cases the particle remains at the centre of the duct, which can be also understood with the migration to the previously observed weak attractor region at similar Ericksen numbers³³ (see Table 2 for conversion of Ericksen numbers). A noticeable difference is that with increasing confinement the defect ring appears compressed in the smallest duct dimension due to the relative proximity of the walls (see image and Movie S2 (ESI†) for 2R/L_x = 0.8, Er = 22.32).

With increasing Ericksen numbers the shape of the vertically oriented defect ring remains largely unchanged for low and medium confinement, as shown in the first and second column, forth and fifth row, of Fig. 4 (2R/L_x = 0.4 and 0.6) for Ericksen numbers Er = 52.72, 69.77 and Er = 51.86, 102, 12, respectively. At Er = 102.12 a slight change occurs such that the defect close to the mid-plane of the duct in the x-direction are distorted and pulled in the upstream direction, i.e. against the flow. This effect is a precursor of the more dramatic elongation of the Saturn ring that will become even more evident as the confinement ratio increases.

At even higher Ericksen numbers Er = 64.85 and Er = 86.06 and confinement 2R/L_x = 0.8, shown in the third column forth and fifth row of Fig. 4, we observe defects that differ substantially from those discussed before. In these cases the particle migrates fully to one of the walls. This has also been previously observed for similar Ericksen numbers in much lower confinement.³³ But there it occurred when the particle was within a distance of one and a half to two diameters from the walls, depending on the Ericksen number (see Table 2 for a conversion of Er). Given the proximity of the walls in the present work with increased confinement, this means that attraction to the walls should occur in practically all situations. This, however, is not the case as we observe attraction to the walls only for the highest Ericksen numbers and the largest confinement. Thus, increased confinement prevents particle migration to the walls and stabilises trajectories around the centre of the duct. The migration to one of the walls results in a different defect shape such that there is a pronounced elongation of the Saturn ring defect in the upstream direction. There is also the indication of a small satellite region of low order, upstream of the particle that never merges up with the rest of the defect (see image and Movie S3 (ESI†) for 2R/L_x = 0.8 and Er = 86.06).

Before focusing on the director structure at high Ericksen numbers and large confinement in more detail (see Fig. 7), we present briefly a synopsis of the defect rings at different confinement ratios and Ericksen numbers. Fig. 6 shows superimposed, vertically oriented defect rings as they occur after the bend-to-splay transition has taken place. At the lowest confinement ratio 2R/L_x = 0.4, shown in Fig. 6(a), the defect ring remains relatively undistorted across a range of medium to high Ericksen numbers. However, comparing the images at the top with the view along the y-direction across the narrowest gap to those at the bottom with the view along the x-direction across the widest gap gives evidence that the shape of the Saturn ring defects is sensitive to confinement. When confined, the defect rings are located slightly downstream from the particle's equator, whereas they remain situated along the equator in the dimension of no or very small confinement (2R/L_y = 0.075). This feature becomes more pronounced as the confinement increases, discernible through the green defect rings at ratios 2R/L_x = 0.6 in Fig. 6(b), and more so at 2R/L_x = 0.8 in Fig. 6(c) where it results in the compressed appearance (Fig. 6(c) top row). This applies to lower (green isosurfaces) and medium (orange isosurfaces) Ericksen numbers. Increasing both confinement and Ericksen numbers leads to the aforementioned different appearance of the defect rings (purple and magenta isosurfaces).

Fig. 6 Saturn ring disclination lines around the particle for various Ericksen numbers after the bend-to-spay transition has taken place. The confinement ratios are (a) 2R/L_x = 0.4, (b) 2R/L_x = 0.6 and (c) 2R/L_x = 0.8. In the top row the view is along the negative y-direction, the widest duct dimension, with the walls in the narrowest duct dimension at the x-boundaries situated above and below the particle. The bottom row is the view in the positive x-direction. The flow is in the horizontal z-direction from left to right.

It is worth mentioning that the confinement ratios we studied are larger than those in similar studies^31,32 (2R/L_x = 0.25 and 2R/L_x = 0.19, respectively), where the lower confinement has been chosen to eliminate possible effects on the results. However, our case of 2R/L_x = 0.4 is obviously already low enough to feature defect rings that appear undeformed and occur at unaltered relative positions to the particle.

The director field, defect structure and magnitude of the fluid velocity at high Ericksen numbers and the largest confinement ratio 2R/L_x = 0.8 are shown in Fig. 7. At this confinement ratio the walls at the x-boundaries are close to the colloidal particle. The two top rows contain slices in the xz-plane (narrowest duct dimension and flow direction) at y = L_y/2 with walls at the x-boundaries at the top and bottom, whereas the two bottom rows show slices in the yz-plane (widest duct dimension and flow direction) at x = L_x/2 cropped to the vicinity of the colloidal particle. The director field in the first and third row is colour-coded with red indicating an orientation parallel to the flow or z-direction and blue indicating an orientation perpendicular to the flow direction or in xy-plane. The left column shows the situation at moderately high Ericksen numbers Er = 22.32.

Fig. 7 Director field, defect structure and fluid velocity profiles for confinement ratio 2R/L_x = 0.8 and anchoring parameter ω = 48 after the bend-to-splay transition for increasing Ericksen numbers Er = 22.32, 64.85 and 86.06, respectively. The first and third row show the director field d with the magnitude d_z of its z-component indicated through the colour code. The second and fourth row show the magnitude of the fluid velocity u(x,z) and u(y,z) through the centre of the particle, normalised to the maximum velocity u_c at the centre line of the duct, where arrows give a sense of the vectorial dependence of the fluid velocity field. The images in the two top rows represent slices through the middle of the channel in the xz-plane (narrowest duct dimension and flow direction) and have the view along the negative y-dimension. Those in the two bottom rows show slices in the yz-plane (widest and narrowest duct dimension) and have the view in positive x-direction. The flow direction is from left to right in positive z-direction. The opacity of the defect rings (green isosurfaces) has been slightly reduced to enhance the visibility of the local director field.

The defect ring is vertically oriented, noticeably displaced downstream close to the walls at the boundary in x-direction (see Fig. 7 first row first column), and situated at the equatorial region of the particle in the non-confined y-dimension (see Fig. 7 third row first column). The director field structure in yz-plane shows that flow alignment occurs in a short distance from the particle and entails a defect in the equatorial region. Focusing again on the director field in xz-plane reveals that the situation is different in the confined x-dimension. Here, the homeotropic anchoring conditions at the wall and particle surfaces prevent any kind of flow alignment in the narrow gap between the particle and the walls. Considering the left-hand upstream side of the particle it becomes evident that both the normal anchoring conditions on the surface and the flow-alignment close to the surface work in the same sense and promote the same director orientation. This is different on the right-hand downstream side. While downstream directly right from the particle's centre flow-alignment and anchoring are also working in the same sense, this is not the case downstream above right and below right from the centre where flow-alignment invokes a northwest–southeast orientation of the director field, while surface anchoring promotes a northeast–southwest orientation. This leads to the slight downstream displacement of the defect ring.

At higher Ericksen numbers Er = 64.85 and Er = 86.06 the particle migrates readily to one of the walls³³ and the shape of the defect changes markedly (see Fig. 7 first and third row, second and third column). The asymmetric positioning of the particle in the duct is only partly responsible for this. In fact, we observe large differential velocity between the particle and the fluid, which means the particle acts now increasingly as obstacle. Therefore, it is instructive to look again at fluid velocity profiles.

The second and fourth row in Fig. 7 show the magnitude of the fluid velocity u(x,z) = |u(x,z)| in the xz-plane and u(y,z) = |u(y,z)| in the yz-plane normalised to the maximum velocity u_c at the centre line of the duct. The profiles in xz-plane (second row) show that compared to Fig. 3 where the confinement ratio is 2R/L_x = 0.6, the now larger confinement ratio of 2R/L_x = 0.8 leads to much lower relative fluid velocities upstream and downstream on the left and right of the particle. With increasing Ericksen number a region with enhanced flow velocities emerges immediately above the particle where the fluid is forced upwards (see Fig. 7 second row third column). The fluid velocity profiles in yz-plane (see Fig. 7 fourth row) demonstrate even further how the relative fluid velocity drops around the particle with increasing Ericksen number. However, what the colour code and normalisation to the peak flow velocity u_c hide is that the velocity gradients in absolute terms are even larger for larger pressure gradient, a direct consequence of higher absolute values of the peak velocity u_c. In view of the director field and defect structure, it becomes evident that the regions with large fluid velocity gradients are also the regions where the director structure becomes noticeably distorted. This effect in combination with local flow-alignment and surface anchoring leads to local regions of low order, for instance the satellite region of very low order slightly upstream on the left of the particle (see Fig. 7 second and third column), and causes the defect ring to become extended further upstream, albeit never completely engulfing the particle.

We conclude our study with an analysis of the advection behaviour of the colloidal particle at different Ericksen numbers and confinement ratios and compare it to that in a simple Newtonian fluid. For this purpose we draw on the theoretical results obtained by Staben et al.,¹⁵ which have been reproduced in a number of studies. While our Reynolds numbers are typically between and and therefore larger than those in ref. 15, it is worth emphasising that the latter results form still a suitable reference as both regimes can be classed as low-Reynolds number.

A suitable measure to characterise the advection behaviour is the retardation ratio v/u_c of particle velocity v to fluid velocity u_c at the centreline of the duct. In an isotropic Newtonian fluid under Poiseuille flow this ratio is constant and depends only on the distance of the particle from the walls of the duct and the confinement ratio. In particular, v/u_c is independent of the Reynolds number. Without confinement the retardation ratio v/u_c is unity as the particle acts as a tracer and is simply advected with the fluid. At finite confinement ratios below 2R/L_x = 1 the movement of the particle is slowed down in the parabolic Poiseuille flow due to the no-slip boundary conditions on the walls of the duct.

Fig. 8 shows the retardation ratio v/u_c for different confinement ratios 2R/L_x, particle anchoring strengths and Ericksen numbers Er. Using the Ericksen number as abscissa has the advantage that the bend-to-splay transition occurs at similar values aiding the comparison across different confinement ratios. The straight, horizontal lines represent the results for a particle in a simple Newtonian fluid. Dashed-dotted lines give the results from Staben et al.¹⁵ for particles at the centre of the duct. We measure retardation ratios of v/u_c = 0.946, 0.876 and 0.759 for confinement ratios 2R/L_x = 0.4, 0.6 and 0.8, respectively, shown in Fig. 8 with solid lines. These results compare well with those of Staben et al., which are v/u_c = 0.945, 0.871 and 0.746 for the same confinement ratios and particles positioned at the centre of the duct. It is worth emphasising that in our setup the largest confinement ratio 2R/L_x = 0.8 has less than 3 lattice sites between the particle surface and the walls on either side. Nevertheless, the relative deviation between ours and Staben's results for Newtonian host phases is less than 1.8% in the worst case, which means our method is remarkably accurate given the relatively sparse discretisation. However, it should be borne in mind that when modelling a liquid crystalline host phases the sparse discretisation affects also the tensor order parameter Q in addition to the fluid–solid interaction in a Newtonian host phase. While these limitations affect the results in Fig. 8 to a certain extent, there are nevertheless clear and robust trends that we will now discuss.

Fig. 8 Comparison of retardation ratios v/u_c of particle velocity v to fluid velocity u_c at the centre of the rectangular duct for confinement ratios 2R/L_x = 0.4 (blue squares), 0.6 (green triangles) and 0.8 (red circles) and different particle anchoring strengths. Horizontal lines show results in a Newtonian fluid from Staben et al.¹⁵ (dashed-dotted lines) and our approach in the isotropic phase (stars). Open symbols indicate cases where the colloidal particle has been fixed in x-direction for comparison as it would normally migrate away from the centre of the duct to either an off-centre position or to the walls. The vertical lines indicate the approximate position of the bend-to-splay transition.

At low Ericksen numbers we observe retardation ratios v/u_c that are close or identical to their corresponding values in Newtonian fluids. Interestingly, and primarily for no or low particle anchoring strength and low confinement ratios 2R/L_x = 0.4 and 0.6, the retardation ratio can be slightly larger in the nematic host phase than in the Newtonian host phase (see light and medium grey and blue data points in Fig. 8). We postulate this occurs because for a particular pressure gradient the peak flow velocity is lower in the flowing nematic than in the Newtonian fluid. But as the same pressure gradients acts across the particle, the latter does not slow down to the same degree in the flowing nematic, leading to comparably higher retardation ratios. For larger anchoring strengths or in higher confinement both additional elastic forces are exerted on the particle and the effective viscosity in the vicinity of the particle increases, both to the effect of slowing down the particle, resulting in smaller retardation ratios.

As the Ericksen number increases, the nematic host phase undergoes a transition from the bend to the splay phase. This occurs at Ericksen numbers 8.30 < Er < 18.10 (2R/L_x = 0.4), 9.84 < Er < 21.25 (2R/L_x = 0.6) and 10.37 < Er < 17.95 (2R/L_x = 0.8), respectively and is indicated by the vertical green dashed lines in Fig. 8. The transition is accompanied by a noticeable drop in the retardation ratio, which reaches a minimum around Ericksen numbers Er ≃ 20, so just beyond the bend-to-splay transition. The minimum is smaller the larger the particle anchoring strength is, but only for medium and large confinement (see medium grey and green curves as well as black and red curves in Fig. 8) and not so for small confinement (see light grey and blue curves in Fig. 8).

Beyond Ericksen numbers in the range of Er ≃ 20 the retardation ratio v/u_c begins to increase again, giving rise to an overall non-monotonous dependence on the Ericksen number. This is the case across all confinement ratios, and the retardation ratios begin to flatten out towards higher Ericksen numbers, approaching or reaching the values of Newtonian fluids again. This non-monotonous behaviour is therefore a consequence of the decreasing importance of liquid crystalline elasticity and consistent with the idea that with higher the Ericksen numbers the liquid crystal behaves rheologically more like a simple fluid.

Regarding how the retardation ratio v/u_c depends on the particle anchoring strength the same tends as for the minima prevail. Higher anchoring strengths entail smaller retardation ratios unless the confinement is small. For our largest confinement ratio 2R/L_x = 0.8 and strongest particle anchoring strength w_part = 0.05 we observe a very strong decrease. This, however, originates also from the migration of the particles to the walls. The two empty circles in Fig. 8 (and similar empty symbols at the two other confinement ratios) permit us to estimate how the trend would continue if the particles had been prevented from leaving the region of maximum flow velocity at the centre of the duct.

In order to explain these findings, we have to look at several separate mechanisms: first of all, there is the transition from the bend to the splay state, which all particles regardless of their anchoring conditions are subject to. The data points for vanishing particle anchoring strength w_part = 0 (light, medium and dark grey in Fig. 8) are indicative of this. The transition causes the general reduction of the retardation ratios from their initially approximately Newtonian values at low Ericksen numbers to their minima around Er ≃ 20. The reason for this decrease is the drop in apparent viscosity and increase in flow velocity u_c around the centre of the duct, whereas the regions of the particle closer to the walls act as anchor and do not allow the particle to pick up velocity v at the same proportion.

The second mechanism at work is the reorientation of the defect ring at the bend-to-splay transition, provided the particle anchoring strength is large enough for a defect to emerge. The vertical orientation of the defect ring with its ring plane perpendicular to the flow direction and walls increases the effective particle radius in the narrowest duct dimension and therefore the effective confinement ratio. This leads to lower retardation ratio v/u_c the larger the anchoring strength is. However, our results suggest this is only the case provided the confinement is not too large. For instance, at 2R/L_x = 0.4 there is very little difference between vanishing and very strong particle anchoring up to Ericksen numbers Er ≃ 60, while at 2R/L_x = 0.6 and 0.8 differences are clearly visible at all Ericksen numbers. This subtlety can be understood by realising that at the different confinement ratios and flow velocities both velocity and order parameter gradients differ across the particle diameter. At a given flow velocity the gradients are largest in large confinement and vice versa. At a given confinement ratio the velocity gradient is largest at large flow velocities and Ericksen numbers. It is precisely this nonlinear order-flow coupling and the interactions between flow and order structure in the vicinity of the particle that cause the observed minor variations in the retardation ratio.

Finally, there is also the possibility of a direct interaction with the wall anchoring when Ericksen numbers and confinement ratios are large. In these situations the colloidal particle shows a tendency to leave the centre of the duct and migrate to the wall regions. There, the advection velocity and therefore the retardation ratio are reduced as a result of the no-slip boundary conditions at the walls.

5 Conclusions

In microfluidic setups of particle suspensions, confinement is often necessary as it allows a certain degree of lateral control over the particle positions, for instance when techniques like confocal or polarised microscopy are used. Our present study has primarily the goal to address some knowledge gaps as to how defects influence and alter the advection behaviour of colloidal particles in moderate and large confinement.

In homeotropic anchoring conditions at the walls and surface of the particle the director field is in the H- or bend state at low Ericksen numbers and has a Saturn ring defect which is oriented parallel to the walls. Increasing the confinement changes the appearance of the defect ring downstream. It can either thicken the defect ring, or invoke a migration to off-centre off-wall positions which we identify with the weak attractor region in our previous work.³³ The latter entails a slight distortion and tilt away from the centre plane.

At moderately high Ericksen numbers around Er ≃ 20 we observe the transition from the bend or H-state to the splay or V-state. This leads generally to a reorientation of the defect ring with a ring plane perpendicular to the flow direction and walls. The defect ring is slightly peeled off downstream in the confined dimension, but sits at the particle's equatorial region in the unconfined dimension, giving it the appearance of an open mouth when viewed from the flow direction. These features are retained at higher Ericksen numbers and in lower confinement. Highly confined particles show a strong tendency to migrate to the walls, a behaviour we observed also in our previous work.³³ This leads to a highly asymmetric defect and induces a satellite region of low order upstream of the particle, which acts partly as an obstacle and forces the flow to slow down and divert around it. Compared to our previous work we do not observe migration to the walls for all but the highest Ericksen numbers and confinement ratios that we tested. Therefore, increased confinement entails stabilisation of trajectories at the centre of the duct.

The interaction between nematic order and flow on one hand, and the fluid–solid interaction on the other hand results in a non-monotonous dependence of the retardation ratio, the ratio of particle advection velocity to the maximum velocity at the centre of the duct, on the Ericksen number. When the Ericksen number is low, the retardation ratio is close to values observed in a Newtonian host phase in all confinement ratios and particle anchoring conditions. This is also the case for vanishing or low anchoring strength and at high Ericksen numbers, where the nematic liquid crystal behaves increasingly like a simple Newtonian fluid as the relative importance of elastic effects decreases. Intermediate Ericksen numbers, however, are characterised by a pronounced minimum in the retardation ratio. We attribute this to a combination of two effects: firstly there is the bend-to-splay transition, to which particles in all anchoring conditions are subject. Secondly, the defect ring undergoes a reorientation from horizontal alignment with the ring plane parallel to the walls to a vertical orientation, which has the ring plane perpendicular to the flow direction and the walls. This increases the effective particle radius and therefore the confinement. The second effect is only present when the defect ring is properly formed, i.e. for stronger particle anchoring strengths, and when the confinement is lower. This is because the increased retardation that the particle experiences is a consequence of the interaction of the defect with the gradients of the flow velocity and liquid crystalline order.

The present study leaves some questions untouched, for example how planar degenerate or hybrid anchoring conditions affects the defect morphology and advection behaviour of the particles. Planar degenerate wall anchoring is fully compatible with the flow alignment that takes place at higher Ericksen numbers. Hence, there is no bend-to-splay transition, rather a more gradual transition to a state where the director field is flow-aligned at the Leslie angle. Furthermore, planar degenerate anchoring conditions on the particle surface invoke topologically different boojum defects that occur at low Ericksen numbers symmetrically upstream and downstream of the particle on an axis that goes through the centre of the particle. Similarly, particles with homeotropic anchoring conditions as used in the present study, but larger diameters, have also topologically different defects. Instead of the half-integer bulk defect loops with topological charge −1/2 they have dipolar full integer satellite defects with topological charge −1 that sit in a distance from the particle surface. It is not clear how these topological differences would affect the advection and migration behaviour.

Data availability

Data underpinning this publication are openly available from the University of Strathclyde KnowledgeBase research information portal.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank Timm Krüger and Uroš Tkalec for very helpful discussions. We acknowledge support from EPSRC under Grant no. EP/R513349/1. M. L. acknowledges funding from the Mac Robertson Postgraduate Travel Scholarship. For the purpose of complying with UKRI's open access policy, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission. This work used the ARCHIE-WeSt High Performance Computer (https://www.archie-west.ac.uk) based at the University of Strathclyde.

Notes and references

1. C. D. Chin, T. Laksanasopin, Y. K. Cheung, D. Steinmiller, V. Linder, H. Parsa, J. Wang, H. Moore, R. Rouse, G. Umviligihozo, E. Karita, L. Mwambarangwe, S. L. Braunstein, J. van de Wijgert, R. Sahabo, J. E. Justman, W. El-Sadr and S. K. Sia, Nat. Med., 2011, 17, 1015–U138.
2. J. Hu, S. Wang, L. Wang, F. Li, B. Pingguan-Murphy, T. J. Lu and F. Xu, Biosens. Bioelectron., 2014, 54, 585–597.
3. Y. Zhang, H. F. Chan and K. W. Leong, Adv. Drug Delivery Rev., 2013, 65, 104–120.
4. S. K. Golombek, J.-N. May, B. Theek, L. Appold, N. Drude, F. Kiessling and T. Lammers, Adv. Drug Delivery Rev., 2018, 130, 17–38.
5. R. Seemann, M. Brinkmann, T. Pfohl and S. Herminghaus, Rep. Prog. Phys., 2012, 75, 016601.
6. D. Mark, S. Haeberle, G. Roth, F. von Stetten and R. Zengerle, Chem. Soc. Rev., 2010, 39, 1153–1182.
7. D. Qin, Y. Xia and G. M. Whitesides, Nat. Protoc., 2010, 5, 491–502.
8. E. K. Sackmann, A. L. Fulton and D. J. Beebe, Nature, 2014, 507, 181–189.
9. T. M. Squires and S. R. Quake, Rev. Mod. Phys., 2005, 77, 977–1026.
10. D. Di Carlo, D. Irimia, R. G. Tompkins and M. Toner, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 18892–18897.
11. H. Amini, W. Lee and D. Di Carlo, Lab Chip, 2014, 14, 2739.
12. J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. Ebrahimi Warkiani and W. Li, Lab Chip, 2016, 16, 10–34.
13. B. Ho and L. Leal, J. Fluid Mech., 1974, 65, 365–400.
14. P. Ganatos, R. Pfeffer and S. Weinbaum, J. Fluid Mech., 1980, 99, 755–783.
15. M. E. Staben, A. Z. Zinchenko and R. H. Davis, Phys. Fluids, 2003, 15, 1711–1733.
16. M. E. Staben and R. H. Davis, Int. J. Multiph. Flow, 2005, 31, 529–547.
17. A. Nikoubashman, C. N. Likos and G. Kahl, Soft Matter, 2013, 9, 2603–2613.
18. H. Stark and D. Ventzki, Europhys. Lett., 2002, 57, 60–66.
19. S. Grollau, N. L. Abbott and J. J. de Pablo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 011702.
20. I. C. Gârlea and B. M. Mulder, Soft Matter, 2015, 11, 608–614.
21. I. C. Gârlea, P. Mulder, J. Alvarado, O. Dammone, D. G. Aarts, M. P. Lettinga, G. H. Koenderink and B. M. Mulder, Nat. Commun., 2016, 7, 12112.
22. U. Tkalec and I. Muševic, Soft Matter, 2013, 9, 8140–8150.
23. S. Copar, M. Ravnik and S. Žumer, Crystals, 2021, 11, 956.
24. N. Sulaiman, D. Marenduzzo and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 041708.
25. M. Urbanski, C. G. Reyes, J. Noh, A. Sharma, Y. Geng, V. S. R. Jampani and J. P. F. Lagerwall, J. Phys.: Condens. Matter, 2017, 29, 133003.
26. O. Wiese, D. Marenduzzo and O. Henrich, Soft Matter, 2016, 12, 9223–9237.
27. S. Mondal, I. M. Griffiths, F. Charlet and A. Majumdar, Fluids, 2018, 3, 39.
28. Ž. Kos and M. Ravnik, Sci. Rep., 2020, 10, 1446.
29. K. Fedorowicz, R. Prosser and A. Sengupta, Soft Matter, 2023, 19, 7084–7092.
30. S. Mondal, A. Majumdar and I. M. Griffiths, J. Colloid Interface Sci., 2018, 528, 431–442.
31. T. Stieger, M. Schoen and M. G. Mazza, J. Chem. Phys., 2014, 140, 054905.
32. T. Stieger, S. Püschel-Schlotthauer, M. Schoen and M. G. Mazza, Mol. Phys., 2016, 114, 259–275.
33. M. Lesniewska, N. Mottram and O. Henrich, Soft Matter, 2022, 18, 6942–6953.
34. H. Híjar, Phys. Rev. E, 2020, 102, 062705.
35. P.-G. de Gennes and J. Prost, The Physics of Liquid Crystals, Oxford University Press, New York, 1993.
36. D. C. Wright and N. D. Mermin, Rev. Mod. Phys., 1989, 61, 385–432.
37. A. N. Beris and B. J. Edwards, Thermodynamics of flowing systems, Oxford Univ. Press, 1994.
38. M. Škarabot, M. Ravnik, S. Žumer, U. Tkalec, I. Poberaj, D. Babic, N. Osterman and I. Musevic, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 1–8.
39. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 1–18.
40. N. Q. Nguyen and A. J. Ladd, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 046708.
41. A. Ladd, J. Fluid Mech., 1994, 271, 285.
42. A. Ladd, J. Fluid Mech., 1994, 271, 311.
43. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, Switzerland, 2017.
44. O. Parodi, J. Phys., 1970, 31, 581–584.
45. Ludwig GitHub repository, https://github.com/ludwig-cf/ludwig.
46. J. C. Desplat, I. Pagonabarraga and P. Bladon, Comput. Phys. Commun., 2001, 134, 273–290.
47. R. Adhikari, K. Stratford, M. E. Cates and A. J. Wagner, Europhys. Lett., 2005, 71, 473–479.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm01297b

---