Numerical Simulations of Flow andMicrostructure in Nematic LiquidCrystalline MaterialsbyNader NorooziB.Sc., University of Tehran, 2005M.Sc., University of Tehran, 2007A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)September 2013c© Nader Noroozi 2013AbstractLiquid crystals are known for their anisotropic characteristics, which leadto a preferred orientation of their molecules in the vicinity of solid sur-faces. The ability of liquid crystalline materials to form ordered boundarylayers with good load−carrying capacity and outstanding lubricating prop-erties has been widely demonstrated. In order to study the advantages ofimplementing liquid crystals as lubricants, the steady state / time tran-sient isothermal flow of thermotropic / lyotropic, nematic / chiral nematicliquid crystals between two concentric / eccentric cylinders and in planarCouette geometries were studied numerically. To consider the influence ofthe microstructure formation / evolution on the macro−scale attributes ofthe flow, the Leslie−Ericksen and Landau−de Gennes theories were em-ployed. Simplicity of the Leslie−Ericksen theory in capturing the orienta-tional alignment angle of the molecules makes it a viable candidate for mod-elling the flow of flow−aligning nematic liquid crystals. On the other hand,the Landau−de Gennes nematodynamics equations are well suited for pre-dicting texture formation since defects and disclinations are non−singularsolutions of the governing equations. The Landau−de Gennes theory forthe liquid crystalline microstructure along with continuity and momentumequations were solved simultaneously using General PDE and Laminar Flowmodules of COMSOL Multiphysics. The investigation of flow characteristicsand orientation of liquid crystalline molecules for different rotational veloc-ities / shear rates and anchoring angles at the boundaries were presented.Furthermore, nucleation and evolution of singularities in texture of the liq-uid crystalline materials were tracked over the simulation time. Moreover,alterations in the macro−scale attributes of the flow such as velocity pro-file, pressure distribution and first normal stress difference along with theevolution of defects were studied inside the liquid crystalline domain.The implementation of Landau−de Gennes nematodynamic governing equa-tions for LCs flow simulations offered an insight in application of these ma-terials as lubricants. It was shown the LCs could provide protection againstthe wearing mechanism by forming a shielding layer in the vicinity of solidsurfaces. Three−dimensional simulations of a simplified prosthetic hip jointiiAbstractsuggested that liquid crystalline materials should be considered as potentialbio−lubricants.iiiPrefaceThe following desertion is the result of my graduate study performed atthe Mechanical Engineering Department, University of British Columbiaduring the period of September 2008−September 2013. All of the workpresented henceforth was conducted in the Advanced Numerical SimulationLab (ANS), Institute for Computing, Information and Cognitive Systems atthe University of British Columbia, Point Grey campus.A version of Chapter 2 has been published online [Noroozi N, GrecovD. Flow modelling and rheological characterization of nematic liquid crys-tals between concentric cylinders. Liquid Crystals 40:871-883, 2013] andpresented [Flow modelling and rheological characterization of nematic liq-uid crystals between two concentric cylinders. 23rd Canadian Congress ofApplied Mechanics, 2011, Vancouver, BC, Canada]. I was the lead investi-gator, responsible for all major areas of concept formation, data collectionand analysis, as well as manuscript composition. Dr. Dana Grecov was thesupervisory author on this project and was involved throughout the projectin concept formation and manuscript composition.A version of Chapter 3 has been accepted to the Liquid Crystal journal[Noroozi N, Grecov D, Shafiei-Sabet S. Estimation of Viscosity Coefficientsand Rheological Functions of NanoCrystalline Cellulose Aqueous Suspen-sions, June 2013]. I was the lead investigator, responsible for all majorareas of concept formation, data collection and analysis, as well as the ma-jority of manuscript composition. Shafiei-Sabet S. was involved in the earlystages of concept formation, experimental data collection and contributedto manuscript edits. Dr. Dana Grecov was the supervisory author on thisproject and was involved throughout the project in concept formation andmanuscript edits.A version of Chapter 4 was presented as [Numerical simulation of flowinduced structure in nematic liquid crystalline materials between two cylin-ders, 24th International Liquid Crystal Conference, Mainz, Germany, ILCC2012, 5551-0386]. I was the lead investigator, responsible for all major areasof concept formation, data collection and analysis, as well as manuscriptcomposition. Dr. Dana Grecov was the supervisory author on this projectivPrefaceand was involved throughout the project in concept formation and manuscriptcomposition.I was the lead investigator for the projects located in Chapters 2 to4 where I was responsible for all major areas of concept formation, datacollection and analysis, as well as the majority of manuscript composition.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . 61.3.1 The Leslie−Ericksen Theory . . . . . . . . . . . . . . 71.3.2 The Leslie Viscosity Coefficients . . . . . . . . . . . . 81.3.3 The Landau−de Gennes Nematodynamics GoverningEquations . . . . . . . . . . . . . . . . . . . . . . . . 101.3.4 Introduction to COMSOL Multiphysics . . . . . . . . 121.3.5 Application of Liquid Crystals as Lubricant . . . . . 131.4 Summary of the Gaps in Existing Knowledge . . . . . . . . . 151.5 Objectives and Approach . . . . . . . . . . . . . . . . . . . . 161.6 Thesis Layout . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Flow Modelling and Rheological Characterization of Ne-matic LCs Using Leslie−Ericksen Theory . . . . . . . . . . . 202.1 Governing Equations of the Leslie-Ericksen Theory . . . . . 21viTable of Contents2.2 Numerical Setup and Validation . . . . . . . . . . . . . . . . 242.3 Experimental Method . . . . . . . . . . . . . . . . . . . . . . 272.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4.1 Multiplicity and Multistability . . . . . . . . . . . . . 272.4.2 Viscosity and Torque . . . . . . . . . . . . . . . . . . 382.5 Summary and Remarks . . . . . . . . . . . . . . . . . . . . . 423 Estimation of Viscosity Coefficients and Rheological Func-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.1 Nematodynamic Governing Equations . . . . . . . . . . . . . 433.2 Estimation of Rheological Functions and Landau ViscosityCoefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.3 Experimental Method . . . . . . . . . . . . . . . . . . . . . . 503.4 Calculation of Rheological Material Functions and LandauViscosity Coefficients for Aqueous Suspensions of NCC . . . 513.5 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . 553.6 Selected Results . . . . . . . . . . . . . . . . . . . . . . . . . 563.7 Summary and Remarks . . . . . . . . . . . . . . . . . . . . . 604 Landau−de Gennes Nematodynamics Theory . . . . . . . . 614.1 Theory and Governing Equations . . . . . . . . . . . . . . . 614.2 Computational Methods . . . . . . . . . . . . . . . . . . . . 644.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . 674.3.1 Start−up Flow . . . . . . . . . . . . . . . . . . . . . . 684.3.2 Effect of Ren on Microstructure−Flow Properties . . 824.3.3 Effect of Energy Ratio on Microstructure−Flow Prop-erties . . . . . . . . . . . . . . . . . . . . . . . . . . . 904.4 3−D Analysis of NLC Flow in a Simplified Prosthetic HipJoint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 944.4.1 Numerical Methods for 3−Dimensional Model . . . . 944.4.2 Selected Results for 3−Dimensional Model . . . . . . 964.5 Conclusion and Summary . . . . . . . . . . . . . . . . . . . . 1035 Conclusion, Contributions and Future Research Directions 1055.1 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . 1055.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 1085.3 Future Research Directions . . . . . . . . . . . . . . . . . . . 109Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111viiTable of ContentsAppendicesA Leslie−Ericksen Derivations . . . . . . . . . . . . . . . . . . . 119B Landau−de Gennes Derivations . . . . . . . . . . . . . . . . . 121viiiList of Tables2.1 List of dimensionless variables . . . . . . . . . . . . . . . . . . 242.2 Alignment angles for the LCs used in this study . . . . . . . . 252.3 List of the Leslie viscosities and tumbling factors for the LCsused in this study . . . . . . . . . . . . . . . . . . . . . . . . . 252.4 Frank elasticity constants . . . . . . . . . . . . . . . . . . . . 252.5 The anchoring angles of molecules on the inner and outercylinders for the boundary conditions used in this study . . . 282.6 Frank elasticity and rate of entropy generation for negative /positive initial distribution of molecules . . . . . . . . . . . . 383.1 Coefficients ai1, ai2 and ai3 used to calculate Leslie viscosities 483.2 List of all characteristic parametric values for 7 wt% NCCaqueous suspension implemented in the Landau−de Gennesmodel and extra−stress tensor for the start up shear flowsimulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56ixList of Figures1.1 Schematic representation of transition from solid crystal toliquid crystalline meso−phases . . . . . . . . . . . . . . . . . 21.2 Schematic representation of various configurations of LC moleculesbased on the molecular arrangement . . . . . . . . . . . . . . 21.3 An experiment by Oleg D. Lavrentovich from Kent Universityon defect patterns in a nematic LCs . . . . . . . . . . . . . . 41.4 Orientation of liquid crystalline molecules in vicinity of defectpoints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.5 Dynamic texture of 5 wt% nanocrystalline cellulose aqueoussuspension under steady shear rate a) 0.01 b) 0.05 by Shafieiet al. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Schematic representation of LC molecules with respect to thedirector n for (a) disk-like (b) rod-like, (c) schematic rep-resentation of director vector of liquid crystalline moleculeswith respect to velocity vector . . . . . . . . . . . . . . . . . . 222.2 Principal deformation of liquid crystal molecules, (a) bend,(b) twist, (c) splay . . . . . . . . . . . . . . . . . . . . . . . . 232.3 Distribution of θ vs. ω˜ for (a) DDA9 and (b) AZA9; initialguess for θ was set to θ = pi/10 and for boundary conditionsan anchoring angle of 0.0 (rad) was implemented . . . . . . . 262.4 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = pi/10 and for boundary conditions case B.C. 1was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 292.5 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = pi/10 and for boundary conditions case B.C. 2was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 30xList of Figures2.6 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = pi/10 and for boundary conditions case B.C. 3was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 312.7 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = pi/10 and for boundary conditions case B.C. 4was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 322.8 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = pi/10 and for boundary conditions case B.C. 5was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 332.9 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = −pi/10 and for boundary conditions case B.C. 1was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 342.10 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = −pi/10 and for boundary conditions case B.C. 2was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 352.11 Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ ofinner cylinder for (a) MBBA and (b) PAA; initial distributionwas set to θ = −pi/10 and for boundary conditions case B.C. 3was implemented . . . . . . . . . . . . . . . . . . . . . . . . . 362.12 Direction of change in orientational angle of liquid crystallinemolecules from unstable alignment angle (UAA) to stablealignment angle (SAA) . . . . . . . . . . . . . . . . . . . . . . 372.13 (a) Apparent viscosity of MBBA for different boundary con-ditions (b) Comparison between apparent viscosity of MBBAderived from the Leslie-Ericksen theory and the measured ap-parent viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . 40xiList of Figures2.14 (a) The torque on the inner cylinder measured for MBBA andfrom numerical simulations for B.C. 1 as a function of shearrate (b) Ratio between the computed torque for MBBA andthe computed torque for the hypothetical Newtonian fluid asa function of shear rate . . . . . . . . . . . . . . . . . . . . . 413.1 The Miesowicz viscosities, corresponding to the orientation ofthe director along one of the three axis ηa, for the directoralong the flow direction, ηb for the director along the velocitygradient direction, and ηc for the vorticity axis . . . . . . . . 493.2 Dimensionless Leslie viscosities α˜i as a function of dimension-less concentration cc? for NCC aqueous suspension . . . . . . 513.3 Dimensionless Miesowicz viscosities η˜i, and alignment viscos-ity, η˜al vs. dimensionless concentration cc? for NCC aqueoussuspension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.4 Dimensionless Landau viscosity coefficients ν˜i vs. dimension-less concentration cc? for NCC aqueous suspension . . . . . . . 523.5 Reactive parameter, rotational diffusivity vs. C of NCC aque-ous suspension for different values of the numerical correctionparameter, B . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.6 Characteristic viscosity η¯ vs. C of NCC aqueous suspensionfor different values of the numerical correction parameter, B . 533.7 Alignment viscosity of 5, 7 and 10 wt% NCC aqueous suspen-sion, from experimental and analytical calculations for differ-ent values of the numerical correction factor, B . . . . . . . . 543.8 Shear stress tensor and apparent viscosity as a function ofshear rate, from numerical simulations and from rheologicalexperiments for 7 wt% NCC aqueous suspension . . . . . . . 573.9 Scalar order parameter S at t = 100 [s] for shear rate γ˙ a) 1b) 5 and c) 10[s−1]. . . . . . . . . . . . . . . . . . . . . . . 583.10 Schematic representation of NCC molecules at t = 100 [s] forshear rate γ˙ (top) 1 (middle) 5 and (bottom) 10[s−1]. . . . 594.1 Numerical domain and boundary conditions . . . . . . . . . . 654.2 Orientation angle of MBBA as a function of the dimensionlessgap between the inner and outer cylinder at the centreline,Ren=1.4 × 10−4, De=1.3 and Er=2.4 × 104 for different di-mensionless times . . . . . . . . . . . . . . . . . . . . . . . . . 664.3 Molecular orientation of MBBA, Ren=1.4 × 10−4, De=1.3,Er=2.4× 104, t=105 at the centreline for t=105 . . . . . . . . 67xiiList of Figures4.4 Flow streamlines for De=10, Er=108 and Ren=10−2 at di-mensionless times: a) t=50, b) t=100, c) t=500 and d) t=1000 694.5 Flow streamlines for De=10, Er=108 and Ren=102 at dimen-sionless times: a) t=50, b) t=100, c) t=500 and d) t=1000 . . 694.6 Dimensionless pressure contours forDe=10, Er=108 andRen=10−2at dimensionless times: a) t=50, b) t=100, c) t=500 and d)t=1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.7 Dimensionless pressure contours forDe=10, Er=108 andRen=102at dimensionless times: a) t=50, b) t=100, c) t=500 and d)t=1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.8 Dimensionless first normal stress difference N?1 vs. dimen-sionless rotational angle along the inner cylinder for variousdimensionless times, De=10, Er=108 and Ren=10−2 (left),for Newtonian fluid, Re=10−3 (right) . . . . . . . . . . . . . . 724.9 Dimensionless first normal stress difference N?1 vs. dimen-sionless rotational angle along the inner cylinder for variousdimensionless times, De=10, Er=108 and Ren=102 (left), forNewtonian fluid, Re=10 (right) . . . . . . . . . . . . . . . . . 724.10 Dimensionless pressure p? vs. dimensionless rotational an-gle along the inner cylinder for various dimensionless times,De=10, Er=108 and Ren=10−2 (left), for Newtonian fluid,Re=10−3 (right) . . . . . . . . . . . . . . . . . . . . . . . . . 734.11 Dimensionless pressure p? vs. dimensionless rotational an-gle along the inner cylinder for various dimensionless times,De=10, Er=108 and Ren=102 (left), for Newtonian fluid,Re=10 (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.12 Dimensionless wall shear stress τ ?w vs. dimensionless rota-tional angle along the inner cylinder for various dimensionlesstimes, De=10, Er=108 and Ren=10−2 (left), for Newtonianfluid, Re=10−3 (right) . . . . . . . . . . . . . . . . . . . . . . 744.13 Dimensionless wall shear stress τ ?w vs. dimensionless rota-tional angle along the inner cylinder for various dimensionlesstimes, De=10, Er=108 and Ren=102 (left), for Newtonianfluid, Re=10 (right) . . . . . . . . . . . . . . . . . . . . . . . 744.14 Distribution of scalar order parameter S between two eccen-tric cylinders for De=10, Er=108 and Ren=10−2 for dimen-sionless time s) t=1, b) t=10, c) t=50, d) t=100, e) t=500,and f) t=1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . 76xiiiList of Figures4.15 Distribution of scalar order parameter S between two eccen-tric cylinders for De=10, Er=108 and Ren=102 for dimen-sionless times a) t=1, b) t=10, c) t=50, d) t=100, e) t=500,and f) t=1000 . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.16 Magnified snapshots of the ellipsoidal representation of theLC molecules for De=10, Er=108 and Ren=10−2 for dimen-sionless times a) t=50, b) t=100, c) t=500, d) t=1000 andschematic representation of defect points with strength Sd=± 12 804.17 Magnified snapshots of the ellipsoidal representation of theLC molecules for De=10, Er=108 and Ren=102 for dimen-sionless times a) t=50, b) t=100, c) t=500, d) t=1000 andschematic representation of defect points with strength Sd=± 12 814.18 Dimensionless first normal stress difference N?1 vs. dimen-sionless rotational angle along the inner cylinder for variousdimensionless Reynolds numbers, De=10, Er=108 at t=1000 824.19 Dimensionless pressure p? vs. dimensionless rotational anglealong the inner cylinder for various dimensionless Reynoldsnumbers, Ren for De=10, Er=108 at t=1000 . . . . . . . . . 834.20 Dimensionless wall shear stress τ ?w vs. dimensionless rota-tional angle along the inner cylinder for various dimensionlessReynolds numbers, Ren for De=10, Er=108 at t=1000 . . . . 844.21 Scalar order parameter contours forDe=10, Er=108 at t=100for various Reynolds numbers a) Ren=10−4, b) Ren=10−2, c)Ren=10−1, d) Ren=1, e) Ren=10, f) Ren=102 . . . . . . . . 864.22 Scalar order parameter contours forDe=10, Er=108 at t=1000for various Reynolds numbers a) Ren=10−4, b) Ren=10−2, c)Ren=10−1, d) Ren=1, e) Ren=10, f) Ren=102 . . . . . . . . 874.23 Magnified snapshots of the ellipsoidal representation of theLC molecules, De=10, Er=108 at t=100 vs. Reynolds num-bers a) Ren=10−4, b) Ren=10−2, c) Ren=10−1, d) Ren=1,e) Ren=10, f) Ren=102 . . . . . . . . . . . . . . . . . . . . . 884.24 Magnified snapshots of the ellipsoidal representation of theLC molecules, De=10, Er=108 at t=1000 vs. Reynolds num-bers a) Ren=10−4, b) Ren=10−2, c) Ren=10−1, d) Ren=1,e) Ren=10, f) Ren=102 . . . . . . . . . . . . . . . . . . . . . 894.25 Flow streamlines at t=150, a) R=107, Ren=10−2, b) R=105,Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1 . . . . . 904.26 Dimensionless pressure contours at t=150, a)R=107, Ren=10−2,b) R=105, Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1 91xivList of Figures4.27 Scalar order parameter at t=150, a) R=107, Ren=10−2, b)R=105, Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1 914.28 Schematic representation of molecules at t=150, a) R=107,Ren=10−2, b) R=105, Ren=10−2, c) R=107, Ren=1 and d)R=105, Ren=1 . . . . . . . . . . . . . . . . . . . . . . . . . . 934.29 Schematic representation of the capsular space in prosthetichip joint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 954.30 Contour slices of the dimensionless pressure for a) lyotropicliquid crystal with Er=9.756 × 105, De=10 and Ren=1 b)synovial fluid diagnosed with RA, Re=0.1 . . . . . . . . . . . 974.31 Contour slices of the dimensionless first normal stress dif-ference for a) lyotropic liquid crystal with Er=9.756 × 105,De=10 and Ren=1 b) synovial fluid diagnosed with RA, Re=0.1 984.32 Contours of the dimensionless wall shear stress on the head ofartificial femur for a) lyotropic liquid crystal with Er=9.756×105, De=10 and Ren=1 b) synovial fluid affected by RA,Re=0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 994.33 Contours slices of the scalar order parameter, Er=9.756×105,De=10 and Ren=1 at dimensionless times a) t=10 b) t=1000 1004.34 Schematic of three−axes ellipsoids representing liquid crys-talline molecules . . . . . . . . . . . . . . . . . . . . . . . . . 1024.35 A magnified slice of schematic representation of liquid crys-talline molecules, a) (left) narrow gap (right) wide gap att=10; b) (left) narrow gap (right) wide gap at t=1000; Ren=1,Er=9.756× 105 and De=10 . . . . . . . . . . . . . . . . . . . 102xvNomenclatureG External director body force [N ]N Angular velocity of director vector with respect to velocity field [s−1]A Rate of deformation tensor [s−1]W Vorticity tensor [s−1]Q Order parameter tensorF Flow contributionsHsr Short-range elasticity contributionsHlr Long-range elasticity contributionsv Velocity [m/s]f Body force [N ]t Leslie-Ericksen total stress tensor [Pa]g Intrinsic director body force [N ]n Eigenvector of Q corresponds to largest eigenvalue, director vectorm Eigenvector of Q corresponds to second largest eigenvaluel Eigenvector of Q corresponds to smallest eigenvalueI Second order unity tensorp Pressure [Pa]s˙ Rate of entropy generation [J/K.kg.s]Sd Defect strengthFd Frank elasticity free energy function [Pa]kii Frank elasticity constants [N ]T Resistant torque [N.m]S First scalar order parameterP Second scalar order parameterU Nematic potentialkB Boltzmann constant [J/K]T Temperature [K]H Flow characteristics length scale [m]L1 First Landau elastic coefficient [N ]L Isotropic long-range elasticity length scale [m]L? Non-dimensional Landau elastic coefficientP2 Second order Legendre’s polynomialP4 Fourth order Legendre’s polynomialxviNomenclaturel Length of LC molecules [m]d Diameter of LC molecules [m]r Variable radius [m]R1 Inner radius [m]R2 Outer radius [m]ε Eccentricity [m]C Dimensionless concentrationc Molecular concentration [m−3]Seq Equilibrium scalar order parameterDr Rotational diffusivity [s−1]De Deborah numberEr Ericksen numberR Energy RatioGreek LettersP Chirality contributionsB Numerical correction parameterpi Intrinsic director surface stress tensor [Pa]τ t Landau-de Gennes’ total stress tensor [Pa]τ a Landau-de Gennes’ asymmetric stress tensor [Pa]τ s Landau-de Gennes’ symmetric stress tensor [Pa]τ v Landau-de Gennes’ viscous stress tensor [Pa]τ e Landau-de Gennes’ elastic stress tensor [Pa]τEr Landau-de Gennes’ Ericksen stress tensor [Pa]τw Wall shear stressλ Reactive parameterρ Density [kg/m3]ρ1 Material constant [kg/m]µi Leslie viscosity coefficients [Pa.s]θ Orientational angle of molecular directorθ0 Anchoring angle of the molecules on the wall boundariesω Rotational velocity [s−1]β Material shape factorΘ Ratio of anisotropic distortional elasticity to the isotropic one Levi-Civita symbolγ˙ Shear rate [s−1]L Anisotropic long-range elasticity length scale [m]νi Landau viscosity coefficients [Pa.s]ηs Solvant viscosity [Pa.s]xviiNomenclatureγ Molecular aspect ratio correction parameterηal Alignement viscosity [Pa.s]ηa Miesowicz viscosity along the flow direction [Pa.s]ηb Miesowicz viscosity along velocity gradient [Pa.s]ηc Miesowicz viscosity along vorticity axis [Pa.s]AbbreviationsMUMPS MUltifrontal Massively Parallel sparse direct SolverFEA Finite Element AnalysisCOF Coefficient Of FrictionSEM Scanning Electron MicroscopyROII Relative Optical Interference IntensityRCD Reciprocating Cylinder on DiscCFD Computational Fluid DynamicsRA Rheumatoid ArthritisPAA Para−AzoxyanisoleDDA9 Poly(4, 4′−dioxy−2, 2′−dimethylazoxybenzene−dodeccanediyl)AZA9 22−dimethyl−3−hydroxy−8−methyl−azaspiracidMBBA N−(4−Methoxybenzylidene)− 4−butylanilinexviiiAcknowledgementsI would like to express my eternal gratitude to my supervisor, AssociateProf. Dana Grecov for her influential guidance, encouragement, patienceand tremendous concern for the success of this research. Additionally, Iwould like to acknowledge my supervisory committee members, ProfessorIan Frigaard and Professor Savvas G. Hatzikiriakos. Also, I’d like to expressmy gratitude to Dr. Hatzikiriakos for many helpful discussions and guidancethroughout the preparation of the thesis.xixDedicationI dedicate my dissertation work to my family, friends. A special feeling of gratitudeto my loving parents, Saeid and Homeira Noroozi and my sister Nazbanoo forbelieving in me. I also dedicate this dissertation to my many friends who havesupported me throughout this process especially Farhang, Arash, Meysam, Farzad,Mehdi, Faramarz and Shervin. I dedicate this work and give special thanks to myfianc´e Azin Jahanafrooz for nonstop support and for the many hours of proofreading.Without the moral support of any of you this work would never reach the end.xxChapter 1Introduction1.1 BackgroundIn 1888, Austrian botanist Friedrich Reinitzer encountered a new phaseof matter while performing physicochemical experiments on derivatives ofcholesterol. Although at that time he was not confident about his discov-ery, he reported his observations to one of his colleagues, Otto Lehmann.Reinitzer observed a reversible phenomenon during the process of warmingcholesterol benzoate; a new phase with distinct cloudy colour had occurredbetween 145.5◦C and 178.5◦C. This intermediate phase was also observedin other materials; for example, although his samples drifted as liquids, theyexhibited crystallinity. It was the peculiar properties of these meso−phasesthat led to the term Liquid Crystals [1]. Although not much progress wasoriginally done following the discovery of the liquid crystalline state andof materials with liquid crystalline characteristics, more recent times haveseen an exponential increase in the interest in studying static and dynamicbehaviour of the liquid crystalline materials.In short, the liquid crystal (LC) phase is a state of matter that hasproperties between those of a conventional liquid and a solid crystal. LCsdo not inherently possess positional order but they keep the orientationalorder of solids and the flowing properties of liquids while exhibiting stronganisotropic properties. A schematic representation of the transition fromthe solid crystal phase to the liquid crystalline meso−phase is presented infigure 1.1.LC molecules can be classified into several different categories based onorigin, shape, molecular arrangement and also their formation process. Theyare classified as organic and non−organic in terms of origin, and rod−likeand disk−like in terms of shape. Furthermore, based on their arrangement,LCs are further categorized into nematic, cholesteric, smectic A, smectic C,columnar etc. [2]. A schematic representation of the arrangement of LCmolecules is presented in figure 1.2.11.1. BackgroundCrystal Liquid Crystal (Mesophases) LiquidTemperature or Concentration- Anisotropic- 3D Lattice- Positional Order- Orientational Order- Anisotropic- No Lattice- No Positional Order- Orientational Order- Isotropic- No Lattice- No Positional Order- No Orientational Order- Anisotropic- 1 or 2D Lattice- No Positional Order- Orientational OrderFigure 1.1: Schematic representation of transition from solid crystal to liquidcrystalline meso−phasesNematic Smectic A Smectic CColumnar CholestericFigure 1.2: Schematic representation of various configurations of LCmolecules based on the molecular arrangementNematic liquid crystal (NLCs) is the simplest liquid crystalline phase andalthough it was generally assumed that its properties were uniaxial [3], ne-matic LCs with biaxial properties have been observed [4]. In fact, cholesteric21.1. Backgroundmeso−phases can be considered a specific type of nematic LCs that have atwist in their molecular orientation about an axis normal to the preferredorientation of the molecules. As mentioned above, nematic and cholestericliquid crystalline materials posses orientational order but no translationalorder, however, as smectic LCs maintain a degree of translational order,they stratify a vast number of structures.Smectic C demonstrates a tilted molecular orientation at different layers,contrary to the upright molecular orientation of smectic A.Several different classes of discotic liquid crystalline materials have been in-troduced based on the orientation of the molecular director and degree oftranslational order (columnar LCs is only one example of these classifica-tions). However, the main focus of this thesis was limited to the rod−likeLC molecules and as such, discotic / disk−like liquid crystalline moleculesare not further discussed in this document.In addition to the above classes, liquid crystalline materials are furthercategorized according to the process that leads to the formation of the liq-uid crystalline meso−phase, namely thermotropics and lyotropics. Withthermotropic liquid crystalline materials, changes in temperature affect themolecular motion and the level of free energy, leading to a phase transition.By melting a crystalline solid or cooling an isotropic liquid to a specifictemperature, a thermotropic phase that exists only within a limited tem-perature range can be reached. In contrast, lyotropic LCs are obtained bychanging the concentration of a particular amphiphilic substance in a sol-vent such as water. Although generally, particle shape, size and charge playthe primary role in the phase transition from an isotropic mixture to a ly-otropic suspension / solution, temperature remains an important factor inthis transition.Frank [5] proposed the term ”disinclination” to describe unique patternsof liquid crystalline materials under polarized light microscopy in analogywith dislocations in solid crystalline materials; this term was later changedto ”disclination”. Any local breakage in symmetric patterns of the orienta-tion of liquid crystalline molecules can be described as a defect. In staticliquid crystalline samples between two glass plates under polarized light mi-croscopy, defect points and brushes around them can be easily spotted sincethe position of defects remain constant while the orientation of disclinationlines (brushes) changes with respect to the polarized light angle. Figure1.3 depicts experimental results on defect patterns in a nematic LC sample(Oleg D. Lavrentovich, Kent University).According to Chandrasekhar [6], defect strength can be calculated byknowing the number of disclination brushes attached to the defect point,31.1. BackgroundFigure 1.3: An experiment by Oleg D. Lavrentovich from Kent Universityon defect patterns in a nematic LCsSd = 14 (# of brushes). Defects with Sd = ±12 and Sd = ±1 are mostly ob-served in the literature [6]. The importance of these patterns has been knownsince the discovery of liquid crystalline states by Reinitzer and Lehmann [7],however, the mathematical study and formulation of these textures wereconducted many years later by Oseen [8] and Frank [5] for static samplesof liquid crystalline molecules. Frank suggested various types of defects anddisclination lines according to the definition of defect strength, presented infigure 1.4.Additional studies were conducted on the dynamic patterns of LCs flow-ing in encapsulated domains. Figure 1.5 depicts the liquid crystalline pat-terns under the flow−induced structure for 5 wt% nanocrystalline celluloseaqueous suspension under steady shear rate of 0.01 and 0.05 [s−1] [9]. In thisfigure, cross polarized light microscopy has been implemented to capture theisotropic and anisotropic subdomains.41.2. Motivation21!"s 21" 4,1 #"" c2,1 #"" c 23" 2"1!"d sd sd sd sdsd sd sdFigure 1.4: Orientation of liquid crystalline molecules in vicinity of defectpoints, by Frank [5] (figure extracted from [6])Figure 1.5: Dynamic texture of 5 wt% nanocrystalline cellulose aqueoussuspension under steady shear rate a) 0.01 b) 0.05 [s−1] by Shafiei et al. [9]1.2 MotivationLiquid crystalline materials have a potential to be used in different appli-cations. Structurally, liquid crystalline materials can be used as precursorsin the fabrication of films, fibers and foams and their transformation intohigh-performance materials using conventional processing operations. Dif-ferent studies have been conducted on the application of liquid crystallinemicrostructure as precursors and the influence of rheological properties suchas viscoelastic modules, anisotropic viscosity and flow−induced structure on51.3. Literature Reviewthe processing of liquid crystalline materials [10]. However, this study was fo-cused on the functionality of LC meso−phases, based on their flow−inducedtexture evolution.In terms of function, LC’s have also been used as lubricating agents [11];for example, one of the functional applications of rod−like liquid crystallinematerials is implementing them as lubricants or as performance boosting ad-ditives to standard lubricants, which take advantage of the unique anisotropicrheological characteristics of LCs. The intrinsic properties of LCs allow themto enhance their orientation−dependent viscosity and customize it for anyapplication. Anisotropic characteristics of liquid crystalline materials alsolead to complex flow properties; thus, several different attempts have beenmade to propose a general continuum model in order to predict the flowingcharacteristics of LC materials.The subject of employing liquid crystalline materials as lubricants, or as ad-ditives to other standard lubricants to enhance their performance in journalbearings, requires an examination of their flowing characteristics, which wasthe main motivation for this study. Moreover, overcoming the challengesraised by simulating the flow between two eccentric cylinders (as a simpli-fied model for journal bearings) due to the existence of shear, rotational andextensional forces, provided additional incentives for this study.1.3 Literature ReviewIn the previous section the role of liquid crystalline meso−phases as lu-bricating agent have been stated as the main motivation for this study;unfortunately, the complex anisotropic characteristics of liquid crystallinematerials made their flowing properties intensely complicated, which limitsthe range of predictions for applying LCs as lubricants.Flow of liquid crystalline materials contains a unique micro−scale texturethat evolves in time and space and interacts with macro−scale flow charac-teristics. Investigating the formation and evolution of LC’s textures (liquidcrystalline finger−prints) was the objective of many studies; several differentcontinuum models were suggested trying to capture all different aspects ofthis phenomenon.This study was focused on illustrating various aspects of the complicatedflow characteristics of liquid crystalline materials based on the different con-tinuum models to have a better understanding from flow of LC materialsas lubricants. In the following sub−sections, details of other researchers’accomplishments are presented; starting by the achievements on the contin-61.3. Literature Reviewuum model, which was proposed by Frank Matthews Leslie and Jerald Erick-sen; followed by the contributions of other researchers on the more complexcontinuum model that was proposed by Nobel prize winner Pierre−Gillesde Gennes and Lev Davidovich Landau. At the end, application of liquidcrystalline meso−phases as lubricants were discussed.1.3.1 The Leslie−Ericksen TheoryLCs are anisotropic viscoelastic materials; the combination of fluid−like flowwith crystal−like anisotropy makes their phases interesting as modifiers ofinterfacial behaviour when applied as lubricants. The anisotropy of theviscosity coefficient, with respect to different flow directions, is a uniqueproperty of the liquid crystalline phase. The ability of liquid crystalline ma-terials to form ordered boundary layers with good load−carrying capacity,and to decrease the friction coefficients, wear rates and contact temperaturesof sliding surfaces (contributing to increase the components service life andsave energy) has been widely demonstrated [11–13]. Several authors havepresented experimental, theoretical and numerical results for the tribologyof LCs [11]. The effect of using liquid crystal as additive for mineral oilon the coefficient of friction, which resulted in major improvements, wasinvestigated in [14]. Liquid crystals used as smart lubricants with an activecontrol of friction coefficient, based on electrical and / or hydrodynamic fieldvariation has been studied in [15].Leslie and Ericksen [16, 17] formulated the conservation laws of mass,linear and angular momentum for nematic LCs using a director vector n.Atkin and Leslie [18] simplified the governing equations in the case of simpleCouette flow of nematic LCs based on the assumption of constant Frankelasticity coefficients. Measurements of splay (k11) and bend (k33) Frankelasticity coefficients have been performed on nematic LCs by Bradshaw etal. [19].Rey [20] revealed that in two−dimensional modelling of liquid crystalsbased on the Leslie−Ericksen model, the director escapes from the shearplane under particular conditions. Rey studied transition of the director ofnematic LCs in radial flow and discovered that the out−plane orientationof the director was due to the high Ericksen number (the elastic and shearforces were incapable of keeping the director inside the shear plane). In 2000,Van Horn et al. [21] implemented a simplified version of the Leslie−Ericksenmodel to solve the start−up simple shear flow; Amar and Cummings [22]used the same Leslie−Ericksen model to study instabilities in thin layers ofnematic LCs on rigid surfaces; however, classification of results based on the71.3. Literature ReviewLeslie−Ericksen model was performed by Lima and Rey [23]. They studiedmultiplicity and multistability of solutions of the Leslie−Ericksen model forPoiseuille flow of discotic LCs.1.3.2 The Leslie Viscosity CoefficientsLiquid crystals are artificial and biological anisotropic viscoelastic texturedmaterial. Due to their properties, they can be used as functional materi-als, sensors, actuators, lubricants or as precursor materials in the manufac-turing of fibers, films, foams or composites. Their flow behaviour can bedescribed using the Leslie−Ericksen director vector (n) theory [8, 16, 17],or the Landau−de Gennes tensor order parameter (Q) theory. The uni-axial director n represents the average orientation of the unit axis of therod−like molecules. The degree of order for the unit normals along nis given by the scalar order parameter (S). The orientation variable is aslow long wave−length mode and the scalar order parameter is a fast shortwave−length mode; thus, the selection of the theory depends on the purposeand use of the model. For fast flows and / or descriptions of short lengthscale phenomena such as defect nucleation, the Landau−de Gennes modelmust be used. The Leslie−Ericksen model is adequate for sufficiently slowflows, when the flow time scale slower than the internal one, in which theorientation dominates the rheology and the scalar order parameter reachesclose to its equilibrium value, and long length scale phenomena. In the limitof slow flows the Landau−de Gennes theory converges to the Leslie−Ericksentheory [24], and thus the viscoelastic material parameters of the latter canbe given in terms of parameters of the former.The Leslie−Ericksen model incorporates a modified stress tensor in theNavier−Stokes equations to couple the microstructure with macro−scaleattributes of the flow of Liquid crystalline materials [18]; therefore, six newcoefficients have been introduced in the stress tensor.The evolution of order parameter tensor Q in time has been the principalfocus of the Landau−de Gennes theory in predicting complex patterns in theflow of LCs [25, 26]. In comparable manner to the Leslie−Ericksen model,the Landau−de Gennes nematodynamic theory introduces an asymmetricstress tensor to capture the complicated phenomenological aspects of theflow of liquid crystals. Total stress tensor for the Landau−de Gennes modelhas a viscous component that implements three viscosity coefficients (νi)to model the anisotropy of rheological characteristics of LCs [25]. Grecovproposed a method for calculation of νi as a function of the Leslie viscositiesand scalar order parameter by mapping the two stress tensors [24].81.3. Literature ReviewAnisotropic rheological characteristics of the liquid crystalline materialsare well described by the Leslie viscosity coefficients αis, i = 1, . . . , 6; whereonly five of these are independent due to the Parodi’s relation [27]. Marrucci[28] proposed a formulation to predict the values of Leslie coefficients assum-ing very thin, long and rigid rod−like molecules, and considering α6 = 0.Moreover, he correlated the Leslie coefficients only to the scalar order param-eter. Doi [29] implemented Smoluchowski equation and equilibrium distri-bution of molecular unit vectors around the director. Furthermore, secondand fourth order Legendre polynomials have been applied as 2nd and 4thmoment of the bulk of molecules to express the Leslie viscosity coefficientsof lyotropic liquid crystals with rod−like polymeric molecules under weakvelocity gradients [30]. Larson proposed additional viscous contributions tothe Leslie viscosity coefficients using a rod−solvent friction factor ξstr whichis a function of rotational diffusivity of dilute suspension and temperature[31].Nanocrystalline cellulose (NCC) is composed of nano−fibres with highaspect ratio cellulose and unique rheological characteristics [9, 32]. Due tothe rod−like shape, size and surface charge of the particles, NCC aqueoussuspension exhibits liquid crystalline state in its phase diagram [32], whichdisplay the influence of concentration and temperature on the phase changeof NCC aqueous suspensions. Colloidal suspension of NCC molecules inwater are known to possess chirality in their nematic phase [33]. Dynamicsof chiral nematic liquid crystalline materials have been studied using theLandau−de Gennes theory by introducing an extra term in the evolutionof order parameter tensor; the extra term regarding the chiral pitch of thecholesteric liquid crystals is related to the ratio of elasticity to the helicaldistortions of the molecules [34, 35].As been discussed by many authors before, the Landau−de Gennes ne-matodynamic model is well suited to simulate texture formation since de-fects are non singular solutions to the governing equations; however, findingmeaningful values for various parameters and variables in the model is prob-lematic. In this study different methods that have been proposed by otherresearchers were gathered to shed light on the values of different variables ofthe Landau−de Gennes model with focus on capturing trends of the Landaucoefficients in viscous component of the total shear stress (νi) for lyotropicLCs e.g. NCC aqueous suspension.91.3. Literature Review1.3.3 The Landau−de Gennes Nematodynamics GoverningEquationsNematic liquid crystals are orientationally ordered, textured, anisotropicand viscoelastic materials. They are used as structural or functional ma-terials. One example of their various functionalities is their application aslubricants, which is related to the anisotropy of the viscosity coefficient withrespect to different flow directions; this is a unique property of the liquidcrystalline phase. The remarkable rheological properties of these materialsare due to the flow−induced evolution of molecular configurations. Inde-pendent studies indicate that LCs change their orientational state in thevicinity of a solid surface. The ability of liquid crystalline materials to formordered boundary layers with good load−carrying capacity, and to lower thefriction coefficients, wear rates and contact temperatures of sliding surfaces(contributing to increase the components service life and save energy) hasbeen widely demonstrated [12, 36–38].Modelling the flow of liquid crystals is a dynamic multidisciplinary areaof research for experimental and numerical fluid flow researchers as well asapplied mathematicians. Since 1888 that liquid crystals were examined forthe first time by Friedrich Reinitzer, several continuum theories had beendeveloped to capture different aspects of the rheology and structure for-mation in liquid crystals. The flow of liquid crystals has been investigatedthrough numerical simulations, in several studies, using Landau−de Gennestheory, Doi−Marrucci−Greco theory and extensions of Doi theory for ne-matic polymers. In Doi theory, at any given point in space, the dynamics ofrod−like molecules of nematic liquid crystalline polymers (NLCPs) can befound using a diffusion equation of their orientational distribution functionby implementing the Maier−Saupe nematic potential function [39].Start−up flow of lyotropic LCs between two eccentric cylinders has beenstudied using the Doi theory by [40]; defects generation and the couplingbetween the flow and the molecular structure has been investigated usingquadratic closure approximation to estimate the evolution of molecules [41].Six closure models considered for the Doi theory have been compared witheach other to predict three regimes of director motion: steady alignment,wagging and tumbling. All of the closure models foresaw the same dis-tinct attributes of the liquid crystals behaviour. Complex flow simulation ofNLCPs has been performed along with comparison among different closureapproximations of the Doi theory for rigid rod−like molecules by [42].Essential elements in the pattern formation of the LCP was solved usingthe Doi−Marrucci−Greco and moment−averaged theory by [43]. Slow Cou-101.3. Literature Reviewette flow of NLCPs was modelled using the Doi−Marrucci−Greco tensorialtheory by [44], extending the structural scaling laws which was developedbased on the Leslie−Ericksen theory to incorporate the flow of NLCPs andmolecular elasticity. Furthermore, the mesoscopic Doi−Marrucci−Grecotheory and the transient Stokes flow of NLCPs were simulated by [45] tostudy unsteady structural and orientational transition in unsteady shearflow. The texture formation of discotic nematic carbon fibers during thethermal relaxation after the termination of extensional flow was studied byHong and Chan [46]. The Landau−de Gennes free energy of molecules wasused to predict the difference between molecular reorientation time scaleand thermal relaxation [47].The isothermal, incompressible flow of liquid crystalline polymers wasformulated by Farhudi and Rey using tensors for the microstructure; thetransient and steady simple shear flow of uniaxial NLCs was studied. Fur-thermore, the elastic contribution for the molecular field and the elasticstress tensors were neglected to simplify the model and a spatially homoge-nous director vector for the molecules were assumed. Two types of orienta-tional modes were predicted according to the nematic potential, (a) simplealigning mode, (b) complex modes. Three different subdivision regimes ofcomplex modes were further investigated via the change in shear rates, whichresulted in describing two distinct critical shear rates that characterize thetransition between tumbling, oscillatory and stationary regime [48, 49].Tsuji and Rey solved constitutive equations for evolution of LCs in rec-tilinear simple shear flow. Both, short−range elasticity and long−rangeelasticity were considered in their model. The compatibility of this tensorialtheory was tested with the Leslie−Ericksen and Doi theory for asymptoticcases of De→∞ and Er →∞ respectively. A mechanism for compatibilitybetween the Doi’s tumbling effect and fixed anchoring angle of the moleculeson the boundaries was suggested [50] .A generalized version of tensorial governing equation was proposed byTsuji; in these governing equations, short−range elasticity, flow contribu-tions and long−range elasticity were included. Shear flow of liquid crys-talline materials as a function of the Ericksen and Deborah numbers andfor a fixed anchoring angle of the molecules on the boundaries was inves-tigated. In their study, the counter impact of the microstructure on flowcharacteristics of LC was neglected [51].Grecov and Rey adapted the Landau−de Gennes theory to describe theflow of flow−aligning thermotropic LCs. Well known thermodynamic re-strictions along with ordering of Miesowicz viscosities and flow−aligningangle of the director emerging from the Leslie−Ericksen theory were imple-111.3. Literature Reviewmented to characterize shear viscosity coefficients of NLCs [24]. Start−upshear responses for different anchoring angles of the director at the wallboundaries and various Ericksen and Deborah numbers were investigated[25]. Grecov and Rey implemented the Landau−de Gennes equations alongwith decoupling simplification and constant shear rate, to achieve an insighton defects nucleation in shear−induced flow of flow−aligning liquid crys-talline polymers (LCPs). The understanding of defect−defect annihilationand defect−wall interactions resulted from this study was needed to achievea defect−free monodomain melts of LCPs [52, 53].In summary, the Landau−de Gennes theory is well suited to simulatetexture formation since defects are non−singular solutions to the governingequations. However, solving the Landau−de Gennes model for realistic ap-plications is computationally expensive; thus, numerous simplifications ofthis theory have been done in the past to make it simpler to solve. Forexample, neglecting the counter impact of the microstructure on the flowor neglecting the long−range contribution on the shear stress tensors forsolving the fluid flow have been used as simplifying assumptions.1.3.4 Introduction to COMSOL MultiphysicsCOMSOL Multiphysics is a commercial FEA1 simulation package that hasbeen developed to handle complex entangled multiphysics physical phe-nomenon. The very first version of the software was developed by SvanteLittmarck and Farhad Saeidi as FEMLAB in KTH2. COMSOL uses Galerkinfinite element method to discretize the continuum domain into several finitepieces. In Galerkin finite element, the solution is calculated at nodal pointsand interpolated between nodes using the test functions v(x). In this casethe solution for an arbitrary variable u(x, t) can be written as,u (x, t) =∞∑i=0ωi (x)ui (1.1)For the case of solving the Landau-de Gennes equations for the LC’s mi-crostructure coupled with the Navier-Stokes equations for the LC’s flow,the system of partial differential equations can be discretized by the FEMmethod that results in a linear system of equations. For example, in Pois-son’s equation with a Dirichlet boundary condition, weak formulation of the1Finite Element Analysis2Royal Institute of Technology, Stockholm, Sweden121.3. Literature Reviewtest functions can be written as,∫Ω∇u (x) · ∇v (x) ds =∫Ωfv (x) ds (1.2)Which can be transformed to Au = R where A is a large, usually sparsematrix. COMSOL solves this system of equations using MUMPS3. MUMPSuses multifrontal LU decomposition to generate upper triangle and lowertriangle matrices, A = LU that leads to calculating the inverse of matrix A.By knowing A−1, the unknown nodal vectors can be calculated as u = A·R.Moreover, COMSOL employs multi order Backward Differential Formula tocalculate the time steps in any transient analysis and variable damping factorbase on CFL numbers for fluid mechanic related problems.1.3.5 Application of Liquid Crystals as LubricantNature of liquid crystals leads to a preferred direction of molecules in thevicinity of solid surfaces that gives them outstanding tribological proper-ties. It has been proved [11, 12, 37, 38, 54–60] that liquid crystals can beused as lubricants or as additives to industrial lubricants. Studies performedby [14, 55] on the influence of the molecular structure on the tribologicalperformance of liquid crystalline meso−phase showed that not only the ori-entation of rigid part of the molecules but also the alignment of the flexiblecomponents of the molecular structure have contributions on the tribolog-ical characteristics of LCs; which also demonstrated a decrease in COF4by increasing molecular mass of the mixture (mineral oil lubricant (SAE10W/40) + LCs).Influence of the implementation of LCs as additives to other types oflubricants was further investigated by [54]. In this study, structural for-mation in standard lubricants due to the presence of additives with liquidcrystalline nature such as nematics, cholesteric and smectic were resultedin four regimes of relative COF based on increasing the load for differenttypes of additives. Further investigations have been proved that lubricantsenriched with liquid crystalline particles, capable of configuring planar mi-crostructure close to a contact surface, exhibit anti−friction characteristics.After observing the self−aligning pattern of liquid crystalline molecules closeto the solid surfaces using imaging techniques by [60]; optimum value of 1%weight addition of LCs to reach minimum COF and maximum wear resis-tance between two contact surfaces were reported.3MUltifrontal Massively Parallel sparse direct Solver4Coefficient Of Friction131.3. Literature ReviewTribological properties of liquid crystals (LC 1− 4, LC 5− 8, LC 9− 12,LC 13 and LC 14) as additives to common oil based lubricants (HN−95−10and SAE 15 W/40) were explored by Bermu´dez et al. [56]. Tribology studiesusing pin and disk friction tester were performed on aluminum−steel metalcouple which is more difficult to lubricate using common lubricants. Weightloss during the wearing process was measured with an electronic balanceand worn surfaces were separately analyzed by SEM5. Using SEM imagingtechnique also showed the tangential alignment of the molecules’ orientationwith respect to the solid surface.Additionally, in the case of implementing lyotropic LCs as lubricants onaluminum alloy surfaces, higher load carrying capacity was observed com-pared to the standard oil based lubricants [57]. This outcome was achievedby comparing the widths of the wear scars measured under microscope forlyotropic LCs and liquid crystalline−enriched lubricants.Nano−tribological properties of liquid crystals as an additive were stud-ied with ROII6 technique [13]. As been suggested before [61], orienta-tional order of liquid crystalline meso−phase grants them long−range, quasilong−range and short−range elasticity that was discussed in chapter 4.The influence of long−range order of the molecules on film forming mech-anism and thickness of the lubricating film was studied by [13]. Outstand-ing anisotropic rheological characteristics of liquid crystalline materials givethem outstanding advantages over the isotropic liquids while being used aslubricants. Considering the studies performed by [15], not only the rheo-logical properties of LCs make them a suited candidate to be used as lu-bricants but also the influence of electromagnetic fields on liquid crystallinemeso−phases initiate a possibility of producing smart lubricants capable ofbeing actively controlled to modify their viscosity and adjust any changesin operating parameters.Electro−rheological characteristics of the mixture containing liquid crys-talline molecules as additives were investigated, which resulted in a raise inthe magnitude of viscosity and film thickness by increasing the strength ofthe electric field; which led to rearrangement of the molecular orientation.In 2005, dynamic characteristics of liquid crystals as lubricants in slidingbearings were measured under applied electric field [62]. Experimental dataproved that the orientation of liquid crystalline molecules became tangentialto direction of the electric field; the fact that the main axis of liquid crys-talline molecules coincide with the shear direction on lubricated surfaces5Scanning Electron Microscopy6Relative Optical Interference Intensity141.4. Summary of the Gaps in Existing Knowledgeled to a decrease in the value of COF even in the absence of the electricfield. Controlling the viscosity of liquid crystalline meso−phase by alteringthe voltage in an electromagnetic field was achieved in order to reach lessfriction forces.Lately, ultra low friction fluid mixtures using mesomorphic materialswere studied by [12]; SAE 10W−40 and paraffin oil were used as solvents,D1 and 07/10 mesogenics, which are discotic and rod−like LCs, as addi-tives. Time dependent friction coefficient of such mixtures with variouscombination of different concentration of LCs in oil was examined usingRCD7 tribometer.In conclusion, the role of liquid crystalline materials as lubricants or asadditives to other type of lubricants has been established. It also shouldbe taken into account that most LCs have intrinsic presence in living or-ganisms [63], with lowest energy dissipation associated with their orderedlayers; which make them a viable candidate as bio−lubricants. Waters etal. [38] performed experiments on wear of the artificial hip−joint materialsby using normal aqueous saline, aqueous sodium azide, pure ethanol, aque-ous hyaluronan acid and cholesterol palmitate as lubricants in ball and disktribometer. After performing the examination on worn disks with SEM, sig-nificant decrease in wear scars and formation of liquid crystalline protectionshield on the counter−face surfaces were detected.1.4 Summary of the Gaps in Existing KnowledgeThe flow between two eccentric cylinders, used as the basis of the flow modelin journal bearings, is particularly complex in nature since shear, extensionaland rotational flow features are present. When used in combination with liq-uid crystals, the model becomes even more elaborate. In fact, to the best ofour knowledge no numerical studies have been able to tackle the lubricationof liquid crystalline materials. Despite this, a number of phenomenologi-cal continuum models have been proposed and several studies have beenperformed on the simple shear flow of NLCs.As such, the flow of liquid crystalline materials in the complicated ge-ometry of two eccentric cylinders was one of the challenges that needed tobe addressed. Moreover, the Landau−de Gennes theory, which is one ofthe most sophisticated theories in predicting the flow properties of liquidcrystalline materials has an additional intrinsic complexity for numericalsimulation due to the presence of second and fourth order gradients of order7Reciprocating Cylinder on Disc151.5. Objectives and Approachparameter tensor. As the literature review suggests, many studies conductedin the past on the numerical analysis of LC materials using the Landau−deGennes theory present a variety of shortcomings, namely:• The counter−impact of the microstructure on the flow properties havebeen commonly neglected and by consequence, microstructure and flowfield were separately modelled. In other words, the evolution equationand the Navier− Stokes equations were segregated in the solving pro-cess.• Little attention has been paid to the influence of third dimensions ofLC molecular shape and flow. Most of the studies in this field wereperformed on 2−dimensional models, which lead to lower accuracycompared to 3−dimensional models. This is particularly evident intrying to capture specific aspects of liquid crystalline behaviour suchas tumbling or a director escaping from shear plane.• The Landau−de Gennes theory has not been used to address the flowof liquid crystalline materials between two eccentric cylinders, which isinduced by the rotation of the inner cylinder and therefore has inherentshear, rotational and extensional features that result in more complexflow patterns.• The idea of implementing liquid crystalline materials as lubricants /lubricating agents has been explored by many experimental studiesbut has not been investigated by numerical simulations.In order to accurately predict the nucleation and evolution of the LC’stexture, and understand the implications of LC patterns on its rheolog-ical and tribological characteristics (by using LCs as lubricating agents),the unabridged Landau−de Gennes theory coupled with the Navier−Stokesequations has to be implemented to model the start−up flow of NLCs be-tween two eccentric cylinders (as a simplified model for journal bearings).1.5 Objectives and ApproachOur understanding of the macro−scale attributes of the flow of liquid crys-talline materials as lubricants is relatively lacking compared to that of New-tonian fluids. As a result, this study focused mainly on the flow of liquidcrystalline materials. The following list presents the primary objectives ofthis study, in an order of noteworthiness.161.5. Objectives and ApproachTexture and Macro−Scale Flow AttributesThe nucleation and evolution of the liquid crystalline texture (fingerprints) / defects under flow−induced microstructure due to the pres-ence of the orientational order of LC molecules has to be modelled;any interactions between the liquid crystalline microstructure and themacro−scale flow attributes have to be taken into account to achieve abetter understanding of the nature of the rheological and tribologicalfeatures as lubricating agents.Molecular Orientation vs. RheologyThe alteration of the molecular orientation under simple shear flowand its impact on the apparent viscosity and resistant torque has to bestudied, so that the intrinsic rheological characteristics of LCs can beeffectively employed in the development of high performance / smartlubricating agents.Rheological FunctionsA semi−empirical methodology to predict the anisotropic viscositycoefficients and the rheological functions of lyotropic liquid crystallinematerials with cholesteric nature has to be proposed.When exploring the implementation of liquid crystalline materials as lu-bricating agents, a number of complexities need to be taken into account.They are associated with the flow of liquid crystalline materials, as well asthe interaction between microstructure / texture / defects and macro−scaleflow characteristics. The following is a listing of detailed objectives, each ofwhich corresponds to a dedicated chapter in the remainder of this text.• Initially, the Leslie−Ericksen theory was implemented to model thesteady flow−induced microstructure of several flow−aligning liquidcrystals between two concentric cylinders with a very small gap. Thechief objective in this part of the study was to capture and inves-tigate the multiplicity and multistability of results produced by theLeslie−Ericksen theory. An additional goal was set to explore how theanchoring angle of molecules on the wall boundaries influence the rhe-ological characteristics of LCs, and to compare the resistance torquegenerated by Newtonian fluid between two rubbing surfaces and theone of LCs.• Analyzing the connection between the Leslie−Ericksen and Landau−deGennes theories was the objective of the third chapter. By takinginto account the significance of several viscosity coefficients in the171.6. Thesis LayoutLeslie−Ericksen and Landau−de Gennes governing equations, it waspossible to define a methodology to calculate the values of Leslie /Landau viscosity coefficients by combination of theoretical / experi-mental methods. To test the accuracy of the proposed method, therheological data on the shear stress and the apparent viscosity of cer-tain lyotropic LCs was compared with numerical simulations of thestart−up simple shear flow using the Landau−de Gennes theory andthe proposed Landau viscosity coefficients.• In the final stage of the study, the leading objective was to understandthe formation and evolution of LC microstructures in the vicinity ofsolid surfaces as a protective layer against wearing processes, as well asits implications on the macro−scale flow attributes between two rub-bing surfaces of eccentric cylinders. In order to achieve this goal, tran-sient flow−induced microstructures and microstructure−induced flowwere numerically simulated by implementing the Landau−de Gennestheory fully−coupled with the Navier−Stokes equations. The flow oflubricants in journal bearings was modelled by the flow between twoeccentric cylinders.The final and most critical objective of this thesis was a numericalstudy of a 3−dimensional start−up flow of liquid crystalline materialsbetween two eccentric semi−spheres (as a simplified model of pros-thetic hip joint), which provided a comprehensive understanding ofthe intricacies involved in the application of LCs as bio−lubricants.1.6 Thesis LayoutIn chapter 2, governing equations of the Leslie−Ericksen theory was pre-sented for steady flow of liquid crystalline materials between two concentriccylinders; alignment angle, entropy generation, free energy per unit volume,apparent viscosity and resistant torque was presented and also comparedwith rheological characteristics obtained from rheometry tests.In chapter 3, rheological functions and the Leslie / Landau viscosity coef-ficients were theoretically predicted for various concentration of NCC aque-ous suspensions; chiral patterns of the molecular orientation under start−upsimple shear flow were explored. Shear stress and apparent shear viscositywere numerically calculated and compared with rheological data.In chapter 4, governing equations of the Landau−de Gennes theory werepresented for 2−dimensional (3−dimensional) transient flow of liquid crys-talline materials between two eccentric cylinders (two eccentric semi−spheres)181.6. Thesis Layoutunder start−up flow by rotation of inner cylinder. Flow−induced microstruc-ture and microstructure−induced flow were investigated; schematic repre-sentation of LC molecules, scalar order parameter as a measure for distri-bution of LC molecules, pressure, flow streamlines and first normal stressdifference were calculated and compared with the related one of the hypo-thetical similar Newtonian fluids.In chapter 5, highlights and summary of the previous sections were il-lustrated, followed by a short discussion on the potential future researchdirections; at the end references based on the order of appearance in thetext and the appendix were arrayed.19Chapter 2Flow Modelling andRheological Characterizationof Nematic LCs UsingLeslie−Ericksen TheoryThe present chapter aims at simulating the steady flow of thermotropic flow-aligning nematic liquid crystalline materials (NLCs) between two concentriccylinders with small gap size, which is the preliminary geometry for journalbearings (modelled usually by eccentric cylinders). Using concentric cylin-ders rather than eccentric cylinders significantly reduces the computationalcost of numerical simulations due its symmetric nature. Furthermore, usingrheological characterization and simulations, the viscous torque correspond-ing to a NLC is compared to the torque for a Newtonian fluid.This chapter is restricted to orientation flow processes where the flowtime scale is always slower than the order parameter time scale. In thischapter we characterize the flow aligning rod-like NLCs and show that therheology is solely driven by orientation processes. The orientation variableis a slow, long wave-length mode. Thus the selection of the theory dependson the purpose and use of the model. For sufficiently slow flows, whenthe flow time scale is slower than the internal time scale such that theorientation dominates the rheology and long length scale phenomena, theLeslie-Ericksen theory is sufficient. In this chapter we use the Leslie-Ericksenvector theory to extract a clear rheological characterization of NLCs.The main objective of the present chapter is to describe the rheologyand orientation of the molecules of representative thermotropic rod-like ne-matic liquid crystals that aligns under weak flow with the average moleculesunit normal close to the velocity directions. Generalizations to other flowconditions will require special considerations to the material parameters ap-pearing in the constitutive equations.The particular objectives of this chapter are:202.1. Governing Equations of the Leslie-Ericksen Theory• characterize the orientation of flow-aligning NLCs and how this affectsthe rheology• compare the viscous torque for a liquid crystal and an equivalent New-tonian fluidThe organization of this chapter is as follow: In the next section we presentthe governing equations that describe the microstructure for rod-like ne-matic liquid crystals under arbitrary shear flow. Numerical and experi-mental setup for various NLCs are explained in the following sections. Inthe results section, multiplicity and multistability of numerical results fordifferent boundary conditions and initial guesses are presented. Compari-son between numerical and experimental results for viscosity and resistanttorque for different liquid crystals are shown in the following subsection.2.1 Governing Equations of the Leslie-EricksenTheoryThe shear flow behaviour and rheology of nematic liquid crystals (NLCs)depends on the sign and magnitude of the reactive parameter λ = f (µi)(considered as one of the material properties), which is the ratio of the flowaligning effect of the deformation rate and the tumbling (rotational) effectof the vorticity. For rod-like NLCs it is known that λ > 0 [64]. When λ > 1,the material flow-aligns close to the velocity direction since the rotationaleffect of vorticity is overcome by deformation. When 0 < λ < 1, the directordoes not align close to the velocity direction because the rotational effectof vorticity dominates over the aligning effect of deformation. Materialswith λ > 1 display the flow-aligning mode. The flow-alignment angle isknown as the Leslie angle and for rod-like NLCs exists for λ > 1. Forshear-aligning rod-like molecules, flow tends to align the average molecularunit normals along the velocity direction. The theoretical alignment angleis 2θ0 = cos−1(1λ).The governing equations based on the Leslie-Ericksen theory for theisothermal flow of liquid crystals, including the conservation of mass, andlinear momentum and Oseen equations, are presented below [6, 8].∂ρ∂t+∇· (ρv) = 0 (2.1)ρDvDt= f +∇· t (2.2)212.1. Governing Equations of the Leslie-Ericksen Theory(a) (b) (c)!v nxyR1R2Figure 2.1: Schematic representation of LC molecules with respect to the di-rector n for (a) disk-like (b) rod-like, (c) schematic representation of directorvector of liquid crystalline molecules with respect to velocity vectorρ1D2nDt2= G + g +∇·pi (2.3)where ρ is density, v is velocity, f is the body force, ρ1 is a material constantthat has dimension of inertia per unit volume[ML−1], G is the externaldirector body force, g is the intrinsic director body force, pi is the intrinsicdirector surface stress tensor [6], n is a dimensionless vector, which repre-sents the direction of the preferred orientation for the LC molecules, also tis the total stress tensor. Figure 2.1 shows a schematic of the uniaxial rod-like nematic liquid crystalline phase. It consists of flat, rod-like moleculesmore or less aligned along a common direction, represented by the uniaxialdirector n.The major difference between isotropic liquids and LCs is the effect of ori-entational movement of the director n in the equations represented by theFrank elasticity function Fd [5, 6], which is the free energy per unit volume,that can be defined as,2Fd = k11 (∇·n)2 + k22 (n · ∇ × n)2 + k33 | n×∇× n |2 (2.4)where kii for i = 1, 2, 3 are Frank elasticity constants. Schematic represen-tations of major deformations related to each one of the Frank elasticityconstants [6, 19] are shown in figure 2.2. As the free energy has to be posi-tive definite, it can be proved that kii ≥ 0 for i = 1, 2, 3.Thus, by knowing the asymmetric total stress tensor t, the model is com-plete. The total stress tensor is t = t0 +t′ where the superscript ”0” denotesthe isothermal static deformation and the prime one denotes the hydrody-namic component.t0 = −p · I−(∂Fd∂∇n)· ∇nT (2.5)222.1. Governing Equations of the Leslie-Ericksen Theory(a) (b) c)(Figure 2.2: Principal deformation of liquid crystal molecules, (a) bend, (b)twist, (c) splayt′ = µ1 (nn : A) nn + µ2nN + µ3Nn + µ4A + µ5nn ·A + µ6A · nn (2.6)Here, p is an arbitrary (indeterminate) constant, µi for i = 1, 2, · · · , 6 arethe Leslie viscosities, A is the rate of deformation tensor, 2A = ∇v+(∇v)T,and N is the angular velocity of the director relative to the velocity of thefluid, Nj = n˙j −Wjknk and W is the vorticity tensor, 2W = ∇v− (∇v)T.The flow of rod-like nematic liquid crystals between two concentric cylinderscan be modelled using the equations 2.1−2.3. The director n was assumedto be a function of θ, which represents the angle between the molecules’director and streamlines and can be considered as a function of the radius,’r’.n =[sinθ, cosθ, 0](2.7)Therefore, equations 2.1−2.6 can be simplified and velocity vector can bewritten as,v =[0, rω (r) , 0](2.8)By introducing n and v in equations 2.1−2.6, the governing equations indimensionless form can be defined as,f˜(θ)[1r˜∂θ∂r˜+∂2θ∂r˜2]+12df˜(θ)dθ[1r˜2+ (∂θ∂r˜)2]+(γ˜1 + γ˜2cos2θ) r˜2∂ω˜∂r˜= 0 (2.9)1r˜∂∂r˜[r˜3g˜ (θ)∂ω˜∂r˜]= 0 (2.10)where the definitions of f˜ (θ) and g˜ (θ) are as follow,f˜ (θ) = cos2θ + ε · sin2θ (2.11)232.2. Numerical Setup and Validation2g˜ (θ) = 2µ˜1sin2θ · cos2θ + (µ˜5 − µ˜2) sin2θ + (µ˜6 − µ˜3) cos2θ + µ˜4 (2.12)The dimensionless parameters in this set of equations are presented in table2.1.ω˜ r˜ γ˜i µ˜iωL2η˜k11rLγiη˜µiη˜Table 2.1: List of dimensionless variablesHere, γ˜1 = µ˜2−µ˜3 and γ˜2 = µ˜5−µ˜6; the gap size between two cylinders isL; µ˜i for i = 1, 2, · · · , 6 are the dimensionless Leslie viscosity coefficients; η˜ isthe average of the three Miesowicz shear viscosities that describe anisotropicliquids and ε = k33/k11 [6, 64].2.2 Numerical Setup and ValidationEquations 2.9 and 2.10 are a set of fully non-linear, non-homogeneous andsecond order uni-dimensional ODEs that should be solved numerically. Twoof the most general methods of solving ODEs are shooting and relaxationmethods. Due to the nature of these equations and the high probabilityof reaching multiplicity in the results, the relaxation method was selectedto solve these equations. The continuous domain between two concentriccylinders was discretized using finite difference schemes in order to trans-form the set of ODEs to finite−difference equations or FDEs, which can beexpressed as a block diagonal mass matrix for discretized domain. In re-laxation method the final solution is determined by starting from a randompredication and improving and evolving it, iteratively [65]. In this study,variant of Newton−Raphson method was used for improving the predica-tion. By knowing the fact that, relaxation method allow us to use not onlythe dependent variables but also any algebraic correlations of the dependentvariables; thus, gradients of the main variables such as θ and ω can be de-termined as part of the solution, which speed up the calculation of apparentviscosity and torque on the inner cylinder.Four rod-like NLCs have been used in this study: DDA98, AZA99, PAA10and MBBA11. Alignment angles of all liquid crystalline materials that have8poly(4,4′−dioxy−2,2′−dimethylazoxybenzene−dodeccanediyl)922−dimethyl−3−hydroxy−8−methyl−azaspiracid10para−Azoxyanisole11N−(4−Methoxybenzylidene)− 4−butylaniline242.2. Numerical Setup and ValidationDDA9 AZA9 PAA MBBAθ ◦0 6.1980 7.1291 8.4175 5.7003Table 2.2: Alignment angles for the LCs used in this studyDDA9 AZA9 PAA MBBAµ1 −162 −1320 0.043 −0.018µ2 −170 −1595 −0.069 −0.1104µ3 −2 −25 −0.002 −0.0011µ4 16.01 208.9 0.068 0.0926µ5 162 1470 0.047 0.0779µ6 −10.1 −149.9 −0.023 −0.0336λ 1.238 1.032 1.059 1.02Table 2.3: List of the Leslie viscosities [Pa · s] and tumbling factors for theLCs used in this study [66]been employed in this study were provided in table 2.2. The Leslie viscositiesof all liquid crystalline materials implemented in this study were presentedin table 2.3 and Frank elasticity constants of selected LCs have been shownin table 2.4.In order to validate the simulation code we compared the numerical resultsfor the bulk angle at high rotational velocities with the theoretical align-ment angle of liquid crystalline molecules. For high rotational velocities ofinner cylinder the orientational angle of the bulk of molecules reached toits alignment value, 2θ0 = cos−1(1λ), that verified the fact that numericalsimulation is accurate.Two liquid crystals have been chosen to test our model, DDA9 and AZA9.For the numerical simulations, the gap size was chosen to be 10−4 m. Gridk11 × 107(dyn) k22 × 107(dyn) k33 × 107(dyn)MBBA 10.1 7.5 22PAA 8.5 7.2 17.8Table 2.4: Frank elasticity constants [6, 67]252.2. Numerical Setup and Validationindependency study was performed for 160, 200, 300, 400 and 550 nodesbetween the two cylinders on distribution of the orientational angle of themolecules.r0 0.2 0.4 0.6 0.8 10246!"= 100= 10 1= 102= 103= 104(deg) !!!!(a)r0 0.2 0.4 0.6 0.8 10246!"= 100= 10 1= 102= 103= 104(deg) !!!!(b)Figure 2.3: Distribution of θ vs. ω˜ for (a) DDA9 and (b) AZA9; initial guessfor θ was set to θ = pi/10 and for boundary conditions an anchoring angleof 0.0 (rad) was implementedSince the 300 node grid exhibited high enough accuracy in capturing the262.3. Experimental Methodgradient of orientational angle distribution close to the wall boundaries, itwas chosen to present the results. Infinity-norm of error, L∞ was used asconvergence criterion, where ‖Ei‖∞ = maxi |Ei| and threshold of acceptablerelative error was set to 10−8. Figure 2.3 illustrates the distribution of theorientation angle (θ) for anchoring angle of 0.0 [rad] on the wall boundariesand an initial distribution was set to θ = pi/10 [rad] for DDA9 and AZA9.The results are presented for different dimensionless rotational velocities, ω˜.By increasing the dimensionless rotational velocity, the orientation of themolecules reached the exact alignment angle that was predicted theoretically(see table 2.2).2.3 Experimental MethodA commercial MBBA (98%) supplied by Sigma−Aldrich (Oakville, Canada),a translucent yellow thermotropic LC with rod-like molecules was usedto validate numerical simulations. A rotational rheometer (Kinexus UltraRheometer Malvern Instruments Ltd.,Worcestershire, UK) was implementedto perform the rheological experiments in this study. Steady shear flow wasused to characterize the rheological behaviour of the liquid crystalline meso-phase. A fixture consisting of two stainless steel concentric cylinders with25 mm diameter and a gap between cylinders of 1 mm, at a controlled tem-perature of 23◦C was used. The rheometer was first calibrated with CannonCertified Viscosity Standard oil (Bellefonte, USA). In the steady shear test,shear rates from 10−1 to 103 [s−1] were applied to the samples. It is impor-tant to note that in some industrial applications (e.g. journal bearing) theshear rates in a normal regime are sometimes well over 103 [s−1].2.4 Results2.4.1 Multiplicity and MultistabilityPreviously, the accuracy of the model for capturing the alignment angle ofdifferent liquid crystals was assessed. Results were calculated for numer-ous types of boundary conditions, only few were presented here. Differentstudies confirmed the fact that LC molecules orient themselves completelyparallel to the solid surface or with a certain angle, which depends on theseveral different factors. Study of different mechanism for obtaining diverseanchoring angle of molecules near solid surfaces was not in the scope ofthis study. However, five distinct orientational angles were implemented on272.4. Resultsthe inner and outer cylinders to study the impact of anchoring angle of themolecules on the bulk of fluid that have been presented in table 2.5. FiguresB.Ca. 1 B.C. 2 B.C. 3 B.C. 4 B.C. 5θ ◦ (rad) 0 (0) 5.73 (0.1) -5.73 (-0.1) 11.46 (0.2) -11.46 (-0.2)a Boundary ConditionTable 2.5: The anchoring angles of molecules on the inner and outer cylindersfor the boundary conditions used in this study2.4, 2.5, 2.6, 2.7 and 2.8 present the orientational angles of MBBA and PAAmolecules vs. different dimensionless rotational velocities and θ = pi/10 asinitial guess for the distribution of molecules, and corresponding to bound-ary conditions B.C. 1, B.C. 2, B.C. 3, B.C. 4 and B.C. 5 respectively. Thepositive initial guesses for distribution of the molecules were chosen close tothe stable solutions.At very low rotational velocities elasticity is dominant and the orienta-tional angle of the molecules in the bulk of fluid depends on the anchoringangle at the wall boundaries rather than the viscous terms. On the otherhand, as the rotational velocity increases there is a competition betweenthe contributions of elastic and viscous terms to dominate the orientation ofmolecules. By increasing the dimensionless rotational velocity of the innercylinder a considerable part of the domain between the cylinders has beenaffected and a larger number of molecules reached the alignment angle.At high rotational velocities, the viscous terms become dominant and thedirector reaches the alignment angle, exception a thin boundary layer closeto the boundaries. The alignment angles obtained from numerical resultsare consistent with the theoretical values.Figures 2.9, 2.10 and 2.11 present the orientational angles for MBBA andPAA for different dimensionless rotational velocities and initial guess for thedistribution of molecules θ = −pi/10, for boundary conditions correspondingto B.C. 1, B.C. 2 and B.C. 3 respectively.For both liquid crystals, it was observed that changing the sign of theinitial guess for the orientational angle from positive to negative would causeinstabilities in results, and even divergence for high dimensionless rotationalvelocities. This phenomenon has been reported by other authors as well[23].The negative initial guesses for distribution of the molecules were chosenclose to the unstable solutions. As illustrated in figure 2.12 positive ini-282.4. Resultsr0 0.2 0.4 0.6 0.8 10246= 10-4= 10-2= 100= 102= 104!"(deg) """"(a)r0 0.2 0.4 0.6 0.8 102468= 10-4= 10-2= 100= 102= 104!"(deg) """"(b)Figure 2.4: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = pi/10 and forboundary conditions case B.C. 1 was implemented292.4. Resultsr0 0.2 0.4 0.6 0.8 15.75.7055.715.7155.725.7255.73! = 10-4= 10-2= 100= 102= 104"(deg)""""(a)r0 0.2 0.4 0.6 0.8 166.577.588.5= 10-4= 10-2= 100= 102= 104!"(deg) """"(b)Figure 2.5: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = pi/10 and forboundary conditions case B.C. 2 was implemented302.4. Resultsr0 0.2 0.4 0.6 0.8 1-6-4-20246!= 10-4= 10-2= 100= 102= 104"(deg) """"(a)r0 0.2 0.4 0.6 0.8 1-6-4-202468= 10-4= 10-2= 100= 102= 104!"(deg) """"(b)Figure 2.6: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = pi/10 and forboundary conditions case B.C. 3 was implemented312.4. Resultsr0 0.2 0.4 0.6 0.8 16810= 10-4= 10-2= 100= 102= 104!"(deg) """"(a)r0 0.2 0.4 0.6 0.8 18.599.51010.51111.5= 10-4= 10-2= 100= 102= 104!"(deg) """"(b)Figure 2.7: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = pi/10 and forboundary conditions case B.C. 4 was implemented322.4. Resultsr0 0.2 0.4 0.6 0.8 1-12-10-8-6-4-20246= 10-4= 10-2= 100= 102= 104!"(deg) """"(a)r0 0.2 0.4 0.6 0.8 1-10-505= 10-4= 10-2= 100= 102= 104!"(deg) """"(b)Figure 2.8: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = pi/10 and forboundary conditions case B.C. 5 was implemented332.4. Resultsr0 0.2 0.4 0.6 0.8 1-10-8-6-4-202!= 10-4= 10-3= 10-2= 10-1= 100(deg)"""""(a)r0 0.2 0.4 0.6 0.8 1-18-16-14-12-10-8-6-4-202468= 10-4= 10-3= 10-2= 10-1= 100!(deg)"""""(b)Figure 2.9: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = −pi/10 andfor boundary conditions case B.C. 1 was implemented342.4. Resultsr0 0.2 0.4 0.6 0.8 1-15-10-505= 10-4= 10-3= 10-2= 10-1= 100!!!!!"(deg)(a)r0 0.2 0.4 0.6 0.8 1-20-15-10-50510= 10-4= 10-3= 10-2= 10-1= 100!(deg)"""""(b)Figure 2.10: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = −pi/10 andfor boundary conditions case B.C. 2 was implemented 352.4. Resultsr0 0.2 0.4 0.6 0.8 1-10-8-6-4-20246!(deg)= 10-2= 10-1= 100= 101= 102"""""(a)r0 0.2 0.4 0.6 0.8 1-18-16-14-12-10-8-6-4-202468= 10-4= 10-3= 10-2= 10-1= 100!(deg)"""""(b)Figure 2.11: Distribution of orientational angle θ between two concentriccylinders for different dimensionless rotational velocities ω˜ of inner cylinderfor (a) MBBA and (b) PAA; initial distribution was set to θ = −pi/10 andfor boundary conditions case B.C. 3 was implemented362.4. Resultsv(SAA)(UAA)i>0i<0(a)Figure 2.12: Direction of change in orientational angle of liquid crystallinemolecules from unstable alignment angle (UAA) to stable alignment angle(SAA)tial distribution of molecules has better chance of reaching stable alignmentangle while negative initial distribution of molecules are more probable toreach unstable alignment angle.At very low rotational velocities the elastic term was dominant and thecontribution of the anchoring angle of the molecules at the boundaries wasmore significant. Increasing the dimensionless rotational velocity, the vis-cous terms became more significant and unstable alignment angles werereached. This was an anomaly that observed for all different liquid crys-talline materials and every boundary condition that had been tested usingthe model. In all of these figures orientational director of the moleculesof MBBA and PAA oscillate between their intrinsic stable alignment angleand the unstable angle. By further increase of the dimensionless rotationalvelocities ω˜ ≥ 10−1 these fluctuations propagated inside the system withhigher frequencies.Unstable solutions are also the answer to the Leslie-Ericksen equationsbut do not possess the minimum amount of Frank elasticity energy nor themaximum amount of entropy generation through the process. The rate ofentropy generation in the process of flowing a liquid crystalline materialbetween two concentric cylinders can be calculated by equation 2.13, [6, 27].Temp · s˙ = t′ji Aij − g′i Ni (2.13)Where Temp is the absolute temperature, t′ can be evaluated using equation2.6 and g is the intrinsic director body force. Equation 2.13 in dimensionless372.4. Resultsform can be simplified to,Temp · S˙ =12(r˜ ·dω˜dr˜)2·[2 g˜ (θ)](2.14)MBBA MBBA PAA PAAN. I.a P. I.b N. I. P. I.Fd(Jm3)0.00125 1.6417× 10−4 0.00120 4.2638× 10−4s˙(JkgKm3s)14.573 14.602 0.2476 0.2512a Positive initial distribution of moleculesb Negative initial distribution of moleculesTable 2.6: Frank elasticity and rate of entropy generation for negative /positive initial distribution of moleculesAmount of Frank elasticity energy and rate of entropy generation werepresented in table 2.6 for dimensionless rotational velocity ω˜ = 1 and dif-ferent initial guesses for the distribution of molecules of MBBA and PAAbased on the equations 2.4 and 2.14 respectively.It can be concluded that stable solutions always possess minimum amountof Frank elasticity energy while generating the higher amount of entropy.2.4.2 Viscosity and TorqueThe theoretical apparent viscosity [18] of selected liquid crystalline materialscan be calculated using equation 2.15;η =r2 − r1∫ r2r1dsg[θ (s)](2.15)Where g (θ) is the dimensional form of equation 2.12. Different types ofsurface treatments achieved various anchoring angles of the molecules onthe solid boundaries, which led to different flow curves. According to thenumerical experimentation that had been done on MBBA samples, it re-vealed that shear thinning and shear thickening behaviour of the materialwere functions of anchoring angles of the molecules on the solid surface.In order to verify the accuracy of the model in capturing the apparent viscos-ity of different liquid crystalline materials, set of experiments was performed382.4. Resultsin order to measure the apparent viscosity of MBBA. Two different rheolog-ical characterization tests with two concentric cylinders as geometry havebeen performed as described before. The non-isotropic behaviour of liquidcrystals lead to non-repeatability in results of viscosity.In figure 2.13 the apparent viscosity of MBBA for different anchoringangles of molecules on the wall boundaries along with the measured apparentviscosity were presented. A good agreement between the apparent viscosityobtained from numerical and experimental approaches with the error of lessthan %6 for relatively high shear rates was demonstrated. It was observedthat at lower shear rates the amount of torque produced due to the resistancein the material is within the range of the rheometer’s error (10−7 N.m). Itshould be noted that steady shear measurements are only reported for shearrates higher than 10−1 s−1, given that at lower shear rates the signals weretoo noisy due to torque limitations and inertia of the rheometer. The torqueper unit of length in the z−direction applied on the inner cylinder can becalculated from the contributions of shear stress and director surface stressas follow, [18].T = 2pi[r2trφ + r (zpq np pirq)](2.16)where t is the total shear stress tensor, pi is the director stress tensor and eijkare the components of the altering tensor. By using equation 2.6, appliedtorque on the inner cylinder can be calculated by,T = 2pir3(dωdr)g (θ) (2.17)Comparison between the torque on the inner cylinder estimated by numer-ical simulation and that obtained from the experiment was illustrated infigure 2.14. Good agreement between numerical and experimental resultswere observed for high shear rates. For low shear rates, fluctuations inthe experimental data due to physical limitations of the rotary rheometerimposed difficulty on obtaining accurate comparison between results.In order to compare the behaviour of liquid crystal material to Newto-nian fluid, simulations had been performed on MBBA and a hypotheticalNewtonian fluid with a viscosity having the same value of µ4. According tothe Leslie Ericksen theory, µ4 is the Leslie viscosity equivalent to the New-tonian viscosity of a fluid. Resistance torque exerted on the inner cylinderby the Newtonian fluid was compared with that exerted by MBBA, as pre-sented in figure 2.14. Up to %60 reduction on the resistance can be obtainedusing a liquid crystal in comparison with a hypothetical Newtonian fluid.392.4. Results10-4 10-3 10-2 10-1 100 101 102 103 104 1050.0280.0290.030.0310.0320.033![Pa.s]". [1/s]B.C.B.C.B.C.B.C.B.C.14523(a)10-4 10-3 10-2 10-1 100 101 102 103 104 1050.0250.030.0350.040.0450.05B.C.B.C.B.C.B.C.B.C.Cup & Bob - Test 1Cup & Bob - Test 2![Pa.s]". [1/s]14523(b)Figure 2.13: (a) Apparent viscosity of MBBA for different boundary con-ditions (b) Comparison between apparent viscosity of MBBA derived fromthe Leslie-Ericksen theory and the measured apparent viscosity402.4. ResultsT[N.m]!. [1/s]10-4 10-3 10-2 10-1 100 101 102 103 104 10510-1010-810-610-410-2100Num. B.C.1Cup & Bob - Test 1Cup & Bob - Test 2(a)10-4 10-3 10-2 10-1 100 101 102 103 1040.40.420.440.46!LCNewB.C. 1B.C. 4B.C. 5B.C. 2B.C. 3TT". [1/s](b)Figure 2.14: (a) The torque on the inner cylinder measured for MBBA andfrom numerical simulations for B.C. 1 as a function of shear rate (b) Ratiobetween the computed torque for MBBA and the computed torque for thehypothetical Newtonian fluid as a function of shear rate 412.5. Summary and Remarks2.5 Summary and RemarksThe Leslie-Ericksen theory was implemented to predict the evolution ofalignment angle of PAA, AZA9, DDA9 and MBBA between two concentriccylinders. Multiplicity and multistability of results for the Leslie-Ericksentheory was observed as a function of the initial state of LC molecules. Un-stable state for alignment angle of molecules were detected, which did notpossess the minimum free energy nor the maximum rate of entropy genera-tion in comparison with stable solutions.Investigation of the influence of rotational velocity of the inner cylinderon evolution of the alignment angle of molecules were performed which re-sulted in following conclusions. At high rotational velocities, viscous termsbecame dominant and the director reached the alignment angle for stablesolutions, for the case of unstable solution fluctuations in the alignment an-gle was observed that propagated inside the system with larger frequenciesat rotational velocities ω˜ ≥ 10−1.Finally, viscosity and resistance torque exerted on the inner cylinder wereestimated using the Leslie-Ericksen theory and compared with the data fromrheological experiments, that has been indicating a good agreement.42Chapter 3Estimation of ViscosityCoefficients and RheologicalFunctionsThe objective of the chapter is to propose a methodology to calculate differ-ent rheological functions and viscosity coefficients for lyotropic liquid crys-tals using analytical and experimental rheological data. The validation ofthe procedure has been done using numerical simulations of the Landau−deGennes equations for start up shear flow of 7 wt% NCC aqueous suspensionsand experimental rheological data. In addition, based on concentration andsize of molecules, different characteristic parametric values and dimension-less numbers from the Landau−de Gennes theory have been calculated forNCC aqueous suspension.The organization of this chapter is as follows: In the section 3.1 wepresent an overview of Landau−de Gennes model for cholesteric nematicsunder arbitrary flow and introduce various parameters and variables in themodel. In section 3.2 we present the estimation of Landau viscosity coef-ficients and rheological functions based on the Leslie−Ericksen theory. Insection 3.3 we present the experimental methods. The calculation of differentrheological parametric values for aqueous suspensions of NCC is presentedin section 3.4. In section 3.5 we present the numerical methods for solvingthe Landau−de Gennes microstructure equation. In section 3.6, we presentthe numerical results for 7 wt% NCC aqueous suspension and the validationof the calculated Landau viscosity coefficients using experimental data.3.1 Nematodynamic Governing EquationsIn the Landau−de Gennes theory the microstructure of cholesteric liquidcrystals is described conveniently in terms of a second order, symmetric andtraceless tensor order parameter (Q). The tensor order parameter Q can be433.1. Nematodynamic Governing Equationsexpressed as a function of its main eigenvalues and eigenvectors,Q = S(nn− 13I)+ 13P (mm− ll) (3.1)In equation 3.1, I is the unit tensor, S = 32µn and P = 3µn +32 , whereµn, µm and µl are the largest, second largest and smallest eigenvalues of Q,corresponding to the eigenvectors n (uniaxial director), m and l respectively.The uniaxial director, n is a dimensionless vector, which represents thedirection of the preferred orientation of the liquid crystals’ molecules.The evolution equation describes the dynamics of cholesteric liquid crystalsthrough the following equation 3.1 [25, 35].Q̂ = F (∇v,Q) + Hsr(Dr,Q)+ Hlr(Dr,∇Q)+P (∇Q) (3.2)where Q̂ is the Jaumann derivative of the tensor order parameter,Q̂ =∂Q∂t+ (v˜ · ∇) Q− W˜ ·Q + Q · W˜ (3.3)F is the flow contribution,F= 23βA˜ + β[A˜ ·Q + Q · A˜− 23(A˜ : Q)I]− 12β{(A˜ : Q)Q + A˜ ·Q ·Q +Q · A˜ ·Q + Q ·Q · A˜−[(Q ·Q) : A˜]I}(3.4)Short−range elasticity (Hsr) arises from the inter−molecular forces andlong−range elasticity(Hlr)transmits the surface anchoring effects from thewall boundaries to the bulk of the fluid.Hsr= 1De[ (U3 − 1)Q + UQ ·Q− U (Q : Q) ·(Q + 13I) ]Hlr= 1Er{∇2Q + L?2[∇(∇ ·Q) +[∇(∇ ·Q)]T−23tr(∇(∇ ·Q))]}(3.5)P is the chiral term and is related to the ratio of elasticity to the helicaldistortions of molecules tending to increase the twist pattern of the directors[35].P = ΘEr(mikQmj,k + mjkQmi,k)(3.6)In the above equations, A˜ is the dimensionless rate of deformation tensor,W˜ is the dimensionless vorticity tensor, β is the molecular shape factor,U is the nematic potential, is the Levi−Civita symbol, L? is the ratioof Landau elastic coefficients, L? = L2/L1 where, L1 = K22/2S2 and L2 =(K11 −K22) /S2; Kii are Frank elasticity constants that depend on principal443.1. Nematodynamic Governing Equationsdeformation of molecules and S is the scalar order parameter. Moreover,four dimensionless numbers are defined in this set of governing equationsand can be described as follows: Ericksen number (Er), represents the ratioof viscous flow effects to long−range order elasticity, Deborah number (De),indicates the ratio of viscous flow effects to short−range order elasticity,energy ratio (R), measures the ratio of short−range order to long−rangeorder elasticity and Θ quantify the ratio of anisotropic distortional elasticityto the isotropic one.Er =ckBT γ˙H26L1Dr; De =γ˙6Dr; R =ckBTH2L1; Θ =L2L2; (3.7)In equation 3.7, c is the concentration of molecules per unit volume, kB isthe Boltzmann constant, T is the absolute temperature, γ˙ is the shear rate,H is the flow characteristic length scale, L (L) is the isotropic (anisotropic)long−range length scale [34] and Dr is the pre−averaged rotational diffusiv-ity [68] that is a process by which the equilibrium statistical distribution ofthe overall orientation of molecules is preserved.The Landau−de Gennes model incorporates the influence of microstruc-ture on macro−scale attributes of the flow by modifying the total stresstensor. Evolution of tensor order parameter along with mass balance andmodified momentum balance (using total stress tensor) are the governingequations for the flow of liquid crystals.The total stress tensor has two main components, symmetric and asym-metric. The symmetric component contains viscous and elastic terms. Di-mensionless viscous, component of total stress tensor can be expressed as afunction of tensor order parameter Q and dimensionless rate of deformationtensor A˜ as follows,τ˜ v= ν˜1A˜ + ν˜2[Q · A˜ + A˜ ·Q− 23(Q : A˜)I]+ ν˜4{(A˜ : Q)Q + A˜ ·Q ·Q+Q · A˜ ·Q + Q ·Q · A˜ +[(Q : Q) A˜]I}(3.8)In equation 3.8, ν˜i for i = 1, 2 and 3 are dimensionless Landau viscosity co-efficients. The other components of the total stress tensor were implementedthe same as [25].An important challenge related to the Landau−de Gennes model is thedifficulty to find the characteristic parametric values including the viscositycoefficients for a specific liquid crystalline material, connect the experimentaldata to the model and provide meaningful values for various parametersand coefficients. In the next section we use the Leslie−Ericksen theory453.2. Estimation of Rheological Functions and Landau Viscosity Coefficientsand Doi and Larson models for reactive parameter for lyotropic nematicliquid crystals, to calculate analytically different rheological functions andthe Landau viscosity coefficients.3.2 Estimation of Rheological Functions andLandau Viscosity CoefficientsThe Leslie−Ericksen theory applies to uniaxial nematic flows neglecting theshort−range order elasticity, and hence being unable to describe the changesof the scalar order parameter due to the imposition of sufficiently strong flow.Consequently, in this theory, the microstructure is described by the directorvector n and the scalar order parameter S is assumed to remain constant,that is, unaffected by the flow, and always equal to its value at equilibrium:S = Seq, while the biaxial order parameter P is equal to zero. The totalextra−stress tensor is t = t0 + t′ where the superscript ”0” indicates theisothermal static deformation and the prime one contains hydrodynamicterms [6, 17].t′ = α1 (nn : A) nn + α2nN + α3Nn + α4A + α5nn ·A + α6A · nn (3.9)In correlations 3.9, n is the director of LC molecules, αi for i = 1, 2, · · · , 6are the Leslie viscosities, N is the angular velocity of the director relative tothe velocity of the fluid [6], Nj = n˙j −Wjk · nk, where W is the vorticitytensor, 2W = ∇v − (∇v)T and A is the rate of deformation tensor, 2A =∇v + (∇v)T.The anisotropy of liquid crystals leads to numerous viscosity coefficientsin order to describe the flow properties. These coefficients strongly dependon the concentration of the LC molecules in the solvent and are sensitiveto molecular structure. The qualitative flow behaviour of a nematic liquidcrystal depends on the sign and magnitude of the reactive parameter orequivalently on the Leslie viscosity coefficients α2 and α3. When shearinga nematic liquid crystal two different types of flow behaviour are possible,depending on the signs of α2 and α3. For rod−like molecules, α2 is alwaysnegative, while α3 can be negative for flow alignment systems (λ > 1) orpositive for non−alignment systems (1 > λ > 0) [6].Leslie viscosities can be expressed in terms of two main components,elastic and viscous, αi = αEi + αVi . Viscous components of Leslie viscosities463.2. Estimation of Rheological Functions and Landau Viscosity Coefficientsare given by [1, 31].αV1 =10Cη¯β?V S4 αV4 =421Cη¯β?V (7− 10S2 + 3S4)αV2 =αV3 = 0 αV5 = αV6 =207Cη¯β?V (S2 − S4) (3.10)Elastic components of Leslie viscosities have been formulated by [29, 30] asexpressed in equations 3.11.αE1 =-2η¯β2S4 αE4 = η¯β2(235)(7− 5S2 − 2S4)αE2 =-η¯β(1 +1λ)S2 αE5 = η¯β[17β (3S2 + 4S4) + S2]αE3 =-η¯β(1−1λ)S2 αE6 = η¯β[17β (3S2 + 4S4)− S2](3.11)According to [1], C is the dimensionless concentration, C = cc? where c is thenumber of molecules per unit volume and c? is a characteristics concentra-tion proposed by Doi et al. [30] to be a function of the molecular dimensions,c? = 16/pidl2 , where l and d are length and diameter of the LC moleculesrespectively. Moreover, β?V is a dimensionless group equal to the ratio ofviscous to elastic contributions, β?V = (ξstrD?r) /kBT , where ξstr has a na-ture of drag coefficient, ξstr = (kBT ) /2Dr0; η¯ is the characteristic viscositythat can be expressed as, η¯ = c?kBT/(10D¯?r), where kB is the Boltzmannconstant, T is the temperature and D¯?r represents the microstructural rota-tional diffusivity at concentration equal to c?. According to [69] for rod likemolecules, D¯?r is a function of the rotational diffusivity in dilute solutions,Dr0,D¯?r = BDr0(c?l3)−2where Dr0 =3kBT[ln (p)− γ]piηsl3(3.12)In equation 3.12, B is the numerical correction parameter, ηs is the apparentviscosity of the solvent, γ is a constant to correct the edge effect which is afunction of the molecular aspect ratio (p) and is equal to 0.8 through thisstudy [68].In equations 3.10 and 3.11, S2 and S4 can be obtained using second andforth order Legendre polynomials, P2 and P4 respectively [1].S2≡ 〈P2 (x)〉 P2 (x) = 12(3x2 − 1)S4≡ 〈P4 (x)〉 P4 (x) = 18(35x4 − 30x2 + 3)(3.13)473.2. Estimation of Rheological Functions and Landau Viscosity CoefficientsVarious correlations have been proposed to estimate the value of reactiveparameter as a function of dimensionless concentration. According to [30],reactive parameter can be expressed as, λ = 1−0.091875C−2−0.02765C−4−0.00129C−6and according to [1] the value of reactive parameter as a func-tion of molecular aspect ratio, 2nd and 4th order Legendre polynomials wasreported as, λ = β (5S2 + 16S4 + 14) / (35S2), where β =(p2 − 1)/(p2 + 1)and p is the effective aspect ratio of LC molecules; for rod−like moleculesp = l/d. Both of them converge to one for high concentration of moleculesin solvent. In this study, the first reactive parameter was referred as λDand second one was referred as λL for the sake of simplicity. For the caseof lyotropic LCs S2 and S4 can be simplified and expressed as a function ofdimensionless concentration [29, 30].S2 =1− 0.14726C−2 − 3.344× 10−1C−4 − 1.49× 10−1C−6 + . . .S4 =1− 0.49087C−2 − 2.713× 10−1C−4 − 1.56× 10−1C−6 + . . .(3.14)Doi [30] implemented equation 3.12 to 3.14, to simplify correlations 3.11 asfollows,αEiη¯= ai1C + ai2C−1 + ai3C−3 for i = 1, 2, . . . , 6 (3.15)where, ai1, ai2 and ai3 from equation 3.15 were presented in table 3.1.i ai1 ai2 ai31 −52 0.767 0.1822 −52 −0.215 −0.0113 0 0.123 0.0514 0 0.123 0.0385 52 −0.153 −0.0816 0 −0.245 −0.041Table 3.1: Coefficients ai1, ai2 and ai3 used to calculate Leslie viscosities[30]According to [17, 70], requiring the entropy production to be positive,483.2. Estimation of Rheological Functions and Landau Viscosity CoefficientsVV Vn n na b cFigure 3.1: The Miesowicz viscosities, corresponding to the orientation ofthe director along one of the three axis ηa, for the director along the flowdirection, ηb for the director along the velocity gradient direction, and ηc forthe vorticity axisthe following thermodynamic inequalities must be satisfied:α4≥02α1 + 3α4 + 2α5 + 2α6≥02α4 + α5 + α6≥0−4γ1 (2α4 + α5 + α6)≥(α2 + α3 − γ2)2 (3.16)In equation 3.16 γ1 is rotational viscosity, γ1 = α2−α3 and γ2 is irrotationaltorque coefficient, γ2 = α5 − α6 [6].Beside Leslie viscosities, Miesowicz viscosities are important for describ-ing the anisotropic behaviour and the rheology of liquid crystalline states.[71]. The three Miesowicz viscosities correspond to three primary orienta-tions of the director of LC molecules fixed with an external force field asshown in figure 3.1.• Director parallel to the flow direction, ηa• Director parallel to the velocity gradient direction, ηb• Director parallel to the vorticity axis, ηcηa = 12 (α3 + α4 + α6) ηb =12 (α4 + α5 − α2) ηc =12α4 (3.17)Theoretically, in high enough shear rates, viscous forces are dominant andalign all the molecules approximately parallel to the flow direction and thevalue of apparent viscosities is equal to the alignment viscosity, ηal [24, 72,73]. Dimensionless alignment viscosity was presented in equation 3.18,η˜al = η˜a + 12 (η˜b − η˜a)(1−1λ)+α˜14(1−1λ2)(3.18)493.3. Experimental MethodWhere, α˜i (i = 1, . . . 6) are the dimensionless Leslie viscosity equal to αi/η¯and η˜i are the dimensionless Miesowicz viscosity equal to ηi/η¯.Finally, the values of viscosity coefficients from the Landau−de Gennesmodel can be obtained by mapping the total stress tensor from the Landau−deGennes model by the one from the Leslie−Ericksen model [24]. Based onGrecov’s calculations,α˜1=2ν˜4S2 − β2S2(89 −89S +S212 )α˜2=-13βS(2 + S −S22 )− S2α˜3=-13βS(2 + S −S22 ) + S2α˜4=ν˜1 − 23 ν˜2S +13 ν˜4S2 + 49β2(1− S − S24 )α˜5=ν˜2 + 13βS[2 + S − S22 + (4− S − S2)]α˜6=ν˜2 − 13βS[2 + S − S22 − (4− S − S2)](3.19)By rearranging equations 3.19, values of ν˜1, ν˜2 and ν˜4 can be obtained asa function of Leslie viscosity coefficients (αis), scalar order parameter (S)and β.3.3 Experimental MethodFreeze−dried NCC, produced according to [9], has been dispersed in DI wa-ter to prepare 5, 7 and 10 wt% concentrations of NCC aqueous suspensions.To ensure proper dispersion, the samples were sonicated (Ultrasonic proces-sor VCX-130, Sonics and Materials Inc.) for 1000 [Jg−1 of NCC]. Ultrasonictreatment was carried out in an ice bath to avoid sample overheating.The rheological measurements were performed using a rotational rheome-ter (MCR-501 Anton Paar Physica) with a parallel plate glass geometry of43 [mm] in diameter. The steady state shear viscosity versus shear ratecurves were generated for all samples for shear rates in the range from 0.01to 100 [s−1]. For each sample, the time required to reach steady state at 0.01[s−1] was determined by a transient test; the sampling time used to generatethe flow curves was decreased with increasing shear rate. To avoid sampleevaporation, a small amount of silicone oil was placed on the periphery ofsamples together with an evaporation blocker. All rheological measurementswere performed at temperature of 25◦C.503.4. Calculation of Rheological Material Functions and Landau Viscosity Coefficients for Aqueous Suspensions of NCC3.4 Calculation of Rheological MaterialFunctions and Landau Viscosity Coefficientsfor Aqueous Suspensions of NCCIn this section we present the calculation of rheological material functionsfrom the Leslie−Ericksen theory and Landau viscosity coefficients for aque-ous suspensions of NCC for a range of concentrations, 4 to 22 wt%. Be-side the concentrations, the most important input data is the size of themolecules. Aspect ratio and size of the particles were obtained by [9] throughTEM imaging techniques.Uren˜a et al. [32] reported the concentration of the NCC aqueous sus-pension in liquid crystalline phase consistent with the results of Shafiei et.al [9]. As discussed in previous section Leslie viscosities can be expressedas a function of dimensionless concentration of molecules per unit volumeof their suspensions. In figure 3.2, 3.3 and 3.4 values of dimensionless Leslieviscosities α˜i, Miesowicz viscosities η˜i, alignment η˜al and Landau ν˜i viscositycoefficients as a function of dimensionless concentration have been reportedfor NCC aqueous suspension.Dimensionless Concentration, CDimensionlesseslieViscosity,! i0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3-8-6-4-202468!1!2!3!4!5!6LFigure 3.2: Dimensionless Leslie viscosities α˜i as a function of dimensionlessconcentration cc? for NCC aqueous suspension513.4. Calculation of Rheological Material Functions and Landau Viscosity Coefficients for Aqueous Suspensions of NCCDimensionless Concentration, CDimensionlessMiesowiczViscosity! iandAlignmentViscosity! al0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 310-310-210-1100101! a! b! c!al : L!al : D~ ~~~~~~Figure 3.3: Dimensionless Miesowicz viscosities η˜i, and alignment viscosity,η˜al vs. dimensionless concentration cc? for NCC aqueous suspensionDimensionless Concentration, CDimensionlessLandauViscosity,! i1 1.5 2 2.5 3-2024!1!2!4~~~~Figure 3.4: Dimensionless Landau viscosity coefficients ν˜i vs. dimensionlessconcentration cc? for NCC aqueous suspension523.4. Calculation of Rheological Material Functions and Landau Viscosity Coefficients for Aqueous Suspensions of NCCDimensionless Concentration, C0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 310-210-1100101102103104!Larson!DoiDr " = 1.0 104Dr " = 2.5 104Dr " = 5.0 104Dr " = 7.5 104Dr " = 1.0 105Reactive Parameter and Rotational DiffusivityFigure 3.5: Reactive parameter, rotational diffusivity vs. C of NCC aqueoussuspension for different values of the numerical correction parameter, BDimensionless Concentration, C0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 310-210-1100101! " = 1.0 104! " = 2.5 104! " = 5.0 104! " = 7.5 104! " = 1.0 105CharacteristicViscosity [Pa.s]Figure 3.6: Characteristic viscosity η¯ vs. C of NCC aqueous suspension fordifferent values of the numerical correction parameter, B533.4. Calculation of Rheological Material Functions and Landau Viscosity Coefficients for Aqueous Suspensions of NCCBeside the aspect ratio and concentration of NCC molecules that hasa huge impact on the magnitude of characteristic viscosities, the value ofnumerical correction parameter B is equally important. In order to estimatethe value of B, a series of trial and error numerical experiments have beenperformed; in figures 3.5 and 3.6 values of reactive parameter (λ), rotationaldiffusivity(D¯r)and characteristic viscosity coefficient (η¯) for different val-ues of numerical correction parameter B have been reported. Viscosity ofAlignmentViscosity,! alandApparentViscosity! ap[Pa.s]0.04 0.06 0.08 0.110-310-210-1100101102!ap Rheometry!al " = 1.0 103!al " = 1.0 104!al " = 2.5 104!al " = 5.0 104!al " = 7.5 104!al " = 1.0 105Concentration [ ]wtFigure 3.7: Alignment viscosity of 5, 7 and 10 wt% NCC aqueous suspension,from experimental [9] and analytical calculations for different values of thenumerical correction factor, Blyotropic liquid crystalline state as a function of shear rate has three regions[31], starting with a shear thinning followed by plateau and second shearthinning regions. For aqueous suspensions of NCC, the experimental valueof apparent viscosity in the plateau region represents the flow−alignmentmode; the qualitative alignment state of liquid crystalline phase of the NCCaqueous suspension captured by polarized light microscopy imaging tech-nique employed by Shafiei et al. [9]. The alignment viscosity of aqueoussuspensions of NCC based on Doi model for reactive parameter has beenmatched with the experimental apparent viscosity of the same suspensionsin high shear rate (γ˙ ≈ 10) for concentrations of 5, 7 and 10 wt% NCC sus-pensions and the results have been reported in figure 3.7.543.5. Numerical MethodsAlignment viscosity with B = 5 × 104 had shown good agreement with ex-perimental results and for this value of B dimensionless ratio of viscous toelastic contributions β?V was equal to 0.003.3.5 Numerical MethodsIn order to verify the procedure of obtaining Leslie viscosities, Landau vis-cosity coefficients and the other calculated parametric characteristic val-ues, series of numerical experimentation were performed for shear flow of 7wt% NCC aqueous suspension. A set of partial differential equations (3.1)consisting of five equations and five unknowns was solved simultaneouslyfor transient simple shear flow between two parallel plates using COMSOLMultiphysics 4.3b. A linear velocity profile between the two plates was im-plemented (distance between the plates is one tenth of the length of thedomain), where the upper plate was moving and the bottom plate was fixed.COMSOL is a commercial package that uses finite element methods todiscretize continuously the liquid crystalline material domain. General PDE(G) module was used to predict evolution of flow−induced microstructureinside the domain. Parallel sparse direct solver (PARDISO) with adaptiverelaxation factor and second order backward differentiation formula (BDF )for time stepping were used to solve these set of equations. Acceptablerelative error was set to be equal to 10−6, mesh independency of the finalsolutions were accepted based on the standard methods. The final meshconsists of 4000 quadrilateral elements with average element quality of 0.963and 21105 number of degree of freedom (DOFs).In each one of these tests, start−up simple shear flow was numericallysimulated under different shear rates; initial condition was set to be equal toabsolute randomness between the molecules and periodic condition with zeroflux was set for the inlet and outlet of the domain. Moreover, the anchoringangle of NCC molecules was set to be parallel to the wall boundaries whichwas reflected in order parameter tensor as follows,QB = Seq ·23 0 00 -13 00 0 -13 (3.20)In equation 3.20, Seq represented the equilibrium state of the scalar orderparameter, Seq = 14 +34√1− 83U .In table 3.2, the list of all the characteristic parametric values and viscositycoefficients that has been used for 7 wt% NCC aqueous suspension was given.553.6. Selected ResultsThe values were calculated based on the size of particles and concentration.All the results from next section pertain to insignificant De numbers whereflow only affects macroscopic orientation.Er De R Θ102 − 108 10−7 − 10−1 109 2.2l [nm] d [nm] L1 [pN ] L2 [pN ]100 5 43.29 6.687λD c? [#a/m3] Dr? Dr0.8792 1.02× 1023 1041.8 1334ν˜1 ν˜2 ν˜4 U0.5057 0.3138 −0.4383 4a Number of moleculesTable 3.2: List of all characteristic parametric values for 7 wt% NCC aqueoussuspension implemented in the Landau−de Gennes model and extra−stresstensor for the start up shear flow simulations3.6 Selected ResultsThe modified shear stress was calculated using the Landau viscosity coeffi-cients ν˜1 = 0.5057, ν˜2 = 0.3138 and ν˜4 = −0.4383 in viscous contributionof the stress. The shear stress and apparent viscosity were compared withexperimental results in figure 3.8. Acceptable agreement between numericaland experimental results was observed validating the analytical predictionsfor Leslie viscosity coefficients [31] and for Landau viscosity coefficients [24].In figure 3.8, the numerical predictions of shear stress were slightly under-estimating the experimental results. For relatively low viscosity materials,the resistant torque, which is a basis of calculation in rheometry, is verysmall and measurement error due to the sensitivity of torque sensor is aviable source of error. Furthermore, the calculation of the Landau viscositycoefficients was based on analytical correlations and moreover on the numer-ical parameter B, selected to match the experimental alignment viscosity; aslight change in B has considerable influence on the value of Landau viscouscoefficients followed by change in shear stress magnitude, that potentiallycould be another source of the disparity between numerical and experimentalshear stress. Moreover, the counter influence of the microstructure on themacro−scale attributes of the flow (change in velocity profile and pressure563.6. Selected Results! [ s-1 ]"[Pa]# app[Pa. s]10-2 10-1 100 101 10210-210-110010110-210-1100101"Num"Exp#app - Num#app - ExpShear Rate,Shear Stress,Apparent Viscosity,Figure 3.8: Shear stress tensor and apparent viscosity as a function of shearrate, from numerical simulations and from rheological experiments for 7 wt%NCC aqueous suspensiondistribution) was not considered. A better prediction of apparent viscosityand shear stress would be expected by considering this interaction. Finally,another source of error could be the absence in the total stress of any extraterm related to the chiral elasticity.In figure 3.9, scalar order parameter, S = 32µn was presented for threedifferent shear rates, γ˙ = 1, 5 and 10[s−1]. For low shear rates chiralbehaviour of NCC aqueous suspension was observed; although, by increasingthe shear rate, flow aligning model was dominant inside the system. Inorder to have a better understanding of the orientational order of 7 wt%NCC aqueous suspension, magnified snapshots of schematic representationof molecules were presentedin figure 3.10 for different shear rates at t = 100[s]. Once again the transition from chirality to flow aligning by increasingshear rate is evident. Due to the fact that tensor order parameter hasnegative eigenvalues, for the schematic representation of molecules, a newtensor M was introduced; where M = Q + 13I. Triaxial ellipsoids were usedto represent NCC molecules. Major axis of ellipsoids was selected to complythe direction of the eigenvector corresponding to the largest eigenvalue ofM; while, first and second minor axes of the ellipsoids were aligning to the573.6. Selected ResultsS: 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90a)S: 0.683 0.684 0.684 0.685 0.685 0.686 0.686 0.687b)S: 0.684 0.685 0.686 0.687 0.688 0.689 0.700c)Figure 3.9: Scalar order parameter S at t = 100 [s] for shear rate γ˙ a) 1 b)5 and c) 10[s−1]direction of eigenvectors corresponding to the second largest and smallesteigenvalues of M at any given point in space [74]. The aspect ratio of threeaxes of each ellipsoids were obtained relative to the first, second and thirdeigenvalues of the modified tensor order parameter M.According to the polarized light microscopy tests performed by Shafieiet al [9] for 5 wt% concentration of NCC aqueous suspension, cholesteric fin-gerprints of liquid crystalline state was observed at low shear rates followedby flow alignment at higher shear rate. The same behaviour was observed innumerical experiments for 7 wt% NCC suspension in the bulk of the fluid,although close to the wall boundaries, a small layer of molecules was underthe influence of the strong anchoring. The cholesteric helix uncoiled andbecame flow aligned when the viscous flow torque has overcome the chiralelastic energy.583.6. Selected Results4 5 610-30.01.00.50.01.00.50.01.00.5Figure 3.10: Schematic representation of NCC molecules at t = 100 [s] forshear rate γ˙ (top) 1 (middle) 5 and (bottom) 10[s−1]593.7. Summary and Remarks3.7 Summary and RemarksA methodology was proposed to calculate the rheological functions andLeslie / Landau viscosity coefficients in order to solve Landau−de Gennesgoverning equations for NCC aqueous suspension which is a cholesteric ly-otropic liquid crystals.By adjusting the numerical correction parameter B, theoretical align-ment viscosity of NCC aqueous suspension has been matched with the ex-perimental apparent viscosity of the 5, 7 and 10 wt% NCC aqueous sus-pensions. Evolution of the microstructure of the liquid crystalline state hadbeen studied by implementing the Landau−de Gennes theory for transientsimple shear flow between two parallel plates on 7 wt% NCC aqueous sus-pension under known velocity profile. Values of the Landau viscosity coeffi-cients (ν˜i) were calculated based on the Leslie viscosity coefficients and wereimplemented in the total stress tensor in the Landau−de Gennes model.Apparent viscosity and shear stress were calculated and compared withthe experimental results for the same concentration. Scalar order parameterand schematic representation of the molecules were presented for three dif-ferent shear rates that exhibited transition from chirality to flow aligning byincreasing the value of the shear rate. Therefore, due to the chiral nature ofNCC aqueous suspensions, a cholesteric pattern was observed at low shearrates, however by increasing the shear rates flow−alignment regime wereachieved.60Chapter 4Landau−de GennesNematodynamics TheoryThe main objective of this chapter was to numerically simulate the behaviourof nematic LCs in a geometry similar to a journal bearing, in relation to theapplication of LCs as lubricants. Therefore, a complete set of equationsfor the evolution of the microstructure coupled with 2−D Navier−Stokesequations has been solved. In this chapter, the evolution of the microstruc-ture modelled using the Landau−de Gennes theory, was linked with theNavier−Stokes equations using a modified total stress tensor established inthe momentum equations. The flow characteristics and the orientation of theliquid crystal molecules for different rotational velocities of the inner cylin-der were investigated. Flow and rheological properties, defects nucleationand their evolution inside the domain have been studied. The flow−inducedstructural changes and the orientation−induced flow were analyzed for dif-ferent Reynolds numbers and Energy ratios. At the end, 3−d analysis of theflow of liquid crystalline materials vs. synovial fluid diagnosed with rheuma-toid arthritis was chosen in a simplified geometry of a prosthetic hip joint.This chapter was organized as follows. First we presented the governingequations that describe the microstructure and flow for lyotropic nematics.Then the numerical procedure was presented. Following this, the numericalresults for start−up shear flow was investigated, the effect of the Reynoldsnumber and the Energy Ratio on the microstructure and rheological prop-erties are presented. Numerical methods, selected results for 3−dimensionalflow of LC materials in a simplified prosthetic hip joint were presented.Finally conclusions were given.4.1 Theory and Governing EquationsIn this paper we study the flow of lyotropic liquid crystals between eccentriccylinders. The Landau-de Gennes nematodynamics equations are well suitedfor predicting texture formation since defects are non-singular solutions of614.1. Theory and Governing Equationsthe governing equations. The evolution equations based on the Landau-deGennes theory and the Navier Stokes equations have been used to modelthe isothermal flow of lyotropic liquid crystals.The evolution equation describes the microstructure of liquid crystallinestate using a traceless and symmetric second order tensor order parameter, Q. Based on the definition which has been given by de Gennes and ProstQ=∫(uu− 13)ωd2u [61], where u is the unit vector that is either normalto the surface of disk-like molecules or parallel to the long chain of rod-likemolecules of LCs and I represent the identity tensor (see Figure 2.1) and ωis the orientation distribution function.Furthermore, the tensor order parameter Q can also be expressed as a func-tion of its main eigenvalues and eigenvectors,Q = S(nn− 13I)+ 13P (mm− ll) (4.1)In equation 4.1, S=32µn and P=3µn +32µm, where µn, µm and µl are theeigenvalues corresponding to the largest eigenvector n (uniaxial director),the second largest eigenvector m and smallest eigenvector l respectively.The uniaxial director, n is a dimensionless vector, which represents thedirection of the preferred orientation of the liquid crystals’ molecules. Figure2.1 shows the orientation vector of each molecule with respect to n and theorientation angle θ.The evolution equation describes the dynamics of liquid crystals throughthe following equation,Q̂ = F (∇v,Q) + Hsr(Dr,Q)+ Hlr(Dr,∇Q)(4.2)In equation 4.2, F is the flow contributions,F= 23βA + β[A ·Q + Q ·A− 23 (A : Q) I]− 12β{(A : Q) Q + A ·Q ·Q+Q ·A ·Q + Q ·Q ·A−[(Q ·Q) : A]I}(4.3)where β is a shape parameter; A=12[∇v+(∇v)T]is the rate of deformationtensor and v is the velocity profile. Short-range elasticity (Hsr) arises fromthe inter-molecular forces and can be expressed as,Hsr = 6Dr[ (U3 − 1)Q + UQ ·Q− U (Q : Q) ·(Q + 13I) ](4.4)Where U is the nematic potential and Dr represents the microstructuralrotational diffusivity [68], a process by which the equilibrium statistical dis-tribution of the overall orientation of molecules is preserved and it is defined624.1. Theory and Governing Equationsby,Dr (Q) = Dr[1− 32 (Q : Q)]−2(4.5)In equation 4.5, Dr is the pre-averaged rotational diffusivity.Long-range elasticity(Hlr), transmits the surface anchoring effects fromthe wall boundaries to the bulk of the fluid.Hlr =6DrL1ckBT{∇2Q + L?2[∇ (∇ ·Q) +[∇ (∇ ·Q)]T−23tr(∇ (∇ ·Q))]}(4.6)In equation 4.6, L?=L2/L1 where Li are Landau coefficients; c represents theconcentration of molecules in liquid crystalline state; kB is the Boltzmannconstant; T is the absolute temperature. Q̂ is the Jaumann derivatives ofthe order parameter tensor that can be defined as,Q̂ =∂Q∂t+ (v · ∇) Q−W ·Q + Q ·W (4.7)where W is the vorticity tensor, W=12[∇v− (∇v)T].The evolution of the microstructure has a considerable impact on the veloc-ity profile and pressure distribution that brings up the necessity to modifythe liquid crystals’ stress tensor for completing the model and to take intoaccount the complexity of the physical phenomenon. The total stress ten-sor for the nematic liquid crystalline materials has three main components;the symmetric components, containing viscous and elastic terms, τ v and τ erespectively, τ s=τ v + τ e; asymmetric component, τ a and Ericksen stresstensor, τEr.τ t = τa + τ s + τEr (4.8)τ v=ν1A + ν2[Q ·A + A ·Q− 23 (Q : A) I]+ ν4{(A : Q) Q + A ·Q ·Q+Q ·A ·Q + Q ·Q ·A +[(Q : Q) A]I}(4.9)τ e=ckBT{−2β3 H− β[H ·Q + Q ·H− 23 (H : Q) I]+ β2{(H : Q) Q+H ·Q ·Q + Q ·H ·Q + Q ·Q ·H +[(Q : Q) I]}}(4.10)τ a=ckBT{H ·Q−Q ·H}(4.11)τEr=ckBTL1{L? (∇ ·Q) · (∇Q)T −∇Q : (∇Q)T}(4.12)634.2. Computational MethodsIn the above equations, H represents the combination of short-range andlong-range elasticity, H=Hsr+Hlr; νi are Landau viscosity coefficients thatcan be derived from mapping the shear stress of the Landau-de Gennestheory with the one of Leslie-Ericksen theory [24].For solving numerically the evolution equation along with modified Navier-Stokes equations, a set of dimensionless variables was introduced to simplifythe equations; (a) Ericksen number, Er, the ratio of viscous flow effect tolong-range elasticity (b) Deborah number, De, the ratio of viscous flow effectto short-range elasticity and (c) energy ratio, R, the ratio of short-rangeelasticity to long-range elasticity.Er =V hckBT6L1DrDe =V6hDrR =ckBTh2L1(4.13)In equation 4.13, V is the characteristic velocity or sliding velocity of theinner cylinder and h=R2−R1, where R1 is the inner cylinder radius and R2is the outer cylinder radius. By introducing the Ericksen and the Deborahnumbers, the dimensionless form of the evolution equation is,Q̂ = F +(1De)Hsr +(1Er)Hlr (4.14)Moreover, a new Reynolds number was defined as, Ren=ρV 2/ (ckBT ), whereρ is a characteristic density. Thus, the dimensionless modified Navier-Stokesequation for the isothermal flow of nematic liquid crystalline material canbe defined as,Ren ·[∂v?∂t?+ (v? · ∇?) v?]= −∇?p? +∇? · τ ?t (4.15)Where,τ ?t =τ tckBT= τ ?a + τ?e +(ErR)· τ ?v +(1R)· τ ?Er (4.16)Where, τ ?... described dimensionless shear stress for each component of thetotal shear stress.4.2 Computational MethodsContinuity equation along with equations 4.14 and 4.15 are a set of sixcoupled non-linear PDEs that inherited the parabolic nature from evolutionequation and hyperbolic and elliptic nature from Navier-Stokes equations.The governing equations of liquid crystalline state in two dimensions were644.2. Computational Methodssolved using Laminar flow module (SPF ) and PDE general module (G) ofCOMSOL Multiphysics 4.3b. COMSOL uses Galerkin FEM method and aMUltifrontal Massive Parallel sparse direct Solver (MUMPS) with adaptiverelaxation factor and variable backward differentiation formula (BDF ) fortime stepping to solve these set of equations. Acceptable relative error wasset to be equal to 10−8, absolute convergences and mesh independency ofthe final solutions were established based on the standard methods; in thiscase, wall shear stress on the inner cylinder were calculated for fine, finer,extra fine and extremely fine automatic meshing options in the software,since wall shear stress value was exactly the same for both extra fine andextremely fine methods, extra fine option was implemented to generate thefinal mesh. Final mesh has 48042 elements with minimum element qualityequal to 0.7475 and average element quality equal to 0.9862.The initial orientation of the molecules was set to be random in the flowplane, which was set by implementing the scalar order parameter, S = 0.In this case, if U ≤ 83 the state of the fluid is considered isotropic. For83 ≤ U ≤ 3, dual phase of the nematic and isotropic states are predicted.Moreover, for higher values of nematic potential, uniaxial nematic phaseof liquid crystals is expected. In this article, U = 4 was chosen for thenumerical simulations.p = 0V = 0V R .! "#R#$R%-1-0.8-0.6-0.4-0.20.00.20.40.60.81-0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 10.6 0.8Figure 4.1: Numerical domain and boundary conditions654.2. Computational MethodsNumerical domain was presented in figure 4.1, with the outer cylinderfixed and the inner cylinder rotating with a velocity V . The boundary condi-tions for velocity and pressure are shown in figure 4.1. The eccentricity ratiois ε/ (R2 −R1) =0.5, where ε is the eccentricity. The boundary conditionfor the tensor order parameter corresponding to a tangential orientation ofmolecules with respect to the wall boundaries was set using Q=QB,QB = Seq ·y2x2+y2 −13−x·yx2+y2 0−x·yx2+y2x2x2+y2 −13 00 0 −13 (4.17)In equation 4.17, the equilibrium scalar order parameter is, Seq=14+34√1− 83U ;the simulation parameters that have been used in this study are: β=0.99,L?=0.45, ν?1=1, ν?2=− 1 and ν?4=6, [25].h*0 0.2 0.4 0.6 0.8 101020304050!(deg) t = 102t = 103t = 104t = 105t = 101t = 100Figure 4.2: Orientation angle of MBBA as a function of the dimensionlessgap between the inner and outer cylinder at the centreline, Ren=1.4×10−4,De=1.3 and Er=2.4× 104 for different dimensionless timesIn order to evaluate the accuracy of the numerical method, transient twodimensional flow of MBBA12 between two concentric cylinders was solved12N−(4−Methoxybenzylidene)− 4−butylaniline664.3. Results and Discussionusing the Landau-de Gennes theory. The anchoring angle of the moleculeson wall boundaries was set to zero or the molecules were aligned perpendic-ular to the surface normal vector; initial condition was set to be completerandomness [75].For MBBA the Leslie coefficients at 30◦C were tabulated by [66, 76] intable 2.3. In general, for any flow-aligning nematic liquid crystals the align-ment angle is a function of rotational viscosity γ1=µ3 − µ2 and irrotationaltorque coefficient γ2=µ6−µ5 [23] and can be expressed as θ0=12 cos−1 (−γ1/γ2).By calculating rotational viscosity and irrotational torque coefficients basedon the data presented in table 2.3, we found for MBBA an alignment angleθ0=5.7◦.0.0 0.5 0.6 0.8 1.00.20.1 0.3 0.4 0.7 0.9h*Figure 4.3: Molecular orientation of MBBA, Ren=1.4 × 10−4, De=1.3,Er=2.4× 104, t=105 at the centreline for t=105Figure 4.2 shows the orientation angle, θ, of MBBA molecules at thecentreline, from numerical simulations based on the Landau-de Gennes the-ory. The results are consistent with the theoretical value for the alignmentangle, with an error less than 1%. Based on these results, the orientation ofmolecules between two concentric cylinders at the centreline was illustratedin figure 4.3.4.3 Results and DiscussionThe start-up shear flow of nematic liquid crystals between two eccentriccylinders was modelled using the Landau-de Gennes nematodynamics modelfor 10−4 ≤ Ren ≤ 102, Er=106 and 108, De=10 and R=105 and 107.This section is divided into three sub-sections. In the first sub-section,different flow properties and the microstructure are illustrated for differenttime steps for two different Reynolds numbers. In the second sub-sectionthe same properties and microstructure are presented as a function of the674.3. Results and DiscussionReynolds number for two time steps. Finally the impact of the Energy ratio(R) on the microstructure of LCs is discussed.4.3.1 Start−up FlowFlow, rheological properties and microstructure of liquid crystalline materi-als were presented for De=10, Er=108 and Ren=10−2 and 102.Before presenting the results, it is important to mention some significantaspects related to the flow of liquid crystals, already known from previ-ous studies, that will contribute to our understanding of the flow betweeneccentric cylinders:1. The flow between eccentric cylinders is very complex as it is a combina-tion of shear, rotational and extensional flows, which occur in differentregions of the domain [40].2. The flow of liquid crystalline materials is governed by the balance be-tween the inertial forces and viscous forces, and the structure is gov-erned by the balance between long/short-range elasticity and viscouscontributions.3. For high Ericksen numbers (Er) the influence of long−range elasticitycontributions is relatively small compared to viscous contribution.4. For Deborah number De ≥ 1 we have a molecular process and the flowaffects the eigenvalues and eigenvectors of Q.5. Based on the set of experiments performed by Larson et al. [77, 78] forlarge Ericksen numbers and De ≥ 5, shear flow uniform texture−freesamples are expected.6. Extensional flow aligns the molecules along the streamlines [40].In figures 4.4 and 4.5 the streamlines for start-up flow of liquid crystallinematerial between two eccentric cylinders were presented for Ren=10−2 and102 respectively. For Ren=10−2, the magnitude of viscous forces was smallcompared with the elastic contribution. The recirculations region was stretchedfor dimensionless time, t=50. By further marching in time quasi-steady statecondition was reached; the flow streamlines became more and more stableand the impact of microstructure evolution on the streamlines was decayingin time. Moreover, the occurrence of instabilities was limited to the centralpart of the recirculations region that was illustrated at t=500 and 1000.6810.80.60.40.20-0.2-0.4-0.6-0.8-1a) b)c) d)-0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.2 -0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.210.80.60.40.20-0.2-0.4-0.6-0.8-1Figure 4.4: Flow streamlines for De=10, Er=108 and Ren=10−2 at dimen-sionless times: a) t=50, b) t=100, c) t=500 and d) t=100010.80.60.40.20-0.2-0.4-0.6-0.8-1a) b)c) d)-0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.2 -0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.210.80.60.40.20-0.2-0.4-0.6-0.8-1Figure 4.5: Flow streamlines for De=10, Er=108 and Ren=102 at dimen-sionless times: a) t=50, b) t=100, c) t=500 and d) t=100069p0.3000.1800.060-0.060-0.180-0.300a )p0.3000.1800.060-0.060-0.180-0.300b)p0.3000.1800.060-0.060-0.180-0.300c )p0.3000.1800.060-0.060-0.180-0.300d )Figure 4.6: Dimensionless pressure contours for De=10, Er=108 andRen=10−2 at dimensionless times: a) t=50, b) t=100, c) t=500 and d)t=1000p22.00013.0004.000-5.000-14.000-23.000a )p22.00013.2004.400-4.400-13.200-22.000b)p25.00015.0005.000-5.000-15.000-25.000c )p25.00015.0005.000-5.000-15.000-25.000d )Figure 4.7: Dimensionless pressure contours for De=10, Er=108 andRen=102 at dimensionless times: a) t=50, b) t=100, c) t=500 and d) t=1000704.3. Results and DiscussionFor Ren=102, the magnitude of viscous forces was higher and becamedominant. The amplitude and frequency of the instabilities for the recir-culations streamlines were smaller and the quasi-steady state was reachedearlier than for Ren=10−2. It should be noted that for both Reynolds num-bers, the streamlines contours at steady state were qualitatively similar tothe Newtonian case, showing that the impact of microstructure was not verysignificant.In order to further study the impact of microstructure on the macro-scaleattributes of the flow of liquid crystalline material pressure distribution forvarious times was illustrated in figures 4.6 and 4.7. For Ren=10−2, bymarching in time, quasi-steady state was reached after t=100. However, theabsolute maximum and minimum values of pressure did not change withtime. For Ren=102, the quasi-steady state was reached earlier in time.In order to compare further the macro-scale attributes of the flow of liq-uid crystalline materials with Newtonian fluids, a hypothetical Newtonianfluid was chosen as follows: ν?1 from the liquid crystal viscous stress expres-sion from Landau-de Gennes theory represents the isotropic segment of theviscosity coefficients, therefore viscosity of the hypothetical Newtonian fluidwas set to ν1=Deν?1ckBTh/V . For each one of the cases with Ren=10−2and Ren=102, flow of a hypothetical Newtonian fluid was simulated in thesame geometry; for the hypothetical Newtonian fluid, Re = ρV h/ν1.Certainly, the start-up flow regime is important for journal bearings, forexample with respect to the initial load carrying capacity or the overshootof the friction torque. The figures 4.8, 4.9, 4.10, 4.11, 4.12 and 4.13 presentdimensionless pressure, dimensionless first normal stress difference (N?1 =τ ?txx − τ?tyy), and dimensionless wall shear stress along the inner cylinder forvarious times, showing the approach to steady state values. All the abovefigures show the liquid crystalline material behaviour compared with thehypothetical Newtonian fluid behaviour.The maximum value of dimensionless pressure, dimensionless first nor-mal stress difference and dimensionless wall shear stress for LC materialshad an overshoot around t = 100 for the low Reynolds number and t = 5for the high Reynolds number. From solid design point of view, informationregarding the maximum amount of stress or pressure inside the system isvery important.714.3. Results and Discussion!"#N1*-0.5 0 0.5 1-0.15-0.1-0.0500.050.10.15t=100t=250t=500t=750t=1000!"#N1*-0.5 0 0.5 1-0.0500.05t=100t=250t=500t=750t=1000Figure 4.8: Dimensionless first normal stress difference N?1 vs. dimensionlessrotational angle along the inner cylinder for various dimensionless times,De=10, Er=108 and Ren=10−2 (left), for Newtonian fluid, Re=10−3 (right)!"#N1*-0.5 0 0.5 1-10-50510t=5t=10t=25t=50t=1000!"#N1*-0.5 0 0.5 1-0.500.51t=5t=10t=25t=50t=1000Figure 4.9: Dimensionless first normal stress difference N?1 vs. dimensionlessrotational angle along the inner cylinder for various dimensionless times,De=10, Er=108 and Ren=102 (left), for Newtonian fluid, Re=10 (right)724.3. Results and Discussion!"#P*-0.5 0 0.5 1-0.6-0.4-0.200.20.40.6 t=100t=250t=500t=750t=1000!"#P*-0.5 0 0.5 1-0.2-0.15-0.1-0.0500.050.10.150.2t=100t=250t=500t=750t=1000Figure 4.10: Dimensionless pressure p? vs. dimensionless rotational anglealong the inner cylinder for various dimensionless times, De=10, Er=108and Ren=10−2 (left), for Newtonian fluid, Re=10−3 (right)!"#P*-0.5 0 0.5 1-40-2002040t=5t=10t=25t=50t=1000!"#P*-0.5 0 0.5 1-2.5-2-1.5-1-0.500.511.5 t=5t=10t=25t=50t=1000Figure 4.11: Dimensionless pressure p? vs. dimensionless rotational anglealong the inner cylinder for various dimensionless times, De=10, Er=108and Ren=102 (left), for Newtonian fluid, Re=10 (right)734.3. Results and Discussion!"#$* W-0.5 0 0.5 10.040.060.080.10.120.14t=100t=250t=500t=750t=1000!"#$* W-0.5 0 0.5 10.020.0250.030.0350.0450.050.0550.06t=100t=250t=500t=750t=1000Figure 4.12: Dimensionless wall shear stress τ ?w vs. dimensionless rota-tional angle along the inner cylinder for various dimensionless times, De=10,Er=108 and Ren=10−2 (left), for Newtonian fluid, Re=10−3 (right)!"#$* W-0.5 0 0.5 12345678t=5t=10t=25t=50t=1000!"#$* W-0.5 0 0.5 10.20.30.40.50.60.7t=5t=10t=25t=50t=1000Figure 4.13: Dimensionless wall shear stress τ ?w vs. dimensionless rota-tional angle along the inner cylinder for various dimensionless times, De=10,Er=108 and Ren=102 (left), for Newtonian fluid, Re=10 (right)744.3. Results and DiscussionBy comparing the difference between maximum and minimum amountof dimensionless pressure, it can be concluded that the liquid crystallinematerial has a greater pressure difference hence higher load carrying capac-ity. Not only the magnitude of the maximum pressure for LC was higherfor both Reynolds numbers comparing with hypothetical Newtonian fluidunder the same conditions, but the ratio between the maximum pressureof LC material with Ren=102 to Ren=10−2 was also larger comparing withthe same ratio for hypothetical Newtonian fluid.The quasi-steady state condition for the region close to the inner cylin-der was achieved by marching in time. In the case of LC, for Ren=102 flowproperties converged after few simulation time steps (around t=55) and forRen=10−2 they converged close to time step t=550. On the other hand,for the hypothetical Newtonian fluid for the low Reynolds number, the con-vergence to the steady state condition occurred faster compared with theone of high Reynolds number. Indeed, considering the impact of elasticityat Ren=10−2 for LC can explain this difference. Qualitatively, all the flowproperties are similar for LC and Newtonian fluid, although the magnitudeis different. These results are consistent with the results for flow streamlines.As discussed before, by marching in time, quasi-steady state was reachedfor all the flow properties, and a balance between elastic, viscous and in-ertial forces was reached. Therefore, the stability of liquid crystalline mi-crostructure under the same elastic, viscous and inertial forces was foresee-able through the simulation time. For rod-shape molecules of the liquidcrystalline material, stability is accessible through the alignment of directorvector with respect to the local velocity vectors.In order to analyze the order of the liquid crystalline molecules, thecontours of the scalar order parameter (S) were presented for Ren=10−2and 102 in figures 4.14 and 4.15 at different dimensionless time frames.Scalar order parameter corresponds to the largest eigenvalue of the orderparameter tensor and confines between zero and one. Higher values of scalarorder parameter indicates higher degree of orientational order in the domainwith one indicating absolute ideal order between molecules and lower valuesof S indicates lower degree of orientational order with zero corresponds toabsolute randomness between molecules or isotropic state. Furthermore,lower value of scalar order parameter occurs in vicinity of disclination linesand defect points. Nucleation and evolution of defects depends on variousfactors related to intrinsic properties of liquid crystals and the shape ofdomain.754.3. Results and DiscussionFigure 4.14: Distribution of scalar order parameter S between two eccentriccylinders for De=10, Er=108 and Ren=10−2 for dimensionless time s) t=1,b) t=10, c) t=50, d) t=100, e) t=500, and f) t=1000 764.3. Results and DiscussionFigure 4.15: Distribution of scalar order parameter S between two eccentriccylinders for De=10, Er=108 and Ren=102 for dimensionless times a) t=1,b) t=10, c) t=50, d) t=100, e) t=500, and f) t=1000 774.3. Results and DiscussionViscous contributions in evolution of the microstructure tend to alignall the molecules in the direction of flow; on the other hand, the elasticcontributions transmit the surface anchoring effects from the wall to thebulk of the fluid. For the case of lower Reynolds number Ren=10−2, theorientation of the molecules evolved from absolute randomness at t=0 tonearly ordered at t=100 due to the viscous forces obtained from rotation ofthe inner cylinder. Formation of a disclination line close to the outer rim ofthe gap was observed along with the development of the orientational order,that was stretched through the centre of the recirculation zone between theupward flow and downward flow.However, for the higher value of Reynolds number Ren=102 the patternwas more complicated; at t=20 nearly all the molecules in the domain werealigned except the vicinity of the disclination line; however, by marching intime, the disclination line became more and more unstable and small defectzones were generated and evolved through the simulation time.Through the simulation time, the flow in the small gap has been primar-ily shear flow that reached the flow-alignment mode in the early stages ofsimulation; no defect point or stretch of disclination lines have been observedin the small gap region through the simulation time. Despite the fact thatmicrostructure was evolving in time as it had been illustrated in figures 4.14and 4.15, the macro-scale attributes of the flow reached a quasi-steady stateat time step t=1000, even though the quasi-steady state of the microstruc-ture only occurred at the low Reynolds number. At high Reynolds number,the directors were essentially unsteady as it was observed at t=1000 andRen=102; this phenomenon has been observed before by other researchers[79].In order to obtain a broad understanding and exploring the interconnec-tion of flow and microstructure of liquid crystalline materials, molecular rep-resentation of LCs was illustrated in figures 4.16 and 4.17 for the Reynoldsnumbers equal to 10−2 and 102. Due to the complexity of the flow close tothe outer rim of the domain, magnified snapshots have been presented solelyfor the sake of better visualization of the microstructure. The fingerprinttexture of LCs was demonstrated and defects in the texture of the fluid weretracked over the simulation time. Due to the fact that tensor order param-eter has negative eigenvalues, for the schematic representation of molecules,a new tensor M=Q + 13I was introduced. Triaxial ellipsoids were used torepresent LC molecules. Major axis of ellipsoids was selected to complythe direction of the eigenvector corresponding to the largest eigenvalue ofM; while, first and second minor axes of the ellipsoids were aligning to thedirection of eigenvectors corresponding to the second largest and smallest784.3. Results and Discussioneigenvalues of M at any given point in space [74]. The aspect ratio of threeaxes of each ellipsoids were obtained relative to the first, second and thirdeigenvalues of the modified tensor order parameter M.In this study the initial condition was set to be absolute randomness forthe LC molecules; at first few dimensionless time steps the spatial deriva-tive of the eigenvalues of the tensor order parameter was negligible whichresulted in degradation of the triaxial ellipsoids to spheres. For the sake ofvisualization, the time steps with complete randomness were omitted fromthis report. In figure 4.16 for Ren=10−2 at t=50, part of the domain closeto the inner cylinder was in liquid crystalline state (shown by the ellipsoidalshape); however, part of the domain close to the outer rim was in isotropicstate. At t=100, two different domains were observed. One was close to theinner cylinder under the influence of viscous contributions with alignmentof the LC molecules along the flow streamlines. The second domain waslocated close to the outer cylinder under the influence of the large-elasticitycontributions and the curvature of the wall boundary. In the darker region,these two contributions were both present. At time t=500 and particularlyat t=1000 a disclination line and later defect points were nucleated betweenthese two domains.In order to categorize defects, the number of brushes [6], (Nd) connectedto the defect point has to be counted, Sd=14 (Nd). Most of the defects thathave been observed in this study could be categorized with strength, Sd=± 12 .Schematic representation of each category was illustrated in figure 4.16. Forthe high value of Reynolds number, the texture pattern was more complex.LC molecules aligned with smaller alignment angles with respect to thevelocity vector and two distinctive domains under the influence of viscousand elastic contributions occurred earlier in simulation time; the disclinationline developed as a result of their encounter. At t=500 the disclination linewas dissembled into smaller domains containing defects; the interactionsof these domains resulted in several different smaller defect domains thatpropagated inside the gap between the cylinders.794.3. Results and Discussion-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.51.110.90.80.70.60.5 1.110.90.80.70.60.5a) b)c) d)½½-Figure 4.16: Magnified snapshots of the ellipsoidal representation of the LCmolecules for De=10, Er=108 and Ren=10−2 for dimensionless times a)t=50, b) t=100, c) t=500, d) t=1000 and schematic representation of defectpoints with strength Sd=± 12 804.3. Results and Discussion-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.51.110.90.80.70.60.5 1.110.90.80.70.60.5a) b)c) d)½½-Figure 4.17: Magnified snapshots of the ellipsoidal representation of theLC molecules for De=10, Er=108 and Ren=102 for dimensionless times a)t=50, b) t=100, c) t=500, d) t=1000 and schematic representation of defectpoints with strength Sd=± 12 814.3. Results and Discussion4.3.2 Effect of Ren on Microstructure−Flow PropertiesIn this section, two simulation time steps were selected to compare theeffect of the Reynolds number on the flow and microstructure of NLCs.Other dimensionless groups were set as follows, De=10 and Er=108. Thecomparison range for Reynolds number was set to 10−4 ≤ Re ≤ 102; twotime steps (t=100 and t=1000) were chosen to make sure the start-up andquasi-steady states that have been described previously were captured andcompared for different values of the Reynolds numbers.!/"N* 1-0.5 0 0.5 1-8-6-4-202468 Re = 10-4Re = 10-2Re = 10-1Re = 100Re = 101Re = 102Figure 4.18: Dimensionless first normal stress difference N?1 vs. dimen-sionless rotational angle along the inner cylinder for various dimensionlessReynolds numbers, De=10, Er=108 at t=1000The dimensionless first normal stress difference, dimensionless pressureand dimensionless wall shear stress on the inner cylinder vs. dimension-less rotational angle along the inner cylinder for various Reynolds numberswere illustrated in figures 4.18, 4.19 and 4.20 at t = 1000. By increasingthe magnitude of the Reynolds number, value of dimensionless first nor-mal stress difference and dimensionless pressure were increased; the anti-symmetry that was observed for the dimensionless first normal stress differ-ence at low Reynolds numbers was disturbed at higher Reynolds numbers.With respect to the lubrication problem in journal bearings, the region withhigher magnitude of wall shear stress is probable to be more affected by the824.3. Results and Discussion!/"p*-0.5 0 0.5 1-30-20-100102030 Re = 10-4Re = 10-2Re = 10-1Re = 100Re = 101Re = 102Figure 4.19: Dimensionless pressure p? vs. dimensionless rotational anglealong the inner cylinder for various dimensionless Reynolds numbers, Renfor De=10, Er=108 at t=1000wear process.The influence of the Reynolds number on the scalar order parameterwas presented in figures 4.21 and 4.22. At t=100 and Ren=10−4, the weakviscous forces resulted in less orientational order in the most part of thedomain. However, by increasing the Reynolds number, the disclination lineoccurred in the vicinity of the outer rim and the width of the domain thatwas occupied by defects decreased due to the viscous contributions. Att=100 and Ren=102, only a small portion of the domain was not ordered.The t=1000 case was more complicated, since the disclination line becamemore and more unstable by increasing the Reynolds number and it wasdivided into several different groups of defects. Propagation and evolutionof the defects inside the domain was discussed in the previous section.Based on the results for scalar order parameter, the evolution of mi-crostructure was highly dependent on the value of the Reynolds number. Itwas expected to find the same influence for the schematic representation ofthe molecules. In figures 4.23 and 4.24 the schematic representation of LCmolecules between the eccentric cylinders was illustrated using ellipsoids fort=100 and t=1000. The same method that was described in the previous834.3. Results and Discussion!/"#* w-0.5 0 0.5 110-210-1100101102 Re = 10-4Re = 10-2Re = 10-1Re = 100Re = 101Re = 102Figure 4.20: Dimensionless wall shear stress τ ?w vs. dimensionless rotationalangle along the inner cylinder for various dimensionless Reynolds numbers,Ren for De=10, Er=108 at t=1000section was implemented to visualize LC molecules. For the sake of easierdetection of all minor alterations in orientation of the molecules, magnifiedsnapshots of domain were presented.As expected for t=100 and Ren=10−4 most of the domain was not af-fected by the rotation of the inner cylinder. However, by increasing theReynolds number, viscous contribution became predominant; at t=100 andRen=10−2 two distinct regions co-existed for a period of simulation time.One was under the influence of the viscous contributions due to the rotationof the inner cylinder and the other region was located close to the outer rimmore likely to be under the influence of the anchoring angle of the moleculeson the wall boundary due to the long- range elasticity contributions. Thedarker zone in between, exhibited isotropic state, in other words, the bal-ance between elastic and viscous contributions prevented the dark zone toexhibit any alteration from isotropic state.A polydomain structure was observed by increasing the Reynolds num-ber; the darker portion of the domain became thiner and turned into adisclination line. As observed in 4.23, two distinct flow-aligned domains co-existed at the same time inside the gap between the cylinders. One domain844.3. Results and Discussionwhich was predominant, had smaller alignment angle that tends to zero byincreasing the Reynolds number. The second domain kept the initial ori-entational angle between the director and the flow streamlines through thesimulation time. For t=1000 and Ren=102 the groups of defects predictedby scalar order parameter, were detected inside the domain; the directorbecame more unstable and an onset of defect nucleation was observed in theschematic representation of the molecules.Furthermore, the ordered layer of liquid crystalline molecules in vicinity ofsolid surfaces formed a protective layer against wear for the two rubbingsurfaces that makes the liquid crystalline materials potentially a viable can-didate to be used as lubricants.854.3. Results and DiscussionFigure 4.21: Scalar order parameter contours for De=10, Er=108 at t=100for various Reynolds numbers a) Ren=10−4, b) Ren=10−2, c) Ren=10−1,d) Ren=1, e) Ren=10, f) Ren=102 864.3. Results and DiscussionFigure 4.22: Scalar order parameter contours for De=10, Er=108 at t=1000for various Reynolds numbers a) Ren=10−4, b) Ren=10−2, c) Ren=10−1,d) Ren=1, e) Ren=10, f) Ren=102 87-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.51.110.90.80.70.60.5 1.110.90.80.70.60.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5a) b)c) d)e) f)Figure 4.23: Magnified snapshots of the ellipsoidal representation of the LCmolecules, De=10, Er=108 at t=100 vs. Reynolds numbers a) Ren=10−4,b) Ren=10−2, c) Ren=10−1, d) Ren=1, e) Ren=10, f) Ren=10288-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.51.110.90.80.70.60.5 1.110.90.80.70.60.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5a) b)c) d)e) f)Figure 4.24: Magnified snapshots of the ellipsoidal representation of the LCmolecules, De=10, Er=108 at t=1000 vs. Reynolds numbers a) Ren=10−4,b) Ren=10−2, c) Ren=10−1, d) Ren=1, e) Ren=10, f) Ren=102894.3. Results and Discussion4.3.3 Effect of Energy Ratio on Microstructure−FlowPropertiesIn order to investigate the influence of the Energy ratioR on the micro/macro-scale attributes of the flow, the following dimensionless variables were used:Ren=10−2 and 1, De=10 and Er=106and 108 which resulted in R=105 and107. As discussed before, by decreasing the Ericksen number, long-rangeelasticity contribution was increased comparing to viscous contribution.In figures 4.25 and 4.26, the flow streamlines and pressure distributionsbetween the eccentric cylinders were presented at t = 150. The time step waschosen to emphasize on the differences between the two energy ratios. Forlower value of Energy ratio, long-range elasticity contributions became moreimportant in comparison with short-range elasticity; therefore, numerousvortices were occurred inside the system due to the contribution of elasticityover the viscous contribution in transient mode.10.80.60.40.20-0.2-0.4-0.6-0.8-1-0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.210.80.60.40.20-0.2-0.4-0.6-0.8-1-0.8 -0.6 -0.4 -0.2 0 0.4 0.6 0.8 10.2a) b)c) d)Figure 4.25: Flow streamlines at t=150, a) R=107, Ren=10−2, b) R=105,Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=190p0.6000.3000.000-0.200-0.500-0.800p2.2001.2000.200-0.800-1.800-2.800a )c )p0.5000.1000.020-0.020-0.100-0.500p0.7000.300-0.017-0.200-0.600-1.000b)d )Figure 4.26: Dimensionless pressure contours at t=150, a) R=107,Ren=10−2, b) R=105, Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1S0.9500.8000.6500.5000.3500.2000.050S0.9500.8000.6500.5000.3500.2000.050a )c )S0.9500.8000.6500.5000.3500.2000.050S0.9500.8000.6500.5000.3500.2000.050b)d )Figure 4.27: Scalar order parameter at t=150, a) R=107, Ren=10−2, b)R=105, Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1914.3. Results and DiscussionAlthough, these vortices gradually disappeared by marching in time, forboth Reynolds numbers the streamlines exhibited the same type of vorticesfor the lower value of energy ratio; however, for Ren=1 these vortices oc-curred close to the outer rim comparing with the central location of the vor-tices for Ren=10−2. Moreover, in transient mode of R=105 and Ren=10−2the viscous contributions could not stabilize the flow even in the small gapbetween the two cylinders and a group of vortices occurred in this region.That was not observed in any other combination of dimensionless variables.Influence of the flow for two energy ratios (R) on the scalar order pa-rameter was presented in figure 4.27; as observed before, the disclinationline and centre of the recirculation zone were co-located inside the domainthrough the simulation time. At t=150 for Ren=10−2 by decreasing thethe energy ratio, several different defect domains were nucleated inside thesystem. The same polydomain structure was observed in schematic repre-sentation of the molecules in figure 4.28. Four isotropic domains along withordered domains co-existed to generate a remarkable texture. Inside theordered domain, several sub-domains with different orientations of the di-rector were observed for Ren=10−2 and R=105. By increasing the Reynoldsnumber, the ordered sub-domains were unified under the dominant viscousforces and the orientation of the directors coalesced mainly to the direc-tion of the velocity vectors. In other words, for transient mode when theelasticity is dominant, the flow was induced by evolution of microstructure.924.3. Results and Discussion-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.5-0.20.20.1-0.10.00.30.40.5-0.3-0.4-0.51.110.90.80.70.60.5 1.110.90.80.70.60.5a) b)c) d)Figure 4.28: Schematic representation of molecules at t=150, a) R=107,Ren=10−2, b) R=105, Ren=10−2, c) R=107, Ren=1 and d) R=105, Ren=1934.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Joint4.4 3−D Analysis of NLC Flow in a SimplifiedProsthetic Hip JointIsothermal flow of lyotropic nematic LCs under start−up flow between twoeccentric semi−spheres as a model for prosthetic hip joint’s capsular spacehas been numerically simulated by solving mass, momentum and evolutionof order parameter equations (see equations 4.14 and 4.16). In order to studythe flow−induced texture of liquid crystalline microstructure, formation andevolution of defects, Landau−de Gennes model has been implemented.4.4.1 Numerical Methods for 3−Dimensional ModelThe governing model that has been described by equations 4.14 to 4.16 isa set of ten coupled non−linear partial differential equations (PDEs) thatinherited their intrinsic parabolicity from the evolution equation, and the hy-perbolic and elliptic nature from Navier−Stokes equations.Three−dimensionalgoverning equations of an isothermal liquid crystalline material in simplifiedcapsular domain of an artificial hip joint for transient shear flow inducedby rotation of the artificial head of femur was numerically simulated withlaminar flow module (SPF ) and PDE general module (G) embedded inCOMSOL Multiphysics. COMSOL is a commercial CFD13 package thatdiscretize continuous domain using Finite Element methods and implementsMUltifrontal Massive Parallel sparse direct Solver (MUMPS) with adap-tive relaxation factor and variable backward differentiation formula to findtime steps.The space and time evolution of six components of the tensor orderparameter (Q11, Q12, Q13, Q22, Q23, Q33) were incorporated as a set of sixparabolic equations in G module; at each time step, evolution equations havebeen solved, and the components of the total stress tensor were evaluatedand incorporated in SPF module, which resulted in alteration of the velocityprofile and pressure distribution. At the next time step, modified velocityand pressure were implemented in evolution equations to apply the counterimpact of macro-scale attributes of the flow on microstructure.Mesh independency and convergence of the final solution was establishedaccording to the standard procedures. Acceptable relative infinity−norm oferror (L∞) was used as convergence criterion and threshold of acceptablerelative error was set to 10−6. The grid that was implemented in this studyhas 1, 808, 947 tetrahedral elements with minimum element quality of 0.2555,13Computational Fluid Dynamics944.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Jointx-0.5 0 0.5 1y -1-0.500.51z-0.500.51Figure 4.29: Schematic representation of the capsular space in prosthetichip jointaverage element quality of 0.7681 and 4, 728, 250 degrees of freedom (DOF).The numerical domain was represented in figure 4.29, all the measure-ments were in dimensionless format using the radius of the outer concaveof acetabulum as characteristic length; furthermore, geometry of the nu-merical domain was characterized by two dimensionless parameters, µ =(R2 −R1) /R1 ≈ 0.055 and = ε/ (R2 −R1) =0.2, where, R1 is the radiusof the head of femur, R2 is the radius of acetabulum and ε is the eccentricityalong x−axis. All the dimensions were selected based on [80, 81].Initial condition of molecules was set to absolute randomness or anisotropic state and for boundary conditions at the solid walls, rotationalvelocity of the head of artificial femur was selected to be equal to slow for-ward rotation. The boundary condition for tensor Q was set to ensure thetangential orientation of the anchoring angle of the molecules with respectto the wall boundaries as follows,QB=SeqA(y cosψ)2 -A3 yz sinψ cosψ-xy cos2 ψ -y2 sinψ cosψyz sinψ cosψ-xy cos2 ψ (z sinψ-x cosψ)2 -A3 xy sinψ cosψ-yz sin2 ψ-y2 sinψ cosψ xy sinψ cosψ-yz sin2 ψ (y sinψ)2 -A3(4.18)In equation 4.18, dummy variableA=y2+(z sinψ − x cosψ)2; where, x, y and zare spacial coordinates of the boundaries, and ψ is the transverse angle ofthe acetabular inlet plane [82]. Some other parameters that have been used954.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Jointin this study were set to certain values according to [25] as follows, β=0.99,L?2=0.45, ν?1=1, ν?2=− 1 and ν?4=6.4.4.2 Selected Results for 3−Dimensional ModelThe three−dimensional model which is a set of 10 equations and 10 un-knowns has 4, 728, 250 number of degrees of freedom. Therefore, due to thecomputational cost of solving those equations, a parametrical study wouldhave been very long and difficult. Thus, a set of dimensionless variableswas chosen in order to match the isotropic segment of viscosity coefficient(ν1) with the viscosity of synovial fluid from patients diagnosed with RA14[83], which is a Newtonian fluid. The dimensionless variables for nematicliquid crystals (NLCs) were selected as follows Ren=1, Er=9.756× 105 andDe=10. In the case of Newtonian fluids Reynolds number, Re = ρV h/ν1,was implemented. It has been proven that for large Ericksen numbers andDe ≥ 5, shear flow uniform texture−free samples are expected [77, 78].In the following three figures, in order to illustrate the magnitude ofdifferent flow attributes inside the domain, ten slices perpendicular to thex−axis in y−z plane were presented. In figures 4.30 and 4.31, dimensionlesspressure and first normal stress difference (N?1 =τ?txx − τ?tyy) for LC and forsynovial fluid affected by RA at dimensionless time t=1000 were presentedand compared. The increase in the value of pressure across the gap suggestsan increase in the load carrying capacity. In figure 4.31, the raise in thevalue of the first normal stress difference from synovial fluid affected by RAto liquid crystalline material was observed due to the viscoelastic nature ofLC materials.Dimensionless hydrodynamic wall shear stress on the surface of the pros-thetic head of femur was calculated for liquid crystalline material and syn-ovial fluid affected by RA and presented in figure 4.32. According to thisfigure, the magnitude of wall shear stress was higher for the case of liq-uid crystalline material and the region with higher magnitude of wall shearstress is probable to be more affected by the wear process. However, forLC the distribution of shear stress was more uniform, indicating that theprediction of life cycle for the head of femur will be more accurate.14Rheumatoid Arthritis964.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Jointx-0.500.5y-1-0.500.5z-0.500.51P*0.120.090.070.040.020.00-0.02-0.04-0.07-0.09-0.12a)x-0.500.5y-1-0.500.5z-0.500.51P*0.00140.00110.00080.00060.00030.0000-0.0003-0.0006-0.0008-0.0011-0.0014b)Figure 4.30: Contour slices of the dimensionless pressure for a) lyotropicliquid crystal with Er=9.756 × 105, De=10 and Ren=1 b) synovial fluiddiagnosed with RA, Re=0.1974.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Jointx-0.500.5y-1-0.500.5z-0.500.51N*10.150.120.090.060.030.00-0.03-0.06-0.09-0.12-0.15a)x-0.500.5y-1-0.500.5z-0.500.51N*10.00080.00060.00050.00030.0002-0.0000-0.0002-0.0003-0.0005-0.0006-0.0008b)Figure 4.31: Contour slices of the dimensionless first normal stress differencefor a) lyotropic liquid crystal with Er=9.756 × 105, De=10 and Ren=1 b)synovial fluid diagnosed with RA, Re=0.1984.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Joint0.070.060.050.040.030.020.010.001zyxzyx123456 10-4Figure 4.32: Contours of the dimensionless wall shear stress on the head ofartificial femur for a) lyotropic liquid crystal with Er=9.756× 105, De=10and Ren=1 b) synovial fluid affected by RA, Re=0.1994.4. 3−D Analysis of NLC Flow in a Simplified Prosthetic Hip Jointx-0.500.5y-1-0.500.5z-0.500.51Figure 4.33: Contours slices of the scalar order parameter, Er=9.756× 105,De=10 and Ren=1 at dimensionless times a) t=10 b) t=10001005.3. Future Research Directions• Investigate the electro−
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical simulations of flow and microstructure in...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical simulations of flow and microstructure in nematic liquid crystalline materials Noroozi, Nader 2013
pdf
Page Metadata
Item Metadata
Title | Numerical simulations of flow and microstructure in nematic liquid crystalline materials |
Creator |
Noroozi, Nader |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | Liquid crystals are known for their anisotropic characteristics, which lead to a preferred orientation of their molecules in the vicinity of solid surfaces. The ability of liquid crystalline materials to form ordered boundary layers with good load-carrying capacity and outstanding lubricating properties has been widely demonstrated. In order to study the advantages of implementing liquid crystals as lubricants, the steady state/time transient isothermal flow of thermotropic/lyotropic, nematic/chiral nematic liquid crystals between two concentric/eccentric cylinders and in planar Couette geometries were studied numerically. To consider the influence of the microstructure formation/evolution on the macro-scale attributes of the flow, the Leslie-Ericksen and Landau-de Gennes theories were employed. Simplicity of the Leslie-Ericksen theory in capturing the orientational alignment angle of the molecules makes it a viable candidate for modelling the flow of flow-aligning nematic liquid crystals. On the other hand, the Landau-de Gennes nematodynamics equations are well suited for predicting texture formation since defects and disclinations are non-singular solutions of the governing equations. The Landau-de Gennes theory for the liquid crystalline microstructure along with continuity and momentum equations were solved simultaneously using General PDE and Laminar Flow modules of COMSOL Multiphysics. The investigation of flow characteristics and orientation of liquid crystalline molecules for different rotational velocities/shear rates and anchoring angles at the boundaries were presented. Furthermore, nucleation and evolution of singularities in texture of the liquid crystalline materials were tracked over the simulation time. Moreover, alterations in the macro-scale attributes of the flow such as velocity profile, pressure distribution and first normal stress difference along with the evolution of defects were studied inside the liquid crystalline domain. The implementation of Landau-de Gennes nematodynamic governing equations for LCs flow simulations offered an insight in application of these materials as lubricants. It was shown the LCs could provide protection against the wearing mechanism by forming a shielding layer in the vicinity of solid surfaces. Three-dimensional simulations of a simplified prosthetic hip joint suggested that liquid crystalline materials should be considered as potential bio-lubricants. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-09-26 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0165526 |
URI | http://hdl.handle.net/2429/45131 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_fall_noroozi_nader.pdf [ 181.26MB ]
- Metadata
- JSON: 24-1.0165526.json
- JSON-LD: 24-1.0165526-ld.json
- RDF/XML (Pretty): 24-1.0165526-rdf.xml
- RDF/JSON: 24-1.0165526-rdf.json
- Turtle: 24-1.0165526-turtle.txt
- N-Triples: 24-1.0165526-rdf-ntriples.txt
- Original Record: 24-1.0165526-source.json
- Full Text
- 24-1.0165526-fulltext.txt
- Citation
- 24-1.0165526.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:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0165526/manifest