NUMERICAL MODELLING OF LIQUID CRYSTAL FLOWS IN CONFINED DOMAINS by Somesh Bhatia B.Tech., National Institute of Technology Karnataka, Surathkal, 2016 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2019 © Somesh Bhatia, 2019 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled: Numerical Modelling of Liquid Crystal Flows in Confined Domains submitted by Somesh Bhatia in partial fulfillment of the requirements for the degree of Master of Applied Science in Mechanical Engineering Examining Committee: Dr. Dana Grecov, Mechanical Engineering Supervisor Dr. Ian Frigaard, Mechanical Engineering Supervisory Committee Member Dr. Fariborz Taghipour, Chemical and Biological Engineering Supervisory Committee Member iii Abstract Liquid crystals are anisotropic, viscoelastic materials with properties intermediate of solids and liquids. They are useful structural and functional materials and due to their ability to form ordered layers close to the bounding surfaces they are used as lubricants. Under the application of a hydrodynamic field, based on the type of velocity profile, values of non-dimensional numbers and anchoring angles, different orientation profiles are observed. Leslie-Ericksen and Landau-de Gennes theories are used to understand the evolution of the microstructure. Leslie-Ericksen theory, due to its simplicity can be used to obtain the behavior of flow aligning nematic liquid crystals. Landau-de Gennes theory is a mesoscopic model and apart from its ability to capture singular solutions can be employed in a study of lubrication using liquid crystals. This research work contains studies of liquid crystals in different flow conditions, such as the Couette flow and the Channel flow. In the study of Couette flow of Graphene oxide suspensions in water, the Leslie-Ericksen theory was used to obtain the orientation and viscosity profiles at different shear rates, up until flow alignment was observed. In the numerical study of pressure-driven channel flows, a Marker and Cell based solution methodology was implemented to solve the Leslie-Ericksen hydrodynamic theory. The expected flow alignment of liquid crystals was obtained which validated the solver. A preliminary study combining the moving wall and pressure driven flow showed that the orientation profile obtained depends on the local direction of shear and the direction of shear gradient. The scaling analysis was applied to Landau-de Gennes theory to derive the simplified equations of a planar lubrication theory. The theory was then validated by comparison with the solution obtained from numerically solving the full set of equations using the Couette flow profile. The parametric studies conducted showed that the solution was in the iv Elasticity-driven steady state. A discussion of the physical conditions showed its applicability for films of thickness lesser than 1 mm for a 1 m long domain. v Lay Summary Liquid crystalline phase is a phase of matter having properties between those of solids and liquids. Liquid crystals are fluids which have an order and this inherent order causes the presence of elasticity in the system. The interplay of the elasticity with the flow viscous effects leads to different orientations of the liquid crystals. The different orientations in the domain leads to a change in properties like viscosity. The Leslie-Ericksen theory is used to numerically study the flow of Graphene Oxide suspensions in a Couette flow. Further, the theory is numerically solved in the case of a 2-dimensional channel flow. The Landau-de Gennes theory is simplified and a lubrication theory is derived to numerically study the flow and understand the orientation in a thin domain. vi Preface The identification and design of the research project was performed by Dr. Dana Grecov and was conducted at the Industrial and Biological Multiphysics Laboratory, Department of Mechanical Engineering at the University of British Columbia, Vancouver. All the research described in the thesis, numerical implementation, and analysis was performed by me. vii Table of Contents Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ........................................................................................................................ vii List of Tables ................................................................................................................................ xi List of Figures .............................................................................................................................. xii List of Symbols ........................................................................................................................... xvi List of Abbreviations ............................................................................................................... xviii Glossary ...................................................................................................................................... xix Acknowledgements ......................................................................................................................xx Dedication ................................................................................................................................... xxi Chapter 1: Introduction ................................................................................................................1 1.1 Motivation ....................................................................................................................... 1 1.2 Objectives ....................................................................................................................... 3 1.3 Organization .................................................................................................................... 3 Chapter 2: Background Information ...........................................................................................5 2.1 Liquid Crystals ................................................................................................................ 5 2.2 Modeling of Liquid Crystal flows................................................................................... 8 2.2.1 Leslie-Ericksen theory ................................................................................................ 9 2.2.2 Landau-de Gennes theory ......................................................................................... 10 viii 2.3 Behavior of Liquid Crystals in External Fields ............................................................ 11 2.4 Lubrication Theory Applied to Liquid Crystals ............................................................ 13 2.5 Conclusion .................................................................................................................... 13 Chapter 3: Couette Flow of Graphene Oxide Suspensions in Water ......................................15 3.1 Governing Equations and the Leslie-Ericksen Theory ................................................. 15 3.2 Simplified Equations of Couette Flow .......................................................................... 18 3.3 Numerical Method ........................................................................................................ 20 3.4 Results and Discussion ................................................................................................. 22 3.4.1 Orientation Profile of Graphene Oxide ..................................................................... 23 3.4.1.1 Multistability and Multiplicity of Solutions ..................................................... 26 3.4.2 Viscosity Response of Graphene Oxide ................................................................... 28 3.4.3 Sensitivity Analysis of Frank’s Elasticity Coefficients ............................................ 31 3.5 Summary ....................................................................................................................... 33 Chapter 4: Leslie-Ericksen Theory Applied to Liquid Crystal Channel Flows.....................34 4.1 Governing Equations .................................................................................................... 34 4.1.1 Non-dimensionalized Equations for Channel Flow .................................................. 35 4.2 Numerical Implementation ........................................................................................... 37 4.2.1 Discretization Procedure ........................................................................................... 37 4.2.2 Boundary and Initial Conditions ............................................................................... 38 4.2.3 Solution Algorithm ................................................................................................... 38 ix 4.2.4 Time Stepping ........................................................................................................... 40 4.2.5 Computational Domain and Constants Used ............................................................ 41 4.3 Results and Discussion ................................................................................................. 41 4.3.1 Validation of the Solver ............................................................................................ 41 4.3.2 Parametric Studies on the Channel Flow .................................................................. 43 4.3.2.1 Effect of Changing the Ericksen Number ......................................................... 43 4.3.2.2 Effect of Changing the Anchoring Angle ......................................................... 44 4.3.3 Combination of Planar Couette and Channel Flow .................................................. 46 4.3.3.1 Study at Different Reynolds Number and Ericksen Number............................ 46 4.3.3.2 Study at Different Inlet Ratios .......................................................................... 47 4.4 Summary ....................................................................................................................... 50 Chapter 5: Lubrication theory applied to Landau-de Gennes theory ....................................52 5.1 Q-Tensor Theory ........................................................................................................... 53 5.2 Governing Equations and the Landau-de Gennes Theory ............................................ 53 5.2.1 Non-dimensional Full Equations .............................................................................. 57 5.2.2 Selection of Flow Scales for Lubrication.................................................................. 58 5.2.3 Simplified Lubrication Equations ............................................................................. 59 5.3 Numerical Methods ....................................................................................................... 60 5.3.1 Solution Procedure for the Full Equations ................................................................ 60 5.3.2 Solution Procedure for the Simplified Lubrication Equations .................................. 62 x 5.3.3 Simulation Parameters .............................................................................................. 62 5.4 Results and Discussion ................................................................................................. 63 5.4.1 Validation of the Solver ............................................................................................ 63 5.4.1.1 Validation of the Solution Obtained for the Full Equations ............................. 63 5.4.1.2 Validation of the Solution Obtained for the Simplified Equations ................... 64 5.4.2 Parametric Studies on the Simplified Equations ....................................................... 69 5.4.2.1 Effect of the Ericksen number on Orientation .................................................. 69 5.4.2.2 Effect of the Energy ratio on Orientation ......................................................... 70 5.4.3 Study to Determine the Physical Parametric Space .................................................. 71 5.5 Summary ....................................................................................................................... 73 Chapter 6: Conclusion .................................................................................................................74 6.1 Summary ....................................................................................................................... 74 6.2 Contributions................................................................................................................. 75 6.3 Limitations .................................................................................................................... 76 6.4 Future Work .................................................................................................................. 77 Bibliography .................................................................................................................................78 Appendices ....................................................................................................................................88 ............................................................................................................................... 88 ............................................................................................................................... 90 ............................................................................................................................... 98 xi List of Tables Table 3-1: Non dimensionalized variables for Couette flow ........................................................ 20 Table 3-2: Leslie viscosity coefficients of GO suspensions in water at different concentrations 21 Table 3-3: Comparison of values of theoretical and numerically obtained orientation angles at alignment....................................................................................................................................... 26 Table 3-4: Comparison of calculated and theoretical values of η1 at different concentrations ... 29 Table 3-5: Comparison of calculated and theoretical values of η2 at different concentrations ... 29 Table 4-1: Leslie Viscosity coefficients of MBBA ...................................................................... 41 xii List of Figures Figure 2-1: Different phases of liquid crystals ............................................................................... 7 Figure 2-2: Schematic representation of director n for (a) Calamitic (b) Discotic ......................... 9 Figure 2-3: Representation of Miesowicz viscosities corresponding to different orientations of liquid crystals in shear flow along different axes with the director is in the direction of: (a) velocity, (b) velocity gradient, and (c) vorticity .......................................................................................... 12 Figure 3-1: The different types of elastic orientations along with their associated elastic constants....................................................................................................................................................... 17 Figure 3-2: Couette flow geometry with R1 and R2 representing the radii of the inner and outer surfaces respectively. The inner surface is fixed, and the outer surface rotates at an angular velocity Ω. The enlarged image shows the flow direction along with the orientation θ of a Discotic liquid crystal molecule ............................................................................................................................ 19 Figure 3-3: Schematic showing the stable and unstable alignment angles for a nematic discotic liquid crystals based on the direction of the velocity and the velocity gradient ........................... 22 Figure 3-4: Orientation profile for 15 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates ............... 24 Figure 3-5: Orientation profile for 20 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates ............... 24 Figure 3-6: Orientation profile for 25 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates ............... 25 Figure 3-7: Orientation profile for 30 wt% GO across the dimensionless gap at different dimensionless shear rates showing that that flow alignment is obtained at high shear rates ........ 25 xiii Figure 3-8: Stable solution branch for 15 wt% GO showing the increase in non-dimensionless shear rate results in orientation reaching alignment value ............................................................ 27 Figure 3-9: Unstable solution branch for 15 wt% GO at different shear rates showing the increase in non-dimensionless shear rate results in oscillations within the orientation .............................. 27 Figure 3-10: Viscosity variation obtained at different concentrations for ϵ = 0.5 at anchoring 117 deg. The figure shows that the viscosity decreases with an increase in shear rate eventually reaching a constant value. ............................................................................................................. 30 Figure 3-11: Viscosity variation at different anchoring angles for ϵ = 0.5 at 15 wt% concentration. The figure shows that the value of viscosity at lower shear rates increases with the increase in anchoring angle while the value of viscosity at high shear rate remains unchanged. .................. 30 Figure 3-12: Viscosity variation in the sensitivity study of Frank coefficients. The figure shows that the change in ϵ causes no significant change in the overall viscosity profile ........................ 32 Figure 3-13: Variation in orientation near the boundary for different coefficients of elasticity. The figure shows that variation in orientation profile is minimal with a change in ϵ.......................... 33 Figure 4-1: Representation of Channel which is W units wide and L units long. The flow is in x direction and the inlet flow condition s represented with a parabolic velocity profile ................. 35 Figure 4-2: Representation of cell with staggered arrangement ................................................... 37 Figure 4-3: Comparison of velocity profiles obtained from analytical and numerical solutions. 42 Figure 4-4: Comparison of orientation profiles obtained from analytical and numerical solutions........................................................................................................................................................ 43 Figure 4-5: Variation in orientation obtained for different Er at Re=25 at anchoring angle 0. It is observed that increase in Er takes the orientation to the value of alignment angle ...................... 44 xiv Figure 4-6: Orientation profile obtained for positive anchoring angles at Re=13.2 and Er=10000. The value of θ in the domain reaches the flow alignment value .................................................. 45 Figure 4-7: Orientation profile obtained for the combination of moving wall and pressure driven flow for different parameters at anchoring angle 0. The figure shows that the flow alignment is observed in the lower half of the domain...................................................................................... 47 Figure 4-8: Velocity profiles obtained for different values of IR across the non-dimensional gap showing increasing convexity of the velocity profile with increase in IR .................................... 48 Figure 4-9: Non-dimensional shear rate obtained at different IR across the non-dimensional gap showing - (a) low IR causes the shear rate move towards a constant value and (b) higher IR has a point where the direction of shear rate reverses............................................................................ 49 Figure 4-10: Orientation profile obtained for different IR showing the flow alignment angle obtained in the lower half of the domain and for IR=0.5, the stable orientation angle is negative....................................................................................................................................................... 50 Figure 5-1: Representation of length scales of the system. The length scales along x and z is L and the scale along y is 𝛿L .................................................................................................................. 58 Figure 5-2: Planar Couette flow velocity profile obtained from solution of full equations for the case Re=0.0001, Er=1, R=0.1. The image shows the typical linear variation .............................. 64 Figure 5-3: Pressure profile obtained from solution of full equations for the Re=0.0001, Er=1, R=0.1 ............................................................................................................................................. 65 Figure 5-4: Comparison of solutions of full and simplified equations at H=1, L=10, Er=1, R=0.1, dudy=1. The orientation profiles are significantly different for different solutions ..................... 66 Figure 5-5: Comparison of solutions of full and simplified equations at H=1, L=10, Er=10, R=1, dudy=1. The orientation profiles are significantly different for different solutions ..................... 66 xv Figure 5-6: Comparison of solutions of full and simplified equations at H=0.1, L=10, Er=50, R=1, dudy=1. The orientation profiles obtained for both the solutions do not coincide ....................... 67 Figure 5-7: Comparison of solutions of full and simplified equations at H=0.1, L=10, Er=50, R=1, dudy=10. The orientation profiles obtained for both the solutions do not coincide ..................... 67 Figure 5-8: Comparison of solutions of full and simplified equations at H=0.01, L=10, Er=10, R=1, dudy=10. The orientation profiles overlap in this case. ............................................................... 68 Figure 5-9: Comparison of solutions of full and simplified equations at H=0.01, L=10, Er=10, R=1, dudy=100. The orientation profiles overlap in this case. .............................................................. 68 Figure 5-10: Variation of θ at R=0.001 with increase in Er showing an increase in deviation from the anchoring value ....................................................................................................................... 70 Figure 5-11: Variation of θ at Er=100000 for different R at unit shear rate and H=0.01. The orientation profile obtained shows that the value of orientation across the domain is constant until a transition occurs at R=100 ......................................................................................................... 71 Figure 5-12: Representation of the non-dimensional numbers for different film thickness showing that the present theory is applicable to δL=0.001 (L=1) ............................................................... 72 xvi List of Symbols 𝑐 Concentration of molecules 𝐷𝑒 Deborah number 𝐷𝑟 Pre-averaged rotational diffusivity 𝐸𝑟 Ericksen number 𝑘 Boltzmann constant 𝑘𝑖𝑖 Frank Elasticity coefficients 𝐿 Length 𝐿𝑖 Landau elasticity coefficients 𝑝 Pressure 𝐐 Orientation variable tensor 𝑟 Radial distance 𝑅 Energy ratio 𝑅𝑒 Reynolds number 𝑆𝑒𝑞 Equilibrium scalar order parameter 𝑡 Time 𝑇 Absolute Temperature 𝑢 Velocity in x direction 𝑈 Nematic potential 𝑣 Velocity in y direction 𝑤𝑓 Elastic energy xvii α𝑖 Leslie viscosity coefficients β Thermodynamic parameter 𝛾𝑖 Rotational viscosities ?̇? Shear rate 𝜂 Viscosity θ Orientation of the liquid crystals ρ Density 𝜏 Stress 𝜓 Potential function ω Angular velocity xviii List of Abbreviations GO Graphene Oxide IR Inlet ratio LC Liquid crystal LdG Landau-de Gennes LE Leslie-Ericksen xix Glossary Mesomorphic: Of intermediate form xx Acknowledgements First and foremost, I would like to thank my supervisor, Dr. Dana Grecov for her continued support and guidance. During the project there were instances when things were not going as planned. However, she remained patient and made me realize that research requires a lot of patience along with dedication. I am grateful for the financial support provided by Mitacs through the Globalink Graduate Fellowship. I would like to acknowledge the support that I was given by my parents and my sister. They have patiently listened to me when I was frustrated while constantly providing reassurance that things would eventually work out. I am grateful to my labmates Ali, Javad, Arian, Nick, Arash, Behzad, and Anish who have made working in the office more pleasurable and have prevented any frustration from creeping in. I would especially like to thank Nick for the multiple insightful conversations which helped me greatly with my project. Last, but not least are my friends who have been an invaluable part of my journey at UBC. Suraj, Vinamra, Pesh, my fellow Marauders, Pradeep, Shreya, Devadrita, Kiran, and Eric have always given me immense confidence and moral support. At UBC, Prashant, Ratul, Aishwarya, and Claire have positively impacted my physical and mental health. Harsh, Fagun, Jafar, Shahzaib, Miguel, Adriana, Faisal, Ketaki, Mohammad, Jasmeet, Nirmal, Sarah, Payel, and Juuso were always around and have been a part of the journey and made my program one to remember. I have learnt a lot from them and am thankful that they are a part of my life. xxi Dedication To my family, who has always believed in me. 1 Chapter 1: Introduction This chapter contains the motivation behind the research on the numerical modelling of nematic liquid crystals (LC) in confined flows. It also contains the objectives of the research and finally concludes with the organization of the thesis. 1.1 Motivation LCs have a variety of applications and this has triggered a lot of industrial and academic research to understand their behavior and there has been an interest to characterize their physical properties. There have been applications utilizing the flows of LCs. In the industry, LCs have been used in applications as a structural material. For example, they have been used in development of electronic parts made using injection molding [1]. Moreover, they are used as precursor elements in the manufacture of films, foams, blends and composites [2]. The 2-dimensional shear and pressure driven flow of LCs is seen in industrially important processes like blade coating and cavity filling [3]. Further, the application of liquid crystalline substances in functionalized microchannels is an area of interest for researchers in the growing field of microfluidics [4]. Apart from the examples listed above, LCs are also used as functional materials and their application as lubricants is considered in this thesis. Tribology gained importance through the 1960s as it became apparent that there is a need to study the field and understand its economic impact closely. Since then, research has taken place to not only understand the phenomena of friction and wear but also to develop effective lubricants [5]. Recently, research has been conducted to find naturally occurring environment-friendly lubricants. Research in tribology is conducted to reduce the impact of wear and losses of energy caused by friction. This would lead to better plant efficiency and consequently lead to savings [6]. Apart from industrial use there is 2 increasing interest in studying the lubrication that occurs in the joints of the human body [7]–[9]. Lubricants have typically been studied with certain additives and the tribological performance was studied. LCs have also been used as additives and also individually as lubricants [9]–[14]. There has been an increasing interest in carbon based LCs as they can be used in the synthesis of high-performance carbon materials [15]. Liquid crystalline processing has been used for carbon fiber spinning which is highly oriented [16]. There has been evidence that Graphene Oxide (GO) suspensions form a LC at concentrations in the range 15-40 wt.% [15], [17], [18] and studies show that GO can form a LC even in organic solvents [19]. GO has strong mechanical properties and can be mass produced from graphite [20]. There is a growing interest in understanding the properties of GO and recently its chemical properties and rheological properties [21]–[23]. The applications involving the use of liquid crystallinity of GO are increasing, for example in the case of development of graphene based nanofiltration membranes [24]. The tribological properties of GO as an additive in lubricants have been studied recently indicating its applicability as a lubricant [25], [26]. Further, GO is bio-compatible and its application as a lubricant in joint implants has been studied [27]. The procedure of obtaining the flow solution of LC involves solving the coupled system of Navier-Stokes equations and additional constitutive equations. Since the solution procedure can be complex, many studies have been conducted utilizing simplifying assumptions such as using an imposed velocity profile but not just limited to that. Numerical simulations give unique insights about the LCs in terms of their orientation and phase change transitions [28]. The literature studying the 3-dimensional evolution of LCs is scarce and it is important to develop numerical models to better understand the complex phenomena. The theories used to develop these models will be introduced in the forthcoming chapters. 3 1.2 Objectives The objectives of this thesis were to perform the following studies: 1. Study the Couette flow of GO LC suspensions in water and understand the behavior at different shear rates. In particular, the objective is to study the orientation and the viscosity profiles and understand the behavior at different concentrations and elasticity coefficients. 2. Develop a solver for 2-dimensional channel flow to obtain the flow solution of nematic LCs and study the behavior at different flow conditions. 3. Develop a lubrication theory by application of the Reynolds scaling approach on the LdG theory. Apply the developed theory to understand the effects of different parameters in thin confined flows of LCs. 1.3 Organization Chapter 2 consists of the Background Information pertaining to LC theory, modelling of LCs and a brief introduction of the theories used in the thesis. During the thesis work, two different theories were used to study the flows of LCs. Therefore, to allow for continuity in the flow of ideas, the thesis consists of chapters each of which contain all the aspects namely the Theory, Governing Equations, Numerical Methods involved, and Results and the subsequent discussion. Chapter 3 consists of the steps involved in the study of Couette flow of GO suspensions in water and the discussion of the results obtained. Chapter 4 contains the steps involved in the implementation of a 2-dimensional solver in a channel flow for nematic LCs using LE theory and the discussion of the results obtained. Chapter 5 consists of the steps involved in developing a Lubrication theory for nematic LCs based on the LdG theory and the requisite information of the techniques involved. It contains all the major results and the relevant discussion. Chapter 6 4 contains concluding comments about the presented research work, limitations, and directions for future work. 5 Chapter 2: Background Information The chapter contains the relevant information to understand the liquid crystalline phase and provide the literature review of the research developments made in the numerical modelling of LCs. It also contains the background information in the field of lubrication and its previous implementation in the field of LCs. 2.1 Liquid Crystals Liquid crystalline phase was discovered by Freidrich Reinitzer, an Austrian botanist and chemist in 1888 while he was studying the properties of cholesteryl benzoate [29]. Since then, research has been conducted to understand the liquid crystalline phase in greater detail. The origin of liquid crystallinity is molecular in nature [30]. Typically, it is observed in molecules which have a shape anisotropy. However, it is not the only factor that determines whether a system would have a liquid crystalline phase. Other factors such as molecular flexibility, electric charge distribution and polarizable anisotropy are important too. Liquid crystalline phase can more appropriately be termed as a mesomorphic phase. Mesophases as the word suggests, a phase of matter which has properties intermediate of solids and liquids. The mesophases have long range ordering of orientation but also can have long range ordering of the positional degrees of freedom. LC molecules can diffuse and viscous flow occurs like liquids [28]. Classification of Liquid Crystals: There are multiple ways of classifying LCs and are listed as follows [31]–[33]– 1. Based on the shape of the LC molecules: a. Calamitic LCs are materials which have a Rod-like shaped molecules. b. Discotic LCs are materials which have a Disk-like shaped molecules. 6 2. Based on the different cases of existence of order the mesophases: a. Nematic LCs are achiral substances which have a high degree of long-range orientational order but have no long-range translational order. Locally, the molecules are aligned in a preferred direction. b. Smectic LCs are layered structures containing a well-defined interlayer and have a translational periodicity with the period equal to the spacing between the layers. c. Cholesteric LCs are chiral materials which have a twist in their molecular orientation about an axis normal to the preferred direction. 3. Based on how the transition to liquid crystalline phase occurs: a. Thermotropic LCs are the materials where the phase transition occurs due to a change in temperature. Thermotropic LCs are stable only in a certain range of temperatures b. Lyotropic LCs which are multi component systems typically consisting of amphiphilic molecule and a polar solvent. The phase transition occurs due to the concentration of the amphiphilic molecule. Other factors such as molecule shape, size and charge affect the transition. In the present work, the studies are restricted to nematic LCs. Both thermotropic and lyotropic LCs are studied. Figure 2-1 represents the pictorial representation of the different types of LCs. 7 a. Calamitic Nematic b. Smectic A c. Smectic C d. Discotic Nematic e. Columnar f. Cholesteric Figure 2-1: Different phases of liquid crystals 0.5𝑃𝑜 8 2.2 Modeling of Liquid Crystal flows The models describing the behavior of LCs are of interest to researchers working with their flow related applications. However, studies also interest mathematicians and even recent attempts are made to understand the existence and properties of solutions 3-dimensional models of LCs [34]–[36]. Apart from the theoretical and experimental studies of LCs, simulations have been performed to understand their behavior. Numerical simulations give unique insights about the LCs in terms of their orientation and phase change transitions. Further, molecular level models can be used to develop links between micro and macro level properties [28]. However, as the molecular level approaches are not considered in the thesis they are not discussed in detail. As mentioned in Chapter 1, flows of LCs are seen in various industrial processes such as injection molding, blade coating and the coating of displays (LCDs). In order to ensure that the resulting part or coating is free of defects, the flow of LCs is modeled. In other cases, like the scattering devices developed using Smectic-A LC, the defect sizes need to be controlled. These defects are not necessarily endemic in nature and are affected by the boundaries and other external fields [37]. Extensive research has been performed on the theoretical and computational front to understand the morphology of LCs. The modeling of LCs has been performed using two different approaches. One is the continuum-based approach and the Leslie-Ericksen theory (LE theory) and the mesoscopic Landau-de Gennes theory (LdG theory) fall in this category. The other is the statistically based molecular level approach and Doi-Marucci-Greco theory is one of the prominent members from this category and it has been used for various case studies [38]–[40]. In this thesis, the continuum-based approach is chosen and both the LE and LdG theories have been implemented. A more detailed description about the development of these theories and 9 the related research is presented in the following sub-sections. Also, it is important to mention that the simulations in the field of LCs are not limited just to the development of models to predict the behavior under the application of external fields. Models have been developed to obtain the multiple material coefficients required and understand the phase transitions in LCs. An extensive review about these models is beyond the scope of this thesis and the reader is directed towards the following review articles [28], [41]. 2.2.1 Leslie-Ericksen theory The LE theory was proposed by Leslie [42] who derived the constitutive equations. Leslie proposed a vector called the director (n) which represents the average orientation of the LC molecules. Figure 2-2 represents the director along with the molecular orientations (which are represented by u) for both the Calamitic and Discotic LCs. Considering the symmetry of the director and taking it to be of unit magnitude, Leslie simplified the equations based on Ericksen’s theory of anisotropic fluids [43]. The constitutive equations were developed incorporating the Frank's theory to model the elasticity of the LCs [44]. The Frank elasticity coefficients were based on the idea that LCs have the tendency to resist and recover from distortion analogous to the tendency of elastic solids to resist strain [45]. LE theory u n nFigure 2-2: Schematic representation of director n for (a) Calamitic (b) Discotic a b u u u u u u u 10 is applicable for cases of slow-moving flows in which the LC orientation does not deviate by a large amount from the equilibrium position [46]. The shear flow behavior of nematic LCs depends on the magnitude of the reactive parameter λ which is the ratio of the flow-aligning effect of shear and the tumbling effect of vorticity. The value of the reactive parameter governs whether a LC is flow aligning (|𝜆| > 1) or tumbling in nature (|𝜆| < 1). The theory has been used to study both the flow-aligning LCs and the tumbling LCs to understand the behavior at different shear rates. Moreover, the LE theory is known to admit multiple solutions based on the anchoring angle (the orientation at the boundaries in case of a confined flows) and the flow conditions. In [47], the Poiseuille flow of Discotic nematic LCs was studied and the existence of multiple solutions based on the anchoring angle and the flow conditions was reported. The Couette flow of nematic LCs using LE theory for a fixed velocity profile has been previously implemented [11], [48]. Similar implementations employing a fixed velocity profile have been used to study Poiseuille flow [47], [49] and Jeffrey-Hamel flow [50] as well. The 2-dimensional channel flow has been studied for tumbling nematic LCs [51] and for nematic LC flows under the application of magnetic field by [52], [53]. Further, [54] studied the 3-dimensional pressure-driven flow under the effect of a magnetic field. Relatively recently, [55] studied the 2- dimensional channel flow taking the assumption of unidirectional flow. They explained the existence of multiple solutions and the energy requirements for the stability of the flow solutions. Further, shear and pressure driven flows were studied by [3] under the thin film approximation. 2.2.2 Landau-de Gennes theory LdG theory is relatively a recent development in the modeling of LCs. The initial work on the development of the theory occurred in the last decade of the 20th century [56], [57]. The theory 11 has since then been used to study nematic LCs and progress was made in modelling both the calamitic and discotic nematic LCs [58], [59]. The microstructure of the nematic LCs can be described based on the second order, symmetric, and traceless tensor order parameter Q [31]. The applicability of the LE theory is limited to sufficiently slow flows where the orientation dominates the rheology. However, the LdG theory can be used to model fast flows and/or description of short length scale phenomenon such as defect nucleation [59]. The theory can be used to describe both uniaxial and biaxial LCs. The theory is complete as it accounts for both the short-range and long-range order elastic effects. This allows the theory to model complex phenomena and it has led to the detection of banded textures, defect generation and coarsening phenomena [60]. Similar to the LE theory, the constitutive equations of the LdG theory have material parameters – Landau viscosity coefficients and Landau coefficient for elasticity. Since in the slow flow limit, LdG theory converges to LE theory, these material parameters can be obtained using the material parameters of LE theory [58], [59]. The theory has been applied to simulate simple shear flow and flow in eccentric cylinders of nematic LCs [9], [59], [61], [62]. Researchers have also used the LdG theory to simulate the flow of chiral nematic LCs [63]–[65]. 2.3 Behavior of Liquid Crystals in External Fields LCs are anisotropic in nature, and the anisotropy can be attributed to the presence of orientational order in the system. Besides imparting different optical properties, the anisotropy also induces a change in the physical properties such as the viscosity. Rheological studies on different LCs have shown that the viscosity of LCs changes on the application of an external field. The Miesowicz viscosities represent the viscosity of the LCs in a shear flow when the LCs are aligned 12 along the flow direction, along the velocity gradient direction and along the vorticity direction respectively. Figure 2-3 represents these different arrangements. Recent developments have been made to develop models to obtain these anisotropic viscosity coefficients based on the rheological measurements without the application of magnetic or electric fields [64]. The measurements of viscosity coefficients were based on the measurements conducted in a simple shear flow for different orientations by application of Electric/Magnetic fields. The application of such an external field leads to orientation in three perpendicular directions of the LCs. The technique was introduced by Miesowicz and thereby the viscosities measured are termed as Miesowicz viscosities [66]. The orientation of the nematic LCs changes under the application of any external field. This change is also observed on the application of a hydrodynamic field. Based on the applications considered as a motivation for this thesis, the focus is to study the behavior of flows of LCs. The stable planar solution branches of nematic LCs under the application of shear flow are the in-plane elastic state, in-plane tumbling-wagging state, in-plane viscous driven state, and in-plane wagging state [60]. There are stable out of plane mode solution branches as well [60]. However, in this thesis the focus is on in-plane elastic state and the in-plane viscous driven state and the other solution branches are not discussed here. The a b c Figure 2-3: Representation of Miesowicz viscosities corresponding to different orientations of liquid crystals in shear flow along different axes with the director is in the direction of: (a) velocity, (b) velocity gradient, and (c) vorticity 13 microstructure of the LCs is affected both by the viscous effects and the inherent elasticity. The elasticity is classified as long-range order and short-range order effects respectively. The isotropic-nematic transition is dependent on the short range order elasticity and the spatial orientation inhomogeneity is a result of the long range order elasticity [60]. 2.4 Lubrication Theory Applied to Liquid Crystals The development of lubrication theory or a thin-film theory is achieved by means of the scaling analysis. The analysis of fluid flow in narrow gaps was introduced by Reynolds and he used a scaling analysis to develop a theory which later came to be known as “Reynolds Lubrication Theory” and has been widely applied in the study of various thin film flows [67]. The scaling analysis has been applied to many other situations like porous media, filtration, adhesion and biological flows [68]. The basis of scaling analysis is that the derivatives in the direction of the film are significant compared to those in the other directions. The requirement for lubrication approximation to hold true is that the characteristic length across the flow should be much lesser than the characteristic length along the flow. As mentioned in Chapter 1, LCs have been used as lubricants and there is a research interest to study the thin films of LCs. There have been previous studies in thin film of LCs using the LE theory [69]–[72]. Moreover, recently the full set of equations of LdG theory were used to simulate the flow of LCs in the simplified hip-joint considering the nematic LCs as lubricant [9]. 2.5 Conclusion In this chapter, the background information regarding LCs was presented. This included the classification and the properties of LCs. Further, the developments in the modelling and 14 theoretical studies of LCs were covered. Finally, the lubrication theory and the developments of study of thin film of LCs relevant to lubrication were presented. 15 Chapter 3: Couette Flow of Graphene Oxide Suspensions in Water In this chapter, the LE theory is used to study the Couette flow of GO suspensions in water which is a Discotic LC. As described in Chapter 2, LE theory was one of the initial breakthroughs in the modelling of LC. The theory is primarily applicable to slow moving flows and the constitutive equations combine the effects of anisotropic viscous properties and elasticity present in the LCs by means of modification of the stress tensor in the Navier Stokes equations. The constitutive equations were derived by Leslie [42] and will be described in this chapter. The assumptions for the work described in this chapter are: 1. The flow is unidirectional 2. The flow is steady 3. The flow is incompressible 4. The flow is isothermal 5. Gravity is neglected 6. There is no external field applied to the system The chapter consists of the studies performed to understand the behavior of GO suspensions at different shear rates and the obtained orientation profile and the viscosity behavior are presented. The effect of concentration and elasticity on these profiles are studied. 3.1 Governing Equations and the Leslie-Ericksen Theory The LE theory consists of the equations (3.1) - (3.6) which govern the flow and development of the orientation of LCs. These equations include the equations for mass 16 conservation, linear momentum conservation, angular momentum conservation equations and the required constitutive equations. The mass conservation equation is as shown below- 𝜕𝜌𝜕𝑡+ ∇ ∙ (𝜌𝐮) = 0 (3. 1) The conservation of linear momentum equation is as shown below- 𝜌𝐷𝐮𝐷𝑡= −∇ ∙ ((𝑝 + 𝑤𝐹)𝐈) + ∇ ∙ 𝝉 + 𝐠 ∙ ∇𝐧 + 𝐆 ∙ ∇𝐧 + 𝜌𝐅 (3. 2) The angular momentum conservation equation is represented as shown below- ∇ ∙ (𝜕𝑤𝐹𝜕∇𝐧) −𝜕𝑤𝐹𝜕𝐧+ 𝐠 + 𝐆 = 𝜆𝐧 (3. 3) In equations (3.1) – (3.3), the terms are- • 𝐮: Velocity vector (𝑢, 𝑣, 0) • 𝐧: Unit vector representing the orientation • ρ: Density • 𝑝: Pressure • 𝑤𝐹: Elastic energy • 𝐠: Internal body force vector • 𝐆: Force vector due to an external body moment • 𝐅: External body force per unit mass • 𝝉: Viscous stress tensor • λ: Lagrange multiplier1 1 The Lagrange multiplier in the equation 3.3 occurs on simplifying the balance of elastic and viscous torques [53]. 17 As mentioned in section 3.1, the elastic effects of the LCs are modeled using the Frank’s theory and the term representing the elastic energy is as shown below- 𝑤𝐹 =12𝑘11(∇ ∙ 𝐧)2 +12𝑘22(𝐧 ∙ ∇ × 𝐧)2 +12𝑘33(‖𝐧 × ∇ × 𝐧‖)2 (3. 4) Here 𝑘11, 𝑘22, and 𝑘33 are known as Frank’s Elasticity coefficients for splay, twist, and bend respectively [44]. These orientations for the case of a discotic LC are as shown in Figure 3-1. In the present case, the LC orientation is affected only by the hydrodynamic field. The stress tensor due to viscous forces is as shown below- 𝛕 = α1(𝐧𝐧:𝐀)𝐧𝐧 + α2𝐧𝐍 + α3𝐍𝐧 + α4𝐀 + α5𝐧𝐧 ∙ 𝐀 + α6𝐀 ∙ 𝐧𝐧 (3. 5) The internal body force can be evaluated using the equation shown below- 𝐠 = −𝛾1𝐍− 𝛾2𝐀 ∙ 𝐧 (3. 6) In equation (3.5), 𝛼𝑖 represents the Leslie viscosity coefficients (𝑖 = 1. .6) and are unique to each liquid crystalline material. The Leslie coefficients are subject to certain thermodynamic constraints and these constraints can be found in any introductory material on LC [31], [42]. The equation (3.6) contains the rotational viscosities (𝛾𝑖) which are defined as 𝛾1 = 𝛼2 − 𝛼3 and 𝛾2 =𝛼5 − 𝛼6 respectively. In equations (3.5) and (3.6), the term 𝐍 = ?̇? − 𝐖 ∙ 𝐧 and it describes the Figure 3-1: The different types of elastic orientations along with their associated elastic constants 18 angular velocity of the director relative to the fluid velocity. 𝐖 = 𝟎. 𝟓 ∗ (∇𝐮 − (∇𝐮)𝑇) is the vorticity tensor. The rate of deformation tensor is defined as 𝐀 = 𝟎. 𝟓 ∗ (∇𝐮 + (∇𝐮)𝑇). Flow alignment is a property of the LCs which means that the LC molecules in the bulk of the domain orient at a fixed angle compared to the direction of the velocity and it is dependent on the Leslie viscosities and therefore is a material property [31]. The flow alignment angle can be calculated using θ𝑎𝑙 = 0.5 ∗ cos−1 (1𝜆). Here, λ is the reactive parameter and is calculated as 𝜆 =−𝛾2𝛾1. The value of the reactive parameter governs whether a given LC is flow aligning (|𝜆| > 1) or tumbling in nature (|𝜆| < 1). For discotic LCs, the value of λ < 0 . 3.2 Simplified Equations of Couette Flow The governing equations mentioned in section 3.1 are in the Cartesian coordinate system. However, for the case of Couette flow it is simpler to describe the equations in Polar coordinates and were derived by Leslie and Atkins [48]. As per convention, 𝑟 is the radial distance, and θ represents the angle of a given point from the horizontal in the anti-clockwise sense. Further, these equations were used by [11] to study different nematic LCs. They performed a comparison of the obtained viscosity with the results of rheological experiments. The Figure 3-2 represents the typical Couette flow geometry. The inner and outer surfaces are at radii R1 and R2 respectively. The outer surface is provided with an anti-clockwise angular velocity ω. Further, the orientation of a discotic LC molecule is shown along with the orientation angle. The equations (3.1) – (3.6) can be simplified by choosing the velocity vector as 𝐮: (0, 𝑟𝜔(𝑟), 0) and the director 𝐧: (sin(𝜃) , cos(𝜃) , 0). The resulting equations in the dimensionless form are as shown in equations (3.7) - (3.8). 19 𝑓(𝜃) [1?̃?𝜕𝜃𝜕?̃?+𝜕2𝜃𝜕?̃?2] +12𝑑𝑓(𝜃)𝑑𝜃[1?̃?2+ (𝜕𝜃𝜕?̃?)2] + (𝛾1̃ + 𝛾2̃ cos 2𝜃)?̃?2𝜕?̃?𝜕?̃?= 0 (3. 7) 1?̃?𝜕𝜕?̃?[?̃?3?̃?(𝜃)𝜕?̃?𝜕?̃?] = 0 (3. 8) The terms 𝑓(𝜃) and ?̃?(𝜃) are defined as shown in equations (3.9) and (3.10) respectively. 𝑓(𝜃) = 𝑐𝑜𝑠2 𝜃 + ε𝑠𝑖𝑛2 𝜃 (3. 9) 2?̃?(𝜃) = 2𝛼1̃ sin2 𝜃 cos2 𝜃 + (𝛼6̃ − 𝛼3̃) cos2 𝜃 + (𝛼5̃ − 𝛼2̃)𝑠𝑖𝑛2 𝜃 (3. 10) In the equations (3.7) – (3.10) the tilde represents that the parameters are non-dimensionalized. The parameter ε =𝑘33𝑘11⁄ and the terms are non-dimensionalized as per Table R2 R1 v xy θ Ω Figure 3-2: Couette flow geometry with R1 and R2 representing the radii of the inner and outer surfaces respectively. The inner surface is fixed, and the outer surface rotates at an angular velocity 𝛀. The enlarged image shows the flow direction along with the orientation 𝛉 of a Discotic liquid crystal molecule 20 3-1. In Table 3-1, 𝐿 represents the difference between the radii 𝑅1 and 𝑅2; η represents the average of the Miesowicz viscosities which are calculated as- 𝜂1 =12(α3 + α4 + α6) ; 𝜂2 =12(α4 + α5 − α2); 𝜂3 =12α4 (3. 11) Table 3-1: Non dimensionalized variables for Couette flow ?̃? ?̃? 𝛾?̃? 𝛼?̃? 𝜔𝐿2𝜂𝑘11 𝑟𝐿 𝛾𝑖𝜂 𝛼𝑖𝜂 3.3 Numerical Method In this section, the scheme for the development of the solver is described. The equations (3.7) and (3.8) are non-linear, non-homogeneous and second order uni-dimensional ODEs. They can be solved numerically using either shooting or relaxation methods. Since the LC system is prone to multiple solutions it is better to use the relaxation method to numerically solve these equations. Recently, this method was implemented by [11] to model the flow of nematic LCs between concentric cylinders and the solution was computed on a mesh of 500 elements representing a gap 10−4 𝑚. The mesh independency study was conducted by comparing the change in the orientation profile with mesh refinement. The relaxation method is used to obtain the solution of the boundary value problem of a system of 1st order differential equations. Therefore, the equations (3.7) and (3.8) are re-written into a system of 4 equations and the variables that are solved are 𝜃, 𝜔, 𝜕𝜃𝜕𝑟 and 𝜕𝜔𝜕𝑟. The consequence of this method is that the spatial derivative of 𝜃 and 𝜔 are calculated during the solution procedure. The methodology involves the 21 selection of an initial guess of the solution, based on which the solution branch gets selected. The finite-difference form of the equations is obtained to numerically solve the system of equations. The infinity norm of the change in the solution was selected as the convergence criterion and the value of tolerance was set as 10−6. Validation and constants used: The solution methodology has previously been implemented by [11] and the scheme was validated using nematic liquid crystalline materials DDA9 and AZA9. Further, there was a study performed comparing the obtained viscosity results with the rheological data of MBBA. Table 3-2: Leslie viscosity coefficients of GO suspensions in water at different concentrations Conc. (wt. %) 𝛼1(𝑃𝑎. 𝑠) 𝛼2 (𝑃𝑎. 𝑠) 𝛼3 (𝑃𝑎. 𝑠) 𝛼4 (𝑃𝑎. 𝑠) 𝛼5 (𝑃𝑎. 𝑠) 𝛼6 (𝑃𝑎. 𝑠) 15 -0.0138 0.0025 0.3611 0.0093 -5.75E-04 0.3631 20 -0.17651 0.007183 0.660202 0.012297 -0.00377 0.663616 25 -0.43438 0.011609 1.048167 0.01803 -0.00569 1.054085 30 -0.78526 0.016969 1.525367 0.027032 -0.00697 1.535361 In the governing equations (3.7)- (3.10) two sets of material coefficients namely - Leslie Viscosity and Frank’s Elasticity coefficients are present. The viscosity coefficients of GO were obtained using a model based on the rheological data of Graphene oxide suspensions in water (Private communication) and are as shown in Table 3-2. As shown, the coefficients are different for different concentrations of GO. The elastic coefficients of GO suspensions in water are not available in literature and therefore a sensitivity study was conducted to determine the impact of the elastic coefficients on the orientation of the LCs. The values of the elastic coefficients were 22 chosen to have the order 10 𝑝𝑁 which is the value typically observed for Frank elasticity coefficients of nematic LCs. The results of the sensitivity analysis are presented in section 3.4.3. 3.4 Results and Discussion This section presents the results obtained from the study of Couette flow of GO suspensions in water. Prior to the in-depth discussion of the results obtained, it is important to look at the stability diagram of the LC orientation for Discotic LCs which is shown in Figure 3-3 [59]. 𝜃𝑎𝑙 Stable Stable Unstable Unstable Figure 3-3: Schematic showing the stable and unstable alignment angles for a nematic discotic liquid crystals based on the direction of the velocity and the velocity gradient Velocity gradient Velocity 23 The figure represents the stable and unstable alignment angles of the director. The alignment angle 𝜃𝑎𝑙 is the angle that the director makes with respect to the primary direction of flow. All the nematic LCs are uniaxial and the director states n and -n cannot be distinguished [31]. Therefore, a 180° rotation results in the same orientation which explains both the stable and unstable alignment angles being 180° apart from each other. 3.4.1 Orientation Profile of Graphene Oxide This orientation of LCs depends on the balance of elastic and viscous forces. The elastic forces are dominant at lower shear rates and then as the shear rate increases, the viscous forces dominate the solution. In the numerical studies using LE theory, it is possible to determine if the solution is computed satisfactorily by observing the orientation in the bulk of the flow. Flow alignment of LCs is observed at higher shear rates and if the orientation in the bulk reaches the flow-alignment angle, then it can be concluded that the numerical solution is accurate. The orientation profiles shown in Figure 3-4 – Figure 3-7 were obtained by selecting ϵ = 0.5 at concentrations 15% - 30%. The anchoring angle in this case is 0.65π = 117°. The choice of anchoring angle is not arbitrary and was selected based on the study performed for viscosity and the discussion is provided in the upcoming sections. The initial guess of the orientation was selected as a value between the orientation angle and the anchoring angle. This value was set as a constant throughout the domain and it was closer to the anchoring angle at low shear rates and closer to the alignment value at higher shear rates. The initial guess of the velocity and the gradient was set as the typical Couette flow profile. From the Figure 3-4 – Figure 3-7, it is evident that the alignment value of orientation which is mentioned at the top of each figure is obtained for all the concentrations. The value of alignment angle can be calculated as θ𝑎𝑙 = 0.5 ∗ cos−1 (1𝜆). 24 Figure 3-4: Orientation profile for 15 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates Figure 3-5: Orientation profile for 20 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates 25 Figure 3-6: Orientation profile for 25 wt% GO across the dimensionless gap at different dimensionless shear rates showing that flow alignment is obtained at high shear rates Figure 3-7: Orientation profile for 30 wt% GO across the dimensionless gap at different dimensionless shear rates showing that that flow alignment is obtained at high shear rates 26 The Table 3-3 represents the comparison of alignment value obtained from the numerical solution with the theoretical alignment value of orientation angle for each concentration. The table confirms that the value of flow alignment is obtained for all the concentrations. Table 3-3: Comparison of values of theoretical and numerically obtained orientation angles at alignment Conc. (wt %) Theoretical angle (deg) Numerical angle (deg) Difference 15 94.756 94.792 -0.035 20 95.955 95.955 0.000 25 96.008 96.008 0.000 30 96.021 96.020 0.001 3.4.1.1 Multistability and Multiplicity of Solutions The LCs are prone to having multiple solutions of orientation and the solution branch selected is dependent on the initial guess given to the solver [11]. In this section, studies have been performed to show the stable and unstable branches of the solution. To highlight the difference between the solution branch selection the anchoring angle was chosen as 0.5π = 90°. The study was performed using ϵ = 0.5. To obtain the stable branch of the solution the initial guess was chosen above 90° and to get the unstable branch it was chosen lower than 90°. The stark contrast between the stable and unstable branches of the solution can be highlighted from the Figure 3-8 and Figure 3-9. The stable branch solutions are similar to the cases studied previously. In the case of the unstable branch, the solution develops oscillations and the orientation does not reach the value of alignment on the increase of shear rate. 27 Figure 3-8: Stable solution branch for 15 wt% GO showing the increase in non-dimensionless shear rate results in orientation reaching alignment value Figure 3-9: Unstable solution branch for 15 wt% GO at different shear rates showing the increase in non-dimensionless shear rate results in oscillations within the orientation 28 3.4.2 Viscosity Response of Graphene Oxide In section 3.4.1 the orientation profile was studied, and it was observed that orientation angle reaches the alignment value at higher shear rates. The next step is to examine the viscosity obtained at different shear rates. In the study, the apparent viscosity is calculated, and it is defined as the viscosity of a Newtonian fluid under the same conditions. The apparent viscosity can be evaluated using the equation (3.12) [48]. 𝜂 =𝑟2 − 𝑟1∫𝑑𝑠𝑔[𝜃(𝑠)]𝑟2𝑟1(3. 12) In equation (3.12), the function 𝑔(θ) is based on the definition shown in equation (3.10). However, to calculate 𝜂 with dimensions it is necessary to use the dimensional form of the equation (3.10). The equation (3.12) is based on a function for nematic LCs and therefore is applicable to a discotic LC such as GO. However, its applicability can be easily tested by ensuring that the Miesowicz viscosities calculated match with the values obtained using equation (3.11). Since the scheme is planar, only the values of η1 and η2 can be tested for. Table 3-4 and Table 3-5 represent the theoretical and calculated values of η1 and η2. To obtain the theoretical value, equation (3.11) was used and the value was compared with those obtained from equation (3.12). η1 was calculated using the value of orientation set to 90° in the entire domain and to calculate η2 the value was set to 180°. It is evident that the difference in values is minimal and thereby the use of equation (3.11) is justified. 29 Table 3-4: Comparison of calculated and theoretical values of 𝛈𝟏 at different concentrations Conc. (wt %) Theoretical η1(𝑃𝑎. 𝑠) Calculated η1(𝑃𝑎. 𝑠) Difference (𝑃𝑎. 𝑠) 15 3.668E-01 3.667E-01 -5.000E-05 20 6.681E-01 6.68E-01 4.237E-05 25 1.060E+00 1.06E+00 -4.118E-05 30 1.544E+00 1.5439E+00 2.032E-05 Table 3-5: Comparison of calculated and theoretical values of 𝛈𝟐 at different concentrations Conc. (wt %) Theoretical 𝜂2(𝑃𝑎. 𝑠) Calculated 𝜂2(𝑃𝑎. 𝑠) Difference (𝑃𝑎. 𝑠) 15 3.11E-03 3.10E-03 -1.25E-05 20 6.72E-04 6.72E-04 -2.68E-07 25 3.64E-04 3.66E-04 1.14E-06 30 1.54E-03 1.50E-03 -4.37E-05 Figure 3-10 represents the variation of viscosity for GO at an anchoring angle of 117°. It is evident that the increase in concentration results in a higher viscosity at flow alignment. The choice of anchoring angle may seem arbitrary, but it is selected based on trial and error and the final selection was made based on the viscosity at low shear rates. The objective was to find the viscosity value that would be close to the experimental value made available [private communication]. Further, as shown in Figure 3-11, different anchoring angles do not affect the value at flow alignment. However, at low shear rates, the anchoring angle affects the value of viscosity. 30 Figure 3-10: Viscosity variation obtained at different concentrations for 𝛜 = 𝟎. 𝟓 at anchoring 117 deg. The figure shows that the viscosity decreases with an increase in shear rate eventually reaching a constant value. Figure 3-11: Viscosity variation at different anchoring angles for 𝛜 = 𝟎. 𝟓 at 15 wt% concentration. The figure shows that the value of viscosity at lower shear rates increases with the increase in anchoring angle while the value of viscosity at high shear rate remains unchanged. 31 From Figure 3-10, it can be observed that as the concentration increases the viscosity increases and that GO is a shear thinning fluid. The shear thinning nature is a characteristic of flow aligning nematic LCs [1], [73] and also has been reported for GO previously [17], [21]. Both these observations can be explained based on the Miesowicz viscosities and the orientation angle in the domain. Revisiting Figure 3.2, it is understood that when the orientation angle is 90°, the value of viscosity is expected to be the least as it is equal to 𝜂1 and when the value of orientation angle increases towards 180°, the viscosity can expected to be the highest as it is equal to 𝜂2. From Table 3-3 and Table 3-5 it can be understood that value of 𝜂2 along with the value of alignment angle results in increasing values of viscosity when concentration increases. Similarly, the trend of higher anchoring angle resulting in a higher viscosity can be explained by the orientation in the domain becoming closer to the orientation representing 𝜂2. 3.4.3 Sensitivity Analysis of Frank’s Elasticity Coefficients As mentioned previously, there is no literary source available containing the elasticity coefficients of GO suspensions. To remedy this, a sensitivity analysis was performed using different ratios of ε =𝑘33𝑘11⁄ and this sub-section consists of the results obtained. The values that were tested are ϵ = 0.5, 1, and 2 which lie in the range of values that is typically expected in nematic LCs [50], [74], [75]. In this case GO of 15 wt.% concentration was studied, and the anchoring angle was set as 0.65𝜋 = 117°. It is important to note that the value of 𝐾22 doesn’t affect the simulation except for affecting the value of the average of the elastic coefficients and is maintained constant across the test cases. This value only affects the dimensional value of the Shear rate. Figure 3-12 represents the variation of viscosity with shear rate for the different cases of ϵ. It is evident that the changes in ϵ do not change the overall trend of variation in viscosity. 32 The next step is to examine the changes that the changing ratio causes in the orientation. It is known that at higher shear rates the effect of the elasticity is limited to a boundary layer while the orientation in the domain is primarily dominated by viscous forces [50]. Figure 3-13 shows the close-up of the orientation of the LC close to the boundary. The orientation profile is compared at different values of ϵ at anchoring angle 0.65𝜋 = 117° and at the same non-dimensional shear rate ω = 103. The difference in the profile suggests that the ratio indeed influences the transition from the anchoring angle to the alignment angle. However, Figure 3-13 which has data points collected at the same dimensionless numbers suggests that the different values of ϵ cause the same dimensionless shear rate to have different dimensional values. This was confirmed by dimensionalizing the shear rate values and it was observed that the variation was due to different dimensional values of shear rate. Figure 3-12: Viscosity variation in the sensitivity study of Frank coefficients. The figure shows that the change in 𝛜 causes no significant change in the overall viscosity profile 33 3.5 Summary To summarize, GO suspension in water which is a discotic nematic LC was the material that was studied in this work. The numerical simulation was performed for the Couette flow and the flow alignment behavior of LC molecules at high shear rates was observed. Consequently, the shear thinning behavior of the flow aligned LCs was observed. The effect of concentration and the anchoring angle on viscosity was explained and it was seen that the value of viscosity at flow alignment depends both on the value of alignment angle and η2. However, the effect of η2 is more dominant. Moreover, it was observed that the stability of the solution was dependent on the initial guess made to the solution which indicates the presence of multiple branches of solutions that the LC could take. A sensitivity analysis was performed to examine the effect of different ratios of elastic coefficients as the exact values of the coefficients are not available in literature. Figure 3-13: Variation in orientation near the boundary for different coefficients of elasticity. The figure shows that variation in orientation profile is minimal with a change in 𝛜 34 Chapter 4: Leslie-Ericksen Theory Applied to Liquid Crystal Channel Flows The LE theory was covered in detail in Chapter 3 and in this chapter the emphasis is on implementing a 2-dimensional solver for the channel flow of nematic LCs. The solver is then validated and used to perform parametric studies to determine its applicability at different flow conditions. The assumptions for the work described in this chapter are: 1. The flow is 2-dimensional 2. The flow is incompressible 3. The flow is isothermal 4. Gravity is neglected 5. There is no external field applied to the system 4.1 Governing Equations The governing equations of LE theory are as shown in equations (3.1)– (3.6). As mentioned in Chapter 3, the elastic effects of the LCs are modeled using the Frank’s elasticity coefficients and the equation is as shown in equation (3.4). In the 2-dimensional case of nematic LCs, the term due to the twist becomes equal to zero. Further the expression can be simplified by assuming 𝑘11 = 𝑘33 = 𝑘 [31], [53]. Then, the equation (3.4) takes the form shown below- 𝑤𝐹 = 12𝑘((∇𝐧): (∇𝐧)𝑇)2(4. 1) The constitutive equations of the theory are as shown in equations (3.5) and (3.6). 35 4.1.1 Non-dimensionalized Equations for Channel Flow The equations (3.1) – (3.6) and (4.1) are in the dimensional form. For the implementation of a numerical solution, the equations are first expressed in the non-dimensional form and are as shown in equations (4.2) - (4.8). Further, the director in this case is taken as 𝐧: (𝑐𝑜𝑠(𝜃) , 𝑠𝑖𝑛(𝜃) , 0) and the velocity profile is taken as 𝐯: (𝑢(𝑥, 𝑦, 𝑡) , 𝑣(𝑥, 𝑦, 𝑡) , 0). It is important to note that 𝜃 is measured in the anti-clockwise sense from the horizontal in this case. A representative image of the Channel flow geometry is as shown in Figure 4-1. The channel is of width 𝑊 and the length of the channel is 𝐿. The inlet and outlet of the channel is also marked. The other two boundaries are treated as walls. Since, the channel flow is considered it is convenient to express the equations in the Cartesian co-ordinate system. The change of variables was performed as follows and the methodology was adapted from [53] and is shown below- ?̃? =𝑥𝐿, ?̃? =𝑦𝐿, ?̃? =𝐯𝑈, ?̃? =𝑈𝐿𝑡, 𝑝 =𝑝𝜌𝑈2, ?̃? =𝛕𝜂, 𝑤?̃? =𝑤𝐹𝜌𝑈2, 𝐑 =𝐠𝜌𝑈2, 𝐆 =𝐆𝜌𝑈2(4. 2) Here, the terms 𝐿 and 𝑈 refer to the typical length and velocity scales of the flow. It should be noted that the length scale 𝐿 should not be misinterpreted as the geometry scale. 𝜂 = 0.5𝛼4 x y L W Parabolic inlet Figure 4-1: Representation of Channel which is W units wide and L units long. The flow is in x direction and the inlet flow condition s represented with a parabolic velocity profile Outlet 36 refers to the characteristic viscosity. In the present study, 𝐆 and 𝐅 are zero. Finally, the equations in the dimensionless form are as shown in equations below (with the tilde dropped for the ease of representation)- ∇ ∙ (𝐮) = 0 (4. 3) 𝑤𝐹 =12𝑅𝑒𝐸𝑟[(∂𝜃∂𝑥)2+ (∂𝜃∂𝑦)2] (4. 4) 𝐷𝐮𝐷𝑡= −∇ ∙ ((𝑝 + 𝑤𝐹)𝐈) + (1𝑅𝑒)∇ ∙ 𝐒 + 𝐑 ∙ ∇𝐧 (4. 5) 𝐷𝜃𝐷𝑡=1𝐸𝑟 𝛾1∇2𝜃 −12𝛾2𝛾1[(𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) cos(2𝜃) + (𝜕𝑣𝜕𝑦−𝜕𝑢𝜕𝑥) sin(2𝜃)] −12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥) (4. 6) In the above equations, the dimensionless numbers are the Reynolds number 𝑅𝑒 =ρ𝑈𝐿η and the Ericksen number 𝐸𝑟 =𝑈𝐿η𝐾. Note that the director handles the constraint 𝐧 ∙ 𝐧 = 1. The terms 𝐒 and 𝐑 on the RHS of equation (4.5) are obtained due to the constitutive equations. The stress tensor 𝐒 = 2𝐀 +𝚽 and the momentum equations in 𝑋 and 𝑌 can be written as shown below- 𝜕𝑢𝜕𝑡+ 𝑢𝜕𝑢𝜕𝑥+ 𝑣𝜕𝑢𝜕𝑦= −𝜕𝑝𝜕𝑥−𝜕𝑤𝐹𝜕𝑥+ 𝑅𝑥𝜕𝑛𝑥𝜕𝑥+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑥+ (1𝑅𝑒) (𝜕2𝑢𝜕𝑥2+𝜕2𝑢𝜕𝑦2) + (1𝑅𝑒) (𝜕2𝚽𝐱𝐱𝜕𝑥2+𝜕2𝚽𝐱𝐲𝜕𝑦2) (4. 7) 𝜕𝑣𝜕𝑡+ 𝑢𝜕𝑣𝜕𝑥+ 𝑣𝜕𝑣𝜕𝑦= −𝜕𝑝𝜕𝑦−𝜕𝑤𝐹𝜕𝑦+ 𝑅𝑥𝜕𝑛𝑥𝜕𝑦+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑦+ (1𝑅𝑒) (𝜕2𝑣𝜕𝑥2+𝜕2𝑣𝜕𝑦2) + (1𝑅𝑒) (𝜕2𝚽y𝐱𝜕𝑥2+𝜕2𝚽𝐲𝐲𝜕𝑦2) (4. 8) The terms - 𝑅𝑥𝜕𝑛𝑥𝜕𝑥, 𝑅𝑦𝜕𝑛𝑦𝜕𝑥, 𝑅𝑥𝜕𝑛𝑥𝜕𝑦, 𝑅𝑦𝜕𝑛𝑦𝜕𝑦 and the components of 𝚽 are listed in the Appendix A. 37 4.2 Numerical Implementation The objective was to implement a solver to numerically solve the equations (4.1) and (4.6)- (4.8). The variables that are being solved for are 𝑢, 𝑣, 𝑝, and θ. The solver was based on the work of [53] who worked on obtaining the flow solution of a 2D flow in the presence of magnetic fields. The scheme uses the GENSMAC method which was initially developed for free surface flows by [76]. GENSMAC method stands for Generalized Marker and Cell method which is an advancement of the Marker-and-Cell method which was introduced in the 1960s [77]. The solver was implemented on MATLAB to obtain the solution. This section describes the discretization, time stepping, and the type of grid used to implement the solver. 4.2.1 Discretization Procedure The governing equations were discretized using the Finite Difference Method. The spatial discretization was of second order. Following the work of [53] the staggered grid was chosen. Besides, the advantages of staggered arrangement are already well documented [78]. Therefore, 𝑢𝑖−12,𝑗 𝑢𝑖+12,𝑗 𝑣𝑖,𝑗−12 𝑣𝑖,𝑗+12 𝑝𝑖,𝑗 , 𝜃𝑖,𝑗 Figure 4-2: Representation of cell with staggered arrangement 38 the scalar variables 𝑝 and 𝜃 are evaluated at the center of the cells and the vector variables 𝑢 and 𝑣 are stored at the cell faces. The arrangement of variable storage can be seen in Figure 4-2. The convection term in the simulation of fluid flows needs special attention to ensure that the resulting solution methodology is stable and there is no numerical diffusion produced. In the present work, the convective terms are discretized using the third-order CUBISTA (Convergent and Universally Bounded Interpolation Scheme for the Treatment of Advection) scheme [79]. The scheme has previously been used in the simulation of Non-Newtonian fluids [80], [81] apart from being successfully implemented in the numerical simulation of nematic LCs [53]. 4.2.2 Boundary and Initial Conditions The initial condition was chosen to have zero velocity field and the orientation field was set to a constant value equal to the boundary value of orientation. The walls and inlet were set with Dirichlet boundary conditions for both the velocity and the orientation angle. Meanwhile, the outlet was chosen to have homogeneous Neumann boundary conditions for both the velocity and the orientation angle. At the walls, no-slip boundary condition was employed for the velocity and at the inlet parabolic velocity profile was chosen. 4.2.3 Solution Algorithm The solver is an explicit time stepping scheme and the methodology is adapted from [53] and was first introduced by [76]. Consider that the variables at a given time 𝑡𝑛 are already known and the values are to be obtained at time 𝑡𝑛+1 and the respective position 𝑥𝑘. The steps of the algorithm are repeated march the solution in time and are listed as follows- 39 Step 1: Calculation of intermediate velocity field (?̃? and ?̃?) at time 𝒕𝒏+𝟏 After discretization, the equations (4.9) and (4.10) can be solved using an explicit time stepping scheme. The values of the R.H.S of equations (4.9) and (4.10) are to be calculated using the actual velocity field (𝑢(𝑥𝑘, 𝑡𝑛) and 𝑣(𝑥𝑘, 𝑡𝑛)) and the orientation variable (θ(𝑥𝑘, 𝑡𝑛)) at time 𝑡𝑛. The boundary conditions are same as mentioned in section 4.2.2. It is to be noted that the velocity field thus obtained will have the correct vorticity at time 𝑡𝑛+1 [76]. 𝜕?̃?𝜕𝑡= − 𝑢𝜕𝑢𝜕𝑥− 𝑣𝜕𝑢𝜕𝑦−𝜕𝑤𝐹𝜕𝑥+𝑅𝑥𝜕𝑛𝑥𝜕𝑥+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑥+ (1𝑅𝑒) (𝜕2𝑢𝜕𝑥2+𝜕2𝑢𝜕𝑦2) + (1𝑅𝑒)(𝜕2𝚽𝐱𝐱𝜕𝑥2+𝜕2𝚽𝐱𝐲𝜕𝑦2) (4. 9) 𝜕?̃?𝜕𝑡= −𝑢𝜕𝑣𝜕𝑥− 𝑣𝜕𝑣𝜕𝑦−𝜕𝑤𝐹𝜕𝑦 +𝑅𝑥𝜕𝑛𝑥𝜕𝑦+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑦+ (1𝑅𝑒)(𝜕2𝑣𝜕𝑥2+𝜕2𝑣𝜕𝑦2) + (1𝑅𝑒) (𝜕2𝚽y𝐱𝜕𝑥2+𝜕2𝚽𝐲𝐲𝜕𝑦2) (4. 10) Step 2: Obtaining the potential function 𝛙 The Poisson equation shown below is solved to obtain the potential field at time 𝑡𝑛+1. It is to be noted that 𝜓 is a scalar and is evaluated at the cell centers. 𝜕2𝜓𝜕𝑥2+𝜕2𝜓𝜕𝑦2=𝜕?̃?𝜕𝑥+𝜕?̃?𝜕𝑦(4. 11) For the solution, homogeneous Neumann boundary conditions were implemented at the inflow and walls. Homogeneous Dirichlet boundary condition was implemented at the outflow [76]. Step 3: Obtaining the final velocity field The final velocity field at time 𝑡𝑛+1 can be obtained using- 40 𝑢 = ?̃? − 𝜕𝜓𝜕𝑥; 𝑣 = ?̃? − 𝜕𝜓𝜕𝑦 (4. 12) Step 4: Obtaining the orientation angle In this step the orientation θ at time 𝑡𝑛+1 is obtained by solving the discretized version of equation (4.6) using explicit time stepping. Step 5: Evaluation of Pressure and components of Non-Newtonian tensor The pressure field can be calculated using the equation shown below and the terms of the Non-Newtonian tensor are evaluated using the equations listed in Appendix A. 𝑝 =𝑑𝜓𝑑𝑡(4. 13) 4.2.4 Time Stepping It is known that the explicit time stepping results in divergence if the time-step selected is not small enough. According to [53], for the simulation of LC flows using the scheme the timestep can be selected as follows. There is a spatial restriction which occurs as a fluid particle cannot travel the distance larger than the grid size which results in the relations shown below- 𝑡𝑥 ≤∆𝑥|𝑢|𝑚𝑎𝑥; 𝑡𝑦 ≤∆𝑦|𝑣|𝑚𝑎𝑥(4. 14) Here, |𝑢|𝑚𝑎𝑥 and |𝑣|𝑚𝑎𝑥 represent the maximum modulus of velocities in the 𝑥 and 𝑦 direction respectively. ∆𝑥 and ∆𝑦 represent the grid sizes in 𝑥 and 𝑦. Further, there is a restriction due to the explicit discretization of the momentum equations and a von Neumann stability analysis of the linearized momentum equation shown below - 𝑡𝑣 =≤𝑅𝑒2(1∆𝑥2+1∆𝑦2)−1(4. 15) 41 The time step is then selected as the minimum out of 𝑡𝑥, 𝑡𝑦 and 𝑡𝑣. However, in the present work it was found that the solver did not always converge when the time step was selected in the aforementioned manner. Further literature review suggested the use of a multiplicative factor ε ∈(0,1] and this was incorporated in the solver [76]. 4.2.5 Computational Domain and Constants Used The numerical simulation was performed using MBBA as the LC and the material constants were obtained from [82]. The Leslie viscosities are as shown in Table 4-1. Table 4-1: Leslie Viscosity coefficients of MBBA In (𝑃𝑎. 𝑠) α1 α2 α3 α4 α5 α6 MBBA (near 25 ℃) -0.0181 -0.1104 -0.001104 0.0826 0.0779 -0.0336 The density of MBBA ρ = 1088𝑘𝑔/𝑚3. The domain was selected to have width 𝑊 = 0.001 𝑚 and 𝐿 = 0.02 𝑚 unless specified otherwise. In the numerical studies the convergence criterion was set as the infinity norm of the change in solution should be lower than 10−6. 4.3 Results and Discussion This subsection first provides the results of the validation. Followed by this, the results of the parametric studies are presented. Finally, the discussion with the combination of pressure driven and moving wall case is presented. 4.3.1 Validation of the Solver The solver was validated based on the analytical solution obtained from [53]. The analytical solution was obtained assuming fully developed flow solution in the 𝑥 direction and 42 setting the boundary orientation as θ = 0 at 𝑦 = 0,1. Further, the velocity in the 𝑦 direction, 𝑣 =0 in the entire domain. The equations of the variables 𝑢 and θ are shown in equations below- 𝑢(𝑦) = 𝑝,𝑥𝑦(𝑦 − 1)𝑅𝑒2(1 + 0.5(𝛼3 + 𝛼6))(4. 16) 𝜃(𝑦) = 0.5𝐸𝑟(𝛾1 + 𝛾2)𝑅𝑒2(1 + 0.5(𝛼3 + 𝛼6))𝑝,𝑥 (𝑦(𝑦 − 1)(1 − 2𝑦)6) (4. 17) Figure 4-3 and Figure 4-4 respectively represent the velocity profile and orientation profile comparing the solution obtained from analytical and numerical solutions. The figures show that there is a good agreement of the computed solution with the analytical solution. The other results such as the pressure contour and normal stresses were obtained similar to that shown in [53] and are not presented here. The mesh of size 32x640 was used in the parametric studies and the mesh independency study was performed by comparing the change in orientation profiles with mesh refinement. Figure 4-3: Comparison of velocity profiles obtained from analytical and numerical solutions. 43 4.3.2 Parametric Studies on the Channel Flow Once the solver was validated, parametric studies were conducted to explore the behavior of the solver at different values of the Ericksen number and anchoring angles. The interest was in obtaining solutions at higher values of the parameters to test the applicability of the solver at those values and to understand the behavior of LCs in the channel flow case. In this section, the fully developed orientation profiles for the parametric studies are presented. 4.3.2.1 Effect of Changing the Ericksen Number In this study, the effect of changing the Ericksen number on the flow solution is studied. The Reynolds number is set as 25. The study is conducted with the channel length 𝐿 = 0.04 𝑚 and mesh 32x1280. 𝐸𝑟 is the non-dimensional parameter that compares the viscous effect with the Figure 4-4: Comparison of orientation profiles obtained from analytical and numerical solutions. 44 elastic effect. A higher 𝐸𝑟 represents the dominance of the viscous effect over elasticity. In the case of flow aligning LCs, it is known that the LCs reach flow alignment when the viscous forces are dominant. From Figure 4-5, it is evident that the deviation in orientation profile increases with an increase in 𝐸𝑟. The value of θ increases until the value of alignment angle (5.7°) is obtained within the domain. This study provides the secondary validation of the flow solver. 4.3.2.2 Effect of Changing the Anchoring Angle The anchoring angle (the boundary condition for orientation at wall) in the present work has been set to the Dirichlet boundary condition. This type of anchoring is known as “strong anchoring”. In earlier studies of LCs, Robin boundary conditions have been applied and this is termed as “weak anchoring” [55], [72]. In [55], the two types of solutions are discussed – a strong flow in which the anchoring angles are set such that the director takes a full rotation in the domain Figure 4-5: Variation in orientation obtained for different Er at Re=25 at anchoring angle 0. It is observed that increase in Er takes the orientation to the value of alignment angle 45 from one wall to the other and the weak flow where the anchoring angles are equal in value. In the study reported in the present section, the anchoring angles are set to be equal and therefore according to [55], the weak flow solution is studied. The study is also motivated by the earlier finding reported for Poiseuille flow where the anchoring angle along with the Ericksen number was varied to find multiple branches of the solution [47]. In this study, positive values of anchoring angles were studied at 𝑅𝑒 = 13.2 and 𝐸𝑟 = 104. The choice of high values of the 𝐸𝑟-𝑅𝑒 pair was motivated as the interest was to study the behavior at flow alignment. It can be seen from Figure 4-6 that the anchoring angles higher than 0° took a different solution branch when compared to the cases where the anchoring angle was 0°. However, the value of orientation reached at flow alignment for anchoring angles higher than 0° is equal to 180° − 𝜃𝑎𝑙. Due to the symmetry of nematic LC, it is seen that these obtained orientations in the domain are the same values. However, since the study is not extensive like the one in [47], no other conclusion besides the existence of multiple solutions can be made. Figure 4-6: Orientation profile obtained for positive anchoring angles at Re=13.2 and Er=10000. The value of 𝛉 in the domain reaches the flow alignment value 46 4.3.3 Combination of Planar Couette and Channel Flow In the previous section, the solver was used to perform parametric studies using MBBA which is well characterized. In this study, the top wall is given a velocity and the flow profile turns into a combination of Couette flow and channel flow. However, the study is preliminary, and as per my knowledge there is no available literature for LCs in this area. To implement this in the solver, the moving wall velocity was set as 1 and an input velocity profile was prescribed. The inlet velocity boundary condition is the combination of the planar Couette flow velocity profile and a parabolic velocity profile similar to the previous cases. However, the parabolic velocity profile was scaled using a parameter - Inlet ratio (IR). IR is defined as the ratio of the magnitude of the maximum value of the input parabolic velocity profile w.r.t the velocity of the moving wall. 4.3.3.1 Study at Different Reynolds Number and Ericksen Number In this study, the effect of increasing the 𝑅𝑒 and 𝐸𝑟 is studied. The study was conducted at IR = 0.1 which indicates that the flow is primarily driven by the moving wall. Figure 4-7 shows the change in orientation obtained with a simultaneous increase in 𝑅𝑒 and 𝐸𝑟. The obtained orientation profile reaches the flow alignment value but only in the lower half of the domain. Further increase in the values of parameters suggested a regime change and the results are not presented. Therefore, for the given IR, the behavior of the LC in the domain when the viscous forces are dominant is as shown in Figure 4-7. However, as this is not the typically observed behavior of LCs further discussion is needed. In the following section, an attempt is made to explain these observations. 47 4.3.3.2 Study at Different Inlet Ratios The results from the studies in the previous sub-section lead to questions regarding the phenomenon that is taking place. Although, the difference in the orientation profile due to the combination of moving wall with a pressure driven channel flow can be anticipated, in this section an attempt is made to explain the findings by means of studying different IR at a fixed 𝑅𝑒 = 25 and 𝐸𝑟 = 1.03 ∗ 106. Once again, the existing knowledge is that in the case of: • A pressure driven channel flow, the direction of shear gradient switches at the center of the channel. The orientation profile is anti-symmetric about the midpoint [55]. Figure 4-7: Orientation profile obtained for the combination of moving wall and pressure driven flow for different parameters at anchoring angle 0. The figure shows that the flow alignment is observed in the lower half of the domain 48 • A planar Couette flow, a single value of orientation is obtained in the bulk of the domain at higher shear rates [11]. Firstly, the velocity profiles at different IR across the non-dimensional gap are as shown in Figure 4-8. It is observed that as the value of IR is increased, the velocity profile becomes convex in nature. The plot is supplemented by Figure 4-9 which shows the variation in non-dimensionless shear rate across the non-dimensional gap. The key observations are: • The non-dimensional shear rate at high IR shows a change in direction of shear direction • At low IR, the velocity profile resembles a Couette flow profile. Figure 4-8: Velocity profiles obtained for different values of IR across the non-dimensional gap showing increasing convexity of the velocity profile with increase in IR 49 Considering these observations and combining the information from Figure 4-10, it can be concluded that the LC orientation at flow alignment is dependent on the local shear rate. Especially, for the case of IR=0.5 the point at which the shear rate changes direction correlates to the point where the orientation reaches the same value as the anchoring angle. This observation holds value as the IR can essentially control the orientation at a certain location in the domain. Figure 4-9: Non-dimensional shear rate obtained at different IR across the non-dimensional gap showing - (a) low IR causes the shear rate move towards a constant value and (b) higher IR has a point where the direction of shear rate reverses 50 4.4 Summary In this chapter, the 2-dimensional solver for the channel flow of nematic LCs using the LE theory was implemented on MATLAB. The governing equations were discretized using the Finite Difference method and an explicit time-stepping scheme was implemented. The solver was validated using an analytical solution. The capability of the solver was studied by performing various parametric studies. On performing a study of the change in the orientation profile with the change in 𝐸𝑟 it was observed that the increase in 𝐸𝑟 causes the orientation profile progressing towards the flow alignment value. This observation shows the dominance of viscous effects over elastic effects at high values of 𝐸𝑟 and thereby providing a secondary validation of the solver. A study conducted at different anchoring angles showed the presence of multiple stable solution branches. Lastly, a preliminary study of the pressure driven flow with a moving wall was Figure 4-10: Orientation profile obtained for different IR showing the flow alignment angle obtained in the lower half of the domain and for IR=0.5, the stable orientation angle is negative 51 conducted. Through the study at different IR, it was observed that the orientation profile obtained is dependent on the local shear rate and the direction of shear gradient. Further, the study revealed that the value where the anchoring angle is obtained in the domain at steady state can be controlled using the value of IR. 52 Chapter 5: Lubrication theory applied to Landau-de Gennes theory In this research work, the application of LCs as lubricants has been a major motivation which led to an interest in the study of thin film flow of LCs. As mentioned in Chapter 2, the LdG theory is a complete model used to describe the flow of LCs. The theory incorporates the effects of short-range and long-range elasticity along with the viscous effects. The model is complex and has been previously solved by applying simplifying assumptions such as selection of simplified velocity profiles [59], [62], [83]. However, recently the full set of equations was implemented to simulate the flow of LCs in the arrangement of eccentric cylinders and also a 3-dimensional simulation between the geometry of a hip joint [9], [61]. The approach to studying thin films is well known and developed. The simplification is performed using the scaling analysis and the elimination of terms is based on their respective orders [68]. However, no such study has been conducted employing scaling analysis to the LdG theory. In this chapter, the scaling analysis is applied to the LdG theory to simplify the equations. In this chapter, the LdG theory is described first and then the procedure involved in the development of a Lubrication theory applicable to LCs is described. The developed equations are numerically solved for the case of planar Couette flow and the obtained results are compared with the simulation performed with the full theory. After the validation of the theory, parametric studies are conducted to understand the behavior of the LCs in the confined geometry. The major assumptions for the work described in this chapter are: 1. The flow is 2-dimensional 2. The flow is incompressible 3. The flow is isothermal 53 4. Gravity is neglected 5. There is no external magnetic or electric field applied to the system 5.1 Q-Tensor Theory The microstructure of the nematic LCs can be described based on the second order, symmetric, and traceless tensor order parameter Q [31] - 𝐐 = ∫(𝐮𝐮 −13𝐈) 𝑓 d2𝐮 (5. 1) Here 𝐮 is the unit vector showing the orientation of the molecules of the nematic LC (see Figure 2-2), I represents the unit tensor and 𝑓 is the orientation distribution function. Alternatively, 𝐐 can be represented based on the eigen vectors (𝐧, 𝐦, 𝐥) and eigen values (𝜆𝑛, 𝜆𝒎, 𝜆𝒍). 𝐐 = 𝜆𝑛𝐧𝐧 + 𝜆𝒎𝐦𝐦 + 𝜆𝒍𝐥𝐥 (5. 2) The eigenvalues can be defined as functions of the scalar order parameter 𝑆 and the biaxial order parameter 𝑃 [57] as shown- 𝜆𝑛 =23𝑆; 𝜆𝒎 = −13(𝑆 − 𝑃) ; 𝜆𝒍 = −13(𝑆 + 𝑃) (5. 3) These eigen values have the restriction 𝜆𝑛 + 𝜆𝒎 + 𝜆𝒍 = 𝟎. In the following sections 𝐮 will represent the velocity vector and not the unit vector of the LC molecules unless specified. 5.2 Governing Equations and the Landau-de Gennes Theory In this sub-section, the governing equations and the LdG theory is described. LdG theory is a tensorial theory and the Q-tensor described in Section 5.1 is used to describe the orientation variable. The governing equations of the theory follows from the dissipation function 𝚫 [59] – 𝚫 = 𝝉𝒔 ∶ 𝐀 + 𝑐𝑘𝑇𝐇 ∙ ?̂? (5. 4) 54 Here, 𝝉𝒔 is the symmetric viscoelastic stress tensor, 𝑐 represents the concentration of molecules in liquid crystalline state, 𝑘 is the Boltzmann constant, and 𝑇 is the absolute temperature. A is the rate of strain tensor, 𝐇 is the molecular field and ?̂? is the Jaumann derivative of the order parameter tensor. The tensors are as follows- 𝐀 =12(∇𝐮 + ∇𝐮𝑇) (5. 5) ?̂? =𝜕𝐐𝜕𝑡+ (𝐮 ∙ ∇)𝐐 −𝐖 ∙ 𝐐 + 𝐐 ∙ 𝐖 (5. 6) 𝑐𝑘𝑇𝐇 = −[𝛿𝐹𝛿𝐐][𝒔]= [𝜕𝑓𝜕𝐐 − ∇ ∙ (𝜕𝑓𝜕∇𝐐)][𝒔](5. 7) Here, 𝐖 is the vorticity tensor, 𝐹 is the total free energy, and 𝑓 is the free energy density. 𝐖 and 𝑓 are defined as - 𝐖 =12(∇𝐮 − ∇𝐮𝑇) (5. 8) 𝑓 = (𝑐𝑘𝑇) [12(1 −𝑈3) (𝐐 ∶ 𝐐) −13𝑈(𝐐 ∶ (𝐐 ∙ 𝐐)) + 14𝑈(𝐐 ∶ 𝐐)2+𝐿12𝑐𝑘𝑇{∇𝐐 ∶ (∇𝐐)𝑻} +𝐿22𝑐𝑘𝑇(∇ ∙ 𝐐) ∙ (∇ ∙ 𝐐)] (5. 9) Here 𝑈 = 3𝑇∗𝑇⁄ is the nematic potential where 𝑇∗ is the isotropic-nematic transition temperature, 𝐿𝑖 refer to the Landau coefficients and the superscript [s] means symmetric and traceless [59]. The molecular field 𝐇 can be calculated as- (𝑐𝑘𝑇)𝐇 = (𝑐𝑘𝑇) [(𝑈3− 1)𝐐 + 𝑈𝐐 ∙ 𝐐 − 𝑈(𝐐 ∶ 𝐐) ∙ (𝐐 +13𝐈)+𝐿1𝑐𝑘𝑇∇2𝐐 +𝐿22𝑐𝑘𝑇{∇(∇ ∙ 𝐐) + [∇(∇ ∙ 𝐐)]𝑻 −23(𝑡𝑟∇(∇ ∙ 𝐐))𝐈}] (5. 10) The above expression contains both short-range and long-range contributions. The dynamics of the tensor order parameter are given by – 55 ?̂? = 𝐅(∇𝐮,𝐐) + 𝐇𝒔𝒓(𝐷𝑅̅̅̅̅ , 𝐐) + 𝐇𝒍𝒓(𝐷𝑅̅̅̅̅ , ∇𝐐) (5. 11) Here, F represents the flow contribution, 𝐇𝒔𝒓 represents the short-range contribution, and 𝐇𝒍𝒓 represents the long-range contribution [83]. The expressions are – 𝐅 =23𝛽𝐀 + 𝛽 [𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐀 −23(𝐀 ∶ 𝐐)𝐈]−12𝛽{(𝐀 ∶ 𝐐)𝐐 + 𝐀 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐀 − [(𝐐 ∙ 𝐐) ∶ 𝐀]𝐈} (5. 12) 𝐇𝐬𝐫 = 6𝐷𝑟̅̅ ̅ [(𝑈3− 1)𝐐 + 𝑈𝐐 ∙ 𝐐 − 𝑈(𝐐 ∶ 𝐐) ∙ (𝐐 +13𝐈)] (5. 13) 𝐇𝐥𝐫 =𝐿1𝑐𝑘𝑇∇2𝐐 +𝐿22𝑐𝑘𝑇{∇(∇ ∙ 𝐐) + [∇(∇ ∙ 𝐐)]𝑻 −23(𝑡𝑟∇(∇ ∙ 𝐐))𝐈} (5. 14) Here, β is a thermodynamic parameter [59]. 𝐷𝑟̅̅ ̅ is the microstructural rotational diffusivity which is defined as shown in equation (5.15). [59] 𝐷𝑟̅̅ ̅(𝐐) = 𝐷𝑟 [1 −32(𝐐 ∶ 𝐐)]−2(5. 15) In equation (5.15), 𝐷𝑟 is the pre-averaged rotational diffusivity. It is to be noted that taking 𝐷𝑟̅̅ ̅ independent of Q results in spurious results for nematic LCs as the tumbling and wagging regimes are not predicted [57]. In the above equations, 𝐿𝑖 refer to the Landau coefficients and they are related to the Frank’s elasticity coefficients through the relations shown in equation (5.16). However, it is to be noted that these coefficients were obtained for the uniaxial phase [59]. 𝐿1 = 3 ∗ 𝑘22 − 𝑘11 + 𝑘336 ∗ 𝑆𝑒𝑞2 ; 𝐿2 = 𝑘22 − 𝑘11𝑆𝑒𝑞2 (5. 16) Here- 𝑘𝑖𝑖 refers to the Frank’s elasticity coefficients and 𝑆𝑒𝑞 is the equilibrium scalar order parameter which can be calculated as shown [83] - 56 𝑆𝑒𝑞 = 0.25 + 0.75√1 −83𝑈(5. 17) Apart from the equations describing the evolution of the orientation, the mass and momentum conservation equations are a part of the governing equations and are as shown - 𝜕𝜌𝜕𝑡+ ∇ ∙ (𝜌𝐮) = 0 (5. 18) 𝜌𝐷𝐮𝐷𝑡= ∇ ∙ 𝛔 + 𝜌𝐟 (5. 19) Here, • 𝐮 represents the velocity vector • 𝛔 represents the Total stress tensor • f represents the external body forces. The LdG theory has a set of constitutive equations that modify the stress tensor in the momentum equations. The stress tensor 𝛔 is defined as the sum of the symmetric component (𝝉𝒔), anti-symmetric component (𝝉𝒂) and Ericksen stress term (𝝉𝑬𝒓). The symmetric term 𝝉𝒔 = 𝝉𝒗 + 𝝉𝒆 has the viscous and elastic terms. The total stress tensor takes the form shown - 𝛔 = 𝝉𝒗 + 𝝉𝒆 + 𝝉𝒂 + 𝝉𝑬𝒓 (5. 20) Individual terms of equation (5.15) are listed in equations (5.21) - (5.24). 𝝉𝒗 = 𝜈1𝐀 + 𝜈2 [𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐀 −23(𝐀:𝐐)𝐈]+𝜈4{(𝐀:𝐐)𝐐 + 𝐀 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐀 + [(𝐐:𝐐)𝐀]𝐈} (5. 21) 𝝉𝒆 = 𝑐𝑘𝑇 {−2𝛽3𝐇 − 𝛽 [𝐇 ∙ 𝐐 + 𝐐 ∙ 𝐇 −23(𝐇:𝐐)𝐈] + 𝛽𝟐{(𝐇:𝐐)𝐐 + 𝐇 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐇 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐇 + [(𝐐:𝐐)𝐈]} } (5. 22) 57 𝝉𝒂 = 𝑐𝑘𝑇{𝐇 ∙ 𝐐 − 𝐐 ∙ 𝐇} (5. 23) 𝝉𝑬𝒓 = 𝑐𝑘𝑇{𝐿2(∇ ∙ 𝐐 ∙ (∇𝐐)𝑇) − 𝐿1(∇𝐐: (∇𝐐)𝑇)} (5. 24) In the above equations, the coefficients 𝜐𝑖 are termed as the Landau viscosities and can be related to the Leslie viscosity coefficients using the relations provided in [59]. 5.2.1 Non-dimensional Full Equations The non-dimensional equations of the full theory have been presented multiple times previously [9], [59], [61] and therefore the derivation is not included in this work. However, since the equations are numerically solved, they are shown below for the sake of completeness. Considering the flow to have a velocity scale 𝑈 and a length scale 𝐿 the non-dimensional form is as shown in equations (5.25) – (5.27). The evolution equation – ?̂? = 𝐅 + (𝟏𝐷𝑒)𝐇𝒔𝒓 + (𝟏𝐸𝑟)𝐇𝑙𝑟 (5. 25) In the equations, the following non-dimensional numbers are selected- Ericksen number: 𝐸𝑟 =𝑈𝐿𝑐𝐾𝑇6𝐿1𝐷𝑟 ; Deborah number: 𝐷𝑒 =𝑈6𝐿𝐷𝑟 ; Energy ratio: 𝑅 =𝐸𝑟𝐷𝑒 Reynolds number as 𝑅𝑒 =𝜌𝑈2𝑐𝐾𝑇 Here 𝜌 as the characteristic density apart from the previously introduced material constants. The momentum equations – 𝑅𝑒 {𝜕𝐮𝜕𝑡+ (𝐮 ∙ ∇)𝐮} = −∇𝑝 + 𝐷𝑒(∇ ∙ 𝝉𝒗) + (∇ ∙ 𝝉𝒆) + (∇ ∙ 𝝉𝒂) +1𝑅(∇ ∙ 𝝉𝑬𝒓) (5. 26) In the equations (5.26), the total stress tensor is non-dimensionalized using- 𝛔∗ = 𝛔 (𝑐𝐾𝑇)⁄ , 𝝂𝟏∗ =(𝝂𝟏6𝐷𝑟)(𝑐𝐾𝑇)⁄ , 𝝂𝟐∗ =(𝝂𝟐6𝐷𝑟)(𝑐𝐾𝑇)⁄ , 𝝂𝟒∗ =(𝝂𝟒6𝐷𝑟)(𝑐𝐾𝑇)⁄(5. 27) 58 5.2.2 Selection of Flow Scales for Lubrication Liquid crystalline phase is a multi-scale system and as the lubrication theory consists of simplification of equations based on the length scale the scales of the LCs are discussed. LCs have three length scales defined – the external (𝑙𝑒), the flow (𝑙𝑓), and the internal length scales (𝑙𝑖) respectively. The time-scales of the system are the external time scale (𝑡𝑒), the flow time scale (𝑡𝑓), and the internal time scale (𝑡𝑖) [62], [84]. The different length and time scales are shown below- 𝑙𝑒 = 𝐿 ; 𝑙𝑓 = √2𝐿1𝐷𝑟𝑐𝑘𝑇∗?̇? ; 𝑙𝑖 = √𝐿13𝑐𝑘𝑇∗ ; 𝑙𝑖 ≪ 𝑙𝑒 (5. 28) 𝑡𝑒 =𝑐𝑘𝑇∗𝐿22𝐷𝑟𝐿1 ; 𝑡𝑓 = 1?̇?; 𝑡𝑖 = 1𝐷𝑟; 𝑡𝑖 ≪ 𝑡𝑒 (5. 29) In the above equations, ?̇? is the shear rate. As shown in Figure 5.1, the 𝑥 axis is chosen to have an order of 𝐿 and as is the case usually in the lubrication one of the other axes has an order of δ𝐿 which in the present case is the 𝑦 axis. The primary assumption in the theory is that 𝛿 ≪ 1. In the simplified lubrication theory, δ𝐿 is chosen as the length scale while defining the non-dimensional numbers. x y z δ𝐿 𝐿 Figure 5-1: Representation of length scales of the system. The length scales along x and z is L and the scale along y is 𝛿L 59 5.2.3 Simplified Lubrication Equations This section consists of the simplified equations obtained after non-dimensionalization and performing scaling analysis. The simplification process is described in Appendix B. The resulting equations are as shown in equations (5.30) – (5.34). (−𝜕𝑝𝜕𝑥) +𝜕𝜏12𝐸𝜕𝑦+𝜕𝜏12𝑎𝜕𝑦+1𝑅𝜕𝜏12𝐸𝑟𝜕𝑦= 0 (5. 30) (−𝜕𝑝𝜕y) = 0 (5. 31) (1𝐸𝑟)(𝜕2𝜕𝑦2(𝑄11) −𝐿∗3𝜕2𝜕𝑦2(𝑄22)) +(1𝐷𝑒)((𝑈3− 1)𝑄11 + 𝑈(𝑄112 + 𝑄122) −𝑈(𝑄112 + 2𝑄122 + 𝑄222 + (𝑄11+𝑄22)2) (𝑄11 +13)) =−𝜕𝑢𝜕𝑦[1 − 3(𝑄112 + 𝑄122 +𝑄222 + 𝑄11𝑄22)]2[𝑄12 (1 +𝛽3) − 𝛽𝑄11𝑄12] (5. 32) (1𝐸𝑟)((1 +𝐿∗2)𝜕2𝜕𝑦2(𝑄12)) +(𝟏𝐷𝑒)((𝑈3− 1)𝑄12 + 𝑈(𝑄11𝑄12 +𝑄12𝑄22) −𝑈𝑄12(𝑄112 + 2𝑄122 +𝑄222 + (𝑄11+𝑄22)2)) =𝜕𝑢𝜕𝑦[1 − 3(𝑄112 +𝑄122 + 𝑄222 + 𝑄11𝑄22)]2[𝑄11 (1 − 𝛽2) − 𝑄22 (𝛽 + 12) +𝛽4𝑄112 +5𝛽4𝑄122 +𝛽4𝑄222 −𝛽3+𝛽4𝑄11𝑄22] (5. 33) 60 (1𝐸𝑟) ((1 +2𝐿∗3)𝜕2𝜕𝑦2(𝑄22)) +(𝟏𝐷𝑒)((𝑈3− 1)𝑄22 + 𝑈(𝑄112 + 𝑄122) −𝑈(𝑄112 + 2𝑄122 + 𝑄222 + (𝑄11+𝑄22)2) (𝑄22 +13)) =𝜕𝑢𝜕𝑦[1 − 3(𝑄112 + 𝑄122 + 𝑄222 + 𝑄11𝑄22)]2[𝑄12 (1 −𝛽3) + 𝛽𝑄12𝑄22] (5. 34) 5.3 Numerical Methods In this section the implementation of the solver is described. There are two aspects to the solution methodology- (1) solution of the full set of equations and (2) solution of the simplified equations from scaling analysis. The first stage in the implementation of the simplified theory is to only solve the simplified evolution equations using a fixed velocity profile. The complete methodology for each solution is provided in the respective sub-sections. 5.3.1 Solution Procedure for the Full Equations In order to validate the solutions obtained from lubrication theory they are compared with the results obtained from the full equations. In this section, the setup of the model to obtain the solution of the full set of equations on COMSOL Multiphysics® 5.3a is described. The system of equations (5.25) – (5.27) is nonlinear and is fully coupled. COMSOL is a commercial Finite Element package and it uses the Galerkin Finite Element Method to solve the system of PDEs. The equations shown were implemented using the combination of Laminar flow (spf module) and the General Form Partial Differential equation solver (General Form pde module). In this sub-section, the details of the model such as the flow geometry, Mesh, Boundary conditions, Initial conditions, and Solver are described. 61 Geometry and Meshing: The objective of the present study is to validate the lubrication theory developed in Cartesian coordinate system. Therefore, a rectangular domain is selected as it is the simplest choice. Along the similar lines, quadrilateral mesh was selected, and the mesh independency study was performed using the standard method by comparing the orientation profiles obtained. The final domain was 10 units in length and 1 unit wide. The final mesh consisted of 8000 elements and the element quality was uniformly 0.97 using Volume versus circumradius as the quality measure. Initial condition: The flow field was set to be at rest and the orientation tensor was set to random noise and it is said to represent a homogeneous orientation field with superposed infinitesimal fluctuations [60]. The initial field for the Q - tensor is as shown in equation 𝐐 = 𝑆𝑒𝑞 (𝐧𝐧 −13𝐈) ; 𝐧 = (cos(𝜀) , sin(𝜀) , 0) (5. 35) Where ϵ = (π/180)10−4 represents the noise and 𝑆𝑒𝑞 is the equilibrium scalar order parameter which can be calculated as shown in equation (5.17) Boundary conditions: In the flow module, the inlet and outlet were set to Periodic boundary conditions. The top wall was set to sliding wall with no-slip condition and was excited with a known velocity. The lower wall was stationary and no-slip conditions were used. The boundary condition in the pde module were set as Periodic boundary conditions at the inlet and outlet. The top wall and lower wall were given Dirichlet boundary conditions and the value was set to planar orientation. Solver implementation: The time-stepping was performed using the Backward Differentiation Formula (BDF) and the Multifrontal Massive Parallel sparse direct solver (MUMPS) was used to 62 invert the generated matrix. The tolerance of the solver was set as 10−6 and the solution was computed till steady state was reached. 5.3.2 Solution Procedure for the Simplified Lubrication Equations The simplified non-dimensional evolution equations of lubrication are as shown in (5.32) - (5.34). It is evident that the equations are stationary. Since these equations are simplified by the omission of multiple terms the results require validation. The validation can be performed by comparison for the case of a fixed velocity profile. In the study, Planar Couette flow was chosen as the velocity profile. The system of equations along with the boundary conditions is a Boundary Value Problem (BVP) and the solution of the set of equations can be obtained using the MATLAB solver bvp5c. bvp5c is a finite-difference implementation of the four-stage Lobatto Illa formula and it gives fifth-order accurate solutions. To solve the equations through bvp5c, they are re-written as a system of first order differential equations. The solver requires an initial guess and it was selected equal to the value of the orientation variable at the boundary. Meanwhile the derivative of the orientation variables in y was initialized as 0 throughout the domain. Dirichlet boundary conditions were implemented and the values are the same as those mentioned in section 5.3.1. The domain was set to have 100000 elements and the absolute tolerance was set as 10−8. The selection of the number of elements was based on the minimum number of elements required to satisfy the required tolerance. 5.3.3 Simulation Parameters In the simulations of both the simplified equations and the full equations the values chosen are β = 1 and 𝑈 = 4. The value of 𝐿∗ = −1 in the full simulations. In the simplified theory a study 63 was conducted to study different values of 𝐿∗. The full simulations required the additional parameters – Landau viscosities that were chosen as ν1∗ = 1, ν2∗ = −1 and ν4∗ = 6 [59]. 5.4 Results and Discussion Using the methodology described in section 5.3, the numerical computations were performed. In this section first the results from the validation are presented. It is followed by the results of the parametric studies performed on the simplified equations. 5.4.1 Validation of the Solver The validation of the solver is critical in any numerical study and as the thin film approximation causes the elimination of terms it is imperative to ensure that the resulting equations are still valid and result in physical solutions. In order to validate these simplified equations, the results of the orientation obtained are compared with the orientation profile obtained from the solution of the complete set of equations. For the sake of completeness, the validation of the solution obtained from the complete set of equations is also provided in the following sub-section. 5.4.1.1 Validation of the Solution Obtained for the Full Equations The COMSOL model has been previously implemented and validated in [9], [61] to obtain the flow solution of LCs using the LdG theory. In the present study, Planar Couette flow is the simulated. From Figure 5-2, it is seen that the velocity profile in the domain resembles the Planar Couette flow velocity profile. Figure 5-3 represent the Pressure profile in the domain and it is observed that the value of pressure obtained in the domain is of the order 10−3 and the profile is independent of x. 64 5.4.1.2 Validation of the Solution Obtained for the Simplified Equations The solution from the simplified planar lubrication theory was compared with the solutions obtained from the full theory. To confirm that the equations are valid, only the simplified evolution equations are simulated by prescribing a velocity field and then by comparing the solution with the solution obtained from the simulation of the complete equations. As mentioned previously, Planar Couette flow was the velocity profile used in these simulations. The orientation of the LCs in the domain is the physical variable of importance and was used for comparison. Figure 5-2: Planar Couette flow velocity profile obtained from solution of full equations for the case Re=0.0001, Er=1, R=0.1. The image shows the typical linear variation 65 Figure 5-3: Pressure profile obtained from solution of full equations for the Re=0.0001, Er=1, R=0.1 66 Figure 5-4: Comparison of solutions of full and simplified equations at H=1, L=10, Er=1, R=0.1, dudy=1. The orientation profiles are significantly different for different solutions Figure 5-5: Comparison of solutions of full and simplified equations at H=1, L=10, Er=10, R=1, dudy=1. The orientation profiles are significantly different for different solutions 67 Figure 5-6: Comparison of solutions of full and simplified equations at H=0.1, L=10, Er=50, R=1, dudy=1. The orientation profiles obtained for both the solutions do not coincide Figure 5-7: Comparison of solutions of full and simplified equations at H=0.1, L=10, Er=50, R=1, dudy=10. The orientation profiles obtained for both the solutions do not coincide 68 Figure 5-8: Comparison of solutions of full and simplified equations at H=0.01, L=10, Er=10, R=1, dudy=10. The orientation profiles overlap in this case. Figure 5-9: Comparison of solutions of full and simplified equations at H=0.01, L=10, Er=10, R=1, dudy=100. The orientation profiles overlap in this case. 69 5.4.2 Parametric Studies on the Simplified Equations Through the discussion in section 5.4.1.2 it was observed that reduction in the value of H results in the orientation profile (θ) from the simplified equations becoming closer to the value obtained from the solution of the complete set of equations. It was observed that the value of H=0.01 was a suitable gap length to apply the simplified equations. Therefore, parametric studies were conducted using the gap length H=0.01 to understand the effect of changes in parameters and this subsection contains the obtained results along with the discussion. Before presenting the results of the parametric studies it is important to mention the types of solutions that we can expect. Since the studies are planar, only the in-plane orientation modes can be observed. The in-plane orientation modes are – the elastic state in which the orientation profile is affected largely by the anchoring value and the viscous state in which the viscous effects dominate and the orientation profile deviates from the value at the boundary [60]. 5.4.2.1 Effect of the Ericksen number on Orientation In this study, the variation of the orientation profile was studied by varying the Ericksen number. The value of the Energy ratio remained unchanged throughout the study. Figure 5-10 shows the results of the study and it was observed that the increase in the value of the Ericksen number resulted in an increasing deviation from the boundary value. This observation was consistent for all the cases that were attempted while ensuring that the value of 𝑅 was low enough. The reason for selection of a low 𝑅 was based on the discussion in [60] and the in-plane viscous state was desired. However, the present results show the parabolic profile which is a feature of the elasticity driven steady state [83]. 70 5.4.2.2 Effect of the Energy ratio on Orientation In this study, the variation of the orientation profile was studied by varying the Energy ratio. The Energy ratio is the ratio of short-range and long-range elasticity. The value of the Ericksen number remain unchanged throughout the study. Figure 5-11 shows the results of the study and it was observed that the increase in the value of the Energy ratio does not change the deviation from the boundary value. This observation was consistent for all the cases that were tested. The constant profile obtained for different R was previously seen for the elasticity driven steady state [83]. However, through the simulations it was observed that increasing the value of 𝑅 beyond a certain value resulted in divergence of the solver and therefore the case where the orientation profile dominated by the equilibrium value was not studied. The jump in the deviation Figure 5-10: Variation of 𝛉 at R=0.001 with increase in Er showing an increase in deviation from the anchoring value 71 for 𝑅 = 102 could be due to the transition from in-plane to out of plane mode. Further increase in the value causes solver divergence and it could be due to orientation going to an out of plane mode. However, these results are not conclusive and there is a need to perform a study of the order parameters to understand the observed behavior. 5.4.3 Study to Determine the Physical Parametric Space In the previous set of studies, the effect of various parameters on the variation in the orientation profile were studied. However, those studies were performed based on non-dimensional equations with non-dimensional numbers. Therefore, it is of interest to understand how these studies can be linked to physical parameters and constants. In order to keep this section brief, the procedure for the analysis is described in the Appendix C. To summarize, the values of Figure 5-11: Variation of 𝛉 at Er=100000 for different R at unit shear rate and H=0.01. The orientation profile obtained shows that the value of orientation across the domain is constant until a transition occurs at R=100 72 physical constants were changed until the non-dimensional values in the regime of current interest were obtained. The findings are summarized in this section. The study was performed using the parameters of a lyotropic LC [64]. In this study, attempt is made to link the rotational diffusivity, mass fraction and the obtained non-dimensional numbers with the previously conducted parametric studies. Figure 5-12 represents the plot of data from Table i, Table ii and Table iii in Appendix C. However, the range of 𝐸𝑟 and 𝑅 is restricted which causes no result for δ𝐿 = 10−2 in the scatter plot. The key observation is that the above set of parameters will result in the in-plane elastic state. The values of rotational diffusivity that were tested was the range seen in literature. From the results it can be concluded that δ𝐿 = 10−3 seems to be the optimum scale at which in-plane Figure 5-12: Representation of the non-dimensional numbers for different film thickness showing that the present theory is applicable to 𝛅𝐋=0.001 (L=1) 73 solutions occur. An increase in thickness is likely to cause out-of-plane solutions and their study is beyond the scope of this thesis. 5.5 Summary In this chapter, the LdG theory was simplified using the scaling analysis and the equations for lubrication were derived. The derived equations were validated by comparison of the orientation profile obtained from the solution obtained from solving the full set of equations representing the LdG theory. The study of the validation cases showed that the results from the simplified theory approximate the results of the full theory better with decrease in the gap. Once the simplified theory is validated it was used to study the effect of different parameters on the orientation profile of the LCs. The study of variation in the Ericksen number shows that the deviation increases with an increase in the Ericksen number until a maximum value is reached. Further increase caused a reduction in deviation and on further increase, solver instability. Through the study of the Energy ratio, it was determined that the orientation profile was not significantly affected by the value of the Energy ratio. Both the above results also showed that the results lie in the Elastic-driven in plane solutions regime. Finally, an attempt was made to link the physical constants with the non-dimensional numbers, and it was seen that the most overlapping cases with the present regime were obtained for the film thickness of 1 𝑚𝑚 for a 1 𝑚 long domain. 74 Chapter 6: Conclusion In this thesis, the work has been directed towards studying LCs in different flow conditions. In this chapter, the highlights of each chapter are summarized. It also includes the contributions of the research, limitations, and directions for future work. 6.1 Summary In Chapter 3, the goal of studying the Couette flow of Graphene Oxide suspensions in water was accomplished. Graphene Oxide suspensions form a discotic nematic LCs and the successful implementation of a solver allows the further exploration of its applicability as a lubricant. The results from the numerical studies showed the expected shear thinning behavior. The orientation profiles obtained at different shear rates showed that at lower shear rates the elastic forces dominate the orientation profile of the director and as the shear rate is increased the viscous forces takeover and the orientation in bulk of the domain reaches the alignment angle. Further, it was observed that there are multiple solution branches that can be obtained depending on the initial guess of the solution. Since the Frank’s elastic constants for GO are not available in literature a sensitivity analysis was performed to determine the effect of different ratios of the elasticity coefficients. In Chapter 4, the goal of implementing an in-house 2-dimensional solver for the flow of nematic LCs using LE theory was achieved. The solver was validated using an analytical solution. Parametric studies were conducted to further understand the capability of the solver. A study to understand the effect of the Ericksen number on the orientation of the LCs showed that the increase in 𝐸𝑟 causes the orientation profile to move closer to the flow alignment value. The study served as a secondary validation case for the solver. The study at different anchoring angles showed the presence of multiple solution branches. The preliminary study of the pressure driven flow with a 75 moving wall was conducted. It was observed that the orientation profile obtained is dependent on the local shear rate and the direction of shear gradient. In Chapter 5, the goal of simplification of the Landau-de Gennes theory using scaling analysis to obtain equations for thin film flows was achieved. The obtained simplified set of equations were validated by comparison of with the solution obtained from the complete set of equations. Parametric studies were conducted to understand their effect on the orientation profiles obtained. These studies involved changes in the balance of short-range elasticity, long range elasticity, and viscous forces on the orientation tensor and therefore the orientation. The solution from the simplified equations suggested that the results from the currently studied parameters were in the Elasticity-driven steady state regime. A study of the physical constants used with the non-dimensional numbers showed that the most overlapping cases with the present regime were obtained for the film thickness of 1 𝑚𝑚 for a 1 𝑚 long domain. 6.2 Contributions The numerical study of Couette flow of GO suspensions in water was performed. GO suspension is a discotic nematic LC and as per my knowledge this was the first such implementation of a numerical study of GO. In the study, the numerically obtained values of alignment viscosity showed agreement with the experimental values. Besides gaining an understanding of the GO suspension, the successful implementation on a simple model is a stepping stone towards studies of GO in complex geometries and under different flow conditions. The development of an in-house solver for 2-dimensional flows of nematic LCs will aid the further studies and the solver is capable of handling combination of flow conditions. The numerical study of a combination of Couette and pressure driven flow of LCs was performed. 76 Further, the solver opens an avenue to implementing in-house numerical solvers to study the 3-dimensional flows of nematic LCs. The lubrication theory for the tensorial LdG theory was derived and it was shown that the results obtained from the simplified theory match those obtained from the solution of the complete set of equations. The developed theory is a novel work and there have been no previous implementation of scaling analysis with the LdG theory. The new theory allows the quicker implementation of studies of thin confined flows of nematic LCs. 6.3 Limitations The limitations of the present work are reported in this subsection. Since there are three different chapters with separate methodologies, the limitations are listed under different paragraphs. • The solution methodology involves an initial guess to be provided. During the studies it was observed that obtaining a suitable initial guess for the orientation angle in the domain is a challenge. Selection of a guess away from the solution branches results in divergence. Further, due to a lack of knowledge about the Frank’s elasticity coefficients of GO, obtaining the value of the dimensional value of shear rate is arbitrary. • The solver for channel flow of LCs considers the assumption of equal Frank’s elasticity coefficients which is not the case in reality. Although, there are studies performed previously stating that this assumption is justified, it is still a limitation. Further, in the solution procedure setting the value of pressure gradient is not trivial. • The main limitation of the work in derivation of simplified equations of lubrication theory is that the study is performed for planar cases which results in lack of completeness when 77 it comes to studying the various flow conditions. Also, the lack of coupling of the evolution equations with the Navier Stokes equation is a limitation. 6.4 Future Work The future work in the study of LCs can be take multiple directions. The potential options for further research in numerical modelling of LCs in confined domains are as follows- • Further exploration of GO suspensions and studying its load handling capabilities for application as a lubricant. • Studying GO suspension LC in different geometries and flow conditions. • Using the 2-dimensional Leslie - Ericksen solver to better understand the inertial effects on the flows of LC. • Extending the 2-dimensional Leslie - Ericksen solver to 3 dimensions for studying complex flow profiles and also the flows of cholesteric LCs. • Extending the 2-dimensional Leslie - Ericksen solver to simulate free surface flows of LCs. A further step could be to explore multiphase flows of LCs. • Using the developed Lubrication theory to study the flow of nematic LCs in typical lubrication geometries. • Extending the Lubrication theory to 3 dimensions and conducting studies for different geometries. Also, studies can be performed to gain understanding of the behavior of stress, order parameters in thin films of LCs. 78 Bibliography [1] A. D. Rey and M. M. Denn, “Dynamical phenomena in liquid-crystalline materials,” Annual Review of Fluid Mechanics, vol. 34, no. 1, pp. 233–266, 2002. [2] A. D. Rey, “Flow and texture modeling of liquid-crystalline materials,” Rheology Reviews, vol. 6, pp. 71–135, 2009. [3] J. Q. Carou, B. R. Duffy, N. J. Mottram, and S. K. Wilson, “Shear-driven and pressure-driven flow of a nematic liquid crystal in a slowly varying channel,” Physics of Fluids, vol. 18, no. 2, p. 027105, 2006. [4] A. Sengupta, “Flow of Nematic Liquid Crystals in a Microfluidic Environment,” in Topological Microfluidics, Cham: Springer International Publishing, 2013, pp. 83–135. [5] H. P. Jost, “Tribology — origin and future,” Wear, vol. 136, no. 1, pp. 1–17, 1990. [6] B. Bhushan, Introduction to tribology. John Wiley & Sons, 2013. [7] Z. M. Jin and D. Dowson, “A full numerical analysis of hydrodynamic lubrication in artificial hip joint replacements constructed from hard materials,” Proceedings of the Institution of Mechanical Engineers. Part C, Journal of Mechanical Engineering Science, vol. 213, no. 4, pp. 355–370, 1999. [8] E. P. J. Waters, P. L. Spedding, A. P. Doherty, and R. L. Spedding, “Wear of the artificial hip joint material under lubrication,” Asia-Pacific Journal of Chemical Engineering, vol. 4, no. 1, pp. 80–89, 2009. [9] N. Noroozi and D. Grecov, “Numerical simulation of three-dimensional flow-induced microstructure in a simplified prosthetic hip joint with nematic liquid crystal lubricant,” Rheologica Acta, vol. 53, no. 5–6, pp. 457–465, 2014. 79 [10] M. J. Shariatzadeh and D. Grecov, “Cellulose Nanocrystals Suspensions as Water-Based Lubricants for Slurry Pump Gland Seals,” International Journal of Aerospace and Mechanical Engineering, vol. 12, no. 6, pp. 586–590, 2018. [11] N. Noroozi and D. Grecov, “Flow modelling and rheological characterization of nematic liquid crystals between concentric cylinders,” Liquid Crystals, vol. 40, no. 7, pp. 871–883, 2013. [12] F. J. Carrión, G. Martínez-Nicolás, P. Iglesias, J. Sanes, and M. D. Bermúdez, “Liquid crystals in tribology,” International Journal of Molecular Sciences, vol. 10, no. 9. pp. 4102–4115, 2009. [13] B. I. Kupchinov, V. G. Rodnenkov, S. F. Ermakov, and V. P. Parkalov, “A study of lubrication by liquid crystals,” Tribology International, vol. 24, no. 1, pp. 25–28, 1991. [14] H. A. Elemsimit and D. Grecov, “Impact of liquid crystal additives on a canola oil-based bio-lubricant,” Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 229, no. 2, pp. 126–135, 2015. [15] J. E. Kim et al., “Graphene oxide liquid crystals,” Angewandte Chemie - International Edition, vol. 50, no. 13, pp. 3043–3047, 2011. [16] W. Song, I. A. Kinloch, and A. H. Windle, “Nematic liquid crystallinity of multiwall carbon nanotubes,” Science, vol. 302, no. 5649, p. 1363, 2003. [17] P. Kumar, U. N. Maiti, K. E. Lee, and S. O. Kim, “Rheological properties of graphene oxide liquid crystal,” Carbon, vol. 80, no. 1, pp. 453–461, 2014. [18] S. Al-Zangana, M. Iliut, M. Turner, A. Vijayaraghavan, and I. Dierking, “Confinement effects on lyotropic nematic liquid crystal phases of graphene oxide dispersions,” 2D Materials, vol. 4, no. 4, p. 041004, 2017. 80 [19] R. Jalili et al., “Organic solvent-based graphene oxide liquid crystals: a facile route toward the next generation of self-assembled layer-by-layer multifunctional 3D architectures,” ACS Nano, vol. 7, no. 5, pp. 3981–3990, 2013. [20] S. Chaiyakun et al., “Preparation and characterization of graphene oxide nanosheets,” Procedia Engineering, vol. 32, no. 7152, pp. 759–764, 2012. [21] F. Del Giudice, B. V. Cunning, R. S. Ruoff, and A. Q. Shen, “Filling the gap between transient and steady shear rheology of aqueous graphene oxide dispersions,” Rheologica Acta, vol. 57, no. 4, pp. 293–306, 2018. [22] F. Lin, X. Tong, Y. Wang, J. Bao, and Z. M. Wang, “Graphene oxide liquid crystals: synthesis, phase transition, rheological property, and applications in optoelectronics and display,” Nanoscale Research Letters, vol. 10, no. 1, pp. 1–16, 2015. [23] M. Javadi, S. Naficy, S. Beirne, S. Sayyar, R. Jalili, and S. E. Moulton, “Ionic interactions to tune mechanical and electrical properties of hydrated liquid crystal graphene oxide films,” Materials Chemistry and Physics, vol. 186, pp. 90–97, 2017. [24] A. Akbari et al., “Large-area graphene-based nanofiltration membranes by shear alignment of discotic nematic liquid crystals of graphene oxide,” Nature Communications, vol. 7, pp. 1–12, 2016. [25] H. Kinoshita, Y. Nishina, A. A. Alias, and M. Fujii, “Tribological properties of monolayer graphene oxide sheets as water-based lubricant additives,” Carbon, vol. 66, pp. 720–723, 2014. [26] Z. Chen, Y. Liu, and J. Luo, “Tribological properties of few-layer graphene oxide sheets as oil-based lubricant additives,” Chinese Journal of Mechanical Engineering, vol. 29, no. 2, pp. 439–444, 2016. 81 [27] J. Li, X. Zeng, T. Ren, and E. Van Der Heide, “The preparation of graphene oxide and its derivatives and their application in bio-tribological systems,” Lubricants, vol. 2, no. 3, pp. 137–161, 2014. [28] M. P. Allen and M. R. Wilson, “Computer simulation of liquid crystals,” Journal of Computer-Aided Molecular Design, vol. 3, pp. 335–353, 1989. [29] R. G. Larson, The Structure and Rheology of Complex Fluids. Oxford Univ. Press, 1999. [30] G. Goldbeck-Wood, P. Coulter, J. R. Hobdell, M. S. Lavine, K. Yonetake, and A. H. Windle, “Modelling of liquid crystal polymers at different length scales,” Molecular Simulation, vol. 21, no. 2–3, pp. 143–160, 1998. [31] P. De Gennes and Prost, The Physics of Liquid Crystals, vol. 3rd. Oxford university press, 1996. [32] A. Saupe, “Liquid crystals,” Annual Review of Physical Chemistry, vol. 24, no. 1, pp. 441–471, 1973. [33] S. Chandrasekhar, Liquid Crystals, 2nd ed. Cambridge: Cambridge University Press, 1992. [34] M. Hieber and J. Prüss, “Dynamics of the Ericksen–Leslie equations with general Leslie stress I: the incompressible isotropic case,” Mathematische Annalen, vol. 369, no. 3–4, pp. 977–996, 2017. [35] C. Jin, “On the Existence of Weak Solutions to the Steady Compressible Flow of Nematic Liquid Crystals,” Journal of Mathematical Fluid Mechanics, pp. 1–21, 2018. [36] R. Lasarzik, “Dissipative solution to the Ericksen–Leslie system equipped with the Oseen–Frank energy,” Zeitschrift für angewandte Mathematik und Physik, vol. 70, no. 1, p. 8, 2019. [37] J. W. Goodby, “Introduction to Defect Textures in Liquid Crystals,” in Handbook of Visual Display Technology, J. Chen, W. Cranton, and M. Fihn, Eds. Berlin, Heidelberg: Springer 82 Berlin Heidelberg, 2012, pp. 1289–1314. [38] M. G. Forest, R. Zhou, and Q. Wang, “Chaotic boundaries of nematic polymers in mixed shear and extensional flows,” Physical Review Letters, vol. 93, no. 8, pp. 4–7, 2004. [39] H. Zhou, M. G. Forest, and H. Wang, “Mathematical studies and simulations of nematic liquid crystal polymers and nanocomposites,” Journal of Computational and Theoretical Nanoscience, vol. 7, no. 4, pp. 645–660, 2010. [40] M. G. Forest, Q. Wang, and H. Zhou, “Exact banded patterns from a Doi-Marrucci-Greco model of nematic liquid crystal polymers,” Physical Review E, vol. 61, no. 6, p. 6655, 2000. [41] C. M. Care and D. J. Cleaver, “Computer simulation of liquid crystals,” Reports on progress in physics, vol. 68, no. 11, p. 2665, 2005. [42] F. M. Leslie, “Some constitutive equations for liquid crystals,” Archive for Rational Mechanics and Analysis, vol. 28, no. 4, pp. 265–283, 1968. [43] J. L. Ericksen, “Anisotropic Fluids,” Archive for Rational Mechanics and Analysis, vol. 4, no. 1, p. 231, 1959. [44] F. C. Frank, “Liquid Crystals: On the theory of liquid crystals,” Disc. Faraday. Soc., vol. 25, no. I, p. 19, 1958. [45] J. J. Feng, “Theoretical aspects of liquid crystals and liquid crystalline polymers,” Encyclopedia of Chemical Processing, pp. 2955–2964, 2006. [46] J. Feng and L. G. Leal, “Pressure-driven channel flows of a model liquid-crystalline polymer,” Physics of Fluids, vol. 11, no. 10, pp. 2821–2835, 1999. [47] L. R. P. de Andrade Lima and A. D. Rey, “Poiseuille flow of Leslie-Ericksen discotic liquid crystal: solution multiplicity, multistability, and non-Newtonian rheology,” Journal of Non-Newtonian Fluid Mechanics, vol. 110, no. 2–3, pp. 103–142, 2003. 83 [48] R. J. Atkin and F. M. Leslie, “Couette flow of nematic liquid crystals,” Quarterly Journal of Mechanics and Applied Mathematics, vol. 23, no. 2, pp. 3–24, 1970. [49] R. J. Atkin, “Poiseuille flow of liquid crystals of the nematic type,” Archive for Rational Mechanics and Analysis, vol. 38, no. 3, pp. 224–240, 1970. [50] A. S. K. Ho and A. D. Rey, “Orienting properties of discotic nematic liquid crystals in Jeffrey-Hamel flows,” Rheologica Acta, vol. 30, no. 1, pp. 77–88, 1991. [51] S. Chono, T. Tsuji, and M. M. Denn, “Spatial development of director orientation of tumbling nematic liquid crystals in pressure-driven channel flow,” Journal of non-newtonian fluid mechanics, vol. 79, no. 2–3, pp. 515–527, 1998. [52] P. A. Cruz, M. F. Tomé, I. W. Stewart, and S. Mckee, “Numerical Investigation of director orientation and flow of nematic liquid crystals in a planar 1:4 expansion,” Journal of Mechanics of Materials and Structures, vol. 6, no. 7, pp. 1007–1030, 2011. [53] P. A. Cruz, M. F. Tomé, I. W. Stewart, and S. McKee, “Numerical solution of the Ericksen-Leslie dynamic equations for two-dimensional nematic liquid crystal flows,” Journal of Computational Physics, vol. 247, pp. 109–136, 2013. [54] P. A. Cruz, M. F. Tomé, I. W. Stewart, and S. McKee, “A numerical method for solving the dynamic three-dimensional Ericksen–Leslie equations for nematic liquid crystals subject to a strong magnetic field,” Journal of Non-Newtonian Fluid Mechanics, vol. 165, no. 3–4, pp. 143–157, 2010. [55] T. G. Anderson, E. Mema, L. Kondic, and L. J. Cummings, “Transitions in Poiseuille flow of nematic liquid crystal,” International Journal of Non-Linear Mechanics, vol. 75, pp. 15–21, 2015. [56] Y. Farhoudi and A. D. Rey, “Shear flows of nematic polymers. I. Orienting modes, 84 bifurcations, and steady state rheological predictions,” Journal of Rheology, vol. 37, no. 2, pp. 289–314, 1993. [57] A. P. Singh and A. D. Rey, “Microstructure constitutive equation for discotic nematic liquid crystalline materials. Part II. Microstructure-rheology relations,” Rheologica Acta, vol. 37, no. 4, pp. 374–386, 1998. [58] T. Tsuji and A. D. Rey, “Effect of long range order on sheared liquid crystalline materials: Part 1: Compatibility between tumbling behavior and fixed anchoring,” Journal of Non-Newtonian Fluid Mechanics, vol. 73, no. 1–2, pp. 127–152, 1997. [59] D. Grecov and A. D. Rey, “Theoretical and computational rheology for discotic nematic liquid crystals,” Molecular Crystals and Liquid Crystals, vol. 391, no. February, pp. 57–94, 2002. [60] A. D. Rey and T. Tsuji, “Recent advances in theoretical liquid crystal rheology,” Macromolecular Theory and Simulations, vol. 7, no. 6. pp. 623–639, 1998. [61] N. Noroozi and D. Grecov, “Numerical simulation of flow and structure in nematic liquid crystalline materials between eccentric cylinders,” Journal of Non-Newtonian Fluid Mechanics, vol. 208–209, pp. 1–17, 2014. [62] D. Grecov and A. D. Rey, “Transient rheology of discotic mesophases,” Rheologica Acta, vol. 42, no. 6, pp. 590–604, 2003. [63] D. G. Venhaus, K. S. Conatser, and M. J. Green, “Dynamics of chiral liquid crystals under applied shear,” Liquid Crystals, vol. 40, no. 6, pp. 846–853, 2013. [64] N. Noroozi, D. Grecov, and S. Shafiei-Sabet, “Estimation of viscosity coefficients and rheological functions of nanocrystalline cellulose aqueous suspensions,” Liquid Crystals, vol. 41, no. 1, pp. 56–66, 2014. 85 [65] M. J. Pospisil, P. Saha, S. Abdulquddos, M. M. Noor, V. A. Davis, and M. J. Green, “Orientation relaxation dynamics in cellulose nanocrystal dispersions in the chiral liquid crystalline phase,” Langmuir, vol. 34, no. 44, pp. 13274–13282, 2018. [66] M. Miesowicz, “The three coefficients of viscosity of anisotropic liquids,” Nature, vol. 158, no. 4001, p. 27, 1946. [67] L. G. Leal, Advanced transport phenomena: fluid mechanics and convective transport processes, vol. 7. Cambridge University Press, 2007. [68] R. L. Panton, Incompressible flow. John Wiley & Sons, 2013. [69] J. A. Tichy and Y. Rhim, “A Theory for the lubrication of layered liquid crystals,” Jounal of Tribology, vol. 1, no. January, 1989. [70] J. A. Tichy, “Lubrication theory for nematic liquid crystals,” Tribology Transactions, vol. 33, no. 3, pp. 363–370, 1990. [71] M. Ben Amar and L. J. Cummings, “Fingering instabilities in driven thin nematic films,” Physics of Fluids, vol. 13, no. 5, pp. 1160–1166, 2001. [72] L. J. Cummings, “Evolution of a thin film of nematic liquid crystal with anisotropic surface energy,” European Journal of Applied Mathematics, vol. 15, no. 6, pp. 651–677, 2004. [73] V. M. Ugaz, O. K. Cinader, and W. R. Burghardt, “Origins of region I shear thinning in model lyotropic liquid crystalline polymers,” Macromolecules, vol. 30, no. 5, pp. 1527–1530, 1997. [74] J. P. Straley, “Frank elastic constants of the hard-rod liquid crystal,” Physical Review A, vol. 8, no. 4, p. 2181, 1973. [75] M. J. Bradshaw, E. P. Raynes, J. D. Bunning, and T. E. Faber, “The Frank constants of some nematic liquid crystals,” Journal de Physique, vol. 46, no. 9, pp. 1513–1520, 1985. 86 [76] M. F. Tome and S. McKee, “GENSMAC: A computational marker and cell method for free surface flows in general domains,” Journal of Computational Physics, vol. 110, no. 1, pp. 171–186, 1994. [77] S. McKee, M. F. Tomé, J. A. Cuminato, A. Castelo, and V. G. Ferreira, “Recent advances in the marker and cell method,” Archives of Computational Methods in Engineering, vol. 11, no. 2. pp. 107–142, 2004. [78] S. Patankar, Numerical heat transfer and fluid flow. CRC press, 1980. [79] M. A. Alves, P. J. Oliveira, and F. T. Pinho, “A convergent and universally bounded interpolation scheme for the treatment of advection,” International journal for numerical methods in fluids, vol. 41, no. 1, pp. 47–75, 2003. [80] P. J. Oliveira, “Asymmetric flows of viscoelastic fluids in symmetric planar expansion geometries,” Journal of non-newtonian fluid mechanics, vol. 114, no. 1, pp. 33–63, 2003. [81] M. A. Alves, P. J. Oliveira, and F. T. Pinho, “Benchmark solutions for the flow of Oldroyd-B and PTT fluids in planar contractions,” Journal of Non-Newtonian Fluid Mechanics, vol. 110, no. 1, pp. 45–75, 2003. [82] M. J. Stephen and J. P. Straley, “Physics of liquid crystals,” Reviews of Modern Physics, vol. 46, no. 4, p. 617, 1974. [83] A. P. Singh and A. D. Rey, “Effect of long-range elasticity and boundary conditions on microstructural response of sheared discotic mesophases,” Journal of Non-Newtonian Fluid Mechanics, vol. 94, no. 2–3, pp. 87–111, 2000. [84] D. Grecov and A. D. Rey, “Shear-induced textural transitions in flow-aligning liquid crystal polymers,” Physical Review E, vol. 68, no. 6, p. 61704, 2003. [85] M. W. Johnson and S. Mangkoesoebroto, “Analysis of lubrication theory for the power law 87 fluid,” Journal of tribology, vol. 115, no. 1, pp. 71–77, 1993. [86] A. Acevedo, P. M. Cotts, and A. D. Shine, “Molecular weight dependence of the rotational diffusivity of rodlike polymers in concentrated nematic solutions,” Macromolecules, vol. 38, no. 15, pp. 6648–6655, 2005. 88 Appendices Appendix A has the terms of the momentum equations used in the Chapter 4. 𝑅𝑥𝜕𝑛𝑥𝜕𝑥+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑥=1𝑅𝑒{ −𝛾1𝜕𝜃𝜕𝑥[𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)] − 12𝛾2𝜕𝜃𝜕𝑥[𝑐𝑜𝑠(2𝜃) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) + 𝑠𝑖𝑛(2𝜃) (𝜕𝑢𝜕𝑥−𝜕𝑣𝜕𝑦)] } (A. 1) 𝑅𝑥𝜕𝑛𝑥𝜕𝑦+ 𝑅𝑦𝜕𝑛𝑦𝜕𝑦=1𝑅𝑒{ −𝛾1𝜕𝜃𝜕𝑦[𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)] − 12𝛾2𝜕𝜃𝜕𝑦[𝑐𝑜𝑠(2𝜃) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) + 𝑠𝑖𝑛(2𝜃) (𝜕𝑢𝜕𝑥−𝜕𝑣𝜕𝑦)] } (A. 2) 𝚽𝐱𝐱 = 𝛼1 𝑐𝑜𝑠2(𝜃) [𝑑𝑢𝑑𝑥𝑐𝑜𝑠2(𝜃) +𝜕𝑣𝜕𝑦𝑠𝑖𝑛2(𝜃) +12(𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) 𝑠𝑖𝑛(2𝜃)]−(𝛼2 + 𝛼3) 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) [𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)]+(𝛼5 + 𝛼6) [𝑑𝑢𝑑𝑥𝑐𝑜𝑠2(𝜃) +12𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥)] (A. 3) 𝚽𝐱𝐲 = 𝛼1 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) [𝑑𝑢𝑑𝑥𝑐𝑜𝑠2(𝜃) +𝜕𝑣𝜕𝑦𝑠𝑖𝑛2(𝜃) +12(𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) 𝑠𝑖𝑛(2𝜃)]+(𝛼3 𝑐𝑜𝑠2(𝜃) − 𝛼2 𝑠𝑖𝑛2(𝜃)) [𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)]+12(𝛼5 𝑠𝑖𝑛2(𝜃) + 𝛼6 𝑐𝑜𝑠2(𝜃)) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) + (𝛼5𝑑𝑢𝑑𝑥+ 𝛼6𝜕𝑣𝜕𝑦) 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) (A. 4) 𝚽𝐲𝐱 = 𝛼1 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) [𝑑𝑢𝑑𝑥𝑐𝑜𝑠2(𝜃) +𝜕𝑣𝜕𝑦𝑠𝑖𝑛2(𝜃) +12(𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) 𝑠𝑖𝑛(2𝜃)]+(𝛼3 𝑐𝑜𝑠2(𝜃) − 𝛼2 𝑠𝑖𝑛2(𝜃)) [𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)]+12(𝛼5 𝑐𝑜𝑠2(𝜃) + 𝛼6 𝑠𝑖𝑛2(𝜃)) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) + (𝛼5𝑑𝑢𝑑𝑥+ 𝛼6𝜕𝑣𝜕𝑦) 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) (A. 5) 89 𝚽𝐲𝐲 = 𝛼1 𝑠𝑖𝑛2(𝜃) [𝑑𝑢𝑑𝑥𝑐𝑜𝑠2(𝜃) +𝜕𝑣𝜕𝑦𝑠𝑖𝑛2(𝜃) +12(𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥) 𝑠𝑖𝑛(2𝜃)]+(𝛼2 + 𝛼3) 𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) [𝜕𝜃𝜕𝑡+ 𝑢𝜕𝜃𝜕𝑥+ 𝑣𝜕𝜃𝜕𝑦+12(𝜕𝑢𝜕𝑦−𝜕𝑣𝜕𝑥)]+(𝛼5 + 𝛼6) [𝜕𝑣𝜕𝑦𝑠𝑖𝑛2(𝜃) +12𝑠𝑖𝑛(𝜃) 𝑐𝑜𝑠(𝜃) (𝜕𝑢𝜕𝑦+𝜕𝑣𝜕𝑥)] (A. 6) 90 Appendix B contains the step by step procedure involved in the simplification of the Landau-de Gennes theory to obtain the Lubrication theory. The governing equations for the Navier Stokes equations and the evolution equation are already described in the Chapter 5 and will not be repeated in this section. It is important to note that the simplification is performed using the scaling analysis and under the assumption that the flow is planar and is in the X-Y plane. Understanding the scaling analysis and non-dimensionalization is crucial. As mentioned in Chapter 5, the length scale along 𝑥 is taken as 𝐿 and the one along 𝑦 is taken as δ𝐿 (knowing that 𝛿 ≪ 1). During the process of non-dimensionalization the terms along 𝑥 and 𝑦 get non-dimensionalized differently as will be shown in the following few steps. Since the dominant value of velocity will be in the 𝑥 direction let us consider the value 𝑈 (say the maximum velocity of shear) as the velocity scale. Considering the mass conservation equation first and non dimensionalizing it based on the length scales 𝐿 and δ𝐿. Also considering the scales 𝑈 and 𝑉 for non-dimensionalization of velocities we obtain- 𝑈𝐿𝜕𝑢∗𝜕𝑥∗+𝑉𝛿𝐿𝜕𝑣∗𝜕𝑦∗= 0 (B. 1) The terms with asterisk are non-dimensional. Also, on balance of scales we get 𝑉 = δ𝑈. Therefore, the final non-dimensionalized mass conservation equation takes the form (asterisk is dropped)- 𝜕𝑢𝜕𝑥+𝜕𝑣𝜕𝑦= 0 (B. 2) 91 The time scale was selected as 𝑇 = 𝐿/𝑈. Similar to the non-dimensionalization of the mass conservation equation other terms are non- dimensionalized. The implementation of this idea on the time derivative and the convection terms in the momentum equation is straightforward and is not discussed any further. However, the momentum equation has the divergence of the stress tensor and the treatment is non-trivial. Non-dimensionalization of terms of the momentum and the evolution equations- Starting first with the simplification of the evolution equations. Please note that only during this appendix 𝑼 represents the nematic potential. It is not a tensor. Recalling that the evolution equation is given by – ?̂? = 𝐅 + 𝐇𝐬𝐫 + 𝐇𝐥𝐫 (B. 3) With the flow contribution given by 𝐅 =2𝛽3𝐀 + 𝛽 [𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐀 −23(𝐀:𝐐)𝐈] −𝛽2{(𝐀:𝐐)𝐐 + 𝐀 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐀 + [(𝐐 ∙ 𝐐): 𝐀]𝐈} (B. 4) The short-range elasticity is given by- 𝐇𝐬𝐫 = 6𝐷𝑟̅̅ ̅ [(𝑼3− 1)𝐐 + 𝑼𝐐 ∙ 𝐐 − 𝑼(𝐐 ∶ 𝐐) ∙ (𝐐 +13𝐈)] (B. 5) The long-range elasticity is given by- 𝐇𝐥𝐫 =𝐿1𝜂∇2𝐐 +𝐿22𝜂{∇(∇ ∙ 𝐐) + [∇(∇ ∙ 𝐐)]𝑻 −23𝑡𝑟∇(∇ ∙ 𝐐)} (B. 6) The process of non-dimensionalization is similar to that of [59]. However, in this case, the derivatives w.r.t y will bring the additional factor of 𝛿 which then would help in the scaling analysis and subsequent elimination of terms leading to the simplified equation. 92 The non-dimensionalization terms used in [59] are- 𝑡∗ = ?̇?𝑡, 𝐀∗ = 𝐀 ?̇?⁄ ,𝐖∗ = 𝐖 ?̇?⁄ , ∇∗= 𝐿∇, 𝐿∗ =𝐿2𝐿1⁄ (B. 7) In the present case ?̇? = 𝑈 𝐿⁄ . However, as mentioned before the derivatives and the terms containing both the horizontal and vertical components of velocity are the terms that require attention. We know that the terms of A and W have gradients of velocity. To give an example, let us consider the term 𝐀12 = 0.5 (𝜕𝑣𝜕𝑥+𝜕𝑢𝜕𝑦). Now as per the proposed non-dimensionalization, the term becomes 𝐀12 = (𝑈𝐿⁄ )𝐀∗𝟏𝟐 Looking at the complete term and realizing that the velocity scale 𝑉 = 𝛿𝑈 and the length scale in y is 𝛿𝐿. Therefore, we can conclude that order of magnitude of 𝜕𝑢𝜕𝑦 is 𝑈 𝛿𝐿⁄ and the order of magnitude of 𝜕𝑣𝜕𝑥 is 𝛿𝑈 𝐿⁄ . Clearly, 𝜕𝑢𝜕𝑦 is of a higher magnitude and therefore the simplified term after non-dimensionalization can be written as 𝐀12 =12𝑈𝛿𝐿(𝜕𝑢𝜕𝑦)∗. The reason for the retention of 𝛿 is that it will later allow for the ease of comparison of terms to make the analysis easier. Now, individual simplified terms are listed as follows- 𝐴12 = 𝐴21 = 0.5 ∗𝜕𝑢𝜕y, 𝑊12 = −𝑊21 = 0.5 ∗𝜕𝑢𝜕y (𝐀 ∙ 𝐐)11 = 𝐴12𝑄21, (𝐀 ∙ 𝐐)12 = 𝐴12𝑄22, (𝐀 ∙ 𝐐)21 = 𝐴21𝑄11, (𝐀 ∙ 𝐐)22 = 𝐴21𝑄12 (𝐐 ∙ 𝐀)11 = 𝑄12𝐴21, (𝐐 ∙ 𝐀)12 = 𝑄11𝐴12, (𝐐 ∙ 𝐀)21 = 𝑄22𝐴21, (𝐐 ∙ 𝐀)22 = 𝑄21𝐴12 (𝐖 ∙ 𝐐)11 = 𝑊12𝑄21, (𝐖 ∙ 𝐐)12 = 𝑊12𝑄22, (𝐖 ∙ 𝐐)21 = 𝑊21𝑄11, (𝐖 ∙ 𝐐)22 = 𝑊21𝑄12 (𝐐 ∙ 𝐖)11 = 𝑄12𝑊21, (𝐐 ∙ 𝐖)12 = 𝑄11𝑊12, (𝐐 ∙ 𝐖)21 = 𝑄22𝑊21, (𝐐 ∙ 𝐖)22 = 𝑄21𝑊12 (𝐀 ∙ 𝐐 ∙ 𝐐)11 = 𝐴12(𝑄21𝑄11 + 𝑄22𝑄21), (𝐀 ∙ 𝐐 ∙ 𝐐)12 = 𝐴12(𝑄21𝑄12 + 𝑄22𝑄22), (𝐀 ∙ 𝐐 ∙ 𝐐)21 =𝐴21(𝑄21𝑄11 + 𝑄22𝑄21), (𝐀 ∙ 𝐐 ∙ 𝐐)22 = 𝐴21(𝑄21𝑄12 + 𝑄22𝑄22) 93 (𝐐 ∙ 𝐀 ∙ 𝐐)11 = 𝑄11𝐴12𝑄21 + 𝑄12𝐴21𝑄11, (𝐐 ∙ 𝐀 ∙ 𝐐)12 = 𝑄11𝐴12𝑄22 + 𝑄12𝐴21𝑄12, (𝐐 ∙ 𝐀 ∙𝐐)21 = 𝑄21𝐴12𝑄21 + 𝑄22𝐴21𝑄11, (𝐐 ∙ 𝐀 ∙ 𝐐)22 = 𝑄21𝐴12𝑄22 + 𝑄22𝐴21𝑄12 (𝐐 ∙ 𝐐 ∙ 𝐀)11 = (𝑄11𝑄12 + 𝑄12𝑄22)𝐴21, (𝐐 ∙ 𝐐 ∙ 𝐀)12 = (𝑄11𝑄11 + 𝑄12𝑄21)𝐴12, (𝐐 ∙ 𝐐 ∙ 𝐀)21 =(𝑄21𝑄12 + 𝑄22𝑄22)𝐴21, (𝐐 ∙ 𝐐 ∙ 𝐀)22 = (𝑄21𝑄11 + 𝑄22𝑄21)𝐴12 (𝐀 ∶ 𝐐) = 𝑄12𝐴21 + 𝑄21𝐴12 ((𝐐 ∙ 𝐐) ∶ 𝐀) = (𝑄11𝑄12 + 𝑄12𝑄22)𝐴21 + (𝑄21𝑄11 + 𝑄22𝑄21)𝐴12 (B. 8) The terms of 𝐇𝐥𝐫 are due to the combination of gradient and/or divergence of 𝐐 and similar process as above is used to simplify them. These terms are listed below- [∇(∇ ∙ 𝐐)]11 =𝜕2𝜕𝑥𝜕𝑦(𝑄12), [∇(∇ ∙ 𝐐)]12 =𝜕2𝜕𝑦2(𝑄12), [∇(∇ ∙ 𝐐)]21 =𝜕2𝜕𝑥𝜕𝑦(𝑄22)[∇(∇ ∙ 𝐐)]22 =𝜕2𝜕𝑦2(𝑄22)[∇(∇ ∙ 𝐐)]𝑇11 =𝜕2𝜕𝑥𝜕𝑦(𝑄12), [∇(∇ ∙ 𝐐)]𝑇12 =𝜕2𝜕𝑥𝜕𝑦(𝑄22), [∇(∇ ∙ 𝐐)]𝑇21 =𝜕2𝜕𝑦2(𝑄12) [∇(∇ ∙ 𝐐)]𝑇22 =𝜕2𝜕𝑦2(𝑄22)[∇2𝐐]11 =𝜕2𝜕𝑦2(𝑄11), [∇2𝐐]12 =𝜕2𝜕𝑦2(𝑄12), [∇2𝐐]21 =𝜕2𝜕𝑦2(𝑄21), [∇2𝐐]22 =𝜕2𝜕𝑦2(𝑄22)𝑡𝑟∇(∇ ∙ 𝐐) =𝜕2𝜕𝑦2(𝑄22) (B. 9) The terms of 𝐇𝐬𝐫 are a combination of different tensor products which are trivial and are not listed here. In [59] the non-dimensionalized evolution equation takes the form (asterisk sign indicating the non-dimensional terms is omitted in this equation) – ?̂? = 𝐅 + (𝟏𝐷𝑒)𝐇𝒔𝒓 + (𝟏𝐸𝑟)𝐇𝑙𝑟 (B. 10) For simplification, the following non-dimensional numbers are selected: • Ericksen number: 𝐸𝑟 =𝑈δ𝐿𝑐𝐾𝑇6𝐿1𝐷𝑟 94 • Deborah number: 𝐷𝑒 =𝑈6δ𝐿𝐷𝑟 • Energy ratio: 𝑅 =𝐸𝑟𝐷𝑒 Like the non-dimensional numbers of the usual cases, in the present study the length scale to determine the non-dimensional numbers is the length important in the calculation of shear rate. Note that the Jaumann derivative has the gradient term w.r.t 𝑦. It is also assumed that the orientation tensor is steady. Using the simplified terms that were listed earlier the final equations for the three evolution variables 𝑄11, 𝑄12, and 𝑄22 can be obtained and are written as follows- (1𝐸𝑟)(𝜕2𝜕𝑦2(𝑄11) −𝐿∗3𝜕2𝜕𝑦2(𝑄22)) +(1𝐷𝑒)((𝑼3− 1)𝑄11 + 𝑼(𝑄112 + 𝑄122) −𝑼(𝑄112 + 2𝑄122 + 𝑄222 + (𝑄11+𝑄22)2) (𝑄11 +13)) =−𝜕𝑢𝜕𝑦[1 − 3(𝑄112 + 𝑄122 + 𝑄222 +𝑄11𝑄22)]2[𝑄12 (1 +𝛽3) − 𝛽𝑄11𝑄12] (B. 11) (1𝐸𝑟)((1 +𝐿∗2)𝜕2𝜕𝑦2(𝑄12)) +(1𝐷𝑒)((𝑼3− 1)𝑄12 + 𝑼(𝑄11𝑄12 + 𝑄12𝑄22) −𝑼𝑄12(𝑄112 + 2𝑄122 +𝑄222 + (𝑄11+𝑄22)2)) =𝜕𝑢𝜕𝑦[1 − 3(𝑄112 + 𝑄122 + 𝑄222 + 𝑄11𝑄22)]2[𝑄11 (1 − 𝛽2) − 𝑄22 (𝛽 + 12) +𝛽4𝑄112 +5𝛽4𝑄122 +𝛽4𝑄222 −𝛽3+𝛽4𝑄11𝑄22] (B. 12) 95 (1𝐸𝑟)((1 +2𝐿∗3)𝜕2𝜕𝑦2(𝑄22)) +(1𝐷𝑒)((𝑼3− 1)𝑄22 +𝑼(𝑄112 + 𝑄122) −𝑼(𝑄112 + 2𝑄122 + 𝑄222 + (𝑄11+𝑄22)2) (𝑄22 +13)) =𝜕𝑢𝜕𝑦[1 − 3(𝑄112 + 𝑄122 + 𝑄222 + 𝑄11𝑄22)]2[𝑄12 (1 −𝛽3) + 𝛽𝑄12𝑄22] (B. 13) Coming to the momentum equations, let us consider the terms of the stress tensor first. The viscous, elastic, antisymmetric and the Ericksen stress terms are as follows- 𝝉𝒗 = 𝜈1𝐀 + 𝜈2 [𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐀 −23(𝐀:𝐐)𝐈] +𝜈4{(𝐀:𝐐)𝐐 + 𝐀 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐀 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐀 + [(𝐐:𝐐)𝐀]𝐈} (B. 14) 𝝉𝒆 = 𝑐𝑘𝑇 {−2𝛽3𝐇 − 𝛽 [𝐇 ∙ 𝐐 + 𝐐 ∙ 𝐇 −23(𝐇:𝐐)𝐈] + 𝛽𝟐{(𝐇:𝐐)𝐐 + 𝐇 ∙ 𝐐 ∙ 𝐐 + 𝐐 ∙ 𝐇 ∙ 𝐐 + 𝐐 ∙ 𝐐 ∙ 𝐇 + [(𝐐:𝐐)𝐈]} } (B. 15) 𝝉𝒂 = 𝑐𝑘𝑇{𝐇 ∙ 𝐐 − 𝐐 ∙ 𝐇} (B. 16) 𝝉𝑬𝒓 = 𝑐𝑘𝑇{𝐿2(∇ ∙ 𝐐 ∙ (∇𝐐)𝑇) − 𝐿1(∇𝐐: (∇𝐐)𝑇)} (B. 17) The total stress tensor is 𝛔 = 𝝉𝒗 + 𝝉𝒆 + 𝝉𝒂 + 𝝉𝑬𝒓 (B. 18) The term due to elasticity 𝐇 = 𝐇𝒔𝒓 +𝐇𝑙𝑟 is already non-dimensional. The total stress tensor can be non-dimensionalized similar methodology as [59]. However, in the simplification we have eliminated terms as per equations (B.8) and (B.9). Also, adapting the following terms- 𝛔∗ = 𝛔 (𝑐𝐾𝑇)⁄ , 𝝂𝟏∗ =(𝝂𝟏6𝐷𝑟)(𝑐𝐾𝑇)⁄ , 𝝂𝟐∗ =(𝝂𝟐6𝐷𝑟)(𝑐𝐾𝑇)⁄ , 𝝂𝟒∗ =(𝝂𝟒6𝐷𝑟)(𝑐𝐾𝑇)⁄(B. 19) At this point, let us consider the momentum equations again. They contain the divergence of the stress tensor. Similar to the treatment thus far we can neglect the terms in the divergence of 96 the stress tensor. Since we know that the derivative w.r.t 𝑦 would be much higher than the derivative w.r.t 𝑥 we can neglect the derivative terms w.r.t 𝑥. Also δ is retained for the analysis going forward. For the time being let us assume that the pressure is non-dimensionalized using a scale 𝑃0. This will be revisited after writing the non-dimensionalized momentum equations. Choosing the Reynolds number as 𝑅𝑒 =𝜌𝑈2𝑐𝐾𝑇 with 𝜌 as the characteristic density the non-dimensional form of the momentum equations is as shown (asterisk sign indicating the non-dimensional terms is omitted in this equation) - 𝛿3𝑅𝑒 {𝜕𝑢𝜕𝑡+ 𝑢𝜕𝑢𝜕𝑥+ 𝑣𝜕𝑢𝜕𝑦} =𝛿3𝑃0𝑐𝐾𝑇(−𝜕𝑝𝜕𝑥) + 𝛿𝐷𝑒𝜕𝜏12𝑣𝜕𝑦+𝜕𝜏12𝐸𝜕𝑦+𝜕𝜏12𝑎𝜕𝑦+1𝑅𝜕𝜏12𝐸𝑟𝜕𝑦(B. 20) 𝛿4𝑅𝑒 {𝜕v𝜕𝑡+ u𝜕v𝜕𝑥+ 𝑣𝜕v𝜕𝑦} =𝛿2𝑃0𝑐𝐾𝑇(−𝜕𝑝𝜕y) + 𝛿𝐷𝑒𝜕𝜏22𝑣𝜕𝑦+𝜕𝜏22𝐸𝜕𝑦+𝜕𝜏22𝑎𝜕𝑦+1𝑅𝜕𝜏22𝐸𝑟𝜕𝑦(B. 21) The terms can now be neglected based on their orders. However, before the neglection of terms the scale for Pressure 𝑃0 has to be selected and based on the previous work done in lubrication theory the appropriate choice is 𝑃0 =𝑐𝐾𝑇𝛿3 [67], [68], [85]. The terms on the LHS can be clearly neglected and after dropping the viscous terms (as they are multiplied by 𝛿) the final form of the equations are as follows- (−𝜕𝑝𝜕𝑥) +𝜕𝜏12𝐸𝜕𝑦+𝜕𝜏12𝑎𝜕𝑦+1𝑅𝜕𝜏12𝐸𝑟𝜕𝑦= 0 (B. 22) (−𝜕𝑝𝜕y) = 0 (B. 23) This shows that the elastic, anisotropic and Ericksen stress terms are the terms remaining from the stress tensor that work in balance with the pressure gradient. Moreover, it is clear that 𝑝 97 is independent of 𝑦. The equation (B. 20) will take a form as shown in (B. 24) only in the case where 𝛿𝐷𝑒 ≈ 1. (−𝜕𝑝𝜕𝑥) +𝜕𝜏12𝑣𝜕𝑦+𝜕𝜏12𝐸𝜕𝑦+𝜕𝜏12𝑎𝜕𝑦+1𝑅𝜕𝜏12𝐸𝑟𝜕𝑦= 0 (B. 25) To conclude, equations (B. 26) – (B. 27), (B. 282), and (B. 29) form the equations of the lubrication theory applied to Landau-de Gennes theory. 98 The Appendix C describes the procedure to understand the approach to study the trend in the non-dimensional numbers and the resulting regimes that are found. In the numerical studies of LCs, it is usually a challenge to link the values of the non-dimensional numbers with the physical constants used in the description of those values. Therefore, it is of interest to link the physical constants with the obtained values. Firstly, the geometry of lubrication needs to be revisited. From Figure 5-1 it is understood that the gap is of dimension δ𝐿 as opposed to the large dimension 𝐿. It is assumed that 𝐿 = 1 𝑚. This results in δ representing the thickness of the LC film. The non-dimensional numbers are - 𝐸𝑟 =𝑈δ𝐿𝑐𝐾𝑇6𝐿1𝐷𝑟 , 𝐷𝑒 =𝑈6δ𝐿𝐷𝑟 and 𝑅 =δ2𝐿2𝑐𝑘𝑇𝐿1 Apart from the Velocity scale (𝑈) and Length scale (δ𝐿), the physical constants needed in the description are 𝑐, 𝑘, 𝑇, 𝐷𝑟, 𝐿1. For this study, 𝑘 = 1.38 ∗ 10−23𝑚2 𝑘𝑔𝑠−2 𝐾−1, 𝑇 = 300 𝐾 and 𝐿1 = 2 ∗ 10−11. From Doi’s theory, it is known that 𝐷𝑟 ∝1𝑐2 and this relation is maintained. However, the constant of proportionality is required to determine the values in order to conduct the study. Following the work of [86], the values of the rotational diffusivity for a lyotropic calamitic LC can be calculated. For the present calculations the Volume fraction is 16%. The length of the LC molecule was chosen as 200 𝑛𝑚 and the thickness was taken as 5 𝑛𝑚. Having decided the physical constants or at least their possible magnitudes the next step is to compare different thicknesses and velocities. The tables consist of the results which lie in the reasonable ranges of values in terms of 𝐸𝑟 and 𝑅. The values also consider the eventual pressure scale that would determine the order of magnitude of the pressure in the domain and the values higher than 1011 were neglected. Finally, the tables are separated based on the film thickness. 99 Table i: Magnitudes of non-dimensional numbers at 𝛅 = 𝟏𝟎−𝟐 Dr 𝛅𝑳 U U/𝛅𝑳 Er R Pressure scale 1.00E+00 1.00E-02 1.00E-02 1.00E+00 8.03E+05 4.82E+06 9.64E+05 1.00E+01 1.00E-02 1.00E-01 1.00E+01 2.54E+05 1.52E+06 3.05E+05 1.00E+00 1.00E-02 1.00E-03 1.00E-01 8.03E+04 4.82E+06 9.64E+05 1.00E+01 1.00E-02 1.00E-02 1.00E+00 2.54E+04 1.52E+06 3.05E+05 1.00E+02 1.00E-02 1.00E-02 1.00E+00 8.03E+02 4.82E+05 9.64E+04 1.00E+02 1.00E-02 1.00E-03 1.00E-01 8.03E+01 4.82E+05 9.64E+04 1.00E+03 1.00E-02 1.00E-02 1.00E+00 2.54E+01 1.52E+05 3.05E+04 1.00E+03 1.00E-02 1.00E-03 1.00E-01 2.54E+00 1.52E+05 3.05E+04 Table ii: Magnitudes of non-dimensional numbers at 𝜹 = 𝟏𝟎−𝟑 Dr 𝛅𝑳 U U/𝛅𝑳 Er R Pressure scale 1.00E+00 1.00E-03 1.00E-01 1.00E+02 8.03E+05 4.82E+04 9.64E+08 1.00E-01 1.00E-03 1.00E-03 1.00E+00 2.54E+05 1.52E+05 3.05E+09 1.00E+00 1.00E-03 1.00E-02 1.00E+01 8.03E+04 4.82E+04 9.64E+08 1.00E+02 1.00E-03 1.00E+00 1.00E+03 8.03E+03 4.82E+03 9.64E+07 1.00E+00 1.00E-03 1.00E-03 1.00E+00 8.03E+03 4.82E+04 9.64E+08 1.00E+01 1.00E-03 1.00E-02 1.00E+01 2.54E+03 1.52E+04 3.05E+08 1.00E+02 1.00E-03 1.00E-01 1.00E+02 8.03E+02 4.82E+03 9.64E+07 1.00E+03 1.00E-03 1.00E+00 1.00E+03 2.54E+02 1.52E+03 3.05E+07 1.00E+01 1.00E-03 1.00E-03 1.00E+00 2.54E+02 1.52E+04 3.05E+08 1.00E+02 1.00E-03 1.00E-02 1.00E+01 8.03E+01 4.82E+03 9.64E+07 1.00E+03 1.00E-03 1.00E-01 1.00E+02 2.54E+01 1.52E+03 3.05E+07 1.00E+02 1.00E-03 1.00E-03 1.00E+00 8.03E+00 4.82E+03 9.64E+07 1.00E+03 1.00E-03 1.00E-02 1.00E+01 2.54E+00 1.52E+03 3.05E+07 1.00E+03 1.00E-03 1.00E-03 1.00E+00 2.54E-01 1.52E+03 3.05E+07 100 Table iii: Magnitudes of non-dimensional numbers at 𝜹 = 𝟏𝟎−𝟒 Dr 𝛅𝑳 U U/𝛅𝑳 Er R Pressure scale 1.00E+02 1.00E-04 1.00E+00 1.00E+04 8.03E+02 4.82E+01 9.64E+10 1.00E+02 1.00E-04 1.00E-01 1.00E+03 8.03E+01 4.82E+01 9.64E+10 1.00E+03 1.00E-04 1.00E+00 1.00E+04 2.54E+01 1.52E+01 3.05E+10 1.00E+02 1.00E-04 1.00E-02 1.00E+02 8.03E+00 4.82E+01 9.64E+10 1.00E+03 1.00E-04 1.00E-01 1.00E+03 2.54E+00 1.52E+01 3.05E+10 1.00E+02 1.00E-04 1.00E-03 1.00E+01 8.03E-01 4.82E+01 9.64E+10 1.00E+03 1.00E-04 1.00E-02 1.00E+02 2.54E-01 1.52E+01 3.05E+10 1.00E+03 1.00E-04 1.00E-03 1.00E+01 2.54E-02 1.52E+01 3.05E+10 Notes for interpreting the study: The value of 𝐷𝑟 is changed which in return changes the number of molecules of the LC. The change in number of molecules on the larger scale translates to a change in mass fraction. It is known that the reduction of 𝐷𝑟 increases the mass fraction. These changes are carried out at a constant aspect ratio.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical modelling of liquid crystal flows in confined...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical modelling of liquid crystal flows in confined domains Bhatia, Somesh 2019
pdf
Page Metadata
Item Metadata
Title | Numerical modelling of liquid crystal flows in confined domains |
Creator |
Bhatia, Somesh |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | Liquid crystals are anisotropic, viscoelastic materials with properties intermediate of solids and liquids. They are useful structural and functional materials and due to their ability to form ordered layers close to the bounding surfaces they are used as lubricants. Under the application of a hydrodynamic field, based on the type of velocity profile, values of non-dimensional numbers and anchoring angles, different orientation profiles are observed. Leslie-Ericksen and Landau-de Gennes theories are used to understand the evolution of the microstructure. Leslie-Ericksen theory, due to its simplicity can be used to obtain the behavior of flow aligning nematic liquid crystals. Landau-de Gennes theory is a mesoscopic model and apart from its ability to capture singular solutions can be employed in a study of lubrication using liquid crystals. This research work contains studies of liquid crystals in different flow conditions, such as the Couette flow and the Channel flow. In the study of Couette flow of Graphene oxide suspensions in water, the Leslie-Ericksen theory was used to obtain the orientation and viscosity profiles at different shear rates, up until flow alignment was observed. In the numerical study of pressure-driven channel flows, a Marker and Cell based solution methodology was implemented to solve the Leslie-Ericksen hydrodynamic theory. The expected flow alignment of liquid crystals was obtained which validated the solver. A preliminary study combining the moving wall and pressure driven flow showed that the orientation profile obtained depends on the local direction of shear and the direction of shear gradient. The scaling analysis was applied to Landau-de Gennes theory to derive the simplified equations of a planar lubrication theory. The theory was then validated by comparison with the solution obtained from numerically solving the full set of equations using the Couette flow profile. The parametric studies conducted showed that the solution was in the Elasticity-driven steady state. A discussion of the physical conditions showed its applicability for films of thickness lesser than 1 mm for a 1 m long domain. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-04-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0378278 |
URI | http://hdl.handle.net/2429/69751 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_may_bhatia_somesh.pdf [ 3.34MB ]
- Metadata
- JSON: 24-1.0378278.json
- JSON-LD: 24-1.0378278-ld.json
- RDF/XML (Pretty): 24-1.0378278-rdf.xml
- RDF/JSON: 24-1.0378278-rdf.json
- Turtle: 24-1.0378278-turtle.txt
- N-Triples: 24-1.0378278-rdf-ntriples.txt
- Original Record: 24-1.0378278-source.json
- Full Text
- 24-1.0378278-fulltext.txt
- Citation
- 24-1.0378278.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0378278/manifest