.4...
3.2:. new.
. .

 

 

 

 

 

1.. ....,._.:..=. .3
#33? .

 

«mam? éfimwufif : ., . .‘ #15:...

. Prim: , ..:...i («4.3;

 

5&7 ‘ UBRARY
MiCl‘ligdi'. State
University

 

 

This is to certify that the
dissertation entitled

DESIGN OF MATERIALS WITH SPECIAL DYNAMIC
PROPERTIES USING NEGATIVE STIFFNESS
COMPONENTS

presented by

Jitendra Prasad

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Mechanical Engineering

 

 

ill}

Major ProTé‘ssor’s Signature
“I H W ,L 200}

Date

 

MSU is an affirmative-action, equal-opportunity employer

-..-.-.--—.-.—.—.-.-----.—.-—.-.—---—.-.—.—.—.—.-o--..--.--..

g.-.-._.-.-.-.—.—.---.-.-.-.—.‘.-.-.-.-.—.-.-.-.-._._._.-.-._

 

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

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p'ICIRC/DateDue.indd-p.1

 

DESIGN OF MATERIALS WITH SPECIAL DYNAMIC PROPERTIES USING
NEGATIVE STIFFNESS COMPONENTS

By

J itendra Prasad

A DISSERTATION

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

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2007

ABSTRACT

DESIGN OF MATERIALS WITH SPECIAL DYNAMIC PROPERTIES USING
NEGATIVE STIFFNESS COMPONENTS

By

J itendra Prasad

This work presents design concepts to synthesize composite materials with special
dynamic properties, namely, materials that soften at high frequencies. A typical rubber-
like material hardens with frequency and a material which reverses this behavior will find
application in product design for vibration absorption such as automobile engine mounts.
Such dynamic properties are achieved through the use of a two-phase material that has
inclusions of a viscoelastic material of negative elastic modulus in a typical matrix phase
that has a positive elastic modulus. Possible realizations of the negative stiffness
inclusion phase are presented. One way to realize the negative stiffness phase is by using
a lattice containing bistable structures. A numerical homogenization technique is used to
compute the average viscoelastic properties of such composites. A methodology is
presented for the automatic design of such special materials using topology optimization
techniques. The method and the vibration-isolation properties of a composite material

designed with it are demonstrated through examples.

Neither the thief can steal it, nor can the king take it
Neither divided amongst brothers, nor too heavy to carry
The more you expend it, the more it increases
Knowledge is the prime wealth

- Subhashitani (The wise sayings in Sanskrit)

iii

Dedicated to all the teachers in my life

iv

ACKNOWLEDGMENTS

My foremost appreciation and gratitude are due to my advisor Professor Alejandro Diaz
for his invaluable guidance, teaching, help and support throughout my doctoral studies. I
will be ever grateful to him for the academic confidence and self-reliance that he
inculcated in me through his seamless belief in my abilities and his encouragement to me

to perform better than others.

Deepest appreciation and thanks are also extended to the other members of my PhD
guidance committee -— Professor Chang-Yi Wang, Professor Alan Haddow and Professor
Brian Feeny — for the time off their busy schedule to examine this work and provide their
invaluable suggestions. Special thanks to Professor Haddow for the discussions and

brainstorming sessions during the initial phase of my research.

I thank the department staff and the fellow students who have extended their helping
hands in coping with my day-to-day life as a graduate student — especially to the graduate
secretary Aida Montalvo for taking care of various paperworks, my advisor’s professorial
assistant Sara Murawa for helping me with the development of demos and one of my lab-

mates Amir Zamiri for the interesting discussions on the material science.

This work was sponsored by Toyota Motor Company. I gratefully acknowledge this

support.

TABLE OF CONTENTS

LIST OF FIGURES ........................................................................... viii
CHAPTER 1
INTRODUCTION ............................................................................ 1
1.1 Motivation and Background ................................................... l
1.2 Organization of the Dissertation ............................................... 5
CHAPTER 2
TWO PHASE COMPOSITE MATERIAL ................................................ 6
2.1 Mechanical Models of Viscoelastic Composites ............................ 7
2.2 Conditions for Stable Softening of Dynamic Modulus ..................... 15
2.2.1 Conditions for softening ............................................. 15
2.2.2 Conditions for stability ............................................... 17
2.2.3 Combined conditions for softening and stability ............... 20
2.2.4 Numerical example ................................................ 23
2.3 Estimation of Composite Material Parameters ............................. 30
2.3.1 Example phase properties .......................................... 32
2.4 Stability of Negative Stiffness Phase B ..................................... 34
2.5 Implementation of Negative Stiffness ....................................... 36
2.6 A Lumped System Realization of Phase B1 ................................. 40
2.7 Phase B as a Lumped System ................................................. 48
2.7.1 Stability of the lumped system .................................... 51
2.7.2 Example of phase B ................................................ 53
2.8 Realization of Phase B Using Tileable Bistable Structure ............... 54
CHAPTER 3
EXAMPLES ILLUSTRATIN G MATERIAL SOFT ENING AND IMPROVED
VIBRATION ISOLATION ............................................................... 56
3.1 Homogenization of Viscoelastic Properties ................................ 57
3.2 Example 1 ....................................................................... 58
3.2.1 Transmissibility analysis .......................................... 61
3.3 Example 2 ....................................................................... 66
3.4 Example 3 ....................................................................... 69
3.5 Example 4 ....................................................................... 73
CHAPTER 4
TOPOLOGY OPTIMIZATION OF THE COMPOSITE ............................... 76
4.1 Background ..................................................................... 77
4.2 The Optimization Problem .................................................... 79
4.3 Sensitivity Analysis ............................................................ 87
4.4 Solution Strategy ............................................................... 89
4.5 Examples ......................................................................... 92
4.5.] Optimization example 1 ............................................. 92

vi

4.5.2 Optimization example 2 ....... . ..................................... 99

4.5.3 Optimization example 3 ........................................... 105

4.6 Summary and Discussions ..................................................... 114
CHAPTER 5

SYNTHESIS OF TILEABLE BISTABLE STRUCTURES ........................ 116

5.1 Introduction .................................................................... 1 18

5.2 An Example of a Bistable Periodic Structure ............................... 121

5.3 The Analysis Model ............................................................ 124

5.4 Optimization of the Topology ................................................ 125

5.4.1 The ground structure and design variables ..................... 125

5.4.2 The objective function .............................................. 127

5.4.3 The optimization problem ......................................... 130

5.4.4 The solution scheme ............................................... 131

5.5 Examples ........................................................................ 134

5.5.1 Example 1 ............................................................ 134

5.5.2 Example 2 ............................................................ 140

5.6 Rubber Model ................................................................... 144

5.6.1 New optimization problem ........................................ 147

5.6.2 Example ............................................................ 148

CONCLUSIONS AND FUTURE WORK ................................................. 151

REFERENCES ............................................................................... 156

vii

LIST OF FIGURES

Figure 2.1. Two-phase composite material (dashed box shows the fundamental

cell) .............................................................................................. 7
Figure 2.2. Two types of mechanical models corresponding to a simple mechanical
mixture .......................................................................................... 8
Figure 2.3. Standard linear solid model of viscoelasticity ............................... 10

Figure 2.4. Spring-dashpot system corresponding to the series system in Fig. 2.2(e) 11
Figure 2.5. Spring-dashpot system with boundary conditions .......................... 12

Figure 2.6. Spring-dashpot system with boundary conditions .......................... 23

Figure 2.7. x1(t) , x2(t) and x3 (I) at the reference parameters and forcing
frequency 1 Hz .................................................................................

Figure 2.8. x1(t) , x2(t) and x3(t) at the reference parameters and forcing

frequency 2000Hz ............................................................................ 25
Figure 2.9. x3(t) with K3 =—0.24 and w=2000 Hz .................................... 26
Figure 2.10. x3 (t) with K3 = —0.39 and forcing frequency le ...................... 27
Figure 2.11. x3 (2) with K3 = —o.39 and forcing frequency 200on .................. 27
Figure 2.12. Absolute value (a), real part (b) and imaginary part (c) of the effective
stiffness (k; ) ................................................................................. 29
30

Figure 2.13. Dynamic stiffness vs. the forcing frequency with c2 / k2 = 0.02 .......

Figure 2.14. A potential arrangement of constituents B1 and 132 to form material B 36

Figure 2.15. A simple bistable structure ................................................... 36
Figure 2.16. A typical force-displacement diagram for a bistable structure .......... 38
Figure 2.17. A typical stiffness-displacement diagram for a bistable structure ....... 38

viii

Figure 2.18. Configuration with the maximum negative stiffness ..................... 39

Figure 2.19. Negative stiffness implementation .. ........................................ 40
Figure 2.20. Two-dimensional lattice of phase B1 ....................................... 41
Figure 2.21. Components of the effective elastic tensor of material B ................ 47
Figure 2.22. Two-dimensional lattice of phase B ......................................... 48
Figure 2.23 Inclusion of material B in matrix of material A ............................ 51
Figure 2.24. A lattice of material B with three nodes fixed .............................. 52
Figure 2.25. Material B as a periodic structure ........................................... 55

Figure 3.1. Representative cell characterizing the (periodic) mixture of materials A

and B ........................................................................................... 57
Figure 3.2. Components of the effective elastic tensor after mixing A and B ........ 60
Figure 3.3. Phase B as layered material ................................................... 61
Figure 3.4. Engine mount system ........................................................... 62
Figure 3.5. The effective complex modulus EC ((0) .................................... 63
Figure 3.6. Engine mount as spring-damper system .................................... 63

Figure 3.7. Spring stiffness (k(a))) and damping coefficient (5(0)) of the engine
mount .......................................................................................... 65

Figure 3.8. Transmissibility as a function of frequency ................................. 66

Figure 3.9. Components of the effective elastic tensor mixing materials A and B in
Example 2 ..................................................................................... 68

Figure 3.10. Components of the effective elastic tensor of material B in Example 3 70

Figure 3.11. Components of the effective elastic tensor mixing materials A and B
in Example 3 .................................................................................. 72

Figure 3.12. A simple square fundamental cell having square inclusion of phase B
in the matrix of phase A .............................. . ...................................... 73

ix

Figure 3.13. Plot of the entries of CE for example 4 ....................................

Figure 3.14. Transmissibility of the composite (solid line) as compared to phase A
(dashed line) for example 4 ..................................................................

Figure 4.1. Element densities as design variables .........................................
Figure 4.2. Angle of orientation of material B in an element as a design variable

Figure 4.3. The composite designed in example 1 (initial angle 90
degrees) ..........................................................................................

Figure 4.4 Plot of the target and achieved “can“ for the density problem of
example 1 (initial angle 90 degrees)... .. .. .

Figure 4.5. Iteration history for example 1 for the density problem of example 1
(initial angle 90 degrees) .....................................................................

Figure 4.6. Orientation of material B for example 1 (initial angle 90 degrees) .......

Figure 4.7. Plot of the target and achieved ||c1212|| for the angle problem of example
1 (initial angle 90 degrees) ...................................................................

Figure 4.8. Iteration history for the density problem of example 1 (initial angle 90
degrees) .........................................................................................

Figure 4.9. New initial angle for example 1 ...............................................
Figure 4.10. Optimal base cell for example 1 (initial angle 45 degrees) ...............

Figure 4.11. Plot of the target and achieved ”(31212“ for example 1 (initial angle 45
degrees)” .. .. . .. . .

Figure 4.12. Iteration history for example 1 (initial angle 45 degrees) ................

Figure 4.13. Geometry of the base cell designed in example 2 (initial angle 90
degrees) .........................................................................................

Figure 4.14. Plot of the target and achieved dynamic moduli for the density
problem of example 2 (initial angle 90 degrees) ...........................................

Figure 4.15. Iteration history for the density problem of example 2 (initial angle 90
degrees) .........................................................................................

Figure 4.16. Optimal angle of orientation for material B for example 2 (initial
angle 90 degrees) ..............................................................................

74

75

80

82

94

94

94

95

96

96

97

97

98

98

100

100

100

101

Figure 4.17. Plot of the target and achieved dynamic moduli for the angle problem
of example 2 (initial angle 90 degrees) ............................... . ......... . ...........

Figure 4.18. Iteration history for the angle problem of example 2 (initial angle 90
degrees) .........................................................................................

Figure 4.19. Geometry of the base cell designed in example 2 (initial angle 45
degrees) .........................................................................................

Figure 4.20. Plot of the target and achieved dynamic moduli for the density
problem of example 2 (initial angle 45 degrees) ..........................................

Figure 4.21. Iteration history for the density problem of example 2 (initial angle 45
degrees) .........................................................................................

Figure 4.22. Optimal orientation of material B for example 2 (initial angle 45
degrees) .........................................................................................

Figure 4.23. Plot of the target and achieved dynamic moduli for the angle problem
of example 2 (initial angle 45 degrees) .....................................................

Figure 4.24. Iteration history for the angle problem of example 2 (initial angle 45
degrees) .........................................................................................

Figure 4.25. Geometry of the base cell designed in example 3 (initial angle 90
degrees) .........................................................................................

Figure 4.26. Plot of the target and achieved dynamic moduli for the density
problem of example 3 (initial angle 90 degrees) ..........................................

Figure 4.27. Iteration history for the density problem of example 3 (initial angle 90
degrees) .........................................................................................

Figure 4.28. Optimal orientation of material B for example 3 (initial angle 90
degrees) .........................................................................................

Figure 4.29. Plot of the target and achieved dynamic moduli for the angle problem
of example 3 (initial angle 90 degrees) .....................................................

Figure 4.30. Iteration history for the angle problem of example 5 (initial angle 90
degrees) .........................................................................................

Figure 4.31. Geometry of the base cell designed in example 3 (initial angle 45
degrees) .........................................................................................

xi

102

102

103

103

103

104

105

105

106

107

107

108

109

109

110

Figure 4.32. Plot of the target and achieved dynamic moduli for the density
problem of example 3 (initial angle 45 degrees) ..........................................

Figure 4.33. Iteration history for the density problem of example 6 (initial angle 45
degrees) .........................................................................................

Figure 4.34. Optimal orientation of material B for example 3 (initial angle 45
degrees) .........................................................................................

Figure 4.35. Plot of the target and achieved dynamic moduli for the angle problem
of example 3 (initial angle 45 degrees) .....................................................

Figure 4.36. Iteration history for the angle problem of example 3 (initial angle 45
degrees) .........................................................................................

Figure 5.1. Strain energy versus strain in a typical bistable structure ................
Figure 5.2. Simplified “double curved beam” bistable mechanism ...................

Figure 5.3. A bistable structure (first stable configuration) based on four “double
curved beam” sub-structures ................................................................

Figure 5.4. A sub-structure of Fig. 5.3 that is a bistable mechanism ..................

Figure 5.5. A bistable structure (second stable configuration) based on four
“double curved beam” sub-structures ......................................................

Figure 5.6. A 3x3 periodic arrangement of bistable structures (a) first equilibrium
configuration (b) second equilibrium configuration ......................................

Figure 5.7. A member (bar) of the ground structure .....................................
Figure 5.8. Typical force-displacement diagram for a bistable structure .............
Figure 5.9. The external load factor as a function of time ...............................

Figure 5.10. Example 1. The package space (a) and corresponding ground
structure (b) with dimensions and boundary conditions .................................

Figure 5.11. The force-displacement diagram for the reference structure (Fig. 5.3)

Figure 5.12. Example 1. Bistable structure obtained by the GA (a) first stable
configuration (b) second stable configuration ............................................

Figure 5.13. The force-displacement diagram for the structure in Fig. 5.12 ..........

xii

111

111

112

113

113

120

122

123

123

123

123

127

128

130

134

135

137

137

Figure 5.14. Example 1. Another bistable structure (a) first stable configuration
(b) second stable configuration ..............................................................

Figure 5.15. Force-displacement diagram for the bistable structure obtained by the
GA ...............................................................................................

Figure 5.16. History of the best GA merit function value f ..............................

Figure 5.17. Example 2. The package space (a) and corresponding ground
structure (b) with dimensions and boundary conditions .................................

Figure 5.18. Example 2. Bistable structure obtained by the GA (a) first stable
configuration, (b) second stable configuration ............................................

Figure 5.19. Example 2. Plot of the internal force at the input port A versus the
displacement at the output port B ...........................................................

Figure 5.20. Detail of the solution of Example 2. (a) First stable configuration
(b) Second stable configuration .............................................................

Figure 5.21. Rubber model: A member (bar) of the ground structure .................

Figure 5.22. Prescribed normal stress as a fimction of normal strain in the soft
material ..........................................................................................

Figure 5.23. Bistable structure obtained by the GA (a) first stable configuration
(b) second stable configuration ..............................................................

Figure 5.24. The force-displacement diagram for the structure in Fig. 5.23 ..........

xiii

138

139

140

140

142

143

144

145

146

149

149

CHAPTER 1

INTRODUCTION

1.1 Motivation and Background

Special dynamic properties such as a softening of the dynamic modulus of a viscoelastic
material at high frequencies are desirable for some applications. For example, a material
that softens at high frequencies may be used to make an ideal automobile engine mount,
as an ideal engine mount should maintain high stiffness at a low fi'equency and low
stiffness at higher frequencies (Y 11 et a1 (2000)). Such material may be called a smart
material, as smart materials typically respond to environmental stimuli (e.g., forcing
frequency) with particular changes in some variables (e. g., dynamic modulus). This work
explores concepts for the design of materials with these unusual dynamic properties, with

applications in design of engine mounts and similar vibration-isolation components.

The function of an engine mount is to support the weight of the engine while isolating
unbalanced engine disturbance forces from the vehicle structure. The force transmitted to
the vehicle structure increases with the dynamic stiffiress of the mount and therefore, the
mount should have low dynamic stiffness in order to isolate the vibration caused by the
engine. However, large static and quasi-static engine displacements resulting from low
dynamic stiffness and shock excitation (e.g., low-frequency excitation forces caused by
an uneven road or by sudden acceleration, deceleration or braking) may cause damage to

the engine and vehicle structure. The dynamic stiffness of an ideal engine mount,

therefore, should be a function of frequency of excitation. In other words, the dynamic
modulus of a material ideal for an engine mount should decrease with the forcing
frequency. However, a typical single-phase material does not soften at high frequencies.
In fact, a composite mixture of such typical materials does not soften either. For example,
a Voigt composite (i.e. a two-phase layered composite material in which laminae are
aligned in the direction of the force) will soften if and only if the elastic modulus of at
least one of the constituent phases decreases with fiequency. This is because the effective
elastic modulus of this composite is a convex combination of the elastic moduli of the
individual phases. Similarly, a Reuss composite (i.e. a two-phase layered composite
material in which laminae are aligned perpendicular to the direction of the force) cannot
display frequency-induced softening unless at least one of the constituent phases too
displays such softening. However, it is possible to achieve frequency-induced softening
in a two-phase composite if one of the constituent phases has a negative elastic modulus.
For example, in Wang and Lakes (2004a and 2004b) an inclusion of a ‘negative—stifi‘ness’
phase in a matrix of a typical material is shown to lead to a softening of the dynamic
stiffness at high forcing frequencies. This behavior is the basis for the concepts explored
here. The objective of this work is to generate design concepts that will lead to the
synthesis of two-phase composite materials that soften monotonically with forcing
frequency. The goal is to expose novel concepts and algorithms that could be used as
guidelines by material scientists in future work, rather than to provide a specific recipe to
synthesize the material. The work of Wang and Lakes (2004a and 2004b) serves as an

inspiration to use negative stiffiress inclusions.

A negative stiffness material has a negative elastic modulus. Such materials are not stable
and therefore, do not permanently exist in the negative stiffness state. A block of negative
stiffness material, however, can be made stable by constraining it from all sides by
surrounding the block with a typical (i.e. stable) material, as shown in Lakes and Drugan
(2002). Negative stiffness materials are realizable. For example, Lakes and Drugan
(2002) and Wang and Lakes (2004b, 2005b) proposed lumped lattices to implement
negative sfiffiress inclusions. In such lumped lattices, negative stiffness is basically
achieved by a bistable structure. A bistable structure has two stable configurations under
no external loads and is known to have ‘negative stiffiress’ in the neighborhood of third,

unstable configuration.

The effect of a negative stiffness inclusion in a composite material may be studied
analytically or numerically by approximating the composite as a lumped system. For
example, Wang and Lakes (2004a, 2004b) studied a one-dimensional spring-damper
system as an approximation of a two—phase composite material. In that work, numerical
analysis of the lumped system revealed atypical properties of a composite material having
a negative stiffness inclusion (the negative stiffness inclusion was demonstrated also to
lead to extremely high stiffness and damping). Analysis of the equivalent 1D lumped
model helps in choosing properties of matrix and inclusion phases that can lead to the

desired softening. This approach will be used here too.

With a negative stiffness material at hand and the knowledge that the inclusions of a

negative stiffness phase in a matrix of a typical material phase can lead to frequency-

induced softening, it is possible to tailor a composite material to desirable properties. The
so-called ‘inverse homogenization problem’ in which the effective properties of a
composite material is prescribed and the goal is to find topology that gives the prescribed
material properties, is well-known, for example in Sigmund (1995) and Diaz and Benard
(2003). The material design methodology has also been extended to design a viscoelastic
material, for example in Y1 et al. (2000). In this work, the existing knowledge in the
material design problem has been extended to design the composite materials that exhibit

frequency-induced softening.

The application of the present work is not limited to the frequency-induced softening or
the vibration isolation only; this work may be easily extended to design extremal
materials, where an extremal material refers to one having an elastic / viscoelasic
property equal to the maximum or minimum value allowed by the rigorous theoretical
bounds such as Hashin-Shtrikman bound (originated from Hashin and Shtrikman (1963)),
Cherkaev-Gibiansky bound (detailed in Cherkaev and Gibiansky (1993)) and the
corresponding viscoelastic bounds (given in Gibiansky and Milton (1993), Milton and
Berryman (1997) and Gibiansky et al. (1999)). The inverse homogenization problem has
already been used to design elastic materials with unusual or etxremal properties such as
the materials with zero or negative Poisson’s ratio as in Sigmund (1995) and the materials
having extremal bulk modulus as in Sigmund (2000). Unususal or extremal materials in
the context of viscoelasticty have also been made. For example, using materials with
negative bulk modulus, Jaglinski et al (2007) have built composite materials that have

viscoelastic stiffness greater than diamond. The negative stiffiress phase and the topology

optimization method used in this work may be extended to design materials that can

exhibit such extreme viscoelastic behaviors.

1.2 Organization of the Dissertation

The rest of the dissertation is organized as follows. A methodology to conceptually
design the desired material is described in chapter 2. It introduces a model of a two-phase
composite material and discusses a negative stiffness material, which is used to build the
inclusion phase. The stability and realization of the inclusion phase are discussed in this
chapter. Chapter 3 describes the homogenization method used here to compute the
effective properties of the composite material. The application of the design methodology
and the performance of the designed material are demonstrated through a few examples
in that chapter. The automatic design of the softening materials using topology
optimization is described in chapter 4 and examples are presented to evaluate the
performance of the design method. Chapter 5 gives details on designing tileable bistable
structures which may be used to realize the negative stiffness material. Conclusions are

drawn at the end and references are given.

CHAPTER 2

TWO PHASE COMPOSITE MATERIAL

With the ultimate goal of achieving a composite material that exhibits frequency-induced
softening, this chapter simplifies a two-dimensional composite material into a one
dimensional mechanical model that is made of springs and dampers. The mechanical
model is analyzed for a stable frequency-induced softening. It is found that one of the
springs is required to have a negative stiffness in order to get frequency-induced
softening. It is therefore important to allow the spring stiffness to assume a negative
value. A structure having a negative stiffness component, however, may be unstable and
therefore conditions for the stability of the mechanical model are derived to make sure
the model has a stable frequency-induced softening. Once the parameters (i.e. spring
stiffness and damper coefficients) of the 1D model leading to a stable frequency—induced
softening are determined, these parameters are extended back to construct the
corresponding two-dimensional composite model. The 1D system studies imply that a
composite material softens with frequency if the composite has inclusion of a negative
stiffness material B in a matrix of a typical elastic material A. The stability of the
negative stiffness inclusion is reviewed as the stability criterion of the 1D model may not
be sufficient in two-dimensions. Material A, being a typical elastic material with positive
Young’s modulus, is easily available. In contrast, material B has negative Young’s
modulus and is not available. A negative stiffness material that can be used as material B

is, therefore, built. Material B is proposed as a lumped lattice that has negative-stiffness

springs. The negative stiffness springs realizes the negative Young’s Modulus. The

negative stiffness is in turn realized by bistable structures.

2.1 Mechanical Models of Viscoelastic Composites

The strategy to design a material exhibiting frequency induced softening relies on the
study of non-homogenous materials with a periodic micro-structure designed to reverse
the behavior of typical materials, which stiffen with frequency. In particular, this work
studies a two-phase composite material composed of periodic inclusions of a viscoelastic
phase B in matrix of an elastic/viscoelastic phase A. A schematic arrangement of such
composites is shown in Fig. 2.1. The dashed square in Fig. 2.1 shows a firndamental cell
of the periodic arrangement. The properties of the constituent phases and their shape and

topology are the variables to be determined in order to obtain the desired results.

 

aoooo

oofipo

0.000

Figure 2.1. Two-phase composite material (dashed box shows the fimdamental cell)

 

 

 

Following Fuj ino et al. (1964) and Marinov (1978), the two phase composite material can
be approximated as two types of mechanical models combining the two phases as shown

in Figures 2.2(b,c) and 2.2(d,e). The former mechanical model represents a simple

additivity of contribution of partial stress of each element sliced vertically to the total
stress, and the latter that of partial strain of each element sliced horizontally to the total
strain. In these models, rigid adhesion between the two phases and no interference
between the sliced elements are assumed, i.e., actual stress or strain distribution along the

spherical surface of each phase is much simplified.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(e)

 

Figure 2.2. Two types of mechanical models corresponding to a simple mechanical
mixture, B sphere in A unit cubic lattice - (a) Unit or fundamental cell, (b) & (c)

equivalent parallel model, ((1) & (e) equivalent series model

In the model in Fig. 2.2(e), 21.1 and (I); are the length fractions of phase B in horizontal and
vertical directions, respectively. If both the two phases A and B are perfectly elastic

materials, the average Young’s modulus is obtained (as a function of the Young’s moduli

of the individual phases) by solving equations of equilibrium. This expression for the
average Young’s modulus can be extended to obtain the average complex modulus of the
viscoelatic material by using the correspondence principle. In the present context, the
correspondence principle states that the expression for the average complex modulus of a
viscoelastic composite may be obtained by replacing the phase Young’s moduli by phase

complex moduli in that expression of the average Young’s modulus.

Haddad (1995) presents a general definition of the correspondence principle as follows:
For a large number of technical viscoelastic problems, it is possible to relate
mathematically the solution of a linear, viscoelastic boundary value problem to an
analogous problem of an elastic body of the same geometry and under the same initial
and boundary conditions. This is carried out by transforming the governing equations of
the viscoelastic problem to be mathematically equivalent to those governing a
corresponding elastic problem. In this, both Laplace and Fourier transforms are often
used. Accordingly, one would be able to employ the tools of the theory of elasticity to
solve different boundary value problems in linear viscoelasticity. This analogy is referred
to as the ‘correspondence principle’ and implies the elastic procedures may be utilized to

derive transformed viscoelastic solutions.

Applying the correspondence principle, if E; and E; are complex moduli of phases A

and B, respectively, then the average complex modulus of the two-phase material may be

given by

 

—1
Efll =(1-41)E:1+41[(1—?‘)+flr] (2-1-1)
EA EB

Similarly for the model in Fig. 2.2(e), 7t; and q): are the length fractions of phase B in
horizontal and vertical directions, respectively. The average complex modulus based on

the series model is given by

—l

. (1-¢2) e
E = 2.1.2
”2 £11; +<1—22)EZ+22£}; ( )

Viscoelastic materials are often represented by mechanical models consisting of elastic
springs and viscous dashpots, where the elastic springs describe for the elastic behavior
and the dashpots describe viscous behavior. The mechanical model shown in Fig. 2.2(e)
will be studied. In this example, Phase A is assumed perfectly elastic and, therefore, will
be represented by an elastic spring. The inclusion phase B is modeled as a standard linear
solid. A standard linear solid model of viscoelasticity is sketched in Fig. 2.3, where E, E2

and 77 are the three parameters of the standard linear solid.

Figure 2.3. Standard linear solid model of viscoelasticity

10

The complex modulus corresponding to the viscoelastic model shown in Fig. 2.3 is given
by

:1: E E +i (0
153(0)): 1( 2 77 )
(E1+E2+ma))

 

(2.1.3)

The mechanical model corresponding to Fig. 2.2(e) is shown in Fig. 2.4. This mechanical

model is same as that studied in Wang and Lakes (2004a, 2004b, 2005a).

Figure 2.4. Spring-dashpot system corresponding to the series system in Fig. 2.2(e)

The parameters in Figs. 2.2(e) and 2.3 are related as follows:

EA: k1¢2 =k4ll‘fz) (2.1.4)
(Hod d

and

Eat: _ E1(E2 +i7760) _ [(3 (k2 +iC20) $2
B (E1 + E2 +i77w) (k3 + k2 + iCzW) 112d

 

(2.1.5)

since the viscoelastic material parameters are

11

E1 = k. 32— (2.1.6)

- 22d
52 = k2 2% (2.1.7)
I] = CZ 122% (2~1-8)
where

d 2 ¥ (2.1.9)

W is the width of the fundamental cell ( Fig. 2.2(a) ).
H is the height of the fundamental cell ( Fig. 2.2(a) ).

D is the thickness of the fundamental cell.

F, X3
k4

_IX1
CZ
I
_1X2
k3
Figure 2.5. Spring-dashpot system with boundary conditions

The lower end of the system will be fixed as shown in Fig. 2.5 and a force will be applied

at the upper end. The equation of motion of the system in Fig. 2.5 is given by

12

C2 —Cz 0 1"] k1 + k2 + k4 —k2 —k4 x1 0
-Cz 02 0 i2 + —k2 k2 + [(3 0 1'2 2 O
O 0 0 13 —k4 0 [(4 X3 F

The applied force F is periodic. F is given by

F=de

0

The steady-state solution to the equations of motion is given by

3‘1 x10 .
)52 — X20 81
x3 1‘30

Substituting (2.1.12) in the equation of motion (2.1.10), we get

[(1 + k2 + k4 +11“? -k2 -i0X’2 —k4 x10 0
—k2 —i((X'2 [(2 +113 +iwc2 0 X20 = 0
-k4 0 k4 3‘30 F0

Solving the new equation ofmotion (2.1.13), we get

0 = k4[k3(k2 + i62w)+k1(k2 +k3 +iC2€0)]

X30 [(3 (k2 '1' i620) + (k1 + [(4 )(k2 + k3 + iCzCU)

F

 

Eliminating k1, k2 , k3, k4 and c; in (2.1.14) by using (2.1.3)-(2.1.6), we get

—1
fa: (1—¢2)+ “)2 dzEgzd

x30 E2 (1— I12)E:1+ 4251;

13

(21w)

01H)

(21w)

019)

min)

01w)

where E22 is the average complex modulus for the two-phase material (equation 2.1.2).

The quantity Ezrzd is actually the average complex stiffness (kg) of the firndamental

cell in Fig. 2.2(a). So, we define

k;,=-5-=E;,2d (2.1.16)
x30
or,
F
k}; = 0 _ 1334.1 (2.1.17)
lesoll

 

 

 

 

 

 

 

"kg" and ”E1122" will be referred to as the dynamic stiffness and dynamic modulus,

respectively.

It should be noted that the parameters k1, k2 , k4 and C2 in the example system (Fig. 2.5)
will be treated as positive numbers, while k3 could assume a negative value, when desired.
Bistable structures are known to have negative stiffiiess (e.g. in Wang and Lakes (2004a,
2004b, 2005b)). A negative k3 can be implemented via a bistable structure. The reason to
allow k3 to be negative is to achieve a dynamic stiffness/modulus that softens with

frequency.

14

2.2. Conditions for Stable Softening of Dynamic Modulus

2.2.1 Conditions for softening

From (2.1.17), it can be seen that dynamic stiffness and dynamic modulus are
proportional to each other. For convenience, softening of dynamic stiffiiess of the
example system in Section 2.1.1 will be investigated. We have (from equations 2.1.14

and 2.1.16),

*

=|| k4[k3(k2 +ic2w)+kl(k2 +k3 Hogan] ||
"k3(k2 41.6260) '1' (k1 + k4 )(k2 + k3 + iCzw)"

 

 

 

 

 

(2.2.1)

8 k},
By analyzing (2.2.1), we find that —B_ < 0 if and only if
(0

 

 

 

 

k4[2k12(k2 +k3)+k3 {k3k4 +2k2(k3 +k4)}+2k1{k3(k3 +k4)+k2(2k3 +k4)}]<0.

Since k2 > 0,

 

*
a In II
is negative provided that
(0

a

K4[2K12(1+K3)+K3 {K3K4 +2(1<3 +1<4)}+21<1 {K3(K3 +K4)+(2K3 +K4)}]<O,

where we define

Kl =k— (2.2.2)
2

K3 =l’:_3 (2.2.3)
2

K4 = k—4 (2.2.4)
k2

15

II!

8k

 

 

 

 

m

Since K4 > 0 , we get that a necessary condition for < 0 is

 

8a)

2K12(1+K3)+K3[K3K4 +2(1<3 +K4)]+2K1 [1r<3(1<3 +1<4)+21<3 +K4]<0

We now define the firnction w as

wtKr,K3,Kr>=2K3<1+K3)+K31K3K4 +2<K3 +K4)l+2K1[K3(K3 +K4)+2K3 +K41

(2.2.5)

8k

 

 

*
H
80)

 

 

 

It follows that is negative if and only if t/I(K1,K3 , K 4) < 0.

Since k1 and k4 are made from same material, K4 > 0 also implies K1 > 0. Since both

a:
alkflll
K] and K4 are positive, the requirement ——<0 is satisfied whenever

I
an)

K3min < K3 < K3max , where

 

—21<1 — K,2 -— K4 - K1K4 -\/K14+ 2K13K4 + K} + KEK}
K3min = (22-6)
2+2K1 +K4

 

-2K1 — K12 -— K4 - K1K4 +\/K14+ 2K13K4 + K} + [(1sz
K3max = (2.2.7)
2 + 2K1 + K4

 

K311“,x may be written as

16

 

~2K1—K12—K4 -—K1K4 +\[(K12+K1K4)2 +K}
2+2K1 +K4

 

Km,1x = (2.2.8)

Since both K1 and K4 are positive, K3min is negative and since

 

\/(K12 “(11(4)2 +K} < (K12 +1<11<4)+1<4 for any positive K1 and K4, K3max is also

81"?! ll
negative. Thus, the condition for the frequency-induced softening (i.e., -a—-<0) is
(l)

K3min < K3 < K3max < 0 -

Thus a negative stiffness spring is needed to achieve frequency softening.

2.2.2 Conditions for stability

In this section, the stability of the system in Fig.2.5 will be examined for a static force F.

Equation (2.1.10) can be rewritten as the following three equations

C2X.l -C2X2 +(kl +k2 +k4)xl —k2x2 — [(413 = 0 (2.2.9)
—62X1+62X2 - k2x1+(k2 + [(3 )X2 = 0 (2.2.10)
—k4xl +k4X3 = F (2.2.11)

Adding (2.2.9) and (2.2.10), we get

(k1 + k4)x1+k3x2 -k4X3 = 0 (2.2.12)

17

From (2.2.11) and (2.2.12) for k3 ¢ 0 and k4 at 0,

 

k —F
x] zL (2.2.13)
k4
k +k —k k
x :F(1 4) 14x3 (2.2.14)
k3k4
Differentiating (2.2.13) and (2.2.14),
5,1: 5,, (2.2.15)
2'2 = -k—li’3 (2.2.16)
k3 ‘

Substituting (2.2.15) and (2.2.16) in (2.2. 10), we get

02 (k1 + k31k45c3 + [k2k3 + k1(k2 + k3llk4x3 - Flk1(k2 + k3)+ k3k4 + k2 (k3 + k4 )1 =

0
WM

(2.2.17)

or,

02(k1 + k3)k45‘3 + [k2k3 + k1(l‘2 + k3)lk4x3 — Flk1(k2 + k3)+ k3k4 + k2(k3 + k4)] = 0

 

 

(2.2.18)
Solving the differential equation (2.2.18) we get
x3 = (1-eA’) Flk1(k2 + k3) +k3k4 “(20% + k4)] (2.2.19)
[k2k3 + k1(kz + k3)]k4
where
A:_k2k3+kl(k2+k3) (2220)

02(k1+k3)

18

We define

k : [k2k3 + k1(k2 +k3)]k4
eff [k1(k2 + k3)+ k3k4 +k2(k3 + k4 )]

 

(2.2.21)

keff is actually the static stiffness of the system. From (2.2.19) and (2.2.20), we find that
x3 is bounded if A < 0. Furthermore the effective stiffness of the system (i.e. kefl) must

be positive so that the external force will do positive work. A negative work done by the

external work is an indication of instability. In other words, the system is stable if

A < O and kef > 0 (2.2.22)

Since c2 is positive (and real), the conditions for the stability of the system is

K3+K1(1+K3)
K1+K3

 

> 0 and keff > 0. The first condition (i.e. A < 0) is satisfied if any one

of the following two conditions is satisfied:

(1)K3 +K1(I+K3)>O and K1+K3 >0, (2223)

(11) K3 +K1(I+K3)<O and K1+K3 < 0 (2.224)

Since K1 > 0 and K3 < 0 , the system is stable (i.e. x3 is bounded) if and only if kefl > 0

and any of the following two conditions is satisfied:

(i) K3 +K1(1+K3) > 0, (2.2.25)

19

2.2.3 Combined conditions for softening and stability
From section 2.2.1 we know that the dynamic stiffness of the example system decreases

with increase in the forcing frequency if and only if (11 < 0, whereas in section 2.2.2 we

learned that the example system is stable if and only if A < 0 and keff > 0. The goal of

this section is to investigate whether both the softening and stability conditions can be
satisfied simultaneously. Although there are two stability conditions to be satisfied

simultaneously (viz. A< 0 and keff >0), only the first stability condition (i.e. A<0)

will be applied to test whether the example system in Fig. 2.5 is stable in the frequency-

induced softening regime. The second stability condition (i.e. keff >0) can be

numerically checked later, once the two conditions w< O and A < 0 are simultaneously

satisfied.

As given in section 2.2.2, A < 0 if and only if any of the following two conditions is

satisfied:
(1)K3 +K1(I+K3) > 0,

(ll) K1+K3 <0

We will check whether one of these two conditions ( (i) and (ii) ) and til< 0 can be

satisfied simultaneously.

(i) Satisfying K3 +K1(1+K3)>0 and ul<0

We define Y as

20

 

Y =K3+ K1 (2.2.27)

 

(1+ K1)
01'
K3 =— K1 (2.2.28)
(1+ Kl)

Inequality (2.2.25), i.e., K3 + K1(1+K3) > 0 is true if and only if Y > 0.

Substituting (2.2.28) in (2.2.5), we get

W:[2K14Y+2K13Y(2+K4 +Y)+Y(2Y+K4(2+Y))+2Kl Y(3Y+K4(2+Y))+

K12(2Y0+3Y)+K4(1+41’+Y2)):'/(1+K1)2

(2.2.29)

We find from (2.2.29) that Y > 0 implies I/l>0 and (from Section 2.2.1) we know that

all'érll
I//>O implies -a—->0. It follows that the condition K3+Kl(1+K3)>0 results in
a) -

hardening of the dynamic stiffness with the increase in the forcing frequency. Thus

stability condition K3+Kl(1+K3)>0 and the softening condition V<0 cannot be

satisfied simultaneously.

(ii) Satisfying K1+ K3 < 0 and I// < 0

We will now check whether K] + K 3 < O and (1! < 0 can be simultaneously satisfied.

21

We define Z as

Z = —Kl — K3 (2.2.30)

01'

K3 =-Kl —Z (2.2.31)

The inequality K1+K3 < O is equivalent to Z > 0. Substituting (2.2.31) in (2.2.5), we

get

W = —1<12(1r<4 —2Z)+ 21(122 +Z[1<4 (—2 +2) +22] (2.2.32)

w < 0 implies

 

 

—K12+K4-\/(K12+K1K4)2+K§<Z<—K12+K4+\/(K12+K1K4)2+K}
2+2K1+K4 2+2K1+K4

The conditions (a < 0 and Z > 0 are satisfied if and only if

 

< —1<12+K4+\M<12+K1 [(4)2 +K}
2+2K1+K4

O<Z

 

OI’

 

K12‘K4-\/(K12'I'Kl[(4)2'I'K42
2+2K1+K4

 

0>K1+K3 >

01'

 

2 2
Klz-K4-\[(K12+K1K4) +K4
2+2K1+K4

 

“K1<K3 <—K1

22

In summary, we have found that the example system in Fig. 2.5 exhibit a stable
frequency-induced softening if and only if the following two conditions are
simultaneously satisfied:

(a) keff > 0 and

 

2 2 2
K12—K4—\/(K1+K1K4) +K4
2+2K1+K4

 

(b) -K1<K3<—K1.

2.2.4 Numerical example
Here, the example system (shown in Fig. 2.5 and again in Fig. 2.6) will be studied with

k k . .
K1 = k—l- = 0.3 and K 4 = f = 0.1 . With these values, A < 0 requires
2 2

k
—O.3616 < K3 < —o.3. We choose K3 = k—3 = —0.304 and also ii = 0.0002. The system
2 2

is loaded with a periodic force, F = F0 sin(a1).

F, X3
k4

1X1

k2 02
k 1 __I_x2

k3

 

 

ml'r

Figure 2.6. Spring-dashpot system with boundary conditions

23

 

 

 

 

 

 

 

 

_20 _ .r
O 2 4 6 8 10
time (t)
(a)
20 - 4
2?:
E: 0
xN
_20 - .r
0 2 4 6 8 10
time (t)
(b)

 

 

 

 

time (t)
(C)

Figure 2.7. x] (t) , x2 (I) and x3 (I) at the reference parameters and forcing frequency 1 Hz

24

 

 

 

 

 

 

 

 

 

20 - -
3o.
0
E; O - .
x"
-20
0 1 2 3 4 5
time (t) x 10-3
(3)
20 - -
AN
.5:
C
e. 0 - -
:01
0 1 2 3 4 5
time (t) x 10-3
0))
20 - ‘

 

 

 

time (t) x 10-3

(C)
Figure 2.8. x1 (I) , x2 (t) and x3 (I) at the reference parameters and forcing frequency

2000Hz

25

Figure 2.7 shows plots of (non-dimensionalized) x1 (t), x2(t) and x3 (I) , when the
forcing frequency is 1 Hz. Figure 2.8 shows the plot of x1 (I) , x2 (t) and x3 (I) , when the

forcing frequency is 2000 Hz. As can be seen in the figures, the x,- ’s are bounded,

showing stability of the system.

When the condition —0.3616< K3 <—O.3 is violated, the system is not simultaneously

stable and in the dynamic-stiffness softening regime. For example, when K3 in the

reference parameters is changed to -0.24, the system is not stable about the initial

conditions (i.e. x,- (0)=0 ). Figure 2.9 shows the plot of x3 (t) with K3 =-0.24 and
f =w/(27r)= 2000 Hz. As can be seen in Fig. 2.9, x3 (I) is departs from the initial

condition x3 (0) = 0 , signifying local instability of the system.

 

 

 

 

 

time (t) —3

Figure 2.9. x3 (I) with K 3 = —0.24 and frequency = 2000 Hz

26

When K3 in the reference parameters is changed to -0.39, the dynamic stiffiress increases

with the increase in the forcing frequency. Figures 2.10 and 2.11 show the plot of x3 (t)

when the forcing frequency is 1 Hz and 2000 Hz, respectively. As can be seen, the

amplitude of x3 (I) is smaller when the forcing frequency is higher, indicating hardening.

 

 

 

 

20 r -
:2»:
0
NM '
..20 - .
0 2 4 6 8 10

time (t)

Figure 2.10. x3 (I) with K3 = —0.39 and forcing frequency le

 

I
N
O

r

r

 

 

 

time (t) x 10-3

Figure 2.11. x3 (I) with K3 = —0.39 and forcing frequency 2000Hz

27

. k k k
At the reference parameters (vrz., K1 =-k—]—= 0.3 , K3 =—‘l=—3=—0.304

2 k2 k2

K4 2%:01 and%-=0.0002 ), the corresponding plots of the absolute value of the
2 2

effective stiffness (i.e. "le ), the real part of the effective stiffness (i.e. Re(k;{)) and the

imaginary part of the effective stiffness as functions of the forcing frequency are shown
in Figs. 2.12 (a), (b) and (c) respectively. In the stiffness plots, the stiffness will be scaled

(and made non—dimensional) by dividing it by k2. As can be seen in Fig. 2.12(a), the

static scaled dynamic stiffness is approximately 0.38. With increase in the forcing
frequency, the stiffiress monotonically decreases. The stiffness value is approximately
0.05 at the frequency 2000 Hz. It can also be seen in Fig. 2.12 (b) that the real part of the
effective stiffness is positive for the whole frequency range shown, i.e. in the range 0-

2000 Hz.

Now, the effect of changing £2— will be studied. When i2- is changed to 0.02, we
2 2

observe increase in the curvature of the dynamic stiffness vs. frequency curve, as in Fig.
c . . . . .
2.13. —2— thus seems to be an important parameter for tuning or achrevrng a desrred

2

dynamic stiffness vs. frequency curve.

28

 

 

 

 

 

0 500 1000 1500 2000
Frequency (Hz)

(a)

 

 

 

 

 

 

 

 

 

0 500 1000 1500 2000
Frequency (Hz) '
(b)
0.4 r
x” 0.3 i
.2: o 2 .
5 0.1 .
O 1 1 1
0 500 1000 1500 2000
Frequency (Hz)
(C)

Figure 2.12. Absolute value (a), real part (b) and imaginary part (c) of the effective

stiffness (k2, ).

29

 

 

0.1E

0 500 1000 1500 2000
Frequency (Hz)

 

 

Figure 2.13. Dynamic stiffness vs. the forcing frequency with c2 / k2 = 0.02

2.3. Estimation of Composite Material Parameters

The two-dimensional composite material in Fig. 2.1 was simplified into the one
dimensional model in Fig. 2.5 to analyze the material for frequency-induced softening
and stability. The parameters of the one-dimensional model (i.e. k1, k2, k3, k4 and 0;) will

now be used to construct a two-dimensional composite model. The parameters for the
. . . , :1: :1-
two-drmensronal composrte model are the complex moduli (E A = E A; and E B) of the

two phases A and B, the length fractions of phase B in horizontal and vertical directions

(i.e. x12 and as in Fig. 2.2 ) and (1:522 (equation (2.1.9)). There are 5 unknown

parameters (E A , E2, 22 , (D; and d) for the two-dimensional composite model and 3

equations, namely (2.1.4) (2 equations) and (2.1.5). Thus there are two free choices to be
made to satisfy the three equations relating the parameters of the 1-D model to the

parameters of the composite material.

30

From (2.1.4),

 

 

EAd = (2‘22) =k4(1—¢>2) (2.3.1)

or

(1-¢2)(1-12) zfl (2.3.2)
692 k4

where

o < 2.2 <1 (2.3.3)

and

0<¢2 <1 (2.3.4)

In equation (2.3.2), k, and k4 are given data. 112 and (D2 are chosen such that they satisfy

(2.3.2)—( 2.3.4).

From (2.1.4) and (2.1.5),

513;: k3(k2 How) (He)
EA (k3 +k2 'I' 1.020)) 142k]

 

(2.3.5)

If d is prescribed, E A is computed using (2.3.1) and E; is computed using (2.3.5).

Alternatively, E A may be prescribed and the corresponding d is computed using

31

d = k4(1-¢2) (2.3.6)
EA

2.3.1 Example phase properties
A one-dimensional example which shows frequency-induced softening as well as

. . . k k k
stability has the followrng parameters: K1 = —1 = 0.3 , K4 = —4 = 0.1, K3 = —3— = —0.304

k2 1‘2 k2

and _c_2_= 0.0002. Using equations in this section, the corresponding parameters for the
2

two-dimensional composite model (viz., E A , E}; , 22 , (D; and d) are found.

A square fundamental cell is assumed and thus 22 = Q. From (2.3.2)-(2.3.4),

(1—c>(1—¢2)_ k1_0-3_3

o2 k4 0.1
01'

22 = o; = 0.2087 .. 0.2

EL k3(k2+ic2w) (1—22)__4.053+i0.0008(0

 

 

_ 2.3.7
E A (k3 +k2 +ic2a2) 22k] 0.696+i0.0002w ( )
At a): 0,
11:
£3- = -5.824
EA
(0:0

and thus the Young’s modulus for material B is negative.

32

The elastic tensor for material A is given by

I VA 0
CA: VA 1 0 (2.3.8)
‘VA 0 0 (l—vA)/2

 

where v A is the prescribed Poisson’s ratio for phase A.

The elastic tensor for material B is given by

,3 1 VB 0
CB = B VB I 0 (2.3.9)
‘VB 0 o (l—vB)/2

 

where V B is the prescribed Poisson’s ratio for phase B.

The complex shear modulus of Material B is given by

t
*_ EB

Note that the shear modulus of material B is negative as the material is isotropic and the

Young’s modulus is negative. A negative shear modulus indicates the elastic moduli are

not strongly elliptic (Lakes and Drugan (2002)). If strong ellipticity is violated, the

material may exhibit an instability associated with the formation of bands of

heterogeneous deformation. However, the violation of strong ellipticity does not

guarantee the loss of stability of the inclusions: experiments show that the energy penalty

of band formation suppresses banding when particles of the material are sufficiently

small. Thus an instability criterion based purely on elasticity theory may not contain

33

enough of the physics of the actual material behavior, and hence may predict instabilities

in regimes where such do not occur in reality (Lakes and Drugan (2002)).

In essence, small-size inclusions of Material B may have potential to achieve the
frequency-induced softening; however, elasticity theories predict instability of the
composite material. While the physics suitable for the small-size inclusions of material B
in the matrix of material A needs to be studied further, inclusions of lumped lattices of
material B (made of springs and dampers) are used in the matrix of material A to achieve
frequency-induced softening in this work. Note that negative shear modulus can also
occur in the lumped elements; however the lumped elements cannot form bands and the
continuum conditions of ellipticity do not apply to the lumped elements (Lakes and

Drugan (2002)).

2.4 Stability of Negative Stiffness Phase B

The inclusions of phase B with negative shear modulus will be unstable and therefore
cannot be used directly. Instead, an anisotropic phase B is needed that has a negative
Young’s modulus, but also a positive shear modulus enough to maintain stability under

prescribed displacement at the phase boundary.

In order to avoid loss of stability, a model for an anisotropic phase B is proposed,

whereby B is itself a mixture of two isotropic constituents, Br and B2. In this model, a

constituent with negative stiffness (131) is always surrounded by a second phase of stable

34

material (B2). One such arrangement is shown in Fig. 2.14. The resulting effective
properties of phase B correspond to a “rank 1” layering of B 1 and B2 along the horizontal
direction (direction 1). The effective elastic tensor for phase B is computed using the

following standard layering formulas (given e. g. in Hassani and Hinton (1999)):

 

 

b11((0) (l-g)'11+gqi1 ( )
42(w)=[g11+(1—g)‘“—2)bn<w) (2.4.2)
r11 6111
2 2 2
= 1— — 12— 1— ‘11—? ———(b‘2(w)) 2.4.3
C2200) (gr22+( g)6122) [8,11% g)q11 + 141(4)) ( )
b33(w) = ‘133'33 (2.4.4)

(1- g)'33 + gqss

where qij and ry- are the ij-components of the elastic tensors (C B, and C 32) for phases B1

and B2, respectively. g is the volume fraction of B2 in B. The effective elastic tensor is

given by

b11((0) 512(4)) 0
CB= 1212(2)) 1222(0) 0 (2.4.5)
0 0 (233(0))

With suitable values for the free parameters g, CB, and C 32, the rank-l set up can result in
an elastic tensor for B where the 2-2 entry (C2222) is negative while 3-3 entry (C 1212) is

positive, i.e., B has positive shear modulus. The choice of B2 is somewhat arbitrary (as

35

long as it is a standard material). Both A and 82 are chosen to be made of the same

material. An idealization of B1 using a lumped parameter model is discussed next.

Stable 82 Unstable B1 Stable B2 Unstable B1

i .1 " ' ‘ .
oooo ' oooo
_ $4.: A 3?"; fr??? 537‘ v I ‘ Q 1 z- 7

Figure 2.14. A potential arrangement of constituents I31 and B2 to form material B

 

 

 

 

2.5. Implementation of Negative Stiffness

As discussed in Section 2.2, negative stiffness is important for the frequency-induced
softening of the dynamic stiffness of the system in Fig. 2.5. Negative stiffness can be
implemented using a bistable structure. Figure 2.15 shows a simple example of a bistable
structure. The hollow circles in the figure denote the hinges in the structure. There are
two linear springs in the structure, each hinged at both ends. The Y-shaped support
structure is assumed to be rigid. The springs are connected such that they form an arch-
shape as shown in the figure. h is the vertical distance between the two ends of either

spring. Similarly, L is the horizontal distance between the two ends of either spring.

 

      

L L

 

 

 

..L
F!"

Figure 2.15. A simple bistable structure

36

A vertical force P is applied at the common joint of the two springs, as shown in the
figure. y is the corresponding vertical displacement. The strain energy of the bistable

structure in Fig. 2.15 is given by,

 

‘I’(y) = 26w) (2.5.1)
where
A = \l(h + y)2 + L2 -\/E + L2 (2.5.2)

The (equilibrium) force corresponding to displacement y is given by

:53: 21.x d—A (2.5.3)

dy dy

P

The stiffness of the structure is given by

kn _—._=_ (2.5.4)

Figure 2.16 shows the typical plot of the applied force (P) versus the resulting
displacement (y). Figure 2.17 shows the corresponding plot of the stiffness (kn) versus
the displacement (y). As can be seen in Fig. 2.17, the bistable structure has negative
stiffiiess for some displacement values. The negative stiffiress corresponds to the negative

slope of the force-displacement diagram (Fig. 2.16).

37

 

Force
O

 

 

 

p......................-

 

Displacement

Figure 2.16. A typical force-displacement diagram for a bistable structure

 

Stiffness

 

 

 

 

Displacement

Figure 2.17. A typical stiffness-displacement diagram for a bistable structure

In the present example, the stiffiress is negative for

y€(y1,.vh)

where

 

 

(2.5.5)

(2.5.6)

(2.5.7)

38

The stiffness is most negative at y = —h and the corresponding value of the stiffness at

that point is

h2
k0=2k 1- -L—2-+1 (2.5.8)

Note that y=—h corresponds to the two springs aligned, as shown in Fig. 2.18. The

structure shown is the figure can be used to achieve negative stiffness. However, the

stiffness is not constant and also the displacement (y ) is to be strictly in the range
(y), y,,) in order to attain negative stiffiress. In contrast, in the examples in section 2.2, k3

is a negative constant and thus it does not depend on displacements. The desired negative

stiffness will be approximately achieved by appropriately constraining y.

 

Figure 2.18. Configuration with the maximum negative stiffness

Replacing the negative-stiffness spring (k3) in the original system (Fig. 2.5) with the

configuration in Fig. 2.18, we get a new system that can be represented schematically as

shown in Fig. 2.19.

39

 

Figure 2.19. Negative stiffness implementation

In Fig. 2.19, x2 = 0 corresponds to y = —h . Thus, x2 and y are related as

x2 = y+ h (2.5.9)

For very small oscillations around x2 = 0 (equivalently y = —h ), i.e. for 'le << lJ’h — J’Il

 

1/3
2
(or xle<<2 [ll—+1] —1), the stiffness kn ( in (2.5.4) ) will be approximately

L2

 

constant at k0 (given in (2.5.8) ) and the system in Fig. 2.19 will behave approximately

same as the linear system in Fig. 2.5.

2.6 A Lumped System Realization of Phase B;

One way to visualize the behavior and inner-workings of phase B; is to use a lumped

lattice of springs and dampers where some of the springs have negative stiffness.

40

Negative stiffness may be realized using e.g., bistable structures as shown in the previous
section. Figure 2.20 shows a potential lumped lattice configuration. As shown in the
figure, the lattice is a four-noded square, with two degrees of freedom at each node. The
nonlinear springs in the figure (shown as springs with an arrow across) correspond to
bistable structures that can be used to provide the desired negative stiffness. Material

phase B1 is obtained by tiling a plane with the two dimensional lattice shown in Fig. 2.20.

 

Figure 2.20. Two-dimensional lattice of phase B1

The 2D lumped lattice in Fig. 2.20 is basically made up of 6 standard linear solid (SLS)
elements interconnecting the 4 nodes. kil- and cy- in the figure denote spring stiffiress and

damper coefficient, respectively. The complex stiffness of each of the four identical SLS

elements forming the four sides of the lattice is given by

s = k31(k21 +110021)
k31+k21+icac21

 

(2.6.1)

41

Similarly, the complex stiffness of either of the two identical SLS elements forming the

two diagonals of the lattice is given by

= [(320622 +i0X’22) (2.62)

S
[(32 + [(22 + iwczz

The complex stiffness matrix (Kg?) for the 2D element in Fig. 2.20 is given by

 

 

Sl +S—2 fl- —S1 0 0 O -§-2— -£Z-
2 2 2 2
5% s1+—2 0 0 0 —s1 122— 3%
--sl 0 s1 +5; 4’1 31 52- 0 0
2 2 2 2
0 0 —"—2 s1 +2 52— —12- 0 —s1
K; = 2 2 2 2
o 0 43— 3—2 sl +S—2 ——S-2— —sl 0
2 2 2
0 —s1 32— 4’; —s—2 51 +—2 0 0
2 2 2
_S—2 —S_2 0 —SI 0 51+—2- f;
2 2 2
l- —-s—2 —£2— 0 —Sl O 0 :2- SI +£2—
2 2 2 2 _
(2.6.3)

The complex modulus tensor (CQ) corresponding to the lumped lattice in Fig. 2.20 is of

the form:

42

f11 1’12 0
0 0 1'33

Three numerical tests are conducted to find CQ:

(i) Tensile test 1: Degrees of freedom 1, 2, 3, 4, 5 and 7 are constrained (zero prescribed

displacement) and degrees of freedom 6 and 8 are given unit displacements. f“ is the

total reaction force along the degrees of freedom 6 and 8. fl] is given by

f” = 2S1+ 82 (2.65)

(ii) Tensile test 2: Degree of freedoms 1, 2, 3, 4, 5 and 7 are constrained (zero prescribed

displacement) and degrees of freedom 6 and 8 are given unit displacements. flz is the

total reaction force along the degrees of freedom 1 and 5. f12 is given by

f12 = 82 (2.6.6)

(iii) Shear test: Degree of freedoms 1, 2, 3, 4, 6 and 8 are constrained (zero prescribed

displacement) and degrees of freedom 5 and 7 are given unit displacements. f33 is the

total reaction force along the degrees of freedom 5 and 7. f33 is given by

43

fi3=fi

(2.6.7)

The aim is to construct an isotropic material using the lumped lattice in Fig. 2.20. In other

words, we want to express C Q in the form of (2.3.9), i.e.

 

Cm=l Vm 1 0
‘VBI 0 0 (1—v31)/2

This is possible if

S1 =82
* BS 83
E =._1=_2_
3‘ 3 3
and
v =—
81 3

(2.6.8)

(2.6.9)

(2.6.10)

(2.6.11)

As an example (to be used in section 4), the following parameters of the lumped lattice in

Fig. 2.20 are chosen:

k2] = 845EA
[C31 = —0.304k21

621 ”(21 = 0.0002

h2=b1

44

fiz=h1

Q2=Q1

Note that k31/k21 = —0.304 and 621 /k7_1 = 0.0002 here are kept the same as k3 / k2 and
c2 / k2 , respectively, in the example in section 2.3.1. The corresponding elastic tensor is

C31 in equation (2.6.8) with

 

E. _ _ 6.8501(1+i0.0002a))
3‘ 0.696 + i0.0002(0

and VBl=1/3-

The lumped system of the negative stiffness material B1 designed here may be mixed
with a typical positive stiffness elastic material B2 to build rank-l layered material B (as
described in section 2.4). For example, a purely elastic material may be used as phase B2.

The elastic tensor for phase B2 that will be used here is

 

C32 =1 2 V32 1 0 (2.6.12)
—VBZ 0 0 (1_VBZ)/2

where V32 = 0.3 is the prescribed Poisson’s ratio for phase B2 and E 32 is the Young’s

Modulus of phase B2.

45

If the volume fraction of B2 in B is 0.8, the elastic tensor for the rank-l layered material

B may be given by

541(0)) 142(0)) 0
CB: 512(0) b22((0) 0
0 0 b33(€0)

The absolute value, real and imaginary parts of the resulting b“(a)), b12((0) , b22 ((0) ,
and b33(a)) are plotted in Fig. 2.21. bU-(w) values in the plots are normalized by E A-

From Fig. 2.21 we find that only one component (b22) is negative and the shear modulus

of material B is positive, viz. 0.4954 E A- Material B thus can be made stable when
constrained. we also notice that in Fig. 2.21, “1222 ((0)" decreases with forcing frequency,

while "1911(0)" , ||b12(a))|| and "1233 ((0)" are approximately constant. Based on this

information, it may be possible to synthesize a lumped lattice for material B directly,
bypassing the construction of 81 followed by layering. In the next section, we will use

this information and construct such a lumped lattice.

46

 

 

 

0.01 a

1.42 ----------- I- » - . ~< 1.42

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

 

 

 

 

 

: oms . ......... .........
1.4 - 1.4 ' 0 i
0 1000 2000 0 1000 2000 O 1000 2000
Frequency (Hz) Frequency (Hz) Frequency (Hz)
(a) ”bu/EA” (b) Re(bll/EA) (C) Im(bll/EA)
x 10"3
0.44 f 0.44 r 3 T
// / 2 ......... . . .
043 .... 043, .......... .......... ..
' - l . ....................
0.42 2 0.42 4 0 '
0 1000 2000 0 1000 2000 0 1000 2000
Frequency (Hz) Frequency (Hz) Frequency (Hz)
((1) ”bIZ/EAH (e) Re(b12/EA) (f) Im(b12/E,4)
0.4
0.2 ........... .....
0.4 4 i O ;
0 1000 2000 0 1000 2000 0 1000 2000
Frequency (Hz) Frequency (Hz) Frequency (Hz)
(g) Hbzz/EAII (h) Re(b22/EA) (i) I111(1722/E4)
x 10—3
0.5 f 0.5 a 3
2 . ...............
0.495 0.495»-~-
1 . 1
0.49 i 0.49 ' 0 4
0 1000 2000 O 1000 2000 0 1000 2000
Frequency (Hz) Frequency (Hz) Frequency (Hz)
0) ”baa/EA” (k) Re(b33/EA) (1) 1m(b33/EA)

Figure 2.21. Components of the effective elastic tensor of material B

47

2.7 Phase B as a Lumped System

This section illustrates that it is possible to synthesize a lumped lattice model for phase B
directly (bypassing the construction of B1 followed by layering). The aim is to find a
phase B material which is approximately the same as that in section 2.6. Phase B in that
section is a rank-l layered material, whereas phase B here is a lumped system and is
stable enough to be used as an inclusion. Figure 2.22 shows a potential lumped lattice
configuration for phase B. As shown in the figure, the lattice is a four-noded square, with
two degrees of freedom at each node. ks and cs in the figure represent spring stiffnesses
and damper coefficients, respectively. The concentric circles at the four comers represent
rotational springs at the four corners. The nonlinear springs in the figure (shown as
springs with an arrow across) correspond to bistable structures that can be used to provide
the desired negative stiffness. Material phase B is obtained by tiling a plane with this two

dimensional lattice.

 

Figure 2.22. Two-dimensional lattice of phase B

48

The complex stiffness of each of the two identical standard linear solid elements forming

the left and right sides of the lattice is given by

51(0)) =

The complex stiffness matrix (K E) for the 2D element in Fig. 2.22 is given by

 

where

_ 2’0

R—
L2

k30(k20 “60020)
[(30 + [€20 + iaX'20

i k
k10+—%Q+R

k4_0
2

4‘10

a- >.- o
A 4:-
o Nlo

£40
2

31 +£go—+R

-31
k40

1:40

4'10

0 —R 0
-R 0 —s1
[(40 k40 k4_0
2 2 2
“in”; m _k4_0
2 2 2
ho Mo 1'40
fl +———+R -—~———
2 km 2 2
.50 fl Wm”.
2 2 2
0 —k1 0 0
‘51 0 -R

(27.1)

-540.

2

-540
2

(2.7.2)

1240

542+R

O *- 3':-
Na- Na-
lO [0

-S1

 

d

L is length of the side of the square lattice and r0 is the angular stiffness of the rotational

springs at the four comers (as shown in Fig. 2.22).

The complex modulus tensor (C B) corresponding to the lumped lattice in Fig. 4.2 is of

the form:

49

f11 f12 0
CB = f12 f22 0 (2-7-3)
0 0 f33

Four numerical tests are conducted to find C B :

(i) Tensile test 1: Degrees of freedom 1, 2, 4, 5, 6 and 8 are constrained (zero prescribed

displacement) and degrees of freedom 3 and 7 are given unit displacements. f” is the

total reaction force along the degrees of freedom 3 and 7. f“ is given by

f1] = ”‘10 + 1‘40 (2-7-4)

(ii) Tensile test 2: Degrees of freedom 1, 2, 3, 4, 5 and 7 are constrained (zero prescribed

displacement) and degrees of freedom 6 and 8 are given unit displacements. f22 is the

total reaction force along the degrees of freedom 6 and 8. f22 is given by

f22 = 231+k40 (2.15)

(iii) Tensile test 3: Degree of freedoms 1, 2, 3, 4, 5 and 7 are constrained (zero prescribed

displacement) and degrees of freedom 6 and 8 are given unit displacements. f12 is the

total reaction force along the degrees of freedom 1 and 4. f12 is given by

50

f12 = 1‘40 (2-7-6)

(iv) Shear test: Degree of freedoms 1, 2, 3, 4, 6 and 8 are constrained (zero prescribed

displacement) and degrees of freedom 5 and 7 are given unit displacements. f33 is the
total reaction force along the degrees of freedom 5 and 7. f33 is given by

f33 = 1‘40 +2R (2.7.7)

2.7.1 Stability of the lumped system

The stability of the inclusions of material B is tested here. Let us consider a 5x5 tiled
arrangement of the lumped lattice of Fig. 4.2 as an inclusion in the matrix of material A,
as shown in Fig. 2.23. The nodes are shown in the figure as small circles and a few of
them are numbered. The nodes of the lumped system at the phase interface (e.g. nodes 1,
2, 3, 4, 5, 6, 7, 12, etc.) are fixed to material A. We will test whether the lumped system
is stable as an inclusion. In particular we will determine whether displacements of the
interior nodes of the lumped system are bormded if motion of the nodes at the phase

interface is bounded.

Material A

“$112414“. it
* a 9 1o 41 15

7
MJifl“ I15 16 17]“;- MaterialB

 

 

1}

 

 

 

 

 

 

 

as

Figure 2.23 Inclusion of material B in matrix of material A

51

Let us consider the top left lattice in the lumped system separately. The four nodes
corresponding to this lattice are numbered 7, 8, 1 and 2. As indicated in Fig. 2.24, the
corresponding degrees of freedom are l and 2; 3 and 4; 5 and 6; and 7 and 8, respectively.

Nodes l, 2 and 7 are fixed, i.e., degrees of freedom 1, 2, 5, 6, 7 and 8 (Fig. 2.24) are fixed.

 

Figure 2.24. A lattice of material B with three nodes fixed

While the other three nodes fixed, node 8 is free to move along degrees of freedom 3 and
4. Assume that a small disturbance force acts on node 8. If x3 and x4 are the resulting
static displacements along degrees of freedom 3 and 4, respectively, then the

corresponding strain energy is

(I) =§xTKx (2.7.8)
where

52

 

 

 

k
K: k k k k (2.7.9)
_49_ 20 30 + 40 +R
_ 2 kZO +1630 2 -
x={x3} (2.7.10)
X4

Under the given boundary conditions, displacements x3 and x4 (at node 8) are bounded if
stiffness matrix K is positive definite. The positive definiteness of K insures that
displacements at the other interior nodes (such as 9, 10, 11, 14, 15 etc.) are also bounded.
With K positive definite, the lumped system is stable as an inclusion because under
infinitesimal displacements at the phase interface, displacements at the interior nodes are

unique and bounded.

2.7.2 Example of phase B
The following parameters may be used to obtain phase B approximately the same as that

used in section 2.6:

1:20 =1.672EA
[(30 = -0.304k20
C20 /k20 = 0.0002

k10 = 1:40 = 0.575,,

These values of km, km and 020 correspond to

53

0.5083(1+i0.0002(0) E

s a) =—
‘( ) 0.696+i0.0002(o

 

and therefore,

S1(0) = -0.73EA

These parameters suggest that, for stability,

R > 0.505EA

We select R = km = 0.57E A to maintain stability.

The corresponding elastic tensor is given by

 

21610 +1640 [(40 0

CB: [(40 251+k40 0
0 0 k40+2R

01'

"1.71 0.57 0 ‘
CB: 0.57 —1'0166(1+,'0'0002w)+0.57 0 EA (2.7.11)

O.696+10.0002a)
_ 0 0 1.71_

 

 

2.8 Realization of Phase B Using Tileable Bistable Structure

One way to realize the lumped systems for phase B1 (given in section 2.6) or B is to use
tileable bistable structure (described in chapter 5). Recall that material B1 is the negative

stiffness material and is one of the two component phases of the rank-1 layered material

54

B. As an example, Fig. 2.25 shows geometric realization of material B; as a periodic
structure. The rectangular solid structures shown in the figure are made of a Kelvin-Voigt
viscoelastic material. These structures connect 2D bistable structures with each other.
The fundamental cell of the periodic arrangement shown in the figure is the dashed
square region in the figure. This fundamental cell is equivalent to the 2D lattice shown in
Fig. 2.22. Chapter 5 will describe the tileable bistable structures that can be used as

building blocks of material B1 or B.

      
   
 

  

    
 
  

'«r

\v
9%”)!

““Wit‘www

       

Aux.g,4

 

Ah
\I

   

Figure 2.25. Material B as a periodic structure

55

CHAPTER 3
EXAMPLES ILLUSTRATING MATERIAL SOFT ENING AND IMPROVED

VIBRATION ISOLATION

The aim of the present chapter is to illustrate frequency-induced softening of a two-phase
composite and the corresponding improvement in the vibration isolation properties. Four
examples are presented to illustrate frequency-induced softening. The composite material
designed in the first example is illustrated to improve vibration isolation properties as
compared to the typical matrix material A. In particular, the transmissibility of an engine
mount made of the designed composite material is shown to be less than that of an engine
mount made of material A only. Before the examples are presented, the chapter
introduces to the method used to compute the average (or homogenized) viscoelastic
properties of a composite material. In the first three examples, the inclusion material B is
a rank 1 layered material (as introduced in the previous chapter) composed of alternating
layer of the example lumped material B1 in section 2.6 and the matrix phase A used in
that example. The three examples differ only in the matrix phase A used. A typical rubber
material (having a constant complex modulus) is used as a matrix material A in the first
example, whereas material A is purely elastic in example 2. In example 3, material A is a
Kelvin-Voigt viscoelastic material (i.e. imaginary part of the complex modulus increases
proportional to the forcing frequency). In example 4, the lumped system developed in
section 2.7 is used as phase B. In all the examples, the designed composite material

exhibits softening. While in the third example the material first softens and then hardens

56

with the frequency, in the other examples the softening is monotonous with respect to the

forcing frequency.

3.1 Homogenization of Viscoelastic Properties

With phases A and B in hand, the final step is to compute the effective properties of the
mixture of A and B, e. g., characterized by a simple arrangement such as that in Fig. 3.1.
This can be accomplished by numerical homogenization (notice that this suggests that the
mixture of constituents B1 and B2 takes place at a smaller scale than the mixture of A
and B). The representative cell is discretized using standard, 2D quadrilateral finite
elements and the effective properties of the mixture are obtained by exposing the cell to
three states of (unit) pre-strain, as is standard in numerical homogenization methods-
(details can be found in Yi et. al. (1998)). The computed homogenized elastic tensor of
the composite material has the information of whether the composite material softens
with fiequency. In particular, a monotonic decrease in the absolute value of the second
diagonal entry of the effective tensor indicates softening of the composite material in

direction 2.

 

Material A:
standard, elastic material

:_l//

 

 

 

 

Material B: /

viscoelastic material with
“negative stiffness”

 

Figure 3.1. Representative cell characterizing the (periodic) mixture of materials A and B

57

3.2 Example 1

Here the material softening and the improved vibration isolation properties are illustrated
with a numerical example. A typical rubber-like viscoelastic material is selected first for

phases A and B2. The (complex) elastic tensor of phase A is

 

CA=C91(1+l5A) With C2: ’42 VA 1 0
‘VA 0 0 (l—vA)/2

where Young’s modulus E A is arbitrary but real and positive and Poisson’s ratio
V A =0.45. The structural damping coefficient 6:; is 0.07. As indicated before, this

material is used also in phase B2 of B, i.e., CB2 = CA.

Phase B is a rank-1 layered material. Here phase B is constructed by alternating layers of

phase B2 and a negative-stiffness phase B1. The volume fraction of B2 in B is 0.8.

Phase B1 is made of a material such as the lumped lattice (Fig. 2.20). Its elastic tensor is

as in (2.6.8), i.e.

 

:1- 1 V81 0
C -— E31 v 1 0
’31 2 Bl
l—VBI

58

with V31 =1/3 and

. _ 6.8501(1+i0.000260)
3' 0.696+i0.0002a)

 

* .
The values of E31 and V31 used here correspond to the example material B1 computed

in section 2.6.

The final mixture corresponds to the periodic mixture of phases A and B characterized by
the periodic repetition of the cell in Fig. 3.1. The volume fraction of phase B is 16%. The
unit cell is discretized into 50x50 square plane stress elements for numerical analysis.
The resulting homogenized complex modulus tensor of the two-dimensional composite is

given by

* 011(0)) 012(0)) 0
CH: 012(4)) 022(0)) 0
0 0 633(0))

The absolute values, real and imaginary parts of cl-j (w) are plotted versus the forcing

frequency in Fig. 3.2. Ci]- ((0) values are normalized by E A .

59

 

 

 

 

 

 

 

 

 

 

2 a
l
0 _. 4 _.
0.1 l 10 100 1000
Frequency (Hz)
“CH/EA“
2 a
l
0 __n _.A ..4
0.1 l 10 100 1000
Frequency (Hz)
”Cu/E4“

 

 

 

 

0 ._.
0.1 l 10

 

100 1000
Frequency (Hz)

llczz/EAH

 

 

 

 

__A

O _. __. _.
0.1 l 10 100 1000
Frequency (Hz)

 

“033/54”

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 _t _.
l
O A ;A
0.1 l 10 100 1000
Frequency (Hz)
R€(Cll/EA)
2 a
l
0 _.
0.1 l 10 mo 1000
Frequency (Hz)
RC(C12/EA)
4
2 l
O ._4 __. #4 __4
0.1 l 10 100 1000
Frequency (Hz)
R6(sz/EA)
2
l .
O _. “_‘ M... M
0.1 l 10 100 1000
Frequency (Hz)
RC(C33/EA)

0.1 fl

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.05)
O - _. _. _x
0.1 l 10 100 1000
Frequency (Hz)
Im(C]1/EA)
0.1
0.05 1
j
0 A .4
0.1 l 10 100 1000
Frequency (Hz)
Im(cl2/EA)
4
2 \ 4
O ___4 4 _. ___.
0.1 l 10 100 1000
Frequency (Hz)
Im(sz/EA)
0.1
0.05
0 _. _. __.
0.1 l 10 100 1000
Frequency (Hz)
Im(C33/EA)

Figure 3.2. Components of the effective elastic tensor after mixing A and B

60

As can be seen in Fig. 3.2, “022(w)" decreases with forcing frequency, while ||c11(a))",
||c12(a))|| and "(233(0))“ are almost constant. This implies that frequency-induced

softening is achieved if the composite material is loaded in direction 2 while the material
is constrained in direction 1. Direction 1 is along the thickness of the layers (of phase B1
or B2) inside phase B, as shown in Fig. 3.3. Frequency-induced softening of the dynamic
modulus is achieved under unidirectional loading — when the composite material is
loaded in direction 2.

2 Material A:
I standard, elastic material

   

Material B:
rank-1 layered material

Figure 3.3. Phase B as layered material

3.2.1 Transmissibility analysis

Now the vibration isolation performance of the material designed in this example will be
demonstrated. A cylindrical block of the composite material synthesized is used here as
an engine mount and a transmissibility analysis is carried out. Figure 3.4 shows the
engine mount system, where l is the length of the engine mount. The area of cross-section
of the cylindrical engine mount is a (i.e. the radius of the cylinder is W ). The engine
mount supports an engine of mass m and the unbalanced disturbance force acting on the

mass is F. F s is the force transmitted to the automobile structure. The transmissibility of

61

the system is T = ||FS||/ ||F1|. An effective vibration isolation performance corresponds to a

low transmissibility.

TF

m

int-1:1

//////

in

 

 

 

 

 

 

 

Figure 3.4. Engine mount system

The components of the elastic tensor of the composite material used to build the mount
are those plotted in Fig. 3.2. In Fig. 3.4, the engine mount is unidirectionally loaded in
direction 2 and not constrained in the direction 1. As indicated earlier in this section,
direction 2 is perpendicular to the layering direction, as shown in Fig. 3.3. Under the

given boundary condition the strain (822) and the stress (0'22) in direction 2 are related as:

022 = Ec(w)€22

where

ECU”) = 62(0)) - V(0))012(0))
and

_ 02(0))
W60) _ 011(0))

EC(a)) and the effective stiffness of the engine mount kC (w) are related by

62

kC(€0) = ———EC(w)a

The absolute value, real and imaginary parts of EC((0) are plotted in Fig. 3.5, where

values are normalized by E A (i.e., the Young’s modulus of phase A). It can be seen in the

figure that Re(EC(0)) is approximately 3EA. In other words, the effective Young’s

modulus of the composite material is 3E4. Typically, the static stiffness (i.e., kC (0)) of

the engine mount is prescribed. In this case, suppose kC (0) is set to 170 N/mm. This

stiffness is achieved with, for example, E A = 4 N/mmz, I = 30 mm and a = 425 m2.

 

 

 

 

 

 

 

 

 

 

 

 

 

4 4 1 4
2 . 2 < 2 j
0 A 0 A 0 2
0.1 1 10 100 1000 0.1 1 10 1001000 0.1 1 10 1001000
Frequency (Hz) Frequency (Hz) Frequency (Hz)
"ECWVEA" Re(EC(0))/EA) Im<Ec(w)/EA)

Figure 3.5. The effective complex modulus EC(a))

 

JP
m

km): E rim)

la

Figure 3.6. Engine mount as spring-damper system

 

 

 

63

The engine mount in Fig. 3.4 may be represented as a system of a spring and a damper

connected in parallel, as shown in Fig. 3.6. The corresponding spring stiffness (k(a)))

and the structural damping coefficient (6(a)) ) are

16(0)) = Re(kc(w))

5W): Im(kc(w))
Re(kc(w))

k(a)) and 5(a)) are plotted in Fig. 3.7.

Note that the complex stiffness of the engine mount is kC(a)) , which is directly
pr0portional to the complex modulus EC(a)). As can be seen in Fig. 3.5, the absolute
value, real and imaginary parts of EC(a)) decrease monotonically with frequency. This
implies that the dynamic stiffness of the engine mount "kc (0))“, the real and imaginary

parts of the complex stiffness (i.e. k(a)) and 5(a))k(a)) , respectively) also decrease

monotonically with frequency. However, the effective structural damping coefficient

5(a)) increases with frequency as shown in Fig. 3.7. The dashed lines in Fig. 3.7
correspond to k(a)) and 5(0)) of the engine mount if the mount were made of phase A

alone.

64

 

 

 

 

 

 

 

 

 

 

1000 ~ - 4 w
100[ _ — _ _ _ — 7
E
z 10 r
..‘A
1 l 1
0.1““ 2 . 2-2-2. -
0.1 l 10 100 1000
Frequency (Hz)
103
10‘ .
co /
10_l t. ____________
10-3 - --- .. 2-2---. - -- . L _.
0.1 l 10 100 1000
Frequency (Hz)

Figure 3.7. Spring stiffness (k(a))) and damping coefficient (6(a)) ) of the engine mount

With k(a)) and 5(a)) in hand, the next step is to compute the transmissibility of the

engine mount, defined here as

= 11 1+ i6(a)) 1|
"1— ma)2 /k(a))+ 16(0))”

FS

F

 

T(a))=

 

 

 

 

Note that the above expression for T ( (0) may not be ideal for a nonlinear system where

k(a)) and 6(0)) are functions of frequency; however, this expression is used for

simplicity as is standard in engineering practice.

65

 

50 ~ - . 4

 

 

 

 

is?
v 0
.e
.5
.5.) _50.
.9
E
2 —100»
E
l—
_150 -A-A...i - - . ---“ . --i---1
0.1 1 10 100 1000

Frequency (Hz)

Figure 3.8. Transmissibility as a function of frequency

The transmissibility T(a)) is plotted in Fig. 3.8 (measured in decibels). The solid line
corresponds to the composite material, while the dashed line corresponds to phase A
alone. As can be seen in the figure, at high frequencies the transmissibility of the
composite material is significantly less than that of phase A alone. Thus the composite

material synthesized in this example performs better than a typical rubber material (phase

A).

3.3. Example 2

In this example, material A is purely elastic. The elastic tensor for material A is given by

 

1 VA 0
CA=1EA2 VA 1 0 (3.3.1)
"VA 0 0 (1—vA)/2

66

where V A = 0.3 is the prescribed Poisson’s ratio for phase A and E A is the Young’s

Modulus of material A.

As in example 1, material B here is rank-l layered material, whose constituent materials
are an elastic material A and a negative-stiffness material B1. Material B1 used here is the
same as that used in Example 1. The elastic tensor for the rank-1 layered material B is

same as that in section 2.6.

We consider the periodic mixture of materials A and B characterized by the periodic
repetition of the cell shown in Fig. 3.1. The volume fraction of the phase B is 16%. The
discretization of the unit cell is the same as in Example 1. The non-zero terms of the

resulting effective complex modulus are plotted in Fig. 3.9. As before, cy- ((0) values are
normalized by E A . As can be seen in Fig. 3.9, ||c22(a))|l decreases with forcing frequency,

while "012(w)“ increases with frequency. "011(09)" and "633(0)" are approximately

constant.

67

 

 

 

 

 

 

 

 

1_144 .......... .......... ..
1142 ........... - ........... 1
1.14 L
0 1000 2000
Frequency (Hz)
(a) ”Cu/EA”
0.4
0.35
0.3
0.25 1
0 1000 2000
Frequency (Hz)
((1) Ilclz/EAII
6

 

 

 

 

 

1000 2000

Frequency (Hz)

(g) ”022/154“

0.4

 

0.3995

 

 

 

 

0.399
O

1000
Frequency (Hz)

2000

(1) ”633/521”

 

 

1.144

1.142

 

 

1.14

 

1600
Frequency (Hz)

0 2000

(b) R6(C11/EA)

0.4

 

 

 

 

0.25 1
0 1000 2000
Frequency (Hz)
(6) R6(C12/E,4)
6

 

 

 

 

 

1000
Frequency (Hz)

0
0 2000

(h) Re(sz/EA)

0.4

 

0.3995 ' -

 

 

 

 

0.399
0

1000
Frequency (Hz)

2000

(k) Re(c33/EA)

 

 

 

 

 

 

 

 

1.5
1 ...................
05 . . .....................
O l
0 1000 2000
Frequency (Hz)
(C) 1m(Cu/E.4)
0
_001 .......................
_0.02 .....................
—0.03 i
0 1000 2000
Frequency (Hz)
(f) Im(cl2/E,4)
3

 

........ 1

  

 

 

 

1000

0 2000
Frequency (Hz)
(1) Im(sz/EA)

x 10"

 

 

 

 

 

1 000
Frequency (Hz)

0
0 2000

(1) Im(c33/EA)

Figure 3.9. Components of the effective elastic tensor mixing materials A and B in

Example 2

68

3.4. Example 3

The purpose of this example is to study the effects of introducing damping in material A.
Here the unit cell and its discretization are kept as that in Example 2. However material A
here is a Kelvin-Voigt viscoelastic material with damping coefficient 0.0002. The elastic

tensor of material A is

 

CA=1 42 VA 1 0 (3.4.1)
‘VA 0 0 (l—vA)/2

where E74 = EA(1 +i0.0002w) and VA = 0.3. EA is the Young’s Modulus ofmaterial A.

As in Example 2, material B is a rank-l layered material whose constituent materials are
materials A and B1. Material B1 here is same as that in Example 1. Unlike Example 2,
material B has Kelvin-Voigt viscoelastic material A as a constituent and therefore
material B here is different from that in Examples 1 and 2. The volume fraction of
material A in the rank-1 layered material B is again 0.8, the same value in Examples 1

and 2. Entries in elasticity tensor of material B are plotted in Fig. 3.10.

69

 

 

 

 

 

0 1000

 

2000
Frequency (Hz)
(a) llbll/EA”
l .5

 

 

 

 

 

 

 

 

0 .
0 1000 2000
Frequency (Hz)
(d) “bu/EA“
3 4
2 ...............
1 .
0 1000 2000
Frequency (Hz)
(8) llbzz/EA”

 

 

 

 

 

0 1000 2000

Frequency (Hz)

(.1) 11b33/EA“

 

1.5 ----------- I .......... 1

 

 

 

10 1000 2000
Frequency (Hz)
(b) Re(b11/E.4)
0.5

 

 

 

 

 

'3 0 1000

 

 

 

 

 

 

 

 

 

2000
Frequency (Hz)
(e) Re(b12/EA)
-0.5
—1
-l.5 -
0 1000 2000
Frequency (Hz)
0‘) Re(b22/EA)
0.5
0.35 L
0 1000 2000
Frequency (Hz)

00 Re(b33/E,1)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 > 1
0 i
0 1000 2000
Frequency (Hz)
(C) Im(bn/E,1)
1.5 4
l ....................
05 ) ....................
0 .
0 1000 2000
Frequency (Hz)
(1) Im(b12/EA)
3
2 .....................
1 ...................
0 i
O 1000 2000
Frequency (Hz)
(i) Im(b22/E,4)
1.5
l ...............
0.5
0 L
0 1000 2000
Frequency (Hz)
(1) In“(1233/1520

Figure 3.10. Components of the effective elastic tensor of material B in Example 3

70

In Fig. 3.10, we observe that while the real parts of bij(w) are of the same order of those
in section 2.6 (also used in example 2), the imaginary parts of by(w) are considerably

larger (since A is now a Kelvin-Voigt viscoelastic material).

Entries in the effective complex modulus tensor of the resulting composite material
(mixing A and B) are plotted versus the forcing frequency in Fig. 3.11. In Fig. 3.11, we

notice that cij (w) here are, in general, considerably different from that in examples 1 and
2. As expected, the imaginary parts of 011(0) are significantly higher than that in

example 1 and 2 as a result of frequency proportional damping in material A.

We also notice that "011(60)“, "012(60)“ and “C33 (0))" are monotonically increasing, while

||c22((0)|| first decreases with the increase of the forcing frequency until about 500 Hz,

after which it increases. This composite material may still be used as an engine mount
showing frequency-induced softening, since the frequency of vibration in automobile is
typically less than 500 Hz. Moreover, by reducing the structural damping coefficient of

material A, the maximum frequency showing the softening can be increased.

Also notice that "622 ((0)“ at 2000Hz is less than that at OHz. The drop in ”c22(a))|| in the

fiequency range 0-500 Hz is significantly more than the increase in "c22(a))" in the 500-

2000Hz frequency range. Thus the average dynamic modulus over the given frequency

range 0-2000Hz is still considerably less than the static modulus.

71

 

 

 

 

 

0 l 000 2 000

Frequency (Hz)
(a) ”Cu/EA”

 

 

 

 

 

0 4
0 1000 2000
Frequency (Hz)
((0 IICn/EAII
6

 

 

 

 

 

0 1000

2000
Frequency (Hz)
(g) 11c22/EA11
1 .5

 

 

 

 

 

0 1000 2000

Frequency (Hz)
0) “CB/EA”

 

1.2

115 .......... _ ............

 

 

 

1.1
2000

0 1000
Frequency (Hz)

(b) Re(Cl1/EA)

 

0.35

0.3

 

 

 

0.25

0 1000 2000

Frequency (Hz)
(e) Re(CIZ/EA)

 

 

 

 

 

 

0 1000 2000

Frequency (Hz)
0‘) Re(022/EA)

 

0.4

 

 

 

0.38

0 1000 2000

Frequency (Hz)
(1‘) Re(C33/E.4)

 

 

 

 

 

3 a
2 .................... 1
l ......................
0 .
0 1000 2000
Frequency (Hz)
(C) 1111(611/EA)

 

 

 

 

0 1000 2000

Frequency (Hz)
(1) 1111(012/EA)

 

 

 

 

 

 

 

 

 

00 7000 2000
Frequency (Hz)
(0 1m(022/1111)
1.5
1 ....................
()5 ......... ..
0 3
0 1000 2000
Frequency (Hz)
(1) 1m(033/521)

Figure 3.11. Components of the effective elastic tensor mixing materials A and B in

Example 3

72

3.5. Example 4

In this example we will test whether the lumped system for phase B given in section 2.7
leads to frequency-induced softening. We will use the simple topology of the
fundamental cell used in the previous three examples and shown again in Fig. 3.12. The

volume fraction of phase B in the fundamental cell is now 10.24%.

 

Material A:
standard, elastic material

)3 /

 

 

 

Material B: /

viscoelastic material with
“negative stiffness"

 

Figure 3.12. A simple square fundamental cell having square inclusion of phase B in the

matrix of phase A

Phase A is same as that used in example 1, i.e. the elastic tensor of phase A is

1 VA 0
E42 VA 1 0
’VA 0 0 (l—vA)/2

CA =C94(1+i6A) with C9, =

 

where Young’s modulus E A is arbitrary but real and positive and Poisson’s ratio

v A = 0.45. The structural damping coefficient (5,} is 0.07.

73

 

 

 

 

 

0 AA _‘ A
0.1 l 10 100 1000
Frequency (Hz)

(3)11011/5411

 

 

x

 

 

 

0
0.1 l 10 100 1000
Frequency (Hz)

((1)11012/5411

 

 

 

0 __z _A
0.1 l 10 100 1000

 

 

 

 

 

 

Frequency (Hz)
(8) ”622/134”
3 __ _.
2 D
1 l
0 g ..A __. _._m
0.1 l 10 100 1000
Frequency (Hz)
(1) ”033/54”

0 4.
0.1 l 10 100 1000

 

 

 

 

.__A

0 m _x _A
0.1 1 10 100 1000
Frequency (Hz)

 

(b) R6(C11/EA)

 

 

N

 

 

0 _. _.
0.1 l 10 100 1000
Frequency (Hz)

 

(e) Re(C12/EA)

 

 

 

0 _‘ _.
0.1 l 10 100 1000
Frequency (Hz)

 

(1'1) R6(sz/ E A)

 

 

 

 

4A #4._—A

 

Frequency (Hz)

(1‘) Re(C33/E.4)

0.2 —

 

0.1

 

 

 

._.4_‘

0 4
0.1 1 10 100 1000

 

 

 

 

 

Frequency (Hz)
(C) Im(Cii/EA)
0.2
0.1 l______,/\
0 n n
0.1 1 10 100 1000
Frequency (Hz)
(0 Im(c 12/EA)

 

.__/\

0 _.
0.1 1 10 100 1000
Frequency (Hz)

 

 

 

(1) 1111(022/EA)

0.2 ~ ~

 

0.1’ ‘

 

 

 

0 _._. _. __x _A
0. 1 1 10 100 1000
Frequency (Hz)

 

(l) Ifl'l(C33/EA)

Figure 3.13. Plot of the entries of (7;) for example 4

74

. * . . . . .
The average elastic tensor (CH) 1S computed usrng the numerrcal homogenization

method after discretizing the fundamental cell into 50x50 elements. (7;, is given by

.. 011(0)) 012(4)) 0
CH: 012(0)) 022(0)) 0
O 0 633(0))

The absolute values, real parts and imaginary parts of the entries of (7;) are plotted in Fig.

3.13. As shown in Fig. 3.13(g), ||c22 || decreases with frequency and thus we see that the

composite material softens in direction 2.

Transmissibility analysis is carried out for the material designed here using the procedure
followed in section 3.2.1. The corresponding transmissibility plot is given in Fig. 3.14.
The dashed line in the figure corresponds to material A, while the solid line corresponds
to the composite material. As can be seen in the figure, transmissibility for the composite

material is better than that for the typical rubber-like material A.

 

 

 

 

 

 

50 4 . a a
65
3 0
.Q‘
E, —50»
.9,
: -100t
E?
E_.
—150 - - - ..
0.1 1 10 100 1000

Frequency (Hz)

Figure 3.14. Transmissibility of the composite (solid line) as compared to phase A

(dashed line) for example 4

75

CHAPTER 4

TOPOLOGY OPTIMIZATION OF THE COMPOSITE

As seen in chapter 3, inclusions of the negative stiffness phase in the matrix of a typical
elastic/viscoelastic material phase lead to a composite material exhibiting frequency-
induced softening. The fundamental cell used in that section has very simple topology,
where the inclusion is square-shaped. This simple topology shows softening in only one
direction. When softening in the other entries of the elastic tensor is required, then it may
be difficult to easily find the corresponding topology. To achieve softening in any
prescribed direction, one option is to use topology optimization to design the composite.
In this chapter, the material distribution in the fundamental cell is optimized to obtain a
prescribed elastic tensor that exhibits softening. The objective fimction to be minimized
is the least-square error between the actual elastic tensor of the designed composite and a
prescribed elastic tensor. The stability of the composite is ensured by a constraint that
maintains positive definiteness. The design methodology starts with discretization of a
two-dimensional design domain into a number of finite elements. An element is made of
a mixture of a typical rubber material A and a negative stiffness material B; however, in
the final topology (i.e., the result of the optimization) it is desirable that all elements are
made purely of either material A or B. Material B used here is the lumped system built in
section 2.7 (and used in example 4 of chapter 3, i.e., section 3.5). The volume fraction of
material A in every element is treated as a design variable. Material B is anisotropic and
has negative stiffness in only one entry of its elastic tensor. So that the negative stiffness

of material B can be made available in any desired direction within an element, the angle

76

of orientation of material B in every element is also treated as a design variable. There
are upper and lower bounds on the total volume fraction of A in the unit cell. A gradient-
based optimization method (viz. the method of moving asymptotes) is used to solve the
problem. There are three problems solved and the resulting geometries of the base cell

are shown in sections 45.1-45.3.

4.] Background

The problem of material design by topology Optimization, which is also referred to as
‘inverse homogenization problem’, was introduced by Sigmund (1994 and 1995). The
aim is to design the microstructure of a two-phase composite material that has prescribed
effective elastic properties. The materials designed by the author are periodic, i.e., the
materials are made by end-to-end periodic repetition (tiling) of a base or fundamental cell.
Homogenization methods to find effective properties of a periodic material are well-
known, e.g., in Bensoussan et a1 (1978) and Guedes and Kikuchi (1990), where the
effective properties are determined solely by analyzing the fundamental cell. Sigmund
extends the periodic homogenization to an inverse problem where the homogenized
properties are given and the goal is to find the material distribution. The inverse
homogenization problem has been used to design materials with unusual or etxremal
properties such as materials with zero or negative Poisson’s ratio as in Sigmund (1995)
and materials having extremal bulk modulus as in Sigmund (2000), where an extremal

elastic property refers to the maximum or minimum value allowed by the theoretical

77

bounds (such as Hashin-Shtrikman and Cherkaev-Gibiansky bounds) on a composite’s

elastic properties.

Diaz and Benard (2003) extend Sigmund’s work by using a polygonal fundamental cell
and produced a material distribution that may not be obtained using a simple square or
rectangular fundamental cell. The objective function to be minimized is the least square
error between the prescribed elastic modulus and the achieved elastic modulus. The total
weight of the fundamental cell is constrained to be equal to a prescribed value. This
formulation is slightly better in computational performance than that in Sigmund (1995)
where the objective is to minimize the weight of the composite material and the effective
elastic tensor is constrained to be equal to a prescribed value. Diaz and Benard (2003)
also introduce thickness of the fundamental cell as a design variable which makes it

possible to scale all entries of the elastic tensor without changing the topology.

The methodology to design an elastic material can be extended to design a viscoelastic
material as in Yi et al. (2000), where the homogenization and inverse homogenization in
frequency domain are formulated applying the correspondence principle. According to
the correspondence principle, the viscoelastic homogenization (or inverse
homogenization) process becomes identical to that of the elastic case except that the

variables are complex in the viscoelastic case.

In this work, Diaz and Benard (2003) will be extended to design a viscoelastic composite

that has the prescribed viscoelastic properties.

78

4.2 The Optimization Problem

As in Diaz and Benard (2003), the objective function to be minimized is the least square
error between the prescribed elastic modulus and the achieved elastic modulus. Since the
composite material being designed here is viscoelastic, the elastic tensor has real and
imaginary parts that are functions of the forcing frequency. The objective function is

given in (4.2.1)

N
¢=izwl:%(t;E-;E*)T jw(sz-jE*)] (4.2-1)

i=1 j=l

where t is a scaling parameter, which physically may be interpreted as the thickness of

the cell. 1’ El: and j-E , respectively, are the prescribed and the actual elastic tensors in a
i 1 * . , T
6X1 VCCIOI' form. jE and jE are in the form Of {Cll11,62222,01212,C”22,C“12,62212} .

In J'E* and ;E , the superscript i = 1, 2 and 3, represent real part, imaginary part and

absolute value of the elastic tensor. The absolute value of the elastic tensor is also
included in the objective function as for some applications the absolute value may also be
of interest. The subscript j denotes the jth forcing frequency. The elastic tensors are
prescribed at a prescribed number (No) of forcing frequencies ((01,602 , (qua, ). The
objective function is a sum of 18Na, individual least square error terms and has provisions

to prescribe relative importance or ‘weight’ to these terms. The diagonal weight matrix

79

j-W sets such relative importance to the entries of the elastic tensor. i-W , ZJ-W and 3J-W ,

respectively, are the prescribed weight matrices for the real part, imaginary part and

absolute values of the elastic tensor. For example, if only the absolute value of Cu” is of

interest, then )W and ij are zero matrix. All entries of 3jW are zero except the first

 

 

 
   
 

 

 

 

diagonal entry.
”7 = 1 P4 = 0
100% Material A 100% Material B
\ I
/ 1
p18 = 0.5 1 p20 = 0.2
50% Material A \ 20% Material A
50% Material B 80% Material B

 

 

Figure 4.1 Element densities as design variables

The design domain for a fundamental cell is discretized into a number of finite elements,
as shown in Fig. 4.1. The ‘density’ (pg) of each element of the fundamental cell is the
design variable. p, is the volume fraction of phase A in element e and (1- p.) is the
volume fraction of phase B in that element, as shown in the figure. pe thus can take any

value between 0 and 1, i.e.,

OSpeSI fore=1,2,...,N

where N is the total number of elements in the discretized fundamental cell.

80

The phase B used in this chapter is the same as that built in section 2.7 and used in the

example in section 3.5. This phase is anisotropic and its elastic tensor, as in (2.7.11), is

 

"1.71 0.57 0
. 661+'0.0002
CB: 0.57 101 ( _’ ”+0.57 0 EA (4.2.2)
0.696+10.0002a)
_ 0 0 1.714

 

 

As can be seen in (4.2.2), only the C2222 term in the elastic tensor of material B has
negative stiffness, which is crucial for the frequency-induced softening. To achieve
softening in the terms other than c2222, material B needs to be rotated and suitably
oriented. To achieve softening in more than one term simultaneously, the base cell may

need to have material B oriented in more than one direction. Material B will, therefore, be

allowed to be rotated from the prescribed original orientation. If Be is the angle of rotation

of material B in element e, then the rotated elastic tensor is given by

Cijkl = aipajqakralscpqrs (4-2-3)

where a)!- are the entries of the transformation matrix given by

[cos Be — sin 68]
a = . (4.2.4)
srn (98 cos He

where Be in this work is bounded as follows:

—7r/ZSBeSfl/2 fore=l,2,...,N

81

61:1'1'14 65:0

a7=m2

 

Figure 4.2 Angle of orientation of material B in an element as a design variable

There is a constraint used on the volume fraction of phase A in the fundamental cell. The

constraint on the volume fraction is given by

where A = {A 1, A2, AN} T is a vector containing element areas, v“ is a prescribed upper

. . . I . .
bound on the volume fract1on of phase A in the composrte, v rs the corresponding

prescribed lower bound and

N
42:24;

e=l

is the total area of the cell.

82

Stability of the composite material is ensured by implementing a constraint that basically
requires the real part of the elastic tensor to be positive definite. In particular, eigenvalues
of the real part of the effective elastic tensor (in 3x3 matrix form) are computed and are

desired to have positive real part. The stability constraint may be written as

—Rezl£fl<0

where ,6 is a prescribed negative real number, which is the lower bound for the real part

of xi. Here/1 is an eigenvalue of the real part of jEH and has the real part the lowest

among all the eigenvalues. A is obtained by solving

det[Re jEH 413,3] = 0 (4.2.5)

jEH is the average elastic tensor for the composite in a 3x3 matrix form:

01111 01122 01112

H
jE = 01122 02222 62212
01112 02212 61212

and is given by

H l l T
qu,. 7453183 -—e‘1 ) E(p,0,a)j)(£6 -e’)dY (4.2.6)

83

In (4.2.6), Y represents the fundamental cell i.e. the periodic domain over which the

inverse homogenization problem is defined. |Y| = A}; is the total area of the fundamental

cell. E(p,0, (0) is the element elastic tensor given by

Ee(pe,6e,w)=peEA1w>+(1—p.)EB(6e,w) (4.2.7)

EA((0) and EB (1990)) are the elastic tensors for material phases A and B respectively.
The element elastic tensor Ee as a function of the phase elastic tensors (EA and EB )
and element density (pg ) such as in (4.2.7) is typical in a material distribution problem

(as can be seen in Bendsee and Sigmund (2003)). E8 is typically constructed such that

pe is discouraged to take an intermediate value between 0 and 1 while improving the

objective function or satisfying the constraints. Here Ee(pe,6le,a)) may be physically

interpreted as the average elastic tensor of a Voigt composite made of phases A and B

having volume fractions pe and (l — pe) respectively.

In (4.2.6),

63) = {1,0,0}T

£3 = {0,1,0}T
and
£3 = {0,0,1}T

84

are the unit test strains applied for the purpose of determining the average elastic tensor
of the composite material. sq in (4.2.6) is the strain induced by the test strain 88 and is

solution of this problem (which is solved by finite element approximation)

)Y8(v)TE8(uq)dY =1Y£(V)TE£ng’ for all Y-periodic functions v (4.2.8)

where 21" is a Y-periodic solution. If Y is a unit square, then a function (p(yl,y2) is Y-

periodic if ¢(yl + m,y2 + n) = (0(y1,yz) for any (y1,y2) 6 Y and integers m and n.

In summary, the optimization problem may be written as follows:

Findp= {p1,p2, ...,pN}T,0= (6,, 192, 610}de11001
3 ”6121 T. . .
minimize ¢(p,0, I): :2); (tle E) ;W(tJ'-E—J'-E*))
-lj-l
subject to 03,06 SI fore: 1,2, ...,N (4.2.9)
427236851!” fore=l,2,...,N

—Re/1$fl<0

v l<_erpe (Va

A2 8:]

Since (D is an unconstrained, convex function of t for fixed p and 0, following Diaz and
Benard (2003), it is possible to express the optimal value oft as the solution of 001/01 =

which is

85

‘2
where
3 N . . .
cl=zf(}E*)T )wle (4.2.11)
i=1j=l
3 N . . .
c2 =Zf(}E)T )ij (4.2.12)

Incorporating t in the original objective function, the objective functron to be
minimized may be rewritten as
3 Nwl

¢(p,0)= 22B (”E—j E*)Tj-w(z*jE-J'IE’)) (4.2.13)

1- l j: l
The optimization problem may now be written as

Findp={p1,p2,...,pN}Tand0={61, 192,..- 610}Tthat
,. T. . .
minimize ¢(p, 0): 22):; (t iE-‘E j E) }W(t*;E—j'-E*))
1— 1}"- 1
subject to O Spe SI for 6 =1, 2, N (4.2.14)

—7r/2.<_BeS/T/2 fore=l,2,...,N

— Re x1 _<. ,6 < 0 (stability constraint)
l< — 12 Age 6 S v" (voltune constraint)
A2 e=l '

86

4.3 Sensitivity Analysis

The optimization problem (4.2. 14) is solved using a gradient-based optimization method,
viz. the method of moving asyrnptotes ( Svanberg (1987) ). The gradients of the objective
function and the constraints with respect to the design variables (required for the

numerical solution of the optimization problem) are derived in this section.

(i) Objective Function
The gradient of the objective function with respect to a design variable xe (for example,

pe or He) is given by

 

 

 

3 N s i
d¢ 22(3) *i 1 * 1 d, 1 *djE
—= z -E E w — E (4.3.1)
dxe I-lj-l ( J J ) J dxe J dxe
where
* 3 Na) . . T . d 1:E
d’ = 2.1—( 1'1:qu jE) )w 1 (4.3.2)

The derivatives of the real part, the imaginary part and the absolute value of the elastic

moduli, i.e.,

 

 

 

dle djE
= Re—, (4.3.3)
dxe dxe
013.13: d jE
= Irn (4.3.4)
dxe dxe
and

87

 

 

l 2

0113-15,, 1 1 d 1.12,, 2 d 1.15,,

= 3 J. k + jEk— (4.3.5)
00:, )E1 00, dx,

are easily obtained from

H H
dqu, zaJ‘Eqr = 1 1

—j —( g—e‘1)T(EA(wj)-EB(6,,wJ-))(55—e’)dA

 

 

 

dpe ape A0 82
(4.3.6)
and
d4? = 33:5]; = (1:0)) 1.32185] _gq)T 0133330)) (.95 —gr)dA
(4.3.7)

(ii) Stability Constraint
The derivative of the real part of the eigenvalue with respect to a design variable xe is

the same as the real part of the derivative of the eigenvalue, i.e.,

 

d R
( “Linea (4.3.8)
dxe dxe

where the derivative of the eigenvalue is

fl.=[.d_1_b;22___d”’5 “91%) / (3,12 —2IE/t+115) (4.3-9)

dxe dxe dxe xe

I E , II E and III E are the three invariants of Re j E” and are given as

[E = tr(ReJ-EH) (4.3.10)

88

IIE =%[(0(Re jEH ))2 —tr (Re 1.1:” )2) (4.3.11)

1115 = det(Re jEH ) (4.3.12)

H
d 1.15,,

dxe

 

The gradients of the three invariants are easily computed as is already known

from (4.3.6) and (4.3.7).

(iii) Volume Constraint

The gradient of the total volume fraction of phase A is trivially given by

 

 

d 1 N A

— Ap =4 (4.3.13)
61/913142; e e] A:
and
d 1 N

— A =0 4.3.14
44142.24 .0.) < >
4.4. Solution Strategy

To solve the optimization problem (4.2.14), an efficient strategy (in terms of fast
convergence of the solutions) is found by trial and error. The strategy consists of the

following three steps:

Step 1: Propose initial angle of rotation (00) of material B.

89

There is no hard and fast rule for the initial assignment of 66; however, for faster and
better convergence He should be such that negative stiffness is available in all the
directions along which softening is required. He = 0 and 1t/2 correspond to negative
stiffness available in vertical (i.e. direction 2) and horizontal (i.e. direction 1) directions,
respectively. If the shear modulus is desired to soften, then 6e 2 i'Tt/4 may be used as

this angle makes negative stiffness available in the diagonal direction of the element and
softening in the diagonal direction leads to softening in shear. If softening is desirable in

more than one entry of the elastic tensor, then Be =i1r/ 4 is, in general, found to give

good results.

Step 2: Obtain p*, which is a solution of the following density optimization problem:

Find p = {,0}, p2, pN}Tthat

minimize ¢(D 90): 2:432:0[é' (* JiE JI'EI')T ;w(t*lE_ JiE*):|

i-lj-l

subjectto OSpeSl fore=l,2,...,N (4.4.1)
— Re xi 5 ,6 < 0 (stability constraint)
v1 S —— Z Aepe S v" (volume constraint)
2 3:]

To facilitate a “black and white” solution, a morphology-based filter as given in Sigmund
(2007) has been used. The morphology-based filter is an extension of the original density

filter introduced by Bruns and Tortorelli (2001). An “open operator”, in particular, is

90

used in the present work. As described in Sigmund (2007)), an open operator
corresponds to erosion followed by dilation. In image processing, this operator
corresponds to removing dimensional details smaller than a prescribed filter size. A
Heaviside filter (basic version proposed in Guest et a1. (2004)) is used as a dilation
operator, whereas a modified Heaviside filter (proposed in Sigmund (2007)) is used as an
erosion operator. It was observed that the open filter used in this work tremendously
decreased (with iterations) the total volume fraction of the material corresponding to

pe =1 (i.e. material A). To stabilize the volume fraction of material A, the volume

constraint is found to be important. It was also observed that the optimization process

works better and faster with the following constraint:
1 N

—Zpe(1—pe) S6

N e=l

where 5 is a prescribed number. This constraint has been used in solving the density
problem (4.4.1). Note that while solving the density optimization problem, a parameter of
the Heaviside function (viz. ,6 in Sigmund (2007)) is doubled every 20th iteration as a
scheme to gradually reduce number of gray elements in the design. The change in this
parameter leads to a jump in the objective function value, as can be seen in the iteration

history plots presented for the solved examples.

Step 3: Obtain 0*, which is a solution of the following angle optimization problem:

91

Find0= {61, 92, 6N}Tthat
3 Na) T
2, 1 *. . * . *. . *
minimize ¢(p ,0): l:— t I-E— l-E '-W t 1-E— 1E )
g; 21 J J ) J 1 J J )

subject to —Iz'/2 S 6,, S7r/2 for 8 =1, 2, N (4.4.2)

— Re 21 S ,3 < 0 (stability constraint)

4.5 Examples

In all the three examples given in this section, material phases A and B are fixed. The
prescribed lower and upper bounds on the total volume fraction of A in the composite are
0.7 and 0.9 respectively. A and B are as used in the example in section 3.5. The elastic

tensors for A and B are as follows:

 

0.45 0

(1+i0.07)
A=____2- 0.45 1 0

1’0'45 0 0 (1—0.45)/2

"1.71 0.57 0 1
CB= 0.57 _1.0166(1+20.000200)+057 0

0.696+i0.00020)
_ 0 0 1.71_

 

 

4.5.1 Optimization example 1
In this example, the goal is to design a composite that will exhibit softening in shear. The

prescribed ”01212” are 0.73, 0.51 and 0.31 for the forcing frequencies 0, 100 and 2000 Hz

respectively.

92

The base cell is square-shaped and is discretized into 30x30 elements. The initial guess
for the angle of rotation of material B (in Step 1) is 90 degrees for all elements. The
initial guess is arbitrary. After solving the density Optimization problem in Step 2, an
optimal geometry is obtained, which is shown in Fig. 4.3. The solution was assumed to
be symmetric about both the horizontal and vertical centerlines of the base cell. In all the
examples in this work, only a quarter (top-right) of the base cell is designed and the other
three quarters are created as mirror images of the designed quarter. The dark areas in the
designed base cell represent material A and the white areas represent B. A few grey
elements in the cell represent a mixture of materials A and B. As in a typical material
design problem, it is very difficult to completely eliminate the gray areas because gray
elements are mathematically allowed in the optimal designs. The total volume fraction of
material A in the composite is 0.764. The obtained ||c1212|| is 0.4167, 0.4133 and 0.4011,
respectively, at the three frequencies. The target and the achieved ||c1212|| are plotted in
Fig. 4.4, where the target values are shown as crosses while the achieved values are
shown as circles. The iteration history (i.e. the objective fimction plotted versus iteration
number) is shown in Fig. 4.5. The objective function value is almost constant throughout
the iteration history as can be seen in the figure. This is as expected. The negative
stiffness is currently not oriented along the element diagonals, whereas softening of shear

modulus requires negative stiffness to be oriented along the diagonals.

93

 

Figure 4.3. The composite designed in example 1 (initial angle 90 degrees)

 

 

 

 

l .5
X target
0 achieved
1
r:
0.50 6 Q?
0 I .
0 10 l 00 2000

Frequency (HZ)
Figure 4.4. Plot of the target and achieved ||c1212|| for the density problem of example 1

(initial angle 90 degrees)

 

Objective function
3

 

 

 

O 20 40 60 80 100
Iteration no.

Figure 4.5. Iteration history for example 1 for the density problem of example 1 (initial

angle 90 degrees)

94

The angle optimization problem is Step 3 is now solved in order to improve the objective
function value. The resulting angles of rotation for material B in the elements are shown
in Fig. 4.6, where every element is shown as a square and the line across the square
represents the actual orientation of material B. For example, a horizontal line denotes 0
degree and a vertical line represents 90 or -90 degrees. The optimal angles of rotation
range between 45 and 90 degrees in the top-right and bottom-left quarters of the base cell
(i.e. between -45 and -90 degrees in the other two quarters). These angles approximately
orient the negative stiffness along the diagonals of the base cell and help to achieve
softening in shear. The solution is as expected. The resulting dynamic modulus (”c1212”)
is plotted and compared with the target value in Fig. 4.7. As can be seen in the figure, the
dynamic modulus is the same as the target value and thus decreases with frequency.
Figure 4.8 shows the iteration history. The objective function has, in general, improved
with the number of iterations. The humps in the plot may be avoided by using a smaller

step size in the optimization routine.

 

Figure 4.6. Orientation of material B for example 1 (initial angle 90 degrees)

95

1.5 - ~

 

 

 

 

 

 

 

 

x target
0 achieved
1 . .
<19
0.5 - ® -
GD
0 . .2;
0 10 100 2000
Frequency (Hz)

Figure 4.7. Plot of the target and achieved ”01212“ for the angle problem of example 1

(initial angle 90 degrees)

20

 

..................................

.................................

Objective function

 

 

20 40 60 80 100
Iteration no.

 

 

Figure 4.8. Iteration history for the density problem of example 1 (initial angle 90

degrees)

_E3c_qmple 1 with g_different initial angle
The problem is now solved with a different initial angle of rotation for B. The angle is as
shown in Fig. 4.9. The angle is 45 degrees for the upper-right and lower-left quarters of
the base cell and -45 degrees elsewhere. The base cell obtained by solving the density
problem (Step 2) is shown in Fig. 4.10. The volume fraction of material A in the cell is
0.794. The elastic modulus is plotted in Fig. 4.11. As can be seen in the figure, the

achieved modulus is already the same as the target value and therefore angle optimization

96

is not needed. The optimization problem was, thus, easier to solve with 45 degrees as
initial guess. 45 degrees is a good guess as this angle will orient the negative stiffness
along the diagonals and help shear modulus to soften. This may imply that with a good
guess for the initial angle, the target properties may be achieved by solving the density

problem only.

 

Figure 4.10. Optimal base cell for example 1 (initial angle 45 degrees)

97

 

1.5 - a

 

 

 

 

 

 

 

x target
0 achieved
1 . .
69
0.5 . ® «
, - - T
0 10 100 2000
Frequency (Hz)

Figure 4.1 1. Plot of the target and achieved ||c1212|| for example 1 (initial angle 45

degrees)

 

Objective function

 

 

 

 

 

 

 

 

 

O 20 40 60 80 100
Iteration no.

Figure 4.12. Iteration history for example 1 (initial angle 45 degrees)

The iteration history is shown in Fig. 4.12. As can be seen in the figure, the objective
function undergoes a jump a few times. The jumps in the objective function value at
iterations 20, 40, 60 and 80 are because of the doubling of a parameter of the Heaviside
function (viz. ,6 in Sigmund (2007)) as mentioned in section 4.4 (Step 2). The other jumps
may be because the effective elastic tensor is very sensitive to the design variables (p) to

the extent that a small change in the design variables tremendously changes the objective

98

function value. This jump in the objective function value may be avoided by using a

smaller step size in the optimization routine.

4.5.2 Optimization example 2

This example designs a material that softens in two directions — in the cum and C2222
entries of the elastic tensor. The desired “0111111 and ||c2222|| are equal and given to be 3.07,
2.59 and 0.99 at three prescribed forcing frequencies 0, 100 and 2000 Hz respectively.

The initial guess for the angle of material B is 90 degrees for all elements.

The optimal base cell obtained by the density optimization (Step 2) is shown in Fig. 4.13.
The volume fraction of material A is 0.774. ||c1m|| and IIC2222|| for the three frequencies
{0,100,2000} Hz are {3.2368, 1.9447, 0.8793} and {1.4750, 1.4728, 1.4728}
respectively. The dynamic moduli are plotted in Fig. 4.14. As can be seen in the figure,
”01111” is close to the target value (and therefore softens); however, ||c2222|| is considerably
different from the target and is almost constant with respect to the frequency. 90 degree
angle of rotation for material B corresponds to negative stiffness oriented along
horizontal direction (i.e. direction 1), which helps the composite to achieve softening in
the en” entry. 11C111111 is constant here because it is not possible to soften cu” by

orienting negative stiffness only along direction 2. Figure 4.15 shows the iteration history.

99

 

Figure 4.13. Geometry of the base cell designed in example 2 (initial angle 90 degrees)

4 x target 4 x target
9) O achieved 0 achieved
3

 

 

 

 

 

 

 

 

- * 3):
X X
2 2
O O O
1 as 1 u
0 0
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(81) 110111111 (b)|102222||

Figure 4.14. Plot of the target and achieved dynamic moduli for the density problem of

example 2 (initial angle 90 degrees)

U!
Q

 

NW4}
COO

 

Objective function

.—
O

 

 

O

 

40 60 80 100
Iteration no.

0
0

Figure 4.15. Iteration history for the density problem of example 2 (initial angle 90

degrees)

100

Angle optimization (Step 3) is now carried out. The resulting angle is shown in Fig. 4.16.
It may be seen in the figure that the angles corresponding to the white areas of the base
cell (shown in Fig. 4.13) are either 90 degrees or about 45 degrees. While a 90 degree
angle helps the composite to soften in direction 1, a 45 degree angle helps softening in
the both directions 1 and 2. The resulting ||c1m|| is 3.1098, 2.4162 and 0.9631,
respectively, for the three forcing frequencies. The corresponding ”C2222“ values are
2.9720, 1.4681 and 1.2625. The target and achieved 110111111 are plotted versus the forcing
frequency in Fig. 4.17. As can be seen in the figure, the absolute values of cm) and 02222
decrease with frequency and thus the material designed here shows softening in both —
horizontal and vertical directions. However, the achieved “C2222” is a little off the target
value. With a different initial angle, the results could be improved as illustrated next. The
iteration history is shown in Fig. 4.18, which indicates gradual improvement in the

objective function value with iterations.

 

Figure 4.16. Optimal angle of orientation for material B for example 2 (initial angle 90

degrees)

101

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 - x target ' 4 L X target
0 achieved 0 achieved
3‘59 1 365 .
O x
2 ~ . 2 '
O
1 » 00 1 1 i?
0 A .0 O A ..
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a) 110111111 (1)) ”02222“

Figure 4.17. Plot of the target and achieved dynamic moduli for the angle problem of

example 2 (initial angle 90 degrees)

 

40 T T 3

30 ....... g ...... ...... g ...... f. ..... .

Objective function

 

 

 

 

0 20 40 60 80 100
Iteration no.

Figure 4.18. Iteration history for the angle problem of example 2 (initial angle 90

degrees)

Example 2 with a different initial angle

The initial angle is now changed to i45 degrees as shown in Fig. 4.9 in the previous
example. The solution of the density problem is shown in Fig. 4.19. The volume fraction
of A is 0.7 (i.e. equal to the lower bound). The dynamic moduli are plotted in Fig. 4.20.

The iteration history is given in Fig. 4.21.

102

 

Figure 4.19. Geometry of the base cell designed in example 2 (initial angle 45 degrees)

 

 

 

 

 

 

 

 

4 X target 4 x target '
O achieved 0 achieved
3:: 3x .
x x
2 2 -
O 0 o O O o
l 3< 1 )1
0 ‘ 0
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a) “0111111 (b) “(3222211

Figure 4.20. Plot of the target and achieved dynamic moduli for the density problem of

example 2 (initial angle 45 degrees)

 

 

Objective function

 

 

 

0 20 40 60 80 100
Iteration no.

Figure 4.21. Iteration history for the density problem of example 2 (initial angle 45

degrees)

103

1|C1111|| is 1.6887, 1.5127 and 1.3313, respectively, for the three fi'equencies. The
corresponding ||C2222|| is 1.5455, 1.4526 and 1.3641. ”0111111 and ||C2222|| decrease very

slightly with frequency.

The solution of the angle problem is shown in Fig. 4.22. The angles corresponding to the
white areas of the base cell (Fig. 4.19) vary element to element; however, the angles are
around 45 degrees. ||c1m|1 is 3.1006, 2.4704 and 1.0324, whereas ”02222“ is 3.0688,
2.5559 and 0.9919 respectively for the three frequencies. The dynamic moduli are plotted
in Fig. 4.23, where it can be seen that the composite’s dynamic moduli are the same as

the target values. Figure 4.24 shows the iteration history.

 

Figure 4.22. Optimal orientation of material B for example 2 (initial angle 45 degrees)

104

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 . x target ‘ 4 * x target
0 achieved 0 achieved
369 . 34D 1
6 6
2 . - 2 -
l . 9 l ' ®
0 ‘ 7 O ‘ 3
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a)||01111|| (b) ”0222le

Figure 4.23. Plot of the target and achieved dynamic moduli for the angle problem of

example 2 (initial angle 45 degrees)

 

N L») 4:.
O O O

p—a
0

Objective function

 

 

 

O

 

20 40 60 80 100
Iteration no.

0

Figure 4.24. Iteration history for the angle problem of example 2 (initial angle 45

degrees)

4.5.3 Optimization example 3

In this example, the aim is to design a composite that is approximately isotropic. The
absolute values of all entries of the elastic tensor are desired to decrease with frequency.
The absolute values (in 6x1 vector form) are prescribed to be {3.07,3.07,0.85,l .38,0,0}T,
{2.59,2.59,0.71,1.16,0,0}T and {0.99,0.99,0.27,0.45,0,0}T for three frequencies 0, 100

and 2000 Hz respectively.

105

The initial angle for material B is 90 degrees. Figure 4.25 shows the solution of the
density problem. The volume fraction of A is 0.804. The elastic moduli are
{2.8005, 1.5368, 0.5069, 0.7142, 0, 0}T, {1.8973, 1.5345, 0.5069, 0.6682, 0, 0}T and
{0.7868, 1.5299, 0.5092, 0.6060, 0, 0}T, respectively, for the given frequencies. ||cnn|| is
able to soften with frequency because the 90 degree angle orients the negative stiffness
along direction 1. ”c1122” is also able to slightly decrease with frequency. The other two
entries “02222“ and ||c1212|| are constant with respect to frequency as the negative stiffness
oriented only along direction 1 cannot produce softening in direction 2 and along the
diagonals of the base cell. The elastic moduli are plotted Fig. 4.26. “C2222” differs

significantly from the target values. The iteration history is shown in Fig. 4.27.

 

Figure 4.25. Geometry of the base cell designed in example 3 (initial angle 90 degrees)

106

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 ' X target ‘ 4 ’ X target
0 achieved 0 achieved
3JS . 3X 4
C x X
2 - o 1 2 ' ‘
O O O
l ' )5 l t it
0 ‘ A 0 ‘ r
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(3) “01111“ (b) 110222211
4 * x target ‘ 4 ' X target 1
O achieved 0 achieved
3 1 . 3 L 4
2 2 ' i
X X
"E ’ ‘0 ‘
o 6 g; O Q?
0 ‘ 4 0 ‘ r
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(0)1101212II ((1) 110112211

Figure 4.26. Plot of the target and achieved dynamic moduli for the density problem of

example 3 (initial angle 90 degrees)

 

 

 

Objective function

 

 

 

0 20 40 60 80 100
Iteration no.

Figure 4.27. Iteration history for the density problem of example 3 (initial angle 90

degrees)

107

The results are improved further by solving the angle problem. The optimal angles are
shown in Fig. 4.28. The elastic moduli are {3.0780, 2.3906, 0.5839, 1.2802, 0, 0}T,
{2.2849, 1.7438, 0.5321, 0.8471, 0, 0}T and {0.9934, 1.3567, 0.4106, 0.4612, 0, of,
respectively, for the three frequencies. The elastic moduli and iteration history are shown
in Figs. 4.29 and 4.30. Although all non-zero entries soften with frequency, all of them
are a little off the target values. A different initial angle will now be used to check

whether results can be improved.

 

Figure 4.28. Optimal orientation of material B for example 3 (initial angle 90 degrees)

108

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 - x target ‘ 4 ' x target
0 achieved 0 achieved
36D 1 3 { ..
6 o x
2 . . 2 . O 4
CD
1 - GD 1 - 3‘
0 L . 0 . ...
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a) “01111” (b) ”02222”
4 * x target ‘ 4 ' X target
0 achieved 0 achieved
3 . 3 . .
2 * 2 ~
1:1 . 16.5 6 .
O ‘ fl 0 ‘ '
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(C)||01212|| (d) “0112le

Figure 4.29. Plot of the target and achieved dynamic moduli for the angle problem of

example 3 (initial angle 90 degrees)

Objective function

NEH-km
O

u—a
O

 

O

......2_......:......_...... .......

CO

 

 

 

 

O

Iteration no.

0 20 40 60 80 100

Figure 4.30. Iteration history for the angle problem of example 5 (initial angle 90

degrees)

109

Example 3 with different a initial angle

The initial angle for material B is now 1'45 degrees (as shown in Fig. 4.9 in the first
example). The geometry of the base cell obtained after solving the density problem is
shown in Fig. 4.31. The total volume fraction of A in the designed composite is 0.709 (i.e.
close to the lower bound, which is 0.7). For the three forcing frequencies, the elastic
moduli are {I.8132,l.7057,0.5931,0.6245,0,0}T, {l.6670,l.5837,0.3793,0.5557,0,0}T and
{l.3735,l.4097,0.2464,0.3008,0,0}T. The elastic moduli and the iteration history are
plotted in Figs. 4.32 and 4.33 respectively. As can be seen in Fig. 4.32, the elastic moduli

are considerably different from the target values.

 

Figure 4.31. Geometry of the base cell designed in example 3 (initial angle 45 degrees)

llO

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 - x target i 4 r x target
0 achieved 0 achieved

3:: . 3:: .
x x

20 ‘ 2 ' ‘
O O O O O

l a: l )I

0 ‘ 4 0 ‘ *

0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a) “01111“ (13) ”02222“
4 - x target ‘ 4 ' X target l
O achieved 0 achieved

3 . 3 ' 1

2 2 I

1 11 X
X C

0&5 . 9 0 0 . f.) <5

0 10 100 2000 O 10 100 2000
Frequency (Hz) Frequency (Hz)
(C) llclzlzll (d) llanzll

Figure 4.32. Plot of the target and achieved dynamic moduli for the density problem of

example 3 (initial angle 45 degrees)

 

 

Objective function

 

 

 

0 20 40 60 80 100
Iteration no.

Figure 4.33. Iteration history for the density problem of example 6 (initial angle 45

degrees)

111

The optimal angle, i.e., the solution of the angle problem is shown in Fig. 4.34. Figure
4.35 shows the corresponding dynamic moduli. The dynamic moduli are approximately
the same as the target values and therefore soften with frequency. The achieved values
are {3.0654, 3.0856, 0.8579, 1.3966, 0, 0}T, {2.5122, 2.4920, 0.6821, 1.1010, 0,0}T and
{1.0405, 1.0349, 0.2856, 0.3382, 0, 0}T at the three frequencies. The iteration history is

shown in Fig. 4.36.

 

Figure 4.34. Optimal orientation of material B for example 3 (initial angle 45 degrees)

112

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

4 1 x target 1 4 - >< target
O achieved 0 achieved
349 . 36D 1
6 6
2 ' 2 .
l b e l , f
0 ‘ “r 0 ‘ r
0 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(a) ”Gun” (13) “02222”
4 - x target 1 4 1 x target
0 achieved 0 achieved
3 . . 3 . .
2 - 2 '
6D
165 . 1 . 6 .
6 E
0 _ - d 0 . in 45
O 10 100 2000 0 10 100 2000
Frequency (Hz) Frequency (Hz)
(0) llclzlzll (d) “01122“

Figure 4.35. Plot of the target and achieved dynamic moduli for the angle problem of

example 3 (initial angle 45 degrees)

KI}
O

 

NW4}.
COO

Objective function

—0
O

 

O

 

 

 

20 40
Iteration no.

0

80 100

Figure 4.3 6. Iteration history for the angle problem of example 3 (initial angle 45

degrees)

113

4.6. Summary and Discussion

In this chapter, two-phase composites that show frequency-induced softening are
designed. One phase, phase A, is a typical viscoelastic material. The other, phase B, is a
negative-stiffness, lumped system. A typical topology optimization is used as a tool to
design the composite. However, the present problem is considerably more complicated
than a typical material design problem. In a typical material design problem, only the
element densities (p) are design variables. In the present work, the angles of rotation (0)
of material B in the elements are also treated as design variables. As in a typical material
design problem, the objective is to minimize the least square error between the prescribed
and the actual elastic tensors and there is a constraint on the volume fraction of material
A in the base cell. The elastic tensor is also constrained to be positive definite in order to
ensure stability of the material and this constraint is not needed in a usual material design

problem. As the effective properties of a composite having negative stiffness inclusions
are very sensitive to the design variables (0,, and 68), the present problem is found to be

very difficult to solve using a typical topology optimization method. In view of this, the
design problem is proposed to be solved in two stages. In the first stage, initial angles of
rotation of material B are guessed. The angles are kept constant while element densities
are allowed to vary. In the second stage, the element densities obtained from the first
stage are kept constant and the angles are allowed to vary. If the initial angles of rotation
are judicially chosen, it is possible to obtain the desired properties in the first stage itself.
In a complex problem, where the initial guess is difficult, the angles may be chosen as

i45 degrees (as in Fig. 4.9) as this angle is generally found to give good results. This is

114

because this angle orients the negative stiffness along the diagonals of the base cell,
which not only helps softening in shear but is also found to help softening in directions 1
and 2. In all the three examples, the initial angle of i45 degrees resulted in the effective
properties very close to the target values. The results obtained from the initial angle of 90

degrees are generally found to be a little off from the target values.

The presence of the negative stiffness material phase in the composite is found to lead to
difficulties in obtaining purely black and white geometry of the base cell. For this reason,
an advanced filtering method such morphology-based black and white filtering (using
Heaviside functions) has been used in the present work. There are, however, still a few
gray elements in the base cell designed. These gray areas may be eliminated by refining
the discretization of the base cell and running the density optimization routine for more
number of iterations. Refined mesh will also give mesh-independent, smooth geometry of

the base cell. Refining mesh, however, decreases time efficiency.

115

CHAPTER 5

SYNTHESIS OF TILEABLE BISTABLE STRUCTURES

In the previous chapters, two-phase composites that exhibit frequency-induced softening
have been designed. Once layout of the composite is known, the next step is to physically
realize the composite material. The matrix phase may be a typical rubber-like material,
which is easily available. The negative-stiffness inclusion phase, on the other hand, is not
available and therefore needs to be physically realized first as a first step towards
realizing the whole composite material. As introduced in chapter 2, bistable structures
have potential to provide the composite with the required negative stiffness. The lumped
lattices of material 131 or B proposed here (as given in sections 2.6 and 2.7) have bistable
structures as well as typical elastic/viscoelastic structures. One way to synthesize such a
lumped lattice is to build the components separately and then assemble them together.
The lumped system of material B1 or B is proposed to be built in two steps — (i) build
two-dimensional (tileable) bistable structures and (ii) interconnect the tileable bistable
structures with blocks of typical elastic/viscoelastic materials such that the resulting
structures is equivalent to the lumped system proposed in sections 2.6 or 2.7. An example
of the resulting structure is shown in section 2.8. As the typical elastic or viscoelastic
components of the lattice system are easily available, physical realization of the lumped
system boils down to the physical realization of tileable bistable structures. Accordingly,
this chapter concentrates on the synthesis of two-dimensional bistable structures.
However, the problem of designing bistable structures will be treated as independent and

general purpose. As such, this chapter may seem to be a little disconnected from the other

116

chapters because of the independent treatment; however, this work is important for the

physical realization of the negative-stiffness materials.

As pointed out earlier, the goal of this chapter is to design a two-dimensional, tileable
bistable structure that can be used as a building block of the negative stiffness phase
required to achieve frequency-induced softening. The 2D bistable structures are designed
using a structural topology optimization technique, where the objective is to maximize
the spatial difference in the two stable configurations of the bistable structures (recall that
bistable structures have two stable configurations with no external loads). The geometry
of the bistable structure is represented by a prescribed set of interconnected beams. The
set is known as a ground structure. The aim of the optimization is to remove the
unwanted beams from the ground structure so that the remaining interconnected beams
functions as a bistable structure. The middle region and the end regions are allowed to
have different cross-sectional areas (or materials). A thin cross section or a flexible
material at either end (as compared to a thick or stiff material at the middle) of a few
beams in the structure is a key to achieve bistability. Accordingly, cross-sectional areas
of the end region of the beams are treated as design variables. As an alternative, the
materials of the end regions are also treated as design variables. The Optimization
problem is solved using a stochastic optimization method, viz. a genetic algorithm (GA).
A penalty function is devised and incorporated in the objective function in order to avoid
an invalid structure as a solution. A typical bistable structure synthesized in this work is

as shown in Fig. 5.18.

117

5.1 Introduction

A bistable compliant structure has two stable equilibrium configurations when unloaded.
Once such structure reaches one of its stable configurations, it remains there unless it is
provided with enough energy to “climb” out of an energy well that leads into the other
stable configuration. This feature of bistable structures can be advantageous in designing
a variety of mechanical devices, such as switching devices in MEMS, relays, valves, etc.
The term bistable periodic structure is used to refer to an arrangement of interconnected
bistable structures that tile a plane periodically. The load need not be periodic and the
tiling need not be infinite. This chapter discusses a computational strategy to design
bistable periodic structures using topology optimization. Understanding the behavior and
developing a methodology for design of bistable, periodic microstructures is a first step

towards design of bistable materials, a principal motivation of this work.

A typical example of a bistable compliant structure is a bistable compliant mechanism
(e. g., see Jensen et al. (1999), Jensen et a1 (2001), Jensen and Howell (2003), Masters and
Howell (2003), King and Campbell (2004) ). Such mechanisms typically rely on strain
energy storage to gain bistable behavior. They are made of flexible members, which
produce motion and at the same time, store energy. Compliant mechanisms in general
(e.g., see Howell (2001) ) offer several advantages when used in design, specially in
design of MEMS, and have been the subject of extensive research. Particularly relevant

to this effort is previous work where a topology optimization strategy is used to

118

synthesize compliant mechanisms, e.g., as discussed in Ananthasuresh and Kota (1995),

Sigmund (1997) and Frecker et a1. (1997).

A distinguishing feature of any bistable compliant structure is the shape of its strain
energy curve. Figure 5.1 shows a qualitative description of the variation of strain energy
with strain for a typical bistable compliant structure. The firnction is non-convex. The
strain energy curve has three critical points - two minima (points C and G) and one
maximum (point B). When the system is unloaded, and under small loads, the structure
will operate in and around one of the two minima, corresponding to one stable
configuration. Transition from one configuration to the next requires the addition of
sufficient energy to jump over the small maximum and over into the neighborhood of the
other minimum. An external load provides this activation energy in the form of external
work and switches the structure from one stable configuration to the other. Once the new
state is reached, again the structure will operate in and around this configuration, for
small enough inputs. As indicated in Fig. 5.1, there are other characteristic points on the
strain energy curve, which are described in detail in section 4.2. The shape of the strain
energy curve around these points gives quantitative measures of the bistability. These
points will be used to produce a measure of performance that can guide an eventual

methodology for optimization of bistable periodic structures.

119

 

 

 

I
56
5 H
.5
cu
a
m
t 5
G Strain

Figure 5.1. Strain energy versus strain in a typical bistable structure

The design of compliant mechanisms typically follows one of two principal approaches.
The kinematics approach produces a structure composed of small flexible pivots and
relatively rigid links. This approach results in lumped compliance, such as locally
deforming flexural joints. The structural topology optimization approach is based on the
methodology for topology optimization of structures introduced by Bendsee and Kikuchi
(1988) and uses either a continuum modeling based on plane elasticity (e.g., Sigmund
(1997) and Bruns et al (2002) ) or a ground structure approach that relies on slender

members (beams and bars, e. g., F recker et al (1997), Joo and Kota (2004) ).

This chapter investigates the automatic topology synthesis of two-dimensional bistable
compliant, periodic structures using a ground structure topology optimization approach.
Section 5.2 introduces a concept for a bistable structure that can be repeated periodically
to tile the plane. This concept is used simply to highlight basic features that one may
expect to find in the layout of these structures. Section 5.3 describes the finite element

model used for analysis of compliant bistable structures. A possible topology

120

optimization problem formulated to synthesize a compliant bistable structure is described
in section 5.4. Two examples in section 5.5 and one example in section 5.6 illustrate the

methodology.

5.2 An Example of a Bistable Periodic Structure

This section introduces a design concept that can be repeated periodically to tile the plane
and form a bistable periodic structure. This concept is used to highlight basic features of
these structures that may be used later in deciding how to build a ground structure for

topology optimization.

In Qiu et. a1. (2004), the authors present a monolithic mechanically-bistable mechanism
that uses no latches, no hinges and no residual stress. The typical implementation of this
mechanism involves two curved, centrally-clamped parallel beams, labeled “double
curved beams”. Figure 5.2 shows a simplified “double curved beam” bistable mechanism,
where each of the curved beams (shown as dashed lines) is approximated as two straight
beams. This simple concept may be used to build a two-dimensional bistable structure, in
this case, using four such bistable mechanisms, as shown in Fig. 5.3. In the figure, one of
the four component bistable mechanisms is encircled by a dashed line and shown alone in ,
Fig. 5.4. The original central clamp (the dark, vertical line in Fig. 5.2) is replaced by a
stiff triangular structure in Figs. 5.3 and 5.4. This allows the incorporation of two input
points, where loads are applied on each side of the structure. In addition, the beams are

thickness modulated, i.e., the cross-sectional area of either beam is allowed to vary along

121

its length to provide more control on the bistability characteristics of the structure. A
simple way of implementing a thickness modulation considered here is to use beams of

two different cross-sectional areas.

Bistability is achieved by designing the structure such that some of the members buckle
or “snap”. In the structures shown in Figs. 5.3 and 5.4, members which may snap are
shown as light gray lines. Members shown as dark lines are stiffer, stiff enough to
provide clamping support to the snapping beams. The configuration shown in Fig. 5.3
corresponds to the first equilibrium configuration of the structure. When loaded as shown
by a load of sufficient magnitude, the structure moves into its second equilibrium

configuration, shown in Fig. 5.5.

 

Figure 5.2. Simplified “double curved beam” bistable mechanism

 

Figure 5.3. A bistable structure (first stable configuration) based on four “double

curved beam” sub-structures.

122

Q.

Figure 5.4. A sub-structure of Fig. 5.3 that is a bistable mechanism

ll

Figure 5.5. A bistable structure (second stable configuration) based on four “double

curved beam” sub-structures

 

(b)

 

(a)
Figure 5.6. A 3x3 periodic arrangement of bistable structures (a) first equilibrium

configuration (b) second equilibrium configuration

The bistable structure shown in Fig. 5.3 can be repeated periodically to tile 3 plane,

joining each unit together via rigid connectors (shown as solid rectangles). Two stable

123

equilibrium configurations of such structure are shown in Fig. 5.6. The dashed box

shows the fundamental cell of the periodic structure.

The solution proposed in Fig. 5.3 was generated by analogy, using the “double curved
beam” structure as guide. A computational strategy is being sought that may allow to
synthesize other such structures and meet some prescribed performance specification,
e.g., match certain desirable features of a strain energy curve or a load deformation
diagram. One such strategy based on topology optimization is discussed in the following

sections.

5.3. The Analysis Model

The structure is analyzed using non-linear, corotational Timoshenko beam elements and a
total Lagrangian formulation (a detailed formulation is given in Crisfield (1991), p219.)

The analysis is quasi static and all inertial forces are ignored. At equilibrium,

g<x> =fint(x)"fext(’) = 0 (5.3.1)
where
fint (X) = 6% (5.3.2)

is the internal force, A is the total strain energy, x(t)e 91" represents the configuration of

the structure at time tunder an external force fax, (t) and n is the total number of degrees

124

of freedom. A Newton-Raphson-like scheme is used to find a solution to (1). At the i-th

iteration step the scheme produces

Xi+1 = X? -[a—] , 8(Xi) (5.3.3)

where x6“ =x'. In (5.3.3), fig/ax is the tangent stiffness matrix. x’ represents the

solution to the equilibrium equation (1) at timet.

Simulation results were validated using the commercial finite element package ABAQUS
using the element ‘b21’, a beam element. As the model consists only of beam elements,
joint rigidity may be overestimated, when compared to a full 2D elasticity model (e.g.,
one that uses quadrilateral finite elements). This drawback can be overcome by the

modulation of the cross-sectional area, as discussed in the next section.

5.4. Optimization of the Topology

5.4.1 The ground structure and design variables

The formulation used here is similar to a classical ground structure approach in truss
topology optimization (e.g., Dorn et a1 (1964) or Bendsee et al (1994) ). In a typical
ground structure approach, design variables are cross-sectional areas Aa of each member

(i.e. bar) in the ground structure, where a=l,2,...,nb and nb is the total number of

members in the ground structure. The ground structure contains enough members and

125

possible connections to include, as a subset, a reasonable number of design alternatives.

The moment of inertia of the beam Ia, is assumed to be proportional to the square of the

beam’s cross-sectional area, i.e.,

Ia, = 0.43, c>0 (5.4.1)

A ground structure made up of beams of uniform cross section will not have enough
flexibility to produce a bistable structure (Qiu et. a1. (2004)). In order to provide joints
with flexibility, a measure of compliance is achieved by “modulating” the cross section
of elements along their axis, reducing the cross sectional area of elements near a flexible
joint (alternatively, a soft rubber-like material may be used to model hinges, as illustrated
in Section 5.5). In the present work, each member of the ground structure is divided into

three regions, as shown in Fig. 4.1. The two end regions are identical in length (L1) and
cross-sectional area (A1). The middle region has length L2 and cross-sectional area A2. In
a ground structure, A1 is an independent (design) variable; A1 is allowed to take one of
three possible values, corresponding to: member removed (Al=0), member is a snap-
beam (A1 = A1) or member is a support beam (A1 = A" ), where 0 < A] < A“ and A] and
A” are prescribed values. A2 is a dependent variable, set to 0 (member removed)

whenever A1 =0 and set to A” otherwise. This choice is motivated by the example
bistable structure introduced in Section 5.2 which was constructed using only two types
of beams: a more compliant one, which snaps and is the source of the bistability, and a

stiffer one, which provides support to members of the first type.

126

L1 L2 L1

          
     

A1 A A1

Figure 5.7. A member (bar) of the ground structure

Based on the previous discussion, the cross-sectional areas Aa of the end regions of each

member (1 in the ground structure, a =1, 2,..., ”b , are selected as design variables in this

problem and letAa take values in {O,A',A"} .

5.4.2 The objective function

The objective function used here is built around characteristic points of the load
displacement curve. A typical curve is shown in Fig. 5.8. Points D and F correspond to
values where the load just reaches the critical value that causes snap-through. Under
external loading, B and H are points where the structure settles just afier the snap
through. Starting from the first stable configuration (point C), when the bistable structure
is loaded past the critical level for snap-through, the state of structure follows the path C-
D-H-I, by-passing the segment D-E-F-G-H along which the internal force in the structure
cannot equilibrate the external load. When the structure is unloaded after reaching 1, it
settles in the second stable configuration at G, afier following the path I-H-G. Under
similar loading in the reverse direction, the structure follows the path G-F-B-A, now by-

passing the segment F-E-D-C-B.

127

An objective function is now proposed, which can be used to evaluate the quality of
alternative designs. Here a very simple objective function is used, focusing primarily on
the difference between the two stable configurations, roughly, a measure of the
magnitude of the snap-through jump in a simple snap-through structure. When xC and x6

are two stable configurations of a simple snap-through structure, this measure is

¢=|| xC —XG ll”2 (5.4.2)

 

§
[2 I
D
------- > - - - - - - - - H
C E G
>
V Displacement
B - ------ (— -------
F

A

 

Figure 5.8. Typical force-displacement diagram for a bistable structure

In more complex structures, the direction of the snap-through may play a role and to

account for this, the objective function is defined as

 

0 = (xC - xG)T W0C - x6) (5.4.3)

128

where W is a diagonal, square, semidefinite matrix. Finally, most subsets of the ground
structure will not be bistable and in fact, the majority of such subsets will not even be
properly formed structures (e. g., they may be disjoined or not properly supported). For
any such structure, (D is set to zero. Thus, the goal achieved by maximizing (I) is simply
to find any one bistable structure contained in the ground structure. The function ¢ is the
objective fimction to be maximized in the optimization problem, stated in the following

section.

Loading A typical external loading is given by

f,,, (t) = r(t)f0 (5.4.4)

where 1(t) is as shown in Fig. 5.9. f0 is a prescribed maximum external force and t is a
pseudo time. The time discretization is chosen to suit the finite element analysis. The goal
is to obtain a bistable structure with a critical load less than f 0. A displacement driven
problem (one where the displacement at the input port is prescribed) is solved to capture
the segments D-E—F-G-H and F-E—D-C-B of the load-displacement curve (Fig. 5.8),

which are by-passed in a load driven problem. In that case the displacement x," at the
input port(s) is prescribed by

xm(1) = T(’)x0 (545)

where x0 is a reference displacement amplitude. The force in the ordinate of the force-

displacement curve corresponds to the corresponding reaction force.

129

 

Figure 5.9. The external load factor as a function of time

5.4.3 The optimization problem

The optimization problem is formally stated as:

Find A = {A1,A2,...,A,,b } e {0, A’ ,A"}"b that

 

maximize ¢ = \/(xC - xG)TW(xC — xG)

subject to

(5.4.6)
nb
Z Z(Aar) S ”MAX
a=l
where
(A )_ O for Aa=0
I a — 1 for A0, >0

Data for this problem are:
A] cross-sectional area of the thin beams

A" cross-sectional area of the thick beams

n MAX the maximum number of members allowed

130

W prescribed weight factors used to emphasize a given snap—through

direction

5.4.4 The solution scheme

Problem (5.4.6) is treated as a combinatorial optimization problem and solved using a
genetic algorithm (GA). However, since many, if not most of the subsets of the ground
structure are not bistable and may not even be well formed structures, a significant
portion of the design space has objective function values ¢ =0. In order to provide for a
meaningful evaluation of the merit of such structures, the original objective in (5.4.6) is

replaced by the minimization of the following merit function in the GA:

f=-(w¢)2+w (5.4.7)
In (5.4.7)
11’ = W1 +W2 + W3 +1114 (543)

is a penalty function with entries m, {1}; , W3 , and ya designed to provide all sub-sets of
the ground structure with a meaningful merit function value (even sub-sets that are not

well formed structures). The (constant) scaling factor w is scaled so that we) z 1 , as are

W1, 1;!) , tog , and 1%- These fimctions are defined as follows:

131

 

 

(1) W1 penalizes structures that do not have any structural members at the prescribed

loading ports (nodes). {(11 is defined as

n
1
W1 = :5— 2 5,- (5.4.9)
InaX jzl
where 51- is the distance between the structure and the j-th loading port and 5m“ is the

maximum distance between two nodes in the structure.

(2) W2 penalizes structures with too many bars. V2 is defined as

 

If n >"MAX
1V2 = (n: -nMAx) " (5.4.10)
0 if no S ”MAX
where
nb
n, = Z 1046,) (5.4.11)
a=l

is the actual number of members in the structure and nr>nmx is a prescribed number used

 

to adjust the slope of the penalty and selected to keep the value of {112 near 1 (e.g., n, = 2
nMAX). A lower (or higher) 11, increases (or decreases) the relative weight of tyz with

respect to the other (If; s.

132

(3) W3 penalizes disjoint structures. V3 is defined as

W3 = _s (5.4.12)

where N, is the number of disjoint sub-structures and N, is a prescribed scaling factor

(e.g., N52). To measure N, the structure is represented as an undirected graph. N, is the

number of components of the graph.

(4) V4 penalizes structures that are not properly supported i.e. structures that may

undergo rigid-body motions. (114 is defined as

1,114 =Nx+Ny+N9 (5.4.13)

where N, =0 if the structure’s rigid-body motion in the x direction is constrained, Nx=1
otherwise, and similarly for Ny and N61 Only a straight forward evaluation is made,

simple enough to rule out designs with little computation. For instance, Nx is set to 0 if a

node connected to the structure is on an x- constrained boundary.

133

5.5. Examples

 

 

 

 

 

5.5.1 Example 1
0.6mm (3min 0.6mm
9. ~. X: A
a 1
‘0. ‘ |
O ‘\ |
”B : B"
at
0B ’1’ : \ Bu
8 ” I
E , I
c: : ‘.
° .3. 1 ,3. ~
(a) (b)

Figure 5.10. Example 1. The package space (a) and corresponding ground structure (b)

with dimensions and boundary conditions

In this example the goal is to design a two-dimensional bistable structure that will
operate, loaded vertically at port A and horizontally at port B and fit within a 1.6mm x
1.6mm package space, as shown in Fig. 5.10. This space is partitioned into 8 triangles,
whose boundaries are marked with dashed lines in the figure. The structure is assumed to
be symmetric about these lines. The ground structure is laid on one of these triangles and
the appropriate boundary conditions are applied to enforce symmetry (Fig. 5.10(b) ).
Loads F =6‘c(t) in the range i6 mN are applied at positions A and B. The function 1:(t) is
as in Fig. 5.9. All entries in the weight matrix W are zero except Wm=1, where m is the
degree-of-freedom associated with the load F. This means that snap-through in the

direction of the external load is needed to be maximized.

134

The total ntunber of bars in the ground structure (nb) is 62. In accordance with Section

5.4.1, every member is discretized into three regions with lengths L2=6L1 and every
region is further discretized into two finite elements (Fig. 5.7). A solution is allowed a

. . . -4 2 u
maxrmum of nmx=8 bars, each usrng one of two cross section areas: AI=10 mm and A

:103 m2. The Young’s modulus and Poisson’s ratio are 1380 MPa and 0.3,

respectively.

The ground structure depicted in Fig. 5.10(b) admits, as one solution, the topology of the
bistable structure introduced Section 5.2 (Fig. 5.3), the reference structure. That structure
has an objective function ¢=O.3674 (mm) (this value is used to scale the objective
function f in (5.4.7) and defines w=1/0.3674=2.7218). The structure has the force
displacement diagram shown in Fig. 5.11. The switching force (the force that causes
snap-through) for this mechanism is approximately -3.5 mN in the forward direction and

1.8 mN in the backward direction.

 

Force (mN)

 

 

 

-6
-0.5

 

Displacement (mm)

Figure 5.11. The force-displacement diagram for the reference structure (Fig. 5.3)

135

 

 

Now a GA is used to find other solutions. The GA starts with a population of 100
randomly generated designs, which are sorted based on their merit function (f) values.
Following a standard GA procedure, every individual is assigned a scaled merit function
value depending on its position. An elitist strategy with generation gap of 60% is used.
The best 40% of the solutions in the population are carried forward to the next generation
and the remaining 60% is produced through crossover among the parents, selected
probabilistically from within the current generation. The probabilistic selection of parents
for crossover is based on the standard “roulette wheel” mechanism, where the probability
of selection of an individual for the crossover is directly proportional to its fitness value.

The mutation rate used is 10%.

The process of fitness assignment, crossover, mutation etc. was repeated for 500
generations, which resulted in the solution shown in Fig. 5.12. The objective function
value for this structure is ¢=0.3688 (mm), with penalty function {71:0, i.e., only slightly
better than the structure in Fig. 5.3. The corresponding force displacement diagram
(obtained with the prescribed input displacement) is shown in Fig. 5.13. The switching

force is approximately -1.8 mN in one direction (forward) and 0.9 mN in the other

(backward).

136

 

ff? v 4

LA

)

 

 

 

 

 

 

 

 

 

 

 

 

(a)
Figure 5.12. Example 1. Bistable structure obtained by the GA (a) first stable

configuration (b) second stable configuration

 

Force (mN)

 

 

 

  

38.5 -04 71.3 —0.2 -01 0 0.17 0.2
Displacement (mm)

Figure 5.13. The force-displacement diagram for the structure in Fig. 5.12.

 

Another solution with essentially the same objective fimction value (and t/FO) is shown
in Fig. 5.14. In this structure, however, several members overlap. The present modeling

does not account for interference or contact between members, as the structure switches

137

from one configuration to the other. This is acknowledged as a drawback, unfortunately

one that is difficult to overcome. This will be discussed further in the next section.

 

(a) (b)

Figure 5.14. Example 1. Another bistable structure (a) first stable configuration (b)

second stable configuration

Solutions in Figs. 5.12 and 5.14 show that the deformation in the hinges can be quite
large, raising questions regarding the manufacturability of solutions modeled using these
ideas. Some of these concerns may be alleviated by using a rubber-like material
(polyurethane elastomer) at the hinges. The maximum allowable normal stress in this
material is about 45 MPa, which corresponds to about 440% elongation at break For
instance, returning to the structure shown in Fig. 5.14, one can use the soft material in the

two end sections of each hinged member (shown light gray in the figure) and keep the
area constant at A" over the length of the member. This would result in the load-

displacement curve shown in Fig. 5.15. Note that using the rubber-like material does not
change the size of the jump between stable configurations significantly (but it does

change the magnitude of the actuating force required to switch configurations). The

138

maximum stress in the elastomer is 25.5 MPa, within the allowable limit for the material.
This suggests that using a combination of materials that includes softer materials near the

hinges may result in concepts that could be manufactured more easily.

 

Force (mN)

 

 

 

-6 1 1 1 L
-0.5 -0.4 -O.3 -0.2 -0.1 0 0.1 0.2
Displacement (mm)

 

Figure 5.15. Force-displacement diagram for the bistable structure obtained by the GA

All the solutions shown here are obtained from the same run of the GA. Figure 5.16
shows the history of the best GA merit function (I) value for a generation. At generation
0, i.e. in the randomly generated population, the best merit function value is f =0.1250. At
the end of 500th generation, the best merit function value is f =-1 .008. The finite element
analysis and optimization procedures were implemented on a MATLAB program running
on a Pentium 4 machine with 2.8 GHz clock speed, 496 MB RAM and Windows XP
operating system. 93 hours were spent in running the Optimization program for 500

generations.

139

 

0.2 T I I 1

0
é—oz
>
8
'5 —0.4
o
E
o “0.6 '
.2.
73
3-08
0

 

 

 

 

—1.2 ‘ 1 ‘
0 100 200 300 400 500
Generation Number

Figure 5.16. History of the best GA merit function value f

 

 

 

 

 

 

 

 

 

 

 

5.5.2 Example 2
0.4
0 6 mm mm 0.6 mm A

E ‘\\\ A E A ”,I
E \“ ' I”
”’O. \ | ,’
o \\\ I III

B : B

8 E -- 74:, -- -

B : B
E ,/ 1 \\
E 1 \
so / ' \
0' A : A

(a) (b)

Figure 5.17. Example 2. The package space (a) and corresponding ground structure (b)

 

with dimensions and boundary conditions

140

The package space in this example is the same as that in Example 1. The goal is to
maximize the snap-through at ports B in the horizontal direction, when the structure is

loaded (only) at ports A, in the vertical direction.

A load F in the range of i4 mN is applied at each loading port. With the given loading,
the structure is assumed to be symmetric about the horizontal and vertical centerlines,
which are shown as thick dashed lines in Fig. 5. 17 (a). In addition, the structure is desired
to be orthotropic. For this reason, the layout is assumed to be symmetric about the
diagonals, shown as thin dashed lines in Fig. 5.17(a). Accordingly, the package space is
partitioned into 8 equal triangles, whose boundaries are marked with dashed lines in the
figure. The ground structure is laid on one of these triangles. The ground structure is the
same as that in Example 1, as are other details of the discretization. It should be noted,
however, that since the loading is symmetric only about the centerlines, one quarter of the
structure is analyzed ( Fig. 5.17(b) ). Appropriate boundary conditions are applied to

enforce symmetry, as shown in the figure.

In this example all entries in the weight matrix W are zero except pr=l, where p is the

degree-of-freedom associated with the horizontal motion of the output port B. The GA
parameters such as population size, mutation rate and generation gap have been kept
same as that in Example 1. The GA was run for 500 generations, after which a solution is

obtained with merit function value f = -0.10614, ¢= 0.1197 (mm) and t/FO. The solution

is shown in Fig. 5.18(a). Fig. 5.18(b) shows the corresponding second stable

configuration of the structure.

141

 

The external force at input port A is plotted versus the diSplacement at the output port B
in Fig. 5.19 (obtained fiom a displacement-driven analysis). The maximum switching
forces are approximately -1.9 mN in the forward direction and 0.5 mN in the backward
direction. The maximum switching force in the backward direction is small compared to
that in the forward direction. While the present model does not penalize this difference,
control of the relative magnitudes of the switching forces can be accommodated in the

GA formulation by simply adding an additional penalty term to the objective function.

In Fig. 5.19 it can be noticed that the displacement at the output port B does not change
sign through the input cycle, i.e., as t(t) (Fig. 5.9) goes through one full cycle. The
force-displacement curve is nearly vertical in the neighborhood of the first stable
configuration, suggesting that in this region the output port remains nearly stationary for
a small range of inputs, a behavior reminiscent of a Poisson’s ratio zero material. This is
of course not true in a neighborhood of the second stable configuration, where a small

external force applied at the input port makes the output port B move a finite distance.

 

(a) (b)

Figure 5.18. Example 2. Bistable structure obtained by the GA (a) first stable

configuration, (b) second stable configuration

142

1‘17st 137——

 

 

Force (mN) at the Input Port A
O

 

   

 

   

 

 

_4 i 1 1 L J
-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3
Displacement (mm) at the Output Port B

Figure 5 .19. Example 2. Plot of the internal force at the input port A versus the

displacement at the output port B.

Finally, again it can be observed in the solution proposed to Example 2 that several
members overlap as the structure moves from one configuration to the other, signaling
interference. In this problem, however, no solution was found where interference did not
occur. One way to address this problem is to cast the problem in three-dimensions,
allowing at least some members to have a small curvature out of the plane, to avoid
interference. This can be appreciated in Fig. 5.20, which shows quarter sections of the
solution in its two stable configurations. Here, as in all previous solutions, one can
identify three main sub-structures. One sub-structure, formed by stiff members, remains
fixed in space throughout the deformation and acts as support to the rest of the structure.
A second sub-structure moves essentially as a rigid body and is also formed by stiff
members. This sub-structure typically receives input from or transfers output to

neighboring components. The third sub-structure connects the previous two and is

143

T ‘T

-. 7

‘-

 

formed by compliant members. This sub-structure is primarily responsible for the
motion of the mechanism. In most cases interference is avoided if these three
components are in different planes of the structure. This qualitative description provides
hints on how to actually implement these solutions. Naturally, this is not sufficient in
many applications, e.g., in very small scale designs where a layered solution is not
realistic. In such cases a fill] implementation of a contact detection scheme that rejects

solutions where interference is detected would be necessary.

9' “

......

   

   
 

Compliant
Members

\ Fixed

(a) Members 0»

Figure 5.20. Detail of the solution of Example 2. (a) First stable configuration (b)

Second stable configuration.

5.6. Rubber Model

As also discussed in Section 5.5.1, deformation in the hinges can be quite large, raising
questions regarding the manufacturability of solutions modeled using two cross-sectional
areas. Some of these concerns may be alleviated by using a rubber-like material (e.g.,

polyurethane elastomer) at the hinges. Using a combination of materials that includes

144

 

 

 

softer materials near the hinges may result in concepts that could be manufactured more

easily. Some works have been carried out in this direction and presented next.

As earlier, every bar in the ground structure is divided into two kinds of regions — shown
as dark and light gray regions in Fig. 5.21. Material properties for the region shown in
light gray are variable, taking three possible values, namely, values corresponding to ‘no
material’, a ‘soft’ material and a ‘hard’ material. The region shown dark is made from

either ‘no material’ or the ‘hard’ material.

./

Figure 5.21. Rubber model: A member (bar) of the ground structure

Here the soft material is a compliant, rubber-like material, e.g. natural rubber, silicone
rubber or a polyurethane elastomer. The hard material is relatively stiff and

approximately linear material, e. g. polypropylene.

Material properties are characterized by — (a) the normal stress (0') vs. normal strain (8)
relationship and (b) the shear stress (2') vs. shear strain (7) relationship for the material.
Figure 5.22 shows a representative plot of 0' vs. 8 in a soft material, viz. polyurethane

elastomer. The rvs. yrelationship for both hard and soft material is assumed to be linear.

145

 

Stress (M Pa)

 

 

 

,.. .............. ,,,_".,, .... ,........ .

‘40 L 1 i 1 1
-2 -l 0 l 2 3 4 5
Strain

 

Figure 5.22. Prescribed normal stress as a function of normal strain in the soft material

The material and layout of the structure are controlled by two binary design variables for

each bar a, labeled pa, and #0,. Variable pa controls the layout of the structure: pa, =1
means bar a is present, Pa =0 means bar a is removed. When a bar is present, its
material is controlled by variable ya. #0: =1 indicates that the light gray regions of bar
aare made of the soft material, i.e. abar is hinged at both ends. ya = 0 means that bar

a is made of the hard material. Thus, the normal stress in the light gray region of bar a

may be written as

05"”(panzme) = pa(/1a01(8)+(1- #a)02(6)) (5.6.1)

146

where the superscripts (1) and (2) refer to properties of the soft and hard material,
respectively. Similarly, the shear stress in the light gray region of bar amay be written

as
réraymmtza, 7) = pa (4.2ar +0 wanztn) (5.6.2)
The normal stress in the dark region of bar amay be expressed as

oéa’ktpa,e)=pa02<e) (5.6.3)

while the shear stress in the dark region of bar a is

Tgark (pan 7) = pal-2 ( 7) (5.6.4)

A scaled arc-length method ( Al-Rasby (1991), Bruns et. A1. (2002) ) has been used for
the solution of the non-linear equations. Arc length method is suitable for analyzing
structures undergoing instability (snap through or buckling). Such instability is typically
undesirable structures ( e. g. in Neves (1995) ), whereas such instability is desirable for

bistable structures.

5.6.1 New optimization problem

The new optimization problem is formally written as

147

 

Find 0 = {p1,p2,.-.,pnb}e {0,1}"b. u = {41,42,-..,/1nb}e {0,1}"b that

maximize ¢(p, u) = (xC — xG)T W(xC - xG) (5.6.5)

subject to

"b
2 par S ”AJAX
a=l

5.6.2 Example

The example in Section 5.5.1 is solved using the Rubber Model. The cross-sectional area
is kept constant at 2x10"3 mm2 over the length of every member. The hard material is
polypropylene with Young’s modulus 1380 MPa and Poisson’s ratio 0.3. The soft
material is assumed to be linear in shear with shear modulus 5.2 MPa. The 0-8 plot for

the soft material is given in Fig. 5.22.

The solution obtained by the GA is shown along with its second stable configuration in
Fig. 5.23. The objective fimction (05) value for this structure is 0.3629 (mm), or merit
function f is -0.9758 with penalty function W0 The corresponding force displacement
diagram is shown in Fig. 5.24. The critical force required for the snap through is

approximately 5 .5 mN in one direction (forward) and 3.0 mN in the other (backward).

148

 

 

 

 

 

 

 

 

 

 

(a) (b)

Figure 5.23. Bistable structure obtained by the GA (a) first stable configuration (b)

second stable configuration

 

 

6 I T ' ' 1 Y
. 1'
. l
4“ ' 1" a
. . 1
”<3”, “““““ ) "“"T‘/
A 2- .. .//. \ . If ..,
Z / \ / 3
g / \ 1
3 0— / \ ,1 2
s , ( . 1'
”- \ . 1'
-2 ~ 1 \ .> [,1 2
.' . \\ .
/ s \, / 4
\ 1
/ \ /
f < . q
0 0.1 0.2

 

 

 

-6 1
-0.5 —0.4 —0.3 -0.2 -0.1
Displacement (mm)

Figure 5.24. The force-displacement diagram for the structure in Fig. 5.23

It can be noticed that the structures shown in Figs. 5.12, 5.14 and 5.23 are the same type

as they essentially only have the 3-beam segments on peripheral as bistable elements. In

149

fact, essentially all solutions found work on a similar principle: they exploit the bistability
of simple two- or three- bar segments placed strategically around a supporting frame.
This fact can be exploited to develop different and perhaps more efficient design
strategies. For example, with this insight one could devise a different algorithm, perhaps
one in which the design space is represented not by means of a ground structure, but
instead by a vocabulary that involves only two bar sub-structures and structures to
support them (and a grammar that works on connecting each piece to each other and to

the external loads). This new approach is studied in Prasad and Diaz (2005).

150

 

 

CONCLUSIONS AND FUTURE WORK

This work investigates methodologies to design a structure or material whose dynamic
modulus would decrease with forcing frequency. The motivation to synthesize such
structures or materials is their applicability in vibration isolation devices. The works of
Wang and Lakes (2004a and 2004a), which used inclusions of an elastic negative-
stiffness phase in the matrix of a viscoelastic metal to achieve extremal materials, were
found be relevant and useful for the present work. In order to study the stability of the
composite material, the authors analyzed an equivalent lumped system and found that the
lumped system can also soften with frequency for some parameters (spring stiffness,
damping coefficients in the lumped system); however, the structure was not stable for the
parameters leading to the softening. Nevertheless, the works of Wang and Lakes (2004a
and 2004b) gave inspiration to use negative stiffness inclusions. The question was how to
get stable frequency-induced softening in a composite having a negative stiffiress

component and the present work addresses that question.

The contribution of this work starts with the idea that a composite material having a
negative stiffness component can be used as a ‘smart’ vibration isolation device. Chapter
2 presents a model of a material that exhibits frequency-induced softening. The model is
2D and based on a periodic mixture of two material phases: A (matrix phase) and B
(inclusion phase). In contrast to Wang and Lakes (2004a and 2004b) where phases A and
B are viscoelastic and elastic respectively, the two phases here are elastic and viscoelastic

respectively. Later desirable results were obtained even when both the two phases are

151

 

viscoelastic. The present work thus contributes in finding alternatives to the phases used

by Lakes and his coworkers.

The two-phase composite material was approximated as a one dimensional mechanical
model in accordance with F uj ino et al. (1964) and Marinov (1978). The one-dimensional
lumped model was analyzed for frequency-induced softening and stability. The analysis
determined parameters for stable softening and forms a part of the contribution of this
work. The parameters of the 1D mechanical model that showed stable softening were
extended back to construct the corresponding 2D composite model. The stability of the
2D composite model was further reviewed as the stability criteria for the 1D mechanical
model may not be sufficient to ensure stability of the 2D composite. It was found that the
negative stiffness phase itself needed to be made stable as an inclusion before the stability
of the whole composite is investigated. Studies revealed that the simultaneous softening
and stability are achievable if phase B itself is a small scale mixture of two constituents,
one elastic (B2) and the other one with negative stiffness (B1). The idea of stabilizing the
negative stiffness phase by the rank-1 layering of the two constituents is a contribution of
this work. It should be noted, however, that while the stability of the mixture was verified
against certain modes of instability, one cannot say for sure that the material is stable
with regard to all possible modes of instability. Further analysis may determine that
additional conditions may be required. A full answer can be obtained only through

experimentation.

152

 

To realize the negative stiffness material Bl, a lumped system is proposed to be made of
negative stiffness springs, typical positive-stiffness springs and dampers. A layered
mixture of this lumped system (phase B1) and a typical elastic material (phase 32) may be
used as the inclusion phase B; however, physical realization of such inclusion phase may
be difficult. As an alternative to the layered phase B, therefore, a lumped system is built
for phase B which is stable as an inclusion. The design of lumped systems (B1 and B) and
their analysis to determine the average elastic properties and estimate stability are part of

the contribution of this work.

Chapter 3 presented examples of the composites that exhibit frequency induced softening.
Transmissibility analysis for the composites was carried out. The vibration isolation
performance for these composites (A-B mixture) was found to be better than that of the

matrix material (A) alone.

Phase B is not isotropic and therefore mixing phases A and B may result in a material
that is not isotropic either. In practice one may be interested only in isotropic materials
and therefore the micro-geometry of the mixture of A and B must be designed so that the
mixture is isotropic. This has been done in chapter 4 by casting the problem as an inverse
homogenization problem (as discussed e.g. in Sigmund (1995)). This methodology has
been used to tailor the material tensor to desired specifications, by re-adjusting the shape
of the inclusion of B inside A. The designed materials closely match the target elastic
properties. The methodology developed in that chapter differs fi'om a typical material

design problem and may be considered as a major contribution. Using positive

153

 

definiteness of the elastic tensor as a constraint to maintain stability is a novel idea
emerged from the present work. The extension of the objective function in Diaz and

Benard (2003) to suit viscoelastic materials is also new.

Once the layout of the composite material is available, the next step is to actually build
the composite. The first task in that direction would be to build the inclusion phase B. In
chapter 2, the negative stiffness (in phase B) is proposed to be realized by using bistable
structures. The realization of phase B is proposed in two steps — (i) build a two-
dimensional (tileable) bistable structure and (ii) interconnect the 2D bistable structures
using viscoelastic blocks. While chapter 5 addresses the first step, the second step is left
as a future work. The model of tileable bistable structure developed in chapter 5 is quite
simple. It ignores important constraints, such as interference, contact between members
and stress constraints. The simplifications introduced were necessary to reduce the
computational burden associated with the analysis and, in particular, with the use of
genetic algorithms. However, in spite of its simplicity, the approach results in a useful
tool to explore concepts in design of bistable, periodic structures. The main contributions
in this chapter include the idea of making a tileable bistable structure, the idea of
achieving bistability by making beams flexible at the ends (using thinner area or softer
material) and formulation of the whole optimization problem (including the design of the
objective function as well as the constraints to avoid the invalid structures). Further
studies for this part of work should include a better objective function, capable of
providing sensitivity information even for the monostable designs; the computation of

analytical sensitivity information, consistent with more efficient, gradient based

154

optimization algorithms; stress and interference constraints; and more refined analysis
models (e. g., two-dimensional elasticity instead of beam models) that may reflect the true

behavior of compliant structures more accurately.

To sum up all the chapters, it can be said that the concept presented in this work is
attractive and may provide indications on what avenues to pursue in the future. A
direction that seems perhaps more attractive involves the realization of the negative
stiffness inclusion by means of a bistable structure inside the composite. This could be
accomplished perhaps by pre-stressing one of the constituents of the composite to a post-
buckled state. Regardless, the current state of the work is only theoretical; relevant
experimental work should be pursued in future before one can suggest what the most
potentially rewarding direction to pursue should be in order to actually realize a material

with the desired properties.

155

 

REFERENCE

Al-Rasby, SN, 1991, “Solution Techniques in Nonlinear Structural Analysis,”
Computers and Structures, Vol. 40, No. 4, pp. 985-993.

Ananthasuresh, GK. and Kota, S., 1995, “Designing Compliant Mechanisms,”
Mechanical Engineering, Vol. 117, No. 11, pp. 93-96.

Bendsae, M., Ben Tal, A. and Zowe, J., 1994, “Optimization methods for truss
geometry and topology design,” Structural Optimization, Vol. 7, pp. 141-159.

Bendsee, MP. and Kikuchi, N., 1988, “Generating optimal topologies in structural
design using a homogenization method,” Computer Methods in Applied Mechanics and
Engineering 1988, Vol. 71, pp. 197-224.

Bendsee, MP. and Sigmund, 0., 2003, “Topology Optimization: Theory, Methods and
Applications,” Springer, New York.

Bensoussan, A., Lions, J.L., and Papanicolau, G., 1978, “Asymptotic Analysis for
Periodic Structures,” North-Holland, Amsterdam.

Bruns, T. E., Sigmund, O. and Tortorelli, D. A., 2002, “Numerical methods for the
topology optimization of structures that exhibit snap-through,” International Journal for
Numerical Methods in Engineering, Vol. 55, pp. 1215-1237.

Bruns, TE. and Tortorelli, DA, 2001, "Topology optimization of nonlinear elastic
structures and compliant mechanisms,” Computer Methods in Applied Mechanics and
Engineering, Vol. 190(26-27), pp. 3443—3459.

Cherkaev, AV. and Gibiansky, L.V., 1993, “Coupled estimates for the bulk and shear
moduli of a two-dimensional isotropic elastic composite,” Journal of the Mechanics and
Physics of Solids, Vol. 41 (5), pp. 937-980.

Crisfield, M.A., 1991, Non-linear Finite Element Analysis of Solids and Structures, Vol.
1, Wiley, Chichester.

Diaz, A. R. and Benard, A., 2003, “Designing materials with prescribed elastic properties
using polygonal cells,” International Journal for Numerical Methods in Engineering, Vol.

57 (3), pp. 301-314.

Dorn, W. C., Gomory, RE. and Greenberg, H.J., 1964, “Automatic design of optimal
structures,” Journal de Mécanique, Vol. 3, pp. 25-52.

Flugge, W., 1967, “Viscoelasticity,” Blaisdell Publishing Company, Waltham.

156

Frecker, M. I., Ananthasuresh, G. K., Nishiwaki, S., Kikuchi, N. and Kota, S., 1997,
“Topological synthesis of compliant mechanisms using multi-criteria optimization,”
Journal of Mechanical Design, Vol. 119(2), pp. 238-245.

Fujino, K., Ogawa, Y., and Kawai, H., 1964, “Generalized model representation relating
the degree of mixing to the rheological behavior of a mechanical mixture of two polymer
components,” Journal of Applied Polymer Science, Vol. 8(5) , pp. 2147 — 2161.

Guest, 1., Prevost, J. and Belytschko, T., 2004, “Achieving minimtun length scale in
topology optimization using nodal design variables and projection functions,”
International Journal for Numerical Methods in Engineering, Vol. 61(2), pp. 238-254.

Gibiansky, L.V. and Milton, G.W., 1993, “On the effective viscoelastic moduli of two-

phase media. I. Rigorous bounds on the complex bulk modulus,” Proc. R. Soc. London A,
Vol. 440, pp. 163—188.

Gibiansky, L.V., Milton, G.W. and Berryman, J .G., 1999, “On the effective viscoelastic
moduli of two—phase media. 111. Rigorous bounds on the complex shear modulus in two
dimensions,” Proc. R. Soc. London A, Vol. 455, pp. 2117—2149.

Guedes, J.M. and Kikuchi, N., 1990, “Preprocessing and postprocessing for materials
based on the homogenization method with adaptive finite element methods,” Computer

Methods in Applied Mechanics and Engineering, Vol. 83, pp. 143—198.

Haddad, Y.M., 1995, “Viscoelasticity of engineering materials,” Chapman & Hall,
London.

Hassani, B. and Hinton, E., 1999, Homogenization and Structural Topology Optimization,
Springer-Verlag London.

Hashin, Z. and Shtrikman, S., 1963, “A variational approach to the theory of the elastic
behaviour of multiphase materials,” Journal of the Mechanics and Physics of Solids, Vol.
11, pp. 127-140.

Herrera, I. and Gurtin, ME, 1965, “A correspondence principle for viscoelastic wave
propagation,” Quarterly of applied mathematics, Vol. 22(4), pp. 360-364.

Howell, LL, 2001, Compliant Mechanisms, John Wiley and Sons, New York, NY.

Jaglinski, T. and Lakes, R. S., 2004, “Anelastic instability in composites with negative
stiffness inclusions”, Philosophical Magazine Letters, 84, (12) 803 - 810, Dec. (2004).

Jaglinski, T., Kochmann, D., Stone, D. and Lakes, RS, 2007, “Composite Materials with
Viscoelastic Stiffness Greater Than Diamond,” Science, Vol. 315, pp. 620-622.

157

Jensen, B. D., Parkinson, M. B., Kurabayashi, K., Howell, L. L. and Baker, M. S., 2001,
“Design optimization of a fully-compliant bistable micro-mechanism,” in Proceedings of
the ASME International Mechanical Engineering Congress and Exposition,
lMECEZOOl/MEMS-23 852, New York.

Jensen, B.D. and Howell, LL, 2003, “Identification of Compliant Pseudo-Rigid-Body
Four-Link Mechanism Configurations Resulting in Bistable Behavior,” Journal of
Mechanical Design, Vol. 125, pp. 701-708.

Jensen, B.D., Howell, LL. and Salmon, LG, 1999, “Design of Two-Link, In-Plane,
Bistable Compliant Micro-Mechanisms,” ASME Journal of Mechanical Design, Vol.
121, pp. 416-423.

Joo, J.Y. and Kota, S., 2004, “Topological synthesis of compliant mechanisms using
nonlinear beam elements,” Mechanics Based Design of Structures and Machines, Vol. 32
(1), pp. 17-38.

King, C. and Campbell, M.I., 2004, “On the Design Synthesis of Multistable Equilibrium
Systems,” in Proceedings of the International Design Engineering and Technical
Conference, Salt Lake City, Utah.

Lakes, R. S. and Drugan, W. J ., 2002, “Dramatically stiffer elastic composite materials
due to a negative stiffness phase?”, J. Mechanics and Physics of Solids, 50, 979-1009.

Lakes, R. S., Lee, T., Bersie, A., and Wang, Y. C., 2001, “Extreme damping in composite
materials with negative stiffness inclusions”, Nature, 410, 565-567, 29 March (2001).

Li, J. and Weng, G.J., 1997, “A homogenization theory for the overall creep of isotropic
viscoplastic composites,” Acta Mechanica, Vol. 125, pp. 141-153.

Marinov, P., 1978, “Two-phase mixture - viscoelasticity and general equivalent
mechanical model,” Polymer Engineering and Science, Vol. 18 (4) , pp. 335 — 338.

Masters, ND. and Howell, LL, 2003, “A Self-Retracting Fully-Compliant Bistable
Micromechanism,” Journal of Microelectromechanical Systems, Transactions of
[BEE/ASME, Vol. 12, pp. 273-280.

Milton, G.W. and Berryman, J.G., 1997, “On the effective viscoelastic moduli of two-

phase media. 11. Rigorous bounds on the complex shear modulus in three dimensions,”
Proc. R. Soc. London A, Vol. 453, pp. 1849—1880.

Neves, M.M., Rodrigues, H. and Guedes, J .M., 1995, “Generalized topology design of

structures with a buckling load criterion,”, Structural and Multidisciplinary Optimization,
Vol. 10(2), pp. 71-78.

158

Prasad, J., and Diaz, AR, 2005, “Layout of Tileable Multistable Structures Using
Topology Optimization,” in Proceedings of TopoptSYMP2005, IUT AM symposium,
Rungstedgaard, Copenhagen, Denmark, October 26 — 29, 2005.

Qiu, J., Lang, J. and Slocum, A.H., 2004, “A Curved-Beam Bistable Mechanism,”
Journal of Microelectromechanical Systems, Vol. 13(2), pp. 137-146.

Sigmund 0., 1994, “Materials with prescribed constitutive parameters: an inverse
homogenization problem,” International Journal of Solids and Structures, Vol. 31(17), pp.
2313—2339.

Sigmund, O., 1995, “Tailoring materials with prescribed elastic properties,” Mechanics
of Materials, Vol. 20, pp.351-368.

Sigmund, 0., 1997, “On the Design of Compliant Mechanisms Using Topology
Optimization”, Mechanics of Structures and Machines, Vol. 25, pp. 493-524.

Sigmund, O., 2000, “A new class of extremal composites. Journal of Mechanics and
Physics of Solids,” Vol. 48, pp.397—428.

Sigmund, O., 2001, “A 99 line topology optimization code written in MATLAB,”
Structural and Multidisciplinary Optimization, Vol. 21(2), pp. 120-127.

Sigmund, O., 2007, “Morphology-based black and white filters for topology
optimization,” Structural and Multidisciplinary Optimization, Vol. 33, pp. 401-424.

Svanberg, K., 1987, “The method of moving asymptotes,” International Journal for
Numerical Methods in Engineering, Vol. 24, pp. 359 -373.

Tcherniak, D. and Thomsen, J.J., 1998, “Slow Effects of Fast Harmonic Excitation for
Elastic Structures,” Nonlinear Dynamics 17: 227—246.

Wang, Y. C. and Lakes, R. S., 2004a, "Extreme stiffness systems due to negative
stiffness elements", American .I. of Physics, 72, Jan. (2004).

Wang, Y. C. and Lakes, R. S., 2004b, "Negative stiffness induced extreme viscoelastic
mechanical properties: stability and dynamics", Philosophical Magazine, 35, 3785-3801,
Dec. (2004).

Wang, Y. C. and Lakes, R. S., 2005a, "Stability of negative stiffness viscoelastic
systems", Quarterly of Applied Math, Vol.63, pp.34-55, March (2005).

Wang, Y.C. and Lakes, R. S., 2005b, “Composites with inclusions of negative bulk

modulus: extreme damping and negative Poisson’s ratio,” Journal of Composite
Materials, Vol. 39 (18), pp. 1645- 1657.

159

Yi, Y.-M., Park, S.-H. and Youn, S.—K., 1998, “Asymptotic homogenization of

viscoelastic composites with periodic microstructures,” International Journal of Solids
and Structures, v61. 35 (17), pp. 2039-2055.

Yi, Y.-M., Park, S.-H. and Youn, S.-K., 2000, “Design of microstructures of viscoelastic

composites for optimal damping characteristics,” International Journal of Solids and
Structures, Vol. 37, pp. 4791-4810.

Yu, Y., Naganathan, N.G. and Dukkipati, R.V., 2000, “Review of automotive vehicle

engine mounting systems,” International Journal of Vehicle Design, Vol. 24 (4), pp. 299-
31 9.

160

 

1{111111111113111111,111