University of Southern California Dissertations and Theses

Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures

(USC Thesis Other)

Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content HOMOGENIZATION PROCEDURES FOR THE CONSTITUTIVE MATERIAL
MODELING AND ANALYSIS OF APERIODIC MICRO-STRUCTURES
by
Preetham Aghalaya Manjunatha
A Thesis Presented to the
FACULTY OF THE USC VITERBI SCHOOL OF ENGINEERING
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulllment of the
Requirements for the Degree
MASTER OF SCIENCE
(CIVIL ENGINEERING)
December 2012
Copyright 2012 Preetham Aghalaya Manjunatha
To the Welfare of This Great World.
ii
Acknowledgments
I am highly thankful to Professor L. Carter Wellford, adviser, for his tireless patience
and guidance during the computational and theoretical work. Furthermore, his
immense support in the implementation of digital image based meshing techniques
which were instrumental in solving dicult material topologies. Consequently, Dr.
Wellford's patience has motivated and enlightened my perspective in research as,
\Patience is the key in any research and development!". In addition, I would like
to thank the thesis committee members, Professor Sami Masri, who encouraged
the idea of digital image processing and Professor Vincent Lee, who laid the strong
foundation in Solid Mechanics. Also, I am highly obliged to the thesis committee
members for their time and invaluable suggestions, which made this thesis successful
and on time.
Secondly, I am grateful to Shashank J. S., and Karthik T., MS students in
Image Processing, Electrical Engineering, Viterbi School of Engineering, for their
priceless suggestions in understanding the fundamentals of digital image processing
operations in a short period.
Finally, I am indebted to my parents, friends and roommates (Pandavas) for
their endless support during the dicult times.
iii
Table of Contents
Dedication ii
Acknowledgments iii
List of Tables viii
List of Figures x
Abstract xv
Chapter 1: Introduction 1
1.1 What is a Composite? . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Micro-structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.1 Types of Micro-structures . . . . . . . . . . . . . . . . . 5
1.2.2 A Discussion on Functionally Graded Materials
(FGMs) . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.2.1 Types of FG Micro-structures . . . . . . . . . . 7
1.2.2.2 Applications of FGMs . . . . . . . . . . . . . . 8
1.3 A Brief History on Homogenization Theory . . . . . . . . . . . . 9
1.4 A Review of the Literature . . . . . . . . . . . . . . . . . . . . . 10
1.4.1 Micro-Mechanical Analysis of Composites . . . . . . . . 11
1.4.1.1 The Voigt Approximation . . . . . . . . . . . . 11
1.4.1.2 The Reuss Approximation . . . . . . . . . . . . 12
1.4.2 A Self-Consistent Method . . . . . . . . . . . . . . . . . 12
1.4.3 The Mori-Tanaka Method . . . . . . . . . . . . . . . . . 13
1.4.4 The Method of Cells . . . . . . . . . . . . . . . . . . . . 13
1.4.5 Homogenization Methods . . . . . . . . . . . . . . . . . 14
1.5 Scope of This Thesis . . . . . . . . . . . . . . . . . . . . . . . . 18
Chapter 2: Homogenization Strategies for Heterogeneous Materials 21
2.1 A Mathematical Review of the Asymptotic Homogenization
Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.2 An Aperiodic Homogenization Method . . . . . . . . . . . . . . 35
iv
2.2.1 Assumptions and Geometry of a Sub-Patch . . . . . . . . 36
2.3 Determination of the Homogenized Elastic Constants . . . . . . 38
2.3.1 Aperiodic Elastic Constants . . . . . . . . . . . . . . . . 39
2.3.1.1 Computation of Anti-periodic Modes by an Iso-
parametric Method . . . . . . . . . . . . . . . . 48
2.4 Validation of the Aperiodic Elastic Constants . . . . . . . . . . 51
2.4.1 Case 1: One Dimensional Approach . . . . . . . . . . . . 51
2.4.2 Case 2: Two-Dimensional Approach . . . . . . . . . . . . 60
2.4.2.1 Case 2a: Homogeneous and Isotropic Medium
with Various Mesh . . . . . . . . . . . . . . . . 60
2.4.2.2 Case 2b: Checker Board Type Material Distri-
bution . . . . . . . . . . . . . . . . . . . . . . . 62
2.4.2.3 Case 2c: A Rectangular Hole in Homogeneous
and Isotropic Material . . . . . . . . . . . . . . 64
2.4.2.4 Case 2d: Soft and Hard Composite Material . . 65
Chapter 3: The Finite Element Formulation 67
3.1 The Three-Dimensional Variational Problem . . . . . . . . . . . 68
3.2 Formulation of The Multi-Scale Model . . . . . . . . . . . . . . 75
3.2.1 Discretization Procedure for the Multi-Scale Model . . . 75
3.2.1.1 Iso-Parametric Meshing . . . . . . . . . . . . . 77
3.2.2 Relationship of Micro and Macro Solutions . . . . . . . . 81
3.3 Solution of Micro Solution Variables . . . . . . . . . . . . . . . . 84
3.3.1 Anti-Periodic (Deformation) Modes . . . . . . . . . . . . 85
3.3.2 Computation of the Basic Axial Micro Solution Vari-
ables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
3.3.3 Computation of the Basic Shear Micro Solution Vari-
ables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
3.3.4 Signicance of Micro Solution Variables . . . . . . . . . . 97
3.4 Macro Equilibrium Equations . . . . . . . . . . . . . . . . . . . 98
Chapter 4: Micro Meshing Based on Digital Imaging Techniques 99
4.1 A Digital Image . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.2 Importance of the Discretization Procedures through Digital
Image Processing . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.2.1 A Typical Micro-structure Meshing by Imaging Tech-
niques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.3 Micro Meshing by Digital Image Processing . . . . . . . . . . . 104
4.3.1 Cropping and Scaling Processes in Micro Meshing . . . . 106
4.3.2 Discretization of Micro Element . . . . . . . . . . . . . . 109
v
Chapter 5: Numerical Simulations and Results 110
5.1 Geometry of the Two-Dimensional Cases . . . . . . . . . . . . . 111
5.1.1 Macro Domain . . . . . . . . . . . . . . . . . . . . . . . 111
5.1.2 Micro Domain . . . . . . . . . . . . . . . . . . . . . . . . 112
5.2 Case 1: Vertical Fibers Spacing Equally . . . . . . . . . . . . . . 114
5.2.1 Case 1a: Global and Local Deformation Analysis for
Axial Loading . . . . . . . . . . . . . . . . . . . . . . . . 117
5.2.2 Case 1b: Global and Local Deformation Analysis for
Shear Loading . . . . . . . . . . . . . . . . . . . . . . . . 119
5.3 Case 2: Functionally Graded Vertical Fibers in X Direction . . . 123
5.3.1 Case 2a: Global and Local Deformation Analysis for
Axial Loading . . . . . . . . . . . . . . . . . . . . . . . . 126
5.3.2 Case 2b: Global and Local Deformation Analysis for
Shear Loading . . . . . . . . . . . . . . . . . . . . . . . . 129
5.4 Case 3: Horizontal Fibers Spacing Equally . . . . . . . . . . . . 132
5.4.1 Case 3a: Global and Local Deformation Analysis for
Axial Loading . . . . . . . . . . . . . . . . . . . . . . . . 134
5.4.2 Case 3b: Global and Local Deformation Analysis for
Shear Loading . . . . . . . . . . . . . . . . . . . . . . . . 138
5.5 Case 4: Functionally Graded Horizontal Fibers in Y Direction . 141
5.5.1 Case 4: Global and Local Deformation Analysis . . . . . 143
5.6 Case 5: Periodic Fibers in Bi-direction . . . . . . . . . . . . . . 147
5.6.1 Case 5a: Global and Local Deformation Analysis for
Axial Loading . . . . . . . . . . . . . . . . . . . . . . . . 150
5.6.2 Case 5b: Global and Local Deformation Analysis for
Shear Loading . . . . . . . . . . . . . . . . . . . . . . . . 152
Chapter 6: Summary, Conclusions and Recommendations 154
6.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
6.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
6.3 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . 158
6.3.1 Type of Finite Elements . . . . . . . . . . . . . . . . . . 158
6.3.2 Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
6.3.3 Optimization . . . . . . . . . . . . . . . . . . . . . . . . 159
6.3.4 Types of Composites . . . . . . . . . . . . . . . . . . . . 159
6.3.5 Modeling of Interfaces . . . . . . . . . . . . . . . . . . . 159
6.3.6 Failure Analysis . . . . . . . . . . . . . . . . . . . . . . . 160
6.3.7 Non Destructive Testing (NDT) . . . . . . . . . . . . . . 160
Bibliography 161
vi
Appendices
Appendix A:Generalised Micro-Solution Variables for the Distorted
Finite Elements 166
A.1 Basic Axial Micro Solution Variable in X Direction . . . . . . . 169
A.2 Basic Axial Micro Solution Variable in XY Direction . . . . . . 172
A.3 Basic Axial Micro Solution Variable in Y Direction . . . . . . . 175
A.4 Basic Axial Micro Solution Variable in YX Direction . . . . . . 178
Appendix B: Computer Implementation of the Aperiodic Method 185
B.1 Flowchart of the Homogenization Program . . . . . . . . . . . . 186
B.2 Flowchart of the Digital Imaging Program . . . . . . . . . . . . 187
B.3 Flowchart of the Main Aperiodic Program . . . . . . . . . . . . 188
vii
List of Tables
Table 2.1: Unit Anti-Periodic Displacements . . . . . . . . . . . . . . 50
Table 2.2: Comparison of Homogenized Values in Homogeneous and
Isotropic Medium . . . . . . . . . . . . . . . . . . . . . . 62
Table 2.3: Comparison of Homogenized Values for a Checker Board
Pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Table 2.4: Comparison of Homogenized Values for a Rectangular
Hole . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
Table 2.5: Comparison of Homogenized Values for Soft and Hard
Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Table 5.1: Case 1: Horizontal Spacing of the Vertically Periodic
Fibers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
Table 5.2: Case 1: Summary of FEA Details . . . . . . . . . . . . . 116
Table 5.3: Case 2: Horizontal Spacing of the Functionally Graded
Vertical Fibers . . . . . . . . . . . . . . . . . . . . . . . . 125
Table 5.4: Case 2: Summary of FEA Details . . . . . . . . . . . . . 125
Table 5.5: Case 3: Vertical Spacing of the Horizontally Periodic
Fibers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Table 5.6: Case 3a: Summary of FEA Details . . . . . . . . . . . . . 134
Table 5.7: Case 3b: Summary of FEA Details . . . . . . . . . . . . . 138
viii
Table 5.8: Case 4: Horizontal Spacing of the Functionally Graded
Horizontal Fibers . . . . . . . . . . . . . . . . . . . . . . . 142
Table 5.9: Case 4: Summary of FEA Details . . . . . . . . . . . . . 143
Table 5.10: Case 5: Horizontal and Vertical Spacing of the Periodic
Fibers in Two Directions . . . . . . . . . . . . . . . . . . 149
Table 5.11: Case 5: Summary of FEA Details . . . . . . . . . . . . . 149
ix
List of Figures
Figure 1.1: Types of Common Composites . . . . . . . . . . . . . . . 4
Figure 1.2: Periodic Micro-structures . . . . . . . . . . . . . . . . . . 5
Figure 1.3: Aperiodic Micro-structures . . . . . . . . . . . . . . . . . 6
Figure 1.4: Types of Functionally Graded Micro-structures . . . . . 8
Figure 1.5: A Typical Representative Volume Element . . . . . . . . 14
Figure 1.6: A Typical Layered Functionally Graded Materials . . . . 17
Figure 1.7: A Micro-Macro Scale Coupling [34] . . . . . . . . . . . . 19
Figure 2.1: An Arbitrary 2-D Shaped Body . . . . . . . . . . . . . . 22
Figure 2.2: A detailed view of Micro-structures at a point `A' . . . . 22
Figure 2.3: A Typical Base or Unit Cell of a Composite . . . . . . . 23
Figure 2.4: Homogenization of a Composites in Two Scales . . . . . 35
Figure 2.5: A Typical sub-patch with Two Materials . . . . . . . . . 36
Figure 2.6: A Sub-patch with Graded Materials . . . . . . . . . . . . 39
Figure 2.7: A Distorted Iso-parametric Finite Element . . . . . . . . 49
Figure 2.8: A one-dimensional Macro-structure . . . . . . . . . . . . 51
Figure 2.9: A one-dimensional Micro-structures . . . . . . . . . . . . 51
x
Figure 2.10: A One-dimensional Fixed End Macro-structure . . . . . 58
Figure 2.11: Homogeneous and Isotropic Medium with Perfect Uni-
form mesh . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Figure 2.11: Homogeneous and Isotropic Medium with Highly Dis-
torted mesh . . . . . . . . . . . . . . . . . . . . . . . . . 61
Figure 2.12: Checker Board Material Distribution with Mesh . . . . . 62
Figure 2.13: Rate of Convergence of E
H
ijkl
vs. Mesh Size . . . . . . . . 63
Figure 2.14: A Rectangular Hole in Homogeneous and Isotropic Mate-
rial with mesh . . . . . . . . . . . . . . . . . . . . . . . . 64
Figure 2.15: Soft and Hard Composite Material with Mesh . . . . . . 65
Figure 3.1: Discretization of a Typical Sub-Patch . . . . . . . . . . . 75
Figure 3.2: Multi-Dimensional Micro Meshes . . . . . . . . . . . . . 77
Figure 3.3: Physical and Parent Plane . . . . . . . . . . . . . . . . . 78
Figure 3.4: An Example of Mapping Physical and Parent Plane . . . 78
Figure 3.5: Types of Finite Elements . . . . . . . . . . . . . . . . . . 79
Figure 3.6: Micro-Macro Relationship in a 4-Noded Element . . . . . 82
Figure 3.7: Types of Strain Fields . . . . . . . . . . . . . . . . . . . 84
Figure 3.8: Normal Deformation Modes of a Typical Sub-Patch . . . 85
Figure 3.8: Shear Deformation Modes of a Typical Sub-Patch . . . . 86
Figure 3.9: A One-Dimensional Rod with Four Finite Elements . . . 97
Figure 4.1: Representations of Images . . . . . . . . . . . . . . . . . 100
xi
Figure 4.2: A Cross Section of the Trabecular Bone Micro-structure
(Image Courtesy: Liebschner et. al. [21]) . . . . . . . . . 102
Figure 4.3: A Transition from Micro-graph to a Rened Mesh (Image
Courtesy: Langer et. al. [19]) . . . . . . . . . . . . . . . 103
Figure 4.3: A Transition from Micro-graph to a Rened Mesh (Image
Courtesy: Langer et. al. [19]) . . . . . . . . . . . . . . . 104
Figure 4.4: A Pictorial Overview of an Image Based Micro Mesh-
ing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
Figure 4.5: A Cropped View of a Micro Element . . . . . . . . . . . 106
Figure 4.6: A Typical Representation of Coordinate Systems . . . . 107
Figure 4.7: A Transition of Cropping to Meshing of a Micro Ele-
ment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Figure 5.1: A Typical Two-Dimensional Geometry of the Aperiodic
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Figure 5.2: A Typical Two-Dimensional Geometry with Micro/Sub-
Micro Elements . . . . . . . . . . . . . . . . . . . . . . . 112
Figure 5.3: Case 1: Geometric Details of the Vertical Fibers with
Equal Spacing . . . . . . . . . . . . . . . . . . . . . . . . 115
Figure 5.4: Case 1a: An Axial Deformation Plot . . . . . . . . . . . 117
Figure 5.5: Case 1a: Plot of a Global-Local X Displacement at the
Bottom Edge . . . . . . . . . . . . . . . . . . . . . . . . 118
Figure 5.6: Case 1b: A Shear Deformation Plot . . . . . . . . . . . . 119
Figure 5.7: Case 1b: Plot of a Global-Local X Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 120
Figure 5.8: Case 1b: Plot of a Global-Local Y Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 121
xii
Figure 5.9: Case 1b: Plot of a Global-Local X Displacement at Two
Locations (X-Intercepts) . . . . . . . . . . . . . . . . . . 122
Figure 5.10: Case 2: Geometric Details of the Functionally Graded
Vertical Fibers . . . . . . . . . . . . . . . . . . . . . . . . 124
Figure 5.11: Case 2a: An Axial Deformation Plot . . . . . . . . . . . 126
Figure 5.12: Case 2a: Plot of a Global-Local X Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 127
Figure 5.13: Case 2a: Plot of a Global-Local Y Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 128
Figure 5.14: Case 2b: A Shear Deformation Plot . . . . . . . . . . . . 129
Figure 5.15: Case 2a: Plot of a Global-Local X Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 130
Figure 5.16: Case 2a: Plot of a Global-Local Y Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 131
Figure 5.17: Case 3: Geometric Details of the Horizontal Fibers with
Constant Spacing . . . . . . . . . . . . . . . . . . . . . . 133
Figure 5.18: Case 3a: An Axial Deformation Plot . . . . . . . . . . . 135
Figure 5.19: Case 3a: Plot of a Global-Local X Displacement at the
Right Edge (Axial) . . . . . . . . . . . . . . . . . . . . . 136
Figure 5.20: Case 3a: Plot of a Global-Local Y Displacement at the
Right Edge (Axial) . . . . . . . . . . . . . . . . . . . . . 137
Figure 5.21: Case 3b: A Shear Deformation Plot . . . . . . . . . . . . 138
Figure 5.22: Case 3b: Plot of a Global-Local X Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 139
Figure 5.23: Case 3b: Plot of a Global-Local Y Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 140
xiii
Figure 5.24: Case 4: Geometric Details of the Functionally Graded
Horizontal Fibers . . . . . . . . . . . . . . . . . . . . . . 142
Figure 5.25: Case 4: An Axial Deformation Plot . . . . . . . . . . . . 143
Figure 5.26: Case 4: Plot of a Global-Local X Displacement at Two
Locations (X-Intercepts) . . . . . . . . . . . . . . . . . . 144
Figure 5.27: Case 4: Plot of a Global-Local X Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 145
Figure 5.28: Case 4: Plot of a Global-Local Y Displacement at Two
Locations (Y-Intercepts) . . . . . . . . . . . . . . . . . . 146
Figure 5.29: Case 5: Geometric Details of the Periodic Fibers in Two
Directions . . . . . . . . . . . . . . . . . . . . . . . . . . 148
Figure 5.30: Case 5a: An Axial Deformation Plot . . . . . . . . . . . 150
Figure 5.31: Case 5a: Plot of a Global-Local X Displacement at the
Bottom Edge (Axial) . . . . . . . . . . . . . . . . . . . . 151
Figure 5.32: Case 5b: A Shear Deformation Plot . . . . . . . . . . . . 152
Figure 5.33: Case 5b: Plot of a Global-Local X Displacement at the
Bottom Edge (Shear) . . . . . . . . . . . . . . . . . . . . 153
Figure A.1: A Rotated Iso-parametric Finite Element . . . . . . . . . 166
Figure B.1: A Typical Flowchart of the Aperiodic Homogenization
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
Figure B.2: A Typical Flowchart of the Digital Image Processing Tool
of Homogenization Method . . . . . . . . . . . . . . . . . 187
Figure B.3: A Typical Flowchart of the Main Homogenization FE
Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 188
xiv
Abstract
Composite materials are the well-known substitutes for traditional metals in various
industries because of their micro-structural character. Micro-structures provide
a high strength-to-weight ratio, which makes them suitable for manufacturing
large variety of applications ranging from simple toys to complicated space/aircraft
structures. Since, these materials are widely used in high performance structures,
their stress/thermal analysis issues are of major concern. Due to the high degree of
material heterogeneity, it is extremely dicult to analyze such structures.
Homogenization (rigorous averaging) is a process that overcomes the diculty
of modeling each micro-structure. It replaces an individual micro-structure by
an equivalent material model representation (unit cell). Periodic micro-structures
appear in regular intervals throughout the domain, in contrast aperiodic micro-
structures follows an irregular pattern. Further, this method bridges the analysis gap
between micro and macro domain of the structures. In this thesis, Homogenization
procedure based on anti-periodic displacement elds for aperiodic micro-structures
and aperiodic boundary conditions are considered to model the constitutive ma-
terial matrix. This work could be easily implemented with the traditional nite
element packages. In addition, it eventually increases the convergence accuracy and
reduces the high computational expenses. Dierent problems are analyzed by the
xv
implementation of digital image processing schemes for the extraction of a unit cell
around the Gauss quadrature points and the mesh-generation. In the future, this
research denes a new path for the analysis of any random heterogeneous materials
by its ease of implementation and the state-of-the-art micro-structure material
modeling capabilities and digital image based micro-meshing.
xvi
Chapter 1
Introduction
The heterogeneity is an important criterion for achieving high strength and thermally
ecient properties of any materials. As the heterogeneous materials vary constantly
throughout the domain of interest, it makes a design engineer to hurdle with
new complications and constraints for the optimal product development. For
instance, not all applications or optimal designs are associated with periodic
material distributions and boundary conditions. There are certain problems which
require an exponential distribution of ber and matrix volume fraction, random
distribution of whiskers/
akes in matrix and high volume fraction of matrix other
than ber. Apparently, the prior methods developed by researchers in the eld
of Material Sciences, Aerospace, Mechanical and Civil Engineering were limited
to analyzing the structural or thermal properties of periodic composites or by
considering the layer-wise averaged model. An absence of research work on the
analysis of the mechanical behavior of aperiodic micro-structures has motivated
this work to produce usable tools for industry.
This thesis concentrates on the development of material constitutive representa-
tions for engineering materials. The work focuses on the multi-scale approaches for
general micro-structures such as periodic, structured or graded, and aperiodic or
random structures. The rigorous averaging techniques, combined with the explicit
coupling of the micro-macro quantities like displacements, forces plays an important
role in dening the mechanical behavior of a domain. Nowadays, time constraints
1
are increasing with a regular pace, this work is designed to result in improvements
in the accuracy and eciency of modeling and analysis of composite materials. As
a result, it reduces the huge computational cost for analyzing composites materials
without considering large degrees of freedom (DOF). This proposed research work
consider a sub-patch around Gauss quadrature point to dene the homogenized elas-
tic constants. Thereby, reducing upto 90% of the computational time as compared
with the ne grid nite element analysis.
In the following sections the reader will be familiarized with the fundamental
terminologies such as the details about composites, micro-structures, functionally
graded materials and other common denitions. Also, a rich
avor of homogenization
technique will be introduced to the reader, so that it ensures a smooth transition
of concepts developed in the later stages.
1.1 What is a Composite?
\A composite is a structural material that consists of two or more combined
constituents that are combined at a macroscopic level and are not soluble in each
other", Kaw [3]. One constituent is called the reinforcing phase and the one in
which it is embedded is called the matrix. The reinforcing phase material may
be in the form of bers, particles, or
akes. The matrix phase materials are
generally continuous. Examples of composite systems include concrete reinforced
with steel and epoxy reinforced with graphite bers, etc. Composites materials
being lighter, stronger and tougher are used in various applications and have played
an increasing industrial role for the last 4 decades. Composite materials exhibit
the micro-structural properties by which they attains high load carrying capacity,
high yield strength and most importantly they have less weight in comparison to
2
other traditional metals. Composites are lighter because of the low density of bers
in matrices. Further, weight reduction provides the one of the most important
motivations for the use of composites in aerospace industries. Some of the common
classication relating to the composites are given below, Ever [7]:
Reinforcement
{ Continuous long bers
Unidirectional ber orientation
Bidirectional ber orientation
Random orientation
{ Discontinuous long bers
Random orientation
Preferential orientation
{ Particles and whiskers
Random orientation
Preferential orientation
Laminate congurations
{ Unidirectional Lamina
{ Laminate
{ Bulk Composites
Hybrid Structure
3
{ Dierent material in various laminae
{ Dierent reinforcement in a lamina
In composite materials, ber reinforcement is preferred because most materials are
much stronger in ber form than in their bulk form. This is because, there will be
a sharp deduction in the defects as the material size decreases. Further, the main
factors that drive the use of composites are weight reduction, corrosion resistance
and part-count reduction. Also, composite materials have an enhanced fatigue
life, thermal/acoustical insulation, low thermal expansion, low or high thermal
conductivity and others. Figure 1.1 in the following page shows a typical ber
reinforced uni-directional composite,
ake enriched composite and a composite
laminate.
(a) A Fiber Reinforced Com-
posite
(b) A Flake Reinforced Com-
posite
(c) A Fiber Reinforced Com-
posite Laminate
Figure 1.1: Types of Common Composites
1.2 Micro-structures
\Micro-structure is dened as the structure of a prepared surface or thin foil of
material as revealed by a microscope above 25x magnication", George [9]. The
4
material of a micro-structure can be broadly classied as either metallic, polymeric,
ceramic or composite. The material choice has a strong in
uence on physical
properties such as strength, toughness, ductility, hardness, corrosion resistance,
high/low temperature behavior, wear resistance, and others, which in turn govern
the application of these materials in industrial practice.
1.2.1 Types of Micro-structures
Micro-structures have a spectrum of variations within their physical denition
range. For simplicity, the micro-structures are classied under two categories. They
are:
Periodic Micro-structures
{ Periodic micro-structures have a spatial distribution of materials which
are symmetric and follow a regular pattern through out the domain.
Further, these micro-structures satises the sucient and necessary
periodic boundary conditions.
Matrix
Fiber
X
Y
Periodic in X - direction
Periodic in Y - direction
(a) A Typical Periodic Fiber Cross-
section
Matrix
Fiber
X
Y
Periodic in X - direction
Periodic in Y - direction
(b) A Typical Periodic Fiber Longitu-
dinal Cross-section
Figure 1.2: Periodic Micro-structures
5
Aperiodic/Random Micro-structures
{ Aperiodic micro-structures involve a spatial distribution of materials that
is unsymmetrical and follows an irregular pattern throughout the domain.
Further, these micro-structures fail to satisfy the sucient and necessary
periodic boundary conditions associated with classical homogenization.
Figure 1.3 shows the common types of aperiodic micro-structures.
Matrix
Fiber
X
Y
Aperiodic in X - direction
Aperiodic in Y - direction
(a) A Typical Aperiodic Fiber
Cross-section
Matrix
Fiber
X
Y
Aperiodic in X - direction
Aperiodic in Y - direction
(b) A Typical Aperiodic Fiber
Longitudinal Cross-section
(c) A Grained Longitu-
dinal Cross-section (Con-
crete)
Figure 1.3: Aperiodic Micro-structures
In this work, materials with periodic and aperiodic nature are considered to
analyze the strength property of structures. In particular, aperiodic materials, the
functionally graded materials (FGMs) with dierent kinds of spatial distributions
and boundary condition are analyzed. The next section discusses FGMs in detail
focusing on their properties and application.
6
1.2.2 A Discussion on Functionally Graded Materials
(FGMs)
The concept of functionally graded materials (FGMs) was proposed in 1984 by
materials scientists in the Sendai area as a means of preparing thermal barrier
materials. Continuous changes in the composition, micro-structure, porosity, etc.
of these materials results in gradients in such properties as mechanical strength and
thermal conductivity. In 1987, a national project was initiated entitled 'Research
on the Basic Technology for the Development of Functionally Gradient Materials
for Relaxation of Thermal Stress'. In 1992 when the project was nished, samples
of 300mm square shell and 50mm diameter hemispherical bowls for SiC-C FGM
nose cones were prepared, Koizumi [17].
\Functionally graded materials (FGMs) are a new generation of engineered
materials wherein the micro-structural details are spatially varied through non-
uniform distribution of the reinforcement phase(s), by using reinforcement with
dierent properties, sizes and shapes, as well as by interchanging the roles of
reinforcement and matrix phases in a continuous manner", Hirai [14]. The result is
a micro-structure that produces continuously or discretely changing thermal and
mechanical properties at the macroscopic or continuum scale.
1.2.2.1 Types of FG Micro-structures
The classication of functionally graded micro-structures are based on the spatial
variation of reinforcement through out the domain. In Figure 1.4a, the ceramic
and metallic inclusions are varied vertically and horizontally. Figure 1.4b shows
the micro-structural grading through non-uniform reinforcement spacing. Lastly,
7
in Figure 1.4c, dierent types of reinforcement are varied along vertical dimension
of the domain.
Ceramic matrix with
metallic inclusion
T
Hot
T
Cold
Ceramic phase
Metallic phase
Metallic matrix with
ceramic inclusion
Transition Region
(a) Continuously Graded Micro-structure
Metallic
Matrix
Ceramic
Inclusion
(b) Discretely Graded Micro-structure
Ceramic
Matrix
Inclusion
Phase 1
Inclusion
Phase 2
Inclusion
Phase n
(c) Multiphase Graded Micro-structure
Figure 1.4: Types of Functionally Graded Micro-structures
1.2.2.2 Applications of FGMs
Functionally graded materials are ideal candidates for applications involving se-
vere thermal gradients, ranging from thermal structures in advanced aircraft and
aerospace engines to computer circuit boards. In one such application as shown in
Figure 1.4a, a ceramic-rich region of a functionally graded composite is exposed to
hot temperature while a metallic-rich region is exposed to cold temperature, with a
gradual micro-structural transition in the direction of the temperature gradient,
Aboudi et.al. [2]. By adjusting the micro-structural transition appropriately, op-
timum temperature, deformation and stress distributions can be realized. This
8
concept successfully reduces the thermal fatigue resistance and life of ceramic
thermal barrier coatings.
Micro-structural grading through non-uniform reinforcement spacing, Figure
1.4b, or through the use of dierent types of reinforcement, Figure 1.4c, can also
be eectively used to reduce the mismatch in the thermo-mechanical properties
between dierently oriented, adjacent plies in a laminated plate, Aboudi et.al. [2].
Thus, the reduction of thermally induced inter-laminar stresses at the free edge of
a laminate (which result from the large property mismatch between the adjacent
plies) can be realized by using the functional grading concept to smooth out the
transition between dissimilar plies.
Further, applications are in optoelectronic devices, such as p-i-n diodes, hetero-
junction photo-detectors and lasers, Wosko et al. [33]. In bio-medical applications,
such materials are used in the manufacturing of implants in bone replacement,
Pompe et al. [26]. Also in electric power equipment such materials are employed
in designing high voltage gas insulated switch-gears and transformers with higher
electric eld conditions, Okubo et al. [25].
1.3 A Brief History on Homogenization Theory
Most of the important developments of physics during the twentieth century were
concerned with describing relations between dierent scales, and it is quite inter-
esting how researchers have dealt with this extraordinary challenge of discovering
what happens at mesoscopic and microscopic scales, often from observations made
at a macroscopic level. In order to overcome these serious issues, mathematical
theory of homogenization was started around 1960-1970's by Italian pioneers, Sergio
Spagnolo, Ennio De Giorgi, and Antonio Marino [22].
9
Composites are characterized by the fact that they are made up of two or
more materials, therefore the problem of nding eective overall properties of
structures with rapidly varying material properties is of importance. In any
composite, the scale of heterogeneities is small compared to the structure global
dimensions. Consequently, two scales characterize the composites are required, one
the microscopic which describes the micro-structural behavior and heterogeneity
of micro-structures, second, macroscopic scale which describes the global behavior
of the composites. From the macroscopic point of view, the composite looks
like a homogeneous material. The aim of homogenization is precisely to give the
macroscopic properties of the composite by taking into account the properties of
the microscopic structure, Doina and Donato [6].
1.4 A Review of the Literature
Determining the eective material, strength and other physical properties of in-
homogeneous materials dates back to the end of the 18
th
century. Most of the
methods developed in the beginning of 19
th
century were limited to certain types
of variations in ber and matrix volume fraction, such as one ber surrounded by
a semi-innite to innite matrix medium. Although, some of the past methods
were accurate for simple cases, it was cumbersome to solve complicated problems
and not feasible to implement on computers. It is worthwhile to point out that
micro-mechanical approaches are some kind of approximate methods are required to
characterize the behavior of composite materials. All the approaches are constrained
by the geometry of inclusions and the overall domain. Further, it is practically
impossible, and generally undesirable to consider the actual spatial distribution of
the inclusions inside a matrix/whisker.
10
Consequently, during the 1970's the use of averaging techniques or homogeniza-
tion drew the attention of many researchers in the development of new methods
for weighting the overall properties of the composite medium. In addition, modern
researchers started to exploit computer algorithms to reduce computational time
for composite material analysis with greater eciency.
1.4.1 Micro-Mechanical Analysis of Composites
These were the preliminary methods developed to determine the elastic moduli
of composite materials in the early past century. These methods were precise
but unreliable for the higher order problems. They are based on the strength
of the materials, semi-empirical methods, elasticity approach, the last one is the
most cumbersome and computation unfriendly method which involves solving
the analytical methods by dierential equations approach. Some of the simple
micro-mechanical methods developed are briefed in this work, including other
methods.
1.4.1.1 The Voigt Approximation
This was the rst model introduced by Voigt in 1887. It nds the eective material
stiness as the comnination of the individual material stinesses weighted by the
appropriate volume fraction, it assumes that the strain is constant throughout the
composite domain. It is given by,
[E] =
f
[E
f
] + (1
f
) [E
m
] (1.1)
11
where [E] is the eective material stiness matrix of the composite, [E
f
] is the
stiness matrix of the ber, [E
m
] is the stiness coecients of the matrix material,
and
f
is the ber volume fraction.
1.4.1.2 The Reuss Approximation
This method was developed by Reuss in 1929. It is similar to Voigt's method, but
instead of stiness it uses the eective compliance which is a weighted combination
of individual materials compliance,
[S] =
f
[S
f
] + (1
f
) [S
m
] (1.2)
where [S] is the eective compliance matrix of the composite, [S
f
] is the compliance
matrix of the ber, [S
m
] is the compliance coecients of the matrix material, and

f
is the ber volume fraction.
1.4.2 A Self-Consistent Method
Hill [13] proposed a self-consistent method for determining the material properties
of the composites. In this model a single ber is surrounded by the innite spread
of the matrix material. A unit strain in the ber can be produced by applying
a uniform force on the boundary of the continuum. The uniform strain is then
assumed to be the average over all the bers in the composite. Although this
method is easy to implement and provides reliable results for the material properties,
its accuracy diminishes if high ber volume is considered.
12
1.4.3 The Mori-Tanaka Method
Mori and Tanaka [24] discussed a method of calculating the average internal stress
in the matrix of a material containing inclusions with transformation strain. It was
shown that the average stress in the matrix is uniform throughout the material and
independent of the position of the domain where the average treatment is carried
out. It is also shown that the actual stress in the matrix is the average stress plus
the locally
uctuating stress, the average of which vanishes in the matrix. Average
elastic energy is also considered by taking into account the eects of the interaction
among the inclusions and of the presence of free boundary. This method closely
predicts the overall properties and local elds of continuous matrix in discontinuous
bers, but suers with dense ber volume interaction.
1.4.4 The Method of Cells
Aboudi [1], has written numerous papers on the method of cells, wherein he assumes
that the composite materials follows the periodic pattern in small scale and thereby
extracts a representative volume element
1
around a region of interest is as shown in
Figure 1.5 . A periodic region, called as an unit cell, is considered for the modeling
of inclusion and its surrounding matrix. The eective stiness matrix is derived by
relating the average stresses to the average strains inside the micro-cells and then
averaged over the volume of the unit cell. This method closely follows the Hill's
principle [12] and is limited to the periodic structures.
1
\A representative volume element (RVE) of a material is the smallest part (a square in Figure
1.5) of the material that represents the materials as a whole. It could be otherwise intractable to
account for the distribution of the constituents of the material", Kaw [3].
13
Figure 1.5: A Typical Representative Volume Element
1.4.5 Homogenization Methods
Homogenization is a most widely used method in the areas of mathematics, materials
sciences, meteorology and computational mechanics. This method had its impact
since 1970's when the highly challenging and oscillatory problems need to be solved
accurately. There are enormous number of publications and journal articles on this
method produced by the researchers in the eld of engineering and computational
mechanics, and it is dicult to provide a truly comprehensive bibliography. Some
of the past promising and reliable works are presented here.
Bendse and Kikuchi [4], used the method of homogenization for the optimal
shape design of structural elements based on the boundary variations. They used
the method to turn the shape optimization problem into one of nding the optimal
distribution of materials. Also, they constructed composite with a hole, and
determined the eective properties of materials by the regular and adaptive nite
element methods. Guedes and Kikchi [10], presented the mathematical principles
and nite element formulation for the composite materials with periodic micro-
structures. They implemented three-dimensional homogenization and nite element
14
formulation for their work. In addition, they solved problems by utilizing the
adaptive nite element techniques to minimize mesh size related errors.
Hollister and Kikuchi [15], studied the computational procedure combining
representative volume element based homogenization theory with digital imaging
which is proposed to estimate strains at various levels of bone structure. Homoge-
nization analysis of an idealized model was used to estimate strains at one level of
bone structure around osteocyte lacunae (second-level trabecular micro-structure).
The results showed that strain at one level of bone structure is amplied to a
broad range at the next micro-structural level. These results may be used to bridge
experimental studies of bone adaptation at the whole bone and cell culture level.
However, accuracy of homogenization RVE analyses is limited by a number of fac-
tors including RVE boundary condition assumptions, lack of experimental stiness
data on bone micro-structural components, inaccuracies in thresholding digital
images of bone micro-structure, and numerical artifacts resulting from non-smooth
nite element meshes made from digital images. Displacement boundary condition
assumptions will theoretically lead to an under-prediction of both eective stiness
and local strain. Most crucial limiting factor is uncertainty in relationship between
measured stiness and digital image voxel density. Terada et al. [31], presented the
systematic methodologies to derive accurate micro-structural models for studying
the mechanical behaviors of composite materials by imaging techniques. They,
implemented this approach to demonstrate the potential of digital image-based
(DIB) nite modeling of the unit cell and further using the homogenization method
for obtaining the eective properties.
Kouznetsova et al. [18], studied a gradient-enhanced computational homogeniza-
tion procedure, that allows for the modeling of micro-structural size eects within a
15
general non-linear framework. In this approach the macroscopic deformation gradi-
ent tensor and its gradient are imposed on a micro-structural representative volume
element (RVE). This enables the incorporation of the micro-structural size and the
eect of non-uniform macroscopic deformation elds within the micro-structural
cell. Every micro-structural constituent is modeled as a classical continuum and
the RVE problem is formulated in terms of standard equilibrium and boundary
conditions.
Ghosh et al. [11], presented a multiple scale nite element method combining
the asymptotic homogenization theory with Voronoi cell nite element method
(VCFEM) for micro-structural modeling. In this study, a micro-structural Voronoi
cell nite element method is coupled with the homogenization method for performing
multiple scale analysis of heterogeneous structures with arbitrary micro-structures.
The microscopic Voronoi cell nite element model is created by Dirichlet tessellation
of representative material elements or base cells. The coupled model projects a
strong analysis tool for nding eective (homogenized) material properties as well
as distributions of evolving macroscopic and microscopic variables.
Chung and Raju [20], studied materials with periodic as well as non-periodic
micro-structures. The method is validated through direct comparison with a known
exact solution, that being a rod with linearly varying axial rigidity. The method for
this case is shown to be exact. Further application of their method is used to verify
two dierent periodic materials with periodic assumption employed in traditional
homogenization methods.
A higher order theory was proposed by Aboudi et al. [2]. This theory allows
the thermo-inelastic analysis of materials with spatially varying micro-structure
based on volumetric averaging of the various eld quantities and satisfying of the
16
eld equations. Homogenization is implemented on a regular shaped rectangular
unit cell.
Besides all these techniques used to approximate the micro-structural properties,
there are some other researchers who followed the method of dividing the domain. A
common procedure for analysis of functionally graded materials involves separating
the inclusions by their volume and material density and performing the conventional
homogenization with nite element analysis for each separated layer with simplied
boundary conditions. Figure 1.6 shows a typical layered material transformed from
their original spatial distribution of material density.
Figure 1.6: A Typical Layered Functionally Graded Materials
Reiter et al. [29], studied the micro-mechanical model of two-phase graded
materials transformed to a single volume fraction gradient. The overall elastic
response can be obtained from homogeneous layer models, with the layer moduli
estimated in terms of the local volume fractions, and matrix and inclusion moduli,
by certain existing self-consistent micro-mechanical methods.
Similarly, Vemaganti and Deshmukh [32], used a layer-wise approach to estimate
the overall properties of in-homogeneous materials and to minimize the error due
to layering based on adaptive nite element techniques. Santare and Lambros [30],
17
used a varied material property matrix in the nite element with a simple boundary
condition. Lastly, Kim and Paulino [16], used an exponential variation of material
elastic properties to model non-homogeneous, isotropic and orthotropic, materials.
1.5 Scope of This Thesis
In this thesis, problems involving heterogeneity are considered to solve for the
structural/strength behavior. Especially, functionally graded materials/micro-
structures are taken into the account for the structural solutions. The method of
aperiodic homogenization has a unique capability compared to other techniques, the
two scale approximation of the displacements in micro-macro coupling eld could
be extracted precisely by the transformation of macro displacements to the micro
displacements through the highly characterized transformation equations. Also, a
successful attempt has been made to exploit the homogenized material model for
dierent scales. The eect/importance of micro-macro coupling equations could be
visualized by the Figure 1.7.
From the Figure 1.7 one can visualize that there exists a signicant relationship
between quantum to macro scale. During the manufacturing process of any materials
there is a vital transition of billions of atoms together to form molecules, then
micro-structures and lastly, the macro-structures by the assembly of tiny bits and
pieces. It becomes crucial to develop relationship between macroscopic properties
and micro-structures of a composite. \At the smallest/nano scale, an atomistic or
quantum analysis can assess features such as interface fracture energy and crack
de
ection at the bi-material ber/matrix interface" [34]. \Key information on
interfacial debonding and sliding is then passed into a continuum interface model,
such as, a cohesive zone, used in a micro-mechanical unit cell model consisting
18
Figure 1.7: A Micro-Macro Scale Coupling [34]
of matrix and a number of bers to compute the stress redistribution around a
ber break for a particular material system" [34]. In the macro scale, classical
lamination theory, strength of materials method, homogenization method and other
techniques could be used in order to nd the macroscopic properties.
The presentation of the thesis begins with a detailed discussion of asymp-
totic/periodic and aperiodic homogenization procedures are in the Chapter 2. This
chapter is intended to familiarize the reader with the mathematical specications
and theory of asymptotic homogenization method for periodic structures that is of-
ten used to obtain the homogenized material constants. In addition, the mechanics
related principles behind the method of aperiodic homogenization are explained for
the micro-structures. Further, the anti-periodic displacements elds are used to
derive the homogenized material constants relationship with the constant strain
eld. Also a one dimensional problem is considered to provide a detailed analytical
19
solution for use in obtain the homogenized elastic constants. In addition, two-
dimensional solutions for homogenized constants are compared for periodic and
aperidocic micro-structures.
The nite element formulation used in the aperiodic homogenization technique
is outlined in the Chapter 3. This includes the utilization of micro-macro stiness,
transformation and solution variables relations. A general purpose nite element
method is implemented for obtaining the homogenized material constants for
each unit cell, and then the material constants are utilized at Gauss Quadrature
points of each macro elements for gathering the strength parameters. Due to the
heterogeneity its highly cumbersome to mesh the micro elements and to extract
the material constants like Young's modules and Poisson's ratio at each micro
element. Therefore, a meshing technique based on digital image processing is
developed for functionally graded materials, which in turn eciently provides the
geometric properties (meshed coordinates and others) and material parameters
for each meshed elements. Although the aperiodic homogenization technique is
capable of solving problems with any heterogeneous micro-structures, the scope of
this work is limited to the structural analysis of functionally graded materials.
The numerical results and simulations of the aperiodic homogenization method
with the viable digital image tool for the meshing process are presented in the
Chapter 4. The results of two-dimensional problems with variety of ber/matrix
volume fraction are considered in this chapter. The detailed results for the aperiodic
homogenization technique are compared with the Dassault Syst emes, Abaqus,
Inc., Abaqus
R
, a commercial nite element package. Finally, an overview of this
thesis, conclusions deduced from the results and method developed, and promising
techniques in the future for other applications are detailed in the Chapter 5.
20
Chapter 2
Homogenization Strategies for
Heterogeneous Materials
In homogenization theory it is assumed that the composite materials is locally
formed by the spatial repetition of tiny micro-structures, i.e. microscopic cells, when
compared with the overall `macroscopic' dimensions or the complete domain of
interest. Also, it is considered that the material properties are periodic and follows
the periodic boundary conditions. But in reality this assumption may fail to apply
for most of the cases when the material properties near the boundary of the domain
are of concern. Therefore, to overcome the periodic scenario, a non-periodic method
is developed with the background of asymptotic homogenization for obtaining the
equivalent material constants throughout Gauss quadrature points in the domain
of interest.
In this chapter a rigorous mathematical analysis of asymptotic homogenization
is presented and nite element implementations are discussed. For additional
information, the readers could refer to an excellent resource by Guedes and Kikuchi
[10]. In addition to this, aperiodic homogenization theory is discussed with its
formulation and assumptions. Lastly, as a validation, the results obtained by the
aperiodic method are compared with accepted solutions.
21
2.1 A Mathematical Review of the Asymptotic
Homogenization Theory
Consider an arbitrary composite material formed of spatial repetition of a base cell
made of dierent materials as shown in Figure 2.1. For simplicity, assume that the
body is made of two dierent materials, although it could be made of n dierent
materials. The mixture of two materials is represented by a unit/base cell that is
very tiny, of a magnitude % (where % is a very small positive number) in contrast
with the domain dimensions of the structural body, which is shown in Figure 2.2 at
point `A'.
f
T3
T2
T1
Gt
Gd
X
A
W
Figure 2.1: An Arbitrary 2-D Shaped
Body
Figure 2.2: A detailed view of Micro-
structures at a point `A'
If the body is subjected to load f, traction T
1
; T
2
and T
n
and some boundary
conditions
t
and
d
, the deformation and stresses resulting load would vary rapidly
from point to point in the domain because of the repetition of the base cell or
Representative Volume/Area Element (RVE/RAE). Mathematically, with the high
level of heterogeneity within the material, the displacements and stresses vary
22
rapidly within a very small neighborhood % of any given point x. Therefore, it
is valid to claim that all quantities have two explicit dependencies. First is on
the macroscopic domain x and other in microscopic level x=%. Let g be a general
function which is comprised of these two dependencies i.e. g =g (x;x=%). It could
be noted that, due to the periodic nature of the micro-structure within the domain,
the dependency of the function on the microscopic variabley =x=% is also periodic.
Solving this nature of problem using nite element method would result in huge
computational cost and large degrees of freedom and also this would jeopardize
the accuracy of the method for some random material discretization and explicit
representation in an unit cell. Therefore, a less time consuming method would be
necessary to account for the mechanical behavior of the macroscopic composite
body. Further, this mathematical analysis of asymptotic method would dene a
path for aperiodic method in later stages of this chapter.
y
Y
Y
S
f
p
Figure 2.3: A Typical Base or Unit Cell of a Composite
Let
be an open subset ofR
2
(since the derivation is in terms of two-dimensional
material matrix) with a smooth boundary as shown in Figure 2.1. Let Y be an
23
open rectangular parallelepiped inR
2
, as shown in Figure 2.3 and mathematically
dened by
Y =]0; y
0
1
[x]0; y
0
2
[x]0; y
0
3
[; (2.1)
The above equation physically shows that the base cell is surrounded by innite,
continuous and repetitive base cells. Let # be an open subset of Y with boundary
@# = S (2.2)
and let
U =Yn

# (2.3)
Where the U is the solid part of the cell,

# represents the closure of # and Y
represents the base cell of the composite micro-structure. The material properties
like Young's modulus and Poisson's ratio vary inside Y, # represents a hole inside
Y. Now,
(y) =
8
>
<
>
:
1 if y2U
0 if y = 2U
(2.4)
and extend toR
2
by % periodicity, i.e. repeat the base cell in all two directions.
Then dene

%
=fx2
j (x=%) = 1g (2.5)
i.e.,

%
is the solid part of the domain. Also dene S
%
=
all cells
[
=1
S

Consider the following assumptions:

%
is a connected domain
The closure/hole(s) # has suciently smooth boundary(ies) S
24
None of the holes intersects the boundary of

There exists a unique solution of displacements for any loads and boundary
conditions.
Let,
V
%
=

2 (H
1
(

%
))
3
jj

d
= 0

(2.6)
Where j

d
represents the value on the boundary
d
.
To nd the solution of deformation of a body

%
due to the body force (f) and
tractions (T
1
; T
2
and T
n
) on the boundary
t
with the traction p inside the hole,
it could be stated that,
Find u
%
2 V
%
, such that
Z

%
E
%
ijkl
@u
%
k
@x
l
@
i
@x
j
d
=
Z

%
f
%
i

i
d
+
Z
t
T
%
i

i
d
+
Z
S%
p
%
i

i
dS 82V
%
(2.7)
It is assumed that the stress-strain and the strain-displacement relations are:

%
ij
=E
%
ijkl

%
ij
(2.8)

%
kl
=
1
2

@u
%
k
@x
l
+
@u
%
l
@x
k

(2.9)
and the material matrix/constants have the following properties
E
%
ijkl
=E
%
jikl
=E
%
ijlk
=E
%
ijlk
(2.10)
9> 0 :E
%
ijkl

ij

kl
=
ij

kl
8
ij
=
ji
(2.11)
A unique solution u
%
exists for the Eq. 2.7, under the assumptions f, T and p
are smooth, and the boundaries
d
,
t
and S
%
are regular. Due to the fact that
25
the body forces, traction and the elastic constants vary within the base cell of a
composite, they are the functions of both x and x=% and y =x=%:

%
(x) = (x;y); y =x=% (2.12)
therefore the solution u
e
should also depend on both x and x=%, i.e.
u
e
(x) =u(x;y); y =x=% (2.13)
Solution dependency ony =x=% means that the quantity is dierent within a small
sub-patch with dimensions smaller than those of macroscopic scale. In the vicinity
of a macroscopic pointx there is a large repetition of small cells, which are periodic
and two-dimensional translation of a base cell. Although, one should note that for
dierent positions of x the composite materials structure could be dierent, but if
looked at a microscopic scale the periodicity could be observed as in Figure 2.2.
The close relationship of the solution u
%
with x, x=% and y =x=% i.e. in both
macroscopic and microscopic scale makes it reasonable that u
%
can be expressed as
an asymptotic expansion with respect to the small scale parameter % (a measure of
the microscopic/macroscopic dimension ratio). Therefore the asymptotic expansion
follows as:
u
%
(x) =u
0
(x;y) +%u
1
(x;y) +%
2
u
2
(x;y) + ; y =x=% (2.14)
where,
u
j
(x;y) is dened in (x;y)2
xU
y!u
j
(x;y) is Y periodic
(2.15)
26
To satisfy the Eq. 2.14 with u
0
; u
1
;:::; u
j
, it is useful to note that
@
@x
i
((x;y =x=%)) =
@
@x
i
+
1
%
@
@y
i
(2.16)
and for a Y periodic function (y),
lim
%!0
+
Z

%

x
%

d
!
1
jYj
Z

Z
U
(y)dY d
(2.17)
lim
%!0
+
Z
S
%

x
%

dS!
1
jYj
Z

Z
S
(y)dSd
(2.18)
WherejYj stands for the volume (or area, for two-dimensional domain) of the cell.
It is assumed that the function (y) is extended to all volume within each cell,
such that it takes the value zero on the holes. This is assumed to be smooth. Then,
dene,
V

U
=f((x;y) dened for (x;y)2
Uj(:;y)
Y-periodic; j

d
= 0; smooth enoughg; (2.19)
V

=f((x) dened in
jj

d
= 0; smooth enoughg; (2.20)
V
U
=f((y) dened inUj(y); Y-periodicsmooth enoughg; (2.21)
Introducing the asymptotic expansion 2.14 and 2.16 in Eq. 2.7, the functional takes
the form of,
Z
%
!
E
%
ijkl

1
%
2
@u
0
k
@y
l
@
i
@y
j
+
1
%

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
@x
j
+
@u
0
k
@y
l
@
i
@x
j

+

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
@y
l
+

@u
1
k
@x
l
+
@u
2
k
@y
l

@
i
@y
l

+%(:::)gd

27
=
Z
%

f
%
i

i
d
+
Z
t
T
i

i
d +
Z
%
S
p
%
i

i
dS 82V

U
(2.22)
Assuming that the functions are smooth enough so that the limte when %! 0
+
of
all integrals exists. Therefore,
1
%
2
Z

%
E
%
ijkl
@u
0
k
@y
l
@
i
@y
j
d
= 0 82V

U
(2.23)
1
%
Z

%
E
%
ijkl

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
@x
j
+
@u
0
k
@y
l
@
i
@x
j

d
=
Z
%
S
p
%
i

i
dS (2.24)
82V

U
1
%
Z

%
E
%
ijkl

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
@x
j
+

@u
1
k
@x
l
+
@u
2
k
@y
l

@
i
@y
l

d

=
Z

%
f
%
i

i
d
+
Z
t
T
i

i
d 82V

U
(2.25)
multiplying euqation 2.23 by %
2
and taking the limit %! 0
+
, property 2.17 yields,
1
jYj
Z

Z
U
E
%
ijkl
@u
0
k
@y
l
@
i
@y
j
dY d
= 0 82V

U
(2.26)
Since is arbitrary, choose = (y), i.e. 2 V
U
. Then integrating by parts,
applying the divergence theorem to the integral inU, and noting that due to the
periodic condition the opposite sides gets canceled out,
1
jYj
Z

Z
U

@
@y
j

E
%
ijkl
@u
0
k
@y
l

i
dY +
Z
S
E
%
ijkl
@u
0
k
@y
l
n
j

i
dS

d

= 0 82V

U
(2.27)
28
Since (y) is arbitrary, this yields the boundary value problem of the rst term u
0
of the expansion of the original solution u
%
in the basic cell domainU:

@
@y
j

E
%
ijkl
@u
0
k
@y
l

= 0; y2U (2.28)
E
%
ijkl
@u
0
k
@y
l
n
j
= 0 onS (2.29)
Apply the following proposition,
@
@y
j

E
ijkl
(y)
@(y)
@y
i

=F (y); y2U (2.30)
has a solution in
V
U
=f is smooth enough, is Y-periodicg (2.31)
dened up to an additive constant, for a regular F, if and only if
Z
U
F (y)dY =
Z
S
E
ijkl
@(y)
@y
i
n
i
dS (2.32)
then it follows from Eqs. 2.28 and 2.29 that,
u
0
=u
0
(x) (2.33)
that is, the rst term of the expansion ofu
%
in% depend only the macroscopic scale
x.
29
Introducing Eq. 2.33 into Eq. 2.24, multiply by %, taking the limit %! 0
+
)
and using properties Eqs. 2.17 and 2.18, results in,
Z

1
jYj
Z
U
E
ijkl

@u
0
k
(x)
@x
l
+
@u
1
k
@y
l

@
i
@y
j
dY

d

=
Z

1
jYj
Z
S
p
i

i
dS

d
82V

U
(2.34)
since Eq. 2.34 is satised for any , choosing =(y) yields,
Z
U
E
ijkl

@u
0
k
(x)
@x
l
+
@u
1
k
@y
l

@
i
@y
j
dY =
Z
S
p
i

i
(y)dS 82V

U
(2.35)
Integrating by parts, using the divergence theorem, and applying the periodicity
conditions and cancelling the opposite sides of Y , then from Eq. 2.35,

Z
U
@
@y
j

E
ijkl

@u
0
k
(x)
@x
l
+
@u
1
k
@y
l

i
(y)dY +
Z
S
E
ijkl

@u
0
k
(x)
@x
l
+
@u
1
k
@y
l

i
(y)n
j
dS =
Z
S
p
i

i
(y)dS 82V

U
(2.36)
Since is arbitrary,

@
@y
j

E
ijkl
@u
1
k
@y
l

=
@
@y
j

E
ijkl
@u
0
k
(x)
@x
l

onU (2.37)
E
ijkl
@u
1
k
@y
l
n
j
=E
ijkl
@u
0
k
(x)
@x
l
n
j
+p
i
onS (2.38)
If in Eq. 2.34 is chosen to be only a function of x, i.e =nu(x), then
Z

1
jYj
Z
S
p
i
dS

i
(x)d
= 0 82V
!
(2.39)
30
which implies that,
Z
S
p
i
(x;y)dS = 0 (2.40)
Eq. 2.40 restricts the possible traction inside the closure. Therefore, the other
applied traction is in self-equilibrium with respect to the outer surface. From the
proposition 2.30 there exists a solution u
1
2V
U
dened up to a constant.
Thereby introducing the Eq. 2.33 into Eq. 2.25 and taking the limit %! 0
+
in
Eq. 2.25 implies that,
Z

1
jYj
Z
U
E
ijkl

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
@x
j
+

@u
1
k
@x
l
+
@u
2
k
@y
l

@
i
@y
j

dY

d

Z

1
jYj
Z
U
f
i

i
dY

d
+
Z
t
T
i

i
d 82V

U
(2.41)
choosing =(x) yields to,
Z

1
jYj
Z
U
E
ijkl

@u
0
k
@x
l
+
@u
l
k
@y
l

@
i
(x)
@x
j

dY d

Z

1
jYj
Z
U
f
i

i
dY

d
+
Z
t
T
i

i
d 82V

(2.42)
It should be noted that this is a statement of overall equilibrium in the macroscopic
eld. Now if =(y) is assumed in Eq. 2.41, then
Z

1
jYj
Z
U
E
ijkl

@u
1
k
@x
l
+
@u
2
k
@y
l

@
i
(y)
@y
j

dY

d

=
Z

1
jYj
Z
U
f
i

i
(y)dY

d
82V
U
(2.43)
31
This equivalent to,
Z
U
E
ijkl

@u
1
k
@x
l
+
@u
2
k
@y
l

@
i
(y)
@y
j
dY =
Z
U
f
i

i
(y)dY
82V
U
(2.44)
and represents the microscopic equilibrium of the base cell.
The objective is to nd homogenized elastic constants such that the macroscopic
equilibrium can be described by the same sort of Eq. as 2.7. These homogenized
constants should be such that the corresponding equilibrium equation re
ects the
mechanical behavior of the micro-structure of the composite material without
explicitly using the micro-scale parameter %. Consider once again Eq. 2.35 and
recognizing that this equation is linear with respect to u
0
and p. This leads to
separate problems,
Let
kl
2V
U
be the solution of
Z
U
E
ijpm
@
kl
p
@y
m
@
i
(y)
@y
j
dY =
Z
U
E
ijkl
@
i
(y)
@y
j
dY 82V
U
(2.45)
and let 2V
U
be the solution of
Z
U
E
ijkl
@
k
@y
l
@
i
(y)
@y
j
dY =
Z
S
p
i

i
(y)dY 82V
U
(2.46)
Wherex plays the role of a parameter. The solutions to these problems are assured
by the proposition and the function of x, the solution u
1
can be written as.
u
1
i
=
kl
i
(x;y)
@u
0
k
(x)
@x
l

i
(x;y) + ~ u
1
i
(x) (2.47)
32
Where ~ u
1
i
(x) are arbitrary additive constants in y. Introducing the Eq. 2.47 into
Eq. 2.42 and simplifying results,
Z

"
1
jYj
Z
U

E
ijkl
E
ijpm
@
kl
p
@y
m
!
dY
#
@u
0
k
(x)
@x
l
@
i
(x)
@x
j
d

=
Z

1
jYj
Z
U
E
ijkl
@
k
@y
l
dY

@
i
(x)
@x
j
d
+
Z

1
jYj
Z
U
f
i
;dY

i
(x)d
+
Z
t
T
i

i
(x)d 82V

(2.48)
It can be also denoted as follows,

E
H
ijkl
(x) =
1
jYj
Z
U

E
ijkl
E
ijpm
@
kl
p
@y
m
!
dY (2.49)

ij
(x) =
1
jYj
Z
U
E
ijkl
@
k
@y
l
dY (2.50)
b
i
(x) =
1
jYj
Z
U
f
i
dY (2.51)
The nal form of the Eq. 2.48 takes the following form:
Z

E
H
ijkl
(x)
@u
0
k
(x)
@x
l
@
i
(x)
@x
j
d
=
Z

ij
(x)
@
i
(x)
@x
j
d
+
Z

b
i
(x)
i
(x)d
+
Z
t
T
i
(x)
i
(x)d 82V

(2.52)
Eq. 2.52 ensures the macroscopic equilibrium, while the homogenized elastic
constants,

E
H
ijkl
is dened by the Eq. 2.49 and the average residual stress within
the cell due to the tractions p inside the holes is given by
ij
and b
i
is the average
body force.
33
From the above equation it can be noted that the microscopic and macroscopic
problems are not coupled, i.e., the homogenized elastic tensor can be computed
within the basic cell by solving the Eqs. 2.45 and 2.46, which do not depend on the
macroscopic deformationu
0
. Supposedly, if a composite material has a uniform cell
structure throughout the whole domain
, as well as a uniform applied traction p
on the boundaries of the closure of the cells, then the microscopic Eqs. 2.45 and
2.46 need to be solved only once. For an opposite scenario, if they are not uniform
in the domain
, these microscopic equations must be solved for every point x of

at which cell structure and traction p are dierent from others.
Drawbacks of Asymptotic Homogenization Method
Mathematically complicated and cumbersome
Requires higher order terms to be included for more accurate results and
higher convergence
Unfeasible for some of the problems with discontinuous material distribution
Micro-structures must be characterized by the periodic assumption and
restricted boundary conditions
Inapplicable for nonlinear problems
No explicit coupling of multi-scale properties
34
2.2 An Aperiodic Homogenization Method
In this method the eective properties of microscopic structures are considered
in order to dene the macroscopic behavior. Micro-structures are made up of an
hierarchy of two or more dierent materials. The characteristic length of the entire
body is much larger than the characteristic length of the micro-structure.
Non-periodic
Boundary
W
Y
W
Non-periodic
Boundary
Y
Figure 2.4: Homogenization of a Composites in Two Scales
In Figure 2.4, a sub-patch in certain domain
is considered. The sub-patch Y
has relatively small sub-domain with respect to
. It is practically impossible to
determine the mechanical properties at each point x in domain
. To address this
diculty, an eective phase property at each sub-patch are taken into account, and
later homogenized constants are obtained. By considering the specic sub-patch
around a region/point of interest it reasons for the more ecient way of analysis
and minimizing the computational expenses. Further, the microscopic quantities
can be transformed to obtain the macroscopic quantities like displacements, forces,
and stresses.
35
2.2.1 Assumptions and Geometry of a Sub-Patch
Materials
{ It is assumed that the material is linearly elastic and isotropic within
a sub-patch refer Figure2.5. Although the domain is comprised of
composite materials, all materials inside a micro-cell are homogeneous
and isotropic in nature.
E
f
, n
f
E
m
, n
m
Figure 2.5: A Typical sub-patch with Two Materials
This gure shows a sub-patch with two dierent materials randomly
distributed and having E
f
and E
m
for the elastic modulus of bers and
matrix correspondingly. Similarly,
f
and
m
for the Poisson's ratio of
bers and matrix.
36
Geometry
{ It is assumed that the geometry of a sub-patch is a regular quadrilateral
i.e. a square or a rectangle.
{ The sub-patches are aligned with the global element i.e. unrotated with
respect to the global axis (X,Y), but the elements within the sub-patch
can be distorted
{ Geometry of inclusions are random in nature within domain
and
sub-patch domain Y
{ Iso-parametric method is considered for the analysis of a sub-patch,
where the elements are 4 noded bi-linear in nature
{ Geometric non-linearity of the domain is ignored
Boundary Conditions
{ The surrounding boundaries of a sub-patch are discontinuous (non-
periodic)
{ May not necessarily have periodic boundaries
{ Sub-patch or unit cell or base cell or RVE could intersect the boundary
of a domain

{ Displacements and strains vary for each individual sub-patch, even if it
is adjacent
{ Each sub-patch satises the rigid body modes
{ Perfect bonding between matrix and inclusions
37
2.3 Determination of the Homogenized Elastic
Constants
Homogenized elastic constants are purely dependent on the material distribution
and constraints along the boundary of a sub-patch. The general principle is, the
eective properties are represented by the eective stiness, in terms of the average
stress or strain [27]. In other words, the average properties of all the micro elements
within a sub-patch/base cell is equivalent to the property of a macro element or
base cell itself.
The strain energy is dened as,
U =
1
2
Z

ij

kl
d
(2.53)
where
ij
is the stress acting at any point of the body
. According to the
equivalence of the strain energy, one can dene,
1
V
Z

1
2

ij

ij
d
=
1
2

ij

ij
(2.54)
But expanding the stress terms,
ij
and
ij
, leads to
1
V
Z

1
2

E
H
ijkl

ij

kl
d
=
1
2
E
ijkl

kl

ij
(2.55)
Eq. 2.55 is obtained by Hill [12] and is referred to as Hill's principle. In Eq. 2.55,

E
H
ijkl
denotes the undetermined elastic constants, E
ijkl
represents the overall/total
material stiness within a base cell and lastly,
ij
stands for the strains. In this
work, a dierent form of the Eq. 2.55 is used to obtain the homogenized elastic
38
constants. More or less, the indices in the Eq. 2.55 are expanded and simplied,
which provides the required materials stiness constants.
2.3.1 Aperiodic Elastic Constants
As briefed earlier, the eective material properties plays a preliminary role in deter-
mining the macroscopic structural properties. The rigorous averaging techniques
proves to be a viable tool for nding the multi-scale coecients of interest. In this
section, a detailed approach to obtain the elastic cocients are presented with a
strong foundation of the earlier concepts of mechanics.
Consider a sub-patch made up of graded composites as shown in the Figure 2.6.
V
1
V
2
U
2
U
3
V
4
U
4
V
3
U
1
U
n
and V
n
Figure 2.6: A Sub-patch with Graded Materials
The base cell is 4 noded iso-parametric bi-linear element, comprised of 12x13 micro
cells of dierent materials .

U
n
and

V
n
are the microscopic displacements, similarly
39
U
1 4
and V
1 4
are the macro displacements dened only in the corner of an
element.
Now consider the microscopic displacements and their total strain energy,
U =
1
2
Z

ij

E
ijkl

kl
d
(2.56)
But,

ij
=
1
2

U
i;j
+

U
j;i

(2.57)
is the strain energy in tensorial form, where i =x direction and j =y direction.
The multi-scale displacement approximation for the displacement is given by
(details in Chapter 3),
U
i
=

J

U
J
i
(2.58)
and followed by the micro-macro coupling equation (details in Chapter 3)

U
J
i
=T
JKi
U
K
i
(2.59)
Therefore by considering the anti-periodic (transformation) equation,

U
i
=

J
T
JKi
U
K
i
(2.60)
dierentiate Eq. 2.60 w.r.t. y, then

U
i;j
=

J;j
T
JKi
U
K
i
(2.61)
consider

U
j
=

J
T
JKj
U
K
j
(2.62)
Similarly, dierentiate Eq. 2.62 w.r.t. x, then

U
j;i
=

J;i
T
JKj
U
K
i
(2.63)
40
Substituing, Eqs. 2.61 and 2.63 in Eq. 2.57 leads to,

ij
=
1
2

U
i;j
+

U
j;i

=
1
2

J;j
T
JKi
U
K
i
+

J;i
T
JKj
U
K
i

(2.64)
from theorem of tensors,
=
1
2

J;j
T
JKi

im
+

J;i
T
JKj

jm

| {z }

D
ijkm
U
K
m
Where
ij
is Kronecker Delta,

D
ijkm
and

D
klln
are the strain transformation matrices
of micro elements. Therefore,

ij
=

D
ijkm
U
K
m
(2.65)

kl
=

D
klln
U
L
n
(2.66)
using Eqs. 2.65 and 2.66, Eq. 2.56 take the form of:
U
total
=
1
2
Z
A
t

D
ijkm
U
K
m

E
ijkl

D
klln
U
L
n
dA (2.67)
where,

K
klmn
=
Z
A
t

D
ijkm

E
ijkl

D
klln
(2.68)
is the micro stiness Eq. and t is the thickness of the domain
and
U
total
=
1
2
U
K
m

K
klmn
U
L
n
(2.69)
is the total strain energy of the micro displacements.
41
Using same argument, now considering the formulation for a macro element,
Let,
U
i
=
K

U
K
i
(2.70)
Let

U's be the macro displacements and dierentiate Eq. 2.70 w.r.t. x
j
U
i;j
=
K;j

U
K
i
(2.71)
U
j
=
K

U
K
j
(2.72)
Similarly, dierentiate Eq. 2.72 w.r.t. x
U
j;i
=
K;i

U
K
j
(2.73)
strain energy is given by,

ij
=
1
2
(U
i;j
+U
j;i
) (2.74)
Substituing Eqs. 2.71 and 2.73 in Eq. 2.74,
=
1
2

K;j

U
K
i
+
K;i

U
K
j

(2.75)
=
1
2
(
K;j

im
+
K;i

jm
)

U
K
m
(2.76)
Where
ij
is Kronecker Delta and

D
ijkm
and

D
klln
are the strain transformation
matrices of macro element. Therefore,

ij
=D
ijkm

U
K
m
(2.77)

kl
=D
klln

U
L
n
(2.78)
Similar to Eq. 2.56 and 2.69,
U
Hom
=
1
2
Z
A
t
ij
E
ijkl

kl
dA (2.79)
U
Hom
=
1
2
Z
A
tD
ijkm

U
K
m
E
ijkl

D
klln

U
L
n
dA (2.80)
42
where
K
klmn
=
Z
A
tD
ijkm
E
ijkl
D
klln
(2.81)
is the macro stiness Eq. and t is the thickness of the domain
and
U
Hom
=
1
2

U
K
m
K
klmn

U
L
n
(2.82)
is the homogenized strain energy of the macro displacements. The displacemnts

U
K
m
or U
K
m
can be dened in terms of anti-periodic modes, U
K
m
. Therfore,

ij
=D
ijkm
U
K
m
(2.83)

kl
=D
klln
U
L
n
(2.84)
Consider the Eq. 2.79 and expanding on k = 1; 2 and 3 and l = 1; 2 and 3 (i.e. x, y
and z respectively and ignoring z dependency,
U
Hom
=
1
2
Z
A
t
ij
E
ijkl

kl
dA
=
1
2
Z
A
t (
ij
E
ij11

11
+
ij
E
ij12

12
+

ij
E
ij13

13
+
ij
E
ij21

21
+
ij
E
ij22

22
+

ij
E
ij23

23
+

ij
E
ij31

31
+

ij
E
ij32

32
+

ij
E
ij33

33
)dA (2.85)
=
1
2
Z
A
t (
ij
E
ij11

11
+
ij
E
ij12

12
+
+
ij
E
ij21

21
+
ij
E
ij22

22
)dA (2.86)
43
now expanding in terms of i; j=1, 2, and 3
=
1
2
Z
A
t
0
@

11
E
1111

11
| {z }

A
+
12
E
1211

11
| {z }

B
+
21
E
2111

11
| {z }

B
+
22
E
2211

11
| {z }

C
+
11
E
1112

12
| {z }

B
+
12
E
1212

12
| {z }

D
+
21
E
2112

12
| {z }

D
+
22
E
2212

12
| {z }

E
+
11
E
1121

21
| {z }

B
+
12
E
1221

12
| {z }

D
+
21
E
2121

21
| {z }

D
+
22
E
2221

21
| {z }

E
+
11
E
1122

22
| {z }

C
+
12
E
1222

22
| {z }

E
+
21
E
2122

22
| {z }

E
+
22
E
2222

22
| {z }

F
1
A
dA (2.87)
rearranging the terms in Eq. 2.87
=
1
2
Z
A
t
2
4

11
E
1111

11
| {z }

A
+
12
E
1211

11
+
11
E
1112

12
+
21
E
2111

11
| {z }

B
+
11
E
1121

21
| {z }

B
+
22
E
2211

11
+
11
E
1122

22
| {z }

C
+
12
E
1212

12
| {z }

D
+
21
E
2112

12
+
12
E
1221

12
+
21
E
2121

21
| {z }

D
+
22
E
2212

12
| {z }

E
+
22
E
2221

21
+
12
E
1222

22
+
21
E
2122

22
| {z }

E
+
22
E
2222

22
| {z }

F
3
5
(2.88)
rearranging the terms in Eq. 2.88 in a matrix form,
=
1
2
Z
A
t
8
>
>
>
>
<
>
>
>
>
:

11

22

12
9
>
>
>
>
=
>
>
>
>
;
T
2
6
6
6
6
4

A

C

B

C

F

E

B

E

D
3
7
7
7
7
5
| {z }

E
8
>
>
>
>
<
>
>
>
>
:

11

22

12
9
>
>
>
>
=
>
>
>
>
;
dA (2.89)
44
Where

E
3x3
is the homogenized material matrix and,

A = E
1111
(2.90a)

B = E
1112
+E
1121
=E
1211
+E
2111
(2.90b)

C = E
1122
=E
2211
(2.90c)

D = E
1212
+E
2112
=E
1221
+E
2121
(2.90d)

E = E
2212
+E
1222
=E
2221
+E
1222
(2.90e)

F = E
2222
(2.90f)
also, 2.89 can be written as,
U
Hom
=
1
2
Z

tfg
T
1x3
h

E
i
3x3
fg
3x1
dA (2.91)
fg
3x1
=
h

D
i
38
fg
81
(2.92)
h

D
i
38
is the strain transformaation matrix andfg
81
denotes the anti-periodic
displacements
By considering the indices i and j
fg
i
3x1
=
h

D
i
38

i

81
(2.93)
substituting the Eq. 2.93 in 2.91
U
Hom
=
1
2

i

T
81
Z

t
h

D
i
T
38
h

E
i
3x3
h

D
i
38
dA

i

81
(2.94)
=
3
X
i=1
3
X
j=1
1
2

i

T
81
Z

t
h

D
i
T
38
h

E
i
3x3
h

D
i
38
dA
| {z }
[

K]
88
! macro stiness matrix

i

81
(2.95)
45
or else,
U
Hom
=
1
2

i

T
81
h

K
i
88

i

81
(2.96)
Further considering j
U
i;j
Hom
=
1
2
Z

t

i

T
81
h

D
i
T
38
| {z }

i
13
h

E
i
3x3
h

D
i
38

j

81
| {z }

j
31
dA (2.97)
Consider i = 1 and j = 1,
U
1;1
Hom
=
1
2
Z
A
t
8
>
>
>
>
<
>
>
>
>
:

11
0
0
9
>
>
>
>
=
>
>
>
>
;
T
2
6
6
6
6
4

A

C

B

C

F

E

B

E

D
3
7
7
7
7
5
| {z }

E
8
>
>
>
>
<
>
>
>
>
:

11
0
0
9
>
>
>
>
=
>
>
>
>
;
dA (2.98)
Multiplying 2.98 leads to,
U
1;1
Hom
=
1
2
Z
A
t
11

A
11
dA (2.99)
But,
R
A
1dA:t =V , therefore
U
1;1
Hom
=
1
2
V
11

A
11
(2.100)
But for micro-domain,
U
total
=
1
2
fg
T
n1
[K]
nn
fg
n1
(2.101)
Although [K] is, [K]
nn
considering nn micro elements, but for consistency
[K]
nn
is transformed back to the macro domain, and it becomes [K]
88
U
i;j
total
=
1
2

i

T
81
[K]
88

j

81
(2.102)
Now, [K]
88
is the average matrix associated with the micro-structure
46
Considering i = 1 and j = 1,
U
1;1
total
=
1
2

1

T
81
[K]
88

1

81
(2.103)
From the theorem of eective properties [12],
U
1;1
Hom
=U
1;1
total
(2.104)
Substituing Eqs. 2.100 and 2.103 in Eq. 2.104 and simplifying provides,

1
2
V
11

A
11
=

1
2

1

T
81
[K]
88

1

81
(2.105)

A =
1
V
2
11

1

T
81
[K]
88

1

81
(2.106)
Now considering i = 1 and j = 2,
U
1;2
Hom
=
1
2
Z
A
t
8
>
>
>
>
<
>
>
>
>
:

11
0
0
9
>
>
>
>
=
>
>
>
>
;
T
2
6
6
6
6
4

A

C

B

C

F

E

B

E

D
3
7
7
7
7
5
| {z }

E
8
>
>
>
>
<
>
>
>
>
:
0

22
0
9
>
>
>
>
=
>
>
>
>
;
dA (2.107)
Multiplying 2.107 leads to,
U
1;2
Hom
=
1
2
Z
A
t
11

C
22
dA (2.108)
But,
R
A
1dA:t =V , therefore
U
1;2
Hom
=
1
2
V
11

A
22
(2.109)
From micro-domain,
U
1;2
total
=
1
2

1

T
81
[K]
88

2

81
(2.110)
U
1;2
Hom
=U
1;2
total
(2.111)
47
Substituing Eqs. 2.109 and 2.110 in Eq. 2.111 and simplifying provides,

1
2
V
11

C
12
=

1
2

1

T
81
[K]
88

2

81
(2.112)

C =
1
V
11

12

1

T
81
[K]
88

2

81
(2.113)
Similarly, for ij = 13; 33; 32; and 22,

B =
1
V
11

12

1

T
81
[K]
88

3

81
(2.114)

D =
1
V
2
12

3

T
81
[K]
88

3

81
(2.115)

E =
1
V
22

12

3

T
81
[K]
88

2

81
(2.116)

F =
1
V
11

12

1

T
81
[K]
88

2

81
(2.117)
Next it is nescessary to calculate the variables in homogenized constants. In
estimating the anti-periodic displacements
83
there are two approaches, one by
iso-parametric method and the other by assuming some constant displacements
associated with the deformation modes given in Figure 3.8. It should be noted
that, since the shear strains are opposite and equal in nature,
3
is the combination
of strain mode xy! 12 and yx! 21.
2.3.1.1 Computation of Anti-periodic Modes by an Iso-parametric
Method
Although, in this work, it is assumed that the macro-element geometry is regular in
shape (square or rectangle), this method works for both regular and distorted cases.
Consider a distorted macro-element as shown in Figure 2.7, the nodal numbering
system follows the general counter-clockwise convention.
48
X
1
, Y
1
X
Y
X
2
, Y
2
X
4
, Y
4
X
3
, Y
3
X
cen
, Y
cen

Figure 2.7: A Distorted Iso-parametric Finite Element
Compute X
cen
;Y
cen
,
X(; )j
==0
= [ (;)]j
==0
fX
n
g (2.118)
Y (; )j
==0
= [ (;)]j
==0
fY
n
g (2.119)
where, (;) are the 4 noded bi-linear shape functions
Compute 3 modal displacements at each node n as,
u
1
n
=X
n
X
cen
v
1
n
= 0
9
>
=
>
;
Mode 1!
1
(2.120)
u
2
n
= 0
v
1
n
=Y
n
Y
cen
9
>
=
>
;
Mode 2!
2
(2.121)
u
2
n
=Y
n
Y
cen
v
1
n
=X
n
X
cen
9
>
=
>
;
Mode 3!
3
(2.122)
49
For simplicity, constant strains such as -1 for negative direction and +1 for positive
direction could also be used, without losing generality. Table 2.1 shows some general
values for deformation modes.
Nodes
1

2

3
1 -1 0 -1
2 0 -1 -1
3 1 0 -1
4 0 -1 1
5 1 0 1
6 0 1 1
7 -1 0 1
8 0 1 -1
Table 2.1: Unit Anti-Periodic Displacements

11
;
11
and
12
are calculated by the product of transformation matrix

D and
i
,
where i = 1; 2; 3. Thereby considering the terms in normal and shear direction
only.
In summary

A =
1
V
2
11

1

T
81
[K]
88

1

81

B =
1
V
11

12

1

T
81
[K]
88

3

81

C =
1
V
11

12

1

T
81
[K]
88

2

81

D =
1
V
2
12

3

T
81
[K]
88

3

81

E =
1
V
22

12

3

T
81
[K]
88

2

81

F =
1
V
11

12

1

T
81
[K]
88

2

81

E
Hom
=
2
6
6
6
6
4

A

C

B

C

F

E

B

E

D
3
7
7
7
7
5
(2.123)
50
2.4 Validation of the Aperiodic Elastic Con-
stants
In this section homogenized elastic constants, obtained by aperiodic method, are
validated in solving of the general problems like homogeneous media with and
without Poisson's eect and other various periodic patterns of material distribution
provided in the literature. While one-dimensional problem is solved analytically in
this work, emphasis is given to two-dimensional problems with various gradation of
materials.
2.4.1 Case 1: One Dimensional Approach
To begin with, the accuracy of the aperiodic homogenization method is evaulated
in application to a thin rod problem.
Consider a typical one-dimensional bar with one macro-element and two nite
elements in micro-scale as shown in gures 2.8 and 2.9.
1 2
U U
L
1 2
A
1
, E
1
A
2
, E
2
1 2
Figure 2.8: A one-dimensional Macro-
structure
1 2
U U U
L
L
2
1
1 2 3
A
1
, E
1
A
2
, E
2
1 2 3
Figure 2.9: A one-dimensional Micro-
structures
Where A
1
, E
1
, L
1
and A
2
, E
2
, L
2
are the area, Young's modulus, length of micro
51
elements 1 and 2 respectively.

U
1;2;3
andU
1;2
are the micro and macro displacements.
From Eq. 3.31, the transformation matrix follows as below:
8
>
>
>
>
<
>
>
>
>
:

U
1

U
2

U
3
9
>
>
>
>
=
>
>
>
>
;
=
2
6
6
6
6
4
T
111
T
211
T
121
T
221
T
131
T
231
3
7
7
7
7
5
[T ]
32
8
>
<
>
:

U
1

U
2
9
>
=
>
;
(2.124)
Now including the functions,
8
>
>
>
>
<
>
>
>
>
:

U
1

U
2

U
3
9
>
>
>
>
=
>
>
>
>
;
=
2
6
6
6
6
4
1
x
1
1
L
+

1
1
L
x
1
1
L

1
1
L
1
x
2
1
L
+

1
2
L
x
2
1
L

1
2
L
1
x
3
1
L
+

1
3
L
x
3
1
L

1
3
L
3
7
7
7
7
5
[T ]
32
8
>
<
>
:

U
1

U
2
9
>
=
>
;
(2.125)
and micro stifness matrix,
(2.126)

K

=
2
6
6
6
6
4
AE
1
L
1

AE
1
L
1
0

AE
1
L
1
AE
1
L
1
+
AE
2
L
2

AE
2
L
2
0
AE
2
L
2
AE
2
L
2
3
7
7
7
7
5
(2.127)
Consider an anti-periodic eld as follows:
8
>
<
>
:
U
1
U
2
9
>
=
>
;
=
8
>
<
>
:
1
1
9
>
=
>
;
U (2.128)
Substituting Eq. 2.128 in Eq. 2.125,
8
>
>
>
>
<
>
>
>
>
:

U
1

U
2

U
3
9
>
>
>
>
=
>
>
>
>
;
=
2
6
6
6
6
4
1 +
2x
1
1
L

2
1
1
L
1 +
2x
2
1
L

2
1
2
L
1 +
2x
3
1
L

2
1
3
L
3
7
7
7
7
5
[G]
31
U (2.129)
52
Then the micro equation is given by,

K

33

U

31
= 0 (2.130)
Introducing Eq. 2.129 in Eq. 2.130

K

33
[G]
31
U = 0 (2.131)
Since U is arbitrary, assume it is 1
j
then

K

33
[G]
31
= 0 (2.132)
or,
2
6
6
6
6
4
AE
1
L
1

AE
1
L
1
0

AE
1
L
1
AE
1
L
1
+
AE
2
L
2

AE
2
L
2
0
AE
2
L
2
AE
2
L
2
3
7
7
7
7
5
8
>
>
>
>
<
>
>
>
>
:

1
1

1
2

1
3
9
>
>
>
>
=
>
>
>
>
;
=
L
2

K

33
8
>
>
>
>
<
>
>
>
>
:
1 +
2x
1
1
L
1 +
2x
2
1
L
1 +
2x
3
1
L
9
>
>
>
>
=
>
>
>
>
;
| {z }
fHg
31
(2.133)
Now substituting for

K

33
and simplifying the right hand side of Eq. 2.133, it
can be seen that
L
2

K

33
fHg
31
=
8
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
:
L
2
n
2AE
1
x
1
1
L
1
L

2AE
1
x
2
1
L
1
L
o
L
2
n

2AE
1
x
1
1
L
1
L
:::
+2(
AE
1
L
1
+
AE
2
L
2
):::
:::
x
2
1
L

2AE
2
x
3
1
L
2
L
o
L
2
n

2AE
2
x
2
1
L
2
L
+
2AE
2
x
3
1
L
2
L
o
9
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
;
53
Constraining the micro displacements,

1
1
=
1
3
= 0 (2.134)
Thus, Eq. 2.133 becomes,

AE
1
L
1
+
AE
2
L
2

1
2

=
L
2

2AE
1
x
1
1
L
1
L
+ 2

AE
1
L
1
+
AE
2
L
2

(2.135)
LetL
1
=L
2
=
L
2
andx
1
1
= 0; x
2
1
=
L
2
andx
3
1
=L. Substituing all these parameters
and simplifying Eq. 2.135 result in,

2AE
1
+ 2AE
2
L

1
2

=
L
2
(AE
1
AE
2
) (2.136)
thus,

1
2
=
L
2

E
1
E
2
E
1
+E
2

(2.137)
A special case is E
1
=E
2
= 0,
1
2
= 0.
Now introducing the required values and x coordinates (i.e. x
1
1
= 0; x
2
1
=
L
2
andx
3
1
=L) in Eq. 2.125 leads to,
[T ]
32
=
2
6
6
6
6
4
1 0
1
2
+
1
2
E
1
E
2
E
1
+E
2
1
2

1
2
E
1
E
2
E
1
+E
2
0 1
3
7
7
7
7
5
(2.138)
Now transforming the stiness matrix and letting

=
E
1
E
2
E
1
+E
2
,
[K] = [T ]
T
32

K

33
[T ]
32
(2.139)
54

K

33
[T ]
32
=
2
6
6
6
6
4
2AE
1
L

2AE
1
L
0

2AE
1
L
2AE
1
L
+
2AE
2
L

2AE
2
L
0
2AE
2
L
2AE
2
L
3
7
7
7
7
5
2
6
6
6
6
4
1 0
1
2
+
1
2

1
2

1
2

0 1
3
7
7
7
7
5
(2.140)
Let

=

AE
1
L
+
AE
2
L

,

K

33
[T ]
32
=
2
6
6
6
6
4
2AE
1
L

AE
1
L

AE
1
L

AE
1
L
+

AE
1
L

2AE
1
L
+

+

2AE
1
L
+

AE
2
L

AE
2
L

2AE
2
L

AE
2
L
+
AE
2
L
3
7
7
7
7
5
(2.141)
Now considering [T ]
T
32

K

33
[T ]
32
results in,

K

33
[T ]
32
=
2
6
4
1
1
2
+
1
2

0
0
1
2

1
2

1
3
7
5
2
6
6
6
6
4
2AE
1
L

AE
1
L

AE
1
L

AE
1
L
+

AE
1
L

2AE
1
L
+

+

2AE
1
L
+

AE
2
L

AE
2
L

2AE
2
L

AE
2
L
+
AE
2
L
3
7
7
7
7
5
(2.142)
Simplifying Eq. 2.142 becomes,

K

33
[T ]
32
=
2
6
4

2AE
1

L
+

2
+

2

2

+

2

2

AE
2
L

2

2

2

AE
2
L

2

2

+

2
+
2AE
2

L
3
7
5
(2.143)
A special case is E
1
=E
2
=E!

= 0 and

=
2AE
L
then,
[K]
[22]j
E1=E2=E
=
2
6
4
AE
L

AE
L

AE
L
AE
L
3
7
5
(2.144)
55
The homogenized stiness coecients are dened through the following procedure,
The strain energy of the original micro-structural material is given by,
U
total
=
1
2
fg
T
21
[K]
22
fg
21
(2.145)
where is the global displacement vector (anti-periodic eld)
fg =
8
>
<
>
:
1
1
9
>
=
>
;
(2.146)
Then the total strain energy is,
U
total
=
1
2
8
>
<
>
:
1
1
9
>
=
>
;
T
2
6
4
K K
K K
3
7
5
8
>
<
>
:
1
1
9
>
=
>
;
(2.147)
=
1
2
8
>
<
>
:
1
1
9
>
=
>
;
T
8
>
<
>
:
2K
2K
9
>
=
>
;
(2.148)
=
1
2
4K (2.149)
= 2K (2.150)
Where K =

2AE
1

L
+

2
+

2

.
Now computing the homogenized energy using E
Hom
and the associated stiness
matrix [K
Hom
]
22
dened as follows:
[K
Hom
]
22
=
2
6
4
K
Hom
K
Hom
K
Hom
K
Hom
3
7
5
(2.151)
56
Where K
Hom
=
AE
Hom
L
.
The strain energy of the homogenized structure is given by,
U
Hom
=
1
2
fg
T
21
[K
Hom
]
22
fg
21
(2.152)
=
1
2
8
>
<
>
:
1
1
9
>
=
>
;
T
2
6
4
K
Hom
K
Hom
K
Hom
K
Hom
3
7
5
8
>
<
>
:
1
1
9
>
=
>
;
(2.153)
= 2K
Hom
(2.154)
According to the eective properties of multi-scale [12], the strain energies must be
the same to attain the equilibrium,
U
Hom
=U
total
(2.155)
or else,
2K
Hom
= 2K
total
= 2K (2.156)
or,
2AE
Hom
L
= 2K (2.157)
and,
E
Hom
=
L
A
K (2.158)
Substituting for K,

and

. Further, simplifying,
E
Hom
=
L
A

2AE
1

L
+

2
+

2

2
!
(2.159)
E
Hom
= 2
E
1
E
2
E
1
+E
2
(2.160)
57
Special Cases
Consider a one-dimensional bar as shown in Figure 2.10 subjected to an axial load
P .
1 2
d d d
L
L
2
1
1
2 3
A
1
, E
1
A
2
, E
2
1 2 3
P
Figure 2.10: A One-dimensional Fixed End Macro-structure
Let A
1
, E
1
, L
1
and A
2
, E
2
, L
2
be the area, Young's modulus and length of micro
elements 1 and 2 respectively.

1;2;3
are the micro displacements. The total de
ection
of the prismatic bar is given by,
=
1
+
2
+
3
(2.161)

1
= 0 (2.162)

2
=
PL
1
AE
1
(2.163)

3
=
PL
2
AE
2
(2.164)
=
PL
1
AE
1
+
PL
2
AE
2
(2.165)
but L
1
=L
2
=
L
2
=
PL
2AE
1
+
PL
2AE
2
(2.166)
=
PL
A

E
1
+E
2
2E
1
E
2

| {z }
1
E
Hom
(2.167)
58
Thus,
E
Hom
= 2
E
1
E
2
E
1
+E
2
(2.168)
Case 1: One-dimensional prismatic bar with same materials
Consider E
1
= E
2
= E, homogeneous material. Substituting in 2.160 and
simplifying,
E
Hom
= 2
E
2
2E
(2.169)
E
Hom
=E From aperiodic method (2.170)
let E
1
=E
2
=E. Substituting in 2.168,
E
Hom
=E From analytical method (2.171)
Case 2: One-dimensional prismatic bar with two dierent materials
Consider E
1
= 1; E
2
= 2, homogeneous material. Substituting in 2.160 and
simplifying,
E
Hom
= 2
2
3
(2.172)
E
Hom
= 1:33 From aperiodic method (2.173)
let E
1
= 1; E
2
= 2. Substituting in 2.168,
E
Hom
= 1:33 From analytical method (2.174)
Thus, from the above two cases it is evident that the results produced by the
aperiodic method are in perfect correlation with the analytical method.
59
2.4.2 Case 2: Two-Dimensional Approach
In this case, heterogeneous problems in two-dimensions are compared. As a prelim-
inary test, aperiodic method is compared with the analytical material constants of
homogeneous media. Further, veried against the well-known literature of Kikuchi
et al. [4, 10] for three dierent material patterns.
2.4.2.1 Case 2a: Homogeneous and Isotropic Medium with Various
Mesh
In this case the problem with homogeneous media is considered to prove the
accuracy of aperiodic method against theoretical solutions. The Young's modulus
is E
f
= 10 for ber and E
m
= 10 for matrix, and Poisson ratio is zero and 0.3.
Geometry of the macro element are shown in Figure 2.11.
Material Pattern of Micro Elements
(a) A Homogeneous Macro-element
Material Pattern of Micro Elements

E=10.00 Nu=0.00
E=10.00 Nu=0.00
(b) A Macro-element with 10x10 Mesh
Figure 2.11: Homogeneous and Isotropic Media with Perfect Uniform Mesh
A rigorous test method is considered for materials of homogeneous medium
with arious types of mesh. First, uniform mesh as shown in Figure 2.11b. Second,
distorted by advancing front technique in Figure 2.11c. Lastly a highly distorted
60
Material Pattern of Micro Elements

E=10.00 Nu=0.00
E=10.00 Nu=0.00
(c) A Macro-element with Advanceing
Front Mesh
Material Pattern of Micro Elements

E=10.00 Nu=0.00
E=10.00 Nu=0.00
(d) A Macro-element with Highly Dis-
torted Mesh
Figure 2.11: Homogeneous and Isotropic Media with Highly Distorted Mesh
mesh represented in Figure 2.11d for convergence check. Dimension of the macro-
element is 10x10 units.
The analytical material matrix for a plane stress case is given by,
[E]
33
=
E
1
2
2
6
6
6
6
4
1 0
1 0
0 0
1
2
3
7
7
7
7
5
(2.175)
E denotes Young's modulus, is the Poisson's ratio.
Substituting E
f
=E
m
=E = 10 and
f
=
m
= 0:3, for two separate test cases,
the analytical and homogenised values are given in Table 2.2.
61
Poisson's Ratio Mesh Type E
11
E
12
E
22
E
33
= 0
Analytical (exact) 10 0 10 5
10x10 10 0 10 5
Advancing Front 10 0 10 5
Highly Distorted 10 0 10 5
= 0:3
Analytical (exact) 10.989 3.296 10.989 3.846
10x10 10.988 3.296 10.989 3.846
Advancing Front 10.988 3.296 10.988 3.846
Highly Distorted 10.988 3.296 10.988 3.846
Table 2.2: Comparison of Homogenized Values in Homogeneous and Isotropic
Medium
2.4.2.2 Case 2b: Checker Board Type Material Distribution
An exact checker board pattern of material distribution as given in Guedes and
Kikuchi [10, pp. 166] is considered as test problem in this case. The substance ma-
terial is characterized byE
1111
=E
2222
= 30 andE
1122
=E
1212
= 10. This example
was solved by Guedes and Kikuchi [4, pp. 166-168] by asymptotic homogenization
method. The unmeshed and meshed geometry is shown in Figure 2.12.
(a) A Un-meshed Checker Board Mate-
rial Pattern
Material Pattern of Micro Elements
(b) A 8x8 Checker Board Material Pat-
tern
Figure 2.12: Checker Board Material Distribution with Mesh
62
Further, the convergence plot of E
H
ijkl
against the dierent mesh size (2 2, 4 4,
6 6, 8 8, 10 10, 12 12, 16 16 and 100 100) is given in Figure 2.13. For
additional details reader could refer to Guedes and Kikuchi [10, pp. 168].
0 0.05 0.1 0.15 0.2 0.25
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Rate of Convergence of E
H
ijkl
mesh size H
1/E
H
ijkl

1/E
H
1111
1/E
H
1122
1/E
H
1212
Figure 2.13: Rate of Convergence of E
H
ijkl
vs. Mesh Size
The comparison of asymptotic and aperiodic method of homogenization is given
in the Table 2.3. It should to be noted that the aperiodic values are for 100x100
(H 0:005) mesh.
Mesh Type E
1111
E
1122
E
2222
E
1212
Asmptotically @ H=0 [10, pp. 167] 7.7172 3.1508 7.7172 3.2718
Aperiodic @ H 0:005 7.9905 3.0195 7.9905 3.6923
Table 2.3: Comparison of Homogenized Values for a Checker Board Pattern
63
2.4.2.3 Case 2c: A Rectangular Hole in Homogeneous and Isotropic
Material
In this case, there exists a rectangular hole (void) inside the material, the substance
material is characterized by E
1111
= E
2222
= 30 and E
1122
= E
1212
= 10. This
example was solved by Bendoe and Kikuchi [4, pp. 209] with their regular approach
and also by adaptive nite element method to demonstrate the capability of the
method to assess the voids in material distribution. The unmeshed and meshed
geometry is shown in Figure 2.14. Each nite element is of one unit in length and
width, and a total of 20x20 units.
Material Pattern of Micro Elements
(a) A Rectangular Hole in Homogeneous
and Isotropic Material
Material Pattern of Micro Elements

E=26.67 Nu=0.33
E=26.67 Nu=0.33
(b) A Rectangular Hole in Homogeneous
and Isotropic Material with 20x20 Mesh
Figure 2.14: A Rectangular Hole in Homogeneous and Isotropic Material with mesh
The homogenized and adaptive results obtained by Bendsoe and Kikuchi [4] and
aperiodic values are provided in the Table 2.4
64
Mesh Type E
1111
E
1122
E
2222
E
1212
Inital 13.015 3.241 17.552 2.785
1
st
Adaptive 12.910 3.178 17.473 2.714
2
nd
Adaptive 12.865 3.146 17.437 2.683
3
nd
Adaptive 12.844 3.131 17.421 2.668
Aperiodic 14.202 3.502 18.972 2.730
Table 2.4: Comparison of Homogenized Values for a Rectangular Hole
2.4.2.4 Case 2d: Soft and Hard Composite Material
A composite material which consists soft material with elastic modulus E
soft
= 10,
hard materialE
hard
= 1000 and Poisson's ratio of = 0:3 is considered for this case.
This example was solved by Bendoe and Kikuchi [4, pp. 207] using asymptotic
approach and also by implementing the adaptive nite element method. The
unmeshed and meshed geometry is shown in Figure 2.15. Each nite element is of
one unit in length and width, and a total of 16x16 units.
(a) A Soft and Hard Composite Macro-
element
Material Pattern of Micro Elements
(b) A Soft and Hard Composite Macro-
element with 16x16 Mesh
Figure 2.15: Soft and Hard Composite Material with Mesh
65
The homogenized and adaptive results obtained by Bendsoe and Kikuchi [4] and
aperiodic values are provided in the Table 2.5
Mesh Type E
1111
E
1122
E
2222
E
1212
16x16 149.80 71.61 149.80 87.12
1
st
Adaptive 127.12 62.91 127.12 75.90
2
nd
Adaptive 125.79 62.92 125.79 75.28
Aperiodic 165.62 69.12 165.62 65.40
Table 2.5: Comparison of Homogenized Values for Soft and Hard Material
From the four dierent cases of the two-dimensional problems, it is evident that
the homogenized elastic constants obtained from the aperiodic method is almost
undeviated from the analytical and literature result. Although, aperiodic results are
close to references, one should critically note that aperiodic method is very dierent
from the asymptotic periodic method. Until and unless the periodic boundary
conditions are exactly mimicked in the aperiodic method to compare against the
other methods, its dicult to expect the results that are similar. Therefore, this
explains the reason for some of the unmatched elastic constants in the above cases.
66
Chapter 3
The Finite Element Formulation
In this chapter, the procedures to implement the nite element method for the
aperiodic homogenization method are discussed in detail. Homogenization proce-
dures, although they could be solved through the analytical approach for simple
one-dimensional problems, it is always viable to develop computational procedures
and fast algorithms for obtaining the solution for various loading and boundary
conditions.
In order to eectively encounter the diculties of practical problems with various
specications, the aperiodic homogenization method as a nite element package is
implemented. Considering the fact that, the normal (
1
J
,
2
J
) and shear (
3
J
,
4
J
)
micro-solution variables plays an important role in developing a solution procedure
for aperiodic method. A detailed nite element procedure with the elaboration
on obtaining four micro-solutions variables is discussed. Additionally, to confront
the issues with the sub-patch extraction, discretization of the sub-patches and
assigning the materials properties for a micro-cell within a sub-patch, a digital image
based processes were implemented to exactly mesh and to extract the material
properties in each sub-micro element, this is discussed comprehensively in Chapter
4. A regular iso-parametric technique (for 4, 8, 9 noded elements) is employed
for the discretization of the global domain. Lastly, the generalized derivation of
micro-solution for aperiodic method in natural coordinate system and computer
implementation of the aperiodic homogenization, digital image based meshing and
main nite element program are discussed in the appendices A and B.
67
3.1 The Three-Dimensional Variational Problem
Most of the continuum problems, including the solid and structural mechanics
problems are formulated according to the two methods, namely dierential equation
method and variational method. For the purpose of the nite element implementa-
tion of the aperiodic homogenization, the variational method is used to obtain the
quantitative data such as displacements and stresses.
Consider the general problem formulation in three-dimension, although it is
implemented for a two-dimensional models in this work. The potential energy of
an elastic body
p
is dened by,

p
=W
p
(3.1)
Where is the strain energy stored in the system, and W
p
is the work done on
the body due to external forces. According to the prinicple of minimum potential
energy,

p
(u; v; ;w) =(u; v; ;w)W
p
(u; v; ;w) (3.2)
Where u; v; ;w are the displacements in Cartesian co-ordinate system. The strain
energy of the system is given by,
=
1
2
ZZZ
V
~
T
~ dV (3.3)
Where V is the volume of the body, ~ is the strain state and ~ is the stress state.
By using the denition of stress-strain and the work done by the external work we
get,
68
=
1
2
ZZZ
V
~
T
[E] ~ dV
ZZZ
V
~
T
[E] ~
0
dV (3.4)
W
p
=
ZZZ
V
~
f
T
b
~
U dV +
ZZ
S
~
f
t
~
U dS (3.5)
Where ~ and ~
0
are the regular and initial strains, [E] is the material stiness
matrix,
~
f
b
is the body force vector,
~
f
t
is the traction vector and
~
U is the nodal
displacement vector.
Step 1: The solid body or area is divided into n nite elements
Step 2: The displacement model within an element n

is assumed as

U =
8
>
>
>
>
<
>
>
>
>
:
u(x; y; z)
v(x; y; z)
w(x; y; z)
9
>
>
>
>
=
>
>
>
>
;
= [N]
~
Q
n

(3.6)
Where
~
Q
n

is the vector of nodal displacements degrees of freedom of the
element, and [N] is the matrix of shape functions.
Step 3: The stiness matrices and load vectors are to be derived from the principle
of minimum potential energy. The potential energy functional of the body
can be written as,

p
=
n
X
i=1

n
p
(3.7)
Where
n
p
is the potential energy of element n

, and given by,
69

p
n

=
1
2
ZZZ
V
n

~
T
[E] (~ 2~
0
) dV
ZZZ
V
n

~
f
T
b
~
U dV
ZZ
S
n

~
f
t
~
U dS (3.8)
Where V
n

is the volume of the element,
~
f
b
is the vector of body forces per
unit volume and S
n

is the portion of the surface of the element over which
forces and surface tractions are prescribed.
The strain vector ~ in Eq. 3.8 can be expressed in terms of the nodal
displacement vector
~
Q
n

by dierentiating Eq. 3.6 and given by,
~ =
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:

xx

yy

zz

xy

yx

zx
9
>
>
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
>
>
;
=
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
@u
@x
@v
@y
@w
@z
@u
@y
+
@v
@x
@v
@z
+
@w
@y
@w
@x
+
@u
@z
9
>
>
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
>
>
;
~ =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
@
@x
0 0
0
@
@y
0
0 0
@
@z
@
@y
@
@x
0
0
@
@z
@
@y
@
@z
0
@
@x
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
8
>
>
>
>
<
>
>
>
>
:
u
v
w
9
>
>
>
>
=
>
>
>
>
;
= [B]
~
Q
n

(3.9)
70
Where,
[B] =
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
@
@x
0 0
0
@
@y
0
0 0
@
@z
@
@y
@
@x
0
0
@
@z
@
@y
@
@z
0
@
@x
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
[N] (3.10)
The stresses ~ can be obtained from the strains ~ using the following equation,
~ = [E] (~ ~
0
)
= [E] [B]
~
Q
n

[D] ~
0
(3.11)
Now substituting the Eqs. 3.6 and 3.9 into Eq. 3.8 the potential energy of an
element takes the form of,

p
n

=
1
2
ZZZ
V
n

~
Q
n
T
[B]
T
[E] [B]
~
Q
n

dV
ZZZ
V
n

~
Q
n
T
[B]
T
[E] ~
0
dV
ZZ
S
n

~
Q
n
T
[N]
T
~
f
t
dS
ZZZ
V
n

~
Q
n
T
[N]
T
~
f
b
dV (3.12)
In Eqs. 3.8 and 3.12, only the body and surface tractions are considered.
However, in some cases the external concentrated forces will be acting on
some of the nodes. If
~
F denotes the vector of nodal forces acting in the
directions of the nodal displacement vector
~
Q of the total structure or body.
Then the total potential energy of the structure or body can be written as,
71

p
=
n
X
i=1

n

p

~
Q
T
~
F (3.13)
Where
~
Q =
8
>
>
>
>
>
>
>
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
>
>
>
>
>
>
>
:
Q
1
Q
2
Q
3
Q
4
.
.
.
Q
m
9
>
>
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
>
>
;
is the vector of nodal displacements of the entire
structure or body and m is the total number of nodal displacements or
degrees of freedom. The summation of Eq. 3.13 implies the expansion of
element matrices to structure or body size followed by the summation of
overlapping terms. Thus Eqs. 3.12 and 3.12 gives

p
=
1
2
~
Q
T
2
4
n
X
i=1
ZZZ
V
n

[B]
T
[E] [B] dV
3
5 ~
Q
~
Q
T
n
X
i=1
0
@
ZZZ
V
n

[B]
T
[E] ~
0
dV
+
ZZ
S
n

[N]
T
~
f
t
dS
ZZZ
V
n

[N]
T
~
f
b
dV
1
A

~
Q
T
~
F (3.14)
Eq. 3.14 expresses the total potential energy of the structure or body in
terms of the nodal degrees of freedom
~
Q The static equilibrium conguration
of the structure can be found by solving the following necessary conditions of
the minimum potential energy.
@
p
@
~
Q
=
~
0 for all nodal displacements (3.15)
72
From the Eq. 3.14, Eq. 3.15 can be expressed as,
2
6
6
6
6
6
4
n
X
i=1
ZZZ
V
n

[B]
T
[E] [B] dV
| {z }
Element stiness matrix
3
7
7
7
7
7
5
| {z }
Global or overall stiness matrix [K]
~
Q
|{z}
Global vector of displacements
=
~
F
c
+
n
X
i=1
0
@
ZZZ
V
n

[B]
T
[E] ~
0
dV +
ZZ
S
n

[N]
T
~
f
t
dS
ZZZ
V
n

[N]
T
~
f
b
dV
1
A
| {z }
Total vector of nodal forces
~
F
(3.16)
Denoted as,
"
n
X
i=1

K
n

#
~
Q =
~
F
c
+
n
X
i=1

~
F
n

ini
+
~
F
n

s
+
~
F
n

b

=
~
F (3.17)
Where,

K
n

=
ZZZ
V
n

[B]
T
[E] [B] dV
! Element stiness matrix (3.18)
~
F
c
=
ZZZ
V
n

[B]
T
[E] ~
0
dV
! Element load vector due to concentrated
forces on the nodes (3.19)
73
~
F
n

ini
=
ZZZ
V
n

[B]
T
[E] ~
0
dV
! Element load vector due to initial strains (3.20)
~
F
n

s
=
ZZ
S
n

[N]
T
~
f
t
dS
! Element load vector due to surface tractions (3.21)
~
F
n

b
=
ZZZ
V
n

[N]
T
~
f
b
dV
! Element load vector due to body forces (3.22)
Step 4: The equilibrium of the overall structure or body can be expressed by Eq.
3.17
[K]
~
Q =
~
F (3.23)
where,
[K] =
n
X
i=1
= Assembled/Global stiness matrix (3.24)
and
~
F =
~
F
c
+
n
X
i=1
~
F
n

ini
+
n
X
i=1
~
F
n

s
+
n
X
i=1
~
F
n

b
= Assembled/Global nodal load vector (3.25)
Step 5: The stresses and other quantitative data of interests are obtained at the
time of post-processing. It should be noted that the stiness matrices [K]
and

K
n

are symmetry and cannot be inversed if they are singular.
74
3.2 Formulation of The Multi-Scale Model
3.2.1 Discretization Procedure for the Multi-Scale Model
Let

e
be the domain of a typical nite element in
, and
e
be the domain of a
typical sub-patch or unit cell inside

e
as shown in Figure 3.1.
W
We
Global/Macro Finite Element
Y
Integration Point
(Gauss Quadrature
Point
G
e
Micro Finite Element /
Sub-Micro Finite Element
Sub-Patch / Unit cell /
RVE
X
X
X
y
y
1
y
2
2
X
1
Point 2x2)
Figure 3.1: Discretization of a Typical Sub-Patch
Letx represent a position vector associated with a point in

e
,x
i
; i = 1; 2 represent
the components of position vector ~ x. Similarly, y
i
; i = 1; 2 be the components
of position vector ~ y. It is assumed that the macro element around the Gauss
75
quadrature point could be contiguous. In other words, there is a possibility of one
macro element covering a nite element, but not extending over the nearby global
nite element.
By the use of this geometric decomposition, a spectrum of simulations can be
developed. In two-scale simulations, the sub-patch
e
could be employed as a micro-
scale domain used for the computation of the homogenized material constitutive
parameters for the use in macro models on
.
In single scale models, it is assumed that the same coordinate system is used on
the macro and micro domains x =y. In this work to demonstrate the workability
of the analysis, the models are of two scale representations.
The discretization of the sub-patch
e
is achieved by dening an appropriate
C
0
or C
1
nite element model capable of representing the material distribution on
the sub-patch and distinguishing all material boundaries, in this work a 4 noded
bi-linear element is employed to represent the discretized sub-patch from the global
nite element.
The micro-mesh must be generated within each sub-patches that are extracted
around nn Gauss quadrature points of a global nite element. Figure 3.2 shows
a typical one-dimensional and two-dimensional micro-mesh. The outer boundary
of the generated micro-mesh must be aligned with the outer boundary of the
sub-patch (regular polygon such as square or rectangle). Further, there must exist
a micro-node at the location of each macro-node in the sub-patch.
The position of micro-mesh node J in the direction i is denoted as Y
J
i
. The
micro displacement variable at micro node J in direction iis denoted by

U
J
i
. The
micro solution variable is analogous to the cell solution variables in the clasical
76
U
1
U
1
U
2 c
1
1 2
Macro Node Micro Nodes
(a) A Typical One-Dimensional Micro-Mesh
U
1
U
1
U
3
c
1
1
2
Macro Node Micro Nodes
3 4
U
2
U
4
(b) A Typical Two-Dimensional Micro-Mesh
Figure 3.2: Multi-Dimensional Micro Meshes
homogenization method, and is dened as
K
J
, K = 1; 4, where K is the micro-
function number associated with the anti-periodic modes and J is the micro node
number. The micro solution variables must be restrained on the sub-patch boundary
and at the macro-node positions within the sub-patch. Since the elements in micro-
mesh have arbitrary material properties, they must also have arbitrary stiness
matrices.
3.2.1.1 Iso-Parametric Meshing
All the commercial nite element packages are equipped with a mesh generator
for discretizing the domains and their sub-domains. In addition to the image
based meshing system (Appendix B.2), it is a necessary to develop a mini mesh
generator for the computer program Homo.Aperiodic to discretize the global domain.
Therefore, an automated iso-parametric based mesh generator was implemented in
this work for dividing the global domain.
77
Consider a structure or arbitrary shaped body in physical and natural co-
ordinate system as shown in the Figure 3.3
X
Y
1
2
1
2
Physical Plane / Parent Plane /
Natural Co-ordinate Cartesian Co-ordinate
x
h
Figure 3.3: Physical and Parent Plane
In the natural co-ordinate systems each global domain (1 and 2) are represented as a
square, eventhough it is distorted in the physical co-ordinate system. The following
steps describes the process of mesh generaration by iso-parametric method.
X
Y
1
Physical Plane / Parent Plane /
Natural Co-ordinate Cartesian Co-ordinate
x
h
1
X
1
,Y
1 X
2
,Y
2
X
3
,Y
3
X
4
,Y
4
1 2
3
4
x=-1
h=-1
x=-0.5
h=-0.5
x=+1
h=-1
x=0
h=0
x=+1
h=0
x=-1
h=0
x=0
h=+1
x=+1
h=+1
x=-1
h=+1
2
3 4
2
4 3
Figure 3.4: An Example of Mapping Physical and Parent Plane
78
Step 1: Divide the natural coordinates according to the number of elements needed.
For example if 2x2 elements are required in the physical system, then the
natural co-ordiantes of the divided elements in parent plane follows the Figure
3.4
Step 2: From the divided natural coordinates the physical coordinates are mapped
by using the interpolation Eqs. 3.26. During the mapping process each shape
functions corresponding to the element (4, 8, and 9 noded) are evaluated
at each natural coordinates of the vertex in parent plane. Evaluating the
interpolation equations at each natural co-ordinate vertex fetches the physical
coordinates within the pre-dened global vertices values (X
1
; Y
1
::: X
n
; Y
n
).
1
2
3 4
(a) A Typical 4-Noded
Bi-Linear Element
1 2 3
4
5 6 7
8
(b) A Typical 8-Noded
Serendipity Element
1 2 3
4
5 6 7
8 9
(c) A Typical 9-Noded
Lagrangian Element
Figure 3.5: Types of Finite Elements
The mapping equation follows as below,
X
i
=
1
(
i
;
i
)X
1
+
2
(
i
;
i
)X
2
+
3
(
i
;
i
)X
3
+:::
n
(
i
;
i
)X
n
Y
i
=
1
(
i
;
i
)Y
1
+
2
(
i
;
i
)Y
2
+
3
(
i
;
i
)Y
3
+:::
n
(
i
;
i
)Y
n
(3.26)
Where i represents the vertices in parent plane, (;
i

i
) are the natural co-
ordinate values at each node in parent plane.
1
;
2
; :::
n
are the shape
functions corresponding to the element.
79
For 4 noded element
1
;
2
;
3
; and
4
are the shape functions at nodes 1-4.

1
(
i
;
i
) =
1
4
(1
i
)(1
i
)

2
(
i
;
i
) =
1
4
(1 +
i
)(1
i
)

3
(
i
;
i
) =
1
4
(1 +
i
)(1 +
i
)

4
(
i
;
i
) =
1
4
(1
i
)(1 +
i
)
(3.27)
Similarly, for 8 noded element
1
;
2
;
3
; :::
8
are the shape functions at
nodes 1-8.

k
(
i
;
i
)=
1
4
(1 +
i

k
)(1 +
i

k
)(
i

k
+
i

k
1) k= 1; 3; 5; 7

k
(
i
;
i
)=
1
2
(1
2
i
)(1 +
i

k
) k= 2; 6

k
(
i
;
i
)=
1
2
(1
2
i
)(1 +
i

k
) k= 4; 8
(3.28)
Lastly, for 9 noded element
1
;
2
;
3
; :::
9
are the shape functions at nodes
1-9.

k
(
i
;
i
)=
1
4
(1 +
i

k
)(1 +
i

k
)(
i

k
)(
i

k
) k= 1; 3; 5; 7

k
(
i
;
i
)=
1
2
(1
2
k
)(1 +
i

k
)(
i

k
) k= 2; 6

k
(
i
;
i
)=
1
2
(1
2
i
)(1 +
i

k
)(
i

k
) k= 4; 8

k
(
i
;
i
)= (1
2
k
)(1
2
i
) k= 9
(3.29)
Step 3: All the above steps are repeated tilli is the last node in natural coordinate
system.
80
3.2.2 Relationship of Micro and Macro Solutions
The multi-scale displacement approximation for the displacement is similar to the
one derived from the conventional homogenization theory, and it follows as,
U
i
=

J

U
J
i
(3.30)
where

j
denotes the global/local shape function associated with micro-node J, u
J
i
represents the micro displacement of J
th
node.
In order to obtain a multi-scale relationship, it is important to ensure a strong
link between the macroscopic and microscopic variables. This can be achieved by
dening the independent macro and micro displacements and linking them through
explicit constraints. By this approach several objectives could be achieved.
First, it is necessary to insure that the total micro-scale correctional displace-
ment increment is formed by a multiplicative process, involving products of
macro displacements and micro correction terms that, in combination, mimics
the role of cell solutions and macro displacements in the homogenization
method.
Second, by the multiplicative process, it can be assured that the resulting
nite element model has the proper rigid body modes, even when micro
solution variables are included in the model
Further, it should be noted that, in the multiplicative model, the micro-scale
variables need to be computed only once with the formulation automatically scaling
up the micro-scale corrections i.e.
K
J
functions, and this property can be used to
dene the averaged or homogenized elastic constants independent of the macro
displacements values.
81
The constraint equations, involving the discrete displacements at micro-node J
in the direction i are dened through the following micro-macro coupling equation:

U
J
i
=T
JKi
U
K
i
(3.31)
Where T
JKi
is the transformation matrix between the macro and micro displace-
ments, J! micro-nodes and K! macro-nodes (refer Figure 3.6). It should be
noted that the summation on i on the right hand side of Eq. 3.31 is not considered.
1 2
3 4
L
1
K=
K= K=
L
2
K=
J=1 J=2 J=3
Figure 3.6: Micro-Macro Relationship in a 4-Noded Element
The transformation matrix components, associated with the displacements in the
x
1
coordinate direction are dened in terms of the micro solution variable
K
J
at
each micro node and it follows as:
T
J11
=

1
Y
J
1
L
1
+

1
J
L
1

1
Y
J
2
L
2
+

4
J
L
2

(3.32)
T
J21
=

Y
J
1
L
1

1
J
L
1

1
Y
J
2
L
2
+

4
J
L
2

(3.33)
T
J31
=

Y
J
1
L
1

1
J
L
1

Y
J
2
L
2

4
J
L
2

(3.34)
T
J41
=

1
Y
J
1
L
1
+

1
J
L
1

Y
J
2
L
2

4
J
L
2

(3.35)
82
In the similar argument, the transformation matrix components associated with the
displacements in x
2
coordinate direction are dened in terms of the micro solution
variable
K
J
at each micro node and it follows as:
T
J12
=

1
Y
J
1
L
1
+

3
J
L
1

1
Y
J
2
L
2
+

2
J
L
2

(3.36)
T
J22
=

Y
J
1
L
1

3
J
L
1

1
Y
J
2
L
2
+

2
J
L
2

(3.37)
T
J32
=

Y
J
1
L
1

3
J
L
1

Y
J
2
L
2

2
J
L
2

(3.38)
T
J42
=

1
Y
J
1
L
1
+

3
J
L
1

Y
J
2
L
2

2
J
L
2

(3.39)
The above transformation equations are the shape functions dened to capture the
spatial variations (micro displacements) within the sub-patch or unit cell due to
the in
uence of constant strain eld.
In order to dene the micro solution variables
K
J
, the micro level approximation
must be used. From the variational problem dened before after the application of
weighting function,

K
LJik

U
J
k
=

F
Li
(3.40)
where,

K
LJik
=
Z

tE
ijkl

J;l

L;j
dV (3.41)
and,

F
Li
=
Z

tf
i

L
dV (3.42)
83
3.3 Solution of Micro Solution Variables
Micro solution variable
K
J
, K = 1; 4 can be computed using the Eq. 3.40 and the
transformation Eqs. 3.32 through 3.39. For this work it is assumed that there is no
external forces acting on the sub-patch, except that the anti-periodic displacements
are applied to introduce the disturbances in the spatial domain of a sub-patch.
Therefore if the force term is zero in Eq. 3.40. Further, substituting Eq. 3.31
results,

K
LJik
T
JKi
U
K
i
= 0 (3.43)
The four micro solution variables
K
J
, K = 1; 4 can be computed from 3.43
by assuming macro deformation elds U
K
i
, induces the necessary macroscopic
strain elds. These strain elds are responsible for uncoupling the macroscopic
and microscopic problems in order to allow the computation of the microscopic
solutions. In terms of macro deformation elds, the most likely possibilities are the
periodic and aperiodic elds as shown in Figure 3.7
U U
1 2
Macro Node Micro Nodes
(a) A Typical One-Dimensional Periodic
Patch
-U U
1 2
Macro Node Micro Nodes
(b) A Typical One-Dimensional Aperiodic
Patch
Figure 3.7: Types of Strain Fields
If the periodic displacement eld is applied to the element associated with the Eq.
3.31, the equation is satised exactly for arbitrary values of the micro solution
variables
K
J
, if and only if it satises the necessary and sucient rigid body
84
conditions without strain. Therefore, the compatibility of the shape functions T
JKi
goes well with the assumption that these functions are demonstrating the proper
rigid body modes. Thus, it becomes obvious that the anti-periodic displacements
properly represents the mechanism to produce the micro solution variables.
3.3.1 Anti-Periodic (Deformation) Modes
Deformation modes imposes the boundary conditions on a sub-patch. For a two-
dimensional problem there are four deformation modes, two in the normal directions
and the other two in the shear directions as shown in Figure 3.8.
U = -1
U = -1 U = +1
U = +1
(a) Deformation Mode u
11
in x Direction
V = -19 = -1
9 = +1 V = +1
(b) Deformation Mode u
22
in y Direction
Figure 3.8: Normal Deformation Modes of a Typical Sub-Patch
85
V = -1
V = -1
V = +1
V = +1
(c) Deformation Mode u
21
in yx Direction
U = +1
U = -1 U = -1
U = +1
(d) Deformation Mode u
12
in xy Direction
Figure 3.8: Shear Deformation Modes of a Typical Sub-Patch
In-fact these modes are responsible for the macro nodal displacement of the
4 noded bi-linear element, which in-turn implemented to obtain the functions.
The deformation modes are used in the nite element calculation of eective
stiness coecients of a composite by the direct average method [27]. Further, the
properties of these modal patterns are such that they produce unit values of strains
(
11
;
22
;
12
and
21
) for a unit value of U or V's.
86
3.3.2 Computation of the Basic Axial Micro Solution Vari-
ables
Consider the macro-deformation eld associated with the constant strain
11
. Refer
Figure 3.8a. The strain eld is obtained by specifying the following nodal macro
displacements based on a set of nodal axial displacement U follows as:
U
1
1
=U
U
2
1
= +U
U
3
1
= +U
U
4
1
=U
U
K
2
= 0 K = 1; 4
(3.44)
For simplicity, letU = 1, since anti-displacements are scalar values, it is immaterial
if it is 1 or 1000. It has to be noted that the magnitude of anti-periodic displacements
will not aect the micro solutions. From Eq. 3.43,

K
LJik
T
JKi
U
K
i
= 0
let,i = 1; k = 1
For the computation in natural coordinate system, let

1
(Y
J
1
) = 1
Y
J
1
L
1
(3.45)

2
(Y
J
2
) = 1
Y
J
2
L
2
(3.46)
87
and the transformation of natural to physical coordinate is given by,
Y
1
=
L
1
2
(1 +
1
)
1
! (3.47)
Y
2
=
L
1
2
(1 +
2
)
2
! (3.48)
Substituting the Eqs. 3.47 and 3.48 in shape functions 3.45 and 3.46, further
simplifying leads to,

1
(Y
J
1
) =
1
2
(1
J
1
) (3.49)

1
(Y
J
2
) =
1
2
(1
J
2
) (3.50)
now considering the T
JKi
U
K
i
part and expanding on K,
T
JK1
U
K
1
=T
J11
U
1
1
+T
J21
U
2
1
+T
J31
U
3
1
+T
J41
U
4
1
(3.51)
Let
1
J
=

1
J
L
1
and
4
J
=

4
J
L
2
Further introducing the anti-periodic constants U =1 gives,
T
JK1
U
K
1
=T
J11

1
U
1
1
+T
J21

1
U
2
1
+T
J31

1
U
3
1
+T
J41

1
U
4
1
(3.52)
Substituting the Eqs. 3.49 and 3.50 in 3.32 through 3.35
88
=

1
2

1
J
1

+
1
J

1
2

1
J
2

+
4
J

+

1
2

1 +
J
1

1
J

1
2

1
J
2

+
4
J

+

1
2

1 +
J
1

1
J

1
2

1 +
J
2

4
J

1
2

1
J
1

+
1
J

1
2

1 +
J
2

4
J

(3.53)
=

1
2

1
J
1

1
2

1
J
2

+
1
2

1
J
1

4
J
+
1
2

1
J
2

1
J
+

1
J

4
J

+

1
2

1 +
J
1

1
2

1
J
2

+
1
2

1 +
J
1

4
J

1
2

1
J
2

1
J

1
J

4
J

+

1
2

1 +
J
1

1
2

1 +
J
2

1
2

1 +
J
1

4
J

1
2

1 +
J
2

1
J
+
1
J

4
J

1
2

1
J
1

1
2

1 +
J
2

1
2

1
J
1

4
J
+
1
2

1 +
J
2

1
J

1
J

4
J

(3.54)
=

1
2

1
J
1

+

1
2

1 +
J
1

1
2

1 +
J
1

+

1
2

1 +
J
1

4
J
+

1
2

1
J
2

1
2

1
J
2

1
2

1 +
J
2

1
2

1 +
J
2

1
J
+

1
4

1 +
J
1

1
J
2

1
4

1
J
1

1
J
2

+
1
4

1 +
J
1

1 +
J
2

1
4

1
J
1

1 +
J
2

(3.55)
=

1
4

1 +
J
1

J
2

J
1

J
2

1
4

1
J
1

J
2
+
J
1

J
2

+
1
4

1 +
J
1
+
J
2
+
J
1

J
2

1
4

1 +
J
2

J
1

J
1

J
2

(1
J
2
) (1 +
J
2
)

1
J
(3.56)
89
=
h
1 +

J
2
1

J
2
i

1
J
+

1
4

1

J
2
+

J
2
+
J
1
+
J
1

1

J
1

J
2

+

1
4

J
1

J
2
+

J
1

J
2
+

J
1

J
2

J
2
+

J
2
+
J
1
+
J
1
+

1

1

(3.57)
T
JK1
U
K
1
=2
1
J
+
J
1
(3.58)

K
LJ11

2
1
J
+
J
1

= 0 (3.59)

K
LJ11

2
1
J

=

K
LJ11

J
1

| {z }

T
1
J
(3.60)
But,

J
1
=
2Y
J
1
L
1
1 (3.61)
Thus,

1
J
=
1
2

K
LJ11

1

T
1
J
(3.62)
or

1
J
=
1
2

K
LJ11

1

K
LJ11

2Y
J
1
L
1
1

(3.63)
Substituting for
1
J
and rearranging,

1
J
=
L
1
2

K
LJ11

1

K
LJ11

2Y
J
1
L
1
1

(3.64)
or

1
J
=
L
1
2

K
LJ11

1

T
1
J
(3.65)
Eq. 3.65 represents the uncoupled form of the micro solution
1
J
from macro
solution.
90
Consider the macro-deformation eld associated with the constant strain
22
. Refer
Figure 3.8b. The strain eld is obtained by specifying the following nodal macro
displacements based on a set of nodal axial displacement V follows as:
U
1
2
=V
U
2
2
=V
U
3
2
= +V
U
4
2
= +V
U
K
1
= 0 K = 1; 4
(3.66)
Again, to ease out the calculations, let V = 1. From Eq. 3.43,

K
LJik
T
JKi
U
K
i
= 0
let,i = 2; k = 2
Now considering the T
JKi
U
K
i
part and expanding on K,
T
JK2
U
K
2
=T
J12
U
1
2
+T
J22
U
2
2
+T
J32
U
3
2
+T
J42
U
4
2
(3.67)
Let
3
J
=

3
J
L
1
and
2
J
=

2
J
L
2
Further introducing the anti-periodic constants V =1 gives,
T
JK
U
K
2
=T
J12

1
U
1
2
+T
J22

1
U
2
2
+T
J32

1
U
3
2
+T
J42

1
U
4
2
(3.68)
Substituting the Eqs. 3.49 and 3.50 in 3.36 through 3.39,
91
=

1
2

1
J
1

+
3
J

1
2

1
J
2

+
2
J

1
2

1 +
J
1

3
J

1
2

1
J
2

+
2
J

+

1
2

1
J
1

+
3
J

1
2

1
J
2

+
2
J

+

1
2

1 +
J
1

3
J

1
2

1
J
2

+
2
J

(3.69)
T
JK2
U
K
2
=2
2
J
+
J
2
(3.70)

K
LJ22

2
2
J
+
J
2

= 0 (3.71)

K
LJ22

2
2
J

=

K
LJ22

J
2

| {z }

T
2
J
(3.72)
But,

J
2
=
2Y
J
2
L
2
1 (3.73)
Therefore,

2
J
=
1
2

K
LJ22

1

T
2
J
(3.74)
or

1
J
=
1
2

K
LJ22

1

K
LJ22

2Y
J
2
L
2
1

(3.75)
Substituting for
2
J
and rearranging,

2
J
=
L
2
2

K
LJ22

1

K
LJ22

2Y
J
2
L
2
1

(3.76)
or

2
J
=
L
2
2

K
LJ22

1

T
2
J
(3.77)
Eq. 3.77 represents the uncoupled form of the micro solution
2
J
from macro
solution.
92
3.3.3 Computation of the Basic Shear Micro Solution Vari-
ables
Consider the macro-deformation eld associated with the constant strain
21
. Refer
Figure 3.8c. The strain eld is obtained by specifying the following nodal macro
displacements based on a set of nodal axial displacement V follows as:
U
1
2
=V
U
2
2
= +V
U
3
2
= +V
U
4
2
=V
U
K
1
= 0 K = 1; 4
(3.78)
Let V = 1. From Eq. 3.43,

K
LJik
T
JKi
U
K
i
= 0
let,i = 2; k = 1
Now considering the T
JKi
U
K
i
part and expanding on K,
T
JK2
U
K
2
=T
J12
U
1
2
+T
J22
U
2
2
+T
J32
U
3
2
+T
J42
U
4
2
(3.79)
Let
3
J
=

3
J
L
1
and
2
J
=

2
J
L
2
Further introducing the anti-periodic constants V =1 gives,
T
JK
U
K
2
=T
J12

1
U
1
2
+T
J22

1
U
2
2
+T
J32

1
U
3
2
+T
J42

1
U
4
2
(3.80)
93
Substituting the Eqs. 3.49 and 3.50 in 3.36 through 3.39,
=

1
2

1
J
1

+
3
J

1
2

1
J
2

+
2
J

+

1
2

1 +
J
1

3
J

1
2

1
J
2

+
2
J

+

1
2

1 +
J
1

3
J

1
2

1 +
J
2

2
J

1
2

1
J
1

+
3
J

1
2

1 +
J
2

2
J

(3.81)
T
JK2
U
K
2
=2
3
J
+
J
1
(3.82)

K
LJ21

2
2
J
+
J
1

= 0 (3.83)

K
LJ21

2
3
J

=

K
LJ21

J
1

| {z }

T
3
J
(3.84)
But,

J
1
=
2Y
J
1
L
1
1 (3.85)
Thus,

3
J
=
1
2

K
LJ21

1

T
3
J
(3.86)
or

3
J
=
1
2

K
LJ21

1

K
LJ21

2Y
J
1
L
1
1

(3.87)
Substituting for
3
J
and rearranging,

3
J
=
L
1
2

K
LJ21

1

K
LJ21

2Y
J
1
L
1
1

(3.88)
or

3
J
=
L
1
2

K
LJ21

1

T
3
J
(3.89)
Eq. 3.89 represents the uncoupled form of the micro solution
3
J
from macro
solution.
94
Lastly, consider the macro-deformation eld associated with the constant strain
12
.
Refer Figure 3.8d. The strain eld is obtained by specifying the following nodal
macro displacements based on a set of nodal axial displacement U follows as:
U
1
1
=U
U
2
1
=U
U
3
1
= +U
U
4
1
= +U
U
K
2
= 0 K = 1; 4
(3.90)
Let U = 1. From Eq. 3.43,

K
LJik
T
JKi
U
K
i
= 0
let,i = 1; k = 2
Now considering the T
JKi
U
K
i
part and expanding on K,
T
JK1
U
K
1
=T
J11
U
1
1
+T
J21
U
2
1
+T
J31
U
3
1
+T
J41
U
4
1
(3.91)
Let
1
J
=

1
J
L
1
and
4
J
=

4
J
L
2
Further introducing the anti-periodic constants V =1 gives,
T
JK
U
K
1
=T
J11

1
U
1
1
+T
J21

1
U
2
1
+T
J31

1
U
3
1
+T
J41

1
U
4
1
(3.92)
Substituting the Eqs. 3.49 and 3.50 in 3.32 through 3.35,
95
=

1
2

1
J
1

+
1
J

1
2

1
J
2

+
4
J

1
2

1
J
1

1
J

1
2

1
J
2

+
4
J

+

1
2

1 +
J
1

1
J

1
2

1 +
J
2

4
J

+

1
2

1
J
1

+
1
J

1
2

1 +
J
2

4
J

(3.93)
T
JK1
U
K
1
=2
4
J
+
J
2
(3.94)

K
LJ12

2
4
J
+
J
2

= 0 (3.95)

K
LJ12

2
4
J

=

K
LJ12

J
2

| {z }

T
4
J
(3.96)
But,

J
2
=
2Y
J
2
L
2
1 (3.97)
Thus,

4
J
=
1
2

K
LJ12

1

T
4
J
(3.98)
or

4
J
=
1
2

K
LJ12

1

K
LJ12

2Y
J
2
L
2
1

(3.99)
Substituting for
4
J
and rearranging,

4
J
=
L
2
2

K
LJ12

1

K
LJ12

2Y
J
2
L
2
1

(3.100)
or

4
J
=
L
2
2

K
LJ12

1

T
4
J
(3.101)
Eq. 3.101 represents the uncoupled form of the micro solution
4
J
from macro
solution.
96
3.3.4 Signicance of Micro Solution Variables
In homogeneous materials the strain eld are consistent with the material properties
and their spatial distribution. This will be dierent in the case of heterogeneous
materials. For example, consider an one-dimensional rod as shown in Figure
3.9, the rod has four nite elements. Supposedly, if the rod has same mate-
rials (E
1
=E
2
=E
3
=E
4
) in all four elements, then the strains are consistent
(
1
=
2
=
3
=
4
) for a unit displacements in the two ends with opposite nature.
On the second hand, if the materials (E
1
6=E
2
6=E
3
6=E
4
) are dierent in each
element, then the strains have to be dierent (
1
6=
2
6=
3
6=
4
) because of the
heterogeneity, else this would leads to the inconsistency in the strain denition on
each node.
1 2 3 4
U = -1 U = 1
L L L L
Figure 3.9: A One-Dimensional Rod with Four Finite Elements
In order to include the eects of heterogeneity in any given sub-patch, a special
function, namely function are dened for obtaining the micro-scale displacements
and its variation through the spatial coordinates.
=
du
dx
(3.102)
For homogeneous
=
du
dx
+
d
dx
(3.103)
For non-homogeneous
97
3.4 Macro Equilibrium Equations
Given the micro solution variables
K
J
, the macro level nite element problem can
be formulated by introducing the micro-macro coupling equation into the micro
equilibrium equations as follows:

U
J
i
=T
JKi
U
K
i
micro-macro coupling equation

K
LJik

U
J
k
=

F
Li
micro equilibrium equation

K
LJik
T
JKi
U
K
i
=

F
Li
(3.104)
Pre-multiplying 3.104 by the transpose of the transformation matrix T
JKi
K
MKik
=T
LMi

K
LJik
T
JKk
(3.105)
F
Mi
=T
LMi

F
Li
(3.106)
These micro-macro equilibrium equations assist in obtaining the micro-macro
solutions just in one step and by transposing it, whereas in case of classical
homogenization method these micro-macro are uncoupled and solved separately.
98
Chapter 4
Micro Meshing Based on Digital
Imaging Techniques
Discretization of the domain, also known as meshing is the central and tedious part
in the nite element analysis. Taking into account of the multi-scale approaches in
nite element methods, it is even harder to obtain a very well dened mesh for the
purpose of computation. To overcome some of the important issues in modeling and
discretization, digital image based techniques proves to be promising and ecient,
which will be explained in upcoming sections of this chapter.
4.1 A Digital Image
\A digital image is a two-dimensional numerical representation of an image. De-
pending on whether the image resolution is xed, it may be of vector or raster type.
Without qualications, the term "digital image" usually refers to raster images
also called bitmap images. Raster images have a nite set of digital values, called
picture elements or pixels. The digital image contains a xed number of rows and
columns of pixels. Pixels are the smallest individual element in an image, holding
quantized values that represent the brightness of a given color at any specic point"
[35].
Therefore, the raster images are turned to be more powerful and ecient in ex-
hibiting its tiniest hidden details, such that there exists a maximum productivity
99
A 2-D Pixel
(a) A Typical Two-Dimensional Image Repre-
sentation by Pixels (Image Courtesy: Bernd
J ahne [5])
l
m
n
z
x
y
A 3-D Voxel
(b) A Typical Three-Dimensional Image Rep-
resentation by Voxels (Image Courtesy: Bernd
J ahne [5])
(c) A Detailed View of Pixels in a Com-
puter Screen (Image Courtesy: FunctionX
[8])
(d) A Computational Representation of a RGB
Channel (Image Courtesy: Mathworks [23])
Figure 4.1: Representations of Images
in exploiting the smallest particulars of multiple scales. Some of the typical repre-
sentations of images are shown in Figure 4.1.
To elaborate each gure in detail, Figure 4.1a represents a spatial distribution
of intensities in a two-dimensional image, and each square denotes a pixel. Further,
100
total number of rows and columns dependents on the resolution of the image.
Similarly, Figure 4.1b displays a spatial distribution of intensities in a three-
dimensional image, cubic shaped object are the voxels. Volumetric resolution is
dependent on the number of voxels in an image.
Furthermore, Figure 4.1c shows the pixels representation on a computer monitor
screen, it should be noted that in a enlarged part two pixels are red and blue
respectively. This means, if it is completely red then the RGB channel is [R-255,
G-0, B-0] and for blue color it is [R-0, G-0, B-255]. Generally, in a typical color
representation RGB channels are used, each channel is of 8 bits, thus it varies from
0 to 255. Collectively, RGB channels makes 24 bits and it is referred as true colors.
Lastly, in Figure 4.1d, a convenient data structure (a 3-D array) for storing the
RGB channels is shown. Extreme top left pixel value is stored in the extreme top
left corner of the three-dimensional array (in each plane) and subsequent pixels
values are stored in its corresponding array by traversing from extreme top left to
extreme bottom right positions. The size of the array is given by the size of an
image i.e. total number of columns-by-rows.
4.2 Importance of the Discretization Procedures
through Digital Image Processing
In nature not all the structures, textures and representations are well dened. Some
of them are highly crooked, and others are moderately distorted. This makes the
modeling and analysis processes cumbersome. Figure 4.2 shows a cross section of
the trabecular bone micro-structure which is highly complex in nature for modeling
due to its intricacies in geometric denitions.
101
vertebralbody.Figure4illustratessomeofthesekey
architecturaldifferencesforatrabecularportionofa
human vertebral body and iliac crest. The former
has a combination of rods and plates; the latter is
composed of thicker plate sections, devoid of a
truss-like structure. Aside from morphological dif-
ferences, there are clear structural anisotropy dis-
tinctions. The vertebral sample is more aligned in
the superior-inferior direction in accordance with
thedirectionofloadduetoweightinthespine.The
iliaccrestislessorganizedinitsfabricdirections.In
addition, the two samples (although unmeasured)
probably have similar yield and ultimate strains on
the apparent level, but the stresses they experience
are drastically dissimilar. Disuse osteoporosis,
whichissensitivetomechanics,affectsthevertebral
body, causing perforations and reductions in trabe-
culae,althoughtheiliaccrestisnotvulnerable.The
implication for tissue engineering is that constructs
must be designed with a priori information on
the site, its boundary geometry, and the loading it
experiences.
Mechanotransductive Considerations
This previous arguments give rise to further ques-
tions regarding the mechanotransductive nature of
bone. If, in fact, a certain variety of bone cell is
responsive to a stress-derived parameter as opposed
to a strain-derived parameter; shouldn’t the tissue
engineeringdesignstrategybetailored accordingly?
Similarly, won’t the significance of isostress and
isostrain vary with this dependence? On the other
hand, if the scaffold material or bone tissue can be
Figure 4 Site-specific architectures of the trabecular
portion of a human lumbar vertebral body (left) and iliac
crest (right).
Figure 3 Isostrain and isostress of a com-
posite material. In isostrain each component
has a uniform deformation; in isostress each
material has a uniform stress. The stress
and strain, for isostrain and isostress, re-
spectively, are in general additive but de-
pend on the moduli and volume fraction of
each component.
222 SEMINARS IN PLASTIC SURGERY/VOLUME 19, NUMBER 3 2005
Figure 4.2: A Cross Section of the Trabecular Bone Micro-structure (Image Courtesy:
Liebschner et. al. [21])
It can be noted that it is most likely impossible to draw each curve, line , voids and
other geometric details of this highly irregular structure. And further, transforming
it into a three-dimensional model for a strength analysis.
Some of the reasons for this diculties are:
Demanding time and work in modeling multiple scales domain
Micro-structures are complex in crystallography and material properties [19]
Multiple scales meshing needs excessive details of the various domains
Innitesimal scales requires a very complex algorithm to perform the dis-
cretization of its domain
Therefore, it becomes necessary to develop some techniques for modeling and
meshing of the structures with highly irregular shapes. To bolster this statement,
image processing algorithms like morphing techniques, image segmentation, edge
detection and meshing procedures are the reliable tools for obtaining accurate and
rened mesh.
102
4.2.1 A Typical Micro-structure Meshing by Imaging Tech-
niques
Micro-structures are dicult to model and mesh. Usually, a high resolution micro-
graphic picture is obtained from an electronic microscope or by tomography, where
in which each material could be distinguished explicitly and can be meshed almost
accurately. \A picture is worth a thousand words", as a matter of fact the detailed
scenario is displayed in the following Figure 4.3.
MAY/JUNE 2001 19
imation of the image, and the generated mesh is
an approximation to the material image. Forcing
the boundary to be well defined and then approx-
imating the boundary by the boundaries of the fi-
nite elements just introduces another level of ap-
proximation into the existing hierarchy.
A second approach would be to take every
pixel in the material image and create one
quadrilateral or two triangular elements from it.
ppm2oof can do this (with triangular elements
Figure 5. The image from Figure 4c, after despeck-
ling the selection and applying elkcepsed.
Despeckle has filled in small holes, and elkcepsed
has eliminated small islands and peninsulas.
Figure 6. A material image. The two colors corre-
spond to the two material properties assigned to
the micrograph in Figure 1. The few remaining
small light regions inside dark areas are artifacts
of the selection tools used, and user intervention
can easily clean them up.
(a)
(b)
(c)
Figure 4. A portion of Figure 1 after applying (a)
a nonlinear smoothing operation and repeating
the (b) burn and (c) demography selections from
Figure 3.
(a) A Typical Micro-graph of a Micro-
structure
(b) Image Segmentation Mode of a Micro-
graph
Figure 4.3: A Transition from Micro-graph to a Rened Mesh (Image Courtesy:
Langer et. al. [19])
Figure 4.3a shows the a micro-graph of a micro-structure. In Figure 4.3a the
segmented portion of the light material is highlighted by jet black color. The
image segmentation is a non-trivial process which involves complex algorithm for
dierentiating regions based on the geometry or by color.
103
MAY/JUNE 2001 21
Figure 8. Steps in the meshing process. (a) The initial uniform mesh overlaid on the material image. (b) The mesh after
10 annealing steps. (c) The triangles with the largest E are divided in two. (d) The mesh after 10 more annealing steps
and one interface refinement. (e) After further annealing and refinement, the mesh closely follows the material bound-
aries. (f) A final refinement resolves some of the smaller features. In (e) and (f) the triangles are colored according to the
material type that they inherit from the pixels.
(a)
(c)
(e)
(b)
(d)
(f)
(c) First Stage of meshing
MAY/JUNE 2001 21
Figure 8. Steps in the meshing process. (a) The initial uniform mesh overlaid on the material image. (b) The mesh after
10 annealing steps. (c) The triangles with the largest E are divided in two. (d) The mesh after 10 more annealing steps
and one interface refinement. (e) After further annealing and refinement, the mesh closely follows the material bound-
aries. (f) A final refinement resolves some of the smaller features. In (e) and (f) the triangles are colored according to the
material type that they inherit from the pixels.
(a)
(c)
(e)
(b)
(d)
(f)
(d) Final Rened Mesh
Figure 4.3: A Transition from Micro-graph to a Rened Mesh (Image Courtesy:
Langer et. al. [19])
Figure 4.3c represents the an initial uniform mesh grids on a partially smoothed
image of a micro-structure. Finally, in Figure 4.3d rened mesh grid is displayed,
in addition it should be noted that the two dierent materials are separated very
smoothly. Further, the little red triangles are the mesh grids considering material
type at each pixels.
4.3 Micro Meshing by Digital Image Processing
In this work, discretization procedures are implemented by the inspiration of image
processing and its nature of accuracy, eciency and ease of micro sub-division.
The micro elements were crucial in order to obtain the homogenized material
constants around each Gauss integration points for each macro element. Since the
geometry and ber are comprised of regular quadrilateral and straight lines, there
never existed an extra work to implement image segmentation related algorithms.
104
Global Domain
Global/Macro Element
Micro Element Micro Element Sub-Micro Element
Figure 4.4: A Pictorial Overview of an Image Based Micro Meshing
Further in this section, a brief overview of the image based meshing is presented
with an aid of Figure 4.4 which displays the complete process of discretization of
the micro elements to sub-micro elements.
The Figure 4.4 represents the exact model of functionally graded materials
(FGMs) used in this work, o-course it is one of the case solved in the Chapter 5.
Thin `red' strips are the representation of bers and the `black' area corresponds to
the matrix distribution within the macro domain. Its all start with an extraction of
a micro element from a macro element around nn Gauss integration points. In
this thesis, extraction of a micro element is done by the cropping technique, which
is highly advantageous in Mathworks
R
, Matlab Program. Further, the cropped
part is scaled accordingly to the main domain. After cropping and scaling, there
105
exists a procedure to demarcate the meshing boundaries, where the boundary lines
are 1 pixel thick, which is negligible compared to the micro or macro domain.
4.3.1 Cropping and Scaling Processes in Micro Meshing
Cropping is a process in which the the region of interest is separated from the whole
image or to change the aspect ratio. As the computation process proceeds, the
four corner coordinates of the micro elements are obtained by the iso-parametric
method. Later, these coordinates are used by the image processing program
(function). Figure 4.5 shows a cropped view of a micro element by the application
of Mathworks
R
, Matlab's function `imcrop'.
Global/Macro Element Micro Element
Figure 4.5: A Cropped View of a Micro Element
Scaling is one of the most frequent operation in micro meshing. Since there is a
dierence between the image coordinate and world/Cartesian coordinate system.
In image coordinate system the length and depth of the image are expressed by the
total number of pixels along the dimension. But, it should be noted that in image
coordinate system all the coordinates are positive integers represented by pixels
count and
oat values cannot exist. Figure 4.6 shows the image and Cartesian
coordinate system, in Figure 4.6a, it should be noted that P
(0;0)
P
(m1;n1)
represents the pixels, althoughP
(0;0)
covers bers and matrix (just for exaggeration),
its worthwhile to note that pixel will be much tinier than P
(0;0)
.
106
Lastly, Figure 4.6b exhibits the functionally graded materials (FGMs) macro domain
in general Cartesian coordinate system with origin (0,0), length and depth.
Y
X
P
(0,0)
P
(m-1,n-1)
n - Columns / n -pixels
m - rows / m - pixels
(0,0)
origin
(a) Image Coordinate System
Y
X
Length (x - units)
Depth (y-units)
(0,0)
origin
(b) Cartesian Coordinate System
Figure 4.6: A Typical Representation of Coordinate Systems
From the Figure 4.6, it is evident that there should be a conversion between the
image and Cartesian system on each occasion when the values are traversed from
back and forth.
107
Suppose for an illustration, let
Length of a cantilever FGM beam =x units
Depth of a cantilever FGM beam =y units
Similarly,
Number of pixels along horizontal direction of an image =n
Number of pixels along vertical direction of an image =m
To nd the horizontal pixel number at distance
x
2
, it follows by an elementary
equation
H
pixel
=
x
2
n
x
(4.1)
Similarly for vertical pixel number at
y
2
,
V
pixel
=
y
2
m
y
(4.2)
It is vice-versa in order to nd the distance or Cartesian parameters from a known
pixel value along horizontal/vertical sides. Lastly, one should note that there is a
considerable amount of shift in Y values of image coordinate system as compared
to the Cartesian system.
To consider the shift,
(m 1)i =
y
unknown
m
y
(4.3)
i = 0 (m 1) i.e. number of rows
108
4.3.2 Discretization of Micro Element
Meshing of the micro element is the penultimate process in imaging aspect which
precedes the coordinates storing process and succeeds after the cropping of micro
element.
Micro Element
W/ Mesh
Sub-Micro Element Micro Element
Figure 4.7: A Transition of Cropping to Meshing of a Micro Element
In the Figure 4.7, the cropping process which follows the meshing of micro element
is displayed. Number of white demarcation lines depends on the number of micro
elements required. Once the mesh size is decided upon the accuracy, white lines
of a pixel thick are drawn. Next, the pixel pointer traverses from bottom extreme
left to extreme top right (consistent with Cartesian coordinate system). As the
pixel pointer encounters a white pixel, further the pointer validates the pixel color
before and after the current white pixel, when there is a dierence of colors in the
validation process, pointer records the current pixel coordinate then immediately
maps to the world coordinate and stores the same. Suppose, just white pixels are
encountered or there is an intersection of two white lines, again the neighborhood of
current pointer is validated. The repetitive process continues till the top right pixel,
all the Cartesian coordinate values, material properties of each sub-micro element,
number of rows, columns and other details are stored for further processing.
Since the programming style varies from author to author, just a typical
owchart
is provided in the Appendix B.2 for the readers reference purpose.
109
Chapter 5
Numerical Simulations and
Results
In this chapter the performance of aperiodic homogenization assembled with a
nite element solver is quantied against the ne grid nite element. To accomplish
the ne grid solutions, a commercial nite element package, Dassault Syst emes,
Abaqus, Inc., Abaqus
R
is utilized.
To validate the method, rigorous tests are performed for ve independent cases.
They are,
Case 1: Vertical bers spacing equally
Case 2: Functionally graded vertical bers in X direction
Case 3: Horizontal bers spacing equally
Case 4: Functionally graded horizontal bers in Y direction
Case 5: Periodic bers in bi-direction
110
5.1 Geometry of the Two-Dimensional Cases
5.1.1 Macro Domain
In order to verify the results of aperiodic method, graded and periodic two-
dimensional problems were considered. Only 16 global elements are used in all the
cases, except one or two and veried against the ne grid nite element solutions.
h
x1
h
x2
h
x3
h
xn
Width
Height
Y
X
h
yn
h
y1
h
y2
h
y3
Figure 5.1: A Typical Two-Dimensional Geometry of the Aperiodic Method
Figure 5.1 shows a typical geometry of a beam used in the aperiodic nite element
analysis. h
x1
, h
x2
, h
x3
::: h
xn
are the spaces between the vertical bers. Similarly,
h
y1
, h
y2
, h
y3
::: h
yn
are the spacing between horizontal bers.
All global elements are 9 noded Lagrangian type as shown in Figure 3.5c, and
the micro elements are of 4 noded bi-linear type as shown in Figure 3.5a. Each of the
element types have two degrees of freedom, plane stress assumption is considered
for the solution of displacements.
111
5.1.2 Micro Domain
In macro scale, the required solutions are obtained by the conventional nite element
analysis, elastic constants at the integration (Gauss) point are computed by the
aperiodic homogenization method. Once the macro solution is obtained these are
Width
Height
Y
X
Sub-Micro
Elements
9 Noded Lagrangian 4 Noded Sub-Patch
Element for Micro Nodes
(a) Discretized Domain in Macro Scale
9 Noded Lagrangian
4 Noded Sub-Patch
Element
for Micro Nodes
U
9,
V
9
U
4
V
4
U
1
V
1
U
1
V
1
(b) Sub-Patch in Micro Scale
Figure 5.2: A Typical Two-Dimensional Geometry with Micro/Sub-Micro Elements
mapped back to the micro scale. In order to obtain the sub-micro elements and
their nodal values, solution (displacements) has be transformed from macro to
micro scale. To perform this, it is convenient to implement iso-parametric method
112
for 9-noded macro element given in Eq. 3.28. Further, the mapping equation is
given by,

U( x; y)!
i
( x; y)

U (5.1)
where,

U( x; y)! Displacements to be determined at micro element

i
( x; y)! Macro element shape functions evaluated at nodes of micro element

U! Known displacements of macro element
After gathering the micro elemental displacements, again the homogenization
procedure is invoked to obtain the sub-micro elemental stresses and their nodal
displacements. The transformation equation is the same as Eq. 3.31 and given by,

U

J
i
=T
JKi

U
K
i
where,

U

J
i
! Displacements at sub-micro nodes
T
JKi
! Transformation matrix obtained from homogenization

U
K
i
! Known displacements of micro element
By using

U

J
i
, displacements at required nodes and macro (or micro) elemental
stresses can be computed to understand the variation in micro scale. Similarly, if
micro displacements are known, macro solutions can be obtained by the inverse of
transformation matrix T
JKi
.
113
5.2 Case 1: Vertical Fibers Spacing Equally
A preliminary test problem is designed to validate the accuracy of the aperiodic
method of homogenization, periodic spacing bers in X direction were considered
in this case. Case summary follows as below:
Number of bers = 8
Thickness of bers = 2 units
Width of beam = 50 units
Height of beam = 43 units
Thickness of beam = 4 units
Pressure on free edge = 50 units
Young's' Modulus of bers, E
f
= 1000 units
Young's' Modulus of matrix, E
m
= 100 units
Poisson's ratio bers,
f
= 0
Poisson's ratio matrix,
m
= 0
A total of 16 global/macro elements were used to mesh the entire geometry in the
aperiodic nite element method. In addition to this, the micro meshing of each
global nite element is considered along the free end in order to extract the micro
displacements. There were about 40 micro-nodes (data collections points) at the
free end. The ber alignment, meshing details, and the boundary condition are as
shown in Figure 5.3.
114
H=43 Units
P=50 Units
W=50 Units h
xe
=50/4 Units
h
ye
=43/4 Units
(a) Vertical Fibers with Constant Spacing
Under Axial Load
H=43 Units
W=50 Units
P=50 Units
(b) Vertical Fibers with Constant Spacing
Under Shear Load
H=43 Units
P=50 Units
W=50 Units h
xe
=50/4 Units
h
ye
=43/4 Units
(c) 4x4 Meshing Details for Vertical
Fibers
Figure 5.3: Geometric Details of the Vertical Fibers with Equal Spacing
115
The horizontal spacing of the vertically periodic bers are given in the Table 5.1.
Spacing Units
h
x1
3
h
x2
4
h
x3
4
h
x4
4
h
x5
4
h
x6
4
h
x7
4
h
x8
4
h
x9
3
Table 5.1: Horizontal Spacing of the Vertical Periodic Fibers
In the conventional FEA program, the model required a total of 17200
global/macro elements in order to obtain accurate results. In the same model,
the aperiodic method used 16 macro elements to compute the local and global
deformation values. The model summary for the horizontally periodic bers case is
presented in the Table 5.2.
Specications Abaqus Aperiodic
Total Global Elements 17200 16
Integration Points per Element 4 (2x2) 4 (2x2)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 17473/2/34946 81/2/162
Table 5.2: Summary of FEA Details
116
5.2.1 Case 1a: Global and Local Deformation Analysis for
Axial Loading
The deformation due to the axial loading at the free end of the cantilever by
aperiodic method is pictured in Figure 5.4.
Figure 5.4: An Axial Deformation Plot by Aperiodic Method (Scale: 0.1)
The micro and macro displacements along the bottom edge of the beam are
compared between aperiodic and ne grid nite element method and presented in
Figure 5.5.
117
0 2 4 6 8 10 12 14 16 18
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
Figure 5.5: Plot of a Global-Local X Displacement at the Bottom Edge
It can be observed that the micro displacements obtained from the aperiodic method
exactly followed the Abaqus, FEA results and also the average slope of the global
displacements matched the slope of Abaqus, FEA.
118
5.2.2 Case 1b: Global and Local Deformation Analysis for
Shear Loading
In this case, the geometric and FEA details remains the same, except the shear
load along the free end. The deformation due to the shear loading at the free end
of the cantilever obtained from the aperiodic method is pictured in Figure 5.6.
Figure 5.6: A Shear Deformation Plot by Aperiodic Method (Scale: 0.1)
The micro and macro displacements are compared along the extreme edges and
middle intercepts of the beam for aperiodic and ne grid nite element method
and presented in the following gures.
119
−70 −60 −50 −40 −30 −20 −10 0
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at Y=0
0 5 10 15
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at Y=25
Figure 5.7: Plot of a Global-Local X Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
There is a signicant convergence in the aperiodic solution as compared to the
Abaqus, FEA when the nodes of interest are changed. In Figure 5.7a, all the nodes
along Y=0 are considered which is an exterior part of a domain. Further, in Figure
5.7b there a rapid convergence when the Y-intercept is taken near the mid-part of
the domain.
120
In the similar way, the displacements for the Y direction are considered by using
the Y-intercepts Y=0 and Y=25 in gures 5.8a and 5.8b.
−150 −100 −50 0
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local Y Displacement at Y=0
−140 −120 −100 −80 −60 −40 −20 0
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local Y Displacement at Y=25
Figure 5.8: Plot of a Global-Local Y Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
It is evident that there is a subtle change in the slope for the two plots at Y=0 and
Y=25 when compared to the macro and Abaqus, FEA solutions.
121
Next, considering the displacements in X direction for node at X-intercepts at
X=25 and X=50 (i.e. the extreme right edge). We obtain, the gures 5.9a and
5.9b. From the Figure 5.9a, it is observable that there is almost one-to-one match
−50 −40 −30 −20 −10 0 10 20 30 40 50
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at X=25
−80 −60 −40 −20 0 20 40 60 80
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at X=50
Figure 5.9: Plot of a Global-Local X Displacement at Two Locations X=25 and
X=50 (X-Intercepts)
in the displacements for nodes at X=25 against Abaqus, FEA and a marginally
dierence can be seen from 5.9b as we traverse towards the end of the beam.
122
5.3 Case 2: Functionally Graded Vertical Fibers
in X Direction
In this case the accuracy of the aperiodic method is compared for graded vertical
bers in X direction with a shear force at the free end. Case summary follows as
below:
Number of bers = 32
Thickness of bers = 0.5 units
Width of beam = 50 units
Height of beam = 43 units
Thickness of beam = 4 units
Pressure on free edge = 50 units
Young's' Modulus of bers, E
f
= 1000 units
Young's' Modulus of matrix E
m
= 100 units
Poisson's ratio bers,
f
= 0:3
Poisson's ratio matrix,
m
= 0:3
A total of 16 global/macro elements were used to mesh the entire geometry in
the aperiodic nite element method. In addition to this, there were about 40
micro-nodes (data collections points) at the bottom edge of the beam. The ber
alignment, meshing details, and the boundary conditions is as shown in Figure 5.10.
123
W=50 Units
H=43 Units
P=50 Units
(a) Functionally Graded Vertical Fibers in X
Direction Under Axial Load
P=50 Units
W=50 Units
H=43 Units
h
ye
=43/4 Units
h
xe
=50/4 Units
(b) Functionally Graded Vertical Fibers in
X Direction Under Shear Load
P=50 Units
W=50 Units
H=43 Units
h
ye
=43/4 Units
h
xe
=50/4 Units
(c) 4x4 Meshing Details for FG Vertical
Fibers in X Direction
Figure 5.10: Geometric Details of the Functionally Graded Vertical Fibers
124
The summation of 5% increase to the previous spacing of the vertical bers are
given in the Table 5.3
Spacing Units Spacing Units
h
x1
0.300 h
x17
1.039
h
x2
0.500 h
x18
1.091
h
x3
0.525 h
x19
1.146
h
x4
0.551 h
x20
1.203
h
x5
0.578 h
x21
1.263
h
x6
0.607 h
x22
1.326
h
x7
0.638 h
x23
1.393
h
x8
0.670 h
x24
1.462
h
x9
0.703 h
x25
1.535
h
x10
0.738 h
x26
1.612
h
x11
0.775 h
x27
1.693
h
x12
0.814 h
x28
1.777
h
x13
0.855 h
x29
1.866
h
x14
0.897 h
x30
1.960
h
x15
0.942 h
x31
2.058
h
x16
0.990 h
x32
0.980
Table 5.3: Horizontal Spacing of the Functionally Graded Vertical Fibers
In the conventional FEA program, the model required a total of 12528 global/macro
elements in order to obtain accuracy. In the same model, the aperiodic method
used 16 macro elements to compute the local and global deformation values. The
model summary for the horizontally graded bers case is presented in the Table
5.4.
Specications Abaqus Aperiodic
Total Global Elements 17200 16
Integration Points per Element 4 (2x2) 16 (4x4)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 17473/2/34946 81/2/162
Table 5.4: Summary of FEA Details
125
5.3.1 Case 2a: Global and Local Deformation Analysis for
Axial Loading
The deformation due to the axial loading at the free end of the cantilever obtained
from the aperiodic method is pictured in Figure 5.11.
Figure 5.11: An Axial Deformation Plot by Aperiodic Method (Scale: 1)
The micro (local) and macro (global) X displacements along the Y-intercepts, Y=0
and Y=25 are compared between aperiodic and ne grid nite element method and
are presented in gures 5.12.
126
0 2 4 6 8 10 12 14 16 18
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at Y=0
0 2 4 6 8 10 12 14 16 18
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at Y=25
Figure 5.12: Plot of a Global-Local X Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
It is evident that there is a negligible amount of change in the slope of the gures
5.12a and 5.12b.
127
Further, Y displacements for the nodes along the edges at Y=0 and Y=25 are given
in the Figure 5.13.
0 0.2 0.4 0.6 0.8 1 1.2 1.4
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local Y Displacement at Y=0
−0.35 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local Y Displacement at Y=25
Figure 5.13: Plot of a Global-Local Y Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
From gures 5.13a and 5.13b, it can be seen that there is a considerable amount of
improvement in the displacement values as the nodes of interest lies close to the
neutral axis of the beam.
128
5.3.2 Case 2b: Global and Local Deformation Analysis for
Shear Loading
The deformation due to the shear loading at the free end of the cantilever obtained
from the aperiodic method is given in Figure 5.14.
Figure 5.14: A Shear Deformation Plot by Aperiodic Method (Scale: 0.1)
The micro (local) and macro (global) X displacements along the Y-intercepts,
Y=0 and Y=25 are compared between aperiodic and ne grid nite element method
and are presented in Figure 5.15.
129
−70 −60 −50 −40 −30 −20 −10 0
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at Y=0
0 2 4 6 8 10 12
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at Y=25
Figure 5.15: Plot of a Global-Local X Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
It is evident from the gures 5.15a and 5.15b, that there is a sucient amount of
change in the slope and a rapid convergence of the solution at the region close to
the middle part of the beam.
130
Further, Y displacements for the nodes along the edges at Y=0 and Y=25 are given
in the Figure 5.16.
−150 −100 −50 0
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local Y Displacement at Y=0
−140 −120 −100 −80 −60 −40 −20 0
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local Y Displacement at Y=25
Figure 5.16: Plot of a Global-Local Y Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
It can be noticed that there is a little change in the slope or accuracy for Y
displacements when the interior nodes are considered.
131
5.4 Case 3: Horizontal Fibers Spacing Equally
The performance of the aperiodic method is compared for a periodic horizontal
bers in Y direction with a axial and shear force at the free end. Case summary
follows as below:
Number of bers = 8
Thickness of bers = 2 units
Width of beam = 50 units
Height of beam = 43 units
Thickness of beam = 4 units
Pressure on free edge = 50 units
Young's' Modulus of bers, E
f
= 1000 units
Young's' Modulus of matrix, E
m
= 100 units
Poisson's ratio bers,
f
= 0
Poisson's ratio matrix,
m
= 0
A total of 16 global/macro elements were used to mesh the entire geometry in the
aperiodic nite element method, were used to model the shear and again 504 global
elements are used to model the axial case. In addition to this, there were about
40 micro-nodes (data collections points) at the right edge of the beam. The ber
alignment, meshing details, and the boundary conditions is as shown in Figure 5.17.
132
H=43 Units
P=50 Units
W=50 Units h
xe
=50/4 Units
h
ye
=43/4 Units
(a) Horizontal Vertical with Constant Spac-
ing Under Axial Load
W=50 Units
H=43 Units
P=50 Units
(b) Horizontal Vertical with Constant
Spacing Under Shear Load
H=43 Units
P=50 Units
W=50 Units h
xe
=50/4 Units
h
ye
=43/4 Units
(c) 4x4 Meshing Details for Horizontal
Vertical with Constant Spacing
Figure 5.17: Geometric Details of the Horizontal Fibers with Constant Spacing
The vertical spacing of the horizontal periodic bers are given in the Table 5.5
133
Spacing Units
h
x1
3
h
x2
3
h
x3
3
h
x4
3
h
x5
3
h
x6
3
h
x7
3
h
x8
3
h
x9
3
Table 5.5: Vertical Spacing of the Horizontally Periodic Fibers
5.4.1 Case 3a: Global and Local Deformation Analysis for
Axial Loading
In the conventional FEA program, the model required a total of 17800 global/macro
elements in order to obtain accurate results. In the same model, the aperiodic
method used 504 macro elements to compute the local and global deformation
values. The model summary for the vertically periodic bers case is presented in
the Table 5.6.
Specications Abaqus Aperiodic
Total Global Elements 17800 504
Integration Points per Element 4 (2x2) 4 (2x2)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 18090/2/36180 2117/2/4234
Table 5.6: Summary of FEA Details
The deformation due to the axial loading at the free end of the cantilever obtained
from the aperiodic method is pictured in Figure 5.18.
134
Figure 5.18: An Axial Deformation Plot by Aperiodic Method (Scale: 1)
The X and Y micro (local) and macro (global) displacements results at the right
edge of the beam are compared between aperiodic and ne grid nite element
method and are presented in gures 5.19 and 5.20.
135
5 5.5 6 6.5 7 7.5 8 8.5
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
Figure 5.19: Plot of a Global-Local X Displacement at the Right Edge (Axial)
136
−1.5 −1 −0.5 0 0.5 1 1.5
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
Y Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
Figure 5.20: Plot of a Global-Local Y Displacement at the Right Edge (Axial)
It can be seen from the gures 5.19 and 5.20 that the micro displacements obtained
from the aperiodic method closely followed the Abaqus, FEA results in the middle
part along the edge of the beam, also there is a considerable phase dierence at
the top and bottom ends.
137
5.4.2 Case 3b: Global and Local Deformation Analysis for
Shear Loading
In the conventional FEA program, the model required a total of 17800 global/macro
elements in order to obtain accurate results. In the same model, the aperiodic
method used 16 macro elements to compute the local and global deformation values.
The model summary for the vertically periodic bers case is presented in the Table
5.7. The deformation due to the shear loading at the free end of the cantilever
Specications Abaqus Aperiodic
Total Global Elements 17800 16
Integration Points per Element 4 (2x2) 16 (4x4)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 18090/2/36180 81/2/162
Table 5.7: Summary of FEA Details
obtained from the aperiodic method is pictured in Figure 5.21.
Figure 5.21: A Shear Deformation Plot by Aperiodic Method (Scale:0.1)
138
The X displacements along the nodes on X=25 and X=50 intercepts of the beam are
compared between aperiodic and ne grid nite element method and are presented
in Figure 5.22.
−25 −20 −15 −10 −5 0 5 10 15 20 25
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at X=25
−40 −30 −20 −10 0 10 20 30 40
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at X=50
Figure 5.22: Plot of a Global-Local X Displacement at Two Locations X=25 and
X=50 (Y-Intercepts)
It can be observed from the gure 5.22a and 5.22b, the micro displacements are
close to the Abaqus, FEA results, in the middle part of the beam, but again
139
there is a notable dierence at the very end of the beam results. Further, there
is considerable change in displacements as we traverse away from the exterior
edge. From the gures 5.23a and 5.23b it can be noted that there is a large phase
−40 −35 −30 −25 −20 −15 −10
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
Y Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at X=25
−110 −100 −90 −80 −70 −60 −50 −40 −30
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
Y Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at X=50
Figure 5.23: Plot of a Global-Local Y Displacement at Two Locations X=25 and
X=50 (Y-Intercepts)
dierence between Abaqus, FEA and aperiodic method. Further, there is signicant
amount of decrease in error at the ends when interior nodes are taken into account.
140
5.5 Case 4: Functionally Graded Horizontal
Fibers in Y Direction
A vertically graded horizontal bers in Y direction with a axial force at the free
end is considered for this case. Case summary follows as below:
Number of bers = 29
Thickness of bers = 0.5 units
Width of beam = 50 units
Height of beam = 43 units
Thickness of beam = 4 units
Pressure on free edge = 50 units
Young's' Modulus of bers, E
f
= 1000 units
Young's' Modulus of matrix E
m
= 100 units
Poisson's ratio bers,
f
= 0
Poisson's ratio matrix,
m
= 0
A total of 16 global/macro elements were used to mesh the entire geometry in the
aperiodic nite element method, were used. In addition to this, there were about
40 micro-nodes (data collections points) at the bottom and right edge of the beam.
The ber alignment, meshing details, and the boundary conditions is as shown in
Figure 5.24.
141
H=43 Units
P=50 Units
W=50 Units
h
ye
=43/4 Units
h
xe
=50/4 Units
(a) Functionally Graded Horizontal Fibers
in Y Direction
H=43 Units
P=50 Units
W=50 Units
h
ye
=43/4 Units
h
xe
=50/4 Units
(b) 4x4 Meshing Details for FG Hori-
zontal Fibers in X Direction
Figure 5.24: Geometric Details of the Functionally Graded Horizontal Fibers
The summation of 5% increase to the previous spacing of the horizontal bers are
given in the Table 5.8
Spacing Units Spacing Units
h
x1
0.300 h
x17
1.039
h
x2
0.500 h
x18
1.091
h
x3
0.525 h
x19
1.146
h
x4
0.551 h
x20
1.203
h
x5
0.578 h
x21
1.263
h
x6
0.607 h
x22
1.326
h
x7
0.638 h
x23
1.393
h
x8
0.670 h
x24
1.462
h
x9
0.703 h
x25
1.535
h
x10
0.738 h
x26
1.612
h
x11
0.775 h
x27
1.693
h
x12
0.814 h
x28
1.777
h
x13
0.855 h
x29
1.365
h
x14
0.897 - -
h
x15
0.942 - -
h
x16
0.990 - -
Table 5.8: Horizontal Spacing of the Functionally Graded Horizontal Fibers
142
The deformation due to the axial loading at the free end of the cantilever obtained
from the aperiodic method is pictured in Figure 5.25.
Figure 5.25: An Axial Deformation Plot by Aperiodic Method (Scale: 1)
5.5.1 Case 4: Global and Local Deformation Analysis
Fine grid FEA ides 8600 global/macro elements in order to obtain accurate results.
In the same model, the aperiodic method used 16 macro elements to compute
the local and global deformation values. The model summary for the functionally
graded horizontal bers case is presented in the Table 5.9.
Specications Abaqus Aperiodic
Total Global Elements 8600 16
Integration Points per Element 4 (2x2) 16 (4x4)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 8787/2/17574 81/2/162
Table 5.9: Summary of FEA Details
143
The X nodal displacements (global and local) for the right edge i.e. at X=50 and
X=25 are compared between aperiodic and ne grid nite element method and are
presented in Figure 5.26.
−0.68 −0.66 −0.64 −0.62
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
Y Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at X=25
4 5 6 7 8 9 10
0
5
10
15
20
25
30
35
40
45
Global−Micro Displacements
X Displacement
Node Position (Y−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at X=50
Figure 5.26: Plot of a Global-Local X Displacement at Two Locations X=25 and
X=50 (X-Intercepts)
It is evident that there is a minimization in the error dierence when compared
to exterior and interior nodes. Further, the high variation of materials is partially
captured by aperiodic method.
144
In the following Figure 5.27 global/local X displacements are determined at Y=0
and Y=25.
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local X Displacement at Y=0
0 1 2 3 4 5 6 7 8
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at Y=25
Figure 5.27: Plot of a Global-Local X Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
From the gures 5.27a and 5.27b, it can be visualized that there is a no signicant
change in the displacements at Y=0. Whereas at Y=25 there is a slight dierence
in the slope of the micro solutions.
145
Lastly, a comparison of the Y displacements are shown for the nodes along Y=0
and Y=25 in Figure 5.28.
−3 −2.5 −2 −1.5 −1 −0.5 0 0.5
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(a) Plot of a Global-Local Y Displacement at Y=0
−3 −2.5 −2 −1.5 −1 −0.5 0 0.5
0
10
20
30
40
50
60
Global−Micro Displacements
Y Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
(b) Plot of a Global-Local X Displacement at Y=25
Figure 5.28: Plot of a Global-Local Y Displacement at Two Locations Y=0 and
Y=25 (Y-Intercepts)
In the Figure 5.28a there exists subtle dierence at the right corner end (bottom),
but in Figure 5.28b it is evident that there is a minimum error.
146
5.6 Case 5: Periodic Fibers in Bi-direction
A vertically and horizontally periodic bers in both X and Y direction is considered
for axial and shear type of loading. Case summary follows as below:
Number of bers in X direction = 15
Number of bers in Y direction = 15
Thickness of bers = 1 units
Width of beam = 50 units
Height of beam = 43 units
Thickness of beam = 4 units
Pressure on free edge = 50 units
Young's' Modulus of bers, E
f
= 1000 units
Young's' Modulus of matrix E
m
= 100 units
Poisson's ratio bers,
f
= 0
Poisson's ratio matrix,
m
= 0
A total of sixteen global/macro elements were used to mesh the entire geometry
in the aperiodic nite element method. In addition to this, there were about 40
micro-nodes (data collections points) at the bottom edge of the beam. The ber
alignment, meshing details, and the boundary conditions is as shown in Figure 5.29.
147
W=50 Units
H=43 Units
P=50 Units
h
xe
=50/4 Units
h
ye
=43/4 Units
(a) Peridic Fibers in Two Directions Under
Axial Load
W=50 Units
h
xe
=50/4 Units
h
ye
=43/4 Units
H=43 Units
P=50 Units
(b) Peridic Fibers in Two Directions
Under Shear Load
W=50 Units
h
xe
=50/4 Units
h
ye
=43/4 Units
H=43 Units
P=50 Units
(c) 4x4 Meshing Details for Peridic
Fibers in Two Directions
Figure 5.29: Geometric Details of the Periodic Fibers in Two Directions
The spacing of the horizontal and vertical bers are given in the Table 5.10
148
Spacing Units Spacing Units
h
x1
1.000 h
y1
1.000
h
x2
2.357 h
y2
1.857
h
x3
2.357 h
y3
1.857
h
x4
2.357 h
y4
1.857
h
x5
2.357 h
y5
1.857
h
x6
2.357 h
y6
1.857
h
x7
2.357 h
y7
1.857
h
x8
2.357 h
y8
1.857
h
x9
2.357 h
y9
1.857
h
x10
2.357 h
y10
1.857
h
x11
2.357 h
y11
1.857
h
x12
2.357 h
y12
1.857
h
x13
2.357 h
y13
1.857
h
x14
2.357 h
y14
1.857
h
x15
2.357 h
y15
1.857
h
x16
1.000 h
y16
1.000
Table 5.10: Horizontal and Vertical Spacing of the Periodic Fibers in Two Directions
In ne grid FE program, the model required a total of 9360 global/macro elements
in order to obtain accurate results, whereas the aperiodic method used 16 macro
elements to compute the local and global deformation values. The model summary
of periodic bers in two directions considering both the loading conditions is given
in the Table 5.11.
Specications Abaqus Aperiodic
Total Global Elements 9360 16
Integration Points per Element 4 (2x2) 16 (4x4)
Element Type 4 Noded Bi-linear 9 Noded Lagrangian
Nodes/DOF/Total DOF 9555/2/19110 81/2/162
Table 5.11: Summary of FEA Details
149
5.6.1 Case 5a: Global and Local Deformation Analysis for
Axial Loading
The deformation due to the axial loading at the free end of the cantilever obtained
by the aperiodic method is pictured in Figure 5.30.
Figure 5.30: An Axial Deformation Plot by Aperiodic Method (Scale: 1)
The X micro and macro displacements results at the bottom edge of the beam are
compared between aperiodic and ne grid nite element method and presented in
the Figure 5.31.
150
0 1 2 3 4 5 6
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
Figure 5.31: Plot of a Global-Local X Displacement at the Bottom Edge (Axial)
The micro displacements Figure 5.31 obtained by the aperiodic method almost
traced the Abaqus, FEA results. However, there is a sudden jump in the displace-
ment value at the corner of the bottom edge as shown by Abqus, FEA result which
in deed is not predicted by the aperiodic method.
151
5.6.2 Case 5b: Global and Local Deformation Analysis for
Shear Loading
The deformation due to the shear loading at the free end of the cantilever obtained by
the aperiodic method is pictured in Figure 5.32. Micro/local and the global/macro
Figure 5.32: A Shear Deformation Plot by Aperiodic Method (Scale: 0.1)
displacements in the direction of X at the bottom edge of the beam retrieve from the
aperiodic approach are compared between against ne grid nite element method
and presented in Figure 5.33.
152
−25 −20 −15 −10 −5 0
0
10
20
30
40
50
60
Global−Micro Displacements
X Displacement
Node Position (X−Coordinates)

Homo−Global
Homo−Micro
Abaqus
Figure 5.33: Plot of a Global-Local X Displacement at the Bottom Edge (Shear)
It can be observed from the Figure 5.33 that the tiny disturbances in the
displacements are captured by the aperiodic method. Despite the method under-
estimated the displacements at the end of bottom edge, further there is a notable
dierence in the phase of micro-displacements as traverse to the extreme right side
of the bottom edge.
153
Chapter 6
Summary, Conclusions and
Recommendations
6.1 Summary
Composite materials are characterized by the micro-structures to attain excellent
strength to weight ratio. Micro-structures are highly heterogeneous in nature and
there by it requires distinguished methods for analysis of its structural or thermal
property. In this thesis, a basic framework is established by using an unique
approach of conventional homogenization theory and casting it into a generalized
technique for analyzing aperidic micro-structures. This new approach is based
on capturing the nest variations in the composite material medium through the
introduction of micro-solution variables. These solution variables mimics the role
of cell solutions in classical homogenization theory.
Multi-scale methods are the ultimate tools for bridging the gap between space
and time continuum. This work is developed on the emphasis to transforming
the quantitative data of micro scale and vice-versa. Thereby including the eects
of multi-scale coupling in obtaining the homogenized elastic tensor. Aperiodic
homogenization technique integrated with a general purpose nite element program
was used to perform several common types of nite element composite structure
analysis for a linear elastic material. The results obtained from the aperiodic
154
method were compared against the solution of Abaqus ne grid models and the
agreement was found to be highly satisfying for most of the cases.
In order to quantify the true performance and accuracy of the aperiodic homog-
enization method, micro and macro/global displacement values were computed for
each independent cases. A simple vertical periodic bers included within matrix
and in
uenced by the axial load at the free end (Case 1a). It was observed that
the aperiodic micro displacements exactly followed the Abaqus result. In addition
to this,for shear loading (Case 1b), displacement solutions at various locations of
the beam (X and Y intercepts) are testied against Abaqus, ne grid results. It
is worth to mention that the relative error of the dierence between the aperiodic
method and Abaqus results started minimizing when solution along the interior
nodes and exterior nodes is compared.
To further fatigue the method, functionally graded vertical bers are considered
with a application of axial (Case 2a) and shearing load (Case 2b) on the free end
of the cantilever beam. In this case it was evident that the results obtained from
aperiodic method are 92% accurate in the axial displacements (X direction) and
20% drifting from the Abaqus results for the transverse displacement (Y direction).
Although, for some solutions there existed a considerable variation in accuracy in
aperiodic method as against Abaqus, ne grid solution. It can be observed that
there is a relatively better results when interior nodes are considered for obtaining
the displacements solution.
For the Case 3a, periodic horizontal bers varying in Y direction under axial
load, it was impossible to obtain the satisfying result without increasing the total
number of macro elements. It is obvious that the pressure is uniformly distributed
in Abaqus, whereas in case of aperiodic method it was hardly acting on ten nodes.
155
Therefore, this forced aperiodic method to demand for more number of global
elements than local. It is quite dierent for the Case 3b, where bending dominated
the problem. The micro displacements in axial direction were almost close to the
Abaqus result except at the extreme nodes in top and bottom of that edge. This
could be justied as there were only 16 global elements, the accuracy within the
global structure itself is limited and micro displacements are the special inheritances
of the macro displacements. Thus, the results at the ends are purely dependent on
the macro displacements. In transverse displacement, there is considerable amount
of dierence in the result against Abaqus solution. Since aperiodic method never
considered the eect of strains due to bending, this led the method to under predict
the bending eects. Furthermore, a 4 noded bi-linear element is weak in predicting
the eects due to bending. Also, in this case, error is minimal as we consider the
solution along the nodes in interior part.
Case 4, which is a functionally graded type consisting the horizontally aligned
bers. Micro displacements at the right edge of beam, hardly traced the result of
Abaqus, this is again due to the insucient number of nodes for the distribution
of the axial load. Nevertheless, the displacements far away from the edges are
consistent with the Abaqus, ne grid FEA. Further, the displacements at the
bottom edge were almost close to the Abaqus result.
Lastly, in Case 5a, the results from the aperiodic method were exact for the
axial case, except there exist a slight dierence in the solution for bottom right
end nodes. Despite the displacements due to the bending (Case 5b) there exists a
thin amount of error corresponding to the Abaqus values. The reasons as discussed
earlier for the bending case holds good for this case too.
156
6.2 Conclusions
Multi-scale techniques are the viable tool for gathering detailed information about
the behavior of the micro-structures of the composites and the composites itself.
Aperiodic method proves to be worthy in transforming the required quantities
from two scales. From the amount of mathematical work involved in producing
the results it is noteworthy to claim that the aperiodic method is highly computer
implementation friendly. Further, it just takes a few extra subroutines to be
integrated with a general purpose nite element program.
Secondly, the theory was veried for just a one-dimensional problem including
some special cases. In addition to this, various complex cases in two-dimensions are
also validated. It is evident from the results of the aperiodic analysis, that there
exists a fair amount of similarity between the Abaqus results. However, only 16
global elements are used to obtain the aperiodic result against 12000 elements on
an average in Abaqus FEA.
Thirdly, it is extremely dicult to obtain the converged result along the edges
because of limited number of solution points. Furthermore, the distribution of
forces are only concentrated at the exterior nodes. This exerts an uneven force to
the ber and matrix nodes, in-turn results in non-smooth solution along the edges.
Therefore, it is important to note that the displacement solutions obtained by the
aperiodic method along the exterior edges of the beam are tend be very close or less
deviating in some cases to the Abaqus, FEA. In contrast to the exterior edges, when
the interior nodes are considered for most of the cases, aperiodic homogenization
method provided some surprising convergence to the Abaqus, FEA.
Lastly, this method reduced the computational time by almost 90% compared
to the conventional nite element program
157
6.3 Recommendations
6.3.1 Type of Finite Elements
It is evident that the 4 noded bi-linear element was ecient in obtaining the axial
deformation values almost precise to Abaqus results. However, this element failed
to work in a full range for the problems predominant in bending. Therefore, its
always advantageous to prefer higher order elements like 8 noded serendipity or 9
noded Lagrangian type that are capable of predicting the bending eect.
6.3.2 Meshing
Discretization of the domain is equally important as compared to the nal result.
Meshing of the structure plays a key role in obtaining the accurate outputs. It is
practically impossible to mesh highly heterogeneous micro-structures by conven-
tional geometric modeling methods without having some restrictions or predening
the type of boundary. As a result, in this work, digital image based meshing
technique is implemented to discretize the graded and periodic micro-structures.
Taking into consideration of
oat to integer conversion and vice-versa in both
image and Cartesian coordinates, it is a considerable fact that the mesh coordinates
obtained were almost accurate upto fourth decimal digit.
Also it should be noted that the use of graded meshing and adaptive approaches
in aperiodic method could have obtained an excellent result [4, 10]. Therefore, on
any day using an adaptive or graded meshes provide a very good result as opposed
to a non-uniform mesh.
Furthermore, image based meshing proves to overcome the diculties in model-
ing heterogeneous materials compared to the conventional technique. As a result,
158
three-dimensional models obtained by the X-ray imaging for metals and magnetic
resonance imaging (MRI) for body tissues and bone fractures can be used to perform
three-dimensional (even two-dimensional) meshing. Thereby aperiodic method with
certain changes in the formulation could be used for analyzing the structural or
thermal aspects of heterogeneity through nite element analysis.
6.3.3 Optimization
On an average aperiodic method is inexpensive as compared with other general
purpose programs (not a commercial one). In a general purpose program it would
be around 8000 seconds to solve a problem with 5000 DOF, whereas aperiodic
method utilized just 75 seconds on an avearage for all the problems presented
in this work, this reduces 99% of the computational expenses without any kind
of optimization. Further, with optimization the results could be obtained in few
seconds.
6.3.4 Types of Composites
In this work, only basic functionally graded and periodic micro-structures were
considered, but the method is capable of solving other micro-structures with some
ne tuning in the change of element type, may be for triangular or tetrahedron
for three-dimensional case. Also, it should be noted that the method developed is
multi-dimensional, but it is implemented for only two-dimensional problems.
6.3.5 Modeling of Interfaces
One of the area in which the most work needed lies in the development of accurate
models for predicting the behavior of the ber-matrix interface. The micro scale
159
modeling is be a handy tool to obtain the tiny disturbances in the interface,
bonding and debonding of two or more materials. Further, opens an opportunities
in developing groundbreaking techniques for nano scale technology that could solve
the inter scale mysteries.
6.3.6 Failure Analysis
Failure of composites is one of the important criterion to be considered in many
industries, especially in Aerospace. This method could be used for obtaining the
stresses and strains in the region of interest with less computational time. The
von-mises stresses obtained from aperiodic method could be compared against the
ultimate failure capacity of the material, and an intuitive prediction could be made
about the failure.
6.3.7 Non Destructive Testing (NDT)
NDTs are the most popular analysis methods to evaluate the material, strength
and other property of a component or system without causing any damage to
the structure. Especially, X-ray, CT and techniques can pierce into the hidden
world of materials. By utilizing these nest techniques, a two-dimensional or three-
dimensional image can be reconstructed. Once the images are reconstructed, digital
image based meshing can be accomplished for nite element analysis. It should
be noted that, since aperiodic method can accommodate non-periodic boundary
conditions and material variations, this method is more
exible for analysis purposes.
160
Bibliography
[1] Aboudi Jacob, 1983, \The Eective Moduli of Short-Fiber Composites," Int. J.
Solid Structures, Pergamon Press Ltd., Great Britain, 19, No. 8, pp. 693-707
[2] Aboudi J., Pindera M.J., and Arnold S.M., 1999, \Higher-order Theory for
Functionally Graded Materials," Elsevier Science Limited, Composites: Part
B Engineering, 30, pp. 777-832
[3] Autar K. Kaw, 2006, \Mechanics of Composites Materials," CRC Press, Taylor
& Francis Group, LLC, Second Edition, ISBN: 978-0-8493-1343-1
[4] Bende, Martin Philip and Kikuchi, Noboru 1988, \Generating Optimal
Topologies in Structural Design Using A Homogenization Method," Computer
MEthods in Applied Mechanics and Engineering, Elsevier Science Publishers
B.V., North-Holland, 71, pp. 197-224
[5] Bernd J ahne, 2005, \Digital Image Processing," Springer-Verlag Berlin
Heidelberg, 6th Revised and Extended Edition), ISBN: 978-3-540-24035-
8 Springer Berlin Heidelberg New York, pp. 32
[6] Doina Cioranescu, and Patrizia Donato, 1999, \An Introduction to Homoge-
nization," Oxford Lecture Series in Mathematics and its Applications, Oxford
University Press Inc., New York, 17, ISBN: 0 19 856554 2
161
[7] Ever J. Barbero, 2011, \Introduction to Composites Materials Design," CRC
Press, Taylor & Francis Group, LLC, Second Edition, ISBN: 978-1-4200-
7915-9
[8] FunctionX Tutorials, 2012, \Colors Fundamentals," http://www.functionx.
com/vccli/gdi+/color.htm, Note: Online; accessed 30-August-2012 @ 08:45
PM, Pacic Standard Time (PST)
[9] George F. Vander Voort, (Editor), 2004, \ASM Handbook
R
- Metallography
and Microstructures," ASM International
R
, 09, ISBN: 0-87170-706-3
[10] Jo se Miranda Guedes and Noboru Kikuchi, 1990, \Preprocessing and post-
processing for Materials Based on the Homogenization Method with Adaptive
Finite Element Methods," Computer Methods in Applied Mechanics and
Engineering, North-Holland, 83, pp. 143-198
[11] Somnath Ghosh, Kyunghoon Lee, and Suresh Moorthy, 1995, \Multiple Scale
Analysis of Heterogeneous Elastic Structures Using Homogenization Theory
and Voronoi Cell Finite Element Method," Int. J. Solids Structures, Elsevier
Science Ltd., 32, No.1, pp. 27-62
[12] Hill R, 1963, \Elastic Properties of Reinforced Solids: Some Theoretical Prin-
ciples," J. Mech. Phys. Solids, Pergamon Press Ltd., 11, pp. 357-372
[13] Hill R, 1965, \A Self-Consistent Mechanics of Composite Materials," J. Mech.
Phys. Solids, Pergamon Press Ltd., Great Britain, 13, pp. 213-222
[14] Hirai T, 1996, \Functional Gradient Materials, In: Brook R J, editor," Wein-
heim, Germany: VCH Verlagsgesellschaft mbH Publishers, Processing of
ceramics, Part 2, pp. 293-341
[15] Hollister S. J., and Kikuchi, N., 1994, \Homogenization Theory and Digital
Imaging: A Basis for Studying the Mechanics and Design Principles of Bone
Tissue," John Wiley & Sons, Inc., Biotechnology and Bioengineering, 43, pp.
586-596
162
[16] Jeong-Ho Kim, and G. H. Paulino, 2002, \Isoparametric Graded Finite El-
ements for Nonhomogeneous Isotropic and Orthotropic Materials," ASME,
Journal of Applied Mechanics, 69, pp. 502-514
[17] Koizumi M., 1973, \FGM Activities in Japan," Elsevier Science Limited,
Composites: Part B Engineering, 28B, pp. 1-4
[18] Kouznetsova V., Geers M. G. D., and Brekelmans W. A. M., 2002, \Multi-
Scale Constitutive Modelling of heterogeneous Materials with a Gradient-
Enhanced Computational Homogenization Scheme," John Wiley & Sons, Ltd.,
International Journal for Numerical Methods in Engineering, 54, pp. 1235-1260
[19] Stephen A. Langer, Edwin R. Fuller, Jr., and W. Craig Carter, 2001, \OOF:
An Image-Based Finite-Element Analysis of Material Microstructures," IEEE,
Materials Science, pp. 15-23
[20] Michael J. Leamy, Peter W. Chung, and Raju Namburu, 2004, \Development
of a Nonperiodic Homogenization Method for One-Dimensional Continua,"
Army Research Laboratory, U.S. Military Academy, pp. 1-13
[21] Michael Liebschner, Ph.D., Brandon Bucklen, B.S., and Matthew Wettergreen,
B.S., 2005, \Mechanical Aspects of Tissue Engineering," Tissue Repair, Re-
generation, and Engineering in Plastic Surgery; Seminars in Plastic Surgery,
Thieme Medical Publishers, Inc., 333 Seventh Avenue, New York, NY 10001,
USA., Number 3, 19, pp. 217-228 (Specic: pp. 222)
[22] Luc Tartar, 2009, \The General Theory of Homogenization, A Personalized
Introduction," Lecture Notes of the Unione Matematica Italiana, Springer-
Verlag Berlin Heidelberg, 7, ISBN: 978-3-642-05194-4
[23] Mathworks, 2012, \Product Documentation," http://www.mathworks.com/
help/techdoc/math/f1-86528.html, Note: Online; accessed 30-August-2012
@ 09:05 PM, Pacic Standard Time (PST)
[24] Mori T., and Tanaka, K., 1997, \Average Stress in Matrix and Average Elastic
Energy of Materials with Mistting Inclusions," Acta Metallurgica, 21, pp.
571-574
163
[25] Okubo H., Kato K., Hayakawa N., Hanai M., and Takei M., 2009, \Function-
ally Graded Materials and Their Application to High Electric Field Power
Equipment," Performance of Conventional and New Materials for High Voltage
Apparatus, Cigre SC D1-Colloquium in Hungary, Budapest
[26] Pompe W., Worch H., Epple M., Friess W., Gelinsky M., Greil P., Hempele
U., Scharnweber D., and Schulte K., 2003, \Functionally Graded Materials
for Biomedical Application," Elsevier Science Limited, Materials Science and
Engineering, A362, pp. 40-60
[27] Qing-Hua Qin, and Qing-Sheng Yang, 2008, \Macro-Micro Theory on Multield
Coupling Behavior of Heterogeneous Materials," Higher Education Press,
Beijing, ISBN: 978-7-04-022350-7 and Springer Science+Business Media,
LLC, New York ISBN: 978-3-540-78258-2, pp. 7-56
[28] Singeresu S. Rao, 2011, \The Finite Element Method in Engineering," Elsevier
Inc., ISBN: 978-1-85617-661-3
[29] Thomas Reiter, George J. Dvorak, 1997, \Micromechanical Models for Graded
Composite Materials," Elsevier Science Ltd., J. Mech. Phys. Solids, 45, No. 8,
pp. 1281-1302
[30] M. H. Santare, and J. Lambros, 2000, \Use of Graded Finite Elements to
Model the Behavior of Nonhomogeneous Materials," ASME, Journal of Applied
Mechanics, 67, pp. 819-822
[31] Terada K., Miura T. and Kikuchi, N., 1997, \Digital Image-Based Modeling
Applied to the Homogenization Analysis of Composite Materials," Springer-
Verlag, Computational Mechanics, 20, pp. 331-346
[32] Kumar Vemaganti, and Pushkaraj Deshmukh, 2006, \An Adaptive Global-
Local Approach to Modeling Functionally Graded Materials," Springer-Verlag,
195, pp. 4230-4243
[33] Wo sko M., Paszkiewicz B., Piasecki T., Szyszka A., Paszkiewicz R., and
T lacza la M., 2005, \Applications of Functionally Graded Materials in Opto-
electronic Devices," Optica Applicata, XXXV, No. 3, pp. 663-667
164
[34] Young W. Kwon, David H. Allen, and Ramesh Talreja (Editors), 2008, \Mul-
tiscale Modeling and Simulation of Composite Materials and Structures,"
Springer Science+Business Media, LLC, ISBN: 978-0-387-36318-9, pp. 40-
41
[35] Wikipedia (Wikimedia Foundation), 2012, \Digital Image," http://en.
wikipedia.org/wiki/Digital_image, Note: Online; accessed 29-August-
2012 @ 10:05 PM, Pacic Standard Time (PST)
165
Appendix A
Generalised Micro-Solution
Variables for the Distorted Finite
Elements
In the preceding Chapter 3, micro-solution variables were derived by consider-
ing constant anti-periodic displacements like1 for both horizontal and vertical
directions. Although, constants like1 works for most of the regular shaped
quadrilaterals, it fails to produce a proper anti-periodic modes when the macro
element (also micro element) is distorted or rotated (w.r.t. any axis). In such a
case, anti-periodic displacements has to be changed according to the orientation of
the micro element and utilize the iso-parametric method.
X
1
, Y
1
X
Y
X
2
, Y
2
X
4
, Y
4
X
3
, Y
3
X
c
, Y
c

Q
Figure A.1: A Rotated Iso-parametric Finite Element
166
Figure A.1 shows an iso-parametric element rotated parallel to the XY plane, also
it should be noted that the element could be distorted.
From Eq. 3.43,

K
LJik
T
JKi
U
K
i
= 0 (A.1)
U
K
1
=X
K
X
c
(A.2)
Where, K = 1 4, number of micro elemental nodes (main exterior nodes).
U
1
1
=X
1
X
c
U
2
1
=X
2
X
c
U
3
1
=X
3
X
c
U
4
1
=X
4
X
c
9
>
>
=
>
>
;
(A.3)
Similarly,
U
K
2
=Y
K
Y
c
(A.4)
U
1
2
=Y
1
Y
c
U
2
2
=Y
2
Y
c
U
3
2
=Y
3
Y
c
U
4
2
=Y
4
Y
c
9
>
>
=
>
>
;
(A.5)
The transformation equation T
J11
is given by,
T
J11
=
0
B
B
B
@
1
Y
J
1
L
1
| {z }

1
(Y
1
)
+

1
J
L
1
1
C
C
C
A
0
B
B
B
@
1
Y
J
2
L
2
| {z }

1
(Y
2
)
+

4
J
L
2
1
C
C
C
A
(A.6)
Y
J
1
=
L
1
2

1 +
J
1

(A.7)
167
But,

J
1
=
2Y
1
L
1
1 (A.8)

J
2
=
2Y
2
L
2
1 (A.9)
Therefore,

1
(Y
J
1
) =
1
2

1
J
1

(A.10)

1
(Y
J
2
) =
1
2

1
J
2

(A.11)
Substituing Eqs. A.10 and A.11,
T
J11
=

1
2

1
J
1

+
1
J

1
2

1
J
2

+
4
J

(A.12)
T
J21
=

1
2

1 +
J
1

1
J

1
2

1
J
2

+
4
J

(A.13)
T
J31
=

1
2

1 +
J
1

1
J

1
2

1 +
J
2

4
J

(A.14)
T
J41
=

1
2

1
J
1

+
1
J

1
2

1 +
J
2

4
J

(A.15)
T
J12
=

1
2

1
J
1

+
3
J

1
2

1
J
2

+
2
J

(A.16)
T
J22
=

1
2

1 +
J
1

3
J

1
2

1
J
2

+
2
J

(A.17)
T
J32
=

1
2

1 +
J
1

3
J

1
2

1 +
J
2

2
J

(A.18)
T
J42
=

1
2

1
J
1

+
3
J

1
2

1 +
J
2

2
J

(A.19)
168
A.1 Basic Axial Micro Solution Variable in X Di-
rection
Consider the Eq. 3.43,

K
LJik
T
JKi
U
K
i
| {z }
i=k=1
= 0 (A.20)
T
JKi
U
K
1
=T
J11
U
1
1
+T
J21
U
2
1
+T
J31
U
3
1
+T
J41
U
4
1
(A.21)
=T
J11
(X
1
X
c
) +T
J21
(X
2
X
c
)
+T
J31
(X
3
X
c
) +T
J41
(X
4
X
c
) (A.22)
=
1
4

2
1
J

J
1
+ 1

2
4
J

J
2
+ 1

(X
1
X
c
)

2
1
J

J
1
1

2
4
J

J
2
+ 1

(X
2
X
c
)
+

2
1
J

J
1
1

2
4
J

J
2
1

(X
3
X
c
)

2
1
J

J
1
+ 1

2
4
J

J
2
1

(X
4
X
c
)

(A.23)
Let,
u
11
=X
1
X
c
u
12
=X
2
X
c
u
13
=X
3
X
c
u
14
=X
4
X
c
9
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
;
(A.24)
169
Substituting A.24 in A.23 and then in Eq. 3.43. Further simplifying,

K
LJ11
1
4

2
1
J
u
11
+ 2
4
J
u
11
+ 4
1
J

4
J
u
11
u
11

J
1

2
4
J

J
1
u
11
u
11

J
2
2
1
J

J
2
u
11
+
J
1

J
2
u
11
+u
11
2
1
J
u
12
+ 2
4
J
u
12
4
1
J

4
J
u
12
+u
12

J
1
+
2
4
J

J
1
u
12
u
12

j
2
+ 2
1
J

J
2
u
12

J
1

J
2
u
12
+u
12
2
1
J
u
13
2
4
J
u
13
+ 4
1
J

4
J
u
13
+u
13

J
1

2
4
J

J
1
u
13
+u
13

J
2
2
1
J

J
2
u
13
+
J
1

J
2
u
13
+u
13
+2
1
J
u
14
2
4
J
u
14
4
1
J

4
J
u
14
u
14

J
1
+
2
4
J

J
1
u
14
+u
14

J
2
+ 2
1
J

J
2
u
14

J
1

J
2
u
14
+u
14

= 0 (A.25)

K
LJ11
1
4

2
1
J
(u
11
u
12
u
13
+u
14
) + 2
4
J
(u
11
+u
12
u
13
u
14
)
+4
1
J

4
J
(u
11
u
12
+u
13
u
14
) +
J
1
(u
12
+u
13
u
11
u
14
)
+2
4
J

J
1
(u
12
+u
14
u
11
u
13
) +
J
2
(u
13
+u
14
u
11
u
14
)
+2
1
J

J
2
(u
12
+u
14
u
11
u
13
) +
J
1

J
2
(u
11
+u
13
u
12
u
14
)
+u
11
+u
12
+u
13
+u
14
] = 0 (A.26)

K
LJ11

2
1
J
(u
11
u
12
u
13
+u
14
) + 2
4
J
(u
11
+u
12
u
13
u
14
)
+4
1
J

4
J
(u
11
u
12
+u
13
u
14
) + 2
4
J

J
1
(u
12
+u
14
u
11
u
13
)+
+2
1
J

J
2
(u
12
+u
14
u
11
u
13
)

=

K
LJ11

J
1
(u
12
+u
13
u
11
u
14
) +
J
2
(u
13
+u
14
u
11
u
14
)
+
J
1

J
2
(u
11
+u
13
u
12
u
14
)u
11
u
12
u
13
u
14

(A.27)
170
Let,
C
1
= 2(u
11
u
12
u
13
+u
14
)
C
2
= 2(u
11
+u
12
u
13
u
14
)
C
3
= 4(u
11
u
12
+u
13
u
14
)
C
4
= 2(u
11
+u
12
u
13
+u
14
)
C
5
= 2(u
11
+u
12
u
13
+u
14
)
9
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
;
(A.28)
Let RHS be,

T
1
J
=

K
LJ11

J
1
(u
11
u
12
u
13
+u
14
) +
j
2
(u
11
+u
12
u
13
u
14
)

J
1

J
2
(u
11
+u
13
u
12
u
14
) (u
11
+u
12
+u
13
+u
14
)

(A.29)
Substituting A.28 in Eq. A.27 leads,
C
1

1
J
+C
2

4
J
+C
3

1
J

4
J
+C
4

4
J

J
1
+C
5

1
J

J
2
=

K
LJ11

1

T
1
J
(A.30)
171
A.2 Basic Axial Micro Solution Variable in XY
Direction
Consider the Eq. 3.43,

K
LJik
T
JKi
U
K
i
| {z }
i=1;k=2
= 0 (A.31)
T
JKi
U
K
1
=T
J11
U
1
1
+T
J21
U
2
1
+T
J31
U
3
1
+T
J41
U
4
1
(A.32)
=T
J11
(Y
1
Y
c
) +T
J21
(Y
2
Y
c
)
+T
J31
(Y
3
Y
c
) +T
J41
(Y
4
Y
c
) (A.33)
=
1
4

2
1
J

J
1
+ 1

2
4
J

J
2
+ 1

(Y
1
Y
c
)

2
1
J

J
1
1

2
4
J

J
2
+ 1

(Y
2
Y
c
)
+

2
1
J

J
1
1

2
4
J

J
2
1

(Y
3
Y
c
)

2
1
J

J
1
+ 1

2
4
J

J
2
1

(Y
4
Y
c
)

(A.34)
Let,
v
41
=Y
1
Y
c
v
42
=Y
2
Y
c
v
43
=Y
3
Y
c
v
44
=Y
4
Y
c
9
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
;
(A.35)
172
Substituting A.35 in A.34 and then in Eq. 3.43. Further simplifying,

K
LJ12
1
4

2
1
J
v
41
+ 2
4
J
v
41
+ 4
1
J

4
J
v
41
v
41

J
1

2
4
J

J
1
v
41
v
41

J
2
2
1
J

J
2
v
41
+
J
1

J
2
v
41
+v
41
2
1
J
v
42
+ 2
4
J
v
42
4
1
J

4
J
v
42
+v
42

J
1
+
+2
4
J

J
1
v
42
v
42

j
2
+ 2
1
J

J
2
v
42

J
1

J
2
v
42
+v
42
2
1
J
v
43
2
4
J
v
43
+ 4
1
J

4
J
v
43
+v
43

J
1

2
4
J

J
1
v
43
+v
43

J
2
2
1
J

J
2
v
43
+
J
1

J
2
v
43
+v
43
+2
1
J
v
44
2
4
J
v
44
4
1
J

4
J
v
44
v
44

J
1
+
2
4
J

J
1
v
44
+v
44

J
2
+ 2
1
J

J
2
v
44

J
1

J
2
v
44
+v
44

= 0 (A.36)

K
LJ12
1
4

2
1
J
(v
41
v
42
v
43
+v
44
) + 2
4
J
(v
41
+v
42
v
43
v
44
)
+4
1
J

4
J
(v
41
v
42
+v
43
v
44
) +
J
1
(v
42
+v
43
v
41
v
44
)
+2
4
J

J
1
(v
42
+v
44
v
41
v
43
) +
J
2
(v
43
+v
44
v
41
v
42
)
+2
1
J

J
2
(v
42
+v
44
v
41
v
43
) +
J
1

J
2
(v
41
+v
43
v
42
v
44
)
+v
41
+v
42
+v
43
+v
44
] = 0 (A.37)

K
LJ12

2
1
J
(v
41
v
42
v
43
+v
44
) + 2
4
J
(v
41
+v
42
v
43
v
44
)
+4
1
J

4
J
(v
41
v
42
+v
43
v
44
) + 2
4
J

J
1
(v
42
+v
44
v
41
v
43
)+
+2
1
J

J
2
(v
42
+v
44
v
41
v
43
)

=

K
LJ12

J
1
(v
42
v
43
+v
41
+v
44
) +
J
2
(v
41
+v
42
v
43
v
44
)
+
J
1

J
2
(v
42
v
43
+v
44
v
41
) (v
41
+v
42
+v
43
+v
44
)

(A.38)
173
Let,
C
6
= 2(v
41
v
42
v
43
+v
24
)
C
7
= 2(v
41
+v
42
v
43
v
44
)
C
8
= 4(v
41
v
42
+v
43
v
44
)
C
9
= 2(v
41
+v
42
v
43
+v
44
)
C
10
= 2(v
41
+v
42
v
43
+v
44
)
9
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
;
(A.39)
Let RHS be,

T
4
J
=

K
LJ12

J
1
(v
41
v
42
v
43
+v
44
) +
j
2
(v
41
+v
42
v
43
v
44
)

J
1

J
2
(v
41
+v
42
v
43
+v
44
) (v
41
+v
42
+v
43
+v
44
)

(A.40)
Substituting A.39 in Eq. A.38 leads,
C
6

1
J
+C
7

4
J
+C
8

1
J

4
J
+C
9

4
J

J
1
+C
10

1
J

J
2
=

K
LJ12

1

T
4
J
(A.41)
174
A.3 Basic Axial Micro Solution Variable in Y Di-
rection
Consider the Eq. 3.43,

K
LJik
T
JKi
U
K
i
| {z }
i=k=2
= 0 (A.42)
T
JKi
U
K
2
=T
J12
U
1
2
+T
J22
U
2
2
+T
J32
U
3
2
+T
J42
U
4
2
(A.43)
=T
J12
(Y
1
Y
c
) +T
J22
(Y
2
Y
c
)
+T
J32
(Y
3
Y
c
) +T
J42
(Y
4
Y
c
) (A.44)
=
1
4

2
2
J

J
2
+ 1

2
3
J

J
1
+ 1

(Y
1
Y
c
)

2
2
J

J
2
+ 1

2
3
J

J
1
1

(Y
2
Y
c
)
+

2
2
J

J
2
1

2
3
J

J
1
1

(Y
3
Y
c
)

2
2
J

J
2
1

2
3
J

J
1
+ 1

(Y
4
Y
c
)

(A.45)
Let,
v
21
=Y
1
Y
c
v
22
=Y
2
Y
c
v
23
=Y
3
Y
c
v
24
=Y
4
Y
c
9
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
;
(A.46)
175
Substituting A.46 in A.45 and then in Eq. 3.43. Further simplifying,

K
LJ22
1
4

2
2
J
v
21
+ 2
3
J
v
21
+ 4
2
J

3
J
v
21
v
21

J
1

2
2
J

J
1
v
21
v
21

J
2
2
3
J

J
2
v
21
+
J
1

J
2
v
21
+v
21
+2
2
J
v
22
2
3
J
v
22
4
2
J

3
J
v
22
+v
22

J
1
+
2
2
J

J
1
v
22
v
22

j
2
+ 2
3
J

J
2
v
22

J
1

J
2
v
22
+v
22
2
2
J
v
23
2
3
J
v
23
+ 4
2
J

3
J
v
23
+v
23

J
1

2
2
J

J
1
v
23
+v
23

J
2
2
3
J

J
2
v
23
+
J
1

J
2
v
23
+v
23
2
2
J
v
24
+ 2
3
J
v
24
4
2
J

3
J
v
24
v
24

J
1
+
2
2
J

J
1
v
24
+v
24

J
2
+ 2
3
J

J
2
v
14

J
1

J
2
v
24
+v
24

= 0 (A.47)

K
LJ22
1
4

2
2
J
(v
21
+v
22
v
23
v
24
) + 2
3
J
(v
21
v
22
v
23
+v
24
)
+4
2
J

3
J
(v
21
v
22
+v
23
v
24
) +
J
1
(v
22
+v
23
v
21
v
24
)
+2
2
J

J
1
(v
22
+v
24
v
21
v
23
) +
J
2
(v
23
+v
24
v
21
v
22
)
+2
3
J

J
2
(v
22
+v
24
v
21
v
23
) +
J
1

J
2
(v
21
+v
23
v
22
v
24
)
+v
21
+v
22
+v
23
+v
24
] = 0 (A.48)

K
LJ22

2
2
J
(v
21
+v
22
v
23
v
24
) + 2
3
J
(v
21
v
22
v
23
+v
24
)
+4
2
J

3
J
(v
21
v
22
+v
23
v
24
) + 2
2
J

J
1
(v
22
+v
24
v
21
v
23
)+
+2
3
J

J
2
(v
22
+v
24
v
21
v
23
)

=

K
LJ22

J
1
(v
22
v
23
+v
21
+v
24
) +
J
2
(v
21
+v
22
v
23
v
24
)
+
J
1

J
2
(v
22
v
23
+v
24
v
21
) (v
21
+v
22
+v
23
+v
24
)

(A.49)
176
Let,
C
11
= 2(v
21
+v
22
v
23
v
24
)
C
12
= 2(v
21
v
22
v
23
+v
24
)
C
13
= 4(v
21
v
22
+v
23
v
24
)
C
14
= 2(v
21
+v
22
v
23
+v
24
)
C
15
= 2(v
21
+v
22
v
23
+v
24
)
9
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
;
(A.50)
Let RHS be,

T
2
J
=

K
LJ22

J
1
(v
21
v
22
v
23
+v
24
) +
j
2
(v
21
+v
22
v
23
v
24
)

J
1

J
2
(v
21
+v
22
v
23
+v
24
) (v
21
+v
22
+v
23
+v
24
)

(A.51)
Substituting A.50 in Eq. A.49 leads,
C
11

2
J
+C
12

3
J
+C
13

2
J

3
J
+C
14

2
J

J
1
+C
15

3
J

J
2
=

K
LJ22

1

T
2
J
(A.52)
177
A.4 Basic Axial Micro Solution Variable in YX
Direction
Consider the Eq. 3.43,

K
LJik
T
JKi
U
K
i
| {z }
i=2;k=1
= 0 (A.53)
T
JKi
U
K
2
=T
J12
U
1
2
+T
J22
U
2
2
+T
J32
U
3
2
+T
J42
U
4
2
(A.54)
=T
J12
(X
1
X
c
) +T
J22
(X
2
X
c
)
+T
J32
(X
3
X
c
) +T
J42
(X
4
X
c
) (A.55)
=
1
4

2
2
J

J
2
+ 1

2
3
J

J
1
+ 1

(X
1
X
c
)

2
2
J

J
2
+ 1

2
3
J

J
1
1

(X
2
X
c
)
+

2
2
J

J
2
1

2
3
J

J
1
1

(X
3
X
c
)

2
2
J

J
2
1

2
3
J

J
1
+ 1

(X
4
X
c
)

(A.56)
Let,
u
31
=X
1
X
c
u
32
=X
2
X
c
u
33
=X
3
X
c
u
34
=X
4
X
c
9
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
;
(A.57)
Substituting A.57 in A.56 and then in Eq. 3.43. Further simplifying,
178

K
LJ21
1
4

2
2
J
u
31
+ 2
3
J
u
31
+ 4
2
J

3
J
u
31
u
31

J
1

2
2
J

J
1
u
31
u
31

J
2
2
3
J

J
2
u
31
+
J
1

J
2
u
31
+u
31
+2
2
J
u
32
2
3
J
u
32
4
2
J

3
J
u
32
+u
32

J
1
+
2
2
J

J
1
u
32
u
32

j
2
+ 2
3
J

J
2
u
32

J
1

J
2
u
32
+u
32
2
2
J
u
33
2
3
J
u
33
+ 4
2
J

3
J
u
33
+u
33

J
1

2
2
J

J
1
u
33
+u
33

J
2
2
3
J

J
2
u
33
+
J
1

J
2
u
33
+u
33
2
2
J
u
34
+ 2
3
J
u
34
4
2
J

3
J
u
34
u
34

J
1
+
2
2
J

J
1
u
34
+u
34

J
2
+ 2
3
J

J
2
u
34

J
1

J
2
u
34
+u
34

= 0 (A.58)

K
LJ21
1
4

2
2
J
(u
31
+u
32
u
33
u
34
) + 2
3
J
(u
31
u
32
u
33
+u
34
)
+4
2
J

3
J
(u
31
u
32
+u
33
u
34
) +
J
1
(u
32
+u
33
u
31
u
34
)
+2
2
J

J
1
(u
32
+u
34
u
31
u
33
) +
J
2
(u
33
+u
34
u
31
u
32
)
+2
3
J

J
2
(u
32
+u
34
u
31
u
33
) +
J
1

J
2
(u
31
+u
33
u
32
u
34
)
+u
31
+u
32
+u
33
+u
34
] = 0 (A.59)

K
LJ21

2
2
J
(u
31
+u
32
u
33
u
34
) + 2
3
J
(u
31
u
32
u
33
+u
34
)
+4
2
J

3
J
(u
31
u
32
+u
33
u
34
) + 2
2
J

J
1
(u
32
+u
34
u
31
u
33
)+
+2
3
J

J
2
(u
32
+u
34
u
31
u
33
)

=

K
LJ21

J
1
(u
32
u
33
+u
31
+u
34
) +
J
2
(u
31
+u
32
u
33
u
34
)
+
J
1

J
2
(u
32
u
33
+u
34
u
31
) (u
31
+u
32
+u
33
+u
34
)

(A.60)
179
Let,
C
16
= 2(u
31
+u
32
u
33
u
34
)
C
17
= 2(u
31
u
32
u
33
+u
34
)
C
18
= 4(u
31
u
32
+u
33
u
34
)
C
19
= 2(u
31
+u
32
u
33
+u
34
)
C
20
= 2(u
31
+u
32
u
33
+u
34
)
9
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
;
(A.61)
Let RHS be,

T
3
J
=

K
LJ21

J
1
(u
31
u
32
u
33
+u
34
) +
j
2
(u
31
+u
32
u
33
u
34
)

J
1

J
2
(u
31
+u
32
u
33
+u
34
) (u
31
+u
32
+u
33
+u
34
)

(A.62)
Substituting A.61 in Eq. A.60 leads to,
C
16

2
J
+C
17

3
J
+C
18

2
J

3
J
+C
19

2
J

J
1
+C
20

3
J

J
2
=

K
LJ21

1

T
3
J
(A.63)
180
Since the Eqs. pair A.30, A.41 and A.52, A.63 are coupled and non-linear in
nature. These equations can be solved by numerical methods like Newton's method,
Newton-Raphson method or any other approximate methods too. In this work,
Newton-Raphson method is utilized to nd the unknowns (
K
J
) in Eqs. A.30, A.52,
A.41 and A.63.
Considering Eqs. A.30 and A.41,
C
1

1
J
+C
2

4
J
+C
3

1
J

4
J
+C
4

4
J

J
1
+C
5

1
J

J
2
=

K
LJ11

1

T
1
J
| {z }

1
J
C
6

1
J
+C
7

4
J
+C
8

1
J

4
J
+C
9

4
J

J
1
+C
10

1
J

J
2
=

K
LJ12

1

T
4
J
| {z }

2
J
Let,

1
J
=
1
J
+
1
J
(A.64)

4
J
=
4
J
+
4
J
(A.65)
Substituting Eqs. A.64, A.65 in Eqs. A.30,
C
1
(
1
J
+
1
J
) +C
2
(
4
J
+
4
J
) +C
3
(
1
J
+
1
J
)(
4
J
+
4
J
)
+C
4
(
4
J
+
4
J
)
J
1
+C
5
(
1
J
+
1
J
)
J
2
=
1
J

C
1
+C
3

4
J
+C
5

J
2

1
J
+

C
2
+C
3

1
J
+C
4

J
1

4
J
=
1
J
C
1

1
J
C
2

4
J
C
3

1
J

4
J
C
4

J
1

4
J
C
5

J
2

1
J
9
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
;
(A.66)
Similarly, considering Eq. A.41 and using Eqs. A.64 and A.65,
181
C
6
(
1
J
+
1
J
) +C
7
(
4
J
+
4
J
) +C
8
(
1
J
+
1
J
)(
4
J
+
4
J
)
+C
9
(
4
J
+
4
J
)
J
1
+C
1
0(
1
J
+
1
J
)
J
2
=
2
J

C
6
+C
8

4
J
+C
10

J
2

1
J
+

C
7
+C
8

1
J
+C
9

J
1

4
J
=
2
J
C
6

1
J
C
7

4
J
C
8

1
J

4
J
C
9

J
1

4
J
C
10

J
2

1
J
9
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
;
(A.67)
rewriting Eqs. A.66 and A.67 in matrix form,
2
6
4
C
1
+C
3

4
J
+C
5

J
2
C
2
+C
3

1
J
+C
4

J
1
C
6
+C
8

4
J
+C
10

J
2
C
7
+C
8

4
J
+C
9

J
1
3
7
5
8
>
<
>
:

1
J

4
J
9
>
=
>
;
=
8
>
<
>
:

1
J
C
1

1
J
C
2

4
J
C
3

1
J

4
J
C
4

J
1

4
J
C
5

J
2

1
J

2
J
C
6

1
J
C
7

4
J
C
8

1
J

4
J
C
9

J
1

4
J
C
10

J
2

1
J
9
>
=
>
;
| {z }
residual
(A.68)
Eq. A.68 can be solved by iterative Newton-Raphson algorithm. Using,

1
(0)
J
=
4
(0)
J
= 0 (Initial conditions) (A.69)

1
(i+1)
J
=
1
(i)
J
+
1
(i+1)
J
(Iterative condition) (A.70)

4
(i+1)
J
=
4
(i)
J
+
4
(i+1)
J
(Iterative condition) (A.71)
Further, considering Eqs. A.52 and A.63,
C
11

2
J
+C
12

3
J
+C
13

2
J

3
J
+C
14

2
J

J
1
+C
15

3
J

J
2
=

K
LJ22

1

T
2
J
| {z }

3
J
C
16

2
J
+C
17

3
J
+C
18

2
J

3
J
+C
19

2
J

J
1
+C
20

3
J

J
2
=

K
LJ21

1

T
3
J
| {z }

4
J
182
let,

2
J
=
2
J
+
2
J
(A.72)

3
J
=
3
J
+
3
J
(A.73)
Substituting Eqs. A.72, A.73 in Eqs. A.52,
C
11
(
2
J
+
2
J
) +C
12
(
3
J
+
3
J
) +C
13
(
2
J
+
2
J
)(
3
J
+
3
J
)
+C
14
(
2
J
+
2
J
)
J
1
+C
15
(
3
J
+
3
J
)
J
2
=
3
J

C
11
+C
13

3
J
+C
14

J
1

2
J
+

C
12
+C
13

2
J
+C
14

J
1

3
J
=
3
J
C
11

2
J
C
12

3
J
C
13

2
J

3
J
C
14

J
2

1
J
C
15

J
2

3
J
9
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
;
(A.74)
Now considering Eq.A.63, using Eqs. A.72 and A.73,
C
16
(
2
J
+
2
J
) +C
17
(
3
J
+
3
J
) +C
18
(
2
J
+
2
J
)(
3
J
+
3
J
)
+C
19
(
2
J
+
2
J
)
J
1
+C
20
(
3
J
+
3
J
)
J
2
=
4
J

C
16
+C
18

3
J
+C
19

J
1

2
J
+

C
17
+C
18

2
J
+C
20

J
2

3
J
=
4
J
C
16

2
J
C
17

3
J
C
18

2
J

3
J
C
19

J
1

2
J
C
20

J
2

3
J
9
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
;
(A.75)
rewriting Eqs. A.74 and A.75 in matrix form,
2
6
4
C
11
+C
13

3
J
+C
14

J
1
C
12
+C
13

2
J
+C
15

J
2
C
16
+C
18

3
J
+C
19

J
1
C
17
+C
18

2
J
+C
20

J
2
3
7
5
8
>
<
>
:

2
J

3
J
9
>
=
>
;
=
8
>
<
>
:

3
J
C
11

2
J
C
12

3
J
C
13

2
J

3
J
C
14

J
1

2
J
C
15

J
2

3
J

4
J
C
16

2
J
C
17

3
J
C
18

2
J

3
J
C
19

J
1

2
J
C
10

J
2

3
J
9
>
=
>
;
| {z }
residual
(A.76)
183
By Newton-Raphson algorithm Eq. A.76 can be solved by iterative procedure,

2
(0)
J
=
3
(0)
J
= 0 (Initial conditions) (A.77)

2
(i+1)
J
=
2
(i)
J
+
2
(i+1)
J
(Iterative condition) (A.78)

3
(i+1)
J
=
3
(i)
J
+
3
(i+1)
J
(Iterative condition) (A.79)
Lastly, algorithm in matrix form follows as,
2
6
4
C
a
1
+C
a
3

m
i+1
J
+C
a
4

J
1
C
a
2
+C
a
3

n
i+1
J
+C
a
5

J
2
C
a
6
+C
a
8

m
J
+C
a
9

J
1
C
a
7
+C
a
8

n
J
+C
a
10

J
2
3
7
5
8
>
<
>
:

n
i+1
J

m
i+1
J
9
>
=
>
;
=
8
>
<
>
:

p
J
C
a
1

n
J
C
a
2

m
J
C
a
3

n
J

m
J
C
a
4

J
1

n
J
C
a
5

J
2

m
J

q
J
C
a
6

n
J
C
a
7

m
J
C
a
8

n
J

m
J
C
a
9

J
1

n
J
C
a
10

J
2

m
J
9
>
=
>
;
| {z }
residual
(A.80)
It should be noted that the residual is also evaluated at (i + 1).
184
Appendix B
Computer Implementation of the
Aperiodic Method
Aperiodic homogenization routine is developed by using the Mathworks
R
, Matlab
Program. In addition to the aperiodic routine, a digital image based meshing
tools has been implemented to carryout the laborious task of discretization of
sub-patches and retrieving the material properties at each micro-cell. Finally, all
these subroutines, combined in a single form produces the Homo.Aperiodic program.
In this appendix it is intended to provide a brief working of the Homo.Aperiodic
computer method. Thus three
owcharts are presented in this appendix for the
homogenization method, digital image processor (meshing tool) andHomo.Aperiodic
main program.
185
B.1 Flowchart of the Homogenization Program
START
Read the input data
(sub-patch coordinates,
Local2Global data, material
properties)
Process the input data and form
(joint restraints, equation number, column
height for banded matrices, loads if any and
others)
Subroutine Stiffness
(forms the stiffness matrix)
For total number of
micro elements in sub-patch
True
False
Loop on
Elements
Assembly of stiffness matrices
Find micro-solutions
(K
LJik
, K i’s and transformation
matrix. finally K
macro
)
Subroutine Homo.Constants
(find the find the strain transformation matrix (STM), next find
the constant strains (STM x anti-periodic displacements) , lastly
homogenized constants (Abar, Bbar and etc...) )
STOP
For total number of
Integration points
Loop on
Gauss Points
Figure B.1: A Typical Flowchart of the Aperiodic Homogenization Method
186
B.2 Flowchart of the Digital Imaging Program
START
Read the input data
(image data, sub-patch
coordinates, scaling
height, width, type of
fibers, size of mesh)
Process the input data
(image scaling, detect fibers, spacing, draw
the meshing lines)
STOP
Pixel Traverse
(traverse through each pixel and
record the image coordinate)
Convert back to physical
coordinates and record
them
Pixel Traverse
(traverse through close to four
corner pixel and the middle pixel of
the micro-cell)
Record the material
property
(assigning a unique
number for each material)
Storing
the
Data
(Hcoord,
matpat)
Figure B.2: A Typical Flowchart of the Digital Image Processing Tool of Homoge-
nization Method
187
B.3 Flowchart of the Main Aperiodic Program
START
Read the input data
(sub-patch coordinates,
Local2Global data, material
properties)
Process the input data and form
(joint restraints, equation number, column height for
banded matrices, loads if any and others)
Subroutine Stiffness
(forms the stiffness matrix)
For total number
of micro elements
in sub-patch
False
Loop on
Elements
Assembly of stiffness matrices
STOP
Digital image processor
(extract sub-patch, mesh, store
coordinates of the mesh and material
pattern)
Homogenization
(provides the homogenized
constants)
Colsol
(to inverse assembled stiffness and multiply by
foce to get displacements)
True
Pre-processor
Deformation
plot
Micro-
displacem
ents
Displacement
plots
For total number of
integration points
Loop on
Gauss
Points
Figure B.3: A Typical Flowchart of the Main Homogenization FE Method
188

Abstract (if available)

Abstract Composite materials are the well-known substitutes for traditional metals in various industries because of their micro-structural character. Micro-structures provide a high strength-to-weight ratio, which makes them suitable for manufacturing large variety of applications ranging from simple toys to complicated space/aircraft structures. Since, these materials are widely used in high performance structures, their stress/thermal analysis issues are of major concern. Due to the high degree of material heterogeneity, it is extremely difficult to analyze such structures. ❧ Homogenization (rigorous averaging) is a process that overcomes the difficulty of modeling each micro-structure. It replaces an individual micro-structure by an equivalent material model representation (unit cell). Periodic micro-structures appear in regular intervals throughout the domain, in contrast aperiodic micro-structures follows an irregular pattern. Further, this method bridges the analysis gap between micro and macro domain of the structures. In this thesis, Homogenization procedure based on anti-periodic displacement fields for aperiodic micro-structures and aperiodic boundary conditions are considered to model the constitutive material matrix. This work could be easily implemented with the traditional finite element packages. In addition, it eventually increases the convergence accuracy and reduces the high computational expenses. Different problems are analyzed by the implementation of digital image processing schemes for the extraction of a unit cell around the Gauss quadrature points and the mesh-generation. In the future, this research defines a new path for the analysis of any random heterogeneous materials by its ease of implementation and the state-of-the-art micro-structure material modeling capabilities and digital image based micro-meshing.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Multi-scale modeling of functionally graded materials (FGMs) using finite element methods

PDF

Vision-based and data-driven analytical and experimental studies into condition assessment and change detection of evolving civil, mechanical and aerospace infrastructures

PDF

Modeling of multiscale continuum-atomistic systems using homogenization theory

PDF

Studies into vibration-signature-based methods for system identification, damage detection and health monitoring of civil infrastructures

PDF

Computational fluid dynamic analysis of highway bridge superstructures exposed to hurricane waves

PDF

Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications

PDF

An energy method for earthquake resistant design of RC structures

PDF

Analytical and experimental methods for regional earthquake spectra and structural health monitoring applications

PDF

Vision-based studies for structural health monitoring and condition assesment

PDF

Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems

PDF

Studies into data-driven approaches for nonlinear system identification, condition assessment, and health monitoring

PDF

Design optimization under uncertainty for rotor blades of horizontal axis wind turbines

PDF

Effect of the air blast on glazing systems safety: structural analysis of window panel responses to air blast for mitigating the injuries from flying glass

PDF

Analysis of embedded software architecture with precedent dependent aperiodic tasks

PDF

Optimal design, nonlinear analysis and shape control of deployable mesh reflectors

PDF

Form finding and shape control of space deployable truss structures

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

Analytical and experimental studies of modeling and monitoring uncertain nonlinear systems

PDF

Development of composite oriented strand board and structures

PDF

Damage detection using substructure identification

Asset Metadata

Creator Aghalaya Manjunatha, Preetham (author)

Core Title Homogenization procedures for the constitutive material modeling and analysis of aperiodic micro-structures

School Viterbi School of Engineering

Degree Master of Science

Degree Program Civil Engineering (Structural Engineering)

Publication Date 10/19/2012

Defense Date 06/28/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag aperiodic structures,composites,finite element method,homogenization,image based meshing,multi scale,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Wellford, L. Carter (committee chair), Lee, Vincent W. (committee member), Masri, Sami F. (committee member)

Creator Email aghalaya@usc.edu,preethamam@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-105649

Unique identifier UC11289010

Identifier usctheses-c3-105649 (legacy record id)

Legacy Identifier etd-AghalayaMa-1256.pdf

Dmrecord 105649

Document Type Thesis

Rights Aghalaya Manjunatha, Preetham

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA