A note on couette flow of nematic crystals according to the Ericksen–Leslie theory

In order to model the flow of nematic crystals, the theoretical framework according to Ericksen and Leslie is applied. The essentials of the theory are compiled and then specialized to Couette flow. The profiles for linear velocity and orientation angle will be computed and, in particular, we shall also study the rise in temperature due to viscous dissipation, which is frequently ignored by mechanicians. Analytical and numerical solutions for the fields are derived for different boundary conditions and will subsequently be discussed.


INTRODUCTION
This short note is geared toward a description of the flow of fluids with internal rotational degrees of freedom, specifically liquid crystals. The flow behavior of such materials can be described by the theory of Ericksen and Leslie ( [1] or [2] ). In a nutshell, they restrict themselves to the description of rigid rod-like molecules, the motion of which is captured by introducing a single vector quantity, the so-called director, . Its motion is determined by a more or less postulated balance for the so-called "intrinsic director force." A detailed overview of the original Ericksen-Leslie theory is given in [2]. The most popular application of the director theory is in context with electromagnetic fields, e.g., in order to analyze optical effects. The aim of this paper, however, is to analyze a classical problem of fluid mechanics with the Ericksen-Leslie theory and to provide solutions, which have not been presented in the literature yet.
In literature the standard boundary value problems of liquid crystals were solved analytically as well as numerically. For example the Hagen-Poiseuille flow problem in [3] and different types of shear flow in [2, pp. 176] or [4, pp. 264]. Furthermore, planar as well as cylindrical Couette flow was investigated in [5] in terms of a phase plane analysis. In [6] Couette flow as well as Hagen-Poiseuille flow are analyzed numerically but only average properties of material parameters, e.g., the viscosity variation with pressure drop, are reported. In short, the Couette flow is extensively studied in literature, cf. [1,2,[5][6][7][8] However, the fact that Couette flow is not isothermal due to internal dissipation is frequently ignored by the mechanics community, see for example [9] or [2], Section 4.2.2 on Balance Laws: "We shall always assume isothermal conditions and therefore ignore thermal effects." In contrast to these previous investigations our aim is to study the Couette flow of such fluids not only from the mechanical point of view, i.e., by calculating the fields of linear velocity and orientation angle of the director, but also from a thermodynamic one, namely by studying the temperature field developing during the flow. This will require us to look at the balances of momentum, intrinsic director force, and internal energy in combination. Furthermore, we discuss different boundary conditions and, in addition, we investigate analytical solutions for a special parameter choice. For brevity's sake our presentation will be limited to the absolute essentials.
Finally note that there are several restrictions on the applicability of the Ericksen-Leslie theory as shown in [10,11] and [12]. Since the Ericksen-Leslie theory makes use of a macroscopic director field it is not able to describe the phase transition from the nematic phase to isotropic phase and vice versa. Starting from the mesoscopic point of view it is shown in [10] that the Ericksen-Leslie theory can only describe materials in which the molecules form a time-independent global planar phase or the molecules are local totally aligned. Therefore, the thermal motion of the individual needle-shaped molecules is ignored. Consequently the contribution of the individual rotation to the internal energy is ignored in the Ericksen-Leslie theory.

THE RELEVANT BALANCE EQUATIONS
Mathematically speaking, our objective is to determine the fields of (a) linear velocity, ( , ), (b) the director field, ( , ), and (c) the temperature field, ( , ), in all points, , and at all times, , within a region of space,  through which incompressible nematic crystal matter of constant mass density, 0 , and constant mass density of the director, d 0 , (in units of kg ∕m) is flowing. The determination of these fields relies on field equations for the primary fields. The field equations are based on balance laws and need to be complemented by suitable constitutive relations later. After introducing the substantial time derivative d ∕d of fields the relevant balances read as follows (see [13, pg. 268 and 276] or [14, pp. 354]): • balance of mass, • balance of linear momentum, • balance of director force, d 0 • and the balance of internal energy, where body force, body couple, and volumetric heat supply have been neglected. denotes the (non-symmetric, second order) Cauchy stress tensor, and is the so-called director stress tensor (in units of kg/s 2 , second order). The symbol ":" denotes a double contraction, ∶ = in Cartesian coordinates. The vector is known as the "director production density" or "intrinsic director body force." It is an additional constitutive quantity beyond the requirements of the ordinary Boltzmann continuum with the unit N/m 2 and subject to the constraint [ where the Gibbsian cross is defined via ( ⊗ ) × = × . Finally is the specific internal energy, and is the heat flux. Note that there is a connection between Equation (3) and the general balance of spin, e.g., [15, pg. 12]. This has been detailed in [14, Sec. 10.1].

RELEVANT CONSTITUTIVE EQUATIONS
Before we can turn to the problem of one-dimensional stationary Taylor-Couette flow, we need constitutive equations for the (non-symmetric) stress tensor, , the director stress tensor, , the director force vector, , and for the heat flux, . The usual techniques of rational thermodynamics can be used to derive the following explicit but lengthy non-linear constitutive relations (see [13] or [14,Chapter 10] for details): where = − is the (objective) difference between the angular velocity of the director, , and the vorticity, is the symmetric part of the velocity gradient, = (∇ ⊗ ) ⊤ is the director gradient and the abbreviations ∇ = (⋅) as well as ∇ ⊤ = (⋅) , ⊗ are used to represent derivatives with respect to the director field and the director gradient, respectively. The symbols ′ and ′ refer to dissipative terms of the stress tensor and the director force vector, respectively. Note that there is no dissipative term present in the director stress tensor, .
Note that these constitutive equations are part of the well established Ericksen-Leslie theory used by theoreticians, e.g., [6,9], and experimentalists. [6,[16][17][18] Our paper intends only to solve particular problems and not to review the theory, because this is done in the references listed. [13,14] It should be mentioned that besides the six different shear viscosities, , the Ericksen-Leslie theory contains a multitude of elastic stiffnesses, , originally introduced by, [19] which shows very clearly that the liquid crystals are at the borderline between fluids and solids. They occur in the elastic part of the free energy density, 0 , and they are related to deformation energy during splaying, bending, and twisting of a gradually changing director field, similarly as a gradient of displacement gives rise to the well known elastic energy in linear elastic Hookean solids (see [2, pp. 16]). The additional coefficients 1 and 2 are expressed in terms of the viscosities in order to satisfy the constraint in Equation (5) identically, see [13]. The quantities and are arbitrary parameters, which are used to ensure additional constraints, e.g., the inextensibility of the director. They play a similar role as the pressure in an incompressible fluid, see also the discussion in [14, pp. 366].
Already note that the Navier-Stokes constitutive relation for the stress tensor is obtained if the elastic energy and all viscosities except for 4 vanish. Then, the ordinary shear viscosity is expressed as = 4∕2.
Using the Onsager reciprocal relations, it was found in [20] that only five of the six viscosities are independent and some of them are related by the so-called Parodi relation In [16] this relation has been verified experimentally.

DISCUSSION OF POSSIBLE BOUNDARY CONDITIONS
In order to analyze specific flow situations, boundary conditions are required. In general there are two types of boundary conditions: essential and natural boundary conditions. Usually, essential boundary (or Dirichlet boundary conditions) are requirements for the sought-for fields themselves and are obtained from experimental observations and experience. The natural boundary conditions, however, arise from so-called jump conditions which in turn are consequences of the balance laws, cf. [14,21].
For ordinary fluids it is customary, although by no means obvious, to employ no-slip boundary conditions for the velocity field at solid walls, i.e., the movement of the fluid coincides with the wall in contact, if viscosity is taken into account. The discussion of possible slip at the interface between fluid and solid has a long history and is still ongoing, see [22]. According to [23], the correct velocity boundary conditions depend strongly on the scale of the system. Tangential slip becomes important if microfluid systems are considered, see [24]. We follow [25], where a nicely written overview on the history of the no-slip condition is given and the conclusion is drawn that in many cases the no-slip condition for the velocity can be assumed to be valid, at least as a good approximation. To this end, note that in [22, pg. 1], it is stated that the possible effect of slip is "unnoticeable" in real experiments (of certain size), because the influence of convection is of a much larger magnitude.
A completely different situation arises if the question is asked as to how to impose the correct boundary conditions for the director, . Since the advent of the Ericksen-Leslie theory in the late 1950s, so-called "strong anchoring conditions" were employed: the director field sticks to the wall at a certain angle which depends on the properties of the liquid crystal and the solid surface and is not effected by the flow situation.
Although this condition has been widely used by theoreticians, e.g., [5,6,9,13] and in molecular dynamics simulations, [26] , other boundary conditions were proposed, see [8], where the director angle with a wall is chosen proportional to the velocity gradient. However, there is strong experimental evidence for the anchoring boundary conditions for moderate flow velocities, see [1, pg. 26] and [27]. According to Leslie, [1, pg. [26][27]: "[...] in the presence of flow the optic axis remains strongly anchored at suitably treated surfaces. Some motivation for this stems from the fact that a torque from a magnetic field apparently has little influence upon alignment immediately adjacent to a solid surface, although occasionally some doubts are expressed (e.g., Saupe [87])." Note that there are situations, e.g., where elastic, magnetic and anchoring effects compete, in which also "weak anchoring" can occur. In [28] such a condition is analyzed: the director is fixed at a solid planar boundary but is allowed to rotate in plane. For these weak anchoring conditions an anchoring energy is usually proposed, see [3,29].
The question arises as to how the boundary conditions are realized experimentally. According to [30] the director aligns with the "easy direction" of the surface if strong anchoring is present. The authors propose several ways to obtain an easy direction, e.g., if the surface of the boundary is formed by a single crystal where a crystallographic plane coincides with the surface the easy direction can be prescribed. Another way of achieving a specific strong anchoring is to treat a surface mechanically, i.e., by carefully rubbing it in one direction, or chemically by etching the surface with suitable detergents. If the easy direction is normal to the surface, which may be achieved by etching glass plates, this situation is referred to as homeotropic anchoring [2, pg. 47].

ONE-DIMENSIONAL ANALYSIS OF COUETTE FLOW
We shall now concentrate on stationary one-dimensional Taylor-Couette flow between two parallel infinite plates separated by a distance ℎ in -direction: The motion can be either displacement or force controlled (see [31]), i.e., the top plate is either subjected to a constant velocity 0 = ( = ℎ) or to a constant shear stress, 0 = ( = ℎ). In both cases the movement will be in the -direction but we will assume that all fields depend only on the -coordinate. Note that this assumption leads to vanishing substantial time derivatives and therefore no caloric equation of state is needed in equation (4). From the discussion of the boundary conditions it is clear that in both cases, i.e., the velocity and the force controlled case, the anchoring boundary condition for director field is employed to describe a realistic flow situation. We will assume that the director sticks to the top and bottom plate in fixed but not necessarily equal angles, t and b , respectively.
In order to specialize the balances of mass, linear momentum, director force, and the corresponding constitutive equations of the Ericksen-Leslie theory to this particular flow, we propose the following semi-inverse ansatz for the velocity, the director field and the temperature: The boundary conditions in the velocity controlled case can now be stated as where (9) 2 may be replaced by ( = ℎ) = 0 in the force controlled case. Bearing in mind that all fields depend only upon the vertical coordinate and that the velocity points in -direction, the balance of linear momentum (2) and the director balance (3) reduce to Then after insertion of the constitutive relations (6) as well as the ansatz, and after some rearrangements of the last equations in (10) (see, e.g., [13, pp. 277] or [2, pp. 180]) a coupled system of (non-linear) second order differential equations is obtained (see [6]): Note that the Parodi relation in Equation (7) was used. Furthermore, is defined with a minus sign in order to be a positive number, as 1 ≤ 0 is required, see [13]. In the velocity controlled case the normalization with the plate velocity, ref = 0 , is the natural choice.
In the force controlled case the system can be simplified by integrating Equation (10) 1 first, which reveals that the shear stress is constant, i.e., ( ) = 0 , in the whole domain. Expressing the shear stress in terms of the primary fields, as in Equation (11) where the velocity is now normalized with ref = 2 0 ℎ ∕ 4 = 0 ℎ ∕ based on a convenient normalization for ordinary fluids. If this result is inserted in Equation (11) 2 , the system decouples Therefore, the angle (̄) can be determined first by solving Equation (14). The velocity is then obtained via integration of Equation (13). Note that (10) 2 reduces to an equation for the pressure gradient which is of no interest here. In order to study the temperature development during stationary one-dimensional Couette flow we start from Equation (4) and insert the constitutive relations from Equation (6). Bearing in mind that due to the special assumptions all substantial time derivatives vanish this leads to Note that the internal friction is linear in the velocity gradient which is in contrast to the quadratic term in the Navier-Stokes case.

ANALYTICAL SOLUTION
In the following we will derive an analytical solution for the director angle, , in a rather special case and proceed to find a series solution for the velocity. Since Equations (13) and (14) are highly non-linear no closed form solution to the system has yet been found, to the best of the authors' knowledge. However, for some special parameters, one is able to find analytical solutions. The first and major simplification arises through the so-called one constant approximation, where all non vanishing elastic constants are assumed to be equal, cf. [5,8,30]. By doing so the elastic parameter is equal to one and the auxiliary function becomes a constant, ( ) = 1, so that Equation (14) simplifies to read d 2 d̄2 = (1 + cos 2 ) 2 sin 2 cos 2 + 2 sin 2 + 2 + 1 .
Even though a first integration of this equation is readily obtained, the second is not. Therefore, we proceed to fix more parameters. In order to get an idea of reasonable scales we investigate the material parameters for MBBA (p'methoxybenzylidene-pn-butylaniline) at 25 • C given in [6, Table 1] and obtain the following set of parameters: = −1.031, = 1.228, = 0.078, = 0.958, = −0.214.
It is therefore reasonable to employ the approximations ≈ −1, ≈ 0 and ≈ 1. In addition, we put = − 1 2 for simplicity even though this does not represent the material MBBA. Note that some of these assumptions were made by [5] and [32] for simplification of a stability analysis. Furthermore, we employ an angle shift, = π 2 − . Hence, the right hand side of Equation (17)

becomes constant and the solution is obtained via integration
Now that the angle is computed the velocity can be calculated. With the aforementioned assumptions Equation (13) reduces to 1 where ∈ [0, π 2 ) was assumed and are the Bernoulli numbers. The solution for the velocity is then obtained via integration where 2 F 1 ( , ; ; ) is the hypergeometric function, (⋅) ′ = d ∕d̄and the auxiliary quantities are given by In the case of homeotropic anchoring, i.e., b = t = π 2 , the solution simplifies to a series in terms of the (incomplete) Euler Beta function, B( ; , ),̄(̄) Although the solution is given by an infinite series the coefficients as well as the Beta function decrease rapidly as the index increases. Therefore, the series can be cut off after a few terms for moderate values of . If, however, becomes large the cut off index needs to be increased. Note that due to the assumption ∈ [0, π 2 ) and Equation (19) we are restricted to < 4π. This implies either a restriction on the material parameters or on the system parameter ref ℎ. 1 Alternatively one could use the expansion, cf. [33], .

F I G U R E 1 Mean values of the series terms of the linear velocity solution in the force controlled case (left) and difference between analytical
and numerical velocity solution (right) for the force controlled case and the parameter set below Equation (18) Regarding the temperature Equation (15) where 2 is an integration constant. In the limiting case when the heat conduction parallel to the directors is much larger than perpendicularly to them, i.e., ≪ 1 or ∥ ≫ ⟂ , the temperature equation is readily solved using Equation (20) wherēt =̄(̄= 1).

RESULTS
The quality of the analytical series solution as well as the numerical solution will be assessed first. Then numerical solutions of the director angle, velocity and temperature profiles for the Couette flow of the material MBBA [6] are presented. Several boundary conditions are discussed and comparisons to the effects of ordinary fluid flow are pointed out and explained. The numerical solutions are obtained with Mathematica. [34] The analytical solution for the velocity in the force controlled case with homeotropic anchoring boundary conditions in Equation (23) is given in terms of an infinite series. In the left inset of Figure 1 the integral of the -th term, is plotted against its index for different values of in a semi-logarithmic chart. The terms decrease monotonically but at different rates as the index increases for different values of . To obtain an error of the order of machine precision for (say) a large value = 3.5π, a cut-off index of approximately = 140 is necessary. Therefore, the analytical solution̄1 40 is used for the verification of the numerical results. In the right inset of Figure 1 the error of the numerical analysis, |̄n um −̄a na 140 |, is shown in a semi-logarithmic chart for the parameter set explained below Equation (18). This shows that the numerical method has approximately an error of order 10 −9 to 10 −6 . Note that the error of the angle solution, which is not shown here, is roughly of the same order of magnitude. The numerical results for the velocity and for the orientation angle are shown in Figure 2 for the force controlled case with homeotropic strong anchoring, t = b = π 2 . In the figure the dimensionless coefficient was varied between values 1 and 25. Note that only the analytical solution is restricted to < 4π. The velocities shown in the left inset of Figure 2 turn from a linear to an S-shape for increasing values of . This is due to the boundary conditions of vanishing velocity at the bottom plate and a constant shear stress at the top, i.e., a constant velocity gradient at the top plate, see Equation (11) 1 . An impact of the additional viscosities in the constitutive equations (6) is also visible. The velocity profiles do not reach the normalized value 1 at the top plate as in the Navier-Stokes case (dashed black line), because the same traction 0 applied to the top plate will result in smaller velocities the higher the effective viscosity. It is, however, impossible to compare the variation of directly to the classical Navier-Stokes case (dashed black line): Indeed we could say that small values of correspond to negligible new viscosities 1 , but it could also mean that the stiffness 11 is very large, or that the plate distance is very small, in which case we expect length scale effects. Given that we fix the material parameters, it could also mean that a different force at the top plate is applied. Regarding the director angle a symmetrical rotation with respect tō= 1 2 into a fixed inclined position can be observed, right inset in Figure 2. Naturally, the characteristic of an angle of orientation does not exist in the case of a Navier-Stokes fluid. Note that in the following presentations is used instead of in order to analyze the variation from the initial state of perpendicular directors that is present if the upper plate does not move, see [2].
Also note that, first, the Ericksen-Leslie theory introduces new viscosities, which in turn leads to non-linear velocity profiles for the Couette flow. In contrast to that, the Navier-Stokes theory yields a linear profile for all material values. Thus, in view of experiments, such a higher theory might be required. Furthermore, we consider only a one-dimensional setting and in higher dimensions the different viscosities may have further impact. Second, the director field can be subsequently used to analyze, e.g., light scattering.
Moreover, the results shown in this and all the following plots are for a real material, MBBA [6] , and the assumptions on the constants required for the analytical solution do not apply. The assumptions, e.g., = 0 and = − 1 2 were only used for the sake of finding an analytical solution.
In Figure 3 the velocity controlled case is depicted for the same variation of . The qualitative behavior of the director angle remains the same when compared to the force controlled case. The velocity, however, is now normalized to 1 at the top plate resulting in S-shaped curves, point-symmetric w.r.t.̄= 1 2 . For the same values of slightly larger angles occur when compared to the force controlled case.
In Figure 4 the temperature profiles stemming from a numerical solution of the heat conduction equation (15) are shown. To this end it was assumed that the top and bottom plate are kept at the same constant reservoir temperature res . For the generation of the plots use was made of heat conduction data provided in [6]  glance it is surprising that the increase in temperature is less pronounced for some choices for the coefficient , because we expect more dissipation due to the various other viscosities. This, however, is due to the boundary condition of a constant traction applied to the top plate. In the non-Newtonian case the same drag force 0 will lead to smaller top plate velocities and smaller velocity gradients, which essentially govern the dissipation (see Figure 2 on the left). This difference to the Navier-Stokes case can be seen from Equation (13): The velocity gradient depends upon the angle auxiliary function ( ) while it is constant for an ordinary fluid. In the right inset of Figure 4 the temperature profiles for the velocity controlled case are shown. In this scenario the temperature rise due to viscous dissipation is larger when compared to the Navier-Stokes case (black dashed line). It is noteworthy, that the temperature decreases as the driving parameter increases. An explanation might be the following. In the case of homeotropic anchoring and small values of the driving parameter the directors do not point in direction of the velocity but are rather perpendicular to it which results in larger velocity gradients, see Figure 3. Thus, the dissipation is larger The dashed line represents a Navier-Stokes fluid. For = 25 the directors are displayed at discrete heights and therefore the rise in temperature as well. As the driving parameter increases the director field aligns more and more with the velocity which reduces the flow resistance and causes a smaller temperature rise.
In Figure 5 the velocity and angle profiles are depicted for the force controlled case but with non-homeotropic anchoring conditions, b = 0 and t = π 4 . Note that the boundary condition t = π 4 is a rather theoretical one and may correspond to the transition between perpendicular and tangential anchoring in a weak-anchoring scenario. For small values of the driving parameter and therefore values of smaller than or close to π 2 in the whole domain the liquid crystal moves faster when compared to an ordinary fluid (black dashed line). As the parameter increases the directors tilt beyond = π 2 in the middle of the channel. This results in lower velocities and higher gradients. However, the flow is still faster than in the Navier-Stokes case. The analysis of this scenario shows that a suitable surface preparation can reduce the flow resistance of the liquid crystal.
Similar comments as for the force controlled case apply to the velocity controlled case in Figure (6). A difference is that the angle profiles are more pronounced, see right inset of Figure (6). As before a drastic change in the velocity profile can also Finally, the temperature profiles for non-homeotropic anchoring with b = 0 and t = π 4 are shown in Figure 7 for the force and velocity controlled case. The temperature profiles also lose their symmetry w.r.t.̄= 1 2 . Additionally, the effect shown in Figure 4 is similar for moderate values of : The temperature rise increases with increasing in the force controlled case but decreases in the velocity controlled situation. However, the increase in the force controlled case in not monotonically, as the maximal temperature rise is larger for = 15 than for = 25. The reason for this can be seen from the velocity profiles in Figure 5 left inset: As a higher mode develops from = 5 to = 25 the velocity gradients for = 15 are larger than for = 25.
For the velocity controlled case in Figure 7, right inset, it is vice versa: From Figure 6 is it clear that the velocity gradient in the neighborhood of̄= 0.6 is maximal for = 25. Note that the dissipation in Equation (15) is not only determined by the velocity gradients but also by the angle of orientation. For = 25 much larger angles occur than for smaller values of . Therefore, the local dissipation and hence the temperature rise is larger in this region, which can be seen from Figure 7, right inset.

CONCLUSIONS AND OUTLOOK
In this short note we provided analytical as well as numerical solutions for the velocity profiles of the Couette flow problem based on the Ericksen-Leslie theory for nematic liquid crystals. We also calculated the angle of orientation for the nematic directors. In addition to that the temperature increase due to internal viscous friction was determined. We analyzed the effects of different boundary conditions, such as traction and velocity controlled upper moving plates as well as different director boundary conditions. The resulting effects were pointed out and explained.
In future works the solutions of the Couette flow problem based on the Ericksen-Leslie theory will be compared to solutions using other theories. In particular solutions based on the micropolar fluid theory from Eringen given in [35] will be considered.