Smart Materials and Structures

PAPER You may also like

- Broadband absorption of infrared dielectric
Topology optimization and fabrication of low resonators for passive radiative cooling
Yanning Liu, Xiaolong Weng, Peng Zhang
frequency vibration energy harvesting et al.

- Modelling third harmonic ion cyclotron

microdevices acceleration of deuterium beams for JET
fusion product studies experiments
M. Schneider, T. Johnson, R. Dumont et
To cite this article: Jiadong Deng et al 2015 Smart Mater. Struct. 24 025005 al.

- Ultracold few fermionic atoms in needle-

shaped double wells: spin chains and
resonating spin clusters from microscopic
Hamiltonians emulated via
View the article online for updates and enhancements. antiferromagnetic Heisenberg and t–J
models
Constantine Yannouleas, Benedikt B
Brandt and Uzi Landman

This content was downloaded from IP address 201.16.190.6 on 16/04/2023 at 13:12

 Smart Materials and Structures

Smart Mater. Struct. 24 (2015) 025005 (16pp) doi:10.1088/0964-1726/24/2/025005

Topology optimization and fabrication of low

frequency vibration energy harvesting
microdevices
Jiadong Deng, Katherine Rorschach, Evan Baker, Cheng Sun and
Wei Chen
Department of Mechanical Engineering, Northwestern University, 2145 Sheridan Road, Evanston, IL, USA

E-mail: [email protected]

Received 15 September 2014, revised 13 November 2014

Accepted for publication 17 November 2014
Published 16 December 2014

Abstract
Topological design of miniaturized resonating structures capable of harvesting electrical energy
from low frequency environmental mechanical vibrations encounters a particular physical
challenge, due to the conflicting design requirements: low resonating frequency and
miniaturization. In this paper structural static stiffness to resist undesired lateral deformation is
included into the objective function, to prevent the structure from degenerating and forcing the
solution to be manufacturable. The rational approximation of material properties interpolation
scheme is introduced to deal with the problems of local vibration and instability of the low
density area induced by the design dependent body forces. Both density and level set based
topology optimization (TO) methods are investigated in their parameterization, sensitivity
analysis, and applicability for low frequency energy harvester TO problems. Continuum based
variation formulations for sensitivity analysis and the material derivative based shape sensitivity
analysis are presented for the density method and the level set method, respectively; and their
similarities and differences are highlighted. An external damper is introduced to simulate the
energy output of the resonator due to electrical damping and the Rayleigh proportional damping
is used for mechanical damping. Optimization results for different scenarios are tested to
illustrate the influences of dynamic and static loads. To demonstrate manufacturability, the
designs are built to scale using a 3D microfabrication method and assembled into vibration
energy harvester prototypes. The fabricated devices based on the optimal results from using
different TO techniques are tested and compared with the simulation results. The structures
obtained by the level set based TO method require less post-processing before fabrication and the
structures obtained by the density based TO method have resonating frequency as low as 100 Hz.
The electrical voltage response in the experiment matches the trend of the simulation data.
Keywords: topology optimization, vibration energy harvesting, microfabrication, low frequency

(Some figures may appear in colour only in the online journal)

1. Introduction from several disadvantages including expensive and harmful

disposal, large size and weight compared to more technolo-
Advances in technology have allowed the decrease in both gically-advanced energy harvesting devices, and limited
size and power consumption of complex electronic systems, lifetimes requiring periodic recharging or replacement. For
such as micro-electro-mechanical systems, implanted elec- example, pacemakers typically need battery replacements
tronic devices, wireless sensors and electronic portable devi- every few years, requiring in-patient surgery that risks health
ces. At present, batteries remain the primary power source for and uses significant hospital resources (Karami and
miniaturized electronic systems. However, batteries suffer Inman 2012).

0964-1726/15/025005+16$33.00 1 © 2015 IOP Publishing Ltd Printed in the UK

Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Figure 2. Systematic design loop.

paper will rely on electromagnetic induction energy harvest-

ing mechanism and the specific energy harvester device
concept is shown in figure 1.
In figure 1, the voltage is generated by the motion of a
magnet in and out of a coil attached to circuitry for energy
storage. The magnet is mounted on a resonator structure. The
micro resonator structure, at the millimeter size scale, is
Figure 1. Vibration energy harvester device concept. anchored to a vibrating base with harmonic oscillating dis-
placement described by ug (ω), where ω is the frequency of
An attractive alternative to batteries is to extract (harvest) the harmonic base movement. Base acceleration would cause
energy from the ambient vibration to directly power minia- relative displacement between the coil and proof mass, and
energy is extracted from the mechanical system by an elec-
turized electronic systems. The electrical energy to power
tromagnetic damping mechanism. The magnet acts as the
microelectronics can be generated from kinetic, solar, elec-
proof mass and the resonator structure is considered as the
tromagnetic, thermal or other energy sources (Paradiso and
design domain. This study will investigate the design, fabri-
Starner 2005). The obtained energy can then be used to
cation, and experimental realization of electromagnetic min-
recharge a secondary battery or, in some cases, to directly
iaturized structures capable of harvesting electrical energy
power the microelectronics. The kinetic energy or environ-
from low frequency environmental vibrations.
mental vibration energy harvesting, pervasive in the current
Realization of miniaturized electromagnetic devices for
sustainable energy research, is the focused problem of this
energy harvesting at low frequencies presents a physical
paper (Beeby et al 2006, Mitcheson et al 2008). Measure- challenge, due to the competing/conflicting requirements: low
ments of readily available vibration sources show that many resonating frequency and miniaturization. According to the
are both low amplitude (less than 1 G) and low frequency scaling law (Wautelet 2001) structural frequencies increase
(less than 150 Hz) (Baker et al 2012). For example, vibrating with the decrease of structural size of similar structures. For
parts of subway and commuter rail cars have peak accelera- micro energy harvesting devices, the space constraints require
tions between 0.2 G and 0.7 G at frequencies between 40 Hz that the power source be on the order of millimeters. A pri-
and 90 Hz. The human heartbeat produces peak accelerations mary consideration to improve the energy harvesting perfor-
inside the chest wall of about 0.3 G at 51 Hz (Kudriavtsev mance of such micro energy harvesting device is to achieve
et al 2007). Walking produces peak amplitudes in carried the low resonant frequencies that take advantage of the
accessories of about 1 G at less than 5 Hz. This paper will environmental frequency spectrum. However, due to the
target harvesting energy from low-frequency (sub-100 Hz) scaling law, as the converters are miniaturized to integrate
environmental vibrations. them on microelectronic devices, the resonance frequency
The principle behind energy harvesting from environ- increases, and it becomes much higher than characteristic
mental vibrations is to convert the mechanical energy from frequencies of many everyday mechanical stimuli. To achieve
relative displacement of a moving part or the mechanical the low resonant frequencies of micro energy harvesting
deformation of some structure inside the energy harvesting device, in this study, a low modulus polymer is chosen for
device into electrical energy. This displacement or deforma- building resonant structures to lower the stiffness and thus the
tion can be converted to electrical energy by three techniques: frequency of the structure, and structural topology optimiza-
piezoelectric, electrostatic and electromagnetic induction tion (TO) methods are explored to find the optimal structural
(Mateu and Moll 2005). Though piezoelectric converters are layout.
robust, they are difficult to miniaturize. For electrostatic TO is a mathematical approach that optimizes material
generators, proper operation of the switches or when the layout within a given design space for a given set of boundary
charges are transferred is critical for a good efficiency. This conditions, such that the resulting layout meets a prescribed

2
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

set of performance targets. TO allows complete freedom in implemented TO design problem using a variable thickness
the final topology of the solution and avails engineers of the method to endow meanings of low-density elements in the
best concept design that meets the design requirements. optimized mechanical amplifier. In their work, the targets of
Among the existing TO methods, such as homogenization first several eigenfrequencies of a micromechanical resonator
methods (Bendsøe and Kikuchi 1988, Bendsøe 1989), density are set to be above 1 GHz, so that the optimum solution of TO
based methods (Zhou and Rozvany 1991, Rozvany will not degenerate.
et al 1992, Bendsøe and Sigmund 1999), evolutionary Another design challenge is due to design dependent
structural optimization methods (Xie and Steven 1993), and body loads. As shown in figure 1, the miniaturized structure is
level set methods (LSM) (Sethian and Wiegmann 2000, mounted on the oscillating base and moves with the base. The
Osher and Santosa 2001, Wang et al 2003, Allaire et al 2004), base movement is transferred into equivalent excitation body
the density method gains large popularity because of its easy forces acting on the structure in the analysis model. Hence the
implementation and efficiency (Bendsøe and Sigmund 2003), equivalent excitation force is design dependent which com-
and LSM also plays a significant role, if combined with X-
plicates the design and optimization process.
FEM, especially in design problems where the precise
In this paper, TO of resonating micro structures for
descriptions of structural boundary or interfaces locations are
energy harvesting considering low and single frequency base
important (Luo et al 2007, Chen and Chen 2011, Sigmund
movement is investigated. A multi-objective optimization
and Maute 2013, van Dijk et al 2013). Therefore, in this
formulation including dynamic output displacement max-
study, both density and LSM based TO methods will be
imization and static lateral stiffness maximization is proposed
investigated and their pros and cons for the problem of
interest will be examined and the results will be discussed. to prevent the structure from degenerating and forcing the
Contrary to decreasing the structure frequencies to solution to be manufacturable. Base movement harmonic
prompt structural resonance for dynamic displacement max- displacements are transferred into equivalent oscillating body
imization, most existing TO works attempt to eliminate forces acting on the structure, which are design dependent
resonance by increasing the structure natural frequencies to (inertial load). Thus, contrary to existing TO designs on
minimize the dynamic displacement. For example, research- resonators (Nishiwaki et al 2000, Tcherniak 2002, Maeda
ers have considered fundamental eigenfrequency maximiza- et al 2006, He et al 2012) which consider design independent
tion (Pedersen 2000), lowest eigenfrequencies maximization load, the rational approximation of material properties
(Ma et al 1993, Du and Olhoff 2007), and dynamic dis- (RAMP) interpolation (Stolpe and Svanberg 2001) approach
placement minimization for steady-state response (Ma is employed in this work to stabilize the low-density area in
et al 1995) and for transient response (Min et al 1999). TO with dependent body load. As suggested by (Bruyneel
Optimization of those objectives is equivalent to stiffening the and Duysinx 2005), the RAMP interpolation can also sup-
structure, which guarantees that the structure should have press local vibration in dynamic problems (Pedersen 2000).
enough stiffness. However, in designing low-frequency External damper is included to simulate the energy output of
resonators, the objective is the opposite, i.e. low-frequency the resonator due to electrical damping. Both density and
resonators need to have low stiffness, which naturally raises level set based TO methods are investigated in their para-
the issue of degenerate/disappearing design solutions (low- meterization, sensitivity analysis, and applicability for low
density areas and hinges) in TO. frequency energy harvester TO problems. Detailed continuum
To prevent the design of low-frequency resonators from based variation formulations for sensitivity analysis and the
degeneration in a TO process, many methods have been material derivative based shape sensitivity analysis are pre-
proposed. For example, a multi-objective optimization pro- sented for the density method and the LSM, respectively.
blem was formulated to prevent the structure from degen-
Their similarities and differences in sensitivity analysis are
erating (Nishiwaki et al 1999, Nishiwaki et al 2000,
highlighted to allow the readers to compare and familiarize
Nishiwaki et al 2001), which is accomplished by maximizing
with the two methods. A systematic design procedure is
the mutual mean compliance in the dynamic case and
proposed that includes fabrication and physical testing for
simultaneously maximizing the dynamic stiffness to resist the
validation. The design loop, shown in figure 2, starts with
applied vibration force and the reaction force. (Maeda
et al 2006) investigated a TO method for designing vibrating structural optimization of the resonator. After post-processing
structures by targeting desired eigenfrequencies and eigen- considering manufacturing feasibility, the resonant structure
mode shapes simultaneously, as resonators resonating at the is fabricated with the projection microstereolithography
desired eigenfrequencies may not have the maximum dis- technique (Sun et al 2005). Finally, physical experiment is
placement at the desired output direction. (Tcherniak 1999) conducted to validate the achieved performance of a manu-
constrained the lower limit of static stiffness to prevent the factured design.
structure from disappearing. (Tcherniak 2002) introduced an The concept of the energy harvester considered in this
external damper at the output point to represent the energy work and its engineering analysis model are presented in
output of the system to reduce degeneration. By adding the section 2. Section 3 describes the optimization formulation
external damper, the structure has sufficient stiffness to and sensitivity calculations with our proposed method. The
transfer the excitation force to the location of the external optimization results are presented in section 4, and the fab-
damper to maximize the output displacement. (He et al 2012) rication and testing of those results are given in section 5.

3
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

structure domain Ω , the weak form of governing equation can

be obtained as finding u ∈ Z such that

∫D ⎡⎣ ε (u): D H : ε ( v̄) + iωC Hu ⋅ v̄ − ρH ω2u ⋅ v̄⎤⎦ dΩ

= ∫ ρ H ω2u0 ⋅ v̄ dΩ , ∀ v ∈ Z , (4)
Ω

where v̄ is the complex conjugate of v, which belongs to Z ,

the complex space of kinematically admissible virtual dis-
placement. Only homogeneous boundary conditions are
considered in this paper, so trial displacement u and virtual
displacement v belong to the same space Z . D H is the
effective elastic modulus tensor.
C H not only includes the mechanical damping which is
caused by internal material friction, but also includes the
electrical damping which represents the energy output of the
structure. It should be noted that the coil is also anchored on
the base, so the electrical damping force applied on the
structure is proportional to ur . The effect of the electrical
damping on the structure is considered by applying external
dampers at the output points as shown in figure 3.
From equation (4), the equivalent external load f (ω) can
be stated as:

Figure 3. Schematic analysis model. f (ω) = ∫Ω ρH ω2u0 ⋅ v̄dΩ, (5)

2. Design concept and analysis model—vibration which suggests that f (ω) now is design dependent because it
considering base movement depends on structural mass or density distribution, and f (ω)
applies on every material point of the structure. The schematic
The energy harvester concept for this study is shown in analysis model is illustrated in figure 3.
figure 1. The harmonic oscillating displacement of the base is
described by ug (ω), where ω is the frequency of the harmonic
base oscillation and ur represents the relative displacement of 3. Optimization model
the resonator to the base. For any material point on the
structure, the strong form of the governing equation can be 3.1. Optimization objective
stated as
The objective of electro-magnetic energy harvester design is
ρ H
( ür + ü g ) + C H
ur + Aur = 0, (1) to maximize the displacement of the magnet, which equals to
maximizing the velocity and damping force of the magnet
where ρ H and C H are the material physical density and under given external excitation with frequency ω. To max-
damping property respectively at the material point. Dot imize the output displacement, the frequency of the structure
denotes total time derivative. A is a linear partial differential should be tuned to be around ω in order for the structure to
operator with Aur = − div σ (ur ), where σ denotes stress resonate at low frequency ω. This kind of energy harvester
tensor and div is the divergence operator. design is similar to ‘resonator’ design in the literature
As the base excitation is assumed to be harmonic, we (Nishiwaki et al 2000, Tcherniak 2002). As mentioned ear-
have ug (ω) = u0 eiωt , where u0 is the amplitude of ug . Then, lier, the optimum topology of the resonator tends to degen-
equation (1) leads to erate, i.e. there are lots of gray areas and hinges which are
cost prohibitive or impossible to be fabricated. To prevent the
ρ H ür + C Hur + Aur = − ρ H üg = ρ H ω2u0 eiωt . (2) structure from degenerating/vanishing, multi-objective opti-
Let ur = ueiωt , where u is the amplitude of ur , mization that includes mutual mean compliance and dynamic
equation (2) is simplified as stiffness (Nishiwaki et al 2000), formulations that includes an
external viscous damper (Tcherniak 2002), and static stiffness
( A + iωC H
)
− ω2ρ H u = ρ H ω2u0 . (3) constraints (Tcherniak 1999) have been proposed by early
researchers. In the existing resonator design works (Nishiwaki
The above equation looks like a static problem, but now et al 2000, Tcherniak 2002, He et al 2012), the location of the
u is a complex quantity. The equivalent analysis model is exciting force and displacement output area are different.
illustrated in figure 3. Multiplying the above equation by the However, as can be seen from figure 3, in our problem con-
complex virtual displacement v̄ and integrating over the sidering base excitation, the excitation body force mainly

4
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

applies on the magnet (as the density of the magnet is much

larger than that of the polymer). Thus the location of the
excitation force and the output displacement points are the
same. Therefore the mutual mean compliance concept pro-
posed in reference (Nishiwaki et al 2000) is inapplicable, as is
the external damper idea proposed in paper (Tcherniak 2002).
In our problem, to suppress very weak or degenerated struc-
tures, we consider the static stiffness of the structure and
formulate a multi-objective optimization problem.
The first objective is to maximize the output displace-
ment of upper surface of the magnet under dynamic load,
which is stated as

f1 (u) = ∫S 1
u ⋅ ūdΓ , (6)

where u is the complex response obtained by solving

equation (4) and ū is the conjugate of u. S1 denotes the output
area for dynamical displacement maximization, which is
upper surface of the magnet as illustrated in figure 3. S1 is
nondesignable.
The second objective is to increase the structural lateral
stiffness to resist the undesired lateral deformation, which
would cause contact and friction in the device. This is done
by minimizing the displacement of the system under lateral
static load applied in the middle of the upper surface, which is
defined as Figure 4. Design parameterization in density based TO.

f2 ( us ) = ∫S 2
us ⋅ us dΓ , (7)
combine the two objective functions f2 and f2 as

where us is solution of equation (8), which is the weak form ( f1 )0 f2

of the governing equations of the structure under static load J=w + (1 − w ) , (9)
fs . S2 is the area for static displacement minimization, which
f1 ( f2 )0
is the middle point of the upper surface of the magnet as
illustrated in figure 3. S2 is nondesignable. The weak form of where w is the weighting coefficient for structural dynamic
the governing equations for static load case can be stated as: performance, (f1 )0 and (f2 )0 are the corresponding values of
the initial design which are used to normalize the two
∫D ε ( us ): D H : ε ( vs )dV = ∫S fS ⋅ vs dΓ , ∀ vs ∈ Y , (8) objectives. It is expected that with more weighting factors
w ∈ [0 1], more Pareto optimum designs could be obtained.
where vs is virtual displacement, belonging to virtual dis- Based on the Pareto designs, engineers can better understand
placement space Y . In the static load case, concentrated lateral the tradeoffs of the different objectives and choose a design
force applied in the middle point of the upper surface of the more rationally.
magnet is considered and body force is ignored. Other
boundary conditions are the same as in the dynamic model.
The purpose of minimizing f2 is twofold. First, it can 3.2. Density based TO formulation and sensitivity analysis
improve the structural static stiffness greatly reducing the
probability the structure will vanish. Secondly, it minimizes 3.2.1. Density based TO formulation. For density based TO
the lateral displacement, which is undesirable and should be method, a nonshape design variable—pseudo density field
suppressed for a stable energy harvester. Thus our goal is to function ρ —is introduced to parameterize the material
design a structure that is flexible in the dynamic sense how- distribution in the design domain. As shown in figure 4, D
ever maintains a specified degree of stiffness in the static represents design domain and x ∈ D denotes a point in D .
sense. Let ρ (x) ∈ [0,1] denote pseudo density field function at point
To construct a single aggregate objective function J and x , which describes the material distribution in design domain
find the Pareto optima of the multi-objective optimization D . Then based on the RAMP interpolation scheme, the
problem, the intuitive weighting method is employed to material effective Hooke tensor D H (x) and physical density

5
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

ρ H (x) at point x can be respectively expressed as, following Lagrange equation.

ρ (x ) J = J + W + U. (14)
D H (x ) = D (x )
1 + R (1 − ρ (x )) Taking the first order variation of equation (14), we have
ρ H (x) = ρ (x), (10) the following:

where D (x) is Hooke tensor of the base material and J′ = J′ + W ′ + U′

parameter R in the RAMP interpolation is set to be 4. If = Ju u + Ju s u
Rayleigh damping is assumed, the material damping C H (x)
+ a ( u , v̄ ) + iωc ( u , v̄ ) − ω2m ( u , v̄ ) + a′ ( u, v̄ )
can be calculated by material stiffness and density.
For the energy harvester design problem considered, the + iωc′ ( u, v̄ ) − ω2m′ ( u, v̄ ) − l′ ( v̄ )
multi-objective optimization formulation using the density + A ( u s, vs ) + A′ ( us , vs )
based TO is stated as
= a′ ( u, v̄ ) + iωc′ ( u, v̄ ) − ω2m′ ( u, v̄ )
find ρ (x ) − l′ ( v̄ ) + A′ ( us , vs ), (15)
min J ( u (ρ), us (ρ) )
s.t. W (ρ , u) = a ( u, v̄ ) + iωc ( u, v̄ ) where the adjoint dynamic complex response v̄s and adjoint
static response vs , both satisfying homogeneous boundary
− ω2m ( u, v̄ ) − l ( v̄ ) = 0 conditions, are determined by the following adjoint equations
U ( ρ , us ) = A ( us , v̄s ) − L ( v̄s ) = 0 respectively:

∫D ρdV ⩽ V¯ a ( u , v̄ ) + iωc ( u , v̄ ) − ω2m ( u , v̄ ) = −Ju u , ∀ u ∈ Z

ρmin ⩽ ρ ⩽ 1, (11) A ( u s, vs ) = −Jus u s, ∀ u s ∈ Y . (16)

In equation (15), the prime on the operators represents
where ρmin = 0.001 and the following definitions are used to the explicit dependence terms on the pseudo density field ρ .
simplify notations for the integrands: As can be seen, the adjoint method is used to cancel the
variations u and u s which cannot be calculated directly
a ( u, v̄ ) = ∫D ε (u): D H : ε ( v̄) dΩ because they are implicitly dependent on ρ through the
governing equations. Finally the real part is extracted and the
c ( u, v̄ ) = ∫ C H u ⋅ v̄ dΩ sensitivity is stated as,
D

m ( u, v̄ ) = ∫D ρH ω2u ⋅ v̄dΩ {
J ′ = Re a′ ( u, v̄ ) + iωc′ ( u, v̄ ) − ω2m′ ( u, v̄ )
− l′ ( v̄ ) + A′ ( us , vs ) }, (17)
l ( v̄ ) = ∫ ρ H ω2u0 ⋅ v̄ dΩ
D
where Re means only the real part is considered.
A ( us , vs ) = ∫ ε ( us ): D H : ε ( vs )dV
D
3.3. LSM based TO formulation and shape sensitivity analysis
L ( vs ) = ∫ fS ⋅ vs dΓ . (12)
S 3.3.1. LSM based TO formulation. In LSM based TO, shape
design variable Ω , which is the region occupied by the
3.2.2. Continuum sensitivity analysis via variation
structure, is introduced to parameterize the design space. Ω is
calculation. The continuum based sensitivity analysis
included in a fixed design domain D , i.e., Ω ⊆ D. The
approach is based on the first order variation in the calculus boundary of the structure, ∂Ω , is implicitly embedded as the
of variations. The prime symbol will be used to denote the zero level set of a one-higher dimensional surface ϕ, which is
first variation. For example, let ρ denote a design field called the level set function, as is shown in figure 5. For any
function, F (ρ) be a function that depends on the current point x ∈ D , if ϕ (x) > 0 means x ∈ (Ω \∂Ω ); then ϕ (x) = 0
design ρ and assume that F (ρ) is continuous with respect to means x ∈ (D \∂Ω ) and ϕ (x) < 0 means x ∈ (D \Ω ). The
design ρ . If the current design is perturbed in the direction of greatest advantage of the implicit representation lies in its
ability to naturally handle topological changes, such as
ρ′ (arbitrary), then the variation of F (ρ) in the direction of ρ′
is defined as splitting and merging of the boundary. Calculating the
material time derivative of the equation ϕ (x) = 0 yields the
d ∂F following Hamilton–Jacobi equation:
F′ ≡ F (ρ + τρ′) = ρ′ , (13)
dτ τ= 0 ∂ρ ∂ϕ
+ Vn ϕ = 0, (18)
where τ is a small parameter that controls the ∂t
perturbation size. where Vn denotes the normal directional design velocity field,
For a given design ρ and equilibrium solutions u and us , that drives changes to the level set function in the Hamilton–
we have W (ρ , u) = 0 , and U (ρ , us ) = 0 , which leads to the Jacobi equation during the optimization process. As it can be

6
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

3.3.2. Shape sensitivity analysis via material derivative

approach. To minimize the objective functional, the
change of the objective functional with respect to a small
variation of the shape design variable Ω needs to be
quantified to provide necessary information for updating the
current design. This process is called shape sensitivity
analysis and the result is called shape derivative
(Sokolowski and Zolesio 1992, Choi and Kim 2005). The
shape derivative will be utilized to construct a design velocity
filed Vn using the steepest decent like method, which will be
used to evolve the boundary by solving the level set function
(Sethian 1999, Osher and Fedkiw 2003). The first work
(Sethian and Wiegmann 2000) that explored LSM base TO
method in structural design constructed the velocity filed Vn
by an evolutionary approach, which removed material in
regions of low stress and added material in regions of high
stress. Later, variation approach (Osher and Santosa 2001,
Wang et al 2003) was used to establish the relation between
the design velocity field Vn and the optimization objective
functional. Now, shape design sensitivity analysis
(Sokolowski and Zolesio 1992, Allaire et al 2004, Choi and
Kim 2005), especially the material derivative approach, is
Figure 5. Design parameterization in LSM based TO.
more popular as it is more convenient. In this paper the
material derivative approach will be employed.
The material derivative has been well studied in
seen, setting up the relationship between Vn and the continuum mechanics, and has many names such as
optimization problem plays a critical role in LSM based TO convective derivative, total derivative, etc (Belytschko
method. et al 1999). The optimization objective performances can be
For the energy harvester design problem considered in formulated either directly as a domain or boundary integral or
the paper, the multi-objective optimization formulation using indirectly as a function of the above integrals. For domain-
LSM based TO method is stated as: based objective functional F (Ω , u) = ∫Ω f (u (Ω ))dΩ , its
material derivative can be stated as:
find Ω
min J ( u (Ω), us (Ω) )
Ḟ = ∫∂Ω f (u) Vn dΩ+ ∫Ω f ′dΩ, (21)
s.t. W (Ω , u) = a ( u, v̄ ) + iωc ( u, v̄ )
− ω2m ( u, v̄ ) − l ( v̄ ) = 0 where super dot means material/total derivative. For bound-
U ( Ω , us ) = A ( us , v̄s ) − L ( v̄s ) = 0 ary-based objective functional F (Ω , u) = ∫∂Ω f (u (Ω ))dΓ ,
its material derivative can be stated as:
∫Ω dV ⩽ V¯ , (19)
⎛ ⎞
where the following definitions are used to simplify notations
Ḟ = ∫∂Ω ⎝ ∂∂nf + fκ⎠ Vn dΓ + ∫∂Ω fu u dΓ ,
⎜ ⎟ (22)
for the integrands:
where ∂f ∂n = f ⋅ n, and κ = div n is the curvature of a
a ( u, v̄ ) = ∫Ω ε (u): D: ε ( v̄) dΩ point on the boundary, and the prime denotes the partial
derivative of a physical quantity with respect to the pseudo
c ( u, v̄ ) = ∫ Cu ⋅ v̄ dΩ time. The terms in equations (21) and (22) involving velocity
Ω
are often referred to as convective derivative in continuum
m ( u, v̄ ) = ∫Ω ρω2u ⋅ v̄dΩ mechanics (Belytschko et al 1999).
For a given design Ω and equilibrium solutions u and us ,
l ( v̄ ) = ∫ ρω2u0 ⋅ v̄ dΩ we have W (Ω , u) = 0 , and U (Ω , us ) = 0 , which leads to the
Ω
following Lagrange equation:
A ( us , vs ) = ∫ ε ( us ): D : ε ( vs )dV
Ω
J = J + W + U. (23)
L ( vs ) = ∫ fS ⋅ vs dΓ . (20)
S
Let the prime symbol denote the first variation in the
calculus of variations, and take material derivative of

7
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

equation (23), the following equation can be obtained. calculation. However, by comparing the sensitivity analysis
of the two methods, we can find that the derivation procedures
J ̇ = J ̇ + Ẇ + U̇ are quite similar, and especially the adjoint problems
= Ju u + Jus us′ (equations (16) and (27)) are exactly the same. The difference
+ a ( u , v̄ ) + iωc ( u , v̄ ) − ω2m ( u , v̄ ) + a ( u, v ′) is that for density based TO method, the sensitivity of the
whole field is obtained (see equation (17)). However, for
+ iωc ( u, v ′) − ω2m ( u, v ′) − l ( v ′)
LSM based TO method only the sensitivity on the boundary
+ ∫∂Ω [ ε (u): D: ε ( v ) + iωCu ⋅ v of the structure is used (see equation (28)). This would
usually lead to a slow convergence rate of the LSM method
− ρω2u ⋅ v − ρω2u0 ⋅ v ⎤⎦ Vn dΓ compared with the density method (Sigmund and
+ A ( us′, vs ) + A ( us , vs′ ) − L ( vs′ ) Maute 2013), though more efficient extended LSM (Wang
et al 2007) may be used to alleviate this problem.
+ ∫∂Ω ε ( us ): D: ε ( vs ) Vn dΓ. (24)

As a (u, v ′) + iωc (u, v ′) − ω2m (u, v ′) = l (v ′) and

A (us , vs′ ) = L (vs′ ), equation (24) can be further reduced to 4. Numerical examples
J ̇ = J ̇ + Ẇ + U̇ Figure 6 illustrates the geometrical model and boundary
= Ju u + Jus u s
conditions of the plane strain numerical example considered
+ a ( u , v ) + iωc ( u , v ) − ω2m ( u , v ) in this paper. The dimension of the design domain is
3.69 mm by 5.82 mm. The magnet is the nondesignable
+ ∫∂Ω ⎡⎣ ε (u): D: ε ( v ) + iωCu ⋅ v − ρω2u ⋅ v domain with dimension 3.69 mm by 0.43 mm, and total mass
− ρω2u0 ⋅ v ⎤⎦ Vn dΓ
1.144 g. The whole structure is meshed into 130 by 220
quadrilateral bilinear finite elements. As discussed in
+ A ( u s, vs ) + ∫∂Ω ε ( us ): D: ε ( vs ) Vn dΓ. (25) section 2, to avoid degeneration in TO, two load cases,
dynamic and static, are considered. For the dynamic load
Finally, we have case, we transfer the vertical harmonic base excitation to the
equivalent design dependent harmonic vibration force acting
J̇= ∫∂Ω ⎡⎣ ε (u): D: ε ( v ) + iωCu ⋅ v − ρω2u ⋅ v on every point of the structure, as shown in figure 6. The
frequency and magnitude of the base excitation ω is 100 Hz
− ρω2u0 ⋅ v + ε ( us ): D : ε ( vs ) ⎤⎦ Vn dΓ , (26) and u0 = 1 mm. External/electrical damping is applied to
the upper surface of the structure, with damping coefficient
where the adjoint dynamic complex response v̄ and adjoint being 0.0139 kg s−1. Mechanical damping is considered by
static response v̄s , satisfying homogeneous boundary condi- prescribing Rayleigh damping coefficients of the material,
tions, are determined by the following adjoint equations with mass proportional Rayleigh damping coefficient equal
respectively:
to 1 × 10−4 and stiffness proportional Rayleigh damping
a ( u , v ) + iωc ( u , v ) − ω2m ( u , v ) = − Ju u , ∀ u ∈ Z coefficient equal to 2 × 10−3. S1, which denotes the output
A ( u s, vs ) = − Ju s u s, ∀ u s ∈ Y . (27) area for dynamical displacement maximization, is the upper
surface of the magnet. For the static load case, a con-
As can be seen, the adjoint method is used to cancel the centrated lateral force is applied at the middle point of the
variations u and u s which cannot be calculated directly upper surface of the magnet. S2 , which denotes the place for
because they are implicitly dependent on Ω through the static displacement minimization, is the middle point of the
governing equations. It can be clearly identified that for the upper surface of the magnet.
objective J to decrease, a feasible velocity field is: In order to lower the device natural frequency into the
range expected from common vibration sources, the resonator
{
Vn = −Re ε (u): D : ε ( v ) + iωCu ⋅ v − ρω2u ⋅ v
material is chosen to be 1, 6-hexanediol diacrylate (HDDA,
− ρω2u0 ⋅ v + ε ( us ): D : ε ( vs ) , } (28) Sigma-Aldrich), a polymer with significantly lower elastic
modulus (1200 MPa fully cured) than metals and ceramics.
where Re means only the real part is considered. The Young’s modulus, Poisson’s ratio and density of the
material used are 530 MPa, 0.3 and 1011 kg m−3, respec-
3.4. Features of density based versus LSM based TO
tively. In this research, the resonator designs are fabricated
In density based TO method, nonshape design variable— with projection microstereolithography (Sun et al 2005), a 3D
pseudo density field is introduced and variation based con- fabrication method with a lateral resolution of 7.1 microns
tinuum approach is utilized for sensitivity analysis. In LSM and a vertical resolution of 20 microns. The same method was
based TO method, a shape design variable is introduced and used to fabricate vibration energy harvesters on the order of
material derivative approach is used for shape derivative millimeters in our prior research (Baker et al 2012).

8
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Figure 6. Geometrical model and boundary conditions of the dynamic load case (left) and static load case (right) for the numerical example.

4.1. Density based TO results and 1 × 10−4 , respectively. Figure 8 shows the results with
different weighting coefficients w . When w = 0 , a com-
The density based TO optimization formulation in
pliance minimization solution for a cantilever is obtained.
equation (9) is solved with several different settings of
When w becomes larger, which indicates that the dynamic
weighting factors w to investigate the effect of the dynamic
force and the static force. For solution search, the method of displacement maximization performance is more emphasized,
moving asymptotes optimization solver (Svanberg 2001) is the optimal structure becomes more compliant/flexible. The
employed, and the volume preserving filter (Xu et al 2010) is critical structure member becomes thinner when w = 0.5
used to suppress numerical instability issues (Sigmund and compared with that of w = 0.25, and finally when w = 0.75,
Petersson 1998) associated with the density based TO more compliant structures are formed. In summary, with the
method. The initial values for (f1 )0 and (f2 )0 are set at 3e-8 and increase of w , the dynamic output displacement measure f1
4e-4, respectively. The material volume fraction is 36%. and the static displacement measure f2 increase, while the
Table 1 shows the results with different weighting factors structural frequency decreases, because the structure has
w , using four-node quadrilateral finite element in analysis. It become more flexible.
is noted that the dynamic displacement objective f1 improves It should be noted that the purpose of providing the
a lot with the increase of w . However this improvement is results from two different TO methods is not to conclude
mainly attributed to the one node or one element connected which method is better. Both methods are sensitive to the
hinges on the upper body of the structure. It is these hinges initial starting points which may lead to local minima and
that contribute to the flexibility of the structures. After fab- there are many optimization parameters to be tuned. In
rication, the hinges appearing in the optimal structure have to addition, the initial values to normalize the two objectives in
be post-processed. As the dynamic performance is dependent the density based method and the LSM are not the same. Thus
on these hinge-like areas, the required post-processing would a fair comparison between the two methods is hard to obtain.
deteriorate the dynamic performance of the optimal design. However, by comparing the trend of data in tables 1 and 2, for
In the meantime, the structural static displacement mea- example the frequency of the structure, we can conclude that,
sure f2 increases and structural fundamental frequency for our specific problem, the density based method has been
decreases with the increase of w , as the structure becomes
able to produce more compliant designs. However this can be
more flexible. It can be seen that when w = 0.75, the structural
mainly attributed to the existence of one node or one element
natural frequency is 237 Hz, which is still more than two
connected hinges and gray element in the final optimal
times of the frequency 100 Hz of the excitation force.
structures. Considering fabrication feasibility, the hinges and
gray elements have to be removed. As the dynamic perfor-
4.2. LSM Based TO results
mance is dependent on these hinge-like areas, removing them
For the LSM based TO optimization formulation in may deteriorate the dynamic performance of the optimal
equation (19), the initial values for (f1 )0 and (f2 )0 are 2 × 10−9 design from density based TO design.

9
 Smart Mater. Struct. 24 (2015) 025005
Table 1. Density based TO results.

w 0.0 0.25 0.5 0.75

Structural topology
10

f1 5.7e-10 9.3e-10 1.0e-9 1.1e-9

f2 3.8e-4 2.4e-2 4.2e-2 5.7e-2
Structural frequency (Hz) 761 374 291 237

J Deng et al
 Smart Mater. Struct. 24 (2015) 025005
Table 2. LSM Based TO results.

w 0.0 0.25 0.5 0.75

Structural topology
11

f1 7.4e-10 8.9e-10 1.1e-9 1.4e-9

f2 3.6e-5 3.1e-4 6.2e-4 1.4e-3
Structural frequency (Hz) 813 778 740 604

J Deng et al
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Figure 7. Post-processed structures for density based TO results with (a)w = 0.0, (b) =0.75 and LSM based TO results (c)w = 0.0, (d)w = 0.75.

4.3. Post-processing of the results prepared, the CAD structure is sliced using a MATLAB code
developed specifically for this system. The UV absorber and
For post-processing, each of the four two-dimensional solu- light intensity concentration is tuned to obtain a curing depth
tions are imported and traced in SolidWorks sketches. The of 20 microns, determining the necessary slicing.
hinge-like areas are thickened iteratively until the structures
The silicon wafer is then aligned with the top of the
can be manufactured. The necessary minimum thickness is
liquid resin layer, and the 160 liter PuSL chamber is filled
around 0.1 mm to remove the hinges and small members. The
with nitrogen gas. This reduces the concentration of oxygen
two-dimensional post-processed designs are shown in
within the chamber and ensures optimal solidification and
figure 7.
resolution of the photo-curable resin. Afterwards, the layer
It is worth mentioning that, for TO design of 2D struc-
building process begins. The first sliced bitmap image is
tures using LSM, a quadratic energy functional can be displayed on the LCD dynamic mask (in this case, a
incorporated in the objective function using a penalty para- 1400 × 1050 pixel array), and the wafer drops by 20 microns.
meter to control the geometrical width of the structure (Chen
The system then waits for 30 s for the liquid resin to settle.
et al 2008) and to suppress hinge formation (Luo et al 2008).
The UV lamp is turned on for 16 s, reflects off a beam
Thus for 2D problems it is not necessary to include post-
splitting mirror, passes through a reduction lens and finally
processing to facilitate manufacturing. However, proper tun-
projects onto the surface of the resin in high resolution, with
ing of the penalty parameter is not straightforward (Chen
each pixel corresponding to 7.1 × 7.1 μm2. This process then
et al 2008) and the quadratic energy functional is not
repeats for each bitmap layer in the fabrication. The micro-
applicable to 3D structures. structure is then removed from the PμSL machine, cleaned
with IPA, dried under a low flow rate nitrogen gun. At this
point, the resin within the structure has not completely soli-
5. Experimental validation dified. To finish the curing process and bring the resin to its
final state, the structure is exposed to a UV post-cure.
5.1. Micro-fabrication with projection microstereolithography The two fabricated resonator structures based on the
optimal solutions from the density and the LSM methods
The sketches in figure 7 are then extruded to 3.69 mm in the
respectively are shown in figure 9. Some minor differences
third dimension and the thick base plates are included for ease
are evident between the CAD designs and the fabricated
of fabrication. The resonator structures are fabricated with
structures. In particular, some of the smaller beams inside of
projection microstereolithography (PμSL). PμSL builds
the density based TO designs (left) were near the limit of the
microstructures from a photo-curable resin in a layer-by-layer
PuSL resolution. These beams are fabricated slightly thicker
fashion directly from a 3D CAD design. Each layer is cured in
than designed to ensure structural integrity.
a single exposure by using a liquid crystal display (LCD)
panel as a dynamic mask for the UV light. This allows for a
drastic reduction in fabrication time compared with micro-
5.2. Testing
stereolithography (μSL), which fabricates 3D structures in a
point-by-point scanning fashion. The experimental setup, shown in figure 10, is designed to
The entire process flow is shown in figure 8. Prior to simulate environmental vibration sources. To represent the
fabrication, a photo-curable resin is mixed. For this paper, the general waveform of the vibration, a control signal is selected
resin is a combination by mass of 97.78% 1,6-Hexanediol and generated by a function generator (Agilent 33120A). The
Diacrylate (HDDA, Sigma-Aldrich) as the low viscosity control signal is then amplified by a power amplifier (Bruel
monomer, 2% photoinitiator (Irgacure 819, Ciba), and 0.22% and Kjaer no. 2718) before entering into an electromagnetic
UV absorber (Sudan I, Sigma-Aldrich). When the resin is shaker (LDS V203). The shaker’s waveform and acceleration

12
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Figure 8. Projection microstereolithography process flow.

are monitored by a LabVIEW program using a uniaxial

accelerometer (PCB Piezotronics no. 333B50) and Force
Sensor (PCB Piezotronics no. 201C01) mounted to the shaker
drive output. The displacement of the electromagnetic system
was monitored using the Keyence LK-G32 laser displacement
sensor. An aluminum stand is used to position the device
under testing (DUT) sufficiently away from the magnetic field
of the shaker, so as to not to cause magnetic interference. The
base of the DUT is fixed to the top of the aluminum stand
with adhesive. The device itself consists of three main com-
ponents: the resonator structure, a 1.143-gram magnetic load,
and a 6 mm inner diameter, 8 mm outer diameter and 3 mm
height coil (48 AWG magnet wire) with a resistance of 590 Ω.
The magnetic load is comprised of three layers: a 0.5-gram
tungsten mass sandwiched between two 0.321-gram NdFeB
(N50) rare earth magnets. The layers are fixed to each other
and the top of the spring structure with epoxy. The magnetic
Figure 9. Fabricated resonators for (left) density based method w coil is positioned such that the motion of the magnet is within
= 0.75, (right) level set method w = 0.75. the coil. Lastly, the voltage output of the device is measured
across a resistive load and recorded with LabVIEW.

13
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Figure 10. (Left) experimental setup, (right) neodymium-ion-boron rare earth magnet(s) on top of the projection micro stereo-lithography
fabricated HDDA polymer structure (lower right, VEH with load mass, top right: inside wire coil).

Separate vibration sweeps are performed to extract the 100 Hz, while the LSM optimized structure LSM1 with a
electrical performance of the device. These sweeps represent a 1.143 g magnet has a major peak at around 510 Hz. D2 and
forced, harmonic vibration, with a sinusoidal input over a LSM2 have more added masses to prove that the peaks
45–550 Hz frequency sweep. The signal is controlled using a obtained are reliable.
custom LabVIEW program with proportional gain feedback The trend of the experiment data matches well that of the
control to maintain a constant acceleration of 0.6 G, so as to optimized results. From the simulation results in table 1 and
decouple the test stand dynamics from the device dynamics. table 2, we can see that density based optimized structure with
The voltage output across varying load resistances is then w = 0.75 has natural frequency 237 Hz, and the LSM struc-
collected in the LabVIEW program at a sampling rate of ture 604 Hz, though the frequency of the excitation force is
10 kHz. At each frequency value 20 wavelengths of data are 100 Hz. For the experiment data, the density based structure
collected, the first 10 wavelengths are ignored and peak to resonates at around 100 Hz, and LSM around 500 Hz. The
peak voltage values of the second 10 wavelengths are aver- trend of experiment data matches the trend of optimized
aged and logged as the final value. To maintain a constant results. The mismatch between the simulation data and the
peak-to-peak acceleration, a similar process was used, if the experiment data is attributed to differences in the material
acceleration was 0.6 G +/−5 % then the voltage reading from properties (Young’s modulus and damping coefficients),
physical models, and electrical damping, especially for the
the accelerometer was stored, if the acceleration was outside
material properties, many factors matter in the manufacturing
of that range, the input to the shaker system was adjusted
process, such as the UV light intensity, exposure time, curing
according to a control theory prediction and subsequent tests
time, and volume fraction of different polymers. For the
were performed until the acceleration was in the desired
damping properties of the material, there are no accurate
range. After implementing control theory the system takes
methods to measure them. How to accurately measure and
between one and six tests to reach the acceleration goal. From
predict the material properties and integrate them in the
these tests, the natural frequency and power output for the design process would need future endeavor.
vibration energy harvester are found experimentally.
The energy harvester prototypes containing the optimized
designs showed clear resonances at 135.5 Hz and 171.5 Hz
for the density based designs and level set designs respec- 6. Conclusion
tively with w = 0.75. The structures were tested multiple
times, at constant accelerations ranging from 0.2 G to 0.8 G, In this paper, a TO based formulation is proposed to max-
which is representative for the vibration sources identified, for imize the frequency response in the design of vibration energy
frequency sweeps from 45 Hz to 550 Hz. Figure 11 shows the harvesting devices. Degenerate, nonmanufacturable optimal
voltage responses of the two structures designed by the structures are suppressed by introducing a secondary, static
density based TO optimization method (named D1 and D2) load case that represents undesired lateral forces in the
and two structures designed using the LSM optimization environment. Simultaneously minimizing the displacement
method (named LSM1 and LSM2) with w = 0.75 and tested response to such a force ensures a connected structure that
at constant accelerations. The density based optimized behaves like a dynamic resonator in the desired direction of
structure D1 with a 1.143 g magnet has a peak at around excitation.

14
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

1.4 References
1.2
D1 with 1.143g Allaire G, Jouve F and Toader A M 2004 Structural optimization
1
D2 with 2.143g using sensitivity analysis and a level-set method J. Comput.
LSM1 with 1.143g Phys. 194 363–93
Pack-peak voltage (V)

0.8 LSM2 with 9.143g Baker E, Reissman T, Zhou F, Wang C, Lynch K and Sun C 2012
Microstereolithography of three-dimensional polymeric
0.6 springs for vibration energy harvesting Smart Mater. Res.
2012 9
0.4 Beeby S P, Tudor M J and White N M 2006 Energy harvesting
vibration sources for microsystems applications Meas. Sci.
0.2 Technol. 17 R175
Belytschko T, Moran B and Liu W K 1999 Nonlinear Finite Element
0 Analysis for Continua and Structures (New York: Wiley)
0 100 200 300 400 500 600
Bendsøe M P 1989 Optimal shape design as a material distribution
Frequency (Hz)
problem Struct. Optim. 1 193–202
Figure 11. Experimental measurements of frequency dependent Bendsøe M P and Kikuchi N 1988 Generating optimal topologies in
voltage output for density based design (w = 0.75) and level set structural design using a homogenization method Comput.
design (w = 0.75) with different added masses. Methods Appl. Mech. Eng. 71 197–224
Bendsøe M P and Sigmund O 1999 Material interpolation schemes
in topology optimization Arch. Appl. Mech. 69 635–54
The success of this strategy is demonstrated system- Bendsøe M P and Sigmund O 2003 Topology optimization: theory
atically by device simulation, 3D micro-fabrication, and Methods and Applications (Berlin: Springer)
Bruyneel M and Duysinx P 2005 Note on topology optimization of
vibration testing of resultant designs. The fabricated and continuum structures including self-weight Struct. Multidisc.
tested resonators show that this optimization formulation, Optim. 29 245–56
used in both the density based method and LSM, can produce Chen S, Wang M Y and Liu A Q 2008 Shape feature control in
structural topology optimization Comput. Aided Des. 40
manufacturable resonators with improved (i.e. lower resonant
951–62
frequency) performance. The testing showed structures with Chen S K and Chen W 2011 A new level-set based approach to
100 Hz and 510 Hz resonating frequencies are obtained using shape and topology optimization under geometric uncertainty
density based TO and LSM based TO methods. Further Struct. Multidisc. Optim. 44 1–18
Choi K K and Kim N H 2005 Structural Sensitivity Analysis and
improvements could be made by weighting the dynamic Optimization 1 Linear System (Berlin: Springer)
response more heavily in the optimization formulation, or by Du J B and Olhoff N 2007 Topological design of freely vibrating
using a 3D design formulation. However, new design con- continuum structures for maximum values of simple and
cepts need to be explored if significant reduction of frequency multiple eigenfrequencies and frequency gaps Struct.
Multidisc. Optim. 34 91–110
is desired. He W, Bindel D and Govindjee S 2012 Topology optimization in
With this formulation, we explore and compare the use of micromechanical resonator design Optim. Eng. 13 271–92
two popular TO methods, density based and level set based, Karami M A and Inman D J 2012 Powering pacemakers from
for the design of low-frequency mechanical resonators. In our heartbeat vibrations using linear and nonlinear energy
harvesters Appl. Phys. Lett. 100 042901
specific numerical example, the density based method resul- Kudriavtsev V, Polyshchuk V and Roy D L 2007 Heart energy
ted in lower-frequency resonators when fabricated. The level signature spectrogram for cardiovascular diagnosis Biomed.
set results had fewer areas of intermediate density elements Eng. Online 6
and needed less post-processing. Luo J Z, Luo Z, Chen S K, Tong L Y and Wang M Y 2008 A new
level set method for systematic design of hinge-free compliant
Challenges remain for both density based and LSM based mechanisms Comput. Methods Appl. Mech. Eng. 198 318–31
TO methods. Post-processing is required to widen these Luo Z, Tong L Y, Wang M Y and Wang S Y 2007 Shape and
regions to ensure structural integrity after manufacturing. topology optimization of compliant mechanisms using a
Future extensions to this work include the addition of geo- parameterization level set method J. Comput. Phys. 227
680–705
metrical constraints to improve manufacturability, and the Ma Z-D, Kikuchi N and Cheng H-C 1995 Topological design for
consideration of the nonlinear frequency-dependent dynamic vibrating structures Comput. Methods Appl. Mech. Eng. 121
behavior of the polymer material in the finite element 259–80
formulations. Ma Z-D, Kikuchi N and Hagiwara I 1993 Structural topology and
shape optimization for a frequency response problem Comput.
Mechanics 13 157–74
Maeda Y, Nishiwaki S, Izui K, Yoshimura M, Matsui K and
Terada K 2006 Structural topology optimization of vibrating
Acknowledgments structures with specified eigenfrequencies and eigenmode
shapes Int. J. Numer. Methods Eng. 67 597–628
Mateu L and Moll F 2005 Review of energy harvesting techniques
Support from the National Science Foundation for the and applications for microelectronics Proc. SPIE 5387 359
Graduate Research Fellowship Grant No. DGE-0824162 and Min S, Kikuchi N, Park Y C, Kim S and Chang S 1999 Optimal
research grants No. CMMI-1130948 and CMMI-0955195 is topology design of structures under dynamic loads Struct.
greatly appreciated. Optim. 17 208–18

15
Smart Mater. Struct. 24 (2015) 025005 J Deng et al

Mitcheson P D, Yeatman E M, Rao G K, Holmes A S and Green T C Sokolowski J and Zolesio J-P 1992 Introduction to shape
2008 Energy harvesting from human and machine motion for optimization Introduction to Shape Optimization 16 (Berlin:
wireless electronic devices Proc. of the IEEE 96 1457–86 Springer) pp 5–12
Nishiwaki S, Min S, Yoo J and Kikuchi N 2001 Optimal structural Stolpe M and Svanberg K 2001 An alternative interpolation scheme
design considering flexibility Comput. Methods Appl. Mech. for minimum compliance topology optimization Struct.
Eng. 190 4457–504 Multidisc. Optim. 22 116–24
Nishiwaki S, Saitou K, Min S and Kikuchi N 2000 Topological Sun C, Fang N, Wu D M and Zhang X 2005 Projection micro-
design considering flexibility under periodic loads Struct. stereolithography using digital micro-mirror dynamic mask
Multidisc. Optim. 19 4–16 Sensors Actuators A 121 113–20
Nishiwaki S, Silva E, Saitou K and Kikuchi N 1999 Topology Svanberg K 2001 A class of globally convergent optimization
optimization of actuators using structural flexibility Proc. of methods based on conservative convex separable
3rd WCSMO’99 pp 4–6 approximations SIAM J. Optim. 12 555–73
Osher S and Fedkiw R 2003 Level Set Methods and Dynamic Tcherniak D 1999 Topology optimization of resonating actuators
Implicit Surfaces (Berlin: Springer) Proc. of the 12th Nordic Seminar on Computational
Osher S J and Santosa F 2001 Level set methods for optimization Mechanics pp 73–6
problems involving geometry and constraints I. Frequencies of Tcherniak D 2002 Topology optimization of resonating structures
a two-density inhomogeneous drum J. Comput. Phys. 171 using SIMP method Int. J. Numer. Methods Eng. 54
272–88 1605–22
Paradiso J A and Starner T 2005 Energy scavenging for mobile and van Dijk N P, Maute K, Langelaar M and van Keulen F 2013 Level-
wireless electronics IEEE Pervasive Comput. 4 18–27 set methods for structural topology optimization: a review
Pedersen N L 2000 Maximization of eigenvalues using topology Struct. Multidisc. Optim. 48 437–72
optimization Struct. Multidisc. Optim. 20 2–11 Wang M Y, Wang X M and Guo D M 2003 A level set method for
Rozvany G, Zhou M and Birker T 1992 Generalized shape structural topology optimization Comput. Methods Appl. Mech.
optimization without homogenization Struct. Optim. 4 250–2 Eng. 192 227–46
Sethian J A 1999 Level Set Methods and Fast Marching Methods: Wang S Y, Lim K M, Khoo B C and Wang M Y 2007 An extended
Evolving Interfaces in Computational Geometry, Fluid level set method for shape and topology optimization
Mechanics, Computer Vision, and Materials Science J. Comput. Phys. 221 395–421
(Cambridge: Cambridge University Press) Wautelet M 2001 Scaling laws in the macro-, micro-and nanoworlds
Sethian J A and Wiegmann A 2000 Structural boundary design via Eur. J. Phys. 22 601
level set and immersed interface methods J. Comput. Phys. 163 Xie Y M and Steven G P 1993 A simple evolutionary procedure for
489–528 structural optimization Comput. Struct. 49 885–96
Sigmund O and Maute K 2013 Topology optimization approaches a Xu S L, Cai Y W and Cheng G D 2010 Volume preserving nonlinear
comparative review Struct. Multidisc. Optim. 48 1031–55 density filter based on heaviside functions Struct. Multidisc.
Sigmund O and Petersson J 1998 Numerical instabilities in topology Optim. 41 495–505
optimization: a survey on procedures dealing with Zhou M and Rozvany G 1991 The COC algorithm: II. Topological,
checkerboards, mesh-dependencies and local minima Struct. geometrical and generalized shape optimization Comput.
Optim. 16 68–75 Methods Appl. Mech. Eng. 89 309–36

16