:nvpr. .

\:

.._ an... M.
sam‘mwnm, ,

>
t r

. _
‘ “H?
x

x i:

.xL... r4 :6. .
«iiimwm.
1 Q.

.2 33m
4%.. .
haw nu

 

-
Mafig‘

. - 2:3. .
. ,1"
w Ed. 2.

 

 

 

‘lll - .Hn Fdaglfiufi mLLwJW»
rt , I n . l

This is to certify that the
dissertation entitled

Simulation of Flow-Induced Fiber Orientation with a New
Closure Model Using the Finite Element Method

presented by
Dilip Kumar Mandal

has been accepted towards fulfillment
of the requirements for the

doctoral degree in Mechanical Engineering

 

 

«f/ MZL/L / (ii/W //

Majori’Professor’s Eignature
/
j /G ,/i? f/

Date

‘\
\

MS U is an Affirmative Action/Equal Opportunity Institution

 

 

LIBRARY
Michigan State
University

 

 

 

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 cJCIRC/DatoDuepss-p. 15

 

 

 

4—‘_ ‘ -..—._-_ _ -. —.. . _t.

SIMULATION OF FLOW-INDUCED FIBER ORIENTATION WITH A NEw
CLOSURE MODEL USING THE FINITE ELEMENT METHOD
By

Dilip Kumar Mandal

A DISSERTATION
Submitted to
Michigan Sate University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2004

ABSTRACT

SHVIULATION OF FLOW-INDUCED FIBER ORIENTATION WITH A NEW
CLOSURE MODEL USING THE FIN IT E ELEMENT METHOD

By

DILIP KUMAR MANDAL

Flow induced alignment of particle suspensions is of practical importance for
various industries related to materials processing. The alignment of non-spherical
particles influences directly the properties of the products and is linked to the flow of the
material while manufacturing. For predicting flow induced alignment, it is possible to use
a distribution function. This approach however is computationally intensive since it
requires the solution of a complicated partial differential equation. Instead, an approach
based on the moments of an orientation distribution function is commonly used. This
approach requires a closure model since the equation governing the second moment of
the orientation distribution function involves an unclosed fourth order tensor. The further
development and validation of a closure model developed at Michigan State University
has been pursued in this work. The performance of the model, and some variations of it,
has been studied for various specified flow fields. Attempts have been made to identify a
constant coefficient that provides satisfactory solutions for these flows. A finite element
code has been developed to solve practical problems using this closure model along with
other popular closure models. This involved solving an equation for the evolution of the
second moment of the distribution function (the second order orientation tensor), coupled
with the momentum and continuity equations, as well as with a constitutive model for the

non-Newtonian stresses generated by the suspension. The stress developed due to

presence of fibers has been accounted in the momentum equation for improving the
accuracy of the predicted pressure and stresses in the computational domain. A variety of
two dimensional problems are solved, including flow in a sudden compression, flow
between rotating eccentric cylinder, and flow in diverging sections with various angles.
The results from the simulations with the closure model are compared with the well-
known quadratic and hybrid closure models. The model can predict quantities such as a
first normal stress difference, unlike the quadratic closure model. The computational cost
of these simulations is quite high, even only for two dimensional problems due to the
large number of unknowns. This led to the development of a new method which
combines finite element and analytical solutions for the reduction of computational time
and memory by using Greens functions. The analytical solution, based on Green's
functions, can be calculated for a simple domain and combined with the finite element
method to obtain a solution in the entire domain. This reduces effectively the
computational time for large problems. This method has been applied to two-dimensional

heat transfer, and fluid mechanics problems.

To my parents

ACKNOWLEDGMENTS

I am grateful to my dissertation advisor Dr. Andre Benard for his valuable
academic and non-academic, continuous guidance, and motivation. He helped me a
great deal in the formulations discussed in this work. The constructive suggestions
that he made during the reading of the manuscript of this dissertation helped in
improving it a great deal.

I would like to thank Dr. Charles Petty for the suggestions he made during the
course of my research work. I want to acknowledge Dr. Alejandro Diaz’s and Dr.
Zhengfang Zhou’s help for serving on my dissertation defense committee.

I would like to thank all my instructors at Michigan State University; Indian In-
stitute of Technology, Madras; and Jadavpur University, who taught and encouraged
me to reach this point. I would like to recognize the support and encouragement
provided by my friends. Finally, a. very special thanks to my parents and brothers for
their love and moral support.

I gratefully acknowledge partial financial support of this work by the National Sci-
ence Foundation through the following education and research programs at Michigan
State University: NSF/ECC/CRCD-9980325; NSF CTS-9986878; and an I/UCRC—

MTP Dissertation Completion Fellowship.

CONTENTS

LIST OF FIGURES XX
LIST OF TABLES xxi
1 INTRODUCTION 1

2 PREDICTION OF FLOW-INDUCED ORIENTATION OF SPHER-

OIDAL PARTICLES 4
2.1 Orientation of a single particle ................... 4

2.1.1 Orientation state of fibers .................. 8
2.2 The fiber distribution function ................... 8

2.2.1 Prediction of flow-induced fiber alignment using a dis-

tribution function ....................... 9

2.3 Average orientation state ...................... 13
2.3.1 Algebraic restrictions on average orientation tensors . 14

2.3.2 Prediction of flow induced fiber alignment ........ 18

2.3.3 Realizable orientation states ................. 20

3 THE DISTRIBUTION FUNCTION 24
3.1 Analytical solution ........................... 24
3.1.1 Elongational flow ........................ 25

3.1.2 Simple shear flow ....................... 27

vi

3.2 Numerical simulation of the distribution function ....... 30

3.2.1 Discretization of the distribution equation and the

mesh used ............................ 30

3.2.2 Discretized equations ..................... 32

4 THE CLOSURE MODEL 38
4.1 Existing closure models ........................ 38
4.1.1 The uncoupling (Quadratic) approximation ....... 39
4.1.2 Linear model of Hand ..................... 39
4.1.3 Hinch and Leal model ..................... 40
4.1.4 Hybrid model of Advani and Tucker ........... 40
4.1.5 Orthotropic closure model .................. 41

4.2 Closure model developed at MSU ................. 44
4.2.1 Construction of the fourth order orientation dyadic . . 44

4.2.2 Reducible form of the fully symmetric closure model . 47

4.2.3 Irreducible form of the fully symmetric closure model 50

4.3 Estimating the coefficients of the parameters .......... 51
4.4 Fitted model parameters ....................... 56
4.4.1 Summary ............................ 69

5 PREDICTION OF FLOW-INDUCED FIBER ORIENTATION US-
ING THE FINITE ELEMENT APPROACH 73
5.1 Finite element formulation ...................... 74

5.1.1 Governing equations ...................... 74

vii

5.1.2 Weak form of the equations ................. 75
5.1.3 Vector form of the discretized equations ......... 76
5.2 Testing the closure model implementation in simple shear flow 79
5.3 Flows in complex geometries .................... 82

5.3.1 Approximations introduced in two dimensional calcu-

lations .............................. 83

5.3.2 Pressure driven flow between two parallel plates . . . . 84
5.3.3 Converging channel ...................... 87
5.3.4 Diverging channels ....................... 89
5.3.5 Flow between eccentric cylinders .............. 95

5.4 Summary ................................. 99
PARAMETRIC STUDY 103
6.1 Effects of varying the fiber aspect ratio .............. 103
6.2 Effects of varying the aspect ratio with diffusion ........ 106
6.3 Effects of varying the concentration of the fibers ........ 106
6.4 Effects of varying inlet velocity ................... 106
6.5 Summary ................................. 111

USE OF GREENS FUNCTION TO INCREASE COMPUTATIONAL

EFFICIENCY 115
7.1 Theory .................................. 116
7.1.1 Viscous flow problem ..................... 118

7.1.2 Heat conduction problem .................. 122

viii

7.1.3 Remarks of the proposed method ............. 123

7.2 Computational results ........................ 123

7.2.1 Heat transfer .......................... 124

7.2.2 Creeping flow .......................... 127

7.3 Summary ................................. 131

8 SUMMARY AND CONCLUSIONS 133
APPENDICES

A. Tensor notation used in the dissertation ............. 137

Double-dot product .......................... 137

Dyadic product ............................. 137

B. Numerical model for distribution function ............ 138

C. Parametric approximation by shape functions .......... 139

D. Stress in a suspension of force-free particles ........... 142

Stress prediction by the proposed closure ............ 145

E. Upwinding ................................ 148

F. Matrix entities in the finite element formulation ........ 151

ix

2.1

2.2

2.3

2.4

2.5

2.6

3.1

3.2

3.3

LIST OF FIGURES

Fiber orientation vector .........................

The fiber rotational vector (p) is perpendicular to the orientation vec-

Fiber orientation balance .........................
Tensor representation of average orientation state (similar to [21]) . .

Diagram of the possible states of the orientation tensor in terms of the

eigenvalues (a1 + (L2 + (I3 = 1) ......................

The shadowed region shows the realizable orientation state as a func-

1

tion of the invariants of the anisotropic orientation tensor (b = a — —6).

3

Three dimensional plot for distribution function for uniaxial flow. The
distribution function evolves with time when fibers are randomly ori-

ented initially (i.e. t/J(0, (z), 0) = 1/47r) ...................

Components of average orientation tensor obtained from analytical for-
mulation of the distribution function in a uniaxial flow. It is assumed

that. all fibers are randomly oriented initially ...............

Three dimensional plot for distribution function for simple shear flow.
The distribution function evolves with time when fibers are randomly

oriented initially (i.e. 1MB, (1), 0) = 1/47r). ................

X

10

18

21

23

26

27

28

3.4

3.5

3.6

3.7

3.8

4.1

4.2

Components of average orientation tensor Obtained from analytical for-
mulation of distribution function for simple shear flow. It is assumed

that all fibers are randomly oriented initially ...............

Mesh used for solving the fiber distribution function over the spherical

surface. ..................................

Numerical error Of the normalization condition for simple shear flow
with C1 = 0.0, A = 1 while solving Smoluchowski equation (equation

2.17) .....................................

Numerical error Of the normalization condition for simple shear flow
with C] = 0.01, A = 1 while solving the Smoluchowski equation (equa-

tion 2.17). .................................

The eigenvalues of average orientation tensor for simple shear flow with
the diffusion coefficient C1 = 0.01. All the fibers are assumed randomly

oriented initially. .............................

Graphical representation of the values of Cis of the reducible fully
symmetry model that match the numerical results and values of the
invariants with dimensionless time (17). These results are for simple

shear flow without fiber interactions and shape factor A = 1. .....

Graphical representation of the values of CS that match the numerical
results for the irreducible fully symmetry model and t is the dimen-
sionless time (77). These results are for simple Shear flow without fiber
interactions and Shape factor A = 1 ....................

xi

29

35

36

37

54

55

4.3

4.4

4.5

4.6

4.7

Sensitivity of invariant. II to the constants Cis (52,- = 811/ 0G,) of the
reducible fully symmetry model. The conditions are the same as Figure
4.1. t is the dimensionless time (7‘7) and these results are for simple

shear flow without fiber interactions and shape factor A = 1. .....

Sensitivity of invariant III to the constants GS (53,- : 8111/80.) of
the reducible fully symmetry model. The conditions are the same as
Figure 4.1. t is the dimensionless time (17) and these results are for

simple shear flow without fiber interactions and shape factor A = 1.

Sensitivity of invariant. II to the constants Cis (S2,: = (911/863,) of the
irreducible fully symmetry model. The conditions are the same as
Figure 4.1. t is the dimensionless time (c'fr) and these results are for

simple shear flow without fiber interactions and shape factor A -— 1.

Sensitivity of invariant III to the constants CS (33,- : GUI/ad) of
the irreducible fully symmetry model. The conditions are the same as
Figure 4.1. t is the dimensionless time (1'7) and these results are for

simple shear flow without. fiber interactions and shape factor A = 1.

Comparison between the results obtained from fully symmetric closure
and the results obtained from the exact solution for simple shear flow
without fiber interactions. The solid line shows the results obtained
by the closure model with parametric coefficients in equations (4.45 -

4.48) while the dotted points are the exact values ............

xii

57

58

59

60

62

4.8

4.9

4.10

4.11

4.12

The F SQ model parameter (Cg) has been back computed at the six

nodes of the realizable diagram by isoparametric quadratic shape func-

tions. The corresponding (C2) values are {—4.76738, 0.269573, 0.0693274,

0.275868, 0.398572, 0.830628} respectively. ...............

This result is for simple shear flow field. This result is for simple shear
flow field. Comparison between the results obtained from the numer-
ical solution and the closure model with 02 computed using shape
functions. The dotted points are the plotted from exact result whereas

the solid line from the closure model using shape functions .......

This result is for uniaxial flow field. By using the shape functions
with quadratic interpolation of C2 values in the realizable diagram,
the results also does not match to the expectation. The dotted points
are the plotted from exact result whereas the solid line from the closure

model using shape functions ........................

The calculations are plotted in the invariant diagram for simple shear
flow without fiber—fiber interaction. The dotted points are for the exact
solution. The dark solid line is for FSQ closure with Cg = 1/ 3 whereas

the faded solid line is for FSQ closure with CQ 2 1 / 2. .........

Comparison of the maximum eigenvalues with the solution obtained
from distribution function and FSQ(C2 = 0.37) for simple shear flow
without fiber interactions. X-axis represents the dimensionless time

w.r.to rate of strain. ...........................

64

65

66

68

4.13 Comparison of F SQ closure with two different parametric coefficients

(C2 = 1/3 and C2 = 0.37) and the exact solution for uniaxial flow. . . 71

4.14 Comparison with the solution obtained from distribution function and
F SQ(C2 = 0.37) for simple shear flow with fiber-fiber interaction(C1 :—

0.01) .................................... 72

5.1 Plot. of the maximum eigenvalues of the orientation tensor vs dimen-
sionless time as obtained from 3 sources: the analytical solution of the
distribution function, the. solution of the evolution equation with FSQ
closure, and the finite element solution of the evolution equation with

FSQ closure coupled with Navier—Stokes equation ............ 81

5.2 Contour plots of the principal eigenvalues(Amax) of the orientation vec-
tor superposed with main eigenvectors of the orientation tensor. The
results are shown for four different times and boundary conditions are
with velocity field (Va: = —y2—2 + 0.43; and Vy = 0.0 at inlet, and
Vy = 0.0 and a—SI—I = 0.0 at outlet) for Re 2 36 with random orienta—
tion state initially. ............................ 85

xiv

5.3

5.4

5.5

5.6

Contour plots of main (maximum) eigenvalues of fiber stress tensor
(01) are shown at different. times. The eigenvectors corresponding to
the main eigenvalues are also plotted with a weighting factor of the
respective main eigenvalues. The results are shown for four different
times and boundary conditions with velocity field (VI = ——y.; + 0.43)
and Vy = 0.0 at inlet, and Vy = 0.0 and 93:5 = 0.0 at outlet) for Re
2 36 with random orientation state initially. ..............
Contour plots of the pressure field are shown for four different times.
The boundary conditions with velocity field (V3: = —=";— + 0.43) and
Vy = 0.0 at. inlet, and Vy = 0.0 and 07% = 0.0 at outlet) for Re 2 36
with random orientation state initially. .................
The top figure shows a contour plot of Cos(x) superposed with main
eigenvectors of the orientation tensor. The contour plot of the main
eigenvalues ((71) of the fiber stress tensor as well as the corresponding
eigenvectors with weighting factor of the main eigenvalues are plotted
in the middle figure. The pressure field is plotted in the bottom figure.
The specified velocity field is (VI = —2(y2 — 0.0842) and Vy = 0.0 at
inlet, and Vy = 0.0 and a—d‘f = 0.0 at outlet) ...............
Contour plots of Cos()() superposed with main eigenvectors of the ori-
entation tensor with velocity field (VI = —2(y2 — 0.152) and Vy = 0.0
at inlet, and Vy = 0.0 and 937‘ = 0.0 at outlet) for different diverging
angles. Number of nodes and elements diverging cones 200, 300, 400

are 6813, 7075, 7183 and 1642, 1707, 1734 respectively. ........

XV

86

88

90

92

5.7

5.8

5.9

5.10

The contour plots of two eigenvalues of the fiber stress tensors (third
component is negligible) are shown. The eigenvectors corresponding
to main eigenvalues are also shown with weighting factor of the main
eigenvalues. The streamlines and pressure field Of the flow field are also
shown subsequently. The velocity field is (VI = ——2(y2 — 0.152) and
Vy :2 0.0 at inlet. and Vy = 0.0 and (13:2; 2 0.0 at. outlet) for different

diverging angle of 200. ..........................

The contour plots Of the main eigenvalue of the stress 01 as well as the
corresponding eigenvectors with weighting factor of the main eigen-
value. The corresponding eigenvectors are similar to the eigenvectors
of the orientation state at that local position. (V17 2 —2(y2 — 0.152)
and Vy = 0.0 at inlet, and Vy = 0.0 and 2i; 2 0.0 at outlet) for

different. diverging angles. ........................

Contour plots of eigenvalue(A) superposed with main eigenvectors of

the orientation tensor. Also the main eigenvectors of the fiber stress

tensor. The streamlines with the pressure field of the flow field are also
2

shown below. Re 2 (’17? = 0.336 in eccentric fiow (outer cylinder is

rotating and the inner one is fixed) ....................

Comparison of the FSQ model with the Advani Tucker hybrid model for
realizable orientation state (i.e. all the eigenvalues are positives) with
Re = 11:3 = 0.336 in eccentric flow field (outer cylinder is rotating

and the inner one is fixed) .........................

xvi

93

94

96

97

5.11

5.12

5.13

6.1

6.2

These test cases for courser mesh with Comparison of the FSQ model
with 2108 nodes and 496 elements. Advani Tucker hybrid model for
realizable orientation state (i.e. all the eigenvalues are positives) with
Re 2 3:3 = 0.336 in eccentric flow field (outer cylinder is rotating

and the inner one is fixed) .........................

The main eigenvalues of the stress tensor for three different models
have been plotted. The FSQ and quadratic models provide results,

which are close to each other unlike the hybrid model ..........

The pressure contours for three different. models have been plotted.

The profiles are almost similar. .....................

The effects of varying the fibers aspect ratio while there is no fiber—fiber
interaction. The eigenvalues of fiber stress tensor are clearly high for
longer fibers. The velocity field is (V13 = —2(y2 — 0.152) and Vy = 0.0
at inlet, and Vy = 0.0 and %£_1 = 0.0 at outlet) for different diverging

angles. ...................................

The effects of varying the fiber aspect ratio on Cos(x) while there is
no fiber-fiber interaction. The fiber orientation state does not. change
much with the aspect ratio of the fiber. The velocity field is (Vzr =
—2(y2 -— 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and 8—3;? = 0.0 at
outlet) for different diverging angles. ..................

xvii

98

100

101

104

 

 

6.3

6.4

6.5

6.6

6.7

The effects of varying the fiber aspect ratio with diffusion on Cos()().
The eigenvalues of the average orientation state is with lower magni—
tude than the corresponding eigenvalues without diffusion. The veloc-
ity field is (Vgr = —2(},/2 — 0.152) and Vy = 0.0 at. inlet, and Vy = 0.0
and 59%;“ = 0.0 at outlet) for different diverging angles ..........
The effects of varying the fiber aspect ratio with diffusion on the prin—
cipal fiber stress. The fiber stress is clearly higher than those without
diffusion. The velocity field is (VI = —2(y2 — 0.152) and Vy = 0.0 at
inlet, and Vy = 0.0 and d—gr—f = 0.0 at outlet) for different diverging
angles ....................................
The effects of varying the concentration of the fibers in the suspension
on C08“). The average orientation state does not defer much with the
concentration. The velocity field is (V11: 2 —2(y2 —0.152) and Vy = 0.0
at inlet, and Vy = 0.0 and % = 0.0 at. outlet) for different diverging
angles ....................................
The effects of varying the concentration of the fibers in the suspension
on the principal fiber stress. The fiber stress does change with the
concentration of the fiber suspension. The velocity field is (VT 2
—2(y2 — 0.152) and Vy = 0.0 at. inlet, and Vy = 0.0 and Q6? 2 0.0 at
outlet) for different diverging angles. ..................
The effects of varying the inlet velocity profile on Cos(x). The average
orientation state with different inlet velocity profile. The norm inlet.
velocity profile is (VI- : —2(y2 —- 0.152) and Vy = 0.0) ........

xviii

107

108

109

110

112

 

 

6.8

7.1

7.2

7.3

7.4

7.5

7.6

7.7

The effects of varying the inlet velocity profile on the principal fiber
stress. The norm inlet velocity profile is (V17 2 —2(y2 — 0.152) and

Vy = 0.0) .................................

Outline of the mesh required for the proposed combined finite element

and analytical method for a complex physical domain ..........

The boundary of the analytical subdomain ((21) has been divided into

N elements. ................................

Mesh and temperature distribution are Shown for full finite element
and combined method when top and bottom sides (y = 3:24) are

insulated, at left side (it? = —2.4) T = 50.0, and at right side (r = 2.4)

Same physical problem as in Figure 7.3 but an eight-node quadratic

element has been considered instead of four-node linear element. . . .

Temperature distribution when all the boundaries with temperature
T = 250 and heat generation, g = 10.0 for full finite element method

(left side) and combined method (right side) ...............

Grid of a complex geometry for full finite element method (left side)

and combined method (right side). ...................

Temperature distribution without heat generation(top two figures) and

113

116

119

126

127

with heat generation( bottom two figures) by full finite element method(left

side) and combined method(right side) .................

xix

129

 

 

7.8

7.9

8.1

8.2

8.3

8.4

Solution for Poiseuille flow by both the methods. Vy = 0.0 at inlet and
outlet; Pressure P = 24.0 at inlet and P = ~24.0 at outlet .......
Solution for fluid flow when Vy = 0.0 at inlet and outlet; pressure
P = 24.0 at inlet and P = —24.0 at outlet, at top and bottom portion

there is a suction of 5.0 unit ........................

Isoparametric quadratic shape function approximation has been taken
for the model parameter in FSQ C2 ....................
The normalized first normal stress difference N1 / [$77612 / (ln(2a.p) — 1.8)]
versus shear rate of the glass fibers of aspect ratios (1,, = 276 in New-
tonian solvents, glucose wheat syrup, and epoxy resin. The volume
fraction 40 is in the range for 0.0002 — 0.0011. FSQ closure model
(C2 = 0.37) provides very good results whereas the quadratic model
falls far short. ................ - ...............
The normalized first. normal stress difference NI / [(157702 (ln(2a.,,) — 1.8)]
versus shear rate of the glass fibers of aspect ratios up = 552 in New-
tonian solvents, glucose wheat syrup and epoxy resin. The volume
fraction ob is in the range for 0.0002 — 0.0011. FSQ closure model
(C2 = 0.37) provides very good results whereas the quadratic model
falls far short ................................
The velocity vector(V) is in XY plane. V is aligned with the axis -s

for streamwise upwinding. ........................

130

130

141

146

147

2.1

3.1

4.1

4.2

4.3

5.1

6.1

7.1

LIST OF TABLES

Boundaries of the realizable orientation states .............
Boundary conditions for mesh ......................

Contracted notation for Orthotropic closure model ...........
Coefficients of the four fully symmetric operators (reducible model)
Coefficients for four fully symmetric Operators with irreducible model

and B : 15 — 50 det(a) — 30(a : a)(1 — (a: a)) .............

The summary of the case studies are presented. A and 01 indicate the
maximum eigenvalue of the average orientation tensor and fiber stress
tensor respectively. The ("tomputational time (1000 see as an unit) is

also shown for each case on the right most column. ..........

Parametric study of the cone 200 diverging angle. A and 01 indicate the
maximum eigenvalue of the average orientation tensor and fiber stress
tensor respectively. The computational time (1000 see as an unit) is

also shown for each case on the right most column. ..........

Computational information for both the finite element and combined
method (NN : number of nodes, NE : number of elements, RM 2
required computational memory, RT : required time for solving, and

ME : maximum error according to analytical solution) .........

xxi

22

34

42

49

52

114

132

CHAPTER 1

INTRODUCTION

Polymers are Often reinforced with fibers to retain the benefits of each constituent
in the final product. Discontinuous fibers are introduced to enhance the strength
and stiffness of the products. Their orientation inside the product determines the
mechanical properties Of the material. Common manufacturing methods for polymer
composites with discontinuous fibers used in large volume production include injection
molding, extrusion and compression molding. In such processes, the alignments of
the fibers are influenced by the flow field, the materials used, the geometry of the

part, and the processing parameters.

During manufacturing of composites, various flow fields are encountered at dif-
ferent locations and time. The forces acting on the fibers cause rotation of the fibers.
AS a result, the fibers will tend to align with the flow direction. One approach to
Characterize and predict the orientation of particles in a flow field is based on using
a. fiber distribution function which provides the probability of orientating particles
in Certain direction. However, even with today’s available computational power, it is
Still very difficult to solve the governing equation for the fiber distribution function
in a. complex flow. A method based on the averaging the moments of the distribution
fuIlction, which is less computationally demanding, is preferred to represent particle

Orientation inside the flow fields. The averaging procedure however gives rise to a

closure problem involving a fourth order average orientation tensor (i.e. there are
more unknown parameters than the number of governing equations). Therefore a
relationship is needed for this fourth order average orientation tensor.
Various closure models, which are based on relating the fourth order orientation
tensor with second order orientation tensor, have been developed previously by a
number of researchers. Few of these however satisfy required algebraic properties of
the fourth order tensor. A new class of closure models has been developed recently
in Michigan State University. This class of models retains all the six symmetry and
projection properties of the fourth order orientation tensor. In this work, simulations
are performed to study the performance of this closure model. Various prescribed
flow fields have been used to check the realizability of the closure model so that the
Inathematical solution of the orientation state is realizable. The solutions with the
selected closure model are compared with the exact solutions for different prescribed
flow fields. These results are also compared with experimental results and results
Obtained by using the quadratic and hybrid closure models.

In order to solve practical problems, one must solve the momentum equations,
the continuity equation, the governing equations for the orientation state, and the
equations of a constitutive model for the fiber stresses. The equations are solved with
the help of the finite element method. The numerical implementation is discussed in
this work. Fiber orientation has been simulated for various geometries and processing
Parameters. The predictions of fiber orientation from the finite element calculations

fOF various geometries and flow fields have been discussed. The performance of this
Closure model is also compared with the pOpular quadratic and hybrid closure models.

2

Thereafter, the influences on the flow parameters; like aspect ratio of the fibers, fiber
concentration in the suspension, diffusion effect, and Reynolds numbers in the How
have been studied.

A new method which combines Green’s functions and the finite element method
is proposed for the efficient solutions of partial differential equations (not related
to flow induced fiber orientation). This approach requires identifying subdomains
where the solution to the governing equations can easily be calculated using Green’s
functions. Combining the analytical solution based on Green’s function and the finite
element calculation, outside of the analytical domain, a solution of the problem is
achieved without requiring a discretization of the analytical subdomain. The proposed
method is compared with the Galerkin finite element method by using examples
Of heat. transfer and fluid mechanics. The computational efficiency is discussed in
the examples. The possible sources of numerical errors, developed by the proposed
method, is discussed thereafter. This method can also be applied to unsteady, non-

linear elliptic and parabolic problems with multiple unknowns.

CHAPTER 2

PREDICTION OF FLOW-INDUCED

ORIENTATION OF SPHEROIDAL PARTICLES

Two most common methods are used to describe the statistics of fiber orientation
inside short fiber composites. The first. approach is based on using a fiber distribution
function, which provides the probability of orienting particle in certain direction. A

second approach, used extensively in this work, is based on using the moments of the

orientation distribution function.

2.1 Orientation of a single particle

The orientation of a single particle can be expressed with a. unit vector p aligned
Wi th the fiber. The unit vector can be characterized with the Euler angles 0 and gt
1 1 , 2] as shown in Figure 2.1. Only two variables (6, qt) are sufficient to describe the

Vector in spherical coordinate system. The vector can be represented in cartesian

coordinates by [3]

\

/ sin 6 cos (1')

P = sindsinqfi (2-1)

( c086 )

 

 

Due to symmetry, the head and tail of a fiber are indistinguishable and the choice

4

 

 

 

Figure 2.1: Fiber orientation vector

of the direction can be both ways. This implies that p —+ —p or alternatively 6 —+ 7r—6
and d) —+ 7r + 95. This property is valid only for rigid straight particles. A different
description can be used for flexible or curved particles[4].

The geometry of a spheroidal particle is Often described with an hydrodynamic

interaction coefficienti[5]

 

(2.2)

where (1,, is the aspect ratio (length / diameter). It is worthwhile to mention that A = 1
represents for long fibers with 00 aspect. ratio, A = 0 indicates for sphere, whereas
A = —1 for disks. In this dissertation only rigid particles have been considered. The
density of the fibers is assumed to be equal to the density of the fluid (buoyancy is
negligible). The suspensions are also assumed to be incompressible.

The motion of a particle p can be characterized by the velocity of the centroid as
in Figure 2.2 and the rotation vector will be p[6]. From Figure 2.2 it is clear that the
velocity vector p is perpendicular to the vector p since there can not be a relative
Velocity inside the rigid particle along the vector p.

The motion caused only by fluid movement, in was derived by Jeffery[7] in 1922.
He derived an equation of motion of a single rigid ellipsoidal particle suspended in
Newtonian fluid in the absence of fiber interactions[8] and no slip boundary condition

On the surface of the particle as follows

In = -w - p+A(S - p — S 2 mm) (2.3)

Where 117 and S are the vorticity tensors and rate of strain respectively.

6

Fiber

P P

X“) Fluid motion

 

 

X1

Figure 2.2: The fiber rotational vector (p) is perpendicular to the orientation vector

p.

2.1.1 Orientation state of fibers

In practical applications, a. very large number of fibers are found inside a product.
Even with today’s computational power, it is still impractical to follow every particle
individually. An average statistical approach is favored to describe the motion of a
group Of particles. A small volume element sufficiently large to contain many fibers,
but small enough to have uniform orientation statistics, is therefore considered.

There are two possible approaches to describe the alignment of many fibers inside
a small volume element. The first one is based on a fiber distribution function IMO, ct),
that describes the fraction of particles in a given orientation direction. The second
method is based on spatially uniform the average orientation state. Both approaches

are discussed below, beginning with the distribution function.

2.2 The fiber distribution function

The fiber distribution function 2,0(6, gt) or ti'(p) represents the fraction of particles
in a given orientation state inside a small domain. For example, if ten percent of
the fibers are orientated in one direction, the distribution function in that direction,
L (P) is one tenth. Since every particle must be oriented in a particular direction, the

Spatial integral of the distribution function equals to unity. The distribution function

Inust have the following relation[9]

27r
/0 fofiii'w )sin dddde'b: fed) p)dp— — 1 (2.4)

The above relation (equation 2.4) is called the normalization property of the orienta-

tion distribution function[2]. As the head or the tail of a fiber can not be distinguished,

8

the following symmetry condition is also obtained as

We) = Iii-p) 01‘ 123(9. <25) = MW — 9, 7? + ¢> (2-5)

The distribution function contains the complete information about the alignment of
fibers because the function provides the fraction of fibers inside the small volume for

each possible orientation state.

2.2.1 Prediction of flow-induced fiber alignment using a dis-

tribution function

The forces exerted on a particle as it moves depend on the fluid properties, the flow
field, as well as the properties of the suspended particles. The forces are mainly di-
vided into three parts: the hydrodynamic force, the interaction forces among the par-
ticles (Brownian motion, affinity among the particles), and external effects (i.e. mag-
netic forces). In this work, it is assumed that. there are no external non—hydrodynamic

effects.

In the presence of hydrodynamic forces, the deviation of the particle vector from
the velocity vector results in an unbalanced force (Figure 2.2). This implies that
the forces exerted on the fiber create a moment, which may cause rotation. Thus
fibers will be influenced by the flow field, and it is necessary to predict the expected

alignment [4] .

For a small control volume, the velocity gradient tensor can be used to describe

9

 
  
     

unit sphere

  

control
volume

Figure 2.3: Fiber orientation balance

a. flow field as[2, 10]

“1.1 111,2 “1,3

VU : ] 11.2.1 U22 U23 I (2.6)
\ U3,1 “3.2 713.3 )
“filere the rate of strain and vorticity tensors are defined respectively as

Vu + VuT

S = —— (2.7)
2
Eind

._ i T

w : \7u 2VIL (2.8)

Where S is symmetric and w is antisymmetric tensor.

10

Assuming a small control volume inside the domain (Figure 2.3), the rotational
BuJ-c , 15¢! at the boundaries must. be in balance with the change of orientation with

respect to time inside the control volume (unit. sphere). Flux through boundaries can

Flux = — {f (be) - wads = — /V// v. - (25¢) W (2.9)

\vhere the divergence theorem is used[1l], and the surface gradient is defined as

be vvritten as

VSZEZeBB 1 8

for 1:11; case of a unit sphere. The change of orientation with respect to time inside

the control volume can be written as % fff i,’)(p)dV. This results in fiber orientation
v

— ff (15w) ads = g,- /// new (211)
A V

allcl Using the divergence theorem, the following equation can be obtained[12]

bal a 11cc given by

De”) (9 .
D, — 71; - (In!) (2.12)

 

Egu ation (2.12) is known as the Fokker-Plank equation[4]. D'l/J/Dt is the material

deriVative. O/Bp is the special derivative in the spherical coordinate system. p, the

ISgt-Ethion vector of the particle p, depends mainly on the convective motion of the
particle caused by fluid movement and the interaction among the fibers.

The motion caused by hydrodynamic movement, 13.; is derived by Jeffery in
1922M. The governing equation as discussed before, is represented by the follow-
ing equation as

R: = -w - p+A(S ' p — S 2 mm) (2.13)

11

 

To account for the interaction forces among the fibers, there are few models avail-
able. The rotational contribution due to interactions among the fibers is represented
by (p — 133). To model fiber rotation caused by the fiber interactions, the following

hypothesis is used in this dissertation

81/;
. —— . 'il,’ : —D ' 2014
(P PJ)lr (P) R—ap ( )

xvllere DR is a diffusion coefficient that. accounts for the interactions between fibers.

Hence, p can be written as[13, l4]
. 31D
p=—w-p+A(S-p—S:ppp)—DR—a; (2.15)

In t he literature, different. approaches for modelling the rotary diffusion coefficient
(D R) are found. For this work, only one model has been used.
Arguing that the rotary diffusion coefficient should be dependent on the local

Strain rate, Folgar and Tucker [15, 16] developed the following expression for DR
D[{:C](Q5)VQSZS (2.16)

\Vllere C1(¢) is the dimensionless particle-particle interaction coefficient. This coeffi—
Qieht depends on the volume fraction (Z) of the suspension and is determined empiri—
Q 8L11y. The fiber volume fraction d), can be defined as 9'9, = %L = fl where v”, =
f1her, fluid and suspension volume respectively. The Folgar-Tticker relation implies
that for Vu = 0, the rotation of fibers caused by fiber interactions becomes zero
instantaneously. The foregoing model does not relax to an isotropic orientation state

after removing the hydrodynamic force like the Brownian diffusion does.

12

Substituting p of equation (2.15) into equation (2.12) and expanding in the spher-

ical coordinate system, the following result can be obtained[13, 17]

 

 

 

 

 

£3? = 3%; + gig—235% + DRcot6%%——
%{(A;1VHTIQTQO) + (A:1nggre6)}+
sifi6%g{A2—1VQT:§§¢ + A31Vg2ere¢}+ 2,9(3AV1LI : grief) (2.17)
vvl'lere g” = {sinBcos 9’9, sianin (1), cos q)},g0 = {cosdcos ¢,cosdsin¢,—sin6},gf’ =
f ‘ sin 95, cos (1‘), 0}, and DR is the rotary diffusion coefficient.

The solution of equation (2.17) provides the evolution of the distribution function
vvi t1’1 time. Hence, it represents the exact and full prediction for fiber orientation inside
the suspension. It is not possible to solve analytically the equation (2.17) for general
Cases. Numerical methods have been used to solve the governing equation, but even

Wi t 11 faster computer it takes days for each particle (see section 3.2).

2 - 3 Average orientation state

A second method to describe fiber alignment is based on the average orientation

tQIISQrB]. An average of N numbers of B,- : i = LN 13 defined as

N
I
B =-— B,- 2.18
< > N; ( )

The average < p > for the orientation vector and all odd number of products of
the orientation vector (i.e. < ppp >) equals zero because of the symmetry implied
by p —> —p. Thus, average orientation can be defined by only even numbers of
vector outer products[11]. For example, a,,- =< p,p,- > is a second order tensor and

13

Ami-.1 =< p,p,pkp, > is a fourth order tensor. (1,,- and Aijkl are defined in terms of

fiber distribution function as[18]

271' 7r
(1,]- =< p,[)j >= f1),p,~i,(i(p)dp = /0 /0 p,pjzii(6,¢)sin ddfidqé (2.19)

and

27r 7r
Aijkl =< ptpjpt-pz >= fprryptpuiipldp = f / prpjpkpz¢(9,¢)sin9d9d¢ (2-20)
0 0

T1] ere actually exists an infinite number of higher-order orientation tensors. The dis-
tri I) ution function can also be constructed from higher order tensors with the following

series[3]

1 315

, 15
1MP) = — + gbzafz‘jlpl + 32—Whjtifrjt—dp) + ° " (221)

V ' - __ 1 ‘ , _ l
Vllel Q bij — at} — 360 and bijkl — aijkl — 7((L1j6k1‘f‘ aik5j1+ (1,153), + (119160 + (1316,}, ‘l'

aj I: (52- 1) + 3%(6U6H + (5,),de + (5,1610. Using additional higher order tensors will increase

the accuracy of the approximation of the distribution function.

2-3 - 1 Algebraic restrictions on average orientation tensors
It can also be noticed that. a change in the indices in < pm, > and < pipjpkp, >
(ecluation 2.19 and 2.20) does not change the average orientation states. Hence, from
the equations (2.19) and (2.20), it is also clear that the second order tensor, a and

EOurth order tensor, A are symmetric. Thus
aw- : a], (2.22)

and
Aijkl = Ajikl = Aijlk = Aklij = Aikjl (2.23)

14

Furthermore, due to normalization condition (equation 2.4), the trace of the second
order tensor is unity.

tr(a) = (1112‘ = 1 (2.24)

where tr(a) : the trace of the square orientation matrix a. Since the dot product of
tvvo orientation vectors p - p equals one, the projection between any two orientation
vectors in the fourth order orientation tensor < pppp > must reduce to the second
order orientation tensor i.e. < pp >. Therefore, the fourth order orientation tensor
should satisfy the following projection property[19] (please see appendix A. for tensor

2111 a] ysis notation)

3 3 3 3 3 3
(lij = Z Assij = Z Asisj = Z Asijs = Z: Aissj = Z Aisjs = ZAijss (225)

3:1 3:1 5:1 3:1 3:1 5:1

and for higher order tensors
3
A,,-,,l : Z Afikmm (2.26)
m=l

Using the symmetry and the normalization conditions (equations 2.22-2.24), five
independent entities of the second order tensor (an, (112, (113, (122, (123) are obtained.
The form of the fourth order tensor can also be studied by applying the symmetry
and projection properties in order to determine the number of independent entries.
Tlle tensor Aijkl has 81 components and using Aijkl = Ajz’kl 2 AU”; reduces the
Possibilities to 36 distinct coefficient at most. Again by using Aijkl = Aklij, the
I1umber of distinct components in Aijk, reduces further from 36 to 21. The symmetry
Condition by interchanging second and third indices (i.e. Aijkl = Aikfl) given in
equation (2.23) leaves only 15 independent entries in the fourth order orientation
tensor[20]. The 15 independent entries are, for example, the components 1111, 1112,

15

1113, 1122, 1123, 1133, 1222, 1223, 1233, 1333, 2222, 2223, 2233, 2333, 3333 of AW.
Since after projecting the Aijkl into the second order tensor aij, six equations can be
obtained the number of unknown in the fourth order orientation tensor will be nine.
However, the equations (2.25) and (2.26) imply that higher order tensors contain

complete information about the lower order tensors.

Since the orientation tensors contain information only about the average of the
fiber distribution inside the small domain, the knowledge of complete fiber distrib—
ut ion function is required to present the exact alignment of fibers. But the main
dis acivantage of using the fiber distribution for computing the orientation of particles

in a. flow field is its complexity and computationally inefficiency (see equation 2.17).
PreC:ise numerical representation of a single control volume with the fiber distribution
funC tion requires a large amount of memory in order to solve the partial differential
eqllation for general case. On the other hand, a description based on the average

OrleIatation tensor just requires five entries.

It is important to mention that one can calculate the average orientation tensors
when the distribution function on the orientation sphere (Figure 2.1) is known as
mentioned in equations (2.19 and 2.20). It is possible to define a fiber distribution
function, 7,0(p) when only the second order orientation tensor is known but it is not
possible, with only the second order tensor, to recover unique information of the
original distribution, since there might have multiple possibilities to arrange fibers

16

with the same average orientation tensor. For example the orientation vector

(1/2 0 0)

a: 0 1/2 0 (2-27)

(0 00)

can be described as planar random, planar biaxial, or all other planar states with an

 

 

average 1 / 2 by 1 / 2 for both axes. Figure 2.4 shows different orientation states and the
appropriate representation with the orientation tensor[21]. Therefore, planar random
and biaxial orientations will have different distribution functions, even though the
average orientation tensors are the same. A three dimensional random orientation
State provides a uniform distribution function of 1/47r because the surface area of the
uni t. sphere is 47r. For two dimensional random orientation, the distribution function
beC Omes a disc of zero thickness and infinite diameter lying in the planar level of
OrieIntation. In the case of uniaxial alignment it is an infinite long cylinder with zero
diaIrleter directed in the direction of alignment. The state of uniaxial alignment is
the only case for which the exact structure of the distribution function is known from
t11e orientation tensor.
Finally, one can notice that for exact planar or uniaxial alignments, the distribu-
tion function in equation (2.17) goes to infinity with an infinitesimal region because in
e<111ation (2.17) there are coefficients with 1/ sin 6 and 9 can be 0 or 7r and the surface
integral integral of the distribution function over the unit sphere will be 1. Thus, it
is impossible to describe the distribution function numerically for those orientation
states. To describe the distribution function for these cases, it is required to have
unusually high resolutions over the orientation sphere.

17

Uniaxial aligned Planar 30 random

 

 

 

 

 

 

   

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

X random X biaxial
1 1 X1 7
\ X 1‘1; \ 9/
'b

Fl er / k 73‘ /$

\Z /

X2 X2 X2 X2

\ / x.

0 0 0 1/2 0 0 1/3 0 0
Tensor at]: 0 1 0 aij: 0 1/2 0 at]: 0 1/3 0
0 0 0 0 o 0 0 0 1/3

Figure 2.4: Tensor representation of average orientation state (similar to [21])

2 - 3 - 2 Prediction of flow induced fiber alignment

As mentioned above, it is very convenient from a computational view point to

118(3 an alternative approach based on the moments of the distribution function. This
is eSpecially true when additional terms such as the influence for fiber interactions
are included. The derivation of equation (2.17) can be rewritten in terms of average
Orient ation tensors. An evolution equation for the moment of the distribution function
is Obtained by multiplying by pip, both sides of the equation (2.17) and integrating
OVer the spherical surface. The time derivative of the second order orientation tensor

Can be expressed as

D . .
E < pp>=<pp>+<pp>

= <131p>+<(b-na)p>+<pf)1>+<p(i>—151)> (228)

18

Where p, in are defined in section 2.2. The dyadics < {up > and < ppJ > can be

written as

< {up >2 w - a+)\(S - a — 2S : A) (2.29)

and

< ppJ >= a - wT+A(a- S — 2S : A) (2.30)

T118 dyadic < (p — pJ)p > for interaction forces among the particle in equation (2.28)
02111 be written as

27r7r

< (p - pap >= / /(p — pnpwsinededcb (2.31)
0 0

Use of the rotary diffusion model in equation (2.14) provides

27r1r

O O 61L, I
< (p — pJ)p >= (—DR—)psm (9d6dq5 (2.32)

Stlbstituting equations (2.13) and (2.14) into equation (2.28) and utilizing the equa-

thhs (2.29 - 2.32), the following equation can be derived as

Da

.5?=_w-a—a-wT+A(S-a+aoS—2S:A)+2DR(6—3a) (2.33)

Where DR is the rotary diffusion coefficient. The solution of equation (2.33) provides
the evolution of the second order fiber orientation tensor, which is dependent on the
Velocity profile and time. It contains a and A, so that a relation for the fourth order
moment tensor (A) must be defined. This is a so—called closure problem and closure
approximation models must be used to describe the fourth order orientation tensor.

19

2.3.3 Realizable orientation states

The average orientation tensor is a symmetric tensor (equation 2.22). Thus,

all the eigenvalues of the orientation tensor are non-negative and the corresponding
eigenvectors are orthogonal to each other[11]. The normalization condition given
111 equation (2.24) implies that the sum of the eigenvalues of the average orientation
tensor is 1. Therefore, only two of the eigenvalues are independent. For exact uniaxial
alignment the largest eigenvalue becomes 1 and both other eigenvalues equals are
zero. If two eigenvalues are equal to 1 / 2, the tensor indicates a two—dimensional
1‘ anclom, planar biaxial, or planar orientation where fibers are equally orientated in
bOtli directions. If all eigenvalues are 1 / 3, the orientation tensor indicates three-
diIl’lensional random, triaxial, or any other orientation that gives this average (see
Figme 2.4). These three orientation tensors are called the basic orientation states.
When the alignment lies in the same direction as the coordinate axis, the diagonal
63163ljients represent the eigenvalues of the tensor and the off-diagonal elements are zero.
T1163 eigenvalues represent the amount of fibers in the direction of its corresponding
ei gen vector. For an orientation state to be realizable, the eigenvalues of the symmetric
a-\"erage orientation tensor must be non-negative. When all diagonal elements of a

IIlatrix are greater than or equal to zero (a,,- 2 O) or the Schwarz’s inequality[22]

laijlz 2 am..- (2.34)

must be fulfilled for a realizable orientation state.

Here, two approaches to represent average orientation states are discussed. The
first consists of possible combinations of the three eigenvalues, which allows realiz-

20

a(2)
1.0

0.5

FigUIe 2.5: Diagram of the possible states of the orientation tensor in terms of the

_4\

 

a(3)=0.0

a(3)=a(2)

Isotropic

a(1)2 a(2) 2 a(3)

a(1)=a(2)

Planar Isotropic

 

 

Realizable states

Uniaxial Alignment

 

eigC‘tnvalues (a1 + a2 + a3 = 1)

a‘bility [23]. Figure 2.5 shows this result. Every orientation tensor must lie in the

triE3LIigle, which consists of three straight lines and the corners represent the three

basic orientation states.

The second method is based on invariants of the average orientation tensor, which
are independent of the principal axis, while the eigenvalues are dependent on the
direction of the principal axis of the tensor. Therefore, possible combination of the

three eigenvalues results in the realizability condition. Table 2.1 shows the possible

orientation states.

’
1.0 a“)

21

 

 

 

 

 

 

 

Orientation Eigenvalues Invariants Remarks
State A1 A2 A3 11 III
Uniaxial alignment 1 0 0 2 / 3 2 / 9

 

Planar anisotropic l-x x 0 11:2/9+2III 0_<_x§1/2

 

 

 

Planar isotropic 1 / 2 1 / 2 0 1 / 6 —1 / 36

 

Planar axisymmetry x x 1-2x 2III=6(II/6)1'5 1/3gx31/2

 

Isotropic 1/3 1/3 1/3 0 0

 

(”[07:03 10 :lfl

Axial axisymmetry x x 1-2x 2III:=6(II/6)1'5 0Sx31/3

 

 

 

 

 

 

Table 2.1: Boundaries of the realizable orientation states

The anisotropic orientation tensor is defined as
b = a — —6 (2.35)

Wllere 6 is the identity matrix. Three invariants of b can now be defined as I:
tr ( b) =2 0, II: tr(b . b), and III: tr(b ' b - b) [23]. Two nontrivial invariants II and
III are used to describe the realizability of an orientation-tensor, which is shown in
F igure 2.6 in the shadowed region.
Let /\1, A2, A3 are the three eigenvalues of an average orientation tensor a. The
eigenvalues A1, A2, A3 are independent of the choice of the coordinate system and the
Choice of coordinate system does not influence the realizability of the orientation-
tensor. The invariants II and III can be expressed as eigenvalues of the matrix a

as

11: (A1 —1/3)2 + (A2 —-1/3)2 + (As - 1/312 (336)

22

  
   

Planar
anisotropic

 
    

W
orientation states

 

Planar isotropic
Axial Axissymmetric

Planar axisymmetric

Isotropic

F i glue 2.6: The shadowed region shows the realizable orientation state as a function

Of t he invariants of the anisotropic orientation tensor (b = a - %6)

8.1—1(1
111: (A1 —1/3)3 + (x2 — 1/3)3 + (A3 — 1/3)3 (2.37)

As the sum of the three eigenvalues is unity, A3 is chosen to be dependent on /\1
a11d A2. A two—dimensional solution space can be obtained from equations (2.36 and

2.37). Table 2.1 also shows the solution space.

23

 

CHAPTER 3

THE DISTRIBUTION FUNCTION

T11e fiber distribution function is governed by a complex equation for a general flow
field. It is only under very specific conditions that the fiber distribution function
earl be solved analytically. Solutions for other cases must be calculated numerically,
especially when the fiber-fiber interactions are included. In this chapter, analytical

Solutions and a numerical method to solve distribution function are presented.

3. 1 Analytical solution

Dinh and Armstrong[16] showed that the fiber distribution function without in-

ter ‘clctions among fibers (A 1) and random alignment as initial condition, is given

by
,( trim A- )‘i (31)
1}) pa' — 471' “pp '
for a three dimensional case, and
/ 1 ‘1

(0X,- / 8:17,) is the displacement tensor, X,

for a two-dimensional case. Where AU-

the position vector at time t = 0 (referential configuration) and 36,- the position vector
at time if (current configuration). The displacement tensor must satisfy the following

24

differential equation [9]

(1A,? ‘
thJ = (LT/’1’]: + AS“) Akj (3.3)

 

A

vvith initial condition (t = 0) Aij = 6,-1- (Identity tensor), Sik is the rate of strain
tensor, and 2:”, is the vorticity tensor. In this section, the analytical solutions of fiber
distribution function with time, '¢v(0,¢,f) for a few specified flow fields have been

presented taking random orientation as initial condition.

3 - 1 .1 Elongational flow

If the velocity gradient of elongational fiow is[24]

Vu=3 0 —1 0 (3-4)

(0 0 4)

T 1 1Q deformation tensor will be found by solving the differential equation (equation

 

3 - 3 ) when there is no interaction among the fibers (CI = 0) and shape factor (A = 1)

as

{€_2t 0 0\

A: 0 et 0 (3.5)

(0 06‘)

Where it = )T is the dimensionless time where T is the actual time. Using equation

 

 

(3.1) the following solution is obtained

_3
2

i ((6—2t + e’)si112 6 cos2 ()5 + et cos2 6)

4W (3.6)

(’(pi) =

25

Figure 3.1 illustrates how the fiber distribution function evolves with the dimen-
sionless time t. It should be mentioned here that in order to obtain the average fiber
orientation tensor from the distribution function, one must integrate equation (3.6)
over the unit spherical surface (equation 2.19).

There is no easy analytical solution for this integration, and a numerical method
must be therefore used. But this becomes very difficult during the later time steps
because the distribution function has very sharp peaks at certain infinitesimal region
as shown in Figure 3.1. This fact requires very small step size for a numerical in—
teg‘ration of equation (2.19) (i.e. the resolution of the grid of the spherical surface
much be very high). Figure 3.2 shows the average orientation tensor obtained from

nIJ-Iiierical integration. After sufficient time, it is almost impossible to calculate the

a""erage orientation tensor by numerical integration.

 

Figure 3.1: Three dimensional plot for distribution function for uniaxial flow. The

distribution function evolves with time when fibers are randomly oriented initially

(i.e. 1109,30) : 1/47r).
26

 

 

 

 

311
03+
0.6 l
0A~
02«
an and ass
0 I r 1 1
o 5 1o 15 20 25

Time

F i.gure 3.2: Components of average orientation tensor obtained from analytical for-
rI11.1]ation of the distribution function in a uniaxial flow. It is assumed that all fibers

are randomly oriented initially.

3 - 1.2 Simple shear flow

An example of a steady simple shear flow for the three dimensional case is pre-

Sented. Here the velocity gradient is

{010\

Vu=1 0 0 0 (3-7)

(000)

27

 

 

            
  
     
  
  
 

 

c. N.
- .

»

~ Io .

 
  
 
    

. .
e at ’1‘
ofi'cfi:

 

   
    
     
 

r
o
’/

   

.
-
*9
’1
r

(11!,
g l
. Li #, I‘¢
x‘ (,1 '1' a], rut: r:
Q ‘.'"z";'~" ’

‘.":.".v

V
\
’ef‘d
\0

 
  
 

~
w

 

      

Time = 0.3

Figure 3.3: Three dimensional plot for distribution function for simple shear flow. The

distribution function evolves with time when fibers are randomly oriented initially (i.e.

2,0(0, a, 0) : 1/47r).

and the distribution function for this flow with no diffusion (CI = 0) and shape factor

(A = 1) is,

1/J(p,t) = 41 (cos2 9 + sin2 ¢c032¢ — tsin 2¢sin20 + (1 + t2) sin2 ¢sin2o)‘% (3.8)

7
Just like for the elongational flow problem, equation (2.19) can not be integrated
analytically with distribution function in equation (3.8). Thus, numerical methods
are used again. Also in this case, the step size for the numerical integrations must be
extremely small in order to achieve sufficient accuracy.

Figure 3.3 shows how the fiber distribution function evolves with the dimensionless
time t in case of simple shear flow. For this case the step size for numerical integration
ShoUld be very small as it has sharp peaks at certain infinitesimal regions. However,
this ecl‘l.1ation can be solved analytically by a different approach. In equation (2.19),

28

 

 

 

 

 

6111
0.8 ~
0.6 «
0.4 ~
a
0.2 ~ 22
a33
0 I I I I l I I I
0 5 10 15 20 25 30 35 40

Time

Figure 3.4: Components of average orientation tensor obtained from analytical for-
mulation of distribution function for simple shear flow. It is assumed that all fibers

are randomly oriented initially.

29

the unit vectors can be presented at reference coordinate system when the fibers are
randomly orientated initially (i.e.‘zgi: = 1/47r)[5] and the average orientation state can

be defined as

 

27r 1r . . ,
a: / / A’kpkAJ‘p’ 112(6, 3, 0) sin9d0d¢ (3.9)
0 0 Amn A quvag

Figure 3.4 shows the results obtained from equation (3.9) when fibers are randomly

oriented initially and there is no fiber interactions.

3.2 Numerical simulation of the distribution function

There is no easy analytical solution for the distribution function when the effects
of diffusion are included in the flow field. The equation (2.17) can not be solved
analytically for a general flow field especially when fiber interactions are included in

the governing equation. This equation is based on the Folgar—Tucker model (equation
2- 16) for the interaction coefficient[15]. Here, the derived algorithm can easily be
changed for other models for fiber interactions. The numerical solution of the Smolu-
chowski equation (equation 2.17) is based on finite difference technique. To obtain
the average orientation tensor < pip]- >, the equation (2.19) is numerically integrated

over the sphere to get the average orientation tensor a at each time step.

3- 2 - 1 Discretization of the distribution equation and the mesh
used

The general equation of the distribution function is expressed in the spherical
coordinate system. TWO degree of freedoms 6, (f) and the dimensionless time t span

30

 

Figure 3.5: Mesh used for solving the fiber distribution function over the spherical

surface.

the solution space of the equation (2.17). Due to the normalization condition, the

spatial integration of the distribution function over the unit spherical surface is 1

(equation 2.4). Thus, a mesh has been constructed over the unit sphere. As the
distribution function, 1/)(p,t) is symmetry, it is sufficient to model only one half of the
sphere to reduce the computational cost. The spherical surface is divided in n, nodes
in 6 direction and 71,-,- nodes in q) direction. That means the spacing in 6 direction is

A9:

for one axis with 6 = A6 and ends at 7r — A6. These conditions allow to avoid the

7r/n, and in (z) direction Agb = 7r/ 71,5. Figure 3.5 shows that the mesh starts

Sing‘Lllarity as some of the denominator of equation (2.17) contain sin6 and sin 6 = 0
When 6 = 0 or 7r. The other axis goes from 0 to 7r — Aqfi. This mesh is similar to the
coordinate of the earth. The nodes are such that, 2' is in the 6 direction and ii in the
Qi‘direction with n, 1,, indicates the grid point. Thus, the difference to the neighbouring
pOint for each point of the mesh is 66 in 6—direction and (id) in (lb-direction.

31

The Crank Nicolson scheme can be used to solve the Smoluchowski equation

(equation 2.17). Therefore the spatial discretization can be written as[25]

 

/ i'i'_1— I'll — 1’ _.
d1»: = _1_ (at. 2.11 12.1 02-1) (3.10)

do 2 260 + 260

and

 

 

 

(3.11)

2,11 “1:1 , ,1—1 ,11—1 .31 .1 t
d ‘ii’i _ 1 ’li‘>’z'+1 + 9/141 “ 2% "(L/1+1 + i61—1 - 2762'
(102 2

_ war (on)?

The discretization error is given as O ((6a)2, (it?)

3.2.2 Discretized equations

To achieve the highest precision with the smallest number of nodes and higher
time step, the Crank-Nicolson method is used because of its smallest truncation error
(better than both implicit and explicit method). Equations (3.10 and 3.11) have been

substituted in equation (2.17) to achieve the discretized equation. In this Crank-
N icolson method there will be mum-,- equations to solve simultaneously; eventually
it: requires a storage space for memory. When the diffusion coefficient Dr, the fiber
aspect ratio (A), and the flow field (Vu) in equation (2.17) are known, the equation

for node (1,12) can be written as

t t ,t 1 t _
Cla‘n'z'li/H‘Ji + C2.i.ii¢’t+1,ii + C3.i.iiwi.ii+1 + C4,i,iii(’1—1,ii + 05331161314 —

,1—1 )1—1 , 1—1 ,1—1 , 1—1
Dual/11,11 + DZ.i.ii7’r’z‘+1,a + D3.i,ii?61,ii+1 + 043321161431 + 135331161314 (3-12)

for 2 =1— 77.,- and a : 1 — n“.
The set of n,.n,«,- couple equations with same number of unknowns for each time

Step Can be written as the following matrix form (please see appendix B. for entries

32

of each component)
( 01.1.1 02.1.1 Cnvi,nii,1 \

01,1,2 02,1,2 Cni,n1~i,2

 

 

 

 

t
\Cldminii CQ,1,ni.m'i Cni,nii,ni.nii/ me'mii)
( D D D. .. \ f \
1.1.1 2.1.1 mmtul $1.1

,1—1

01,1,2 D2,1,2 011131112 2,1
(3.13)

 

 

 

 

t—l
\D1,1,ni,nii D2,1,m',nii Dni,nii,ni,nii} Kwnimii)

Periodic boundary condition has been applied to calculate the distribution func-
tion 0’): J at the mesh boundaries. Because of the symmetry conditions of the fiber
distribution function, the nodes inside the mesh (top half of the mesh of the sphere)
are found with identical values as the outside neighbor nodes (bottom half of the
mesh of the sphere) of the distribution function in Figure 3.5. Table 3.1 shows the
boundary conditions have been used in this problem.

Equations in (3.13) are solved for every time step. For the first time step, t = 6t
an initial values (t = 0) of the distribution function must be specified. For example,
the distribution function of three-dimensional random aligned fibers is 1/47r for each
“Ode. Although the most distanced off diagonal entries of the matrices are zeros, few
dista-nt off diagonal entities are found to be non zero due to the boundary conditions
in Ta ble 3.1. Therefore, the whole matrix has been solved with LU decomposition[26].

A 0+ ~+~ code has been developed to solve for distribution function for equation (2.17)
TO calculate the average orientation tensor, the error accumulates at two stages

33

 

Node number outside the mesh Node number outside the mesh

 

 

 

 

 

 

 

 

i ii i ii
0 x n, x
n, + 1 x 1 x
X 0 ni+1 — 1 mm
x 11,-, + 1 Ill-+1 — z 1

 

 

 

Table 3.1: Boundary conditions for mesh

first solving the distribution equation (2.17) and second for calculating the second
order average orientation tensor (a) from the distribution function by using numerical
integration of the equation (2.19). A convergence study was performed by reducing
t. he time step and increasing the number of nodes. Convergence criterion is to deliver
the normalization condition (equation 2.4) where the trace of the average orientation
dyadic must be 1. In the calculations, all six different values of the symmetric second
order orientation tensor are calculated. Thus, the normalization condition is a terrific
quality criterion for the good numerical solution.

For a simple shear flow, a good convergence is found for 90 nodes in 6 and 90
nodes in ([5 direction for a time step size for the dimensionless time t = 0.025. Figure
3-6 and Figure 3.7 show that the error of the normalization condition increases with
dim ensionless time. At dimensionless time t = 10, the error is 0.385% for simple shear

flow Without fiber interactions for 90 nodes in 6 and 90 nodes in g0 direction for a time
Step Size of 0.025 (i.e. the trace of the average orientation tensor at dimensionless
time Of 10 is 1.00385). The same is 0.328% when fiber interactions is included. The

34

20 nodes

o\°

.5

g 30 nodes

Ill
45 nodes
60 nodes
90 nodes

 

Figure 3.6: Numerical error of the normalization condition for simple shear flow with

C1 = 0.0, A = 1 while solving Smoluchowski equation (equation 2.17).

error is less when fiber interactions is included because fiber interactions tend to align
the fibers to a random orientation state, which makes the distribution function to be
uniformly spaced over the unit sphere. Therefore, the error in numerical integration
is reduced.
The three eigenvalues of the average orientation tensor for simple shear flow with
fiber interactions when interaction coefficient C1 = 0.01, are shown in Figure 3.8.
Since the fiber interactions are included in this problem, equation (3.1) can not be
used to solve this problem. Here, the fibers are not going to align exactly to the
direCition of velocity vector. Instead it will have periodic movement with respect to

the Velocity vector as shown in Figure 3.8 due to interaction among the fibers.

35

 

 

 

 

 

 

20 nodes
6 .
.\°
.5
§ 4 '
E 30 nodes
2 . 80 nodes
‘i : 45 nodes
; 60 nodes
0 . - . --.._ =—> 90 nodes
0 2.5 5 7.5 10
t

Figure 3.7: Numerical error of the normalization condition for simple shear flow with

CI = 0.01, A = 1 while solving the Smoluchowski equation (equation 2.17).

36

 

 

 

 

 

 

 

 

5 10 15 20

Figure 3.8: The eigenvalues of average orientation tensor for simple shear flow with

the diffusion coefficient 01 = 0.01. All the fibers are assumed randomly oriented

ini t. ially.

37

CHAPTER 4

THE CLOSURE MODEL

Even with today’s high computational power, it is still very difficult to keep track with
each particle in a complex flow field using the fiber distribution function. Therefore,

a method based 011 the moments of the distribution function is preferred to represent
particle orientation inside the flow fields. However, this averaging procedure gives rise
to a closure problem of a fourth order average orientation tensor (i.e. there are more
number of unknown than the number of governing equations). Some of the existing
closure models have been laid out followed by a new class of closure models developed
at Michigan State University. Simulations are performed to study the performance

Of this kind of closure models for different prescribed flow fields.

4- 1 Existing closure models

To predict the average orientation state by solving the evolution equation of the
fiber orientation tensor, the fourth order orientation tensor must be known. The
information available of the fourth order tensor are its symmetry and projection
pr Operties (equations 2.22 and 2.25). It may be possible to reconstruct the fourth

or der orientation tensor from combinations of the lower order orientation tensor,
WitIlcmt prior knowledge of the distribution function. The following hypothesis is

38

commonly made

A = F (a) (4.1)
where a and A are second and fourth order average orientation tensors respectively.
There is no known exact general function for the equation (4.1) to represent A and
numerous models were developed to provide a closure model between A and a. They
are based on different approaches and provide different results of equation (2.33).

Some of them are discussed briefly.

4.1.1 The uncoupling (Quadratic) approximation

The simplest model is the uncoupling approximation of Doi[4] known as quadratic
model. It regards the fourth order tensor as the dyadic product of two second order

tensors as

A =< pppp >=< pp >< pp >41) 24,-ij = aijakl (4.2)
The quadratic model provides exact results for perfectly aligned fiber. It can be
Shown that it does not fulfill all the six symmetric conditions. However, it does fulfill
Only two of the six projection properties (see equations 2.25). This model has the

adV'antage of being very simple.

4- 1 .2 Linear model of Hand

Hand proposed[27] in 1962 to use all products of the second order tensor aij and
the unit tensor 6,). According to the symmetry and projection properties, only the
lineal terms are used[28]. This procedure gives

39

1
Aijkl = —% (5ij5k1 + (M531 + 6i16jk)+

1
? (aijdkl + aikéjl + aildjk + akldij + ajldik + 071,56“) (4.3)

This model fulfills all six symmetry and projection properties. It gives precise results

only for randomly aligned fibers.

4.1.3 Hinch and Leal model

Hinch and Leal derived a number of closure models. Here only one model is
shown. Based on the investigations of Advani and Tucker[29, 30], one model seems to
deliver the best results of this group of closure models, the so called “H&L1”-model.

It. is given as

2 1 3
Aijkl = g (5ijakz + (105k!) — Baud/c1 + 3 (aikajl + ailajk) —
3
2
E Z (62'jakmaml + aimamj6kl) (44)
m=1

Hinchs and Leals model satisfies two of the six symmetry properties, but it does not

Satisfy any of the six projection properties.

4- 1.4 Hybrid model of Advani and Tucker

Since the quadratic model provides a correct answer for perfectly aligned fibers
and the linear model is exact for random alignment, Advani and Thcker[29] combined
the linear and the quadratic model to obtain

Aijk, = (1 — 27 det a) aijakz + 27 det a((6ij6kl + (M631 + (Siléjk) +
1

§ (aij5kz + (Ii/4511 + a¢z5jk + akléij + a316,}, + ajk6,,)) (4.5)

40

This model is identical to the linear-model for three-dimensional random aligned
fibers. In case of uniaxial or biaxial alignment it becomes identical to the quadratic-
model because (det a) will be zero. Only two projection properties are fulfilled. This
model is used in most computer simulations concerned with predicting flow-induced

fiber alignment.

4.1.5 Orthotropic closure model

Orthotropic closure models[13] are the latest family of models from the group of
Tucker. It consists of using the orientation tensor in its diagonal form. Because of the
normalization condition, the diagonal second order orientation—tensor has only two

independent entries an and (L22. The fourth order tensor also satisfies the symmetries
Aijkl = Ajikl = Akjil = Aljki = 14::ka = Ailkj = Aijlk (45)

which can be rewritten with the contraction notation of Jones[20] as
an = Aijkl (4.7)

The contracted notation is shown in Table 4.1.

Therefore, a six by six matrix is obtained as
(31181213130 0 ol
321 322 B23 0 0 0

B31 B32 B33 0 0 0 '
an: (4'8)

 

 

 

 

 

 

 

 

 

Contracted notation m or n Tensor notation z' j or [cl
1 11
2 22
3 33
4 23 or 32
5 13 or 31
6 12 or 21

 

 

 

 

Table 4.1: Contracted notation for Orthotropic closure model

The additional symmetry of the fourth order tensor (equation 46) indicates that
an has nine independent entries. Using the fact, that the fourth order tensor is
symmetric with respect. to any pair of indices (equation 2.23). Thereby B12 = B66,
323 = 844, and 813 2 B55. This leaves six independent entries for the contracted

Inatrix an. There is also a normalization property, which implies that

3
aij = E Assij
s=1

therefore
Bu + 355 + 366 = (111 (4-9)
B22 + B44 + 366 = 022 (4-10)
333+B44+Bss 20.33 = 1—011—022 (4-11)

These three equations (4.9-4.11) and six unknown end up with three independent
entities of the fourth order orientation tensor. Finally, there are three independent

42

entities in the fourth order orientation tensor. They are function of two independent

entities of the diagonal second order tensor as

311 = f11(011,a22) (4-12)
B22 2 f22fa'11, 022) (4'13)
B33 = f33((111,0«22) (4-14)

Cintra and Tucker[13] present three different ways of fitting values for B11, B22
and 833. The “orthotropic smooth closure” interpolates linearly between the points
of three dimensional random, planar random and uniaxial orientation (Figure 2.5).

I nterpolating linearly between the states the three unknowns can be represented as

B11 —0.15 +1.15011— 0.10022
322 = —0.15 + 1.15a11 + 0.90a22 (4-15)
B33 0.60 — 0.606111 — 0.60a22

For the quadratic approximation Cintra and Tucker used the data from exact dis-
tribution function for different set of velocity fields and evaluate the exact orientation
terlsor (1,-3- and an at certain time. The data are then fitted with least-square method
to achieve the fourth order orientation tensor as quadratic function of an and (122.

This closure was named the “Orthotropic fitted closure”. To get representative values
for injection molding they decided to set A = 1 and to set the fiber interactions coef-
ficient C1 = 0.01. They used a simple shear flow, two different shearing/stretching,
an I—111iaxia1 elongational and a biaxial elongational flow for the fitting procedure. The

43

result is as follows

{Bll\

B22 2

was)

 

 

{0.06096 0.0371243 0.555301 —0.36916 0.318266 0.371218\

0.12471 —0.389402 0.258844 0.086169 0.796080 0.544992

 

 

(1.22898 —2.054411 0.821548 —2.26057 1.053907 1.819756)

4.2 Closure model developed at MSU

The previous models have flaws in their development. To enhance the precision of
predictions of flow-induced particle alignment based on closure models, a new group
Of closure models has been developed[1, 31]. These models are based on satisfying

the six symmetric and projection properties, which can not be disregarded by closure

Solutions in order to provide realizable results.

 

011
2
all
6122

022

K 0116122 )

4- 2.1 Construction of the fourth order orientation dyadic

The hypothesis for this closure is based on expressing the fourth order tensor
only in terms of the second order tensor (equation 4.1). To represent a closed form

of eCauation (4.1), six tetradic operators Pm are necessary. The various combinations

44

 

are represented as

P1: S(6,6)
P2 2 5(633)
P3=S(a,a)

(4.17a)
(4.17b)
(4.17c)
(4.17d)
(4.17e)
(4.17f)
(4.17g)

(4.17h)

To ensure the symmetry of the fourth order tensor one symmetry operator is con-

S(X7 Y) = XY + XYT23 + XYT24 + YX + YXT23 + YXT24

Structed with multiple terms in the form of

Change of the 01’" and Bthindices of the dyadic product of the X, Y tensors.

(4.18)

Where S(X,Y) is symmetry operator of X and Y tensor and XYTO‘fi indicates a

Satisfying the six-fold symmetry condition implies that the permutation of any

45

indices of a dyadic must give the same result. Thus, S (X, Y) delivers always a

dyadic with the required six-fold symmetry. Permuting all six possible combinations

of indices can prove the symmetry of function S (equation 4.18) as follows

S(X, Y) = XUYH + X,k)/,-, + Xqu, + YUXH + KkXfl + YuXk, (4.19a)
(i H j) = inYkl + XjkYil + XjIYki + 365sz + ijXil + szsz' (4-19b)
(i H k) = XrJ-Kz + XMYJ-z + XHYU- + ijxi, + Yk,X,-, + Yklxij (4.19c)
(2‘ H z) = X,,-Y,..- + X,,.Y,.- + Xquj + YIJ-Xk. + mx,.- + )4.ij (4.19d)
(1H ’1‘) Z Xiijz + Xinkl + Xilek + Yikaz + Yinkl + Yquk (4-196)
(J H1) 2 Xilij + XikYzj + Xinkl + YuX/cj + YilXjk + Yilej (4.19f)

(k H1) 2 Xinlk + Xilek + Xilej + Yinlk + YilXjk + Yilej (4-198)

where (H) is used to indicate the indices of the entries permuted.

All higher order terms of Pm, where m > 8, can be expressed by lower order
terms. P6 and P8 can also be represented by 131.5 and P7. This can be explained
by Cayley-Hamilton theorem, which states that a square matrix a satisfies its own

characteristic equation[22].

C(33) = det(:1:6 — a) = 0

:> a.a.a — Ila.a + 123 — I36 2 0 (4.20)

Where I1 = tr(a), I2 = 0.5(I? — a : a), and I3 = det(a) are the invariants, and 6 is the

identity tensor. Thus < pppp > in equation (2.20) can be written as E aiPi. The
i=1

fully symmetry model has been constructed with linear combination of four different

terms, which consist of first, second, third, and fourth order terms of a consecutively

46

as follows

A1 =< PPPP >1: 0111P1+ 012P2 (4-21)

A2 =< PPPP >2: 021131 + 023133 + 024134 (422)

A3 =< pppp >3: (131131 + 034134 + a35P5 + a36P5 (4.23)
A4 =< PPPP >4: 0141131 + 044134 + 0147107 + a48P8 (4-24)

Consequently, the closure model provides the combination of the four equations

(4.21-4.24), and can be written as

4
A = 2 Cu < pppp >n (4-25)

n=1

C" are to be determined for this closure.

4.2.2 Reducible form of the fully symmetric closure model

Each of the four operators given in equations (4.21-4.24) is symmetric. The
coefficients am” are selected so that, the normalization condition (equation 2.25)
for the first three operators are fulfilled. For the last operator (< pppp >4), the
coefficients are defined such that the dot product of any two vector in the last operator

(< pppp >4) ends up being a null second order matrix i.e.

Assij = Asisj = Aissj = Aisjs = Aijss = 0 (426)

The Cayley-Hamilton theorem implies, that the coefficients anm (equations 4.21-
4.24) are only dependent on the tensor a, the identity-tensor and invariants of a. For
the first operator the derivation of the coefficients 011m are shown in the following

paragraph.

47

Taking a is any arbitrary second order average orientation tensor, < pppp >1

should satisfy all the projection properties in equation (2.25). Therefore

3
Z Assij = aij (4-27)
3:1

and

A1 =< PPPP >1: 011131 + 0112132
= 2011(5ij5k1+ 521:5]! + ($151.7) +

012 (50014 + (Maj: + (Sitar-j + aijdkz + (Ii/c531 + 0115/0) (4-28)

Now projection of A1 for the first two vectors (i.e. < p - ppp >1) will be from

equation (4.28) as [32]

< p-ppp>1=a=2ozn(tr(6)6+6-6+6-6)+

(112(tr(5)a + 6 - a + 6 - a+tr(a)6 + a - 6 + a - 6) (4.29)

01'

a = 100116 + 0112 (7a + 5) (4.30)
as tr(a) = 1.
811 = 7012a11+ 0112 + 10am.
a22 = 70112311+ 012 + 100111-
Hence, ( a33 = 7012811 + an + 10011. indicating that 0112 Z —

_ 1
a12 = a21 = 7012312 0111 — —=,'5

a13 = a31 = 7012313.

 

323 = a32 = 70112313-

48

 

 

 

 

 

 

 

 

 

 

an", n=1 n=2 n=3 n24

m=1 —7i0 éa:a 715+%a:]f:—:] %(a-a):(a-a)+%(a-a)2
m=2 % 0 0 0

77123 0 % 0 0

772.24 0 —% —%(a:a)—1 —%(a:a)

m: 5 0 0 (211a)_1 0

m=6 0 0 —%(a a) 1 0

m =7 0 0 0 §

m=8 0 0 0 —%

 

 

 

 

 

 

 

 

Table 4.2: Coefficients of the four fully symmetric Operators (reducible model)

The coefficients for the other operators can be derived in the same way. Table

4.2 shows the results[31, 33].

Since each of the first three operators satisfies the normalization condition (equa-
tion 2.25), a requirement for the fully symmetric parameters (equation 4.25) to fulfill

the normalization condition together is given by

3
Z Cm; and C4 arbitrary (4.31)

i=1

Therefore, only three of four fully symmetric (FS) parameters are independent. These
closure approximations for the fourth order orientation tensor in terms of the second

order tensor have six fold symmetry and fulfill all six projection properties of the

49

exact solution. The four terms of the closure model from Table 4.2 are

 

 

 

A1=< pppp >1: —71—OS(6, 6) + ;S(6,a) (4.32)
A2 =< pppp >2: §5(a : a)S(6, 5) + %S(a,a)—%S(6, (a - a)) (4.33)
A3 = < pppp >3: (710 + 3—25a :((:;:))) S(6,6) — —S(6, (a-a))+
1 4 1
(aza)S(a,(a-a))—§(a:a)S(5,(a-a-a)) (4.34)

A4 = < pppp >4: (3139 - a) : (a'a)+7—10(a : a)2) S(6,6) —

1 1 2
$(a : a)S(6,(a.a))+§S((a.a),(a-a))—?S(6,(a-a-a-a)) (4.35)

4.2.3 Irreducible form of the fully symmetric closure model

In case of irreducible model, all the terms are reduced using the Cayley-Hamilton
theorem. The coefficients dmn (different notation from reducible model to avoid
confusion) are selected so that, the normalization condition (equation 2.25) for all
four operators are fulfilled. Here the calculations for the coefficients for < pppp >4

are shown below. As all the operators satisfy the normalization property, therefore
3
Z Assij = aij (4-36)
3:1
and
A4 =< PPPP >4: dillpl + d42P2 + (51431032 + (5144134 + d45P5 '1‘ (147197 (437)

Here, P6 and P8 can be reduced to lower order by using Cayley-Hamilton theorem.
Now applying projection property of A4 for the first two vectors =< p - ppp >4 will

50

be from equation (4.37) as

< p - ppp >4: a = 106416 + 0242 (7a + 6) + 2643 (a+2a - a) +
(144 (7a - a + (a : a)6) + (345 (5a - a+ (—2+3(a : a)) a + 4det(a)6) +

26-47((1+ 2(a : a))a - a+(2 det(a) — 1 + (a: a))a + 2det(a)6) (4.38)

Now taking (142 = ($43 = (144 = 0 as equation (4.38) consists of only three inde-

pendent equations, the following equations can be obtained

10d“ + 4 det(a)a’-45 + 4det(a)c’r47 = 0
(—2 + 3(a : a)) (345 + 2(2 det(a) — 1 + (a: a))d47 =1

50,45 + 2(1+ 2(a Z a))o’z47 = 0

The coefficients for the other operators can be derived in the same way as in case
of the reducible model. Table 4.3 shows the results for irreducible model.

Since all the four operators satisfy the normalization condition, a requirement for
the fully symmetry parameters to fulfill the normalization condition together can be

given by

4
i=1

where Cms are the coefficients of the irreducible model.

4.3 Estimating the coefficients of the parameters

For the fully symmetric closure, the coefficients for the parameters (equations 4.31
and 4.39) are not universal constant. These are computed by calculating the average
orientation tensor a from the distribution function (equation 2.17) and comparing it

51

 

 

 

 

 

 

 

 

 

 

 

 

cinm 7221 71:2 7123 n=4
, . aza 4det]a] det(a]]3-—4a:a)
7” 2 1 —i $3 ' a Tlo (% —2+3a:a _ —2+3a:a) B
m = 2 % 0 0 0
m = 3 0 % 0 0
__ 2 5 1 5 10a:a
m _ 4 0 _7 7 2—3a:a + B
.. 1
m — 5 0 0 2-3aza 0
m = 6 0 0 0 0
._. 25
m — 7 0 0 0 _§[_3
m = 8 0 0 0 0

 

 

 

 

 

 

 

 

Table 4.3: Coefficients for four fully symmetric operators with irreducible model and

B = 15 — 50dct(a) — 30(a : a)(1 — (a : a))

52

with the entities of a of equation (2.33) using the closure model at different times.

Equation (2.33) can be written as follows

S:A=%(—-ll))—:l—w-a—a-wT+/\(S.a+a-S)+2DR(6—3a)) =K(t)

(4.40)
After solving for distribution function Mp), the average orientation tensor a can be
calculated by multiplying pip,- and integrating over the sphere (equation 2.19). Hence,
the right hand side of the equation (4.40) (K(t)) is known for a specified flow field.
For the fully symmetric closure model, < pppp >,- are also known from Table 4.2 or
4.3 as functions of a and 6. Therefore, equation (4.40) can be written as

(043 :< pppp >4) = K(t) (4-41)

4
:1

for reducible model and

Z (Cy-s :< pppp >,-) = K(t) (4.42)

=1

N .

for irreducible model. Now comparing each entities of K(t) with the left hand side, a
system of four equations can be obtained to solve the 0,8 [34] or C48. The results of
these computations are presented in Figures 4.1 and 4.2 at different dimensionless time
(t = 1T and "y is the rate of strain) for simple shear flow, without fiber interactions,
and taking isotropic state as initial condition. Invariants II and III as defined in
Section 2.3.3, have also been calculated from a.

It appears that the values of C,s or C,s may be related to the invariants. The
sensitivity of the coefficients with respect to the anisotropic invariants (II: tr(b - b),
III= tr(b - b - b) where b = a — 1/36) allows to find the magnitude of change of

53

 

 

 

 

G
9 999?

 

 

 

 

hoammmm
P 99%?
aoammmm

*7

 

 

246810, 246810

 

 

 

,
'FH- .W-O—O—H-O—O—H‘

 

 

 

p peer
eoeawam
C4
0 OO—K-l
hOAmNON

 

 

 

m 4+ - . W.
3415:7’8910 31156778910
t

 

 

 

 

c xxx-fl

 

 

l

H
P 09??

achmmmm
HI

9 999?

hoammmm

 

 

 

 

 

 

345:789610 23406678910
t

Figure 4.1: Graphical representation of the values of 0,8 of the reducible fully symme-
try model that match the numerical results and values of the invariants with dimen-
sionless time (f'y'r). These results are for simple shear flow without fiber interactions

and shape factor /\ = 1.

54

 

 

 

 

 

—O.4r

 

 

 

 

 

1.6
1.2)

«5- 0.8 »
0.4 .

 

 

—0.4

 

 

1.6

1.2»

0.4.»

 

 

—o.4t

 

 

 

Figure 4.2: Graphical representation of the values of Cis that match the numerical

results for the irreducible fully symmetry model and t is the dimensionless time (’77).

These results are for simple shear flow without fiber interactions and shape factor

Azl.

55

the response of the function. This provides some sense of the importance of each
coefficient in the fully symmetry closure model over the average orientation tensor.
The sensitivity of invariant II and III to the coefficient C,- (C1, C2, C3, C4) for reducible

model and C,- (C1, C2, C3, C4) for irreducible model is defined respectively as

 

811 BIII
521' — 301 and S31 — ‘55; (4.43)
and
. 811 - BIII
S,=—7andSi=——.— 4.44
2 80,- 3 .- ( l

The sensitivities have been calculated with simple shear flow without rotary diffusion
(i.e DR = 0 in equation 2.33) or Brownian motion, starting from an isotropic state.
The sensitivity of the coefficients has been studied for different dimensionless times
for the invariants. Results are presented in figures 4.3—4.6 for both the reducible and
irreducible models. The investigation has been performed with simple shear flow
direction without rotary diffusion or Brownian motion, starting from an isotropic
state. Clearly the value of C4 in the reducible model has almost no influence over
the value of the invariants of a. This indicates that only 02 and C3 are required to

obtain an accurate solution.

4.4 Fitted model parameters

This section summarizes the effort made at identifying the Optimal model para—
meters of the fully symmetric model. It can be observed that in the fully symmetric
model there can be many possibilities depending on the choices of the parametric co-
efficients. An optimal value of the parameters means that the results calculated from

56

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.4 0.4
0.2) 0.2)

‘3' o a? o
_0,2.wf —02 ]
e4; A - . . —o.4r -

2 4 6 8 1o 2 4 6 10
t t
0.4» 0.4-
0.2 0.2

d? 0 J3" o------:::;; ....

—o.2,. ----;=:==e: —o.2-
—o.44 —o.4r
‘ is” 4“3‘7“$””7””é‘ ”9‘40 343:: 7 a 3 10

Figure 4.3: Sensitivity of invariant II to the constants Cis ($2, = 811/590,) of the
reducible fully symmetry model. The conditions are the same as Figure 4.1. t is
the dimensionless time (77) and these results are for simple shear flow without fiber

interactions and shape factor A = 1.

57

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.04 0.04»
0.02 0.02
0f 0 ‘ a? 0 4
—0.02] —0.02 \
—0.04~ D —O.O4L A;:::;::::::
2 4 3 8 10 2 4 0 3 10
t t
0.041 0.04»
0.02] 0.02]
«E 0 ‘ a? 0W)
—0.02~ —0.02]
—0.04« In...“ —0.04)
L 2 '4 - 6 3 10 2 4 0 3 10
t t

Figure 4.4: Sensitivity of invariant III to the constants 013 (S3, = 0111/00.) of the
reducible fully symmetry model. The conditions are the same as Figure 4.1. t is
the dimensionless time (#7) and these results are for simple shear flow without fiber

interactions and shape factor A = 1.

58

of

d?

 

004+
0.02,

 

—0.02*
—0.04>

 

 

 

 

r
l

0.04]

0.02‘.

 

—0.02

 

~0.04r

Vvvvvvvvjvvvv'v'vv

 

 

0.04 >
0.02 .

—0.02
—0.04

0.04
0.02

Q

«5‘ 0

-0.02 *
—0.04 '

 

 

 

 

 

 

 

 

 

 

Figure 4.5: Sensitivity of invariant II to the constants as (82, = BII/ (9C1) of the

irreducible fully symmetry model. The conditions are the same as Figure 4.1. t is

the dimensionless time (AN) and these results are for simple shear flow without fiber

interactions and shape factor /\ = 1.

59

 

 

 

 

 

 

 

 

 

 

 

0.02 r——— — - a 0.02
0.01] 0.01»

«15" 0 . ‘3 0
-0.011 ' -0.01 1 ' ' ' ' ' ' ' ' ""'"" "
*0-02 ‘2 ‘4‘“2“*2‘ ‘10 W 2 4‘4 “2“‘12

t t

0.02;-——-— we ~ 0.02
0.01: 0.01

 

 

25‘ 0. a ,3 0
t W]

—0.01 ] _0_01 .

 

 

 

 

 

l
-002 ‘—-- 2 —0.02
2

Figure 4.6: Sensitivity of invariant III to the constants as (83, = BIN/80,) of the
irreducible fully symmetry model. The conditions are the same as Figure 4.1. t is
the dimensionless time (")T) and these results are for simple shear flow without fiber

interactions and shape factor /\ = 1.

60

equation (233) using the fully symmetric closure model fits to analytical or numerical
results of exact solution with the smallest error. To obtain the best parameters, the
following procedures have been taken.

In case of both reducible and irreducible models, the parametric coefficients are
calculated from equation (4.41) at. different time for simple shear flow. The invariants
of the orientation tensor are also obtained from the orientation tensor at the corre-
sponding time for same velocity field. From Figures 4.1 and 4.2, it can be noticed
that the magnitude of coefficients in the irreducible model (Figure 4.2) are very high
near the initial time (i.e. near the isotropic state as random initial condition has been
taken). Hence, the irreducible model is more susceptible to produce error than that.
of reducible model. Therefore, irreducible model has been rejected.

Following the sensitivity analysis it can be concluded that the coefficients are
dependent on the invariants in the realizable diagram. Assuming that the coefficients

are function of the invariants, the C,s are obtained by the regression analysis with

best fit as
01 = 3.46 — 26.811 + 76.9111 (4.45)
CQ 2 1.85 — 2.2001 (4.46)
C3 = 1 — 01 — 02 (4.47)
C4 = 2.45 — 2.20C1 (4.48)

The average orientation equation (2.33) has been solved with the fully symmetric
Closure taking the parametric coefficients as in equations (4.45 to 4.48). The results

are compared with the data obtained from exact solution for simple shear flow without

61

 

.0
01

 

.0
a

.0
m

p
(.0
'V 'Y'T-Y T'V'T—‘T—‘Y-V—fifi'

 

.0

 

 

r—‘vf—‘v— r "Y—1*"‘7-~ ~w—r

0

Figure 4.7: Comparison between the results obtained from fully symmetric closure
and the results obtained from the exact solution for simple shear flow without fiber
interactions. The solid line shows the results obtained by the closure model with
parametric coefficients in equations (4.45 - 4.48) while the dotted points are the exact

values.

interactions among the fibers[5]. Figure 4.7 shows that the closure model with those
coefficients provides good results near the isotropic region. But the fibers do not
approach to the uniaxial position as it should be for this flow field. Therefore, taking

all the coefficients fully symmetric closure also does not provide good result.

Next approach to achieve better solution is done by neglecting the higher order
terms ( < pppp >3 and < pppp >4) in the fully symmetric model. Hence C3 and
C4 are zero. Therefore C1 = 1 — 02 due to normalization condition of the orientation

62

tensor (equation 4.31). Considering the quadratic model (i.e. taking C3 = C4 = 0
and 01 + C2 = 1 ), the values of 02 can be back computed from equation (4.41) as
discussed, and can be written as function of the invariants. The realizable diagram
(Figure 2.6) can be approximated as triangle with three corners even though two sides
are not straight. Any parameters inside the assumed triangle can be approximated
with the Lagrangian shape functions of the triangle[35]. This shape functions for the
realizable diagram are similar to the shape functions for finite element approximations.
As the three corner angles are quiet different from each other, it is better to construct
the shape functions directly from the triangle to avoid the mapping error]36]. The
calculations of shape functions are in the appendix (C.). Figure 4.8 shows the six nodes
on the boundary of the realizable diagram to construct quadratic shape functions.
Each parameters (Cg), and invariants can be represented through shape functions. For
simple shear flow there will be N numbers of back calculated 02 and corresponding
second and third invariants. These calculations have been done by using the known
exact solution for simple shear flow without fiber interactions and compared with the
equation (4.41). The C2 values at the corresponding six points shown in Figure 4.8
are {—4.76738, 0.269573, 0.0693274, 0.275868, 0.398572, 0.830628} respectively. This
method provide pretty good result for the simple shear flow without fiber interactions
as shown in Figure 4.9. But for another flow field like uniaxial flow it does not predict
even close to the result as shown in Figure 4.10. This flow field has been taken as

(10 0\

Va = 0 _1/2 0 (4.49)

 

 

(0 0 —1/2)

63

 

 

Figure 4.8: The FSQ model parameter (Cg) has been back computed at the six
nodes of the realizable diagram by isoparametric quadratic shape functions. The
corresponding (C2) values are {—4.76738, 0.269573, 0.0693274, 0.275868, 0.398572,

0.830628} respectively.

64

 

ll
9 .o .o
N (A) b
T—T—‘T—T '. 1‘ I T— V 'V' 'T—T_Y_ r—T "T“ Y'_I'—7—T T—_Y' V V V' T Y—Y 'V V I V V V 7‘1

 

 

 

 

0.05 0.1 0.15 0.2
III

Figure 4.9: This result is for simple shear flow field. This result is for simple shear
flow field. Comparison between the results obtained from the numerical solution and
the closure model with Cg computed using shape functions. The dotted points are the
plotted from exact result whereas the solid line from the closure model using shape

functions.

65

 

rw——r’——~——i—~ T——r~v _f"'—‘T_—‘_~——.T v v ‘r r T v r r v Y v r r v Y v i

0.6 ~

2.5 ~ “' i

 

0.4 ‘1 . °'

0.3

0.2

 

0.1

'0

 

fif" fifirfifi'fifi‘fi-fi—V’ w—Pr‘Tfi—V_Y—‘ r V V V T
I
O
O

O 0.05 O 1 0.15 0.2

Figure 4.10: This result is for uniaxial flow field. By using the shape functions with
quadratic interpolation of Cg values in the realizable diagram, the results also does
not match to the expectation. The dotted points are the plotted from exact result

whereas the solid line from the closure model using shape functions.

This closure using interpolation functions provides even unrealizable result. Therefore
this model parameters also have been rejected.

An attempt has been taken to achieve simpler FSQ parametric coefficients. The
following two steady state average orientation states, which indicate top two corners

in the realizable diagram have been considered as

(1/200\1 (100)

a1= 0 0 0 .812: 0 0 0 (4-50)

(001/2) (000)

66

 

 

 

 

Under steady state orientation inside a flow field, Da/Dt term in the equation (4.40)

will be zero and it will result in

41%?=—w-a—a-wT+A(S-a+a-S—2S:A)+2DR(6—3a)=0 (4.51)

With no fiber fiber interactions (DR = 0) and A = 1, equation (4.51) becomes

g:_w.a—a-wT+(S-a+a-S—2S:A)=0 (4.52)

Under steady state conditions, the time derivatives of the orientation tensors in the
evolution equation can be computed. The entities of the derivatives of the biaxial

orientation tensor al are

mZMZMZD—(flazD—(flhzmzo (453)
Dt Dt Dt Dt Dt Dt '

and using equation (4.52) with the velocity field mentioned above and using fully

symmetric quadratic model the following equation is achieved

D(a1)23 _ 3 _
Dt _ 7O( 2+4CQ) _ 0 (4.54)

which results in CQ 2 1 / 2. On the other hand, the steady state orientation tensor for

uniaxial alignment. provides

D(az)11 : D(a2)12 : D(a2)13 : D(32)22 = M : M = 0 (4 55)
Dt Dt Dt Dt Dt Dt '

which results in 02 = 1 / 3. Therefore, it is found that FSQ model parameter Cg must
be equal to 1 / 2 for biaxial alignment whereas at the right top corner of realizable
diagram (Figure 2.6) C2 must be equal to 1 / 3. The solutions for average orientation
equation (equation 2.33) taking FSQ closure model with Cg = 1/3 and 02 = 1/2 are

67

 

0.6 r

0.5 ~

 

0.3

 

0.2

0.1

 

 

fi__"——-—r—T—_,—T—7—Yfifi_—T— 7—1 V r T v 7 v v
b
p
. 2
__'______ _
I
O
X C
.
2
‘2
a.
l
u
.
in
a
n
1
1
2
u
n
u
2
4 . 4

O 0.05

Figure 4.11: The calculations are plotted in the invariant diagram for simple shear
flow without fiber-fiber interaction. The dotted points are for the exact solution. The
dark solid line is for FSQ closure with Cg = 1/3 whereas the faded solid line is for

FSQ closure with 02 = 1/2.

shown in Figure 4.11 when simple shear flow is applied without fiber interactions.
The exact analytical solution for the same flow field is also shown. When CQ 2 1/3
the solution provide good results near the random orientation state but after certain
time the solution provides unrealizable result as it crosses the realizable boundary.
In case of 02 : 1/2 even though the solution provides all realizable result, it fell far

short compare to exact solution.

The parametric estimation obtained by the above procedures fails to provide good

results. Since the steady state FSQ parameter 02 falls between 1 / 3 and 1/2, there

68

might be a value of Cg in between 1 / 3 and 1 / 2, which might provide good result.
Now, by deviating the magnitude of the parameter Cg from 1 / 3, it has been noticed
that Cg = 0.37 provides pretty good results for simple shear flow (see Figure 4.12).
This is valid not only for simple shear flow flow but also with different flow field even
including fiber interactions. For example, in case of uniaxial flow, Figure 4.13 shows
that the parametric coefficients Cg = 0.37 and C3 = C4 = 0.0 provide good results.
As observed from the exact solution for uniaxial flow without fiber interactions, the
fibers reach a. steady state condition in the uniaxial region.

The results obtained from fully symmetric quadratic (FSQ) closure with Cg = 0.37
are shown in Figure 4.14 for simple shear flow when fiber interactions coefficient
C, = 0.01. Fiber distribution function has been solved numerically (see section 3.2)
for this prescribed flow field using equation 2.17. Then the distribution function have
been numerically integrated over the spherical surface to obtain average orientation

tensor. It seems that FSQ closure with Cg = 0.37 provides very good result.

4.4.1 Summary

The fully symmetric closure satisfies all the symmetry and projection properties.
The main difficulty is in finding the parametric coefficients (Cis). Different procedures
have been tested to obtain a better closure. FSQ closure with (Cg = 0.37,C1 =
1 - Cg, C3 = C4 = 0.0) has been selected as a good closure model. In the following
chapters, this FSQ closure has been used to simulate the fiber orientation inside

complex flow field using the finite element methods.

69

 

 

 

 

 

0.8 4
0.8 - . .
Analyn cal Solution
0.4 - ’
0.2 - - - - - Solution using FSQ closure
0 l l l l I l
0 5 10 1s 20 25 30

Time

Figure 4.12: Comparison of the maximum eigenvalues with the solution obtained
from distribution function and FSQ(Cg = 0.37) for simple shear flow without fiber

interactions. X-axis represents the dimensionless time w.r.to rate of strain.

70

,2 _ 02:0.33 unrealizable

 

 

 

 

 

 

1 4 Wi— :“ — — —[ 1 ——————
/,..—
0.8 -
0.6 ~ v '
Analytical solution CFO-37 realizable
0.4 J
0.2 4
0 ar— , 1 f 1 —)

Time

Figure 4.13: Comparison of FSQ closure with two different parametric coefficients

(Cg = 1/3 and Cg = 0.37) and the exact solution for uniaxial flow.

71

 

0.6

V I I 1 V
4
L A 1 4A

 

 

 

 

0.5 » J
i 1
i l
04 » 4
= 0.3 : FSQ closure with Q = 0.37 .
: Solution from distribution function ]
0.2 f
i 1
0.1 ~ 1
0 i i

0 0.05 0.1 0.15 0.2

III

Figure 4.14: Comparison with the solution obtained from distribution function and

F SQ(Cg = 0.37) for simple shear flow with fiber-fiber interaction(C1 = 0.01)

72

CHAPTER 5

PREDICTION OF FLOW-INDUCED FIBER

ORIENTATION USING THE FINITE

ELEMENT APPROACH

In the previous chapter, simulations for the fiber orientation have been performed for
prescribed velocity fields. The results for fiber orientation are with a predetermined
velocity field and the forces exerted on the fluid due to suspension of the fiber (see
appendix D. for detail) are unaccountable. In practical applications for manufacturing
complex parts, the flow fields change with time and location. Furthermore, if the fiber
concentration level is semi—dilute or concentrated [37], the fiber stress is comparable
with the Newtonian stress developed by the fluid flow[38]. The concentrations are

defined by the following

u< 171; dilute
312 << 1/ << d—ig Semi-dilute (5-1)
1/ = 0 (fig) - - - Concentrated

Therefore, fiber stress should be included in the governing equations to achieve accu—
rate orientation unless the fiber concentration is very dilute.

For a complex geometry, one must solve the continuity equation, the momentum
equations, the orientation equations, and the equations for a constitutive model for

73

the fiber stress. There is no analytical solution to account for fiber forces and applied
moments, unless the flow field is really simple. A finite element approach has been
used to solve the governing equations for fiber orientations taking into account the
force due to the suspension of the fiber inside the fluid. A finite element C++ code
has been developed for this problem. Several numerical approaches have been taken

to achieve efficient numerical calculations.

5.1 Finite element formulation

In this section the governing equations have been stated for Newtonian-flow, and

finite element formulations have been developed to predict the fiber orientation.

5. 1.1 Governing equations

Mass and momentum have been balanced inside the suspension. Furthermore,
the constitutive equations for the stress due to fiber suspension and the constitutive
equations for fiber orientation have been used for the finite element calculations.

The continuity equation is

V-V=0 (5.2)

where V is the velocity vector. The momentum equation including the stress due to

rotation of the fibers, is given by[39]

p<aa—Y+(V-V)V)+VP—pV2V—V-a'f20 (5.3)

where p is the density of the fluid, P indicates pressure, p is the dynamic viscosity of
the solvent, and of is the stress on the fluid due to presence of the fibers.

74

The constitutive equation for the stress caused by the fibers[38, 40] is given by

(see appendix D. for the derivation of the equation)

of = k A: S+2kDR(a —— 36) (5.4)

The equation for the average fiber orientation (equation 2.33) is again

%+(V.V)a=_w.a—a-wT+A(S-a+a-S—2A:S)+2DR(6-3a) (5.5)

The unknown quantities in these equations are V, P, apand a. It can be noticed
that the above average orientation state equation and constitutive equation for stress
(5.4) are nonlinear. Therefore nonlinear solver should be used for this problem. Fur-
thermore, a is symmetric and the sum of the diagonal terms of a is equal to 1, which

reduces the number of unknown variables in the average orientation tensor[41, 42].

5.1.2 Weak form of the equations

The flow domain QT is approximated by subdividing the domain (0) into bi-
quadratic, isoparametric Lagrangian elements. The weak forms of these equations
(5.2 to 5.5) are obtained via the usual way [35, 43, 44]. The residual equations are
multiplied by a weighting function, and integrated over the domain 9. Thereafter,
the differentiations are distributed among the various variables equally such that the
variables and the weighted function are differentiable only once with respect to the
independent coordinate variable (23,-). This is achieved with integration by parts. The
third step consists in expressing the boundary integral terms as the function of known
quantities. The formulations are derived below.

75

The continuity and momentum equations lead to

/ v14,” (v . V) 210 = 0 (5.6)
Q ,

f [p (Wt-45:5 + WfUPG (v . V) V) + W,“ (VP — V . of) + u (vvvwr) d0
0

= / tW,“dI‘ (5.7)
I‘

where t is the Newtonian traction force on the boundary = 44%}. For the stress

developed[45] in the flow due to fiber orientation, the corresponding weighted equation
is
f w] [of — k A: S—2kDR(a — 36) ] 20 = 0 (5.8)
o

The weak form for the average orientation state is
e

(9
W,151 [0? +woa+a-wT—A(S-a+a-S—2S :A) —2DR(6—3a)] 010:0 (5.9)
where W, are the weighting function and the superscripts u, p, f, a are used to dis-
tinguish the weighing functions for the velocity, pressure and orientation tensor re-
spectively. W,SUPG is the weighted function due to upwinding (see appendix E..)
The above equations (5.8 and 5.9) are nonlinear due to S:A. Furthermore, numer-

ous unknown parameters involve a multiplication of the velocity components V,- with

orientation state (a,,—).

5.1.3 Vector form of the discretized equations

The parameters are discretized by the finite element shape functions[46]. The
same interpolation functions as the weighted function (isoparametric) are used to ap—

76

proximate the dependent variables, velocity, pressure, stress, and average orientation
tensor by Lagrangian polynomials. Hence V, a" are approximated with higher order
interpolation function than P to satisfy the LBB (Ladyzhenskaya, Babuska, Brezzi)
condition[36, 47]. The average orientation state tensor a has been approximated with
quadratic Lagrange polynomials[48, 49]. The discretization has been formulated in

the following way

Ni
V. = 2 8311/3 .—_ Wu’rv, (5.10)
j=1
N2
P: 2 Wfp, = WPTP (5.11)
N3 m
of]. = Z I/Vrf1 (0,1,) = WuTa'fj (5.12)
m=1
N4
(Ii-7‘: 2: W202; = WaTaij (5.13)
m=1

where W“ WP, IVf, W“ are in C1(f2).
Now substituting the equations (5.10-5.13) into the weak formulations (equations
5.6—5.9) the following equations can be obtained

Continuity(one equation)

 

,druT
/ W108” 0v, = 0 (5.14)
szd

Momentum (three for three dimensional problems, two for two dimensional problems)

()V 8W"T 5W"
u uT [SUPG 11T__WPT d0 P
[/ pW w are] at +[/va1 (W, V,) 8:1: dn] V-— [n— 62: ]

uT u 141‘
{fr/“33” .2].,(,.[/,.(6;V 81;; )22]2=/.Wu.2 2...)
Q :1: Q 1: x F

77

 

 

 

 

 

For stress due to fiber suspension
[/ Ill/’“IIVufdfl] afi, — / k W“ [A: S+2DR(a — 36)] d0 = 0 (5.16)
e 0
Fiber orientation equation
r-u ru aakl u u u
[/ w W Tdfl] E— + U W (W, Tv,-) W Tde 221+
n n
/[w'a+a-w—A(S-a+a-S—ZS:A)]df2=0 (5.17)
e

The construction of the element matrix associated equations, (5.16 and 5.17) one
of the terms (V,- or an) was done in such a way that both V,- and aij are treated as
unknowns in the element matrix. For example, equation (5.16) can be discretized
as Mf ((l(1,,-/dt) + K3114- : 0.0 taking all the 0,,- terms inside K31. But discretizing
the equation (5.16) as (Mffi + K33)afj + K3114- : 0.0 (i.e. keeping some V,- terms
inside K 31 and letting aijs are unknown) will help for a better convergence rate. Same
procedure has been followed for equation (5.17).

Using the Galerkin - Finite element method the non-linear equations [50, 51] are

assembled in the following matrix form

dX
M1E + K(X)X = F (5.18)

where X = (z”,:1:p,x3,x7) and X E R", n = 2n.) + 71,, + 312.8 + 372., is the number
of degree of freedoms (for two dimensional formulation, 2n” = velocity DOFs, 71,, =
pressure DOFs, 3n, = three stresses (DOFs) due to suspension of fibers considered
for the 2D flows, and 3n, = the orientation state (DOFs) components for 2D flows) .

The M 1 is the mass matrix and K (X) is the matrix for nonlinear discretization of the

78

 

 

 

 

 

into partitioned matrices and vectors as[52, 53, 54, 55]

 

 

 

 

steady-state parts of the governing equation. The above equation can be expanded

Aii, 0 0 0 513” K11 K12 0 K14] 513” - f”
0 0 0 0 2 24° K21 0 0 0 cc” f"
0 0 Mf 0 dt 2:8 K31 0 K33 0 2:3 f’
0 0 0 0 ] 27 J K21 0 K43 K44 32’ . f’

 

where the entities of the matrix afe provided in the appendix F..

Solution of the above system requires the boundary conditions for velocity and
pressure at one point and initial condition for the orientation state and velocity. A
first order forward numerical scheme for time dependent calculations.

As mentioned above V, up will be approximated with nine node quadratic ele-
ment whereas P will be four node linear element due to the inf-sup condition[36, 47].
The average orientation tensor, a also has been approximated with a nine node
quadratic element since the orientation equation is nonlinear. Various problems are
solved below with different geometries. The numerical results are first validated by

comparing with known solutions.

5.2 Testing the closure model implementation in simple shear
flow

In validating the finite element implementation with the F SQ closure two cases
must be addressed: validating the closure model and testing the finite element imple-
mentation. The closure model is validated by comparing the results with a solution
of the evolution equation (2.33) with the exact solution where the velocity field is

79

specified. The finite element implementation is tested by comparing with the results
obtained from the evolution equation (equation 2.33).

For simple shear flow (Va, = y and Vy = 0), the analytical solution of the distri—
bution function when fibers are randomly oriented initially (i.e.t/J(p,t = 0) = 1/47r)

and A = 1, can be written as[5]

_§
2

111(p)=Zl7; (AT.A : pp) (5.19)

where A is the deformation tensor with respect to time (t) defined as

(140)

A= 0 1 0 (5-20)

(001)

Since the distribution function can be obtained from equation (5.19) utilizing equation

 

 

(5.20), a can be derived using surface integration over the sphere equation (2.19).

In Figure 5.1 the exact solution is plotted along with the finite difference solution
of the equation (2.33) for simple shear flow. For the solution of equation (2.33) A = 1
and 111(6, q‘),t = 0) = 1/ (4n) has been taken as in case of exact solution. Fourth order
Runge-Kutta scheme with a time step of 6t = 0.0005 was used in the solution of the
equation (2.33) with FSQ closure. The finite element results of the evolution equation
with F SQ closure coupled with Navier—Stokes equation have been superposed in the
Figure 5.1. From Figure 5.1, it can be assumed that the results obtain by the finite
element approximation using FSQ closure is quite good, which validates the C++
finite element code.

80

   

Analytical Solution

 
 

 

 

0.4 J - - - - Solution using FSQ closure
02 - 0 FE Solution using FSQ closure
0 . , , , , .
0 5 10 15 20 25 30
Time

Figure 5.1: Plot of the maximum eigenvalues of the orientation tensor vs dimen-
sionless time as obtained from 3 sources: the analytical solution of the distribution
function, the solution of the evolution equation with FSQ closure, and the finite ele-
ment solution of the evolution equation with FSQ closure coupled with Navier-Stokes

equation.

81

5.3 Flows in complex geometries

The prediction of fiber orientation state inside complex geometries are necessary
for their practical application in manufacturing[56, 57]. The test cases presented
below are (i) pressure driven flow between two parallel plates, (ii) a 3:2 converging
step, (iii) diverging geometry, and (iv) eccentric cylinder. Pressure driven Poiseuille
flow has been considered for simulation of relatively simple geometry. The converging
channel and different angles of diverging channels have been considered because of
their complex flow patterns and practical applications. The eccentric cylinder has
been considered as it provides circulation inside the computational domain and the

behavior of the fiber orientation inside the circulation has been tested[58].

The following calculations have been performed using the parameters associated
with glycerin i.e., density p = 1260kg/m3, and dynamic viscosity u = 1.5Ns/m2.
The volume fraction by volume of the fibers is d) = 0.4%, (a semi dilute concentration

if L/d = 50)[38].

In the following cases the boundary and initial conditions have been specified in
the respective figures. All the fibers in the computational domain have been taken
as initially random (i.e. an = egg = (133 = 1/3 and off diagonal terms are zero).
The calculations are unsteady and need numerous iterations at each time step due
to non-linearity of the finite element formulation. A successive substitution scheme
with a relaxation factor of 0.5 was used. For all the cases, the time step has been

taken as 0.05 with first order backward difference scheme.

82

5.3.1 Approximations introduced in two dimensional calcula-

tions

In case of two-dimensional flow fields, the following approximation has been taken:
the components of 013, egg and 0&3, 053, 033 are set to zero due to their weak con-
tributions. This is justified since initial values of 013,023 = 0.0, the magnitude of
(113,023,053, and 053 are always equal to zero for two dimensional flow field. Even
though the stress component a§3 is not zero for FSQ closure, it has negligible effect
on the velocity field through momentum transfer (equation 5.3) when compared to
the magnitudes of (0111 052, gig). The magnitude of 033 is determined by the normal-
ization condition (033 = 1 —- an — agg). Therefore average orientation tensor a will
have three independent entities (an, egg, alg) as 0.21 = (112 and all + 0.2g + 033 = 1.
Only three variables (0],, 252, Gig) are to be evaluated for the fiber stress tensor in
two dimension. Therefore, the total number of equations for the two dimensional flow
field is nine: one continuity, two momentum equations, three orientation equations,

and three fiber stress equations due to fiber suspension.

The computed microstructure can be interpreted with the help of eigenvalues and
eigenfunctions of a. The eigenvalues of a has the property of 1 / 3 S Amax g 1.0; where
Amax is the maximum eigenvalue in the domain. Also if Amam 2 1 (i.e. Amin 2 0.0), the
fibers are aligned in the direction of corresponding eigen vector. On the other hand,
if Am,x 2 1/3 (i.e. Am,“ 2 1/3), the fibers are oriented randomly in all directions.
Setting emax as the normalized eigen vector associated with the maximum eigenvalue,
the microstructure has been interpreted by plotting (Arm,x — 1/3)€m2x- This will result

83

in a vector of zero length for a random orientation state.

5.3.2 Pressure driven flow between two parallel plates

In this problem the inlet velocity profile has been specified for Reynolds num-
ber(Re) of 36. The results are presented here for 4961 nodes with 1200 nine node
quadratic elements (II/1'2) [36]. It has been assumed that at the inlet (11,.) the fibers
have been injected randomly. Figure 5.2 shows the fiber orientation at different times.
Figure 5.2d can be considered as steady state condition as negligible changes in the
orientation state have been observed after that time.

At the inlet, the dot points indicate that the fibers are randomly oriented. Further
downstream the results show that the fibers become oriented along the flow field
rapidly near the boundary when compared to the center (i.e Amax is higher near the
boundary than the the center region). The velocity gradient (Vu) increase from the
center towards the boundary region. Therefore, the orientation vectors, which are
largely dependent on the velocity gradient, is longer towards the boundary.

The stresses created by the presence of fibers depend on the orientation states of
the fibers. Figure 5.3 indicates that the main eigenvalues of the fiber stress tensor (01)
is higher near the boundary than the center region. In Figure 5.3b, 01 is higher than
at other times because the rate of change of the orientation at that time (7' = 6.0) is
high. The eigenvectors corresponding to the main eigenvalues are also plotted with a
weighting factor of the respective main eigenvalues. The corresponding eigenvectors
are similar to the eigenvectors of the orientation state at that local position. The
pressure fields are also shown in Figure 5.4 at different times. The pressure field is

84

  
     
  
    

_4 , 0.65
- 342 (Mm =6.0 0,60
' ' 0.41
0.39
0.37

(a) Time =2.0

Random orientation at inlet

 

1.00
(C)Time =10.0

Aligned with X axis

 

Figure 5.2: Contour plots of the principal eigenvalues(Am2x) of the orientation vector
superposed with main eigenvectors of the orientation tensor. The results are shown for
four different times and boundary conditions are with velocity field (Vz = 4’; + 0.4y
and Vy = 0.0 at inlet, and Vy = 0.0 and 6%} = 0.0 at outlet) for Re = 36 with

random orientation state initially.

85

(a) Time :20

(b)Time :60 0.22

. 0.1a

Maximum61 =0.l73 0.14
Maximum, =0.224 0.10

. Random orientation at inlet

 

 

 

 

 

 

 

o

o

#1

.0

m

.4 V' ,‘ .
xiv

.2

a)

m

m

A

(c)Time =1o_o 0.22 (d) Time = 15.0 022

~ 0.18 0.18

Maximumct =0.193 ' 0.14 Maximum61=0.159 g: 0.14
0.10 i“: 0.10

Random orientation 0-06 = 0.06

0.01 1.2 _ 0.01

   

Figure 5.3: Contour plots of main (maximum) eigenvalues of fiber stress tensor (01)
are shown at different times. The eigenvectors corresponding to the main eigenvalues
are also plotted with a weighting factor of the respective main eigenvalues. The
results are shown for four different times and boundary conditions with velocity field
(Vac = —y2—2 + 0.4y and Vy = 0.0 at inlet, and Vy = 0.0 and 37V:- 2 0.0 at outlet) for

Re = 36 with random orientation state initially.

86

not same at any downstream cross section. Therefore, the specified parabolic velocity
profile (at inlet) changes in the downstream due to the fiber stress and the fluid no

longer behaves as Newtonian fluid.

5.3.3 Converging channel

A 3:2 converging channel is used to study the evolution of the microstructure in
a more complex flow. The fibers have random orientation at the inlet. The Reynolds
number is 1.3, and 5803 nodes and 1400 nine node quadratic elements are used. Due
to symmetry of the flow field, only the top half of the physical domain has been
considered.

In Figure 5.58 the average orientation states are indicated by the vector (Amam —
1/3)€m2x, where em,x is the normalized eigen vector associated with the maximum
eigenvalue. The solution is shown at dimensionless time 'ymxr = 6.72; where 17mm, 2
maximum rate of strain (0.336 sec—1). At the converging section, the changes in the
velocity field affects the orientation of the fibers as the unbalanced forces tries to align
the fibers with the flow field. Further downstream, they tend to get oriented along
the flow field as the flow progress.

Another parameter of interest in interpreting the microstructure is Cos(x), which
is defined as the dot product between the normalized orientation vector corresponding
to Arm,x and the normalized velocity vector. Cos(X) indicates the local alignment of
the fibers with the velocity vector. A contour of Cos(x) is shown in Figure 5.5a.
At the inlet as fibers are randomly oriented the values of Cos(x) are much less than
unity. In Figure 5.5b the main eigenvalues of fiber stress tensor (01) are plotted. Also

87

. 4.00 b . 4.00

(a) Time :20 32° ( )Tlme =6.0 ‘ 3.20

$33 E 3'28

Random orientation at Inlet 030 1.2 E 0.80

   

P P
400 . 4-00

(C) Time =10.0 3.20 (d) Time = 15-0 3.20

2.40 £2. 2.40

1.60 g; 1.60

1.2 0-80 1.2 22 0.80

 

 

Figure 5.4: Contour plots of the pressure field are shown for four different times. The
boundary conditions with velocity field (VJ: = 4’23 + 0.43; and Vy = 0.0 at inlet,
and Vy = 0.0 and a—g’x—z = 0.0 at outlet) for Re = 36 with random orientation state

initially.

88

the eigenvectors with the weighting factor of their corresponding eigenvalues indicate
that the orientation states have great influences in determining the fiber stress in the
suspension. At the converging section, as the velocity vector changes its orientation,
the fibers also undergo changes in their orientation resulting in the higher fiber stress
at the converging section. The pressure field, which accounts for the fiber stress as
well as the Newtonian stress of the fluid suspension, is also shown in Figure 5.5c.

At the corner of the converging section, the continuity equations does not con-
verge to a satisfactory level, which results in an erroneous velocity gradient. Mesh
refinement study has been performed at the corner. But a couple of elements fail to
satisfy the continuity even with very fine mesh. As evolution equation of the average
orientation state (equation 4.18) is very sensitive to the velocity gradient, 8 small
corner region does not provide a realizable state of orientation tensor in the finite

element calculation (i.e., one of the eigenvalues of a is negative. Figure 5.5a).

5.3.4 Diverging channels

The effect of the geometry of the orientation states is studied by solving the
flows in diverging channels with different diverging angles[59]. At the inlet, a random
orientation state is specified for all different diverging angles (Figure 5.6) and the
inlet velocity profile is same. As the flow progress from the inlet, the fibers tend to
align themselves with the velocity vector. At the diverging section, their orientation
changes rapidly due to change in the velocity profile. Immediately after the diverging
section, some fibers are oriented in the vertical direction (along y axis) because at the

diverging section the y component of the velocity field is prominent. Further down

89

(a)

    

Random orientation Unrealizable orientation state

0 . . 0.3

Max = 1.16
(b) M... ‘3}. .02

 

 

 

 

Figure 5.5: The top figure shows a contour plot of Cos(x) superposed with main
eigenvectors of the orientation tensor. The contour plot of the main eigenvalues (01)
of the fiber stress tensor as well as the corresponding eigenvectors with weighting
factor of the main eigenvalues are plotted in the middle figure. The pressure field is
plotted in the bottom figure. The specified velocity field is (Vr = —2(y2 — 0.0842)

and Vy = 0.0 at inlet, and Vy = 0.0 and 0—8‘2—1 = 0.0 at outlet).

90

stream, the fibers become oriented themselves along the velocity field as shown in

Figure 5.6.

In Figure 5.6 Cos(x) is within 0.00 — 0.5 at the diverging section. Fibers in
that region are undergoing changes in the orientation due to change in the velocity
directions. Downstream, the fibers are aligned with the velocity field, which results in
a Cos(x) magnitude near about 1. Another observation at the converging section in
Figure 5.6 is that an increase in diverging angle tends to align the fibers along y-axis
at the diverging section. In case of diverging angle of 40°, the unrealizable orientation
tensors are found. This is due to the same reason as in case of converging channel

(the continuity does not satisfy at that corner by the finite element analysis).

In Figure 5.7 the contour plots of two eigenvalues of the fiber stress tensor and
the pressure profile are shown for 200 diverging channel. The third eigenvalues will
be negligible as discussed before (section 5.3.1). The eigen vectors corresponding to
main eigenvalues are also shown with weighting factor of the main eigenvalues. The
streamlines of the flow field are also shown. The main eigenvalues of the fiber stress

tensor are much higher than those of the second eigenvalues (01 >> 0g).

Figure 5.8 shows the eigenvectors of the corresponding main eigenvalues of the
average orientation tensors with weighting factor of the main eigenvalues for the
different diverging angles. Also comparing with the orientation vectors in Figure 5.6,
it can be concluded that the main eigenvalues of the orientation tensor have great
influences on the fiber stress tensor. Figure 5.8 shows that a reduction in magnitude
of 01 at the diverging section as the angle increased.

91

(a) e = 20°

 

 

 

 

Figure 5.6: Contour plots of Cos(x) superposed with main eigenvectors of the ori-
entation tensor with velocity field (Vx = —2(y2 — 0.152) and Vy = 0.0 at inlet, and
Vy = 0.0 and 87v; = 0.0 at outlet) for different diverging angles. Number of nodes
and elements diverging cones 20°, 30°, 400 are 6813, 7075, 7183 and 1642, 1707, 1734

respectively.

92

 

3.04
2.39
1.74
1.08
0.43

lllffifiv

 

 

Figure 5.7: The contour plots of two eigenvalues of the fiber stress tensors (third
component is negligible) are shown. The eigenvectors corresponding to main eigen-
values are also shown with weighting factor of the main eigenvalues. The streamlines
and pressure field of the flow field are also shown subsequently. The velocity field is
(V2 = —2(y2 — 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and % = 0.0 at outlet)

for different diverging angle of 20°.

93

 

 

 

 

 

Figure 5.8: The contour plots of the main eigenvalue of the stress 01 as well as
the corresponding eigenvectors with weighting factor of the main eigenvalue. The
corresponding eigenvectors are similar to the eigenvectors of the orientation state at
that local position. (Var = —2(y2 — 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and

66:1 = 0.0 at outlet) for different diverging angles

 

94

5.3.5 Flow between eccentric cylinders

The classic problem of the flow between two rotating eccentric cylinders[60] is
also studied with 1680 nine node quadratic elements and 6880 nodes and the outside
cylinder is rotating. The boundary conditions are shown in the caption and Reynolds
number is 0.336 as shown in Figure 5.9. The fibers have been assumed to be randomly
oriented initially. A recirculating region is generated just outside the inner cylinder in
the large gap as shown where the fibers are randomly oriented (Figure 5.9). Outside
the circulating region, the fibers tend to orient with the velocity field. Since the
velocity gradient is high in the smaller gap between the cylinders, the fibers tend to
be highly aligned. However, with a Reynolds number more than 3.0, the continuity
equation does not converge quite well in that small region, where velocity gradient is
very high, resulting in an unrealizable orientation states.

The performance of the FSQ closure is compared with the quadratic and hybrid
models in this problem. The quadratic approximation of the fourth order orientation

tensor is (see section 4.1)

A =< pp >< pp > (5.21)

where < pp >< pp > is the dyadic product of pp and pp. The approximation of

the hybrid model is defined as (see section 4.1)

A = 27deta < pppp >1+(1— 27deta) < pp >< pp > (5.22)

In Figure 5.10 the orientation vectors of all three models have been shown. The
maximum eigenvalue predicted by FSQ closure is 0.87, which is lower compare to the

95

Fsa MODEL )‘rn - 0-87 Fso MODEL

      

Figure 5.9: Contour plots of eigenvalue(A) superposed with main eigenvectors of
the orientation tensor. Also the main eigenvectors of the fiber stress tensor. The
streamlines with the pressure field of the flow field are also shown below. Re =

”—le = 0.336 in eccentric flow (outer cylinder is rotating and the inner one is fixed).

96

01 FSO MODEL )Vma? 0-87 QUADRATIC MODEL Amay? 0-98
. ' _ 0.1 ’- .
: ‘ km"? 0.33 . Ami”: 0.33

0.05 0.

O
01

>'0

-005 -0.05 i

 

 

   

-0.1 005 0 -. 0.05 0.1 -0.1
X

HYBRID MODEL 4mm? 0-97
A . = 0.33

min

0.11-

     

0.05:

-0.05:

 

Figure 5.10: Comparison of the FSQ model with the Advani Tucker hybrid model for
realizable orientation state (i.e. all the eigenvalues are positives) with R8 = 6% =

0.336 in eccentric flow field (outer cylinder is rotating and the inner one is fixed).

97

FSQ MODEL Xmas; 0-83 QUADRATIC MODEL Amaf 1-40

0.05

 

 

 

Figure 5.11: These test cases for courser mesh with Comparison of the FSQ model
with 2108 nodes and 496 elements. Advani Tucker hybrid model for realizable orien-
2

tation state (i.e. all the eigenvalues are positives) with Re = (”73 = 0.336 in eccentric

flow field (outer cylinder is rotating and the inner one is fixed).

98

other two models (0.98 and 0.97 for quadratic and hybrid model respectively). How-
ever, with the courser mesh of 2108 nodes and 496 nine node quadratic elements the
results show that the orientation state for the quadratic and hybrid models are unreal-
izable in some region (Figure 5.11), whereas FSQ model is still realizable throughout
the region. The FSQ model appears to be more robust but the accuracy of the results
still need to be validated.

The eigenvalues of the stress tensor also have been compared with each other
(Figure 5.12). The eigenvalues of the stress tensors of FSQ and quadratic are very
close to each other unlike the hybrid model. For a complex flow field (unlike simple
shear flow; see appendix D.), the quadratic model provides very good results. The
pressure fields, which are influenced by both hydrodynamic stress as well as fiber
stress, are shown in Figure 5.13 for different closure approximations. The pressure
profile of FSQ and quadratic are very close to each other unlike the hybrid model,

which is 30% higher than the FSQ model.

5.4 Summary

Table 5.1 summarizes the case study in a tabulate form. Table 5.1 shows com-
putational time required for the each case to calculate upto 400 time steps with step
size of 0.05. Increase in node numbers and Reynolds number requires higher compu—
tational time. As the FSQ model is more complicated than other the two models, it

requires more computational time than the quadratic and hybrid model (Table 5.1).

99

F

       

SQ MODEL
0.1 . ..

-0.05 '

 

'- . -0.05

HYBRID MODEL

0.19
0.15
0.11
0.07
0.03

0.1
0.05'
> 0 .

-0.05

 

 

Figure 5.12: The main eigenvalues of the stress tensor for three different models have
been plotted. The FSQ and quadratic models provide results, which are close to each

other unlike the hybrid model.

100

FSQ MODEL QUADRATIC MODEL
pm; 0.95

= 0.95
Pmax

  

0.1
0.83
0.05 0.49
0.15
-o.19
; m; -0.53
-0.05 , -0.87

 

HYBRID MODEL
p = 0.93
max

  
 
 
 
 
 
  

0.83
0.49
.1 0.15
, -0.19
-0.53
-0.87

Figure 5.13: The pressure contours for three different models have been plotted. The

profiles are almost similar.

101

 

 

 

 

 

 

 

 

 

 

Test Case Re A 0'1N / m2 CPU Time
Poiseuille flow. Figure (5.2 and 5.3) 36 0.85 0.224 43.146
Converging channel (Figure 5.5) 1.3 1.45 1.16 696.618
Diverging channel (200). Figure 7.56 0.88 0.032 222.668
(5.6 and 5.8)
Diverging channel (300). Figure 7.56 0.92 0.243 211.005
(5.6 and 5.8)
Diverging channel (300). Figure 7.56 0.93 0.190 284.507
(5.6 and 5.8)
Eccentric cylinder (FSQ model and 0.336 0.83 0.14 1197.094
Figure (5.9)
Eccentric cylinder (Quadratic 0.336 1.4 0.15 1115.969
model and Figure 5.11)
Eccentric cylinder (Hybrid model 0.336 1.03 0.20 1130.018
and Figure 5.11)

 

 

 

 

 

 

 

Table 5.1: The summary of the case studies are presented. A and 01 indicate the
maximum eigenvalue of the average orientation tensor and fiber stress tensor respec-
tively. The computational time ( 1000 sec as an unit) is also shown for each case on

the right most column.

102

CHAPTER 6

PARAMETRIC STUDY

In this chapter the effects on the orientation state and fiber stress have been studied
by varying the concentration of fiber, aspect ratio, and diffusion effect of the fibers.
Only one complex geometry (cone 200 diverging channel) has been considered for this
evaluation. Low diverging angle also provides good convergence. The FSQ model has

been used for all the simulations.

6.1 Effects of varying the fiber aspect ratio

The longer is the length of the fiber, the larger force requires to move it. From
equation (D.19) it is clear that the fiber stress ( of, = k Aijkz SH ) depends on the

l

coefficient k where k is defined as fi% (8). Depending on the fiber length and

fiber concentration[38] and f (e) = 1715—6, where 5 = [ln(2L/d)]"l, the magnitude of
the coefficient (k) will vary. Therefore, the fiber stress depends on the fiber aspect

ratio, concentration of the fibers in the suspension, viscosity of the fluid.

From the Figure 6.1 the fiber stress is highest for an aspect ratio 40. Figure 6.2
shows that the fiber orientation state does not change much with the relatively large
aspect ratio. The aspect parameter A [(A = (0,2, — 1) / (0,2, + 1) )] in the orientation
equation (equation 2.33) depends on the aspect ratio. But the magnitude of A does

not change much when aspect ratio is increased from 20 to 40.

103

0.5 5

(a) ap= 20 Max 6, = 0.032

0.4 -

 

 

 

 

 

 

 

 

Figure 6.1: The effects of varying the fibers aspect ratio while there is no fiber-fiber
interaction. The eigenvalues of fiber stress tensor are clearly high for longer fibers.
The velocity field is (V3: = —2(y2 — 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and

595—: = 0.0 at outlet) for different diverging angles.

104

max - ' Cos(ki)

 

 

 

Cos(ki)
1 .00

 

 

Figure 6.2: The effects of varying the fiber aspect ratio on Cos(x) while there is no
fiber—fiber interaction. The fiber orientation state does not change much with the
aspect ratio of the fiber. The velocity field is (Va: = —2(y2 — 0.152) and Vy = 0.0 at

inlet, and Vy = 0.0 and % = 0.0 at outlet) for different diverging angles.

105

6.2 Effects of varying the aspect ratio with diffusion

Fiber interactions have the effects to orient the fibers randomly (i.e. isotropic
state). Therefore, the eigenvalues of the average orientation state will be less when
fiber interactions are accounted for, than without fiber interactions. Comparing fig-
ures 6.3 and 6.2, the Amax is less by 10% when fiber interactions are effective. From
equation (5.4), it is clear that the diffusion will contribute to the fiber stress. So the
fiber stress is little higher when diffusion is effective. Therefore, the fiber stress in

Figure 6.4 is higher corresponding to Figure 6.1.

6.3 Effects of varying the concentration of the fibers

V

The fiber concentration ¢f is defined as (bf = VL = where V“, = fiber,

VFVL,
fluid and suspension volume respectively. In the average orientation state (equation
2.33) does not depend on the concentration of the fibers. Figure 6.5 shows that the
average orientation state does not change much with different concentrations. But
from equations (5.4), it is clear that the fiber stress depends on the concentration of

the fibers. Figure 6.6 shows that the fiber stress increases with the increase of fiber

concentration in the suspension.

6.4 Effects of varying inlet velocity

The inlet velocity profile in Figure 5.6 is Va: = —2(y2 — 0.152) and Vy = 0.0.
An increase in the inlet velocity profile means an increase in the velocity gradient.
As the rate of strain of velocity field influences to predict the average orientation of

the fibers, the increase in velocity gradient tends to align the fibers more along the

106

0.4 km, = 0.7996

 

 

 

 

Figure 6.3: The effects of varying the fiber aspect ratio with diffusion on Cos(x).
The eigenvalues of the average orientation state is with lower magnitude than the
corresponding eigenvalues without diffusion. The velocity field is (Vac = —2(y2 —

0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and % = 0.0 at outlet) for different

diverging angles.

107

0'5 f (3) ap = 20 Max 0, = 0.038 01

0.4 -
0.12

0.10
l ’i 0.07
0.05
0.02
0.00

 

0.2 -

 

   

0.1 5.

 

 

 

0.12
0.10
0.07
0.05
0.03
0.00

 

 

 

0.6 g—(C) 3,, = 40
Max 0'1 = 0.12 0',

0.12
0.10
0.07
0.05
0.03
0.00

 

 

 

Figure 6.4: The effects of varying the fiber aspect ratio with diffusion on the principal
fiber stress. The fiber stress is clearly higher than those without diffusion. The
velocity field is (Vx = —2(y2 — 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and

% = 0.0 at outlet) for different diverging angles.

108

 

 

 

0-6 (C) (D =0.005

 

Figure 6.5: The effects of varying the concentration of the fibers in the suspension on
Cos(x). The average orientation state does not defer much with the concentration.
The velocity field is (Va: = —2(y2 — 0.152) and Vy = 0.0 at inlet, and Vy = 0.0 and

86—? = 0.0 at outlet) for different diverging angles.

109

0.5 -

0.4 - (a) q) = 0'001 Max 01 = 0.04

 

0.20
0.16

 

0.08
0.04

 

 

0.20

0.12

 

 

 

Figure 6.6: The effects of varying the concentration of the fibers in the suspension on
the principal fiber stress. The fiber stress does change with the concentration of the
fiber suspension. The velocity field is (Vac = —2(y2 — 0.152) and Vy = 0.0 at inlet,

and Vy = 0.0 and M = 0.0 at outlet) for different diverging angles.
01:

110

streamlines. Therefore, the orientation vectors (i.e. the maximum eigenvalues) shown
in Figure 6.7 are longer as the inlet velocities increase. From equation (5.4) it is clear
that the fiber stress is directly proportional to the velocity gradient. Consequently,

the Figure 6.8 shows that the fiber stress increases with the increase in inlet velocity.

6.5 Summary

Table 6.1 summarizes the parametric study in a tabulate form. Table 6.1 shows
the computational time required for the each case to calculate the 400 time steps
with step size of 0.05. Increase in aspect ratio, concentration and Reynolds number
requires higher computational time. Also adding diffusion on the simulation reduces
the computational time because the maximum eigenvalue of the orientation tensor

will be less, which makes a better convergence.

111

(a) V=1.0*norm 7me= 0-88

 

 

0.5 (b) V=1.5*nonn

 

Cos(ki)
1 .00

0.00

 

Figure 6.7: The effects of varying the inlet velocity profile on Cos(x). The average
orientation state with different inlet velocity profile. The norm inlet velocity profile

is (Vac = -2(y2 — 0.152) and Vy = 0.0)

112

0.5 -
04 (a) V=1.0*norm Max 0, =0.032
‘ F 0.35

0.28
0.20
0.13
0.06

 

   

 

 

O ‘ 072 ' 074 die ' ‘ofs i ‘ ”1:2 " '174' I 1.6

(b) V=1.5*norm

0. -

5 Max 0, =0.222 01
0.35
0.28
0.20
0.13
0.06
-0.01

 

 

0.35
0.28
0.20
0.13
0.06
-0.01

 

 

Figure 6.8: The effects of varying the inlet velocity profile on the principal fiber stress.

The norm inlet velocity profile is (VJ: = —2(y2 — 0.152) and Vy = 0.0)

113

 

Varying Parameter A 01 CPU Time

 

 

 

a.p = 20 0.8822 0.032 222.668

 

Different aspect ratios (Figure 6.1 a,D = 30 0.8848 0.064 242.684

and 6.2)

 

a7) = 40 0.8867 0.106 256.717

 

 

up = 20 0.7996 0.038 211.005

 

Different. aspect ratios with fiber- (LP 2 30 0.7999 0.074 229.027
fiber interactions. (Figure 6.3 and

6.4)

 

a.p = 40 0.7993 0.120 252.194

 

 

45 = 0.001 0.885 0.040 225.961

 

Different concentrations. (Figure 6.5 d) = 0.002 0.887 0.083 251.663

and 6.6)

 

¢=0.005 0.890 0.220 284.677

 

 

Re = 7.56 0.88 0.032 222.668

 

Different Reynolds number. (Figure Re = 11.34 0.91 0.222 279.731

6.7 and 6.8)

 

Re=15.12 0.94 0.356 294.485

 

 

 

 

 

 

 

Table 6.1: Parametric study of the cone 200 diverging angle. A and 01 indicate
the maximum eigenvalue of the average orientation tensor and fiber stress tensor
respectively. The computational time (1000 sec as an unit) is also shown for each

case on the right most column.

114

CHAPTER 7

USE OF GREENS FUNCTION TO INCREASE

COMPUTATIONAL EFFICIENCY

This work summarizes on efforts directed at increasing the computational efficiency
for the numerical solution of second order differential equations. Many of the gov—
erning equations in mechanics problems involve with second order partial differential
equations. Analytical solutions to these problems are available only for very simple
geometries. Even with high computational power, it may still require a few days to
solve a complex problem. There is still a need for further development of computa-
tional strategies that will reduce the cost of such computations. A possible approach
toward reaching this goal is based on combining analytical and numerical solutions.
In the proposed method, simple sub domains, where the Green’s function is easy to
construct, have been taken inside the physical domain. Combining the analytical so-
lution using Green’s function with finite element calculation of outside the analytical
domain (Figure 7.1), a unique solution of the problem has been calculated without
requiring functional node inside the circle. The computational times and memory
requirements for both the proposed and Galerkin finite element methods have been
compared.

The theoretical foundation of the proposed method is presented first. The advan-
tages over the Galerkin finite element calculations are then discussed. Heat transfer

115

No mesh
needed

Analytical

solution inside
the circles

 

Figure 7.1: Outline of the mesh required for the proposed combined finite element

and analytical method for a complex physical domain.

and fluid mechanics problems are solved with both the Galerkin finite element method
and the proposed method. The efficiency of the method is discussed in regard to the
required memory and computational time for the solved problems. The possible
source of errors and limitations of the method developed by the proposed method are

also discussed later.

7.1 Theory

The theory of Green’s function is presented to achieve an analytical solution,
which has been discretized to combine with the finite element weak formulation. If
U(r) 6 672(9) is an arbitrary function and r E Q (Q 6 IR" is open, bounded, and U(r)

116

on the boundary (99 is C 1). <I>(:c) is defined as the fundamental solution of Laplace
equation for :1: E R", when :1: 3:4 0. Using Green’s second identity, the following result

can be obtained [61]

[Q My — w)v2U<y)dy (7.1)

where v is the unit outer normal.

Utilizing equation (7.1) and considering the following Poisson problem,

—V2U = f in {2
(7.2)
U = g on an

where g is continuous on boundary, the solution of equation (7.2) is

U<x)=— [angmg (xMHdS / f (7.3)

where G (:12, y) is the corresponding Green’s function. Green’s function is also a funda-
mental solution, but the boundary value of the Green’s function is zero (i.e. <I>(a) = 0
on the surface). The advantage of the zero boundary value of Green’s function has
been used to achieve equation (7.3).

A simple Green’s function is necessary for the solution of equation (7.3) to reduce
the computational time. The construction of a Green’s function for an arbitrary
shape is difficult. For example, the Green’s function of a rectangular geometry can
be expressed as the summation of infinite series whereas the Green’s function of a
sphere is simple to construct. If g on 89 and f in Q are known, the solution of
equation (7.2) can be achieved by using equation (7.3).

117

This approach is tested for a two dimensional problem in this work. Circles
are considered for identifying the subdomains because Green’s functions are easy
to construct for two dimensional case. Taking N number of linear elements on the
boundary of the circle (Figure 7.2), equation (7.3) can be discretized for a subdomain

(Ql)int.o the following equation

N gee) f (—) %%(x,y)d8<y)+
U1“) 2 ‘2 an;
90:...) f (we) ?£($,y)d3(y)

ye+1 ‘I/c

60;
+ a f(y)G(x,y)dy (7-4)
U103) = Zy<x.)K(m.)+ a f(y)G(:c.y)dy (7.5)

where 9(178) is the unknown degree of freedom on the boundary of the analytical
domain. Equation (7.4) is then substituted into the finite element weak formulation

to obtain the solution over the entire domain.

The solution, U1(CE) has been combined with the finite element weak formulation
to achieve a unique solution throughout the computational domain. From equation
(7.4), f (:13) should be known inside the analytical domain to achieve solution to the

problem defined in equation (7.2).

7.1.1 Viscous flow problem

The governing equations for steady state viscous flow are[62]

V2P = V. (g — u.Vu) = M(x) (7.6)

118

 

Figure 7.2: The boundary of the analytical subdomain (521) has been divided into N

elements.

119

and

uv2v = VP + pV.VV (7.7)

where g is the body force and u - Vu is the convective term. Assuming M (x) =
V. (g -— u - Vu), the solution of pressure P(x), and velocity V(x) of equation (7.6)

and (7.7) can be written with the help of equation (7.3) as

 

   

P(X) = - L9 P(y)%%(x,3’)d5(3') — Q M (Y)Gp(x,y)dy (7-8)
and
V(x) =— mafia". yy>dS( )— f(y)G..(x,y)dy (7.9)
anl 5U :21

where f(x) = fi(VP + pV 0 VV).
Considering equations (7.4) and (7.8), the pressure component P(x) on the ana-

lytical boundary can be discretized as[44, 35]

22
M2

P(x.>Ke(x>> — fn M(y)a,(x,y)dy (7.10)

e=1

where Ke(x) = —asjf (My—)6 8U(x, y)dS(y )— f (Ll—1)%;G(x,y)d5(y). Substi-

YC+l—'Ye 601 yc-ye l

tuting equation (7.10) of the analytical domain (Q1) into equation (7.9) and utilizing

equation (7.4), the velocity field can be discretized as follows[44, 35]

N'

V(xemeux»+§j<P<sce>Ke2(x>)— / f<y)c.<x,y>dy (7.11)

91

22
M2

621 e=1
where V(re) and P(re) are the corresponding unknowns on the boundary of the
analytical domain and K e1(x) and / or K e2(x) are correspond to the non-linear terms
(M(x) and f(x)).

120

The derivatives of pressure and velocity can be written as follows

8P(x) 6Ke(x ),,(0G' ,y)
8m. ~Z— 81:, )—/( M(y a: ———yd (7.12)

and

82;) ~i(u-(r )QEAj—(fl)+i(P() ”1:831:20: +>~ /f.-(y)G

6:1 6:1

 

(7.13)

The finite element weak formulation[62] of the governing equations (7.6) and (7.7)

are

[a (VP.VWP)dQ— fa M(x)WPdQ= fr tW“dI‘ (7.14)

8V
/ [p (Wu? + W3”PG (V - V) V) + W“VP + n (VV.VW")] d0
Q

= / tW"dI‘ (7.15)
F

Substituting P(xb), (ff—SS), and ear—(:9 into the equations (7.14) and (7.15), the element
matrix can be constructed.

M (x) and f (x) should be known prior to the solution of equations (7.8) and (7.9)
inside the analytical subdomain ((21). If body force g = 0, the non-linear convective
part (pV - VV) should be known prior to the solution. To solve equations (7.8)
and (7.9), some velocity field is assumed initially and numerically integrated over the
analytical subdomain. The pressure and velocity field are then substituted into the
finite element weak formulation. The solution is achieved until the iterations converge
to satisfy the compatible boundary conditions.

121

If n is the number of functional nodes of an element next to the analytic sub-
domain and N is the total boundary nodes along the analytic subdomain (figure
7.2), the dimension of the element matrix is n x N. For creeping flow approximation

(convective terms can be neglected), Aflx) will consist of only a constant gravity

 

f 8Ke(x) 6K81(x)
61‘"

force. Therefore, after discretization the element matrix (consists o , 6:1:- ,
0

elf—5:99) next to the analytical subdomain will neither depend on boundary values
of the parameters [P(x) and V(x)], nor the corresponding values of inside the ana-
lytical subdomain. Instead, it will be function of some weighted value depending on
the corresponding Green’s functions and analytical boundary, like the Gauss value of
shape function at Gauss points. To include the nonlinear convective terms into the
equations, the Picard or Newton iteration method[26] can be applied as mentioned

above.

7.1.2 Heat conduction problem

The governing equation of heat conduction problem is[63]

g'(x.t) : _1_aT(x, t)

V2T t
(X’ )+ k a fit

(7.16)

where T(x, t) is temperature at x point and time t, thermal conductivity k is constant,
g’(x, t) is the heat generation per unit volume at x point and time t, and a is the
constant thermal diffusivity. For a steady state condition, the governing equation is

similar to a Poisson’s equation as follows(letting g(x, t) = g’ (x, t) / k)

V2T(x, t)+g(x, t) = 0 (7.17)

122

The equation (7.17) is comparable with the pressure equation (7.6) in fluid flow.

Hence, the temperature can be written as follows using equation (7.3)

0G

M) = — A01T<y>55<xy>ds<y>+ [m g(y>G<x.y)dy (7.18)

The transient. solution can be obtained the same manner as finite element calcu-

lations for unsteady process.

7.1.3 Remarks of the proposed method

1. The standard Galerkin finite element method(GF EM) needs to have the mesh
generated throughout the whole computational domain, whereas the proposed

method needs mesh generated only in the finite element subdomain.

2. N 0 extra boundary conditions are needed for the boundary between the analyt-

ical and numerical domains[64, 65].

3. The size of the global matrix will be much less than the GFEM, eventually

reducing the memory required for solving the same problem.

4. This method can deal with non-homogeneous, nonlinear equations.

7.2 Computational results

The proposed method is first tested with simple problems, where the analytic so—
lutions are known. The results obtained by the proposed method have been compared
with standard finite element results. A C++ finite element code has been developed
to achieve solutions for the proposed method and the standard finite element method.

123

7.2.1 Heat transfer

The governing equation for two dimensional steady state condition without any

heat generation is same as Laplace equation[63]

627 1 e: _
0:132 (93/2 —

0 (7.19)

A square region has been taken to test the proposed method. Inside the square
(Figure 7.3), a circle has been considered for analytical calculation using Greens
function. The mesh generated inside the domain is shown in Figure 7.3; left figures
deal with standard finite element calculations, whereas the right side figures deal with
the proposed method calculations. Four-node linear elements have been considered
for this problem. According to the boundary condition, the temperature gradient
will be constant according to analytical solution. In the case of the prOposed method
the temperature distribution inside the circle has been post computed (Figure 7.3).
The computational information is stated in Table 7.1. The error with respect to
the maximum temperature(100°), generated by finite element method might be due
to affine mapping[36] and the round of error[44]. Affine mapping depends on the
determinant and norm of the geometric interpolation matrix and all elements are
not square inside the computational domain. Figure 7.4 shows the same physical
phenomena inside the square region, but with eight-node quadratic elements. Here
the error (shown in Table 7.1) is lower than corresponding linear approximation in
case of proposed method.

When source term is included in the heat equation, which is [63]

WT BQT
. _ + __ + — 7.20
82:2 83/2 g 0 ( )

124

Computational grid Computational grid

 

 

Analytical solution
Inside the circle

 

 

 

 

 

 

 

 

 

Figure 7.3: Mesh and temperature distribution are shown for full finite element and
combined method when top and bottom sides (y = i2.4) are insulated, at left side

(a: = —2.4) T = 50.0, and at right side (r = 2.4) T = 100.0.

125

Finite Element Solution Combined Method Solution

         

 

 

2:- 2:-
>-o_— >0:—
-2 X 1 2 .2

Figure 7.4: Same physical problem as in Figure 7.3 but an eight-node quadratic

element has been considered instead of four-node linear element.

Assuming the heat generation per unit area to be g = 10.0W/m2 and temperature
on all the boundaries to be 25°C, the results are shown in Figure 7.5. The maxi—
mum temperature difference between standard finite element solution and combined
method is 0.42% w.r.to the boundary temperature(25°C).

Complex geometries have also been considered where simple analytic solution
is not available. Figure 7.6 shows the mesh for both the computational methods.
Equations (7.19) and (7.20) have been solved with both standard finite element and
proposed method. Different boundary conditions have been applied (top figures) to
solve the problems. Top two figures of Figure 7.7 show the maximum temperature
difference between standard finite element solution and proposed method solution is
0.33% with respect to maximum temperature(100°C'). The bottom two figures of
Figure 7.7 show if heat generation is considered (equation 7.20) and the boundary
temperature is 25°C. The maximum difference between the standard finite element

126

Finite Element Solution Combined Method Solution

 

 

 

 

 

 

 

 

 

Figure 7.5: Temperature distribution when all the boundaries with temperature T =
25° and heat generation, 9 = 10.0 for full finite element method (left side) and

combined method (right side).

and the proposed method is 0.46%. Therefore, it can be concluded that the proposed

method provide very good results.

7.2.2 Creeping flow

The governing equations of creeping flow in the steady state condition with-

out any body force are[39, 62]

 

(9213 8213
El; + a—y2 — (7.21)
(92V 82V 1 6P
1 .__x _ —— = .22
01:2 + (9342 )1 6r 0 (7 )
2 2
av, av,_1§_13_0 (7.23)

 

 

(9r? 6y2 11 By —
where 11 is the dynamic viscosity. Assuming 11 = 1, the above three equations (7.21-
7.23) have been solved (Figure 7.8) by both methods. Due to stability considerations,

127

       
 
 
 
    

    

    

as:-
Q ~
‘1‘-

 
 
     

2*

  

‘I
as.
::.‘.'ii
{E‘s-Ill

 

 

ti)
o.—
N."
A
a)

 

Figure 7.6: Grid of a complex geometry for full finite element method (left side) and

combined method (right side).

the pressure interpolation function is linear, whereas velocity approximations will be
quadratic[36, 47, 66]. Figure 7.8 shows the velocity profile and streamline for both
methods. Streamlines are straight lines along positive X direction as it should be
according to analytical solution of the same boundary conditions[39]. The maximum
differences by both the methods have been calculated in the computational regions.
The difference between both methods in pressure calculation is 0.36% with respect
to maximum pressure difference. The difference in velocity calculation in X -direction
V1. is 0.63%, and in Y—direction Vy is 0.34% with respect to theoretical maximum

velocity (28.8m/ sec).

Another flow field governed by the equations (7.21-7.23) is shown in Figure 7.9.
The boundary conditions have been altered to illustrate the generality of this method.
The difference in pressure by both methods is 0.36% w.r.to maximum pressure dif-
ference, maximum velocity difference in X direction V, is 0.63%, and in Y direction

128

No heat

       
 

   

All boundaries with Al boundaries with

Temperature 25.0 Temperature 25.0

Source term g=10.0 $002M“ 9=10.0
T

 

 

 

 

 

 

 

 

Figure 7.7: Temperature distribution without heat generation(top two figures) and
with heat generation(bottom two figures) by full finite element method(left side) and

combined method(right side)

129

Finite Element Solution Combined Method Solution

 

 

 

 

 

 

Figure 7.8: Solution for Poiseuille flow by both the methods. Vy = 0.0 at inlet and

outlet; Pressure P = 24.0 at inlet and P = —24.0 at outlet.

Finite Element Solution

 

 

 

 

Figure 7.9: Solution for fluid flow when Vy = 0.0 at inlet and outlet; pressure P = 24.0
at inlet and P = —24.0 at outlet, at top and bottom portion there is a suction of 5.0

unit.

130

Vy is 0.37% with respect to maximum velocity.

7.3 Summary

The error in the proposed method is higher than standard finite element calcula-
tions for test cases. The extra errors arise from interpolation criteria at the analytical
boundary. Furthermore, the floating error will be more because the maximum front
size in the frontal solver will be higher compared to standard finite element calcula-
tions. But there will not be an accumulation of the error from inside the analytical
domain. In fact for more complex geometries, the proposed method might provide
better results as there will not be errors developed by finite element method[44] inside

the analytical domain.

131

 

Case

Galerkin Finite Element

Combined Method

 

7.3

NN=627, NE(linear)=578,

RM=0.4917MB, RT:0.515unit,

ME:0.0001%.

NN2326, NE(linear)=248,

RM=0.2030MB, RT=0.468unit,

ME=0.077%.

 

7.4

N N: 1827, NE(quadratic)———578,

RM=2.9616MB, RT:1.672unit,

ME:0.0001%.

NN:896, N E(quadratic):248,

RM=1.2525MB, RT=0.906unit,

ME=0.037%.

 

NN=1827, N E(quadratic):578,

RM=2.6886MB, RT=1.515unit.

NN=896, N E(quadratic)=248,

RM=1.1062MB, RT=0.875unit.

 

7.7

NN=2469, N E(quadratic) :782,

RM=4.725MB and, RT=2.968unit

first

for case(g=0.0), and

RM=4.5528MB, and RT:2.859unit

for second case(g=10.0).

NN=1538, NE(quadratic)=452,

RM=2.26MB, and RT=1.484unit

for first case(g=0.0), and

RM=2.1519MB, and RT=1.343unit

for second case(g=10.0).

 

7.8

NN=1827, NE(quadratic)=—578,

RM=16.1717MB, RT=11.843unit.

NN=896, NE(quadratic)=248,

RM=6.6631MB, RT=7.438unit.

 

7.9

 

 

NN:1827, NE(quadratic)=578,

RM=16.1717MB, RT=12.016unit.

 

NN=896, NE(quadratic) 2248,

RM=6.6631MB, RT=7.656unit.

 

 

Table 7.1: Computational information for both the finite element and combined
method (NN = number of nodes, NE 2 number of elements, RM 2 required computa-
tional memory, RT = required time for solving, and ME 2 maximum error according

to analytical solution).

132

CHAPTER 8

SUMMARY AND CONCLUSIONS

An approach based on computing the average orientation tensor is used in this thesis
to predict flow-induced alignment of particles. The method is more economical than
a classical approach based on an orientation distribution function. The governing
equation for the second order tensor however requires a closure model for a fourth
order orientation tensor. Such closure model has been developed at Michigan State
University and it retains all the six symmetry and projection prOperties of the fourth
order tensor, unlike many others closure models. Different forms of this model have
been presented in this work. The coefficients (the 0,5) of this model have been
back computed with the help of exact solution of the governing equation for the
orientation distribution function. An irreducible form of the closure model has been
developed and discarded due to the significant changes in the coefficients near a
random orientation state which may lead to large errors. A sensitivity analysis of
the coefficients with respect to the invariants of the anisotropic component of the
orientation tensor shows that the coefficients are dependent on the invariants of the

orientation tensor.

Four coefficients arise in the most complete form of the closure model. These
coefficients have been found as functions of the invariants after a regression analysis.

However, solutions obtained with these coefficients do not compare well to the exact

133

solutions for simple prescribed flow fields. Thereafter only linear and quadratic terms
of in the fully symmetric model have been considered (the so—called FSQ model). The
sole coefficient of the FSQ model (C2) at the two corners of the realizability diagram
has for values 1 / 2 and 1 / 3. C2 = 0.37 was found to provided satisfactory results
for various flow fields. Furthermore with C2 = 0.37, many investigated flow fields
including fiber-fiber interactions are found to be realizable (non-negative eigenvalues)
and solutions of the evolution equation compare well to available exact solutions. A
model parameter C2 = 0.37 has been selected for the subsequent studies.

A finite element code has been developed for predicting the flow induced orienta-
tion of suspensions in complex geometries. This code allows for solving the continuity,
momentum, orientation, and stress equations in a fully coupled manner. Incorporat-
ing the fiber stress in the momentum equation causes the suspension to behave like
non-Newtonian fluid. The orientation and fiber stress equations are nonlinear and
a Picard iterative method has been adopted for obtaining the solution. Quadratic
approximations of the unknowns have been used in the simulations except for pres-
sure. The inclusion of the non-linear convective term in the momentum equations,
when the Reynolds number is greater than unity, is found to slow significantly the
numerical solution convergence rate.

Various two dimensional problems involving complex geometries have been solved
using the FSQ model with a C2 coefficient set at 0.37. These problems include
flow between eccentric rotating cylinders, a contraction flow problem, and flow in a
diverging nozzle.

It is found that the stress tensor must be considered in order to obtain a good

134

approximation of the pressure. In the diverging flow problem, the change in the
diverging angle results in different orientation states. The second eigenvalue of the
fiber stress tensor is much lower than the maximum eigenvalue indicating that the
fiber stress is prominent in the orientation direction. In the eccentric cylinder, the
quadratic and hybrid models have also been used so as to compare the performances
of these closure models with the FSQ closure. The fiber stress predicted by hybrid
model differs by 30% from the FSQ model. In case of a coarse mesh, the FSQ model
provided realizable results whereas the other two popular models showed unrealizable

results a certain region in the computational domain.

A parametric study of the flow parameters, which include the aspect ratio of the
fibers, concentration of the fibers in the suspension, and inlet velocity, has been per-
formed on the diverging flow problem to investigate their effects on flow induced fiber
orientation. An increase in the aspect ratio increases the fiber stress but does not
affect significantly the orientation tensor. This is because a change in aspect ratio
affects the coefficient of fiber stress whereas it marginally affect the dimensional para-
meter (/\) appearing in the orientation equation. Inclusion of the fiber interactions in
the governing equations tends to orient fibers randomly. The maximum eigenvalues
of the orientation tensor are reduced by almost 10% compared to the values with
negligible fiber interactions. The increase in concentration of the fibers also increases
the fiber stress without affecting significantly the orientation tensor. A change in the
magnitude of the velocity field at the inlet affects both the average orientation tensor

as well as the fiber stress tensor.

A new method has been developed to solve second order linear partial differential

135

equations. The method is based on using simple Green’s function inside the domain of
interest and combining the analytical solutions with the computational method. The
geometry for the Green’s function should be carefully chosen in order to result in a
simple expression. Modification of the governing equations may be needed to achieve
a well defined Green’s function. This new method has been applied to the solution of
heat transfer as well as fluid dynamics problems, which involve inter-related second
order partial differential equations. The information of calculating time and required
storage memory space in table provide real advantages of this combined method over
the finite element method. The reduction of the computational time depends on the
possible size of the analytical domains. Even though nonlinear equations have not
been done yet, it can be safely stated that the proposed method has the capability of

solving those problems in less time and memory space.

136

APPENDICES

A. Tensor notation used in the dissertation

Double-dot product

Double-dot product of two 3 x 3 tensors (A, B) is as follows

A I B = 2 21415333: allb11+ a21b12 + a31b13 +

1:1 3:1

0121721 + 022b22 + 032b23 + 0131131 + 023b32 + 033b33

Dyadic product

Dyadic product of two vectors (x, y) is as follows

1151291 3313/2 33193

3 3
x.y=E E x,yje,-ej= 1152311 3323/2 1723/3
121 j=1

$391 115392 3333/3

 

 

137

Dyadic product of two 3 x 3 tensors (A, B) is as follows

011511
011b21
alibai
021511
021b21
a21b31
a31b11
a31b21

(131531

 

b

B. Numerical model for distribution function

(111512
011522
011532
021012
0'21b22
021b32
031512
031b22

031532

S(AB)=Z

a31b23

a31b33

0

1:1 J
a12b11
0121921
0121931
0221311
0221921
a22b31
a32b11
0323221

0321731

3

1

E E aijbklefijeke) 2

i=1 J

=1

012512 012b13

6112522
012532
0221712
022b22
022532
0321312

032 522

012523
0121933
(1221913
a22b23
022533
a32b13

0321723

032 1’32 (132533

0131311
(113521
013531
(1231911
a23b21
023531
0331711
0333321

033b31

aisbiz
013522
013532
(1231312
023522
023532
(133512
0331322

033b32

013b13
0131323
013b33
0231913
0231123
023b33
0335313
033b23

033533

Finite difference formulation for Smoluchowski equation (equation 2.17)is

t t t t t
C 1.i.iiwi,ii + 02..-,.~.-w.-+1,.-.- + Cami/(1,161“ + C4.i.z‘i¢1_1,u + Cs,z‘.u¢1,1i—1 =

/t—l t—l t—I
mem + D2.1.iii/J1+1,a + Dayan/2.).“ + 04m“?

The entities of the matrix are as follows

02,1,1'1' :

C1,i,ii :

—

 

 

L.

138

t—l
i—1,ii

+ Dwarf/J

—((8¢)2(D38t + (86)?) + D38t(06)2(csc 0)2+

1.5(0¢)26t(89)2)\ (Vu) : er(0, ¢)er(0, 76))

t—l
Lit—1

1

 

d

(8¢)20t(0.25DRc90 cot 9 + 0.258601 + 1) (Vu) : er(6, ¢)eg(6, (b)

+0.51)R(a¢)2 — 0.25009 — 1) (Vu)T : e.(0, ¢)e9(9, a))

 

 

C (0¢)20t(—0.25D386 cot 6 + 0.258601 + 1) (V11) : e,(6, ¢)eg(6’, cf)
3,1",11‘ :
+0.5DR(0¢5)28t + 0.2539(A — 1) (Vu)T : e,.(6, ¢)e0(9, a))
' 1
C (636)20t csc 6(0.5DR csc6 — 0.258¢(A + 1) (V11) : er(6, ¢)e¢((9, (0)
4.1.21 =
+0.256¢(/\ — 1) (Vu)T : e.(9, ¢)e,,(a, (1))
C (86)20t csc 6(0.5DR cscO + 0.258¢(A + 1) (V11) : er(9, ¢)e¢(0, QB)
5,1,11‘ =
41.250119 — 1) (V11)T : e.(a, ¢)e,(9, a))

 

 

D _ —((8¢)2(DRat — (86?) + DR8t(66)2(CSC 6V—

1.5(aa)2at(ae)2A (W) : e40, ¢)er(6. ¢))
D2,i,iz‘ = ’C2a'.ii

D3,i,ii = —C3,i,ii
D4,i,ii = _C4,i,ii

Dam = —C5,i,ii

C. Parametric approximation by shape functions

Any value of a parameter inside the realizable diagram can be computed as
quadratic approximation of six nodal values of the parameter. Assuming the FSQ
parameter C2 is a quadratic function in :c and y (III and II indicates as 113-axis and

y-axis respectively), C2 can be written as
C2 = (11 + 022: + 03y + any + a5x2 + (163/2 (0.2)

N ow the six nodal points considered quadratic approximation are {{0.0, 0.0}, {0.2, 0.6},
{—0.0278, 0.167}, {0.1, 0.39}, {0.114, 0.4482}, {—0.0139, 0.1050}} (Figure 8.1). The

139

magnitude of C; at any point inside the domain can be defined as,

C; = Adi (CB)

where C; is the value of C2 at the 1"" node and A is calculated for the above six points

as
1 0 o o 0 o
1 0.2 0.6 0.12 0.04 0.36
1 —0.0278 0.167 —o.00564 0.000773 0.02789
A: (04)
1 0.1 0.39 0.039 0.01 0.1521

1 0.114 0.44822 0.0511 0.013 0.201

 

 

1 —0.0139 0.105049 0.00146 0.000193 0.01103

i—

Again using Lagrange shape functions, the approximation can be done by
02 (3.9) = N.- (23.31) Cf (05)

where N,- (r, y) are the shape functions corresponding to those six points. Now using
the above four equations (equations C.2-C.5)[46], the isoparametric shape functions

can be calculated as function of :r and y.

Several values of C2 are back computed as well as the corresponding invari-
ants (section 4.3) are calculated for simple shear flow without fiber interactions.
Least square approximation has been done to achieve the optimal values of 02 at
the six nodal points. They are {—4.76738, 0.269573, 0.0693274, 0.275868, 0.398572,
0.830628} at the six nodal points in Figure 8.1 respectively.

140

 

 

Figure 8.1: Isoparametric quadratic shape function approximation has been taken for

the model parameter in FSQ C2.

141

D. Stress in a suspension of force-free particles

The particle rotates in a suspension when it experiences a couple by external
means, immersed in a Newtonian fluid[67, 68]. The fluid will experience some stress
due to movement of the particle and vice-versa. The stress can be expressed in
terms of instantaneous particle orientations. When the suspension is statistically
homogeneous inside a small control volume, the bulk stress and bulk velocity gradients
in the suspension are defined as average of an ensemble of realization; these averages
being equal to the integrals over a suitably chosen volume of ambient fluid and particle
together. The following is extracted from Batchelor[40].

The volume and surface of a typical particle in the suspension volume V are

denoted by W) and A0. Then the bulk stress in the suspension is

1
EU 2 —/ (CU — puiu’) dV (13.6)
V v 3
where u’ = u — 11 where u’ is the fluctuation of the mean velocity (ii). Since the

ambient fluid is considered to be Newtonian, with viscosity )1, the bulk stress can be

written[40] as

 

 

1 811- 311-
21 = — — 6. ‘ ’ dV
J V v—zvol p J+#(6$j+axi)} +
1 1 ,
_ 2.. _ _ . ’.dV D.7
V2: (fvoajdv) V./V]: 11,11] ( )

where the summation is over all the particles in V. The volume average velocity

gradient (OM/8:13,) can be written as

CUz' 1 Bu,-
= — — v
8113.1 V A GIL‘jd

 

 

1 011,- 1
: — dV + V Z to uinJ-dA (13.8)

V V—EVo a$j

142

where n,- is the outward unit normal to surface A0. Further the following relation

exists while using divergence theorem

faijdV2/ aikrjnk— g—afixjdV (D.9)
V0 A0 V0 xk

As Bait/(9r), is divergence of stress, flaw/8:13), term will be equal to the body force
per unit volume resulting

if“ = pf: (1110)
55k

 

f,-’ is local acceleration relative to the average acceleration in V and p is assumed to
be uniform throughout the suspension. With the aid of equations (D8), (D9) and

(D.10), the equation (D.7) can be written as

 

8U BH-
2, = —6,-- V ' ——’- (7 DH
.7 ]./V—)3Vopd +#(a$j + 0.731;) +01] ( )

where
1
of, = V E {aikrjnk — M(Uinj +u,-n,~)}dA —
A0
1 / , 1 7 I
— pfix-dV — —/ puiudV (D.12)
v2 .0 . v v .

Assuming the effect of inertia forces associated with the fluctuations about the average
motion negligible, the last two terms of equation (D.12) can be ignored. Further, the

angular momentum of one equation for one particle is

Li Z Eijk/ 0:)de (D.13)
V0

However, the angular momentum depends on the antisymrnetrical part of the bulk

stress. Therefore, equation (D.12) can be modified as

1 1 1
0f]- : ? 2L0 {507,115+ ajkxi)nk — [.1 (mm + 111-7711)} (114 "l" W ZeijkLk (D.14)

143

It can be noticed that, since no—slip condition, the integral of (Uinj + ujni) over
the surface A0 is zero for either translational or rotational motion (or a combining
therefore) of a rigid particle.

The particle inside the suspension causes some disturbances [69]; and the velocity,

pressure, and stress in the fluid can be written as

 

 

U5 2 23%? +712, P=P0 +17,
’ (D.15)
O’ij = —p0(5,:j + 2pm.) + 0'2]-
According to Lamb (1932), the disturbance velocity can be written as
, 3$i$jxk xjéik — $1,50-
u, = DJ)c (— 2T5 + r3 + - ~ (D.16)

where D3), are coefficients, which are determined by the inner boundary condition,
and the condition of zero volume flux across the surface enclosing the particle requires

Di,- 2 0. Therefore, the disturbance stress follows as

Bu.- 6117' )

 

 

6x,- + 8x.-

3$k(xi(5jl + CCjCSu) 15xiccjxkx)

Uij = —P’51‘j + #(

 

([117)

Since equation (D.14) provides the same value when taken over any closed surface
surrounding the particle, it is assumed the center of the local coordinate is same as
the center of the particle. Again Jeffery (1922) shows when a rigid ellipsoid particle
on which a couple (no force) is imposed by external means, the disturbance velocity
inside the suspension depends on the tensor Djk. With unit vectors parallel to the
principle axes of ellipsoid denoted by p, q, r, the relation between the components of

the tensor D,,- for the general and special sets of axes is

Dij = Dkl’hfljz (D18)

144

where ”7’11 = ”711 = pi, ”712 = ’72:“ = (12', and 713 = 731' 2 Ti-

f

After some calculations, 0,]. can be written as[40]

(722- = k Aijkl Ski (D.19)

where k 2% (6) is a factor depending on the fiber length and fiber concentra-

tion [38] and 773 is the Newtonian viscosity of the fluid, vL3 is the volume fraction

of the particle, (1,, = L/d is the aspect ratio of the particle, f (e) = 1713;, and
e = [ln(2L/d)]_1. The concentration of fibers in the suspension
6 = 57121.0 /4 = met/(411;) (D.20)

Therefore, the factor k largely depend on the concentration of the fibers.
To include the diffusion term in the fiber stress equation the fiber stress equation
will be[38]

0';- = k Aijkl Ski + 2kDR(a,-j — 3013') (13.21)

Stress prediction by the proposed closure

The stress model given in equation (D.21) largely depends on the closure model
of fourth order orientation tensor (AU-k1). The experimental results have been re-
ported for fiber suspension by Zirnsak[70]. The experiment was done in Newtonian
solvents, namely glucose wheat syrup and epoxy resin. The first normal stress dif-
ference [70, 71] data (N1) in Figure 8.2 and 8.3) are normalized to a test theory
of Carter (1967), who proposed fiber-fiber collisions occur, yielding a prediction
N1 oc Wmag/ (ln(2ap) — 1.8). Taking n = 0.6Nsm‘2 for this solution, the results
are compared with the proposed FSQ model as well as popular quadratic model. The

145

Aspect ratio = 276

 

 

 

 

100 -.
0 EXPERIMENTAL

10 4
— FSQ CLOSURE

1 « —QUADRATIC
CLOSURE
0.1 4
0.01 . . . .
0.1 1 10 100 1000

Shear rate

Figure 8.2: The normalized first normal stress difference Nl/ [rtnag/ (ln(2ap) — 1.8)]
versus shear rate of the glass fibers of aspect ratios (1p 2 276 in Newtonian solvents,
glucose wheat. syrup, and epoxy resin. The volume fraction (0 is in the range for
0.0002 — 0.0011. FSQ closure model (C2 = 0.37) provides very good results whereas

the quadratic model falls far short.

146

 

Aspect ratio = 552

 

 

100 5
. EXPERIMENTAL
1o 1 — FSQ CLOSURE
—QUADRATIC
1 . 8 CLOSURE
0.1 ]
0.01 . . .
0.1 10 100 1000

Figure 8.3: The normalized first normal stress difference N1/ [qfinag/ (ln(2ap) — 1.8)]
versus shear rate of the glass fibers of aspect ratios 0,, = 552 in Newtonian solvents,

glucose wheat syrup and epoxy resin. The volume fraction 05 is in the range for

0.0002 — 0.0011. FSQ closure model (Cg = 0.37) provides very good results whereas

Shear rate

the quadratic model falls far short.

results show that the FSQ closure provides very good result whereas the quadratic
model does not provide the results, that are close to experimental results. In fact,
it provides almost same result with increase in shear rate. Therefore, it can be con-

cluded that FSQ model will provide very accurate results as it can well predict the

interaction forces between fiber and suspended fluid.

147

E. Upwinding

To achieve better and accurate result the convective term in the momentum (equa—
tion 5.7) is included. With Reynolds number greater than 1, generally upwinding
technique is used to lower the damping effect of the numerical calculations[43]. In
this case, streamwise upwinding Petrov—Glerkin (SUPG) technique has been consid-
ered to solve for the problems[72]. The brief theory for this finite element upwinding

is described below.

The steady state convective term of a parameter can be written as

—VTDV¢ + V-ng = 0 (E22)
g} 1 0
where 06 is a functional parameter, V Xy E , D E , and V = velocity
a
{9—}, 0 1
u
vector 2
v

 

 

Considering the Figure 8.4, it is possible to introduce a new coordinate system
(s, t) of the equation (E22), such that velocity of the element (V) is co-linear with
the direction of axis S, eventually perpendicular to axis t. Thus the equation (E22)

can be rewritten as

—V3.DV..¢ + IVI V754 = 0 (13.23)

where |V| is the magnitude of the velocity vector. The gradient operator V is related

between the two coordinate system by

Vs: = TVXY (E24)

148

 

 

 

 

 

 

 

>X

Figure 8.4: The velocity vector(V) is in XY plane. V is aligned with the axis -s for

streamwise upwinding.

149

cos 6 sin 6
where T = . The equation (E23) can be in one equation because

 

 

—— sin 6 cos 6

there is no Velocity in the t direction. The equation can be written as

0 605 _
‘0? (D WI 2:) + |V| g _ 0 (E25)

A balancing diffusion should be added only in the direction of the flow to reduce

the damping effect for _the convective term. An anisotropic balancing diffusion of the
form D’ = 2 is added in the diffusion term D; where h is the characteristic
0 0

length, whiEh is determined by the cross sectional length of the element along 5

 

 

direction. Hence the equation (E25) has been modified to

<9 , 375 645 _
7); ((D + D) |V| 5;) + WI 5; _ 0 (E26)

N ow the resulting equation (E26) has been converted back to X Y coordinate system.

Therefore, the modified equation of (E22) after upwinding will be

(9 86) au,h 8¢ .2? _
Bag. [138% + 2lV| 2:115] +0.62%. _ 0 (E27)

 

 

 

The weak formulation of the equation (E27) can be obtained after rearranging as

follows

 

666w," ..
/D(5—, ax.)dQ+/..“‘ (W +21V12u’671d")

=/ DW“8¢dI‘ (E28)
avdr

The same upwinding procedure has been incorporated for the momentum equation
to reduce the effect of diffusion error.

150

F. Matrix entities in the finite element formulation

The global matrix of the finite element formulation is

 

. . - - F , - . -
All!) 0 0 0 (EU K11 K12 0 K14 xv f” ]
0 0 O 0 :13” K21 0 O 0 (13p fp
d _
r; + _
0 0 M,8 0 1:8 K3. 0 K33 0 22’ f’
0 0 0 0 2:7 K41 0 K43 K14 1” f7

 

 

 

 

 

 

 

 

 

The entities of the matrix are as follows

For momentum equations:-
m1111 0
A’IIU = where 777.1111 2 Inc quWudQ

0 7711111

killl k1112

 

 

 

 

K11 = where 191111 = fa, 11(66‘? 35’: + 33;" 5&3) cl!) and
klll2 kllll
_ IS U PG 61V“ 3W“
(131112 — file pl/L (116—x + U?) d9
k1211 61V“
K12 = where ((71211 = ch 62:.- WPdQ
k1211
k141i 0 (C1413
[C14 2 where [01413 = — fa: Wuaa%dfl

0 (C1413 131411

For continuity equation:-

K21:

 

[C2111 [C2111 3:1:

where ($72111 = [C1211 = ff) 6“,.“ WPdQ
For orientation equation:-

m3111 0
Mi? = where 7721111: If) WuWudQ

0 m3111

151

/ k3111 k3112 \

 

k3111 k3112 where 163111 = er(—l + C2 + 2C2a§3)W“ag:u (if! and

 

 

\k3lll k3112 )
1:31” : er(_1 + 02 + 202a§3)W“%3-d9
( k3311 [C3312 k3313 \

K33 = where

[C3321 k3322 k332i}

 

 

\k3331 k333i! £53311)

3((—(1+6(32)i‘é+ -j—”+Cg(15a12(j—“- y+d")
k3311= er 7 a” “ Wuwuda,

+(—1 + 7022)j—:))A) + 3—101102(51% —- 8%)

k f %(—d '§+ 02(5012G—Z + 3:) + %))’\+ Wuwudfl
3312 = n, ,
W012C2( _4dv))‘

35

;(_ 35(— —7+/\+6C'2)\)— —(7+A+6C'2 ))+

602A)» + 335012C2(9% _ 83—3»

k33l3 : fag WuWudfla
328:;‘501202( — — Qd—yv)A
k f %(—602(j—: + 2:)6112 +— (g(l + 02(— 1+ Wuwucm
3321 = 9e a
7022))A) -— %a1102(4% — 3%)A
k f %((—1 + C2)d_u +_ d_v + 6023;; —1502(j—;+ wuwudo
3322 = 96 ,
3%)a1g))\ — 32302202(8%;i — 51%)A
l((—@(—7 + A + 602A) — @(7 + A+
[$3323 = er ( 7 dy dz WuWudQ,

1—14((12C'2(5% - 2%)(112A + g("7‘l‘ \

[93331 = fee (—3 — 402 + 2802a22)A) + j—ZW + (—3— Wuwudfl,

 

 

K 402 + 2802022))\))) — 3430.1102(%5 + %)A /

152

( {3((1202(—2gg + 5%)a12A + $7 — 3A— \

k3332 = In, 402A) — j—;(7 + (3 + 4CQ)A)))— W“W“d9,

 

 

\ 773602202(%5 + %)A )
k3333— — er (17(2801202(% + %)A) Wul’VudQ

For fiber stress equation:-

{ k4111 k4112 \

K41 = k where 1:41“ = ch(—1 + 02 + 202ag3)wuig—“dn and

194111 [94112

 

 

\ k4111 [C4112 )
k4112 = fne(_1+ CZ + 202033)WuQ%yI—udfl

{164311 [$1312 k4313\‘

m) vL3

K43 : k ‘1161 (2L/d)

where k— — f ( ) is a factor depending on

k432i M322 194323

 

 

K k4331 194332 194311
the fiber length and fiber concentration [38] and 773 is the Newtonian viscosity of the

fluid, vL3 is the volume fraction of the particle, L / d is the aspect ratio of the particle,
f(5) = ,—_}—5—, and 5 = [ln(2L/d)]"1.

%(6(—1+C2)——1502(%-: +ZZ>212+d—Z(

194311 = fag Wuwucm,
1+ Cg — 7C2022)) + 3—ga11C'2(—51%;-‘ + 8%)
l((—1+C'2)fl +6C2(-d*y-+ dv)012)+
7 d d dx

k4312 = fag y y WuWudQ,

'2'5'(1::C'2(m311i + 45-1-5“)
13,313 = fQ 1+ 02)( +g)+ Ha1202(8“ 433)) W’“W"dQ,

702a22))+ +315a1102(4% — 3%)

%<(—1+02)3—: —3( 2(—1+C2)%;’-)+

wuwudn,

(§(
7
1(602(:-—: —)a12 +d—' d:1:_( 1+ Cz—
m = In ( 7 “‘1 Wuwudn,
k4322 = er (

5CQ(% +da: fl)(l12)+ +3—15a2202(8% — 513—3)

153

33_5((_1 + C2)(Z—: +d—:)+
k4323 = foe d Wuwudo,

3—1501202(9:—:' — 8%)

 

 

{ %((—% + 02%;“ — % + 02% — 1502%a12+ \‘
[C4331 = foe 602%012 — 702%022 — 702%022» WuWudQ’
\ +£01102(:—:' ‘l‘ 332:) )
l(“—1 + C2lég+ +( 1+ Czlg‘d'
k4332 = ch 7 dy (11+ Wuwudo, and
302((2% — 5%)a12)) + 3502202(Cdl—: + %)
k4333= f9 (1—3058C' Z—Zd-fi» Wuwudn

( [C4411 0 0

K44 = 0 [$4411 0 where [£4411 2 fm Wuwudo

 

 

\ 0 0 k4411 )

154

 

BIBLIOGRAPHY

[1] Imhoff A., “Validation of a closure model for flow-induced alignment of fibers”,
Diplomarbeit in Michigan State University (Mechanical Engineering), 2000.

[2] Bird R.B., Armstrong R.C., Hassager 0., “Dynamics of polymeric liquids”, Vol-
2, Fluid Mechanics second edition, 1987.

[3] Advani S.G., Tucker C.L, “The use of tensors to describe and predict fiber ori-
entation in short fiber composites”, Journal of Rheology, Vol-31, ppz751-784,
1987.

[4] Doi M., Edwards SE, “The theory of polymeric dynamics”, Oxford University
Press, 1998.

[5] Rao B. N., “Rheology and flow of fiber suspensions in composite material man-
ufacturing”, Dissertation of Mechanical Engineering Department, at University
of Oklahoma, 1996.

[6] Bretherton F. P., “The motion of rigid particles in a shear flow at low reynolds
number”, Journal of Fluid Mechanics, Vol-14, pp:284-304, 1962.

[7] Jeffery GB, “The motion of ellipsoidal particles immersed in a viscous fluid”,
Proceedings of the Royal Society of London. Vol-102 (Nov 1922).

[8] Ericsson K.A., Toll S., Manson J.-A. E., “The two-way interaction between
anisotropic flow and fiber orientation in squeeze flow”, Journal of Rheology,
Vol-41(3), pp2491-511, 1997.

[9] Lipscomb II G.G., Denn M.M., Hur D.U., Bogar D.V., “The flow of fiber sus-
pensions in complex geometries”, Journal of Non-Newtonian Fluid Mechanics,
Vol-26, pp:297—325, 1988.

[10] Bird, R.B., RC. Armstrong, and O. Hassager, “Dynamics of polymeric liquids”,
Volume 1, Fluid Mechanics, Second Edition, Wiley-Interscience, 1987.

[11] Greenberg M.D, “Advanced engineering mathematics”, Prentice Hall, NewJer-
sey, 1998.

[12] "flicker C.L., “Flow regimes for fiber suspensions in narrow gaps”, Journal of
Non-Newtonian Fluid Mechanics, Vol-39(3), pp:239—268, 1991.

155

 

[13] Cintra J .S., Tucker C.L., “Orthotropic closure approximations for flow-induced
fiber orientation”, Journal of Rheology, Vol—39(6) pp:1095-1122, 1995.

[14] Feng J ., Chaubal C. V., Leal L.G., “A comparison of closure approximation for
the Doi theory of LCPS”, Journal of Rheology, Vol-39(1) pp273-103, 1995.

[15] Folgar F., Tucker C.L., “Orientation behavior of fibers in concentrated suspen-
sions”, Journal of Reinforced Plastic Composites, Vol—3, pp:98-119, 1984.

[16] Dinh. S.M, Armstrong R.C, “A rheological equation of state for semiconcentrated
fiber suspensions”, Journal of Rheology, Vol-28, pp:207—227,1984.

[17] Dinh, S.M., Armstrong R.C., “A rheological equation of state for semiconcen-
trated fiber suspensions”, Journal of Rheology, Vol-28(3), pp:207-227, 1984.

[18] Bay R.S., Tucker C.L., “Fiber orientation in simple injection moldings. Part I:
Theory and numerical methods”, Polymer Composites, Vol-13(4), pp:317—331,
1992.

[19] Nguyen, C. T., “Molecular orientation of liquid crystalline polymers under simple
shear flows”, MS Thesis, Michigan State University, 2001.

[20] Mase G.T, Mase G.E.,“Continuum mechanics for engineers”, CRC Press, 1999.
[21] Carlsson L. A., Gillespie J .W.,“Micro-models for composite materials parti-

cles and discontinuous fibers”, Delaware composites design encyclopedia, Vol-2,
pp:91-142, 1989.

[22] Cullen C.G, “Matrices and linear transformations, second edition”, Dover Pub-
lications, New York, 1990.

[23] Parks SM, Weispfennig K, Petty CA, “An algebraic preclosure theory for the
Reynolds stress”, Physics of Fluids, Vol-10(3) pp:645-653 , 1998.

[24] Altan M.C, Advani S. G, Gueceri S. I, Pipes R.B, “On the description of the ori— -

entation state for fiber suspensions in homogeneous flows”, Journal of Rheology,
Vol-33, pp:1129-1155, 1989.

[25] Tannehill C., Andersor D., Pletcher R., “Computational fluid mechanics and
hear transfer”, second edition, Hemisphere publishing corporation, 1997.

[26] Stoer.J., Bulirsch.R., “An introduction to numerical analysis”, second edition,
Springer-Verlag, New York, 1993.

156

 

[27] Hand G.L., “A Theory of anisotropic fluids”, Journal of Fluid Mechanics, Vol-13,
pp233-46, 1962.

[28] Advani S.G., Tucker 111 CL, “Closure approximations for three—dimensional
structure tensors”, Journal of Rheology, Vol-34(3), ppz367-386, 1990.

[29] Advani S.G., Tucker, C.L., “Closure approximations for three-dimensional struc-
ture tensors”, Journal of Rheology, Vol-34, ppz367-386, 1990.

[30] Leal L.G., Chaubal C.V., “A closure approximation for liquid—crystalline polymer
models based on parametric density approximation”, Journal of Rheology, Vol-
42, pp:177-207, 1998.

[31] Petty C. A., Parks S. M., Shao S. M., “Flow induced alignment of fibers”,
in Proceedings of 12th International Conference on Composite Materials,
ICCM12/TCA, Paris, July 5-9, 1999.

[32] Parks S. M., Petty C. A., “Prediction of low-order orientation statistics for flow-
induced alignment of fibers and platelets, paper presentation”, Fundamental Re-
search in Fluid Mechanics: Particulate and Multiphase Flow I, AIChE Annual
Meeting, Dallas, TX, October 31- November 5, 1999.

[33] Kim Y. C., Nguyen C. T., Parks S.M., Petty C. A., Mandal D. and Benard A.,
“Rheological and alignment properties of liquid crystalline polymers”, Poster
Session on Fluid Mechanics, 2001 Annual AIChE Meeting, November 4—9, Reno
Hilton, Reno, NV, 2001.

[34] Mandal D., Parks 8., Petty C. and Benard A., “Discussion of a closure model
for flow-induced alignment of fibers”, Polymer processing society, 17th annual
meeting, May, 2001.

[35] Reddy J .N., “An introduction to the finite element method”, McCraw-Hill Pub-
lishing Company: New York, 1984.

[36] Ciarlet PC. and Lions J .L., “Handbook of numerical analysis, finite element
methods (Part-1)”, Vol-2, Elsevier Science Publisher, 1991.

[37] Bibbo M. A., Dinh S. M., Armstrong R. C., “Shear flow properties of semicon-
centrated fiber suspensions”, Journal of Rheology, Vol-29(6), pp:905-929, 1985.

[38] Larson R.G., “The structure and rheology of complex fluids”, Oxford University
Press, 1999.

157

[39] Munson B.R, Young D. F, Okiishi T.H, “Fundamentals of fluid mechanics”, John
Wiley & Sons, New York, 1998.

[40] Batchelor G. K., “The stress system in a suspension of force-free particles”,
Journal of fluid Mechanics, Vol—41, ppz545-570, 1970.

[41] Rajagopalan D., Robert C. Armstrong, Robert A. Brown, “Finite element meth-
ods for calculation of steady, viscoelastic flow using constitutive equations with

a newtonian viscosity”, Journal of Non-Newtonian Fluid Mechanics, Vol-36,
pp:159—192, 1990.

[42] Paul J .N., Robert C., Armstrong, Brown. R.A., “Finite element calculations
of time-dependent two dimensional viscoelastic flow with the explicitly elliptic

momentum equation formulation” , Journal of Non-Newtonian Fluid mechanics,
Vol-36, pp:109-133, 1990.

[43] Reddy. J .N ., Gartling D.K., “The finite element method in heat transfer and
fluid mechanics”, CRC Press, 2001.

[44] Bathe K J ., “Finite element procedure in engineering analysis”, Prentice—Hall:
New Jersey, 1982.

[45] Servais C., Luciani A., Manson J. E., “Fiber-fiber interaction in concentrated
suspensions: Dispersed fiber bundles”, Journal of Rheology, Vol-43(4), pp:1005-
1018, 1999.

[46] Hughes T.J.R., “The finite element methods”, Prentice—Hall: New Jersey, 1987.

[47] Brezzi F, Bathe KJ, “A discourse on the stability conditions for mixed finite ele—
ment formulations”, Computer Methods in Applied Mechanics and Engineering,
Vol-82, pp:27-57,1990.

[48] Khomami B., Talwar K.K., “Higher order finite element techniques for vis-
coelastic flow problems with change of type and singulatiries”, Journal of Non-
newtonian Fluid Mechanics, Vol-59 (1), pp:49—72, 1995.

[49] Khomami B., Talwar K.K., Ganpule H.K., “A note on selection of spaces in
computation of visco—elastic flows using the hp—finite element method”, Journal
of Non-newtonian Fluid Mechanics, Vol-52 (3): , pp:293—307, 1994.

[50] Chiba K. Nakamura K., “A numerical solution for the flow of dilute fiber sus—
pensions through an axisymmetric contraction”, Journal of Non-newtonian Fluid
Mechanics, Vol-35, pp:1-14, 1990.

158

[51] Kunji Chiba, Kazunori Yasuda, Kiyoji Nakamura, “Numerical solution of fiber
suspension flow through a parallel plate channel by coupling flow field with fiber

orientation distribution”, Journal of Non—Newtonian Fluid Mechanics, Vol-99,
pp:145—157, 2001.

[52] Yu C.-C., Heinrich J .C., “Petrov-Galerkin methods for the timedependent con-
vective transport equation”, International Journal for Numerical methods in En-
gineering, Vol-23, pp:883-901,1986.

[53] Hood. P, “Frontal solution program for unsymmetry matrices”, International
Journal for Numerical Methods in Engineering Vol—10, pp:379-399, 1976.

[54] Sloan. S. W., “A fortran program for profile and wavefront reduction”, Inter-
national Journal for Numerical Methods in Engineering Vol-28, pp:2651-2679,
1989.

[55] Sloan. S. W., “An algorithm for profile and wavefront reduction of sparse matri-
ces”, International Journal for Numerical Methods in Engineering Vol—23, pp:239-
251, 1986.

[56] Advani, S.G., Tucker III CL, “A numerical simulation of short fiber orientation
in compression molding”, Polymer Composites, Vol-11(3), pp:164—173, 1990.

[57] Advani S.G., “Flow and rheology in polymer composites manufacturing”, Vol-10,
Composite Materials Series, Elsevier, 1994.

[58] Mandal D. K., Benard A., “Numerical simulations of microstructured fluids using
a novel closure model”, WCCM V - Fifth World Congress on Computational
Mechanics, Vienna, Austria, July 7-12, 2002.

[59] Gartling D., “A test problem for outflow boundary conditions - flow over a
backword facing step”, International Journal for Numerical methods in Fluids,
Vol-11, pp:953—967,1990.

[60] F eng J ., Leal L.G., “Simulating complex flows of liquid-crystalline polymers using
the Doi theory”, Journal of Rheology, Vol-41(6), pp:1317—1335, 1997.

[61] Lawrence C. Evans, “Partial differential equations”, American Mathematical
Society: Rhode Island, pp:20-35, 1998.

[62] Gresho, Sani, Engleman, “Incompressible flow and the finite element method”,
John Wiley and Sons: New York, pp:367-375, 2000.

159

[63] Ozisik M. Necati, “Heat conduction”, John Wiley and Sons: New York, 1980.

[64] James H. Kane, “Boundary element analysis in engineering continuum mechan-
ics”, Prentice Hall: New Jersey, 1994.

[65] Raamachandran J ., “Boundary and finite elements”, Narosa Publishing House:
New York, 2000.

[66] Gunzburger M. D., “Finite element methods for viscous incompressible flows”,
Academic Press: San Diego, 1989.

[67] Phillip M. L., Brady J .F., “The hydrodynamic force on a rigid particle under-
going arbitrary time dependent motion at small Reynold number”, Journal of
Fluid Mechanics, Vol-256 pp:561—605, 1993.

[68] Phan-Thien N., “Constitutive equation for concentrated suspensions in New-
tonian liquids”, Journal of Rheology, Vol—39(4), pp:679-695, 1995.

[69] Stover C.A., Koch D.L., Cohen C., “Observations of fibre orientation in sim-
ple shear flow of semi-dilute suspensions”, Journal of Fluid Mechanics, Vol-238,
pp:277-296, 1992. ‘

[70] Zirnsak M. A, Hur D. U, Boger D. V, “Normal stresses in fibre suspension”,
Journal of Non-Newtonian Fluid Mechanics Vol-54, pp:153—193, 1994.

[71] Strand S.R., Kim S. Karrila S.J., “Computation of rheological properties of sus-
pensions of rigid rods: stress growth after inception of steady shear flow”, Journal
of Non-Newtonian Fluid Mechanics Vol-24, pp:311-329, 1987.

[72] Hughes T.J.R., Tezduyar T.E., “Finite element methods for first-order hyperbolic
systems with particular emphasis on the compressible Euler equation”, Computer
Methods in Applied Mechanics and Engineering Vol-45, pp:217—284, 1984.

160