.‘ .,1,.. hn—unt

1 3 8:
. 1 r. v...)(.r .
(11.: :

mu.»

 

 

“SIR

3;“ 33 { LIBRARY

University

This is to certify that the
thesis entitled

EQUILIBRIUM MICROSTRUCTURE FOR LIQUID
CRYSTALLINE POLYMERS

presented by
Hemant Kamalakar Kini

has been accepted towards fulfillment
of the requirements for the

MS. degree in Chemical Engineering &
Materials Science

 

Michigan State I

 

@«kgi Vim

 

Major Professor’ 8 Signature
Q»
A MRJ \11’3 O 3
W /

Date

 

MSU is an Affirmative Action/Equal Opportunity Institution

—t-—‘—-—--~—o_-—-—--oo—.—o-.—.-.-

 

PLACE IN RETURN Box to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJClRC/DatoDuopes-ms

EQUILIBRIUM MICROSTRUCTURE FOR LIQUID CRYSTALLINE
POLYMERS

By

Hemant Kamalakar Kini

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Chemical Engineering and Materials Science

2003

ABSTRACT
EQUILIBRIUM MICROSTRUCTURE FOR LIQUID CRYSTALLINE
POLYMERS
By
Hemant Kamalakar Kini
The relaxation of the microstructure of liquid crystalline polymers subject
to Brownian motion and a Maier-Saupe (MS-) excluded volume potential is
predicted by the following dynamic equation for the second order moment of the

orientation distribution:

d<pp>

1
d1:_ =-(< 22>"3’D+U(<BE>'<EE>‘ <pp>I< 2222 >)

A recently developed algebraic closure for < ppp p > , the fourth order moment of

the orientation distribution function, provides a closure for the abOve equation.
The new theory predicts the existence of multiple equilibrium states for a
MS-coefficient U between 4.656 and 5.000. At steady state, an isotropic
orientation dyadic satisfies the moment equation for all values of U. The theory
shows that the isotropic equilibrium microstructure is stable for U < 4.656 and
unstable for U > 5.000. For 4.656 < U < 5.000, three steady states are predicted: a
conditionally stable isotropic state, an unstable anisotropic state, and a
conditionally stable anisotropic state. The equilibrium semi-positive quadratic
form associated with the orientation dyadic has a prolate shape (i.e., two equal
eigenvalues smaller than the maximum eigenvalue). The temporal relaxation of a
non-equilibrium state to an equilibrium state occurs for time scales comparable to

l / DR, where DR is a phenomenological rotary diffusion coefficient.

To my parents

iii

ACKNOWLEDGEMENTS

My utmost gratitude to my advisor Dr Charles A. Petty for all his help,
guidance and inspiration during the course of this research. It is due to his
unrelenting support and prodding that I was able to complete the research and
write up this thesis in the given time frame. My sincere thanks to Dr André
Bénard and other members of the research group for their valuable suggestions
towards my research.

I am obliged to my family and friends for their love and emotional
support.

I would also like to thank the National Science Foundation (N SF/CTS-
9986878) and the Department of Chemical Engineering and Materials Science for

their financial support to this research.

iv

TABLE OF CONTENTS

LIST OF TABLES
LIST OF FIGURES
LIST OF NOTATION

CHAPTER 1. BACKGROUND

1.1 Introduction

1.2 Characterization of the Microstructure

1.3 Dynamic Equation for the Orientation Dyadic
1.4 Closure Model for the Orientation Tetradic

1 .5 Objectives

1 .6 Methodology

1.7 Scope

1.8 Figures

CHAPTER 2. RESULTS

2.1. Steady States and. Equilibrium States
2.2. Relaxation to Equilibrium
2.3. Figures

CHAPTER 3. CONCLUSIONS

3.1 Tables and Figures

CHAPTER 4. RECOMMENDATIONS

APPENDIX A. INVARIANTS OF THE ORIENTATION DYADIC

APPENDIX B. SCALAR ORDER PARAMETER FOR THE
PROLATE STATES

APPENDIX C. TETRADIC CLOSURE AT THE NEMATIC STATE

APPENDIX D: NUMERICAL DATA with C2 = 1/3

APPENDD( E: NUMERICAL DATA with C2 = 0.31

APPENDIX F: NUMERICAL DATA with C2 = 0.29

APPENDIX G: COMPUTER CODE

LIST OF REFERENCES

Page
vi
viii
ix

10
10
15
15
18
19

25
27
29

40

41

45

47
51

53
54
57
60
63

95

LIST OF TABLES

Table 3.1 Comparison of the Uip, Upi, api , values for the isotropic-prolate
transition predicted by the FSQ-closure (C2=l/3) with the BL],

and the de-coupling approximation

Table 3.2 The Uip, Upi, api , values for the isotropic-prolate transition
predicted by the FSQ-closure with different C2’s.

Table A-1. Invariants of the structure tensor

Table D] The eigenvalues of g , the structure parameter and the invariants

of 2 for U less than Upi and with C2=1/3

Table D2 The eigenvalues of a__ , the structure parameter 0t and the invariants

of 2 for U greater than Uip and with C2=1/3

Table D.3 The eigenvalues of g , the structure parameter and the invariants

of 1:) for Up, S Us Uip and with C2=1/3.

Table E] The eigenvalues of g , the structure parameter and the invariants

of g for U less than Upi and with C2=O.31

Table E2 The eigenvalues of g , the structure parameter or and the invariants

of 2 for U greater than Uip and with C2=O.31

Table E3 The eigenvalues of a , the structure parameter and the invariants

Of 2 for Upi S U S Uip and With C2=O.31

vi

Page

42

43

49

54

55

56

57

58

59

Table F.l The eigenvalues of g , the structure parameter and the invariants

of 1:) for U less than Up; and with C2=0.29 60

Table F .2 The eigenvalues of g , the structure parameter'a and the invariants

of 1:) for U greater than Uip and with C2=0.29 61

Table R3 The eigenvalues of g , the structure parameter and the invariants

of 2 for Uip S Us Upi and with C2=0.29 62

vii

LIST OF FIGURES

Figure 1.1. Progression of order for liquid crystal phases formed by rod

like molecules

Figure 1.2a,b,c Isotropic phase, Nematic phase, Smectic phase

Figure 1.3 Representation of an axisymmetric polymer molecule

Figure 1.4 Orientation vector distribution over the unit sphere

Figure 1.5 Equilibrium Phase Diagram indicating the prolate-isotropic

transition and the isotropic-prolate transition.

Figure 2.1 The isotropic—prolate transition predicted by the FSQ-closure with
C2=1/ 3

Figure 2.2 The isotropic-prolate transition predicted by the F SQ-closure for

C2=1/3 with the region between Uip (5.0) and Up, (4.626) enlarged.

Figure 2.3 The isotropic-prolate transition predicted by the F SQ-closure with
C2=0.3 1

Figure 2.4 The isotropic-prolate transition predicted by the F SQ-closure for

C2=0.31with the region between Uip (4.9) and Upi (4.629) enlarged.

Figure 2.5 The isotropic-prolate transition predicted by the FSQ-closure with
C2=0.29

Figure 2.6 The isotropic-prolate transition predicted by the FSQ-closure for
C2=0.29 with the region between Uip (4.83) and
Upi (4.606) enlarged.

viii

Page

20

21

22

23

24

30

31

32

w

3

34

35

Figure 2.7 or versus 1 (dimensionless time) for the relaxation of LCP

microstructure from an initial nematic state with C2=1/3.

Figure 2.8 or versus 1: (dimensionless time) for the relaxation of LCP

microstructure from an initial nematic state with C2= 1/3

Figure 2.9 or versus 1 ( dimensionless time) for the relaxation of LCP
microstructure from different initial states with C2= 1/3
and U= 4.7.

Figure 2.10 or versus 1: ( dimensionless time) for the relaxation of LCP

microstructure from different initial states with C2= 1/3 and U=6.

Figure 3.1 Comparison of the isotropic-prolate transition predicted by the
F SQ-closure (C2=l/3 and 0.29) with the MU and the

de-coupling approximation.

Figure A.1 Invariant diagram of the structure tensor g. The area and

the boundaries of the hypothetical triangle represent the set

of realizable orientation states

ix

36

37

38

39

44

50

Symbol

2 = < 22 >

an, 822, 333
< 2222 >

at, 32, a3

2

b1, b2, b3
llb

lllb

ID

< 2222 >

Q

DR

LIST OF NOTATION

Orientation dyadic

Diagonal components of the orientation dayadic

Orientation tetradic
Eigenvalues of g

Anisotropic orientation dyadic
Eigenvalues of 2

Second invariant of 1:)

Third invariant of 2

Probability distribution function

Unit normal vector

Orientation tetradic

Scalar order parameter

Nematic strength

Critical U beyond which the isotropic states are unstable
U corresponding to prolate-isotropic transition
or corresponding to isotopic-prolate transition

Boltzman constant

Absolute temperature

Rotary diffusion coefficient

Dimensionless time

Symbol

C1 Coefficient of linear terms in FSQ-closure

C2 Coefficient of quadratic terms in FSQ-closure

xi

CHAPTER 1
BACKGROUND
1.1 Introduction

Liquid crystals are ubiquitous in the world of digital electronics in particular and in
our life in general. Nematic liquid crystals can be found in display devices like laptops,
calculators, watches to name a few. Cholesterics or chiral nematics (a special case of the
nematic state) are used as temperature indicators in medical diagnoses or in paints to
produce angle dependent iridescent colors. Smectic liquid crystals due to their layered
structure, have potential application as lubricants. Liquid crystals by their structure and
properties can be considered as intermediate between classical solids and classical
liquids. Liquid crystals are complex fluids, which can flow like liquids but their
mechanical properties are anisotropic like those of crystals. Most of the historic work on
liquid crystals was with low molecular weight materials.

Liquid crystalline polymers (LCPs) are not only liquid crystalline but also
viscoelastic. A polymeric mesophase (positional and/or orientational long-range order in
one or two dimensions) was first discovered by Bawden and Pirie (see page 503, Larson),
when they observed that, a solution of tobacco mosaic virus formed two phases above a
critical concentration, one of which was birefringent. A synthetic polymer solution,
poly(y-benzyl-L-glutamate) exhibiting liquid crystalline phase, was reported in 1950 by
Elliot and Ambrose (see page 1, Weiss,1990). Generally in solids the molecules possess
both positional and orientational order. This means that the centers of mass of the
molecules lie at specific locations and the molecular axes point in certain directions.

Spherically symmetric molecules can have positional order but no orientational order.

Such positional and orientational order is absent in classical liquids. In solids, thermal
motion may lead to small changes in the positions and orientations of the molecules, but
their positions are in general fixed at specific lattice points and the motion is with respect
to the perfect geometrical lattice. In liquids, the molecules diffuse freely and randomly in
all directions. In LCPs, the molecules diffuse through the sample while still retaining
some positional and orientational order. Liquid crystals exhibit long-range orientational
order. A liquid crystalline phase (mesomorphic phase) can either be the result of
temperature change (therrnotropic), or the addition of a solvent (lyotropic). Figure 1.1 at
the end of this chapter, shows the sequential progression of order for various liquid
crystal phases exhibited by rod-like particles. A very common example of a lyotropic
liquid crystalline polymer is l,4-phenyleneterephthalamide, sold commercially as Kevlar.
This progression wherein the liquid crystal phase is intermediate between an isotropic
liquid and crystalline solid can be described by three types of order: orientational order,
positional order and bond orientational order (see page 3, Collings and Patel, 1997).

The nematic phase is characterized by the long axes of the molecules, tending to
align along a certain direction as they move about in the sample. This is the only long
range order present in the nematic phase and hence its properties are those of an
anisotropic fluid. Figure 1.2b is a schematic representation of the molecular order in a
nematic phase. The molecules are represented by oblate ellipses. The ellipses can be seen
to be exhibiting some average orientation; however there is no long range order in the
relative positions of the ellipses. These phases are optically anisotropic and fluid-like,

without any long range positional order.

The nematic liquid crystals are used in display devices, wherein their direction of
preferred orientation can be switched by the application of an electric field. And with the
optical properties such as birefringence strongly dependent of the orientation, the
nematics can be used as optical switches. A smectic phase exhibits both positional and
orientational order. Smectic phases have layered structure, and the positional order may
be along the layer or between the layers. There can be three types of positional order:
short range positional (SRO) order wherein the order is evident only over a finite distance
as in a simple fluid; long range positional order (LRO) which is seen in a three
dimensional crystal; quasi-long range order wherein the order is between that of SRO and
LRO. Depending on the average orientation angle made by the molecules with the layers
there can be different smectic phases. An isotropic state does not posses long range
orientational order, the structure can be considered to be reminiscent of the proverbial
‘bag of nails’. The liquid crystalline phases can be identified by various techniques like
miscibility and optical observation of textures, the most comprehensive understanding of
the molecular order can be attained by the X-ray diffraction studies on oriented samples.

As in LCPs, the behavior of rigid rod like particles dispersed in a viscous fluid is
of great importance in the processing of composite materials. Due to their low cost, fast
cycle times and good mechanical properties, the use of composites has increased
tremendously over the years. Many composite materials are formed, by combining a
reinforcing fiber with a polymeric or metallic matrix. These materials are processed in
different ways depending on the type of reinforcement used and the fiber orientation in
the final composite, which influences the final composite mechanical properties. As the

composite is being formed, flow induces the fibers to be oriented, this orientation

template is retained in the finished product. This fiber orientation template in the
composite material controls its mechanical properties. The direction in which most of the
fibers are oriented is stronger and stiffer than the direction with the least orientation (see
Advani and Tucker 1987; Parks et a1. 1999; Tucker, 1988).

There are both theoretical and experimental models for correlating the mechanical
properties with the direction of fiber orientation. Hence there is an ample incentive for
understanding the complex link between the composite manufacturing process, the final
fiber orientations and the composite properties. This knowledge would allow one to
design and control the process to obtain the most desirable mechanical properties in the
composite material. Predicting full-scale processing flows, where the flow is
inhomogeneous, necessitates tracking of all the particles, which is a tedious task. In
polymeric liquids, any molecule experiences an average effect of its surroundings and the
problem reduces to calculating the mean properties from the behavior of individual
molecules. Though these assumptions simplify the problem, they are made while
retaining enough information to extract accurate predictions of the rheology and
molecular orientation (see F olgar and Tucker, 1984; Advani and Tucker 1987).

The change in concentration of rigid rod like polymers causes dramatic changes in
the rheological properties of LCPs and other polymeric solutions. Onsager and Flory (see
page 229, Doi, 1981) identified critical concentration for solutions of such molecules,
above which they form liquid crystals with a nematic structure. This theoretical
prediction can be validated by viscosity measurements and polarized light microscopy of
the microstructure of polyisocyanates as done by Aharoni (1979,1980). Both these

studies have shown the existence of a critical concentration of the polymer above which

there is a transition from isotropic state to an anisotropic state, as the polymer
concentration is increased. The phenomenological theory of Erickson, Leslie and Parodi
(see page 229, Doi, 1981), was successful in describing many dynamical phenomena in
nematic, thermotropic liquid crystals. But this theory does not predict the nonlinear
viscoelasticity, which is quite conspicuous in such systems. Doi (1981), proposed a
model based on molecular approach, which predicted a nonlinear viscoelasticity for
LCPs. Also the predicted dependence of the viscosity on molecular weight and
concentration agreed fairly well with experiments.

Using the fluids microstructure numerous fluid mechanical theories have been
formulated. Jeffery (1922), studied the fluid motion in the vicinity of a suspended
ellipsoid to explain the increase in viscosity due to its presence, in an otherwise
Newtonian fluid. For a suspension of non-interacting dumbbell shaped particles, Prager
(see page 33, Hand, 1961) deduced a constitutive equation for stress and equations
determining the preferred direction adopted by the particles. The orientation state of the
fibers can be represented in numerous ways. The probability distribution function for the
fibers is one way (see Advani and Tucker, 1990,1994; Cintra and Tucker, 1995). It is
suitable for problems involving nontrivial flow fields with planar orientation. A more
compact way of representing the fiber orientation is by using various parameters. The
formulation of these parameters depends on assumptions about the direction of the
principal axes of orientation and some symmetry in the probability distribution function.
An ideal description of fiber orientation would be general like the distribution function
and compact like orientation parameters. These requirements are satisfied by a set of

tensors that have been called tensorial order parameters or orientation-moment tensors.

The underlying theory for this research utilizes this form of representation for the
orientation state, as suggested by Hand (1962) and the Doi (seee Doi and Edwards, 1986)
theory for liquid crystalline polymers. The tensorial form of description while providing
an efficient way to compute flow-induced fiber orientation, requires an accurate closure
approximation for the higher-order moments of the orientation distribution function. Over
the years several closure approximations have been proposed, some of which are

discussed later.

1.2 Characterization of the Microstructure

1.2.1 Orientation Characterization and Distribution

The polymer molecule is considered to be a rigid rod uniform in length and diameter.

The two ends of this polymer are identical. The orientation of this molecule can be
represented by an unit vector p , directed along its molecular axis or by angles 0 and (I)
made by the molecular axis. Fig 1.3 shows these forms of representation. The two forms
of representation are related by,
pl=sin0 cos¢ p2=c080 p3=sin03in¢ . (1.1)
pl , p2 , p3 are the components of the unit vector 2 along the co-ordinate axes X,Y and
Z respectively.

As the ‘head’ and ‘tail’ of the molecule are indistinguishable from each other, the

representation of the polymer direction is unchanged by changing p to —p or 0 and (I)

as below

9 —> n-e ¢ —> n+¢ . (1.2)

A probability distribution function u; (0,4)) describes the orientation state of a
particle in space, it is the fraction of orientation states per unit area on a sphere of unit

radius. The probability P, of finding a particle between 01 and 01 + A0; and (I). and

¢1+ Ad) is given by

P(0$039+A9,¢$¢_<_¢+A¢)=w(0,¢)sin0AOA¢. (1.3)

This fimction is normalized by integrating the function over the entire sphere i.e., the

probability of finding a particle at some orientation over the sphere is unity,

fiw(0,¢,t)sin0d0d¢=l . (1.4)

From the definition of p we can write that

ME) = w(-g) ’ (1-5)

W(9a¢) =‘ll(7t-9a¢+7r) - (1-6)

1.2.2 Continuity Equation
Considering an arbitrary area over the unit sphere, as shown in Figure 1.4. This

area is bounded by the contour line C. Relative to a material frame of reference, the flux

of the orientation states across the boundary C is p w. Writing a balance of orientation

states across the surface A,

—l(gw)-gds = illw(e.<1>.t)dA
c th

(1.7)

Using Stokes theorem ( see Deen, 1998), the lefi-hand-side of Eq. (1.7) can be written

as,

-(p_\y)dA . (1.8)

The equation for the evolution of w for the spatially homogeneous LCP mixture can be

written as,

a .
—=-—-——- . 1.9
E (2 v) ( )

The term (2 \V) represents the rotary flux, where B is the angular velocity of the LCP

5
component and .6—p— is the surface gradient operator. The rotary flux can now be written

as the contributions from convective and diffusive fluxes,

9v 5 133v+(9v-93w) . (1.10)

f) 3 is the angular velocity caused by fiber/fluid interactions .

Now Eq. (1.9) can be written in terms of the convective and diffusive fluxes

6w_ 6 . 6 . .
333— (3p (gal!) (3p (2 BOW). (1.11)

1.2.3 Rotary Diffusive Flux
The rotary diffirsive flux has a contribution due to Brownian motion and a contribution

due to an excluded volume potential :

. . 5W 5 AUMS
_ =—D _+ ____
(g 23w {62 val—)(kBT) ,

where DR is the rotary diffusion coefficient and has the units of l/time. (1.12)
1.2.3.1 Maier-Saupe potential

Orientation distribution can be calculated from a nematic potential which signifies
the influence of one molecule or a rod’s orientation on that of its neighbor. Every
polymer molecule has a region around it wherein the stearic effect precludes other
molecules from that region. This is called the excluded volume. Onsager initially
pr0posed a theory to calculate such a potential. This theory is applicable at low
concentrations at which pairwise excluded volume interactions are the dominant ones.
The Doi theory utilizes the Maier-Saupe potential AUMS, which is more suitable for

thermotropic nematics.
3 1 3 l
AUMS= —§UkBT(pp—§I).(pp-§l). (1.13)

U is a dimensionless measure of the polymer number density n and is directly
proportional to the excluded volume between rigid rods. When the concentration of the

polymer is increased, the interaction forces become stronger, this increases the alignment

of polymer molecules. Hence there can be a transition from a random to an ordered phase

even in the absence of flow.

1.3 Dynamic Equation for the Orientation Dyadic

1.3.1 The Smoluchowski Equation

The Smoluchowski equation with p 3: 0 (no external flow field), can be written

as

QE=DRE.[6_‘V+ Pug—"fin , (1.14)

at 5p 5p Wag kg

The Smoluchowski equation gives the evolution of the orientation distribution function.
Multiplying Eq. (1.14) by p p and integrating the equation over the entire unit sphere

gives

d<pp>
d1?

9

1
= -(< 22> ‘§l)+U(< 22> ° < 22> - < 22*2222 >>
(1.15)

where ‘C E 6 DRt is the dimensionless time.

1.4 Closure models for orientation tetradic
1.4.1 Background

As evident from the earlier discussion the distribution function \u(9, <1), t;U)
provides a description of the fiber orientation state. A compact representation of the
orientation for numerical simulations is in the form of orientation tensors. The dyadic

product of the vector p with itself and the integration of this product, with the

10

distribution function over all possible directions, yields one set of orientation tensor. As
the distribution function is even, the odd order integrals are zero. The integrals of
importance are the even-ordered ones (see Irnhoff, 2000)., The second order tensor or

dyad can be written as
<p__p_>= jjppw(0,¢,t)sin0d0d¢ . (1.16)
Similarly the fourth order tensor or tetrad is

<pppp>=fipp pw(0,¢,t)sin0d0d¢ . (1.17)

The fourth order tensor must reduce to the second order tensor by contracting any two

orientation vectors:

<P.P.>=<P PP.P.>=<PP.PP.>=<PP.P.P >
1) s31) SIS] SIJS

=<p.ppp.>=<p.p up >=<p.p.pp >
rs SJ rs _] s 1 JS 8
(1.18)

Clearly, the orientation tetrad also has six fold symmetry inasmuch as

>=<P3P.P P >= <Pil3

>=< >
J k 1 PP PiPPP

le Ijk

<
pipjpkpl

= <pk pl pi p3- >=< pj p1 pi pk >=< p3- pk pi p1 >

(1.19)
An infinite number of even order tensors can be written. The higher the number of
orientation tensor, the better is the description of the orientation state. But the calculation

of higher order tensors is very cumbersome. Hence the second and fourth order tensors

which provide an adequate representation of the orientation state are commonly used.

11

When the equation for < pp > is written, it has the higher order moment < pppp >,
which needs to be determined. Whenever an equation for any finite moment of p is
written, it will always contain a higher order moment of B. So to solve for < pp > , we
need to close the higher order moment < pppp >. Some of the commonly used closures

and the closure used in this research (i.e., the FSQ- closure) are now discussed.

1.4.2 The linear closure model of Hand

This linear closure proposed by Hand (1962) is given by

l 1
<pppp>z<pp>=[(——+—<pp>:<pp>)l + —<pp>+—<pp>-<pp>]

(1.20)
This closure satisfies all the projection and symmetry properties of the exact orientation
dyadic. The predictions made by this model are found to be accurate near the isotropic

state.

1.4.3 The uncoupling approximation

The uncoupling approximation which was first used by Doi (1981), is
<pppp>z<pp>=<pp>(<pp>z<22>). (1.21)
This is one of the simplest closures and has been widely used to describe LCPs. Here the
fourth order moment is expressed as the product of the first order dyads. For simple shear

flows, this approximation fails to show the director tumbling and wagging transitions

predicted by the unapproximated Doi theory (see Feng et al. ,(1998)).

12

1.4.4 The quadratic closure of Hinch and Leal (HL1 closure)
This closure was proposed by Hinch and Leal (1975,1976). Hinch and Leal
developed a closure for the orientation tetrad by interpolating between two regimes: near

equilibrium (weak flows ) flows and strong-flows.

<pppp>z<pp> =§<EB>+%(6<BB>'<BB>'<EB>‘2<EB>‘<EB>
-<pp><pp>:<pp>+§~<pp><pp>z£ +21<pp>z<pp>
2 2
-g£<22>=£-21<22>-<22>=<22>+§£<BB>°<22>=D

(1.22)

The accuracy of the closure varies considerably when applied to the excluded volume

term < pp >:< pppp > and the flow term §2<pppp>, where g is the strain rate

dyadic. The quadratic closure is found to be more accurate than the HL1 closure for the
nematic term, while the HL1 closure is found to be more accurate than the quadratic
closure for the flow term.
1.4.5 The fully symmetric quadratic ( FSQ-) closure

This closure was proposed by Petty et a1. (1999). This closure has a linear and a

quadratic part, and is defined as follows

<BBBB >= C1<EEBB >1+ C2<BEBB >2 (123)
< ppp> < >—[(——+-1-< p>< >)I+—3—< >+fi< > <p >]
E--— 1 BB 35 7 p ‘31—) = 35 BB 7 BB —-9 (1.24)

13

2
<pppp>2<pp>=[(3§<pp><pp>——tr(<pp> <pp> <p_p>))l+
39 2 6
(gs—<pp><pp>)<pp>——7-<pp> <pp>+7<pp> <pp>o<pp>]
(1.25)

The linear part of the closure i.e., < pppp >l is identical to the Hand closure. The closure
retains the six fold symmetry and projection properties of the exact orientation tetradic.
The tetradics < pppp >l and < pppp >2 , individually satisfy the six fold symmetry and

projection properties. The projection property for Eq. (1.23) is recovered by requiring

C1+C2=1. (1.26)

At the nematic state C2= 1/3 ( see Appendix D). Eq. (1.23) can be written as,

< >< >-2[(——+— >< >)I+—§—< >+i< > < >]
2322 BB 3 35 7 W 22 = 35 BE 7 ER BB
2 2
3K? pp><pp>—7tr(<pp> <pp> <pp>))l+
39 2 6
(E<pp><pp>)<pp>—7<pp> <pp>+7<pp> <pp> <pp>]

The FSQ-closure satisfies the six fold symmetry and contraction properties discussed in
Section 1.4.1.

14

1.5 Objectives

Equation (1.15) determines the evolution of the orientation dyadic. Solving the

differential equation requires the evaluation of the term < [32 >:< pppp >. This term, can

be evaluated by the various closure models discussed earlier and the FSQ-closure (see
Section 1.4.5).

In this research the following characteristics of Eq. (1.15) are examined:
1) the critical vales of U for the transitions between the liquid crystalline phases ; 2) the
equilibrium microstructure and phase transitions for lyotropic LCPs using the FSQ-
closure ; and 3) the predictions made by the uncoupling approximation and the HL1

closure.

[.6 Methodology
The second order orientation moment satisfies Eq. (1.15) .The F SQ-closure, which

satisfies the above conditions, is used to calculate the fourth order orientation moment

<pppp >. The goal of the simulations is to predict the critical concentrations for the

transitions between the different prolate and isotropic states, and to get the equilibrium
phase diagram. Figure 1.5 shows a typical equilibrium phase diagram with the different
prolate and isotropic states. The arrows “a” and “b” indicate the transition between the
conditionally stable prolate states to the conditionally stable isotropic and stable prolate
states respectively. The other arrows, as marked, show the prolate-isotropic and isotropic-
prolate transition. The equilibrium orientation states depend on the value of U. Below the
Up), the microstructure relaxes to a stable isotropic state. Beyond Uip, the isotropic states

are unstable, while the prolate states are all stable. Between Up; and Up, the isotropic and

15

prolate states are conditionally stable. The faint curve represents the transition between
conditionally stable isotropic and prolate states (i.e, finite perturbation of a microstructure
on the faint line will change the microstructure to either the conditionally stable prolate
state or the isotropic state). As U is increased beyond Up the microstructure gets
increasingly oriented in one direction (i.e, it tends towards the nematic state). However,
Brownian motion of the solvent molecules prevents the dispersed phase from becoming

nematic for U < 00.

1.6.1 Prolate-Isotropic transition Up)
Eq. (1.15) is solved using a fourth order Runge-Kutta algorithm (see Matlab code
in Appendix H) subject to the following conditions for 0 < 1: S 800 :

< > =8 e .
BB 0 —3—3

The temporal step size Ar= 0.01 .
The simulations are started with U=1. The nematic state (initial state) corresponds to
point A on the invariant diagram ( see Figure A], Appendix A). At the end of the
simulation a33(t) is plotted as a firnction of the dimensionless time, 1:. This gives the
visual conformation that the orientation state is fully relaxed to a steady state. At the

isotropic state the orientation dyadic has the form,

The calculation is repeated for larger values of U. At a particular value of U=U;p ,
the nematic state relaxes to an anisotropic state ( see Figure 1.5). At values of U less than

Up) the nematic state relaxes to a stable isotropic state. This is the region of stable

l6

isotropic states. For U 2 Up, the nematic state relaxes to a stable prolate state. The
prolate state has the form

<BE>= 3119.191 +a22929.2 +a33§3§-3'

where a“ 2 a22 < a”.

1.6.2 Isotropic-prolate transition,Uip

For U > Uip, the isotopic state is unstable. To locate this point for a particular
value of U beyond Uip, the simulation is started with an initial state very close to the
isotropic state, with a step size of 0.001 for a few hundred time steps. If the
microstructure tends to relax back to the isotropic state, then the calculation is repeated
with a slightly higher value of U. This is continued until U=Uip , beyond which the
microstructure tends to relax to the prolate state from an initial state close to the isotropic

state as well as from the nematic state.

1.6.3 Region of multiple steady states Up; < U < Uip

To identify the region of multiple steady states, the simulations are run similar to
the ones developed to identify Uip. An initial state close to the isotropic state is selected.
As illustrated in Figure 1.5, if the initial state is at Point “a” for a value of U between Uip
and Upi , then it relaxes to the isotropic state. If it is at Point “b”, then it relaxes to the
prolate state. The simulations are repeated till the points “a” and “b” merge to yield the

conditionally stable prolate and conditionally stable isptropic states between Up) and Up .

17

l. 7 Scope

C2 is the coefficient of the quadratic terms in the F SQ-closure (see Eq. (1.23) ).
This research examines the influence of C2 on the predictions made by the F SQ-closure.
C2 is equal to 1/3 at the nematic state. Two other values of C2 = 0.29 and 0.31 are used
for this purpose. The relaxation of the microstructure from an initial state is studied by
varying U between 0 and 15. U= 0 corresponds to relaxation due to Brownian motion. As
the value of U increases the contribution of the excluded volume increases. The

relaxation times for different values of U are compared.

18

1.8 Figures

19

 

Isotropic

 

 

 

I Increasing order

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Nematic Molecular Orientational order
I l
Smectic-A
Positional order normal to layers
Hexatic-B
Various tilted Bond Orientational order
smectic phases
C stal-B
ry Positional order within layers
Crystal-E
Asymmetric Axial site symmetry V

 

 

 

 

 

 

Figure 1.1. Progression of order for liquid crystal phases formed by rod like molecules.

20

\I‘

Figure 1.2a Isotropic phase

Figure 1.2c Smectic phase

 

 

Showing order between layers

21

(mi
N
fit“

Figure 1.2b Nematic phase

I'D

\

 

>Y

 

‘----—-----4 -—--—-—------

2:131 E1 +P2 92 +P3 93

p =(sin0cos 4)) e, + (0030) 92 + (sinOsin 4)) 93

Figure 1.3 Representation of an axisymmetric polymer molecule

22

 

Unit sphere

rotary flux of orientation states in

a material frame of reference

 

 

. d
~l(gv)-nds = —llv(9,¢,t)dA
c dt A
\ J \ Vi ‘1
Y
loss of orientation states from area A per accumulation of orientation states
in A unit time due to fiber/fluid and in area A

fiber/fiber interactions

Figure 1.4 Orientation vector distribution over the unit sphere

23

al.—>1 for U —)oo

 

 

  

 

 

 

 

 

f
1 \ ..........
stable prolate states
OI. A
”pi ........................... IP- transition
' PI- _ _
transition unstable isotropic states
Ila.
v I
0 5 1
U

Pi ip

Figure 1.5 Equilibrium phase diagram indicating the prolate-isotropic transition and the

isotropic-prolate transition. For definition of a see Appendix B.

24

CHAPTER 2

RESULTS

2.1 Steady States and Equilibrium States

The evolution equation for the second order moment ( see Eq. (1.15) ) is

D<pp>

1
Dr =_(<EE>-§£)+U(<EB>'<BE>_<BB>KBBBE>

The following properties of the orientation dyadic < pp > must be preserved :

1) The trace of the orientation dyadic predicted must also be 1 for all 1 >0.

2) The orientation dyadic<pp > is symmetric. For this to be true at all times, the

orientation tetradic < pppp > must also have the symmetry properties (see Eq.(1 . 19)).

The F SQ-closure, which satisfies the above conditions, is used to calculate the fourth

order moment <pppp >. This differential equation is solved as explained in the earlier

chapter. The relaxation of the microstructure is predicted with C2= 0.29, 0.31 and 1/3.

2.1.1 Determining the point of prolate-isotropic transition,Uip with C2 = 1/3

The differential equation, Eq. ( 1.15) is solved with C2 = 1/3 and other conditions
described in Section 1.6. The initial run of simulation is started with U = 1. At this value
of U = 1, the microstructure relaxes to an isotropic state. The simulations are repeated
with a higher value of U and other conditions remaining unchanged, if the microstructure
relaxes to an isotropic state. At Up) = 4.656, the microstructure relaxes to a stable prolate

state. At this point, a=0.2257 and

25

0 .2581 0 0
<pp> = 0 0.2581 0
0 0 0.4838

Up; corresponds to the point of prolate-isotropic transition. 1

2.1.2 Determining the point of isotropic-prolate transition,Uip with C2 = 1/3
Eq. (1.15) is solved with C2 = 1/3 and other conditions described in Section 1.6.
To locate this point Uip = 4.6 , the simulation is started with an initial state or = 0.05. This

corresponds to the orientation dyadic,

0.3167 0 O
<pp>: 0 0.3167 0
0 0 0.3666

The components of < p p > are calculated from or as follows,

 

< > =or+l/2
BE 3'3 3/2

_ 1_<BE>M
<BE>I.I'<BB>2.2=——2—

If the microstructure tends (i.e., or becomes less than 0.05) to relax back to the isotropic
state then the calculation is repeated with a slightly higher value of U. This is continued
until U= Uip , beyond which the microstructure tends (i.e., or becomes greater than 0.05)
to relax to the prolate state form an initial state or=0.05. This is the point of stable-

unstable isotropic transition, Uip.

26

2.1.3 Determining the region of multiple steady states Up) < U < Uip with C2 = 1/3

The purpose of the calculations is to determine an or (conditionally stable prolate
state) for the value of U lying between Up) and Uip. This is demonstrated by calculating or
for a particular value U = 4.7. The calculations are done as follows :

Step 1. With U = 4.7 and or=0.1 as the initial condition, the microstructure tends

to relax to the isotropic state.

Step 2. With U = 4.7 and a=0.2 as the initial condition, the microstructure tends

to relax to the prolate state.

Step 3. With U = 4.7 , 0t is varied between 0.1 and 0.2 as limits. At 0L=0.129, the

microstructure does not show a tendency to relax to the isotropic state or to the
stable prolate state. or=0.129 is taken to be the conditionally stable prolate state
corresponding to U=4.7.

Therefore, this method is repeated for discrete values of U between U... and Up), to
determine the conditionally stable prolate region, as indicated by the dashed lines on

Figure 2.1.

2.2 Relaxation to Equilibrium

As evident from Figure 2.7, the time needed for the microstructure to relax from
an initial nematic state to an equilibrium isotropic state (U < 4.656 with C2=1/3) is
considerably longer for values of U which are closer to Up) = 4.656 than for values of U
which are considerably less than Upi=4.656. This can be attributed to the microstructure
trying to reach an equillibrium between the competing Brownian forces (tending to make
the microstructure isotropic) and the excluded volume potential (tending to align the

microstructure). At U=0 there is no excluded volume potential. The relaxation of the

27

microstructure is entirely due to the Brownian motion. Beyond Uip, when the excluded
volume potential dominates the Brownian forces the relaxation time is relatively short for
the values of U close to Uip.

For U lying between Uip and Up) there are multiple steady states. Figure 2.9 shows
the multiple steady states for U=4.7, which lies between Upi= 4.6556 and Uip =5.0.
0t=0. 129 corresponds to the conditionally stable prolate state. So any microstructure with
or higher or lower than a=0.129 will relax to the stable prolate state (0t=0.2869) and
isotropic state ((1:0) respectively.

Figure 2.10 shows the relaxation pathways for different initial states (OF 1,
015,005) with U=6.0 and C2=1/3. or: 0.15 and 0.05 are initial states, which are close to
the isotropic state. Since U=6 is greater than Uip =5 , small perturbation from the
isotropic states will take them to a stable prolate state. This state can also be reached by

relaxing from an initial nematic state with the same value of U.

28

2.3 Figures

29

1.0 -

 

 

0.8 -
0.6 -
0.4 4
0.2 -
0.0

3

 

Figure 2.1 The isotropic-prolate transition predicted by the FSQ-closure with C2=1/3.

see Appendix E for the numerical values for the graph.

30

0.6 -

 

 

 

 

0.4

0'2 - unstable states

0.0 ' .1 ' '
4.6 4.8 5 5.2 5.4

Figure 2.2 The isotropic-prolate transition predicted by the FSQ-closure for C2=1/3

with the region between Uip (5.0) and Up) (4.626) enlarged

31

 

 

 

0.8 -
0.6 - §
a3, =o.2240
0.4 -
------ 1' """""" UIp =4.9

0.2 - /

0 I I l I 1

3 5 7 9 11 13 15

Figure 2.3 The isotropic-prolate transition predicted by the F SQ-closure with C2=0.3 1.

See Appendix F for the numerical values for the graph.

32

 

0 .2 - I unstable states

 

0 l ‘ l 1 I

4.5 4.7 4.9 5.1 5.3 ' 5.5

U

 

Figure 2.4 The isotropic-prolate transition predicted by the F SQ-closure for C2=0.3]

with the region between Up (4.9) and Up) (4.629) enlarged

33

 

 

 

0.8 '
0.6 ‘
(133:0.2413
0.4 ‘
0.2 1.
0,, =4.83
0 I I I I I I
3 5 7 9 11 13 15

Figure 2.5 The isotropic-prolate transition predicted by the FSQ-closure with C2=0.29.

See Appendix G for the numerical values for the graph.

34

 

 

 

0.4 P
0 2 L unstable states
.-
0 l t l l l l
4.50 4.70 4.90 5.10 5.30 5.50

Figure 2.6 The isotropic-prolate transition predicted by the F SQ-closure for C2=0.29

with the region between Uip (4.83) and Upi (4.606) enlarged.

35

1.00 1
0.90 -
0.80 -
0.70

 

0.60 J U: 6.0
0.50
0 .40 . U: 5.0

 

0-30 ‘ U= 4.65
0.20 «

0.10 - 511:0 U=4-0 U=4.5
0.00 -

0 10 20 30 40 50 60 70 80 90 100 110

U: 4.58

 

 

 

Figure 2.7 or versus 1 (dimensionless time) for the relaxation of LCP microstructure

from an initial nematic state with C2=1/3.

36

1.00

 

 

 

 

 

 

0.90
0.80
0.70 q l.
0.60 - ,3 U: 6.0
0.50 .
(X 0 40 . . __ u: 5.0
0.30 -
U= 4.65
0.20 - flu: 4.58
0.10 3 04.5
”=0 U=4.o
0.00 "' l I
0 10 20 30
T

Figure 2.8 a versus 1 (dimensionless time) for the relaxation of LCP microstructure from

an initial nematic state with C2= 1/3.

37

1.00 r
0.90 ‘
0.80 ‘
0.70 r
0.60

0.60 «l

(X. 0.40 ‘
0.30 ‘ a=0.2869

 

0.20 r
0.10

0.001fi I I I I I I I I I I I I
010 20 30 40 50 60 70 80 90100110120130140

 

a: 0.129

 

 

 

 

 

13

Figure 2.9 or versus 1: ( dimensionless time) for the relaxation of LCP microstructure

from different initial states with C2= 1/3 and U= 4.7.

38

1.00

0.90

0.80 -
0.70 4
0.60 r
0.50 r
0.40 ‘
0.30 r
0.20 r
0.10 -
0.00 r r u .

 

 

 

 

Figure 2.10 or versus 1: ( dimensionless time) for the relaxation of LCP microstructure

from different initial states with C2= 1/3 and U=6.

39

CHAPTER 3

CONCLUSIONS

Eq. (1.15) is solved using the F SQ-closure to evaluate the orientation tetradic, for
different values of C2 and U. By analyzing the results obtained, the following conclusion

are made :

1.The Uip =5.0 (see Kini et al. 2003) predicted by the FSQ-closure with C2=1/3 is the
same as that predicted by Chaubal et al.(1995) with the HL1 closure .

2. The values of Up, up), Up) predicted are dependent on the value of C2 chosen. The
region of transition (i.e., region between Up and Upi) reduces as the value of C2
decreases.

3. The magnitude of U relative to Uip influences the relaxation times of the microstructure
from an initial nematic state to an equilibrium state (see Chapter 2, Section 2.2).

4. The F SQ-closure exhibits the isotropic-prolate transition, qualitatively similar to those
shown by the de-coupling approximation and the HL1 closure. A comparison of the Uip,
up), Up; for the de-coupling approximation, HL1 and FSQ-closures are listed in the Table
3.1. From Figure 3.1 it is evident that with C2=1/3, the predictions made by FSQ-closure
are closer to those made by the HL1 closure, and those with C2=0.29 are closer to the

predictions made by the de—coupling approximation closure.

40

3.] Tables and Figures

41

Table 3.1 Comparison of the Uip, Up), api, values for the isotropic-prolate transition

predicted by the F SQ-closure (C2=1/3) with the HL1, and the de-coupling approximation.

 

 

 

 

Closure Up Up) api Reference
de-coupling 3.00 2.667 0.2500 Bhave et
approximation al. 1993
FSQ 5.00 4.656 0.2258 This research
(C2= 1/ 3)
HL1 5.00 4.898 0.125 . Chaubal et
al. 1995

 

 

 

 

 

 

 

42

Table 3.2 The U;,,, Upi, up), values for the isotropic-prolate transition predicted by the

FSQ-closure with different C2’s.

 

 

 

 

C2 Uip Upi up;
0.2900 4.83 4.606 0.2413
0.3100 4.90 4.629 0.2249

1/3 5.00 4.656 0.2258

 

 

 

 

 

43

 

 

 

 

0.8 -
0.6 - , -
OI. /:/ + FSQ(C2=1/3)
0.4 - / ”3'" HL1
De-coupling
+FS C2=0.29
0.2 3 Q( )
0 I I I I I I I I
0 2 4 6 8 10 12 14 16

Figure 3.1 Comparison of the isotropic-prolate transition predicted by the F SQ-

closure(C2=1/3 and 0.29) with the HL1 and the de-coupling approximation.

44

CHAPTER 4

RECOMIVIENDATIONS

Experiments

The F SQ closure predicts an equilibrium phase diagram, which is qualitatively
similar to those made by other closures like the HL1 and the de-coupling approximation.
From Figure 3.1 it can be inferred that, the predictions made by the FSQ-closure are
influenced by the particular value of C2 chosen. These conclusions are based on the
results of the theoretical results. Comparing these results with experimental data on the
relaxation phenomena of lyotropic LCPs, will not only help in judging the predictions
made by the FSQ-closure, but will also provides a basis for the selection of a C2. The
experimental data will also help in evaluating some of the underlying assumptions of the
theory that DR, the rotary diffusion coefficient and C2 are constant.

Aharoni and Walsh (1979) observed the isotropic to anisotropic transition for
liquid crystalline solutions of two, Poly (isocyanates). Beyond a critical concentration
they observed that the solution separated into an isotropic and an anisotropic phase
(liquid crystalline phase). The optical microscopy with cross-polarized light of the
anisotropic phase showed that the isocyanate polymer was being aligned in a particular
direction. Though these experiments were done with lyotropic LCPs that undergo phase
separation, similar experimental studies with LCPs that do not phase separate should be
able to determine the critical concentrations for phase transitions and provide
experimental evidence for fixing a particular value for C2.

The nematic strength U, is dependent on the number density of the molecules or

the concentration of the molecules. So devising experiments to observe and characterize

45

the microstructure with different values of U, should yield an equilibrium phase diagram,
which could then be compared with those deduced from theoretical models. Measuring
the relaxation time for the limit of U=0, would provide an estimate of DR the rotary
diffusion co-efficient. The numerical results indicate that starting from an initial nematic
state beyond Uip, the equilibrium states all lie on the prolate line (see. Appendix .B).
Characterizing the experimentally observed microstructure should help in understanding

this theoretical observation.
Variation of C2

For this research the values of C2 used were 1/3, 0.31 and 0.29. Further

simulations could be carried out with values of C2 greater than 1/3.

46

APPENDIX A : INVARIANTS OF THE ORIENTATION DYADIC

The orientation dyadic g=<p p > has three non-negative eigenvalues a1, a2, and

a3 . Also associated with the tensor are the three invariants. Ia , IIa and 1113 .

13,: 0(2): al +a2 +a3 =1 (A.l)
Ha: tr(gog)= (61.)2 + ((12)2 +013)2 ‘ (A2)
111,: my; .3): (22,)3 + (4,)3 + (a3)3 . (A3)

The structure tensor b is defined as,

1:): g — 1/3 (A.4)
Because g=gT and 1:? it follows that =_b=2T .
Similarly, the invariants of 2 can be written as
Ib= tr(g)=bl +b2 +b3 =0 (A5)
[1,: egg): (6,)2 +(b3)2 +(b,)2 (A.6)
IIIb= tr(2.2.2)= (b,)3 + (b2)3 + (b3)3 . (A.7)

where b1, b2, b3 are the three eigenvalues of I; .

47

The eigenvalues of a__ and 2 are related by

bi = a - 1/3 . (A-8)

Therefore the nontrivial invariants of bin terms of the eigenvalues of a are given by

11.: (al —1/3)2 +(a2 -1/3)2 +(a3 —1/3)2 (A9)

111,: (al —1/3)3 +(a2 —1/3)3 +(a3 -1/3)3 . (A.10)

Figure A] Identifies the set of realizable orientation states (i.e., O S ai $1) .

Table A-1 defines the boundary of the invariant diagram (see Petty, et al. 1999; Nguyen,

20013

48

Table A-1. Invariants of the structure tensor

 

Orientation State

Invariants of g

Eigenvalues of < pp >

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ilb 111;, l 2 3 Notes
A. Nematic 2/3 2/9 0 0 1
B. Planar Anisotropic IIb=2/9 + 2 III), 0 a l-a 0< a<1/2
C. Smectic 1/6 -1/36 0 1/2 1/2
D. Oblate 11.,=6(- rub/6)”3 1-2a a a 1/3< a<1/2
E. Isotropic 0 0 1/3 1/3 1/3
F. Prolate IIb=6( IIIb/6)2’3 a a l-2a o< a<1/3

 

49

 

IIb I

 

1115

Figure A] Invariant diagram of the structure tensor 2.

A. Nematic State

B. Planar Anisotropic States
C. Planar Isotropic State

D. Planar Axisymmetric

E. Isotropic States

F. Axial Axisymmetric

50

APPENDIX B : SCALAR ORDER PARAMETER FOR THE PROLATE STATES
The orientation order parameter on is a scalar measure of the orientation of the

microstructure. or varies between the values 0 (isotropic) and 1 (nematic) and is defined

 

 

 

 

by:
3
or: (--b:b) . (8.1)
2::
From Eq.(A.4), it follows that
3( I)( 1)
or: — a—=—: a-;
2 = 3 = 3
=\/2[a:a—2+ll I]
2 = = 3 9 = =
3 1
)3... 3.
3 1
a=\/§(alz+a22+a32)—§ . (B.2)

At the prolate states the three eigenvalues of a are related as,

al =a2 =(1-a3)/2 .

Therefore Eq. (B.2) can be written in terms of a3 as follows

51

 

 

_ 2 1-..,2 2_1
a—J2[2[ 2 )+a, 3] . (3.3)

0t=—a3 —— . (B.4)

a=1 represents all the molecules perfectly aligned in one direction, while 0t=0

corresponds to an isotropic state.

52

APPENDIX C: TETRADIC CLOSURE AT THE NEMATIC STATE
C2 is the coefficient corresponding to the quadratic term in the FSQ-closure.
The following calculations show that at the nematic state C2= 1/3.

At the nematic state

<m>=(I-C2)<m>l + C2 <fl>2 . (Cl)

1 1
< > s——S II +—S I < > C2
2222 1 35 [2:1 7 I: 22 1 < >

2 2
< pppp >2: E < pp >:< pp > S[ll]+ S[< pp >< pp >]—-7-S[I(< pp > - < pp >)]
(C.3)

S[é,§]=AijBk,+AikBj, +A3IBjk +Aleij+A B +A33Bil

ji ik

(C4)

The FSQ-closure at the nematic state reduces to,

1 1
< pipjpkp1 >=—§(1-3C2)(Iijlkl +likljl +Iillkj)+ 7(1_3C2)(Iij§j§3 +Ii3§j§3 +Ii393-e-j
”1.19393 +I3le39k +1333 9391)+3C2§393§393

(C.5)

Clearly C2 = 1/3 yields the orientation state tetradic g3 g3 93 e3 , which is representative of
the nematic state

53

APPENDIX D: NUMERICAL DATA with C2=1/3
The fourth order tetradic is computed with C2=l/3 for the F SQ-closure and the

differential equation Eq.(1.15) is solved for different values of U.

Table D.1 The eigenvalues of g , the structure parameter and the invariants of 2 for U

less than Up; and with C2=1/3

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a2 a3 0: "0 "1:,
3.00 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.00 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.20 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.40 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000

 

 

For U < Up) , the microstructure relaxes to an isotropic state. The orientation dyadic is a

diagonal matrix. Hence an=a,, a22= a2 and a33= a3 ,Where a,, a2 and a; are the

eigenvalues of a_ . llb, Illb are the second and third invariants of the anisotropic

orientation tensor 0. The isotropic state corresponds to point B on the invariant diagram

(see Appendix A, Figure A-l).

54

Table D2 The eigenvalues of 2:1 , the structure parameter or and the invariants of 1:) for U

greater than Up; and with C2=1/3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a2 a3 a "IL_ 1111,
4.659 0.2550 0.2550 0.4900 0.2350 0.0368 0.0029
4.660 0.2542 0.2542 0.4916 0.2374 0.0376 0.0030
4.700 0.2377 0.2377 0.5246 0.2869 0.0549 0.0052
4.730 0.2304 0.2304 0.5392 0.3088 0.0636 0.0065
5.000 0.1944 0.1944 0.6111 0.4167 0.1158 0.0161
5.500 0.1604 0.1604 0.6792 0.5188 0.1794 0.0310
6.000 0.1389 0.1389 0.7222 0.5833 0.2268 0.0441
6.500 0.1232 0.1232 0.7535 0.6303 0.2649 0.0556
7.000 0.1111 0.1111 0.7778 0.6667 0.2963 0.0659
8.000 0.0932 0.0932 0.8136 0.7204 0.3460 0.0831
9.000 0.0805 0.0805 0.8390 0.7585 0.3835 0.0970
10.000 0.0709 0.0709 0.8581 0.7872 0.4131 0.1084
11.500 0.0603 0.0603 0.8794 0.8191 0.4473 0.1221
13.500 0.0450 0.0450 0.9100 0.8650 0.4988 0.1438
15.000 0.0447 0.0447 0.9105 0.8658 0.4997 0.1442

 

 

 

 

 

 

 

 

For U > Uip, the microstructure relaxes to a prolate state. The orientation dyadic is a

diagonal matrix. The prolate states correspond to the line F on the invariant diagram (see

Appendix A, Figure A-l). At the prolate state a H = a 32 < a 33 .

55

Table D.3 The eigenvalues of g , the structure parameter and the invariants of 2 for

Up; S US Uip and With C2=1/3.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a2 a3 a 11!— "lb
4.70 0.2903 0.2903 0.4193 0.1290 0.01 1 1 0.0005
4.75 0.3013 0.3013 0.3973 0.0960 0.0061 0.0002
4.80 0.3133 0.3133 0.3733 0.0600 0.0024 0.0000
4.85 0.3183 0.3183 0.3633 0.0450 0.0014 0.0000
4.90 0.3217 0.3217 0.3567 0.0350 0.0008 0.0000
5.00 0.3333 0.3333 0.3334 0.0001 0.0000 0.0000 '

 

The values of O. tabulated correspond to the unstable states on Figure 2.2. These are the

unstable states or the conditionally stable prolate states.

56

 

APPENDIX E: NUMERICAL DATA with C2= 0.31
The fourth order tetradic is computed with C2=0.3] for the FSQ-closure and Eq.(l.15) is

solved for different values of U.

Table E.1 The eigenvalues of g , the structure parameter and the invariants of 1:) for U

less than Upi and with C2=0.31

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a} a3 a. “1,—— “lb
4.000 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.500 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.626 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.627 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000

 

For U < Up. , the microstructure relaxes to an isotropic state. The orientation dyadic is a
diagonal matrix. Hence an=a,, an: a; and a33= a3 .Where a,, a2 and a; are the

a . llb, ”ID are the second and third invariants of the anisotropic

eigenvalues of

orientation tensor 2. The isotropic state corresponds to point B on the invariant diagram

(see Appendix A, Figure A-l).

57

Table 13.2 The eigenvalues of g , the structure parameter and the invariants of 2 for U

greater than Up; and with C2=0.3]

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a; a3 a "9._ "lb
4.630 0.2561 0.2561 0.4877 0.2316 0.0358 0.0028
4.631 0.2539 0.2539 0.4922 0.2383 0.0379 0.0030
4.633 0.2512 0.2512 0.4975 0.2463 0.0404 0.0033
4.635 0.2493 0.2493 0.5014 0.2521 0.0424 0.0036
4.650 0.2405 0.2405 0.5190 0.2785 0.0517 0.0048
4.660 0.2365 0.2365 0.5270 0.2905 0.0563 0.0054
4.800 0.2070 0.2070 0.5861 0.3792 0.0958 0.0121
5.200 0.1678 0.1678 0.6644 0.4966 0.1644 0.0272
5.500 0.1495 0.1495 0.7010 0.5515 0.2028 0.0373
6.000 0.1275 0.1275 0.7449 0.6174 0.2541 0.0523
7.000 0.0990 0.0990 0.8020 0.7030 0.3295 0.0772
8.000 0.0806 0.0806 0.8388 0.7582 0.3832 0.0969
9.000 0.0675 0.0675 0.8650 0.7975 0.4240 0.1 127
10.000 0.0576 0.0576 0.8848 0.8272 0.4562 0.1258
15.000 0.0305 0.0305 0.9389 0.9084 0.9389 0.9084

 

 

 

 

 

 

 

 

For U > Uip , the microstructure relaxes to a prolate state. The orientation dyadic is a
diagonal matrix. The prolate states correspond to the line F on the invariant diagram (see

Appendix A, Figure A-l). At the prolate state a I = a 2 < a 3 .

58

Table 13.3 The eigenvalues of g , the structure parameter and the invariants of 1:) for

Upi S U S Uip and with C2=0.31.

 

 

 

 

 

 

 

 

 

4.90

 

 

 

 

 

 

U a1 a; a3 a "b "lb
4.63 0.2667 0.2667 0.4667 0.2000 0.0267 0.0018
4.66 0.2767 0.2767 0.4467 0.1700 0.0193 0.0011
4.70 0.2917 0.2917 0.4167 0.1250 0.0104 0.0004
4.76 0.3077 0.3077 0.3847 0.0770 0.0040 0.0001
4.80 0.3167 0.3167 0.3667 0.0500 0.0017 0.0000

0.3333 0.3333 0.3334 0.0001 0.0000 0.0000

 

The values of on tabulated correspond to the unstable states on Figure 2.4. These are the

unstable states or the conditionally stable prolate states.

59

 

APPENDIX F: NUMERICAL DATA with C2 = 0.29
The fourth order tetradic is computed with C2=0.29 for the FSQ-closure Eq.(l.15) is

solved for different values of U.

Table RI The eigenvalues of g , the structure parameter and the invariants of 1:) for U

less than Up. and with C2=0.29

 

 

 

 

 

 

U a1 3L a3 a "2_ "1:,
2.500 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
3.500 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.200 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.600 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000
4.605 0.3333 0.3333 0.3334 0.0000 0.0000 0.0000

 

 

 

 

 

 

 

 

 

For U < Up, , the microstructure relaxes to an isotropic state. The orientation dyadic is a

diagonal matrix. Hence an=al, 322: a; and a33= a3 _Where a1, a; and a3 are the

eigenvalues of i . Uh, Uh, are the second and third invariants of the anisotropic

orientation tensor b. The isotropic state corresponds to point B on the invariant diagram

(see Appendix A, Figure A-l).

60

Table F2 The eigenvalues of g , the structure parameter and the invariants of 2 for U

greater than Up, and with C2=0.29

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

U a 1 a; a3 a "L... "13,
4.607 0.2495 0.2495 0.5010 0.2515 0.0422 0.0035
4.608 0.2477 0.2477 0.5046 0.2569 0.0440 0.0038
4.610 0.2452 0.2452 0.5096 0.2644 0.0466 0.0041
4.620 0.2379 0.2379 0.5242 0.2863 0.0546 0.0052
4.650 0.2259 0.2259 0.5483 0.3225 . 0.0693 0.0075
4.660 0.2229 0.2229 0.5542 0.3313 0.0732 0.0081
4.670 0.2202 0.2202 0.5595 0.3393 0.0767 0.0087
4.690 0.2155 0.2155 0.5691 0.3537 0.0834 0.0098
4.800 0.1964 0.1964 0.6072 0.4108 0.1125 0.0154
5.500 0.1389 0.1389 0.7222 0.5833 0.2268 0.0441
7.000 0.0872 0.0872 0.8256 0.7384 0.3635 0.0895
8.000 0.0683 0.0683 0.8634 0.7951 0.4215 0.1117
10.000 0.0447 0.0447 0.9107 0.8661 0.5000 0.1444
11.500 0.0333 0.0333 0.9334 0.9001 0.5401 0.1621
13.500 0.0200 0.0200 0.9600 0.9400 0.5891 0.1846
15.000 0.0168 0.0168 0.9664 0.9496 0.6012 0.1903

 

 

 

 

 

 

 

 

For U > Uip , the microstructure relaxes to a prolate state. The orientation dyadic is a

diagonal matrix. The prolate states correspond to the line F on the invariant diagram (see

Appendix A, Figure A-l). At the prolate state a l = a 2 < a 3 .

61

Table F3 The eigenvalues of g , the structure parameter and the invariants of 2 for

Upi 5 Us Up and with C2=0.29

 

 

 

 

 

 

 

 

 

 

 

 

 

U a1 a; 33 a "9— "lb
4.620 0.2700 0.2700 0.4600 0.1900 0.0241 0.0015
4.660 0.2850 0.2850 0.4300 0.1450 0.0140 0.0007
4.680 0.2900 0.2900 0.4200 0.1300 0.0113 0.0005
4.760 0.3133 0.3133 0.3733 0.0600 0.0024 0.0000
4.830 0.3333 0.3333 0.3334 0.0004 0.0000 0.0000

 

The values of a tabulated correspond to the unstable states on Figure 2.4. These are

the unstable states or the conditionally stable prolate states

62

 

APPENDIX G: COMPUTER CODE (see. Nguyen, 2001)

The following is a Matlab code that solves the Eq. (1.15) by using a fourth-order Runge

Kutta algorithm. The fourth order tetradic is computed using the FSQ-closure.

List of notations for computer code
Symbol Definition Corresponding notation

used in the computer code

2 = < pp > Orientation dyadic a

< pppp > Orientation tetradic ppppO
2 Anisotropic orientation dyadic b
llb Second invariant of l__) IIb
Illb Third invariant of 1:) IIIb
U Nematic strength U

T Dimensionless time time
1 Identity matrix ID

C1 Coefficient of linear terms in C]

F SQ-closure
C2 Coefficient of quadratic terms in C2

FSQ-closure
dt step size dt

stop number of time steps stop

63

F low chart for computer code

 

Initially specify
dt , stop , U, a , Cl, n=1

 

 

 

 

le

Calculate
a,pppp0

l

 

 

 

 

 

Solve
Eq. (1.15) using fourth-order R—K algorithm
Calculate a, b

n=n+1

 

 

 

if n < stop

 

 

 

 

V

if n = stop

64

dt=0.01 ;
stop=90000;

%for liquid crystaline polymers;
lam=l .0;
U=4.626

C1=2/3;
C2=l-Cl;

% %Initial input: perfect random;
a=[ 0 O O
O O O
0 O l ] ;
%Unit tensor;
ID=[ l O 0
O l O
O O 1 ];

s=0;
w=0;

%Defined parameter;
alpha=[-1/351/7O 0 000 O
0 0 1-2/7 000 O
O 0 O O 000 O

0 O 00 001 -2/7];

%Symmetric and Antisymmetn'c portion of Velocity Gradient tensor;
for n=1 :stop

%First loop of the R-K method
%Determine SH 1];
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SII(i,j,k,l)=ID(i,j)*ID(k,1)+ID(i,k)*ID(j,l)+ID(i,l)*ID(k,j);
end
end

65

end
end

%Determine S[I a];
for i=1 :3
for j=l :3
for k=1 :3
for l=l :3
SIa(i,j,k,l)=ID(i,j)*a(k,l)+ID(i,k)*a(j,1)+ID(i,l)*a(k,j);
SIa(i,j,k,l)=a(i,j)*ID(k,l)+a(i,k)*IDG,1)+a(i,l)*H)(kJ)+SIa(iJ,k,l);
end
end
end
end

%Determine S[a a];
for i=1 :3
for j=l :3
for k=1 :3
for 1:1 :3
SaafiJ,Kl)=a(i.i)*a(k,1)+a(i,k)*a(i,1)+a(i,l)*a(k.i);
end
end
end
end

%Determine S[I(a*a)];
for i=1 :3
for j=1 :3
aa(i.j)=0;
bb(i,j)=0;
for k=1 :3
aa(i.i) = aa(i.i)+a(i,k)*a(k,j);
bb(i,j) = bb(i,i)+(a(i,k)-ID(i,k)/3)*(a(k,j)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=l :3
for k=1 :3
for l=1 :3
SIaa(i,j,k,l)=ID(i,i)*aa(k,l)+ID(i,k)*aa(j,1)+ID(i,l)*aa(k,j);
SIaa(i,j,k,l)=aa(i,j)*ID(k,l)+aa(i,k)*ID(j,l)+aa(i,l)*ID(k,j)+SIaa(i,i,k,l);
end

66

end
end
end

%Determine S[(a*a)a];
for i=1 :3
for j=l :3
aa(i,j)=0;
for k= 1 :3
aa(i.i)=aa(i.i)+a(i,k)*a(k.i);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for 1=1 :3
Saaa(i,j,k,l)=aa(i,j)*a(k,l)+aa(i,k)*a0,l)+aa(i,l)*a(k,j);
Saaa(i,j,k,l)=a(i,i)*aa(k,l)+a(i,k)*aa(j,l)+a(i,l)*aa(k,j)+Saaa(i,j,k,1);
end
end
end
end

%Determine SIaaa;
for i=1 :3
for j=l :3
aaa(i.i)=0;
bbb(i,j)=0;
for k=1 :3
aaa(i.i)=aaa(iJ)+aa(i,k)*a(k.i);
bbb(i,i)=bbb(i,j)+bb(i,k)*(a(k,j)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
SIaaa(i,j,k,l)=ID(i,j)*aaa(k,l)+H)(i,k)*aaa(j,l)+ID(i,l)*aaa(k,j);
SIaaa(i,j,k,l) = aaa(i,j)*ID(k,l)+aaa(i,k)*ID(j,1)+aaa(i,l)*ID(k,j)+SIaaa(i,j,k,l);
end
end

67

end
end

%Determine S[(a*a)(a*a)];
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
Saaaa(i,j,k,l)=aa(i,j)*aa(k,l)+aa(i,k)*aa(j,l)+aa(i,l)*aa(k,j);
end
end
end
end

%Determine S[I(a*a*a*a)];
for i =1 :3
for j=l :3
aaaa(i,j)=0;
for k=1 :3
aaaa(i.i)=aaaa(i.i)+aaa(i.k)*a(k.i);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for 1:] :3
SIaaaa(i,j,k,1)=ID(i,j)*aaaa(k,l)+ID(i,k)*aaaa(i,l)+ID(i,l)*aaaa(k,j);

SIaaaa(i,j,k,l)=aaaa(i,i)*ID(k,l)+aaaa(i,k)*ID(i,l)+aaaa(i,l)*ID(k,j)+SIaaaa(i,j,k,l);
end
end
end
end

%Determine aza (11a) and alpha(2,1);
IIa=aa( l , l )+aa(2,2)+aa(3,3);
IIb=bb(] ,1 )+bb(2,2)+bb(3,3);
alpha(2,1)=2*IIa/3 5;

%Determine alpha(3,1);

ALP1=aaa(1,1)+aaa(2,2)+aaa(3,3);
alpha(3,1)=1/35+(4/35)*(ALP1/Ila);

68

111a=aaa(1,1)+aaa(2,2)+aaa(3,3);
IIIb=bbb(l,l)+bbb(2,2)+bbb(3,3);

%Determine alpha(3,4);
alpha(3 ,4)=- l / (7 * Ila);

%Determine alpha(3,5);
alpha(3,5)=l/IIa;

%Determine alpha(3,6);
alpha(3,6)=-4/(7*Ha);

%Determine alpha(4,l);
ALP2=aaaa( l , l )+aaaa(2,2)+aaaa(3,3);
alpha(4,1)=2/35*ALP2+(l/35)*IIa*IIa;

%Determine alpha(4,4);
alpha(4,4)=(- l/7)* IIa;

%MSU contribution;
%Determine pppp1,pppp2
for i=1 :3
for j=1 :3
for k=1 :3
for 1:1 :3
pppp5(i.j,k,l)=a(i.i)*a(k,l);
ppppl(i,j,k,l)=alpha(l,l)*SII(i,i,k,l)+alpha(1,2)*SIa(i,j,k,l);

%pppp2(i,i,k,l)=alpha(2,1)*SH(i,j,k,l)+alpha(2,2)*SIa(i,j,k,l)+alpha(2,3)*Saa(i,j,k,l)+alp
ha(2,4)*SIaa(i,j,k,l);

pppp2(i,j,k,l)=alpha(2,1)*SII(i,j,k,l)+alpha(2,3)*Saa(i,j,k,l)+alpha(2,4)*SIaa(i,j,k,l);

%Casel :Michigan State University model;
pppp0(i.i.k,l)=C1*pppp1(iJ,k,l)+C2*pppp2(i.i,k,l);

end
end
end
end

%Determine <pppp>:S

for i=l:3
forj=l:3

69

M(i.i)=0;

for k=1 :3
for l=l :3

M(i.i)=M(i.j)+pppp0(i.j,k,l)*s(l,k);

end

end

end
end

%Determine <pp>:<pppp>
for i=1 :3
for j=l :3
N(i.i)=0;
for k=1 :3
for l=l :3
N(i,j)=N(i,j)+ppppo(i.i,k,l)*a(l,k);
end
end
end
end

%Determine S.a+a.S, W.a+a.W(T);
for i=1 :3
for j=l :3
saas(i,j)=0;
waaw(i,j)=0;
for k=1 :3
saas(i,j)=saas(i,j)+s(i,k)*a(k,j)+a(i,k)*s(k,j);
waaw(i,j)=waaw(i,j)+w(i,k)*a(k,j)+a(i,k)*w(i,k);
end
end
end

%Right handside of the equation
for i=1 :3
for j=1 :3
dkl(i,j)=(ID(i,j)/3-a(i,j))+lam*Pe*(saas(i,i)-2*M(i,j))-Pe*waaw(i,j)+U*(aa(i,j)-
N(iJ));
end
end

%an loop of 4th order R_K
for i=1 :3

for j=l :3
a1(i,j)=a(i,j)+dkl(i,j)*dt*0.5;

70

 

end
end

%Determine S[I a];
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIal(i,j,k,l)=ID(i,j)*a1(k,l)+ID(i,k)*al(j,l)+ID(i,l)*al(k,j);
SIal(i,j,k,l)=al(i,j)*ID(k,l)+al(i,k)*ID(i,l)+al(i,l)*ID(k,j)+SIal(i,j,k,l);
end
end
end
end

%Determine S[a a];
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
Saal(i,j,k,l)=al(iJ)*al(k,l)+al(i,k)*al(j,l)+a1(i,l)*al(k,j);
end
end
end
end

%Determine S[I(a*a)];
for i=1 :3
for j=1 :3
aa1(i,j)=0;
bbl(i,j)=0;
for k=1 :3
aal(i,j) = aa1(i,j)+a1(i,k)*al(k,j);
bbl(i,j) = bb1(i,j)+(al(i,k)-ID(i,l<)/3)*(al(kJ)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for 1=1 :3
SIaal(i,j,k,l)=ID(i,j)*aal(k,l)+ID(i,k)*aa1(j,l)+ID(i,l)*aal(k,j);
SIaal(i,j,k,l)=aal(i,j)*ID(k,l)+aal(i,k)*ID(j,l)+aa1(i,l)*ID(k,j)+SIaal(iJ,k,l);
end
end

71

end
end

%Determine S[(a*a)a];
for i=1 :3
for j=1 :3
aal(iJ)=0;
for k=1 :3
aa1(i,j)=aal(i,j)+al (i,k)*al(k,j);
end
end
end

for i=1 :3
for j=l :3
for k=1 :3
for l=l :3
Saaal(i,j,k,l)=aa1(i,i)*al(k,l)+aa1(i,k)*al(i,l)+aal(i,l)*al(k,j);
Saaal(i,j,k,l)=al(i,j)*aal(k,l)+a1(i,k)*aal(i,l)+al(i,l)*aa1(kJ)+Saaal(i,j,k,l);
end
end
end
end

%Determine SIaaa;
for i=1 :3
for j=1 :3
aaal (i,j)=0;
bbbl(i,j)=0;
for k=1 :3
aaal(i,j)=aaal(iJ)+aa1(i,k)*al(k,i);
bbb1(i,j)=bbbl(iJ)+bb1(i,k)*(a1(kJ)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=l :3
for k=1 :3
for 1=1 :3
SIaaal(i,i,k,l)=ID(i,i)*aaal(k,l)+ID(i,k)*aaa1(j,l)+ID(i,l)*aaa1(k,i);
SIaaal(i,j,k,l) =
aaal(i,j)*ID(k,l)+aaal (i,k)*ID(j,l)+aaal(i,l)*ID(k,j)+SIaaal(i,j,k,1);
end
end

72

end
end

%Determine S[(a*a)(a*a)];
for i=1 :3
for j=l :3
for k=1 :3
for l=1 :3
Saaaal(i,j,k,l)=aal(i,i)*aa1(k,l)+aal (i,k)*aa1(j,l)+aal(i,l)*aa1(k,j);
end
end
end
end

%Determine S[I(a*a*a*a)];
for i =1 :3
for j=1 :3
aaaal(i,j)=0;
for k=1 :3
aaaal(i,j)=aaaa1(i,j)+aaal (i,k)*al (1g);
end
end
end

for i=1 :3
for j=l :3
for k=1 :3
for l=1 :3
SIaaaal(i,j,k,l)=ID(i,j)*aaaal(k,l)+ID(i,k)*aaaal(j,l)+ID(i,l)*aaaal(k,j);

SIaaaal(i,j,k,l)=aaaa1 (i,j)*ID(k,l)+aaaa1(i,k)*ID(j,l)+aaaa1(i,l)*ID(k,j)+SIaaaal(i,j,k,l);
end
end
end
end
%Determine aza (11a) and alpha(2,l);
IIa1=aal(1,1)+aal(2,2)+aal(3,3);
Ilbl=bbl( 1,1)+bb1(2,2)+bb1(3,3);
alpha1(2,1)=2*(IIal)/35;

%Determine alpha(3,1);

ALPl l=aaal(l,1)+aaal(2,2)+aaal(3,3);
alphal(3,l)=1/35+(4/3S)*(ALP1l/IIal);
IIIal=aaal(l ,1)+aaal(2,2)+aaal(3,3);
IIIbl=bbb1(1,1)+bbb1(2,2)+bbbl(3,3);

73

%Determine alpha(3,4);
alphal(3,4)=l/(7*IIal);

%Determine alpha(3,5);
alpha1(3,5)=l/IIal;

%Determine alpha(3,6);
a1phal(3,6)=-4/(7*Ilal);

%Determine alpha(4,l);
ALP21=aaaal(l, l )+aaaal(2,2)+aaaa1(3,3);
alphal(4,1)=2/35*ALP21+(1/35)*Ila1*IIal;

%Determine alpha(4,4);
alpha 1 (4,4)=(- 1/7)*IIa1 ;

%MSU contribution;
%Detennine pppp1,ppp92,pppp3,pppp4;
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
pppp51(i.i,k,l)=a1(i.i)*al(k.l);

ppppl 1(i,j,k,l)=alpha(1,l)* SII(i,j,k,l)+alpha(1,2)*SIa1(i,j,k,l);

%pppp2(i,i ,k,l)=alpha(2, 1 )*SII(i,j ,k,l)+alpha(2,2)* SIa(i,j ,k,1)+alpha(2,3)* Saa(i,j ,k,l)+alp

ha(2,4)*SIaa(i,j,k,l);

pppp2 1 (i,j,k,l)=alphal (2, 1)* SII(i,j ,k,l)+alpha(2,3)* Saa l (ij ,k,l)+alpha(2,4)* SIaa1(i,j,k,l)

%Casel :Michigan State University model;
pppp01(i.i.k.l)=C1*pppp1 l(iJ,k.l)+C2*pppp21(iJ,k,1);

end
end
end
end

%Determine <pppp>:S
for i=1 :3
for j=1 :3
M1(iJ)=0;
for k=1 :3
for l=1 :3
M1(i.i)=M1(iJ)+pppp01(iJ,k,l)*s(l,k)
end

74

end
end
end

%Determine <pp>:<pppp>
for i=1 :3
for j=1 :3
Nl(i,j)=0;
for k=1 :3
for l=1 :3
Nl(i.i)=N1(i.i)+pppp01(i.i,k,l)*al(Lk);
end
end
end
end

%Determine S.a+a.S, W.a+a.W(T);
for i=1 :3
for j=1 :3
saasl(i,j)=0;
waawl(i,j)=0;
for k=1 :3
saasl(i,j)=saasl(i,j)+s(i,k)*a1(k,j)+al(i,k)*s(k,j);
waawl(i,j)=waawl(i,j)+w(i,k)*al(k,j)+al(i,k)*w(j,k);
end
end
end

%Right handside of the equation
for i=1 :3
for j=1 :3
dk2(i,j)=(ID(i,j)/3-al(i,j))+lam*Pe*(saas1(i,j)-2*Ml (i,j))-
Pe*waaw1(i,j)+U*(aal(i,j)-N 1(i,j));
end
end

%3rd loop of 4th order R-K

for i=1 :3
for j=1 :3
a2(i,i)=a(i,i)+dk2(i,j)*dt*0.5;
end
end

%Determine S[I a];

75

for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
SIaZGJ,k,l)=ID(i.j)*aZ(k,1)+ID(i,k)*€120,1)+ID(i,l)*aZ(k.i);
SIa2(i,j,k,l)=a2(i,j)*ID(k,l)+a2(i,k)*ID(j,l)+a2(i,l)*ID(k,j)+SIa2(i,j,k,l);
end
end
end
end

%Determine S[a a];
for i=1 :3
for j=l:3
for k=1 :3
for l=l :3
Saa2(i,j,k,l)=a2(i,j)*a2(k,l)+a2(i,k)*a2(j,1)+a2(i,l)*a2(k,j);
end
end
end
end

%Determine S[I(a*a)];
for i=1 :3
for j=1 :3
, aa2(i,i)=0;
bb2(i,j)=0;
for k=1 :3
3320.0 = 332(iJ)+aZ(i,k)*aZ(kJ);
bb2(i,j) = bb2(i,j)+(a2(i,k)-ID(i,k)/3)*(a2(k,j)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for 1=1 :3
SIaa2(i,j,k,l)=ID(i,j)*aa2(k,l)+ID(i,k)*aa2(j,1)+ID(i,l)*aa2(k,j);
SIaa2(i,j,k,l)=aa2(i,j)*ID(k,l)+aa2(i,k)*ID(j,l)+aa2(i,l)*ID(k,j)+SIaa2(i,j,k,l);
end
end
end
end

%Determine S[(a*a)a];

76

for i=1 :3
for j=1 :3
aa2(iJ)=0;
for k=1 :3
aa2(iJ)=aaZ(iJ)+a2(i,k)*a2(kJ);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
Saaa2(i,j,k,l)=aa2(i,j)*a2(k,l)+aa2(i,k)*a2(j,l)+aaZ(i,l)*a2(k,j);
Saaa2(i,j,k,l)=a2(i,j)*aa2(k,l)+a2(i,k)*aa2(j,l)+a2(i,l)*aa2(k,j)+Saaa2(i,j,k,l);
end
end
end
end

%Determine SIaaa;
for i=1 :3
for j=1 :3
aaa2(i.i)=0;
bbb2(i,j)=0;
for k=1 :3
aaaZGJ)=aaaZ(iJ)+aa2(i,k)*a2(kJ);
bbb2(i,j)=bbb2(i,j)+bb2(i,k)*(a2(k,j)-ID(k,j)/3);
end
end
end .

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIaaa2(i,j,k,l)=ID(i,j)*aaa2(k,l)+ID(i,k)*aaa2(j,1)+ID(i,l)*aaa2(k,j);
SIaaa2(i,j,k,l) =
aaa2(i,j)*ID(k,l)+aaa2(i,k)*ID(j,l)+aaa2(i,l)*ID(k,j)+SIaaa2(i,i,k,l);
end
end
end
end

%Determine S[(a*a)(a*a)];

77

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
Saaaa2(i,j,k,l)=aa2(i,j)*aa2(k,l)+aa2(i,k)*aa2(j,l)+aa2(i,l)*aa2(k,j);
end ‘
end
end
end

%Determine S[I(a* a*a*a)];
for i =1 :3
for j=1 :3
aaaa2(i,j)=0;
for k=1 :3
aaaa2(i,j)=aaaa2(i,j)+aaa2(i,k)*a2(k,j);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIaaaa2(i,j,k,l)=ID(i,j)*aaaa2(k,l)+ID(i,k)*aaaa2(j,l)+ID(i,l)*aaaa2(k,j);

SIaaaa2(i,j,k,l)=aaaa2(i,j)*ID(k,l)+aaaa2(i,k)*ID(i,l)+aaaa2(i,l)*ID(k,i)+SIaaaa2(i,j,k,l);
end
end
end
end
%Determine aza (Ha) and alpha(2,l);
11a2=aa2(1,1)+aa2(2,2)+aa2(3,3);
Ilb2=bb2(l ,1)+bb2(2,2)+bb2(3,3);
alpha2(2,l)=2*(IIa2)/35;

%Determine alpha(3,1);

ALP 12=aaa2( l , 1)+aaa2(2,2)+aaa2(3,3);
alpha2(3,1)=l/35+(4/35)*(ALP12/IIa2);
IIIa2=aaa2(] ,1)+aaa2(2,2)+aaa2(3 ,3);
IIIb2=bbb2( 1 , l )+bbb2(2,2)+bbb2(3 ,3);

%Determine alpha(3,4);
alpha2(3,4)=1/(7*IIa2);

78

 

%Determine alpha(3,5);
alpha2(3,5)=l/Ila2;

%Determine alpha(3,6);
alpha2(3,6)=-4/(7*Ila2);

%Determine alpha(4, 1 );
ALP22=aaaa2(1,1)+aaaa2(2,2)+aaaa2(3,3);
alpha2(4, 1 )=2/35 *ALP22+(1/35)*IIa2*IIa2;

%Determine alpha(4,4);
alpha2(4,4)=(-1/7)*Ila2;

%MSU contribution;
%Determine ppppl
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3

pppp52(iJ,k,l)=a2(i.i)*a2(k,l);
pppp12(i,j,k,l)=alpha(l,1)*SII(i,j,k,l)+alpha(1,2)*SIa2(i,j,k,l);

%pppp2(i,i,k,l)=alpha(2,l)*SII(i,i,k,l)+alpha(2,2)*SIa(i,j,k,l)+alpha(2,3)*Saa(i,j,k,l)+alp
ha(2,4)*SIaa(i,j,k,l); _

pppp22(i,j,k,l)=alpha2(2, 1)*SII(iJ,k,l)+alpha(2,3)* Saa2(i,j,k,l)+alpha(2,4)* SIaa2(i,j ,k,l);

pppp32(iJ,k,l)=alpha2(3,1)*SII(i,i,k,l)+alpha2(3,4)*SIaa2(i,j,k,l)+alpha2(3,5)*Saaa2(i,j,k
,l)+alpha2(3,6)*SIaaa2(i,j,k,l);

pppp42(i,j ,k,l)=alpha2(4, l )*SII(i,j ,k,l)+alph32(4,4)*SIaa2(i,j ,k,l)+a1pha(4,7)* Saaaa2(i,j ,k
,l)+alpha(4,8)*SIaaaaZ(i,i,k,l);

%Casel :Michigan State University model;
pppp02(i.i,k,l)=Cl*ppppl2(i.i,k,1)+C2*pppp22(iJ.191);
end
end
end
end

%Determine <pppp>zS
for i=1 :3
for j=1 :3
M2(iJ)=0;

79

for k=1 :3
for l=l :3
M2(i.i)=M2(i,j)+pppp02(iJ,k,l)*S(l,k);
end
end
end
end

%Determine <pp>:<pppp>
for i=1 :3
for j=1 :3
N2(i.i)=0;
for k=1 :3
for l=l :3
N2(iJ)=N2(iJ)+pPPP02(iJ,k,1)*a2(1,k);
end
end
end
end

%Determine S.a+a.S, W.a+a.W(T);
for i=1 :3
for j=1 :3
saa52(i,i)=0;
waaw2(i,j)=0;
for k=1 :3
saasZ(i,j)=saasZ(i,j)+s(i,k)*a2(k,j)+a2(i,k)*s(k,j);
waaw2(i,j)=waaw2(i,j)+w(i,k)*a2(k,j)+a2(i,k)*w(i,k);
end
end
end

%Right handside of the equation
for i=1 :3
for j=1 :3
dk3(i,i)=(ID(i,j)/3-a2(i,j))+lam*Pe*(saas2(i,j)-2*M2(i,j))-
Pe*waaw2(i,j)+U*(aa2(i,j)-N2(i,j));
end
end

%3nd loop of 4th order R-K
for i=1 :3

for j=1 :3
a3(i.il=a(i.i)+dk3(i.j)*dt;

8O

end
end

%Determine S[I a];
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIa3(i,j,k,l)=ID(i,j)*a3(k,l)+ID(i,k)*a3(j,1)+H)(i,l)*a3(k,j);
SIa3(i,j,k,l)=a3(i,i)*ID(k,l)+a3(i,k)*[DO,1)+a3(i,l)*ID(k,j)+SIa3(i,j,k,1);
end
end
end
end

%Determine S[a a];
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
Saa3(i,j,k,l)=a3(i,j)*a3(k,l)+a3(i,k)*a3(j,l)+a3(i,l)*a3(k,j);
end
end
end
end

%Determine S[I(a*a)];
for i=1 :3
for j=1 :3
aa3(i.i)=0;
bb3(i,j)=0;
for k=1 :3
21330.0 = aa3(i.j)+a3(i,k)*a3(kJ);
bb3(i,j) = bb3(ij)+(a3(i,k)-ID(i,k)/3)*(a3(k,j)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for 1=l :3
SIaa3(i,j,k,l)=ID(i,j)*aa3(k,l)+H)(i,k)*aa3(j,l)+ID(i,l)*aa3(k,j);
SIaa3(i,j,k,l)=aa3(i,j)*ID(k,l)+aa3(i,k)*ID(i,l)+aa3(i,1)*ID(k,j)+SIaa3(iJ,k,l);
end
end

81

end
end

%Determine S[(a*a)a];
for i=1 :3
for j=1 :3
aa3(i,j)=0;
for k=1 :3
aa3(i.j)=aa3(i.i)+a3(i,k)*a3(k.i);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
Saaa3(i,j,k,l)=a33(i,j)*a3(k,l)+aa3(i,k)*a3(j,1)+aa3(i,l)*a3(k,j);
Saaa3(i,j,k,l)=a3(i,j)*aa3(k,l)+a3(i,k)*aa3(j,l)+a3(i,l)*aa3(k,j)+Saaa3(i,j,k,l);
end
end
end
end

%Determine SIaaa;
for i=1 :3
for j=1 :3
aaa3(i,j)=0;
bbb3(i,j)=0;
for k=1 :3
aaa3(i,j)=aaa3(i,i)+aa3(i,l<)*a3(k,j);
bbb3(i,j)=bbb3(i,j)+bb3(i,k)*(a3(k,i)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
SIaaa3(i,j,k,l)=ID(i,j)*aaa3(k,l)+ID(i,k)*aaa3(j,1)+1D(i,l)*aaa3(k,j);
SIaaa3(i,j,k,l) =
aaa3(i,j)*ID(k,l)+aaa3(i,k)*ID(j,l)+aaa3(i,l)*ID(k,j)+SIaaa3(i,j,k,l);
end
end

82

end
end

%Determine S[(a*a)(a*a)];
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3 ,
Saaaa3(i,j,k,l)=aa3(i,j)*aa3(k,l)+aa3(i,k)*aa3(j,l)+aa3(i,l)*aa3(k,j);
end
end
end
end

%Determine S[I(a*a*a*a)];
for i =1 :3
for j=1 :3
aaaa3(iJ)=0;
for k=1 :3
aaaa3(i,j)=aaaa3(i,j)+aaa3(i,k)*a3 (k,j);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIaaaa3(i,j,k,l)=ID(i,j)*aaaa3(k,l)+ID(i,k)*aaaa3(i,l)+HD(i,l)*aaaa3(kJ);

SIaaaa3(i,j,k,l)=aaaa3(iJ)*ID(k,l)+aaaa3(i,k)*1D(j,l)+aaaa3(i,l)*ID(k,j)+SIaaaa3(i,j,k,l);
end
end
end
end
%Determine aza (113) and alpha(2,l);
IIa3=aa3(l,1)+aa3(2,2)+aa3(3,3);
IIb3=bb3( 1 , 1 )+bb3(2,2)+bb3(3,3);
alpha3(2,l)=2*(IIa3)/35;

%Determine alpha(3,1);
ALP13=aaa3(1,1)+aaa3(2,2)+aaa3(3,3);
alpha3(3, 1)=1/35+(4/35)*(ALP1 3/IIa3);
IIIa3=aaa3(l,1)+aaa3(2,2)+aaa3(3,3);
IIIb3=bbb3( 1 , 1 )+bbb3 (2,2)+bbb3 (3,3);

83

%Determine alpha(3,4);
alpha3(3,4)=l/(7*Ila3);

%Determine alpha(3,5);
alpha3(3,5)=l/IIa3;

%Determine alpha(3,6);
alpha3(3,6)=-4/(7*Ila3);

%Determine alpha(4,l);
ALP23=aaaa3(l,1)+aaaa3(2,2)+aaaa3(3,3);
alpha3(4,1)=2/35*ALP23+(1/35)*IIa3*IIa3;

%Determine alpha(4,4);
alpha3(4,4)=(-1/7)*IIa3;

%MSU contribution;
%Determine ppppl,pppp2
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
pppp53(i.i,k,l)=a3(i.i)*a3(k.l);
ppppl 3(i,j,k,l)=alpha(l , l )"' SII(i,i,k,l)+alpha( 1 ,2)* SIa3(i,j ,k,l);

%pppp2(i,j,k,l)=alpha(2, l )*SII(i,j ,k,l)+alpha(2,2)*SIa(iJ,k,l)+alpha(2,3)*Saa(i,j,k,l)+alp
ha(2,4)*SIaa(i,j,k,l);

pppp23(i,j,k,l)=alpha3 (2, 1)*SII(i,j,k,l)+alpha(2,3)* Saa3(i,j,k,l)+alpha(2,4)*SIaa3(i,j,k,l);

%Casel :Michigan State University model;
pppp03(i.i,k.l)=C1*ppppl3(i.i,k,l)+C2*pppp23(i.i.k,l);
end
end
end
end

%Determine <pppp>:S
for i=1 :3
for j=1 :3
M3(i,j)=0;
for k=1 :3
for l=l :3
M3(i.i)=M3(i.i)+pppp03(i.i,k,l)*s(l,k);

84

end
end
end
end

%Determine <pp>:<pppp>
for i=1 :3
for j=1 :3
N3(iJ)=0;
for k=1 :3
for l=1 :3
N3(iJ)=N3(iJ)+PPPP03(iJ,k,l)*a3(l,k);
end
end
end
end

%Determine S.a+a.S, W.a+a.W(T);
for i=1 :3
for j=1 :3
saas3(i,j)=0;
waaw3(i,j)=0;
for k=1 :3
saas3(i,j)=saas3(ij)+s(i,k)*a3(k,j)+a3(i,k)*s(k,j);
waaw3(i,j)=waaw3(i,j)+w(i,k)*a3(k,j)+a3(i,k)*w(j,k);
end
end
end

%Right handside of the equation
for i=1 :3
for j=1 :3
dk4(i,j)=(ID(i,j)/3-a3(i,j))+lam*Pe*(saas3(i,j)-2*M3(i,j))-
Pe*waaw3(i,j)+U*(aa3(i,j)-N3(i,j));
end
end

%4nd loop of 4th order R-K

for i=l:3
forj=1:3

%Overall computation

a(i,j)=a(i,j)+(l/6)*(dk1(iJ)+2*dk2(i,j)+2*dk3(i,j)+dk4(i,i))*dt;
E=[a(1,1),a(1,2),a(1,3);a(2,1),a(2,2),a(2,3);a(3,1),a(3,2),a(3,3)];

85

[v,lambdal]=eig(E);

if i=2
if j=3
A(n)=a(iJ);
time(n)=(n- 1 )*dt;
end
end
if i=2
if j==
B(n)=a(i.i);
time(n)=(n- 1 )*dt;
end
end
if i==1
if j== l
C(n)=a(iJ);
time(n)=(n- 1)*dt;
end
end
if i==3
if j=3
D(n)=a(i.i);
time(n)=(n- l )* (it;
end
end

%d(n)=det(E);
e3(n)=v(3,3);
time(n)=(n- 1 )*dt;
e2(n)=v(2,2);
time(n)=(n-1 )*dt;
en(n)=V(2,3);
time(n)=(n- 1 )*dt;
e1(n)=v(1,1);
time(n)=(n- 1 )*dt;
12(n)=lambdal(2,2);
time(n)=(n- 1 )*dt;
l3(n)=lambdal(3,3);
time(n)=(n- 1 )*dt;

- ll(n)=lambdal(l,1);
time(n)=(n- 1 )*dt;

if i==
iszz

86

BB(n)=bb(i,j);
BBB(n)=bbb(i,j);
end
end
if i=2
if j=2
CC(n)=bb(i.j);
CCC(n)=bbb(i,j);
end
end
if i==
if j=3
DD(n)=bb(i,j);
DDD(n)=bbb(i,j);
end
end
end
end

for i=1 :3
for j=1 :3
II(n)=BB(n)+CC(n)+DD(n);
HI(n)=BBB(n)+CCC(n)+DDD(n);
end
end

a
%Computation of finding tau values

%Determine S[I a];
for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIatau(i,j,k,l)=ID(i,j)*a(k,l)+ID(i,k)*a(i,l)+ID(i,l)*a(k,j);
SIatau(i,j,k,l)=a(i,j)*ID(k,l)+a(i,k)*ID(j,1)+a(i,l)*ID(k,j)+SIatau(i,j,k,l);
end
end
end
end

%Determine S[a a];
for i=1 :3
for j=1 :3
for k=1 :3

87

for l=l :3
Saatau(i,j,k,l)=a(i,j)*a(k,l)+a(i,k)*a(j,1)+a(i,l)*a(k,j);
end
end
end
end

%Determine S[I(a*a)];
for i=1 :3
for j=1 :3
aatau(i,j)=0;
bbtau(i,j)=0;
for k=1 :3
aatau(i,j) = aatau(i,j)+a(i,k)*a(k,j);
bbtau(i,j) = bbtau(iJ)+(a(i,k)-ID(i,l()/3)*(a(k,j)-ID(k,j)/3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=l :3
SIaatau(i,j,k,l)=ID(i,j)*aatau(k,l)+ID(i,k)*aatau(j,l)+ID(i,l)*aatau(k,j);

SIaatau(i,j,k,l)=aatau(i,j)*ID(k,l)+aatau(i,k)*ID(j,l)+aatau(i,l)*ID(k,j)+SIaatau(i,j,k,l);
end
end
end
end

%Determine S[(a*a)a];
for i=1 :3
for j=1 :3
aatau(i,j)=0;
for k=1 :3
aatau(i,j)=aatau(i,j)+a(i,k)* 8091);
end
end
end

for i=l:3
forj=lz3
for k=1:3
for l=l:3

88

Saaatau(i;j,k,l)=aatau(i,j)*a(k,l)+aatau(i,k)*a(j,l)+aatau(i,l)*a(k,j);

Saaatau(i,j,k,l)=a(i,j)*aatau(k,l)+a(i,k)*aatau(j,l)+a(i,l)*aatau(k,j)+Saaatau(i,j,k,l);
end
end
end
end

%Determine SIaaa;
for i=1 :3
for j= 1 :3
aaatau(i,j)=0;
bbbtau(i,j)=0;
for k=1 :3
aaatau(i,j)=aaatau(i,j)+aatau(i,k)*a(k,j);
bbbtau(i,j)=bbbtau(i,j)+bbtau(i,k)*(a(k,j )-ID(k,j )/ 3);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3 v
SIaaatau(i,j,k,l)=ID(i,j)*aaatau(k,l)+ID(i,k)*aaatau(j,1)+ID(i,l)*aaatau(k,j);

SIaaatau(i,j,k,l)=aaatau(i,j)*lD(k,l)+aaatau(i,k)*ID(j,l)+aaatau(i,l)*ID(k,j)+SIaaatau(i,j,k,
1);
end
end
end
end

%Determine S[(a*a)(a*a)];
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
Saaaatau(i,j,k,l)=aatau(i,j)*aatau(k,1)+aatau(i,k)*aatau(i,l)+aatau(i,l)*aatau(k,j);
end
end
end
end

89

 

%Determine S[I(a*a*a*a)];
for i =1 :3
for j=1 :3
aaaatau(i,i)=0;
for k=1 :3
aaaatau(i,j)=aaaatau(i,j)+aaatau(i,k)*a(k,j);
end
end
end

for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
SIaaaatau(i,j,k,l)=ID(i,j)*aaaatau(k,l)+ID(i,k)*aaaatau(i,1)+ID(i,l)*aaaatau(k,j);

SIaaaatau(i,i,k,l)=aaaatau(i,i)*ID(k,l)+aaaatau(i,k)*ID(j,l)+aaaatau(i,l)*ID(k,i)+SIaaaatau
(H.191);
end
end
end
end

%Determine a:a (Ila) and alpha(2, l);
IIatau=aatau(l ,1)+aatau(2,2)+aatau(3,3);
IIbtau=bbtau( 1 , 1)+bbtau(2,2)+bbtau(3 ,3 );
alphatau(2, l )=2*IIatau/ 3 5;

%Determine alpha(3,1);

ALP 1tau=aaatau( l , l )+aaatau(2,2)+aaatau(3,3);
alphatau(3 , 1 )= l / 3 5+(4/ 3 5)*(ALP 1 tau/IIatau);
HIatau=aaatau( 1 , l )+aaatau(2,2)+aaatau(3 ,3 );
IIIbtau=bbbtau( 1 , l )+bbbtau(2,2)+bbbtau(3 ,3);

%Determine alpha(3,4);
alphatau(3,4)=-1/(7*IIatau);

%Determine alpha(3,5);
alphatau(3,5)=l/IIatau;

%Determine alpha(3,6);
alphatau(3,6)=-4/(7*IIatau);

%Determine alpha(4,l);

ALP2tau=aaaatau(l , l)+aaaatau(2,2)+aaaatau(3 ,3);
alphatau(4, l)=2/ 35 *ALP2tau+(1/3 5)*IIatau*IIatau;

90

%Determine alpha(4,4);
alphatau(4,4)=(- 1/7)* IIatau;

%MSU contribution;
%Determine pppp1,pppp2,pppp3,pppp4;
for i=1 :3
for j=1 :3
for k=1 :3
for l=1 :3
pppp5taU(i.i,k,l)=a(i.i)*a(k,l);
ppppltau(i,j,k,l)=alpha(1 , l)*SII(i,j,k,l)+alpha( l ,2)*SIatau(i,j,k,l);

%pppp2(i,j ,k,l)=alpha(2, 1)* SII(i,j,k,l)+alpha(2,2)* SIa(i,j ,k,l)+alpha(2,3)* Saa(i,j,k,l)+alp
ha(2,4)*SIaa(i,j,k,l);

pppp2tau(i,j ,k,l)=alphatau(2, 1)* SII(i,j ,k,l)+alpha(2,3)* Saatau(i,j ,k,l)+a1pha(2,4)* SIaatau(
i.j,k,1);

%Casel :Michigan State University model;
ppppOtaU(i.j.k.l)=C1*ppppltaU(i.i,k,l)+C2*pppt>2taU(i.i,k,1);
end
end
end
end

%Determine <pppp>zS
for i=1 :3
for j=1 :3
Mtau(i,j)=0;
for k=1 :3
for l=1 :3
Mtau(i,j)=Mtau(i,j)+pppp0tau(i,i,k,l)*s(l,k);
end
end
end
end

%Determine <pp>:<pppp>
for i=1 :3
for j= l :3
Ntau(i,j)=0;
for k=1 :3

91

for l=1 :3
Ntau(i,j)=Ntau(i,j)+pppp0tau(i,j,k,l)*a(l,k);
end
end
end
end

%Computation of tau
for i=1 :3
for j=1 :3

t2111(i.j)=(&|(i,.i)-1/ 3*ID(iJ))-U*(aa(iJ)-N(iJ))+Pe*M(iJ);

if i=2
if j=3
W(n)=taU(i.i);
time(n)=(n-1 )*dt;
end
end

X(n)=taU(iJ);
time(n)=(n-1 )*dt;
end
end
if i==1
if j==
Y(n)=taU(iJ);
time(n)=(n- 1 )* dt;
end
end
if i==3
if j==3
Z(n)=taU(iJ);
time(n)=(n- 1 )*dt;
end
end
end
end

%Primary normal stress
for i=1 :3
for j=1 :3
P(n)=Z(n)-X(n);
end

92

end

%Secondary normal stress
for i=1 :3
for j=1 :3
Q(n)=X(n)-Y(n);
end
end

%Time average

f(l)=0;

f(n+ l )=f(n)+1 ;

R(n)=0;

if all(f(n)>limitnl) & all(f(n)<limitn2)
R(n)=P(n);

end

end

averagetau=sum(R)/round(limitn2-limitn1 - 1 )
mean(P(n))

xIII=[-1/36:0.0001:2/9];

yII=2/9+2*xIII;

nyI=[0:0.0001:1/6];
xxIII=-6*(nyI/6)."(3/2);

ynyI=[0:0.0001:2/3];
xxxIII=6*(ynyI/6)."(3/2);

Smin=[-.04 :0.01: 0.05];
Smax=[0.1 0.15 .2];

IIbmin=[.08l6 .0816 .0816 .0816 .0816 .0816 .0816 .0816 .0816 .0816];

IIbmax=[.552067 .552067 .552067];

%plot(xIII,yII,xxHI,nyI,xxxIII,ynyI,III,II,':',Smin,IIbmin,'--',Smax,Hbmax,'--

'),xlabel('IIIb'),ylabe1('IIb');

%plot(xIII,yII,xxIII,nyI,xxxIII,ynyI,III,II,'+'),xlabel('IIIb'),ylabel('IIb'),axis([.2 .24 .6

68]);

%plot(xIII,yII,xx111,nyI,xxxIII,ynyI,III.II,'-.').xlabel('IIIb'),ylabel('IIb');
%curve D

%plot(xIII,yII,xxIII,nyI,xxxIII,ynyI,III,II,':',-

0.01 13,0.0915,'o',0,0,'square'),xlabel('IIIb'),ylabel('IIb'),axis([-.03 .001 -0.05 .2]);
%curve B

93

%plot(xIII,yII,xxHI,yyH,xxxIII,ynyI,III,H,':',-
0.0113,0.0915,'o',0,0,'square'),xlabel('IIIb'),ylabel('IIb'),axis([-.05 .25 0.15 .68]);
%curve F

%plot(xIII,yII,xxHI,nyI,xxxIII,ynyI,III,II,':',-

0.01 13,0.09 l 5,‘o',0,0,'square'),xlabel('IIIb'),ylabel('IIb'),axis([-.01 .25 0 .67]);
%plot(time,ev22)

%plot(time,d,'*'),xlabel('det of b'),ylabel('t')

%a(2,3)

% E

%v

%lambdal(3)

%e3

% Hb

% IIIb

%aij curve

%plot(time,A,time,e3,'-.',time,D,'--',time,en,':'),xlabel('t'),ylabel('aij')

% maxlam=1.5*e3-O.5;

% plot(time,A,time,maxlam,'-.',time,D,'--',time,en,':'),xlabel('Time'),ylabel('Alpha')

%eigan director curves
%plot(time,e3,'--',time,e2,':',time,en,'-.'),xlabel('t'),ylabel('vij'),axis([0 15 -1.2 1.2])
%plot(time,l3,time,e3,'--'),xlabel('t'),ylabel('vij')

%eigan value curves
%plot(time,12,'-.',time,l3,'--',time,l 1),xlabel('t'),ylabel('lambdal')
%plot(time,en,'+')

%stress tensor tau curves
%plot(time,W,time,X,'-.',time,Y,':',time,Z,'--'),xlabel('t'),ylabel('tauij')

%Primary and Secondary stress curves
%plot(time,P,time,Q,':'),xlabe1('t'),ylabel('normal stress');

94

LIST OF REFERENCES

Advani, 8.6. and CL. Tucker III, 1987, “ The Use of Tensors to Describe and Predict
Fiber Orientation in Short Fiber Composites ”, J. Rheol., 31 (8), 751-784.

Advani, S.G., and CL. Tucker III, 1990, “ Closure Approximations for Three-
Dimensional Structure Tensors ” , J. Rheol., 34 (3), 367-386 .

Aharoni, S.M., 1979, “ Rigid Backbone Polymers, 2: Polyisocyanates and Their Liquid —
Crystal Behavior ” Macromolecules, 12, 94-103.

Aharoni, S. M., and E. K. Walsh , 1979, “ Rigid Backbone Polymers. 4. Solution
Properties of Two Lyotropic Mesomorphic Poly(isocyanates)” Am. Chem. Soc, 12(2),
271-276.

Aharoni, S.M., 1980, “ Rigid Backbone Polymers, XVII: Solution Viscosity of
Polydisperse Systems, Polymer, 21, 1413-1422.

Bhave, A.V., R.K. Menon, R.C. Armstrong, and RA. Brown, 1993 “A Constitutive
Equation for Liquid-Crystalline Polymer Solutions”, J.Rheol., 37(3), 413-441.

Chaubal, C.V., L.G. Leal and G.H. Fredrickson, 1995, “A Comparison of Closure
Approximations for the Doi Theory of LCPs”, J.Rheol., 39(1), 73-104.

Cintra, J .S. and CL. Tucker III, 1995, “ Orthotropic Closure Approximations for Flow-
Induced Fiber Orientation ”, J. Rheol., 39 (6), 1095-1122.

Collings, CI. and J.S. Patel, 1997, Handbook of Liquid Crystal Research, Oxford
University Press, London.

Deen, W., 1998, Ananlysis of Transport Phenomena, Oxford University Press, New
York.

Doi, M, 1981, “Molecular Dynamics and Rheologoial Properties of Concentrated
Solutions of Rodlike Polymers in Isotropic and Liquid Crystalline Phases”, J.Poly. Sci
Polym.Phy.Ed., 19, 229-243.

Doi, M. and SF. Edwards, 1986, T heory of Polymer Dynamics, Oxford University Press,
London.

Feng, J., C. V. Chaubal, and LG. Leal, 1998, “ Closure Approximations for the Doi
Theory: Which to Use in Simulating Complex Flows of Liquid-Crystalline Polymers ”, J.
Rheol., 42(5), 1095-1119.

95

Folgar, F. and CL. Tucker III, 1984, “ Orientation Behavior of Fibers in Concentrated
Suspensions ”, J. Reinf Plast. Compos., 98-119.

Hand, G.L., 1962, “ A Theory of Anisotropic Fluids ”, J. Fluid Mech., 13, 33-46.

Hinch, E]. and LG. Leal, 1975, “ Constitutive Equations in Suspension Mechanics. Part
1. General Formulation ”, J. Fluid Mech., Vol. 71, Part 3, 481-495.

Hinch, E]. and LG. Leal, 1976, “ Constitutive Equations in Suspension Mechanics. Part
2. Approximate Forms for a Suspension of Rigid Particles Affected by Brownian
Rotations ”, J. Fluid Mech., Vol. 76, Part 1, 187-208.

Irnhoff, A., 2000, Mechanical Engineering, Diplomarbeit: Validation of Closure Models
for Fiber Induced Alignment of Fibers, Michigan State University/University of Aachen,
2000

Jeffery, GB, 1922, “ The Motion of Ellipsoidal Particles Immersed in Viscous Fluid ”,
Proc. R. Soc. Land. A, 102, 161-179.

Kini, H.K., Y.C.Kim, C.T. Nguyen, A. Bénard, and C. A. Petty, 2003, "Flow Induced
Microstructure in Composite Materials", Proceedings of the 14th International
Conference on Composite Materials, ICCM14 CD-ROM, San Diego, California

July 14-18, 2003

Larson, R.G., 1999, The Structure of Complex Fluids, Oxford University Press, London.

Nguyen, C. T. 2001, Molecular Orientation of Liquid Crystalline Polymers Under Simple
Shear Flows, MS Thesis, Michigan State University

Parks S. M., C. A. Petty, and M. Shafer, 1999, “ Flow-Induced Alignment of Fibers in the
Absence of Fiber-Fiber Interactions ”, paper presentation, Symposium on Suspensions,
APS/DFD, New Orleans, LA, November 21-23.

Petty C. A., S. M. Parks, and S. M. Shao , 1999, “ Flow Induced Alignment of Fibers ”,
in Proceedings of 12” International Conference on Composite Materials, ICCM12/T CA,
Paris, July 5-9.

Rey, AD. and M.M. Denn, 2002, "Dynamical Phenomena in Liquid-Crystalline
Materials", Annual Review of Fluid Mechanics, 34, 233-266.

Tucker, CL, 1988, “ Predicting Fiber Orientations in Short Fiber Composites ”, p.95 in
Proc. Manufacturing Int ’1, AMSE Publication, New York.

Tucker 111, CL. and SO. Advani, 1994, “ Processing of Short-Fiber Systems ”, pp. 147-

202, in Flow and Rheology in Polymer Composites Manufacturing, Volume 10,
Composite Materials Series, Editor: 8.6. Advani, Elsevier

96

Weiss R. A. and C. K. Ober, 1990, “ Current Topics in Liquid-Crystalline Polymers ” ,
pp 1-13, in Liquid Crystalline Polymers , ACS Symposium Series 435, American
Chemical Society, Washington, DC, Edited by : R. A.Weiss and C. K. Ober

97

 

 

IIIIIIIIIIIIIIIIIIIIIIIII

Ill ll ll» Hj‘llllll‘ l ll lllll‘ ill ‘I

     

2470 0456