## Abstract

Orientational director effects in nematic liquid crystals with small Ericksen number are investigated. The director field is disturbed by a semi-infinite plate on *y*=0 and and moving with a constant velocity *U*. Strong anchoring conditions at the plate are assumed. The resulting equations are a system of nonlinear partial differential equations for a nematic in the one elastic constant approximation. These equations are reduced to a coupled set of nonlinear ordinary differential equations by a suitable transformation. No such transformation seems possible for the many elastic constant case. The resulting equations are solved using analytical methods and strict bounding solutions obtained. These strict analytical solutions are compared with Picard iterated solutions.

## 1. Introduction

The hydrodynamic behaviour of a uniaxial nematic liquid crystal depends on both a velocity vector field * v* and a unit director vector field with the spatial orientation defined by the twist and tilt angles. The director field is also considered to be non-polar, such that the states

*and are indistinguishable.*

**n**A nematic liquid crystal with small Ericksen number disturbed by a semi-infinite plate on *y*=0 and and moving with a constant velocity *U* is considered (see figure 1). We investigate the orientational director effects driven by the plate velocity when the flow field does not distort the director, which is valid for flow regimes where (see §3*a*). It was observed in some experiments that besides the polar angle of the molecules the azimuthal angle is also altered inside the director domain walls created by external forces (see Neto *et al*. 1984). The director can escape to the third direction developing an axial component to avoid some prohibitive penalties due to some deformation, e.g. splay when homeotropic anchoring is assumed (see de Gennes 1974). In this work, we include all degrees of freedom of the director rotation. Thus, the continuum theory leads to a pair of coupled nonlinear partial differential equations governing the distortional phenomena of the director field described by the twist and tilt angles, *ψ* and *ϕ*, respectively, which are here functions of the Cartesian coordinates *x* and *y* as well as the time *t*. Note that *ϕ* is the angle between the plane *z*=0 and the director field.

We find that this nonlinear system can be reduced to a nonlinear system of ordinary differential equations using a suitable change of variables, i.e. and . The family of curves represents director front disturbances as parabolas or cross-sections of parabolic cylinders propagating with a constant velocity to the right (positive *x* direction). Thus, the direction of the directors, i.e. the twist and tilt angles, diffuse throughout the nematic sample. The motion of the semi-infinite plate is assumed steady. These ordinary differential equations can be integrated to give a nonlinear system of integral equations. It is shown that the solution of this system lies between the solutions obtained by bounding a kind of generalized diffusion coefficient above and below. This nonlinear diffusion coefficient turns out to be a complicated function of the material properties of the liquid crystal and the tilt angle. In addition, this also motivates the selection of an average diffusion coefficient, which is used to obtain approximate solutions. These analytical solutions are compared with numerical solutions obtained using the method of successive approximations (Picard iteration method).

The plan of this paper is as follows. We begin with the Ericksen–Leslie partial differential equations in Cartesian coordinates (see §2). In §3*a*, the boundary-value problem is formulated and the nonlinear system of differential equations describing the orientational director effects is obtained. Approximate solutions for the twist and tilt angles are also investigated here using bounding integral techniques. These results for the twist and tilt angles are given in §3*b* for problems with various boundary conditions. Moreover, in §3*b*, the method of successive approximations (Picard iteration method) is also considered and solutions are compared. The equations in Cartesian coordinates, governing the hydrodynamic behaviour of nematic liquid crystals, are reformulated in terms of the twist and tilt angles in appendix A. Finally, it is shown in appendix B that the director correction term due to the nematic flow is of the order of the Ericksen number.

## 2. Governing partial differential equations

The equations describing the hydrodynamic behaviour of a nematic liquid crystal are well known (Ericksen 1961; Leslie 1968; Stewart 2004). For completeness, these equations are presented here in Cartesian coordinates. The flow of a nematic liquid crystal is described by the velocity field * v* together with a director field

*, which is considered in this work as a non-polar unit vector giving the orientation of the anisotropic axis in these transversely isotropic liquids. With the assumption of incompressibility, the equations are the constraints(2.1)together with the balance laws:*

**n**the linear momentum equations(2.2)

the angular momentum equations(2.3)

and the constitutive relations: these are

the equations for the stress tensor (2.4)

the equations for the tensor (2.5)

the equations for the vector (2.6)

where the energy function *W* is defined by(2.7)with the viscous stress tensor given by(2.8)and the vector given by(2.9)

The rate of strain tensor is defined by(2.10)and the vector * N* is given by(2.11)where the vorticity tensor is defined by(2.12)

The material time derivatives of the director and velocity fields are given by(2.13)and the rotational and torsional viscosities are defined by(2.14)

The usual summation convention (summation over repeated indices) is used. A comma denotes partial differentiation and the superposed dot a material time derivative. The Cartesian coordinates are denoted by . Throughout the paper, unless otherwise stated, indices take the values 1, 2 and 3. In the conservation of linear momentum equations (2.2), *ρ* is the density, * F* any external body force and the stress tensor (asymmetric). In the conservation of angular momentum equations (2.3),

*σ*is a constant inertial coefficient, a generalized body force arising from a body couple present due to magnetic or electric fields,

*is a generalized intrinsic body force that depends upon the director vector field*

**g***through equation (2.6), and is a generalized stress tensor. The inertial term in these equations, , is usually assumed negligible and hence omitted in most treatments. In the constitutive relations (2.4), the scalar*

**n***p*is the arbitrary hydrostatic pressure following from the assumption of incompressibility, while

*γ*in (2.6) and the vector in (2.5) represent similar indeterminacies due to the constraint on the additional kinematic variable . The vector

*in (2.11) is the co-rotational time flux of*

**N***, which describes the internal motion of the director with respect to the fluid. In (2.10), A is the velocity gradient tensor (symmetric), i.e. the rate of strain tensor, whereas in (2.12) is the vorticity tensor (antisymmetric). In (2.8), is the viscous stress tensor, which is an even function of*

**n***. The vector in (2.9) is the hydrodynamic part of*

**n***and define the viscous torque on the director with coefficients of friction given by and , which are defined by (2.14) and usually called rotational and torsional viscosities, respectively. In (2.4), is the well-known Kronecker symbol and in (2.7), is the Levi-Civita symbol. The elastic torque on the director field is given by , where is the molecular vector field defined by(2.15)where*

**g***W*is the energy function and is the usual elastic coefficient related by equation (2.7), whereas in (2.8) is the viscosity coefficient. The elastic and viscosity coefficients are assumed constant if thermal gradient effects are ignored. Various authors have discussed the possibility of restrictions on these coefficients for physically acceptable behaviour. In particular, the coefficients in (2.8) are to be consistent with the thermodynamic inequality(2.16)

Other restrictions have been suggested, e.g. the Parodi relationship (see Parodi 1970) among the viscosity coefficients(2.17)

If in (2.3) we neglect the director inertial terms we obtain the equation(2.18)which together with the equations (2.4)–(2.14) as well as (2.1) and (2.2), and appropriate boundary and initial conditions enable (in principle) the velocity and director fields to be determined.

A slightly different form of these equations can be deduced. We consider an incompressible isothermal nematic liquid crystal with external body forces given by(2.19)where is some scalar energy density function, e.g. can be the magnetic or electric energy (e.g. Atkinson & Pereira 2005*a*). Note that the case of no external body forces is described by .

Differentiating (2.4) we obtain(2.20)

Inserting (2.5) and (2.6) into the angular momentum equations (2.18) we obtain(2.21)

Assuming that the energy function depends only on , we have(2.22)

Using (2.21), (2.20) becomes(2.23)

Considering a non-polar and unit director field we can write(2.24)

Thus, from (2.2), (2.19), (2.22), (2.23) as well as (2.24) we obtain(2.25)

Thus, the angular and linear momentum equations can be written in the form of (2.21) and (2.25), respectively.

## 3. Nonlinear boundary-value problem

### (a) Governing equations. Analytical solutions

In flow regimes where the velocity field follows instantaneously the director field, the back flow effects are minimal and, consequently, neglected. Here, we are interested in the orientational effects of the director field when the viscous forces due to the velocity of the fluid do not distort the director field. In flow regimes with very small Ericksen number , we can neglect the response of the director to the flow field and just use the director field for (e.g. Atkinson & Pereira 2005*b*). This means that the flow field does not create any distortion to the director, but whenever a director field gradient occurs in the nematic sample, a flow perturbation can be induced. The Ericksen number is defined by the ratio of the viscous forces to the elastic forces, i.e. . In this equation, *U* is the characteristic velocity, *L* is the length, which is of the order of the thickness of the sample (10–100) μm down to 1 μm, *μ* the viscosity, which is of the order (10^{−3}–10^{−2}) Pa s and *K* a representative elastic constant which is of the order (10^{−12}–10^{−11}) N whereas is the stress. Under conditions satisfying this flow regime, one can show in a straightforward manner that the perturbation director term due to the nematic flow is of the order of the Ericksen number, i.e. , where is a small average velocity (see appendix B). Provided the Ericksen number is small , the perturbation of the director field due to the slow motion of the liquid crystal is small.

The angular momentum equations are thus given by (2.21) (see §2) for (e.g. Atkinson & Pereira 2005*b*). For the typical elastic and viscosity parameters, which often take values of order mentioned earlier, the condition amounts to . For nematic samples with thickness *d*=*L* of order 10^{−5} m this corresponds to reasonably high speeds of order millimetres per second. Note that thicker nematic samples require lower velocities to fall into the region of applicability of this flow regime.

We consider a nematic sample disturbed by a semi-infinite plate on *y*=0 and (see figure 1). Furthermore, the semi-infinite plate is moving with a constant velocity *U* in the *x* direction, which represents here the characteristic velocity. The director field has the form(3.1)where and are the general twist and tilt angles, respectively. As the order of magnitude of the elastic constants is nearly the same, we adopt here the one elastic constant approximation for the energy function (see appendix A). Even though the elastic constant can be two or three times larger than and , the use of one elastic constant approximation should give qualitatively the same results as using the actual values. In appendix A, a reformulation of the Ericksen–Leslie equations (see §2) is given in terms of the twist and tilt angles. Thus, from (2.9)–(2.13)_{1} (see §2), (A 1), (A 11), (A 14), (A 16) (see appendix A), (3.1) and assuming flow regimes with no viscous forces due to the velocity of the fluid as well as no external forces, we obtain the following nonlinear system of partial differential equations:(3.2)and(3.3)where and . Note that the formulation of the Ericksen–Leslie equations in the twist and tilt angles has a singularity at . As *ϕ* is the angle between the plane *z*=0 and the director field, solutions for can also be obtained using the method prescribed in this work. Similarly, in this last case the direction fails to be represented.

For steady motion, this nonlinear system of partial differential equations can be expressed in terms of the moving coordinate(3.4)

Also, these partial differential equations reduce to ordinary differential equations if one uses the following change of variables:(3.5)where(3.6)

The curves represent director front disturbances propagating with a constant velocity *U* to the right (positive *x* direction).

Note that for the many elastic constant energy function the resulting system of equations analogous to (3.2) and (3.3) can be deduced. However, for this system no simple substitution such as (3.5) seems possible.

We also assume a prior treatment of the semi-infinite plate in order to obtain a perfect alignment of the director field at the boundary , which means that on *y*=0 and we have(3.7)with and specified constants with and .

We expect that and , so the twist and tilt angles become constant as . This requires further director field conditions, i.e.(3.8)where again and are constants with and .

Using (3.1)–(3.6), the nonlinear system of partial differential equations reduce after some algebra to the following nonlinear system of ordinary differential equations:(3.9)and(3.10)

We can rewrite (3.9) as follows:(3.11)where(3.12)

Thus, we obtain the following expression after integration of (3.11):(3.13)where is an unknown real constant.

Using (3.13), the tilt angle equation (3.10) can be written as follows:(3.14)

We can make these ordinary differential equations dimensionless using the following relations:(3.15)

One can see that the parameter *a* is constrained by the small Ericksen number, i.e. .

Thus, we obtain the following nonlinear ordinary differential equations:(3.16)and(3.17)where the dimensionless constant and .

If we denote,(3.18)we can write (3.17) in the form(3.19)

We can also rewrite (3.19) as follows:(3.20)with given by(3.21)

Thus, using (3.18), (3.20) as well as (3.21), we obtain(3.22)where *B* is an unknown real constant and .

We integrate (3.16) and (3.22) for using the boundary conditions described by (3.7) and (3.8) replacing by *η* to evaluate the arbitrary real constants *A*, *B*, *C* and *D* (where *C* and *D* are the real constants appearing after integration of (3.16) and (3.22), respectively) giving(3.23)and(3.24)where the error function is defined by(3.25)

The unknown constants *A* and *B* can be written in terms of as follows:(3.26)

Note that (3.23) and (3.24) with the unknown constants given by (3.26) satisfy the boundary conditions for the twist and tilt angles given by (3.7) and (3.8) replacing by *η*. Moreover, the twist and tilt angle solutions are defined here as continuously prolongated functions at *η*=0.

Equations (3.23) and (3.24) are exact expressions even though we do not know the exact expression for . Setting equal to , (3.24) becomes a nonlinear integral equation for .

Using maximum and minimum principles for differential equations, lower and upper bounding solutions of the nonlinear ordinary differential equations (3.16) and (3.17) with boundary conditions given by (3.7) and (3.8) replacing by *η* can be obtained. We say that is a lower bounding solution of this boundary-value problem if is continuously differentiable at least until second order on and(3.27)with boundary conditions given in this case by(3.28)

We say that is an upper bounding solution of the boundary-value problem if is continuously differentiable at least until second order on and(3.29)

Similarly, the boundary conditions can be written as follows:(3.30)

The term on the right-hand side of the inequalities (3.27) and (3.29) is given by(3.31)

Note that, in general, the equal signs of the boundary conditions in (3.28) are replaced by less or equal signs, whereas in (3.30) the equal signs of the boundary conditions are replaced by greater or equal signs. Furthermore, the continuously increasing function is bounded since *ϕ* only takes values between and . Similarly, the twist and tilt angle solutions are also bounded since they have horizontal asymptotes defined by and and take values higher or equal to and , respectively ( and ). These solutions are also assumed to be continuously differentiable at least until second order on .

To obtain bounding solutions, we consider approximate boundary-value problems with the same boundary conditions given by (3.28) and (3.30) as follows:(3.32)and(3.33)where(3.34)

Note that the integrals in the denominators of (3.34) are equal to . Thus, the lower and upper bounding tilt angle solutions can be written as follows:(3.35)and(3.36)

In (3.35) as well as in (3.36), the complementary error function is given by .

If is a solution of (3.17) with boundary conditions given by (3.7) as well as (3.8) replacing by *η* and with , then we can show that(3.37)

The proof is given by contradiction. Assuming that for is false, so has a positive maximum in . If the positive maximum occurs at an interior point, say , then attains a positive maximum at , so and . As for all , then this last inequality gives the contradiction . Thus, the positive maximum does not occur in . The positive maximum does not also occur at *η*=0 since or with the most general boundary conditions for , we have . In addition, we have . The boundary conditions also impose that cannot be a non-zero constant. Note that at infinity we have as well as . Moreover, as the solution functions, the variable coefficient as well as , , and are bounded on every closed subinterval of , they satisfy the boundness requirement of these maximum and minimum principles. Similarly, one can show that . This completes the proof. Finally, strict lower and upper bounding solutions give strict inequalities, i.e. .

We can work in a slightly different form. From (3.22) and using the boundary conditions given by (3.7) as well as (3.8) replacing by , we define(3.38)such that(3.39)

Note that and . Furthermore, if is a non-negative real constant . Looking for a positive maximum of we need to consider(3.40)

Using equations (3.38) and (3.40) we obtain(3.41)so . Note that the non-negative function between square brackets in (3.41), say , does not depend on *η*, since it can be shown that . Thus, if is higher, given by (3.41) is also higher. Furthermore, has in general a positive maximum in , since which means that . From these results, one can see that increases with increasing of , i.e. if for all *η*, where and . Note that and, consequently, . Thus, for the strict upper bounding solution we consider in (3.39), whereas for the strict lower bounding solution we use in the same equation. Here, and denote the minimum and maximum values of *ϕ* in *J*, respectively. Using these last two results and (3.39), we define strict lower and upper bounding solutions similar to those given by (3.35) and (3.36).

Bounding solutions for the twist angle can be defined and written in a similar form(3.42)and(3.43)which satisfy the following boundary conditions:(3.44)

Note that on we have by definition strict lower and upper bounding solutions for the twist and tilt angles. Moreover, the quotient of the integrals, say , in (3.42) as well as (3.43) define a non-negative decreasing function since .

As the exact solution lies somewhere between the lower and upper bounding solutions, we consider approximate analytical solutions for the tilt and twist angles defined by the functions and , i.e.(3.45)and(3.46)where(3.47)

The twist and tilt angle solutions and satisfy the boundary conditions given by (3.7) and (3.8) replacing by *η*. These solutions are defined using the average diffusion coefficients given by (3.47).

### (b) Theoretical results. Analytical and Picard iterated solutions

To compare our analytical results, a numerical algorithm based on a combination of the method of successive approximations (Picard iteration method), numerical integration and polynomial interpolation is implemented. To reduce the round-off-errors, all calculations are made using double precision. Basically, the Picard iterations are defined by , where is the standard Picard iterated function of order given by (3.24) and (3.26) replacing *ϕ* by and denotes the associated Picard operator. The initial iterated function is defined by (3.24) and (3.26) replacing *ϕ* by or . Similarly, we define the Picard iterations for the twist angle as , where is calculated replacing *ϕ* by in (3.23) as well as (3.26) and the initial iterated function is evaluated using these last two equations but replacing *ϕ* by or . The Picard operator associated with the twist angle is denoted by .

We only present here some illustrative examples.

The lower and upper bounding solutions (, , and ) as well as , and the approximate solutions given by the Picard iteration method ( and ) for the twist and tilt angles are shown in figures 2 and 3, respectively. The remaining curves are some Picard iterated functions of lower order. Here, we consider that the anchoring conditions are given by and . Moreover, as homogeneous boundary conditions, i.e. , are considered. In figures 2 and 3, the starting iterated function is defined by . From these curves, one can see as usual the very fast convergence of the method of successive approximations after some iterations. These results show that and can approximate well the numerical solutions.

In figures 4 and 5, one can see similar curves to those presented before, but the anchoring and boundary conditions are now given by and and . Here, the initial iteration is defined using .

It is possible to construct increasing and decreasing recursive sequences of functions associated with the lower and upper bounds such that the Picard iteration solution of order *n*, , lies somewhere between these lower and upper bounding iterated functions of order *n*. The recursive lower and upper sequence terms of order *n*=1, and , are defined using (3.35) and (3.36), respectively. The terms of order *n*=2 ( and ) are defined using the same equations but now replacing *ϕ* by and in each equation such that the trigonometric functions and are maximized or minimized. Similarly, one can define iterated functions of higher order. The graph of figure 6 shows that using these recursive sequences the bounding solutions can be improved at each iteration. Anchoring and boundary conditions are considered here to have the form and . The initial iteration associated with is defined using . Note that and .

## 4. Concluding remarks

Orientational director effects in nematic liquid crystals with small Ericksen number are investigated. The director field is disturbed by a semi-infinite plate moving with a constant velocity to the right (*x* positive direction). The continuum theory leads to a pair of coupled nonlinear partial differential equations for a nematic in the one elastic constant approximation. These equations are reduced to a system of nonlinear ordinary differential equations by a suitable transformation. No such transformation seems possible for the many elastic constant case. The resulting equations are integrated giving a nonlinear system of integral equations.

Using analytical methods, strict bounding solutions for the twist and tilt angles are obtained. It is shown using the maximum principle for differential equations that the solution of this system lies somewhere between the lower and upper bounding solutions obtained by bounding a kind of nonlinear diffusion coefficient below and above. It also motivates the use of an average diffusion coefficient to obtain approximate analytical solutions for the twist and tilt angles. From the theoretical results, it is shown that these bounding and approximate solutions and compare well with the numerical solutions obtained using the method of successive approximations (Picard iteration method). Decreasing and increasing recursive sequences are proposed to improve the lower and upper bounding solutions. As the standard Picard iterated function at each iteration has to satisfy the boundary conditions, a sequence of unknown constants *A* and *B* is generated. This means that the isotone (increasing) or antitone (decreasing) behaviour of the Picard operators and in as well as for the case with *A* and *B* fixed during all recursive Picard iterative process can be altered in some boundary-value problems. Thus, the isotone and antitone behaviour of these operators is constrained by the boundary conditions and the unknown constants *A* and *B*.

Experiments using this type of moving plate advancing into a sample of a nematic liquid crystal can be used to obtain viscosity measurements. Thus, these results are relevant, since one can obtain the director profiles driven by the plate velocity as well as information about the rotational viscosity.

## Acknowledgements

The author P.J.S.P. thanks the Instituto Superior de Engenharia de Lisboa and Ministry of Education of Portugal for a scholarship.

## Footnotes

- Received June 15, 2005.
- Accepted December 5, 2005.

- © 2006 The Royal Society