2014 Book MultiphysicsSimulation 2
2014 Book MultiphysicsSimulation 2
2014 Book MultiphysicsSimulation 2
Ercan M. Dede
Jaewook Lee
Tsuyoshi Nomura
Multiphysics
Simulation
Electromechanical System Applications
and Optimization
Simulation Foundations, Methods
and Applications
Tsuyoshi Nomura
Multiphysics Simulation
Electromechanical System Applications
and Optimization
123
Authors
Ercan M. Dede Tsuyoshi Nomura
Toyota Research Institute Toyota Central R&D Labs.
of North America Nagakute, Aichi
Ann Arbor, MI Japan
USA
Jaewook Lee
Korea Aerospace University
Goyang-si, Kyonggi-do
Korea, Republic of (South Korea)
Series editor
Prof. Emeritus Louis G. Birta
University of Ottawa
Ottawa, ON
Canada
v
vi Foreword
References
1. Bendsøe MP, Kikuchi N (1988) Generating optimal topologies in structural design using a
homogenization method. Comput Methods Appl Mech Eng 71(2):197–224
2. Bendsøe MP (1995) Optimization of structural topology, shape, and material. Springer, Berlin
3. Bendsøe MP, Sigmund O (2003) Topology optimization: theory, methods, and applications,
2nd edn. Springer, Berlin
Preface
This book developed out of a collaboration by the authors at the Toyota Research
Institute of North America, where multiphysics simulation and optimization is
used on a daily basis for a variety of engineering studies related to electrome-
chanical systems. Multiphysics simulation is a rapidly growing field, and the term
itself is broad and may be applied to an extremely wide variety of coupled-physics
problems. By nature, multiphysics simulation requires an array of technical skills
in different intersecting disciplines. As such, this book aims to narrow down the
topic by specifically focusing on multiphysics simulation for electromechanical
systems, the original target application investigated by the authors. It is our hope
that the collaborative aspects of such studies become apparent as the various
technical topics throughout the book are presented.
Overview
vii
viii Preface
Target Audiences
It is our hope that the content presented in this book will serve as a reference for
industry and academic researchers and engineers in the field of advanced elec-
tromechanical system design. The topics in this book are appropriate for under-
graduate and graduate level students, although many of the design examples may
be of interest to anyone curious about the unique design solutions that arise when
optimization methods are coupled with multiphysics simulation strategies.
Preface ix
Acknowledgments
We would like to thank our colleagues at the Toyota Research Institute of North
America (TRINA), Korea Aerospace University, and Toyota Central Research and
Development Labs (TCRDL), for their support in the completion of this project.
Professor Emeritus Louis G. Birta, Series Editor for Springer’s Simulation
Foundations, Methods and Applications book series, deserves special thanks for
the invitation to write this manuscript. Professor Gregory M. Hulbert from the
University of Michigan plus Dr. Danil Prokhorov and Koji Shiozaki from TRINA
are also acknowledged for their kind encouragement in undertaking this project.
Our appreciation goes to Satoru Sasaki from Toyota Motor Corporation for his
review of the text and technical feedback, Dr. Yan Liu from Toyota Technical
Center for discussions on thermal-fluid system design, and Prof. Xiaoping Qian
from the University of Wisconsin-Madison for interesting conversations on opti-
mization and CAE.
We further thank Dr. Kazuo Sato and Dr. Kunitoshi Nishikawa of TCRDL for
providing the opportunity and support for this collaborative work, Paul Schma-
lenberg from TRINA, Dr. Hisayoshi Fujikawa and Makoto Ohkado from TCRDL
for prototyping and measurement of the topology optimized microstrip coupler,
Prof. Shintaro Yamasaki of Osaka University for technical feedback on the level
set method and contribution of figures, Prof. Tatsuya Kashiwa of Kitami Institute
of Technology and Prof. Mohamed Bakr of McMaster University for discussions
on adjoint analysis with time domain solvers, Prof. Sang Won Yoon of Hanyang
University for the topology optimized DBC prototyping, Keiichi Shimaoka and
Takashi Ozaki of TCRDL for inviting an author to participate in their MEMS
project, Dr. Tadayoshi Matsumori and Dr. Atsushi Kawamoto of TCRDL for
discussions on filtering methods and vectorial topology optimization, and Prof.
Shinji Nishiwaki of Kyoto University for leading guidance in the field of topology
optimization, especially in Asia.
We also would like to thank Prof. Yoon Young Kim of Seoul National Uni-
versity, Prof. Seonho Cho of Seoul National University, Prof. Jeonghoon Yoo of
Yonsei University, Prof. Seungjae Min of Hanyang University, and Prof. Gil Ho
Yoon of Hanyang University for their support and advice on topology optimiza-
tion. Additionally, we wish to express special thanks to Prof. Joo Ho Choi and all
of the Professors at Korea Aerospace University for their encouragement in
completing this project.
A very special acknowledgment is given to Professor Noboru Kikuchi for his
kind support of this project and pioneering research in the field.
Finally, this book would not be possible without the encouragement and strong
support of our families and friends.
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Design via Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Single Physics Versus Multiphysics Simulation . . . . . . . . . . . . . 3
1.3 Challenges of Multiphysics Simulation. . . . . . . . . . . . . . . . . . . 5
1.4 The Role of Structural Optimization Methods . . . . . . . . . . . . . . 6
1.4.1 Topology Optimization . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4.2 Size and Shape Optimization . . . . . . . . . . . . . . . . . . . . 8
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
xi
xii Contents
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
Acronyms
2-D Two-dimensional
3-D Three-dimensional
ABC Absorbing boundary condition
AC Alternating current
ALE Arbitrary Lagrangian–Eulerian
AVM Adjoint variable method
BC Boundary condition
CAD Computer-aided design
CAE Computer-aided engineering
CFD Computational fluid dynamics
CTE Coefficient of thermal expansion
CVW Coulomb Virtual Work
DARPA Defense Advanced Research Projects Agency
DBC Direct bonded copper
DC Direct current
DOF Degrees-of-freedom
DRA Dielectric resonator antenna
EMC Electromagnetic compatibility
EMI Electromagnetic interference
FDTD Finite difference-time domain
FEA Finite element analysis
FE-BI Finite element-boundary integration
FEM Finite element method
GA Genetic algorithm
HPC High performance computing
HV Hybrid vehicle
IC Integrated circuit
IPMSM Interior permanent magnet synchronous motor
LSM Level set method
MEMS Microelectromechanical systems
MMA Method of Moving Asymptotes
MST Maxwell Stress Tensor
NASA National Aeronautics and Space Administration
OC Optimality criteria
xiii
xiv Acronyms
Scalar Quantities
As Surface area
B Magnitude of magnetic flux density
Br PM Permanent magnet strength
Brx x-direction component of residual
magnetic flux density
Bry y-direction component of residual
magnetic flux density
C Curie constant
Cp Specific heat capacity
D Extended design domain
E Elastic (or Young’s) modulus
F Force
f Frequency
fi Inclusion volume fraction
g Behavior constraint
H Magnitude of magnetic field
h Surface convection coefficient
hc Channel height
ht Heaviside function transition bandwidth
I, i Electric current
J Electric current density
k Thermal conductivity
kh Hysteresis coefficient
L Inductance
l Mechanical compliance
Mw Molecular weight
N Nodal shape function
Nc Number of coils
P Pressure
p Penalization parameter
Pr Number of rotor poles
Ps Number of stator poles
P Power
xv
xvi Symbols
v Poisson’s ratio
X Domain
Xd Design domain
Xm Material domain
x Angular frequency or velocity
U Electric scalar potential
UB Magnetic flux
q Density
qe Electric resistivity
qq Electric charge density
ς Electric conductivity
h Inclusion angle
hr Rotor angle
ζ Channel aspect ratio
Vector, Matrix, or Tensor Quantities
A Magnetic vector potential
B Magnetic flux density
Br Residual magnetic flux density
Be Array of derivatives of element nodal
shape functions
C Stiffness tensor
D Electric displacement
E Electric field
E Electric field variable
ê Target unit vector
F,f Force (global, element)
H Magnetic field
J Electric current density
Je External current density
Jeddy Eddy current density
K,k Stiffness matrix (global, element)
^
n Normal unit vector
Q Magnetic load vector
q Heat flux vector
t Surface load vector
^t Tangential unit vector
U,u Displacement (global, element)
v Velocity
x Position
–u Virtual displacement
e Strain
^
λ Lagrange multipliers (adjoint variables)
Stress
xviii Symbols
Functions
Fo Optimization objective function
Fc Convergence function
H Heaviside function
~
H Regularized Heaviside function
L Lagrangian
` Level set function
ˆ Scalar function
~
ˆ Filtered scalar function
Subscripts
avg Average
d Distributed
F,f Fluid
in Inlet
init Initial
l Lower
M Magnetic
max Maximum
min Minimum
o Reference value
out Outlet
S Structural
s Solid
T Thermal
u Upper
Chapter 1
Introduction
Starting in the mid-1950s and progressing into the following two decades, the
mathematical framework behind the finite element method was developed by
researchers due to the potential of the method as a revolutionary tool in the analysis
and design of civil and aerospace structures [17]. Building off such early work [35],
and under financial support from the US federal government, NASTRAN®1 finite
element analysis (FEA) software was developed by NASA in conjunction with indus-
try for aerospace applications. Today, FEA is a highly developed and commercially
available computational approach that allows for the estimation of the response of
structures that are subjected to multiple physical processes.
In this chapter, the role of FEA in design via simulation is explored and the impor-
tance of this approach for advanced electromechanical design is emphasized. Rep-
resentative differences between single and multiphysics simulations are presented
in the context of an example electromechanical system. Common challenges asso-
ciated with multiphysics simulation are then described to motivate the discussion
of the various computational tools presented throughout this book. An introductory
overview of the selected optimization techniques that are presented later on in the
text is additionally provided.
The use of FEA in the early stages of a product design process has become common
across multiple industries including automotive, aerospace, and consumer electron-
ics. Designing complex systems via simulation can lead to more efficient designs
and significant reduction in product development cycle times as evidenced by recent
aircraft structural design studies. In [14], the authors describe the use of advanced
numerical topology, shape, and size optimization techniques to generate efficient
stress and stability designs for aircraft leading edge droop nose ribs. This example
Physics A
Analysis
Concept User Initial Prototype Final
(CAD) Evaluation Design Validation Design
Physics B
Analysis
Fig. 1.1 A design cycle process flow utilizing FEA for performance evaluation
Fig. 1.2 A design cycle process flow utilizing advanced multiphysics simulation and optimization
technologies
of the synthesized system concept [24]. From there, the designer evaluates perfor-
mance and iterates on the proposed concept either during the analysis or prototype
evaluation stage. It is often the case that several iterations are necessary to satisfy
the performance targets, where major redesigns become increasingly costly once at
the prototype evaluation stage.
While the use of theoretical techniques and engineering experience is invaluable,
the goal in adopting advanced simulation and optimization technologies early in the
design stage of a project, Fig. 1.2, is to avoid costly iterations at later phases of the
project where less design flexibility is readily available. This ideal scenario is even
more imperative as electromechanical systems become more highly integrated and
space constrained, since analytical models may be limited and numerous prototypes
must be fabricated when adopting a purely experimental approach in order to fully
understand interactions arising from multiple physical scales and processes. In view
of the foregoing, the need for advanced multiphysics simulation and optimization
techniques is real and continually growing.
In this book, we aim to describe a set of numerical tools that enables efficient
early stage research and development concepts for product design via multiphysics
simulation and optimization. These tools may then be applied to a broad variety
of electromechanical system applications, and this point is emphasized through the
description of several example studies that are provided in later chapters of the text.
Inverter
/Converter
Assembly Battery
Motor/Generator
Internal
Assembly
Combustion
(Power Split &
Engine
Reduction
Devices)
Mechanical
power path
Electrical
power path
Transmission
Fig. 1.3 Example hybrid vehicle drive-train schematic including traditional internal combustion
engine, transmission (with motor/generator), power control unit (comprising inverter/converter), and
hybrid vehicle battery. The four starred items in the figure indicate multiphysics systems involv-
ing some combination of heat transfer, fluid flow, thermal-stress, or electromagnetic phenomena.
Dynamics, vibration, and control are also critical aspects to such systems and must be considered
in the design process
30]. The availability of new multiphysics FEA software plus advanced workstations
and HPC systems enables accurate coupled-field analysis for electronics, motors,
and actuators. Here, electro-thermal, thermal-fluid, thermal-structural, and electro-
magneto-thermal interactions (to name a few) are common and often lead to greater
system complexity in addition to significant system efficiency and reliability chal-
lenges [15].
As an example of a multiphysics system, a representative hybrid vehicle (HV)
drive-train schematic is shown in Fig. 1.3, where the starred items indicate sub-
systems that involve some combination of heat transfer, fluid flow, thermal stress, or
electromagnetic phenomena [18, 21]. For the inverter electronics, managing steady-
state and transient thermal loads plus thermally induced stresses is of prime impor-
tance, while increasing motor torque and simultaneously minimizing torque ripples
is a key challenge for the electric motor [7].
This drive-train example illustrates the importance of physics coupling in advanced
simulation. Specifically, it is not possible to anticipate and mitigate the stresses
induced by operating an electronics system unless the designer understands the power
dissipation and thermal loads that are generated by applying electrical power to the
device package. Additionally, the introduction of a design variable (as in topology
optimization) may also create coupling between the actual physics-based problem
1.2 Single Physics Versus Multiphysics Simulation 5
As explained in [15], the primary benefits and challenges associated with multi-
physics simulation stem from coupled-field analysis, where multiple physical phe-
nomena may be modeled simultaneously. Prior to the development of multiphysics
FEA software, coupled-field analysis was complicated, where the solution from one
physics model and software program had to be passed to another program to solve
for a second physical situation. This file transfer process required time and accuracy
was not necessarily guaranteed since user errors were easily propagated with such
data exchanges. Fortunately, this challenge has mostly been overcome through the
creation of commercially available software packages; see for example [2, 5, 9].
Further examination of coupled-field analysis reveals that two types of coupling
generally exist for multiphysics problems including direct coupling and sequential
coupling [15]. In a direct coupled system, a single matrix system of equations based
on all of the relevant physics is assembled and then solved. One drawback for direct
coupled systems is that finding a solution can be costly in terms of required process-
ing power and computer memory requirements. Alternatively, a segregated solver
may be used for a sequentially coupled system, where the solution from a first set
of field equations is passed to a second set of field equations, which is solved and
then passed to a third set of equations, etc. This segregated process is repeated until a
final solution is obtained. Many of the commercial multiphysics FEA packages now
automatically ‘suggest’ an appropriate default coupling depending on the physics
involved, although the user may customize these settings, as needed, based on solu-
tion time, available computing power, and numerical convergence of the problem at
hand.
As discussed in the prior section, the introduction of a design variable logically
creates coupling between the physics-based problem and a supplementary optimiza-
tion problem. The efficiency of this coupling scheme again affects overall solution
time. For topology optimization, this coupling is typically handled in a sequential
fashion, where information from the FEA solution is passed to the optimizer and
vice-versa until some convergence criterion is satisfied. Interestingly, there is also
recent work related to transforming the design variable into an equivalent state vari-
able in the solution of a single direct coupled physics-based problem [12]. Such
an approach has inherent challenges but also has significant potential for increased
useability and reduced computational cost.
Perhaps one of the greatest challenges with any multiphysics simulation is under-
standing all of the relevant physics involved in a particular problem and setting up a
geometrically accurate model with appropriate loads and boundary conditions. Incor-
porating multiple physical phenomena into a simulation may require consultation
6 1 Introduction
Structural Optimization
Fig. 1.4 The role of size, shape, and topology in structural optimization; ref. [10]. Reprinted from
[7, Fig. 1], Copyright (2012), with permission from Inderscience Enterprises Limited
with experts in different fields, while modern software tools allow for direct import
of more realistic three-dimensional (3-D) solid model computer-aided design (CAD)
geometry. To assist in the process just described, the multiple physics that are com-
monly considered for electromechanical system applications are outlined in Chap. 2,
while governing equations for such systems are reviewed in Chap. 3. Having estab-
lished the proper physics and geometry, the standard finite element method (FEM)
steps, per [17], then apply and include: (1) pre-processing or model setup, (2) numer-
ically solving the problem, and (3) post-processing or evaluation of computational
results.
The field of structural optimization is extremely broad and generally concerned with
the definition of the shape, size, and topology of a structure, as discussed in [10]
and illustrated in Fig. 1.4. Broadly, the goal in a structural optimization problem
is to find the optimal set of assumed design variables that maximize performance
subject to predefined design constraints (e.g., weight, size, cost, etc.). Appropriately,
a significant amount of research has been published encompassing topics that span
a range of physical scales from the optimization of large truss structures to the
optimal design of material microstructure. Efficient numerical methods for such
structural optimization applications are well established and utilize both gradient-
based searches [3, 4] and evolutionary algorithms [8, 10].
In this manuscript, we focus primarily on the use of gradient-based methods for
the structural topology optimization of electromechanical systems, although some
pertinent parametric size analysis examples are also presented. The main significance
of these numerical studies is that they are representative examples of how automated
methods may be exploited in arriving at optimal structural topologies for further
design verification within a commercial multiphysics CAE software environment.
1.4 The Role of Structural Optimization Methods 7
F F
Fig. 1.5 A standard structural compliance minimization topology optimization problem is shown
on the left for a simply supported beam. The optimal topology for minimum compliance, assuming
a 50 % solid material volume fraction, is shown on the right
The origin of numerical topology optimization can be traced back to the highly
influential work of Bendsøe and Kikuchi [3] published in 1988. Topology optimiza-
tion has since become a rapidly developing research field, and the optimization of
both single and multiphysics systems is well established for a variety of applica-
tions including actuators [26, 27], sensors [13], phononic materials [29], compliant
mechanisms [32], fluid flow devices [22, 33], and fluid-structure interaction prob-
lems [34], to name a few. Given the rate at which computing power is increasing,
such optimization techniques are prime for utilization in the design and development
cycle of advanced electromechanical systems.
As explained in detail in [4], topology optimization is a method in which a math-
ematical representation of the geometry is assumed and gradient information within
a finite element framework is used to guide the topological layout of a design. Dis-
cretization of a structure into many finite elements is typically required for topology
optimization, and each element is assigned a design variable value, γ . For a sin-
gle physics case, these design variables are often associated with a single material
physical parameter. A common example is a static structural problem, and a typical
objective function, Fo , is to minimize the structural compliance (i.e., deflection). In
Fig. 1.5, a simply supported beam is shown with an applied load, F, in the middle.
Here, the elastic modulus, E, is interpolated between γ = 0 (void) and γ = 1 (solid).
Additionally, a 50 % solid material volume constraint is assigned in order to deter-
mine the optimal structural topology that minimizes deflection for a given material
volume fraction; refer to the illustration on the right in Fig. 1.5. As explained in [28],
a power law interpolation function, E = γ p E o , is commonly used to control the
material state, where E o represents the elastic modulus of a given isotropic material
and p is a penalty parameter typically set to 3.
The process defined above is similar for a multiphysics problem except that addi-
tional design variables, constraints, and/or objective function terms may be consid-
ered. As an example, in a two-physics problem, the designer may define a two-term
objective function that is related to two separate state variables within the design
domain. Each term within this objective function may then have a weighting value
associated with it, and a family of ‘optimal’ solutions may then be developed depend-
ing on the relative weighting or priority of each objective function term. This family
8 1 Introduction
w1 > w2
F o_1 w2 > w1
Fo_2
Fig. 1.6 A simple Pareto front example diagram illustrating the tradeoff between different objective
function terms, Fo_1 , Fo_2 , in finding the optimal solution in a multi-objective optimization problem.
The vertical arrow indicates the effect of prioritizing the weighting value, w1 , for the first objective
function term, Fo_1 , while the horizontal arrow shows the effect of prioritizing the weighting value,
w2 , for the second objective function term, Fo_2
Size and shape optimization also play a crucial role in the design of multiphysics
systems. Assuming that the specific component size or shape variable (e.g., part thick-
ness or external boundary curve/spline) that one is interested in can be adequately
parametrized, many of today’s advanced multiphysics finite element software pro-
grams allow for the parametric optimization of multiple variables [1, 2, 5]. Here, the
topology (or layout) of a structure is generally known, and specific components are
optimized in terms of their relative size or shape to achieve a particular objective.
A classic multiphysics example of size optimization for electromechanical sys-
tems can be found in the context of layered electronics packages, where thermally
induced stresses play a key role. Thermal stress (along with vibration) is one of the
leading causes of electronics failures [11]. An underlying cause of thermal stress is
the mismatch in the coefficient of thermal expansion (CTE) of the various materials
comprising an electronics package [31]. For example, an order of magnitude differ-
ence exists between the CTE of silicon (2.6 ppm/◦ C) and tin-lead solder (27 ppm/◦ C).
Additionally, in a laminated structural assembly, the length and thickness of the vari-
ous material layers is important. Thus, these geometric parameters may be efficiently
optimized or ‘tuned’ using standard numerical tools to minimize thermal-structural
effects due to CTE mismatch.
Further details regarding parametric size analysis for electromechanical systems
are provided in Chap. 4, and two example studies are presented in Chap. 5. While
parametric size analysis/optimization is primarily covered in this book, many of the
numerical strategies may be logically extended to parametric shape optimization, as
well.
References
1. Altair HyperWorks, Altair Engineering Inc. 1820 Big Beaver Rd. Troy, MI 48083
2. ANSYS Multiphysics, ANSYS Inc. Southpointe, 275 Technology Drive, Canonsburg, PA
15317
3. Bendsøe MP, Kikuchi N (1988) Generating optimal topologies in structural design using
a homogenization method. Comput Method Appl M 71:197–224. doi:10.1016/0045-
7825(88)90086-2
4. Bendsøe MP, Sigmund O (2003) Topology optimization: theory, methods, and applications,
2nd edn. Springer, Berlin
5. COMSOL Multiphysics, COMSOL, Inc. 1 New England Executive Park, Suite 350, Burlington,
MA 01803
6. DARPA (2014) Tactical technology office: advanced vehicle make (AVM). http://www.darpa.
mil/Our_Work/TTO/Programs/Adaptive_Vehicle_Make__(AVM).aspx. Accessed 12 March
2014
7. Dede EM, Lee J, Liu Y, Robert B, Yönak SH (2012) Computational methods for the optimiza-
tion and design of electromechanical vehicle systems. Int J Vehicle Des 58(2–4):159–180.
doi:10.1504/IJVD.2012.047383
8. Goldberg DE (1989) Genetic algorithms in search, optimization and machine learning.
Addison-Wesley Publishing, Reading
9. Isight, Dassault Systèmes. 10, Rue Marcel Dassault, 78140 Vélizy-Villacoublay, France
10. Jakiela MJ, Chapman C, Duda J, Adewuya A, Saitou K (2000) Continuum structural topology
design with genetic algorithms. Comput Method Appl M 186:339–356. doi:10.1016/S0045-
7825(99)00390-4
11. Jiang ZQ, Huang Y, Chandra A (1997) Thermal stresses in layered electronics assemblies. J
Electron Packaging 119:127–132. doi:10.1115/1.2792218
12. Kawamoto A, Matsumori T, Nomura T, Kondoh T, Yamasaki S, Nishiwaki S (2012) Topology
optimization by a time-dependent diffusion equation. Int J Numer Meth Eng 93:795–817.
doi:10.1002/nme.4407
10 1 Introduction
13. Kim JE, Kim DS, Ma PS, Kim YY (2010) Multi-physics interpolation for the topology opti-
mization of piezoelectric systems. Comput Method Appl M 199:3153–3168. doi:10.1016/j.
cma.2010.06.021
14. Krog L, Tucker A, Rollema G (2002) Application of topology, sizing and shape optimization
methods to optimal design of aircraft components. Altair Engineering, Coventry
15. Lethbridge P (2004) Multiphysics analysis. Ind Phys 12:26–29
16. Littmarck S (2001) Solving differential equations. Ind Phys 2:21–23
17. Logan DL (2002) A first course in the finite element method, 3rd edn. Brooks/Cole, Pacific
Grove
18. Matsubara T, Yaguchi H, Takaoka T, Jinno K (2009) Development of new hybrid system for
compact class vehicles. SAE Technical Paper 2009–01-1332. doi: 10.4271/2009-01-1332
19. MSC Software (2014) MSC Nastran. http://www.mscsoftware.com/product/msc-nastran.
Accessed 12 March 2014
20. Nomura T, Sato K, Nishiwaki S, Yoshimura M (2007) Multi-disciplinary multi-objective topol-
ogy optimization of electromagnetics and structural mechanics: for case of optimal dielectric
resonator antenna designs. T Jpn Soc Mech Eng A 73:1111–1119. doi:10.1299/kikaia.73.1111
21. Nozawa N, Maekawa T, Nozawa S, Asakura K (2009) Development of power control unit for
compact-class vehicle. SAE Int J Passeng Cars Electron Electr Syst 2:376–382. doi:10.4271/
2009-01-1310
22. Olesen LH, Okkels F, Bruus H (2006) A high-level programming-language implementation
of topology optimization applied to steady-state Navier-Stokes flow. Int J Numer Meth Eng
65:975–1001. doi:10.1002/nme.1468
23. Savander BR, Shankara P (2012) Propulsion system performance optimization—design by
analysis. CD-adapco Dynamics 33:47–50
24. Schoofs AJG (1993) Structural optimization history and state-of-the-art. In: Dijksman JF,
Nieuwstadt FTM (eds) Topics in applied mechanics. Kluwer Academic Publishers, Nether-
lands, pp 339–345
25. Shigley JE, Mischke CR, Budynas RG (2004) Mechanical engineering design, 7th edn.
McGraw-Hill, New York
26. Sigmund O (2001) Design of multiphysics actuators using topology optimization—Part
I: one-material structures. Comput Method Appl M 190:6577–6604. doi:10.1016/S0045-
7825(01)00251-1
27. Sigmund O (2001) Design of multiphysics actuators using topology optimization—Part
II: two-material structures. Comput Method Appl M 190:6605–6627. doi:10.1016/S0045-
7825(01)00252-3
28. Sigmund O (2001) A 99 line topology optimization code written in Matlab. Struct Multidiscip
O 21:120–127. doi:10.1007/s001580050176
29. Sigmund O, Jensen JS (2003) Systematic design of phononic band-gap materials and structures
by topology optimization. Philos T Roy Soc A 361:1001–1019. doi:10.1098/rsta.2003.1177
30. Thilmany J (2010) Multiphysics: all at once: physical phenomena-and engineers-rarely work
in isolation, so simulation software is addressing those facts. Mech Eng CIME 132:39–41
31. Timoshenko S (1925) Analysis of bi-metal thermostats. J Opt Soc Am 11:233–255. doi:10.
1364/JOSA.11.000233
32. Yin L, Ananthasuresh GK (2002) A novel topology design scheme for the multi-physics prob-
lems of electro-thermally actuated compliant micromechanisms. Sensor Actuat A-Phys 97–
98:599–609. doi:10.1016/S0924-4247(01)00853-6
33. Yoon GH (2012) Topological layout design of electro-fluid-thermal-compliant actuator. Com-
put Method Appl M 209–212:28–44. doi:10.1016/j.cma.2011.11.005
34. Yoon GH, Jensen JS, Sigmund O (2007) Topology optimization of acoustic-structure interaction
problems using a mixed finite element formulation. Int J Numer Meth Eng 70:1049–1075.
doi:10.1002/nme.1900
35. Zienkiewicz OC, Taylor RL, Zhu JZ (1967) The finite element method: its basis and funda-
mentals, 1st edn. Elsevier Butterworth-Heinemann, Burlington
Chapter 2
Overview of Physics for Electromechanical
Systems
Fig. 2.1 Range of physical scales encountered in the design of electromechanical systems
Electromechanical
Systems
Fig. 2.2 Commonly encountered physics couplings for electromechanical system components
including electronic systems, low frequency magnetics, radio frequency (RF) components, and
motors and actuators
Modern electronic equipment is highly complex, and the effect of electric current
flow and the resulting device power dissipation on the temperature and reliability of
the package is significant. In Fig. 2.3, a representative structure is shown to illustrate
the typical components found in an electronics package. While this package features
a planar integrated circuit (IC) device, the main features of the structure are similar
to those found in packages with discrete devices.
Generally, electric current is applied to the device, which generates an electro-
magnetic field and dissipates power in the form of heat. The device is attached to a
2.1 Electronic System Components 13
Device Electro-magnetics
Substrate
+
Bond Mechanics
Layers
Thermal
Fluidics
Cold Plate
Fig. 2.3 Cross-section schematic of a representative electronics package highlighting the physics
involved with various components. Although a planar type of structure is shown, the illustrated
components (i.e., device, substrate, bond layer/TIM, and a cold plate/heat sink) are common to
most electronics applications
The above descriptions highlight the host of multiphysics challenges that arise in
the design of electronic equipment. The governing equations for electronic systems
involving Joule heating, thermally induced stress, and conjugate heat transfer are
provided in Chap. 3 along with several example design studies in Chap. 5. Note that
other single physics phenomenon such a vibration are also important in assessing the
reliability of an electronics package; however, the reader is referred to the existing
literature [23, 24] since this topic has been extensively studied.
Wire Electric
+
Core Thermal
Heat Magnetic
sink
Fig. 2.4 Side view schematic of a representative low-frequency magnetic device with defined
physics couplings. Although an I-shaped type of inductor structure is shown, the illustrated com-
ponents (i.e., wire, core, and heat sink) are common to most magnetic components
2.3 RF Components
Recently, electric motors and actuators have been extensively used as fundamental
components of future vehicle drive-trains and robotic systems. A possible scenario is
that electric motors will eventually replace traditional internal combustion engines for
future eco-friendly vehicles. Actuators are essential in the development of drive-by-
wire technologies, which aim to replace traditional mechanical control systems with
electronic systems. Figure 2.6 illustrates a schematic of a linear actuator structure
that is presented in a simplified form for ease of understanding. The permanent
magnet (PM) and electro-magnet (composed of yoke and wires) generate a magnetic
field loop that passes within the plunger and through the yoke (or rotor and stator,
respectively, for a motor). The generated magnetic field causes a magnetic force
18 2 Overview of Physics for Electromechanical Systems
Plunger Mechanics
Motion (Rotor)
Permanent
Magnet Thermal
Wire
Magnetic
Yoke
(Stator)
Electric
Heat
sink
Fig. 2.6 Conceptual schematic view of a simple linear actuator. The featured components of the
device (i.e., permanent magnet, wire, etc.) are common to not only actuators, but electric motors
that acts on the plunger. Consequently, the electric energy supplied to the wires
is converted into mechanical energy at the plunger through the magnetic energy
surrounding the device. During this energy conversion processes, thermal energy is
produced due to power loss in the electrical wires and magnetic materials (i.e., yoke
and plunger). Logically, a heat sink (or cooling jacket) is then attached for thermal
management of the device.
The physical interactions that occur in motors and actuators are very similar to
those present in magnetic components with the exception of the mechanics-related
physics. Low-frequency electromagnetic field coupling is again a fundamental phe-
nomenon that must be considered in the performance prediction of motors and actu-
ators. The heat generation effect due to Joule and hysteresis losses (i.e., temperature
rise due to heat dissipation) is identical to the effect explained in relation to magnetic
components. The main difference is that motors and actuators are devices that gener-
ate force (or torque) and thus mechanical power. Therefore, it is necessary to under-
stand the dynamics, structural, and vibration characteristics related to the plunger
(or rotor) movement. Also, this mechanical behavior may have cascading effects and
influence the electromagnetic characteristics. Thus, the coupled electromagnetic-
thermal-mechanical problem needs to be addressed in the analysis and design of
motors and actuators.
An active research area with regard to motors and actuators is the development of
methods for simultaneously increasing power density (i.e., generated power amount
per unit device size) and efficiency (i.e., mechanical energy output per consumed elec-
trical energy input). The enhancement of these two performance metrics is required
in order to improve the mobility of motors and actuators. Many extensive research
studies have been conducted, and one technical approach is the development of new
materials based on nano-technologies. Such materials hold promise to offer better
material properties (e.g., permanent magnets that have higher residual magnetic flux
density or iron that has higher magnetic permeability with low electric conductiv-
ity). Improving electric power conversion devices such as inverters and converters
2.4 Motors and Actuators 19
is a second potential approach to achieve high power density and efficiency, where
engineers seek an optimal electrical energy input waveform. Lastly, design opti-
mization based on multiphysics simulation, as addressed in this book, may also be
an effective technique for motor and actuator performance enhancement.
Simulation methods that consider electromagnetics and force generation, relative
to power density and efficiency, have been extensively studied; see for example [1,
2, 19, 21]. Thus, related governing equations for the multiphysics analysis of motors
and actuators will be provided in Chap. 3 along with several design example studies
in Chap. 5. However, it should be mentioned that understanding thermal and vibration
performance is another critical piece of the puzzle in the analysis and design of motors
and actuators. Fully coupled simulation methods that handle electromagnetics, force
generation, heat transfer, and structural dynamics are not yet mainstream due to their
inherent complexity, although future research into motors and actuators may address
all of these effects.
References
1. Bastos JPA, Sadowski N (2003) Electromagnetic modeling by finite element methods. Marcel
Dekker, New York
2. Bianchi N (2005) Electrical machine analysis using finite elements. CRC Press, Boca Raton
3. Bossche AVD, Valchev VC (2005) Inductors and transformers for power electronics. CRC
Press, Boca Raton
4. Chmielak B, Waldow M, Matheisen C, Ripperda C, Bolten J, Wahlbrink T, Nagel M, Merget F,
Kurz H (2011) Pockels effect based fully integrated, strained silicon electro-optic modulator.
Opt Express 19:17212–17219. doi:10.1364/OE.19.017212
5. Del Vecchio RM, Poulin B, Feghali PT, Shah DM, Ahuja R (2010) Transformer design prin-
ciples. CRC Press, Boca Raton
6. Fay CE, Comstock RL (1965) Operation of the ferrite junction circulator. IEEE Trans Microw
Theory Tech 1:15–27. doi:10.1109/TMTT.1965.1125923
7. Flanagan WM (1992) Handbook of transformer design and applications, 2nd edn. McGraw-
Hill, New York
8. Gad-el-Hak M (2002) The MEMS handbook. CRC Press, Boca Raton
9. Gianchandani Y, Tabata O, Zappe H (2008) Comprehensive microsystems, vol 1–3. Elsevier,
Amsterdam
10. Hache F, Ricard D, Flytzanis C, Kreibig U (1988) The optical kerr effect in small metal
particles and metal colloids: the case of gold. Appl Phys A Mater Sci Process 47:347–357.
doi:10.1007/BF00615498
11. Hsu T-R (2008) MEMS and microsystems: design, manufacture, and nanoscale engineering,
2nd edn. Wiley, Hoboken
12. Incropera FP (1999) Liquid cooling of electronic devices by single-phase convection. Wiley-
Interscience, New York
13. Joannopoulos JD, Meade RD, Winn JN (1995) Photonic crystals: molding the flow of light,
1st edn. Princeton University Press, New Jersey
14. Kays WM, London AL (1998) Compact heat exchangers. Krieger Publishing Company,
Malabar
15. Kim S-J, Lee S-W (1996) Air cooling technology for electronic equipment. CRC Press, Boca
Raton
16. Lin SY, Moreno J, Fleming J (2003) Three-dimensional photonic-crystal emitter for thermal
photovoltaic power generation. Appl Phys Lett 83:380. doi:10.1063/1.1592614
20 2 Overview of Physics for Electromechanical Systems
17. Mikami O, Nakagome H (1984) Waveguided optical switch in InGaAs/InP using free-carrier
plasma dispersion. Electron Lett 20:228–229. doi:10.1049/el:19840153
18. Remsburg R (1998) Advanced thermal design of electronic equipment. International Thomson
Publishing, Florence
19. Salon SJ (1995) Finite element analysis of electrical machines. Kluwer Academic Publishers,
Boston
20. Senturia SR (2001) Microsystem design. Springer, New York
21. Silvester PP, Ferrari RL (1996) Finite elements for electrical engineers, 3rd edn. Cambridge
University Press, New York
22. Soref R (2006) The past, present, and future of silicon photonics. IEEE J Sel Top Quant
12:1678–1687. doi:10.1109/JSTQE.2006.883151
23. Steinberg DS (2000) Vibration analysis for electronic equipment. Wiley, New York
24. Steinberg DS (2001) Preventing thermal cycling and vibration failures in electronic equipment.
Wiley, New York
25. Tlusty T, Meller A, Bar-Ziv R (1998) Optical gradient forces of strongly localized fields. Phys
Rev Lett 81:1738–1741. doi:10.1103/PhysRevLett.81.1738
26. Wedlock BD (1963) Thermo-photo-voltaic energy conversion. Proc IEEE 51:694–698. doi:10.
1109/PROC.1963.2261
Chapter 3
Governing Equations for Electromechanical
Systems
In this chapter, the governing equations for the electromechanical systems described
in Chap. 2 are explained in detail. Joule heating, thermal stress, and conjugate heat
transfer problems are first introduced. Low frequency electromagnetic systems are
then addressed followed by high frequency electromagnetics. Prior to providing
the detailed equations for each physical system, a conventional forward analysis
solution of a single physics structural mechanics example problem is provided as a
point of reference for the steps involved in standard finite element analysis. For each
multiphysics system that then follows, a system diagram, the governing equations,
and typical system loads plus boundary conditions are provided.
The utility of the finite element method is that it may be used to produce an approx-
imate solution to a set of equations that characterize the response of a structure by
discretizing it effectively into many elements or sub-regions. This approach is com-
mon for a broad set of problems that cover a range of physics, and a representative
FEA example is related to the computation of the deflection of a structure due to an
applied static load.
For illustrative purposes, a two-dimensional (2-D) simply supported (i.e., trans-
verse displacements constrained, but rotation and longitudinal displacement per-
mitted [15]) beam comprising isotropic linear-elastic material is assumed in this
introductory case, where a load is applied to the center of the structure, as shown in
Fig. 3.1. Here, the reader is assumed to have some background knowledge of struc-
tural mechanics [5] and the finite element method [7, 12, 19]. As described in [12],
the principle of virtual work is utilized to obtain the governing matrix system of
equations for the displacement of the structure
γU = γW, (3.1)
x
F
Displacement, u
Fig. 3.1 A simply supported beam with an applied load in the center. The beam is assumed to be
made of an isotropic linear-elastic material. The dashed line indicates the approximate deformed
shape of the lower edge of the beam
where γU is the virtual strain energy of a finite element due to internal stress and
γW is the virtual work done by the external forces on the element.
The virtual strain energy of a finite element may be written in matrix form in
terms of the virtual strains, γλ, and internal stresses, σ , as
γU = γλT σ dve , (3.2)
ve
λ = Be u, (3.3)
σ = Cλ. (3.4)
The virtual work done by the nodal structural forces, f S , acting on an element may
also be expressed as
γW = γuT f S , (3.6)
kS u = fS , (3.9)
where
kS = Be T CBe dve , (3.10)
ve
is the element structural stiffness matrix, u is the vector of element nodal displace-
ments, and f S is the vector of applied element structural forces. Once each element
stiffness matrix is computed, a standard assembly process is used to generate a global
structural stiffness matrix, K S [SI unit: N/m], and force vector, F S [SI unit: N], to
solve for the global displacements, U [SI unit: m], of a given structure
KS U = FS ; (3.11)
y
1m
0.1 m
x
Fy = -10 N
Pinned Roller uy = -1.17 x 10-8 m
ρ = 7850 kg/m3
E = 200 GPa
υ = 0.33
Fig. 3.2 The left side image illustrates the definition of the structural geometry, material physical
properties, boundary conditions, and applied loads for a simply supported beam required as part
of the FEA pre-processing; a 2-D plane-strain approximation is assumed for the planar model.
The right side image shows the post-processed deformed shape of the structure with y-direction
displacement contours (blue colored contours = larger displacement)
part of the modeling process, which along with the model BCs and loads defines
a system of algebraic equations.
• Numerical Analysis: Numerical algorithms are used to solve the linear or non-
linear algebraic system of equations defined in the preprocessing step.
• Post-processing: The output obtained through the prior analysis is post-processed
and conveyed through the generation of data and figures that quantify the state vari-
able results (e.g., deformed structural geometry, temperature field, and magnetic
field).
For clarity, the pre and post-processing steps are illustrated in Fig. 3.2. In the image
on the left, the structural geometry, material physical properties (including density,
ρ, elastic modulus, E, and Poisson’s ratio, ν), boundary conditions, and applied loads
for a simply supported beam are provided as part of the FEA pre-processing, where
a 2-D plane-strain approximation is assumed for the planar model [12]. The numer-
ical analysis then proceeds, where a convergence criteria for the selected solution
algorithm is typically satisfied. Finally, the vertical deflection and deformed shape of
the structure may be obtained through post-processing of the computational results,
as shown in the right side image in Fig. 3.2. For simple structural geometries, such
numerical results are easily compared with standard analytical predictions [15].
Note that the sizing of elements (plus associated mesh refinement) and the treat-
ment of concentrated point loads play an important role in obtaining an accurate
numerical solution. The specific physics involved, geometry details, nature of the
loads, boundary conditions, and numerical quantity of interest all affect the level
of discretization used in the numerical model and method of applying the loads.
Regions where large gradients occur in the state variables often require a finer mesh,
and the application of concentrated loads needs to be carefully considered relative
to the particular quantity of interest. In practice, however, every physical situation
is different and one should be wary of broadly applying general rules. The selected
modeling strategy should be based on prior experience (coursework and/or applied)
or appropriate engineering references to related work.
3.1 Single Physics Structural Mechanics Example 25
In summary, while the principle of virtual work was exploited in developing the
governing equations for the structural mechanics beam deflection problem above,
other approximation methods for the derivation of the governing equations for dif-
ferent physical systems may logically be used. Regardless, the overall ‘big-picture’
process explained above may be translated to most physical scenarios including the
multiphysics systems described in detail throughout the remainder of this chapter.
∂T
ρC p = ◦ · (k◦T ) + Q, (3.12)
∂t
where the conductor density, specific heat capacity, and thermal conductivity are
given as ρ [SI unit: kg/m3 ], C p [SI unit: J/(kg K)] and k [SI unit: W/(m K)], respec-
tively, and T [SI unit: K] is the temperature state variable. For response under steady-
state conditions, the term on the left hand side of Eq. (3.12) is zero.
Regarding Joule heating, the heat source term, Q, [SI unit: W/m3 ] on the right
hand side of this expression arises from an electromagnetic heat source
Q = J · E, (3.13)
where E [SI unit: V/m] is the electric field, and J [SI unit: A/m2 ] is the current
density.
More specifically, under steady-state conditions the divergence of the current
density is zero
◦ · J = 0. (3.14)
Furthermore, Ohm’s law provides for the current density in terms of the electric
field as
J = ς E + Je , (3.15)
26 3 Governing Equations for Electromechanical Systems
Fig. 3.3 Loads and boundary conditions for a general Joule heating problem
where ς is the electric conductivity [SI unit: S/m] of the conductive material, and Je
[SI unit: A/m2 ] is an externally generated current density.
Finally, substituting this expression into the continuity equation and generalizing
for current sources we have
◦ · (ς E + Je ) = Q j , (3.16)
where the electric field strength may be expressed in terms of the gradient of the
electric scalar potential, φ [SI unit: V], as E = −◦φ. The heat source due to Joule
heating, Q j , determined via Eq. (3.16) is substituted into Eq. (3.12) as the thermal
load to the system.
Following the comprehensive review of coupled electrothermal effects for related
microactuators [1] and MEMS [2], the solution of the electrothermal problem above
requires the specification of boundary conditions in addition to an initial temper-
ature, Tinit , and electric potential field, φinit . As shown in Fig. 3.3, either assumed
temperature, ∂ΩT , and electric potential, ∂Ωφ , (i.e., Dirichlet) boundary conditions,
or assumed normal heat flux, ∂Ωq⊥∗∗ , and electric current density, ∂Ω J⊥ , (i.e., Neu-
mann) boundary conditions are applied. In the special case of a convection heat
transfer boundary condition, the normal heat flux, q⊥ ∗∗ = h(T − T ), is proportional
o
to the temperature difference between the surface and the reference temperature, To ,
of an assumed surrounding fluid, where h [SI unit: W/(m2 K)] is the assumed sur-
face convection coefficient [8]. Although not specifically addressed here, the reader
is referred to the literature for a discussion of radiation heat transfer boundary
conditions [2, 8, 10].
3.3 Thermal Stress 27
Fig. 3.4 Loads and boundary conditions for a general thermal stress problem
The topic of thermal stress is broad, and in the context of electromechanical systems
there are many applications that could be considered. In particular, the subject of
thermal stress and strain in electronics packaging may be broken down into several
subcategories including the thermoelastic analysis of multilayered devices, determi-
nation of thermal stresses in thin films, evaluation of viscoplastic thermal damage
(i.e., cumulative life-cycle) effects in die attachment materials, and the study of the
thermal-mechanical loading of wirebond structures, to name a few [11]. However, to
narrow the field, the topic of thermal stress in this book is focused on the fundamental
case of the response of a thin laminated (i.e., multilayered) composite structure com-
posed of non-temperature-dependent isotropic linear-elastic materials. The reader is
referred to the literature for more complex topics that then build off of this essential
example.
Considering the above caveats, the loads and boundary conditions for a represen-
tative system are illustrated in Fig. 3.4, and the governing equations for quasi-static
thermoelasticity in a given domain, Ω, follow from [11]. These equations include
the heat conduction equation and the equation of motion for the system, which may,
respectively, be written as
∂T
ρC p = ◦ · (k◦T ) + Q, and (3.17)
∂t
− ◦ · σ = Fbody + Fα . (3.18)
In Eq. (3.17), T [SI unit: K] is the temperature state variable, k [SI unit: W/(m K)]
is the material thermal conductivity, and Q [SI unit: W/m3 ] is the volumetric power
28 3 Governing Equations for Electromechanical Systems
density. Additionally, the term on the left hand side of Eq. (3.17) is the time transient
term, which depends on the material density, ρ, and specific heat capacity, C p , as
discussed above in Sect. 3.2. In the latter expression, Eq. (3.18), σ [SI unit: N/m2 ] is
the matrix of thermally induced stresses, Fbody [SI unit: N/m3 ] is the vector of body
forces per unit volume, and Fα [SI unit: N/m3 ] is the vector of thermally induced
loads (per unit volume) specified as
Fα = Be Cα(T − To ), (3.19)
Similarly to the case of thermal stress, fluid mechanics in conjunction with ther-
mal transport is a very broad topic with many areas of subspecialty and numerous
authoritative texts; see, for example [6, 8, 14]. Here, the topic of conjugate heat trans-
fer is limited to the case of conduction plus convection heat transfer (i.e., radiation
effects are ignored), where all flows are considered to be viscous, incompressible,
and generally laminar. However, notwithstanding these conditions, some discussion
of conjugate heat transfer under turbulent flow will also be provided in relation to
logical extensions of the design studies presented herein.
Typical loads and boundary conditions for a domain, Ω, are shown in Fig. 3.5,
where the domain consists of a moving fluid region, Ω f , and a solid body, Ωs [14].
For steady-state fluid flows, the continuity equation (or incompressibility condi-
tion) is a mathematical statement of the principle of conservation of mass
3.4 Conjugate Heat Transfer 29
Fig. 3.5 Loads and boundary conditions for a general conjugate heat transfer problem
◦ · v = 0, (3.20)
ρ (v · ◦v) = −◦ P + ◦ · η ◦v + (◦v)T + ρf, (3.21)
where ρ is the fluid density [SI unit: kg/m3 ], η is the fluid dynamic viscosity [SI unit:
Pa s], P is the pressure state variable [SI unit: Pa], and f is the body force vector per
unit mass [SI unit: N/kg].
Finally, the First Law of Thermodynamics introduced in Sect. 3.2 above is used to
rewrite the steady-state heat equation with consideration of convective heat transfer
ρC p (v · ◦T ) = ◦ · (k◦T ) + Q, (3.22)
where C p represents the heat capacity [SI unit: J/(kg K)], k is the thermal conductivity
of the fluid [SI unit: W/(m K)], Q is the volumetric power density [SI unit: W/m3 ],
and T is the temperature state variable [SI unit: K]. Note that in the case of the solid
domain, Ωs , in Fig. 3.5, the fluid velocity vector v is zero and the left hand side term
in Eq. (3.22) disappears reducing the equation to the form
30 3 Governing Equations for Electromechanical Systems
− ◦ · (ks ◦T ) = Q s . (3.23)
Initial conditions for a conjugate heat transfer system typically involve an initial
temperature, Tinit , for the solid domain, Ωs , plus an initial temperature and hydrostatic
pressure, Pinit , for the fluid domain, Ω f , as shown in Fig. 3.5. Regarding boundary
conditions for the thermal and transport parts of the system shown, the same Dirichlet
(temperature) and Neumann (heat flux) boundary conditions previously discussed
for Joule heating and thermal stress systems again apply here for both fluid and
solid portions of the domain. Also, convection and radiation conditions, as described
in [14], may be applied. For the fluid dynamic portion of the system illustrated in
Fig. 3.5, either the velocity components, ∂Ωv , (i.e., Dirichlet boundary conditions)
or total surface stress, ∂Ωσ , (i.e., Neumann boundary conditions) must be specified
on the boundary of the fluid region; refer to [14] for additional explanation of these
and other specialized (e.g., periodic, slip, and interfacial stress) boundary conditions.
This section explains three types of analysis methods for low frequency electro-
magnetic systems based on the finite element method [4, 16, 17]. Here, the term
low frequency refers to the case in which the frequency of the electromagnetic field
is lower than several hundred kilohertz. In this low frequency situation, the wave-
length of the electromagnetic field is much longer than the device size, and thus
the wave characteristics of electromagnetic field becomes insignificant. The original
governing equations for electromagnetics, i.e., Maxwell’s equations can be written
as
∂B
◦ ×E=− Faraday’s law (3.24)
∂t
∂D
◦ ×H= +J Ampère’s law (3.25)
∂t
◦ · D = ρq Gauss’s law (3.26)
◦ ·B=0 Gauss’s law of magnetics, (3.27)
D =
E (3.28)
B = μH (3.29)
J = ς E, (3.30)
where E, H, B, D, J and ρq are the electric field [SI unit: V/m], the magnetic field [SI
unit: A/m], the magnetic flux density [SI unit: T], the electric displacement [SI unit:
C/m2 ], the electric current density [SI unit: A/m2 ], and the electric charge density [SI
unit: C/m3 ], respectively. Also,
, μ, and ς are, respectively, the electric permittivity
[SI unit: F/m], magnetic permeability [SI unit: H/m], and electric conductivity [SI
unit: S/m].
In the low frequency regime, the time dependent terms of the above Maxwell’s
equations are of little importance, and thus they may be ignored for simplification
depending on the situation and application. Based on the level of simplification,
it is possible to implement one of the three different analysis approaches for low
frequency applications.
The first analysis method for low frequency applications is electric and magnetic
field analysis. In this analysis method, the electric and magnetic fields are, respec-
tively, derived using the electric scalar potential, φ, and the magnetic vector potential,
A, from Eqs. (3.24), (3.26) and (3.27) without ignoring any time-dependent terms of
the original Maxwell’s equation
∂A
E = −◦φ − , (3.31)
∂t
B = ◦ × A. (3.32)
Then, both electric and magnetic potentials are obtained by solving Eq. (3.25) with
the constitutive laws, Eqs. (3.28)–(3.30). The boundary conditions for the electric
and magnetic field analysis are illustrated in Fig. 3.6. To understand the boundary
conditions, the analysis of a 3-D inductor can be considered. The voltage difference
at the terminals of the copper winding is modeled as the Dirichlet boundary condition
of electric potentials, φ. For the insulation of the magnetic field, the zero Dirichlet
condition of the magnetic potential, A, is applied at the boundary of the analysis
domain.
This first analysis method does not simplify the original Maxwell’s equations
and finds both electric and magnetic fields. It is useful to calculate both the current
density, J, and magnetic flux density, B, when the terminal voltages are specified in
the magnetic component. However, this generally requires three dimensional analysis
in order to obtain a meaningful result, which is computationally expensive. Thus,
this high computational cost may limit the use of this analysis method.
32 3 Governing Equations for Electromechanical Systems
Fig. 3.6 Boundary conditions for electric and magnetic field analysis
where ω is the frequency of the time-harmonic vector potential, A. The above time-
harmonic equation can be solved using the finite element method. In the 2-D case,
the finite element formulation of Eq. (3.34) can be represented in the matrix form as
kM A = fM . (3.35)
The element magnetic stiffness matrix, k M , and the magnetic force vector, f M ,
are, respectively, derived as
3.5 Low Frequency Electromagnetics 33
1 ∂ Ni ∂ N j ∂ Ni ∂ N j
and by solving this equation with Eq. (3.32), the time-harmonic magnetic flux density,
B, is obtained for the whole analysis domain.
The third and final analysis method for low frequency applications is magneto-
static analysis. Using this method, all time varying terms in the governing equations
are ignored by making use of the assumption of a static or quasi-static field, and the
magnetic field is calculated when the external current density, Je is given and/or the
permanent magnet generates the residual magnetic flux density, Br . From Eqs. (3.25)
and (3.32), the equation for the static magnetic vector potential, A, can be derived as
1 1
◦× ◦ × A = Je + ◦ × Br . (3.39)
μ μ
The above static equation may be solved again using the finite element method.
In the 2-D case, the element stiffness matrix, k M , and the force vector, f M , in the
finite element formulation can, respectively, be derived as
1 ∂ Ni ∂ N j ∂ Ni ∂ N j
∂ 1 ∂ 1
f M,i = (Je,z Ni )dydx − Br,y dydx + Br,x dydx.
∂x μ ∂y μ
x y x y x y
(3.41)
The matrix equation with Eqs. (3.40) and (3.41) gives us the static distribution of
the magnetic flux density, B, for the whole analysis domain.
Note that the loads and boundary conditions for both the second and third analysis
methods presented in this section are illustrated in Fig. 3.7. The external current
density, Je , is applied at the winding areas as the given load. Then, the magnetic field
distribution generated by the given external current density is calculated considering
eddy currents or with static conditions. For the insulation of the magnetic field, the
zero Dirichlet condition is applied at the boundary of the whole analysis domain.
34 3 Governing Equations for Electromechanical Systems
Fig. 3.7 Loads and boundary conditions for either the eddy current or magnetostatic analysis
method
In this section, two popular numerical methods are briefly described. One is the
Finite Difference-Time Domain method, FDTD [18], and the other is the FEM [9].
The FDTD method is a very popular method for time transient analysis, while the
FEM is usually combined with the frequency domain formulation for high frequency
electromagnetic simulations.
As mentioned in the previous section, electromagnetic systems are described with
Maxwell’s equations
∂B
◦ ×E=− Faraday’s law (3.42)
∂t
∂D
◦ ×H= +J Ampère’s law (3.43)
∂t
◦ · D = ρq Gauss’s law (3.44)
◦ ·B=0 Gauss’s law of magnetics, (3.45)
D =
E (3.46)
B = μH (3.47)
J = ς E, (3.48)
3.6 High Frequency Electromagnetics 35
where E, H, B, D, J and ρq are the electric field [SI unit: V/m], the magnetic field [SI
unit: A/m], the magnetic flux density [SI unit: T], the electric displacement [SI unit:
C/m2 ], the electric current density [SI unit: A/m2 ], and the electric charge density [SI
unit: C/m3 ], respectively. Also,
, μ and ς are, respectively, the electric permittivity
[SI unit: F/m], magnetic permeability [SI unit: H/m], and electric conductivity [SI
unit: S/m].
In the high frequency regime, the time dependent part of the above formulation
dominates; that is, Faraday’s law, Eq. (3.42), and Ampère’s law, Eq. (3.43). With the
Finite Difference-Time Domain method, these two equations are directly discretized
with finite differences.
With Eqs. (3.42) and (3.43) and the constitutive laws, we obtain
∂H
◦ × E = −μ (3.49)
∂t
∂E
◦ ×H=
+ ς E, (3.50)
∂t
and the FDTD discretization yields the following equations for the field in x-direction
t
E x |i,n+1
j,k = E x |i, j,k +
n
E x |i,n+1 (3.51)
j,k
n+1/2 n+1/2
t Hy |i, j+1/2,k − Hy |i, j−1,k
= E x |i,n j,k + (3.52)
z
n+1/2 n+1/2
Hz |i, j,k+1/2 − Hz |i, j,k−1/2
− (3.53)
y
and
where superscripts and subscripts on the right side of the vertical bars (|) are, respec-
tively, temporal indices and spatial indices. FDTD uses a special finite difference grid
designated as a Yee cell [18], as shown in Fig. 3.8. The magnetic field grid is offset
by a half grid size from the electric field grid. The 1/2 term in the spatial indices
indicates that the value is on the half offset grid. These time dependent equations
are updated with an explicit method called a leap-frog algorithm, which updates one
36 3 Governing Equations for Electromechanical Systems
H grid
Hx
E grid
Ex Hy
Ey Hz
z
y
Ez
x
Fig. 3.8 Schematic of a Yee cell. The solid lines denote a cell on an electric field grid and the
dashed lines indicate a cell on a magnetic field grid
field (the electric field) and the other field (the magnetic field) alternatively, and 1/2
in the temporal indices represents the temporal offset of the algorithm. In a similar
manner, the equations for the y and z-directions are also derived.
In the finite element method, there are several options in terms of the main state
variables to be solved. For example, by following the manner in the magnetic case,
it is possible to use the combination of the magnetic vector potential, A, and scalar
electric potential, φ. Another popular option is to use either the electric field or the
magnetic field. By substitution of Eqs. (3.43), (3.46) (3.47), and (3.48) into Eq. (3.42),
the vectorial wave equation in the frequency domain, that is the vectorial Helmholtz
equation, is derived.
1
◦× ◦ × E − ω2
E = − jωJ (3.57)
μ
1
◦× ◦ × H − ω2 μH = ◦ × 1
J , (3.58)
Fig. 3.9 Analysis settings for a closed high frequency electromagnetic system similar to the inside
of a metallic cavity (such as a microwave oven) or a metallic waveguide
Fig. 3.10 Analysis settings for a high frequency electromagnetic system with an open boundary
References
The field of structural optimization can trace its roots back to Anthony Michell’s
influential 1904 paper [29] that was focused on determining which frame structures
were the most economical in withstanding a balanced system of forces. However,
despite Michell’s early work, the pace of research in the field of structural optimiza-
tion did not gain speed until the 1960s when the use of electronic computers became
more common. As reviewed in [45], research related to mathematical programming
based on sensitivity analysis for structural optimization started to emerge during this
time period. By the late 1980s and early 1990s, the field of structural optimization
had evolved from a limited focus on shape and size optimization to include topol-
ogy optimization, as evidenced by early pioneering work [3, 45]. Since then, many
applications and numerical strategies for structural topology optimization have been
developed, as extensively reviewed in [5, 10, 35, 41].
As mentioned in Chap. 1 of this book, the field of structural topology optimization
has grown exponentially in the last two decades to encompass a wide variety of
technical topics. Relevant target applications for this technology include structural
dynamics [18], compliant mechanisms [30, 48], composite materials [17, 36, 37],
sensors and actuators [21, 23, 24, 38, 39], fluid flow structures [2, 6, 31], microfluidic
mixers [1], phononic materials [40], and fluid-structure interactions [27, 49]. More
recently, research studies related to fully coupled multiphysics systems have also
emerged including complex problems involving conjugate heat transfer [12, 14, 26,
50], piezoelectric and fluid-structure interactions [28], and even as many as four
physical processes [51].
While the reader is referred to the literature for numerous examples, the goal of
this section is to review the general form of the topology optimization problems con-
sidered throughout this book. The majority of the structural topology optimization
examples presented in this text utilize the solid isotropic material with penalization
(SIMP) approach [4, 5] for interpolating the design density variables, although struc-
tural optimization using a level set function approach [46] is also briefly reviewed.
Prior to the discussion of the level set method, two separate single physics optimiza-
tion examples for heat conduction and fluid flow are provided in this section as a
precursor to the more complex multiphysics studies introduced in Chap. 5. Interest-
ingly, these two single physics systems are used as the foundation for several of the
multiphysics investigations that follow.
Optimization Problem Description and Computational Process Overview
In the structural topology optimization design examples presented throughout this
book a non-linear optimization problem of the standard form described in [43] is
considered,
Find γ
Minimize Fo (γ )
Subject to gi (γ ) ≤ 0, for i = 1, . . . , m
γ j_min ≤ γ j ≤ γ j_max , for j = 1, . . . , n,
1 Reprinted from Ref. [15], Copyright (2012), with permission from Inderscience Enterprises
Limited.
4.1 Topology Optimization 43
No
Converged?
Yes
Analyze results
(Post-processing)
any number of design variables using only one additional numerical analysis, i.e.,
the adjoint analysis.
For simplicity, the following discussion is posed in terms of the static structural
system first introduced in Chap. 3. By discretizing the governing equation with some
numerical analysis method such as FEM, the following linear algebra equation is
obtained
K(γ )U = F, (4.1)
where K(γ ), U, and F are the design variable embedded system matrix, state variable
vector (or field vector), and excitation vector, respectively. Note that, for compact-
ness, the subscript, S , indicating a structural problem in Eq. (3.11) has been dropped
in Eq. (4.1).
By introducing the adjoint variable vector, λ̂, the objective function is aug-
mented as
Fo∗ = Fo − λ̂T (K(γ )U − F). (4.2)
σU
By choosing λ̂ to satisfy the following equation, Eq. (4.4), the term σγ can be
eliminated from Eq. (4.3)
σFo
− λ̂T K(γ ) = 0. (4.4)
σU
Once the values of the state variables, U, and adjoint variables, λ̂, are obtained, they
may be substituted into the expression for the design sensitivity,
σFo ∗ σK(γ )
= −λ̂T U , (4.5)
σγ σγ
In Eq. (4.6), T is the temperature state variable, Q is the internal heat generation,
and k(γ ) is the material thermal conductivity, which is a function of the design density
variable, γ . The thermal conductivity is then interpolated using a SIMP approach
1 σ T σ2T
Fo = k(γ ) + dρd , (4.8)
2 σx2 σ y2
ρd
Find γ
Minimize Eq. (4.8)
Subject to Eq. (4.6)
ρd γ dρd − vu ≤ 0
0≤γ ≤1
Given Eq. (4.7),
Adiabatic BC
Thus, the optimization process illustrated in Fig. 4.1 was carried out for this exam-
ple problem using this commercial software [9], and AVM was adopted in computing
the sensitivities directly in the program. The MMA algorithm was additionally cou-
pled to the code via a custom MATLAB®3 [25] script to realize the complete flow
of computations (including additional steps beyond traditional FEA), as described
in Fig. 4.1.
The assumed 2-D single physics heat conduction model geometry with boundary
conditions is shown in Fig. 4.2, where the design domain is 1 m square. A fixed
(zero temperature) condition was enforced at the middle of the left edge to represent
a heat sink. The remaining edges of the domain were considered adiabatic (i.e., zero
heat transfer). Constant uniform heat generation (in the form of a 1 W/m2 heat flux
in 2-D) was applied uniformly to the domain. This model was utilized since it is a
standard benchmark in the field of structural topology optimization.
The design domain shown in Fig. 4.2 was meshed using approximately 10,000
quadrilateral elements with 41,000 degrees-of-freedom (DOF). A 40 % solid material
volume constraint (vu = 0.4 m2 ) was used to obtain the optimal topology shown on
the left in Fig. 4.3. Observe that a branching type of topology is constructed which
minimizes the conductive thermal resistance between the heat source and sink.
The single physics heat conduction problem is easily extended to 3-D, and a
related optimization study was undertaken in [11]. The model geometry and optimal
structural topology results for a 40 % solid material volume constraint are shown
in Fig. 4.4, for completeness, where yellow colored regions indicate solid mater-
ial. These results are consistent with those for pure heat conduction found in the
literature [7]. However, note that a relatively coarse mesh with 8,000 hexahedral ele-
ments and approximately 69,000 DOF was used to generate the structure shown in
Fig. 4.4 to avoid incurring a significant penalty in terms of increased computational
cost. Thus, the expansion of the optimization approach into the third dimension can
quickly lead to massively large-scale systems of equations that must be solved, and
this numerical problem is further exacerbated by the inclusion of additional physics.
At this point, it is appropriate to briefly discuss the relationship between mesh size
and the structural designs that are obtained via the topology optimization process
Fig. 4.3 Topology optimization results for the single physics pure heat conduction problem. The
optimal topology shown on the left was obtained using a 40 % solid material volume constraint. Note
dark regions = solid (i.e., thermally conductive) material; light regions = void (i.e., thermally insu-
lative) material. The image on the right shows the computed normalized temperature distribution,
where red regions indicate a higher temperature
Fig. 4.4 3-D model geometry (left) and topology optimization results (right) for the single physics
pure heat conduction problem [11]. The optimal topology was obtained using a 40 % solid material
volume constraint
described thus far. The issue of mesh dependency is explored extensively in the
literature [5], and the reader is referred to it for a deep exploration of the topic.
However, the simplicity of the 2-D pure heat conduction problem allows for efficient
investigation of the basic concept. Generally, finer meshes produce structural designs
that exhibit finer features, and this trend continues without limit in the absence of
some form of filtering.
Here, a useful Helmholtz PDE filtering routine is adopted [20, 22] to illustrate the
manner in which mesh independency and a minimum structural length scale can be
efficiently enforced. Specifically, for pure heat conduction, the optimization problem
is re-formulated as follows
48 4 Optimization Methods for Electromechanical Systems
Find γ
Minimize Fo in Eq. (4.8)
Subject to Eq. (4.6)
γ = H̃(ν̃)
−R 2f ∇ 2 ν̃ + ν̃ = ν
−1
≤ν ≤1
ρd γ dρd − vu ≤ 0
Given Eq. (4.7),
where R f is the filter radius for Helmholtz filtering and H̃(ν̃) is the regularized
Heaviside function of the filtered scalar function, ν̃, with a transition bandwidth
of h t . Note that the design variable has now become ν, and the original design
variable, γ , is now a function of the intermediate projected field, ν̃. Additionally,
the regularized Heaviside function is defined as follows
⎧
⎪
⎪0 (ν̃ < −h t )
⎨ 3 5
H̃(ν̃) = ν̃ ν̃ ν̃
1
+ 15
− 5
+ 3
(−h t ≤ ν̃ ≤ h t ) , (4.9)
⎪
⎪ 2 16 ht 8 ht 16 ht
⎩
1 (h t < ν̃)
Fig. 4.5 Topology optimization results for the single physics pure heat conduction problem
obtained using ∼10,000 finite elements, a solid volume of 40 %, and a Helmholtz PDE filter-
ing routine. The results for a filter radius value, R f , of 2 × 10−5 , 5 × 10−3 , and 2 × 10−2 with a
Heaviside function bandwidth of h f = 1 are shown, respectively, at left, center, and right
convective heat transfer frequently plays an important role in the design of many
electromechanical devices and components. A method for minimizing the dissipated
power in a fluid [6] passing through a flow structure or fluid manifold is adopted.
The implementation of the method is based on [31], where fluid flow in a porous
medium is assumed, and the final design of a manifold is determined by interpo-
lating the inverse permeability between semi-permeable solid and fluid states. Note
that laminar (i.e., deterministic) flow fields with Reynolds (Re) number values below
2100 are assumed in adoption of the method presented below.
The governing equations for steady-state Navier–Stokes laminar fluid flow [34]
in the absence of body forces are
∇ · v = 0, and (4.10)
Ω (v · ∇v) = −∇ P + ∇ · ∂ ∇v + (∇v)T − ς̃(γ )v, (4.11)
where the fluid density and dynamic viscosity are, respectively, Ω and ∂, and the
inverse permeability of the porous medium is ς̃(γ ), which is assumed to be approx-
imately valid for an actual porous medium [31]. The state variables in Eqs. (4.10)
and (4.11) include the fluid pressure, P, and velocity vector, v.
The effective inverse permeability, ς̃(γ ), in Eq. (4.11) is a function that is inter-
polated to define either the solid or fluid regions in the design domain, see [31],
where
q(1 − γ )
ς̃(γ ) = ς̃min + (ς̃max − ς̃min ) . (4.12)
q +γ
In this convex interpolation scheme, q is a tuning parameter set equal to 0.1, while
ς̃min and ς̃max are the minimum and maximum values of the inverse permeability,
respectively.
50 4 Optimization Methods for Electromechanical Systems
6X Outlet Walls
(Uniform out flow) (No slip)
Following [31], a useful objective function for minimizing the fluid power dissi-
pated (i.e., flow resistance) in an assigned design domain is
⎡ ⎤
1 σvi
σv j 2
Fo = ⎣ ∂ + + ς̃(γ )vi2 ⎦dρd , (4.13)
2 σx j σ xi
ρd i, j i
Find γ
Minimize Eq. (4.13)
Subject to (4.10 − 4.11)
Eqs.
γ dρd − vu ≤ 0
ρd
0≤γ ≤1
Given Eq. (4.12).
A simplified 2-D design domain (of size 65 × 40 mm) was assumed for the
example single physics fluid flow design optimization problem, as shown in Fig. 4.6.
A zero pressure side fluid inlet boundary condition was set on the right side of the
design domain, while six circular fluid outlets were distributed uniformly throughout
the domain. The flow rate out of each nozzle was assigned the same value in an effort
to achieve a relatively uniform flow distribution as opposed to simply assigning a
fluid inlet velocity or flow rate condition. A no-slip condition was enforced on the
remaining side wall boundaries. The design domain was meshed using ∼4,400 trian-
gular finite elements for a total of approximately 29,200 DOF including the x and y
components of the fluid velocity vector plus the pressure state and design optimiza-
tion variables. For reference, greater details about a related numerical example may
be found in [15].
The optimal fluid manifold topology that balances the fluid flow rate out of all six
nozzles (with Re ≈ 10) is shown in Fig. 4.7 along with the normalized fluid velocity
4.1 Topology Optimization 51
Fig. 4.7 Topology optimization results for the single physics fluid flow manifold design problem.
The optimal topology was obtained using a 50 % solid material volume constraint. Note dark regions
= solid; light regions = fluid. The image includes the superimposed fluid velocity contours with
larger fluid velocities shown in red at higher elevations
Fig. 4.8 Fabricated manifold prototype based on topology optimization results for the single
physics fluid flow design problem. Both the manifold (shown on bottom) and cap (shown on top)
were machined out of aluminum material
contours. Note that larger fluid velocities are shown in red at higher elevations and
a 50 % solid material volume constraint was used. The curvy side wall manifold has
been compared with traditional manifold designs and was shown to provide reduced
flow resistance (pressure drop) and greater outlet flow uniformity, as expected [15].
Additionally, a related prototype structure (designed based on the optimization con-
52 4 Optimization Methods for Electromechanical Systems
cept) has been fabricated, refer to Fig. 4.8, to further verify flow performance in an
electronics cooling application [13, 16].
The above design optimization method may again be readily extended to 3-D, and
further discussion of this is provided in Chap. 5 in relation to the design optimization
of thermal-fluid (i.e., conjugate heat transfer) multiphysics systems.
The level set method is a technique to represent the surface of a structure implicitly
with a scalar field or so-called level set function. The method was originally proposed
by Osher and Sethian [32], and it has been applied to various problems that involve
time dependent varying free surfaces, such as two-phase flow problems. Also, the
method may be applied to various static problems and structural optimization is one
of them, that is, what we call level set-based topology optimization. Hereafter, the
term, level set method or LSM, refers to level set-based topology optimization in this
section. While both topology optimization and the level set method were proposed
independent of each other, both techniques have a large part in common.
Both methods use a scalar field, the density field in topology optimization and
the level set function in the level set method, to represent the designed structure, and
the field value is updated after each numerical iteration based on design sensitivity
information. Usually, the level set function is updated by solving a time-dependent
equation, a kind of Hamilton–Jacobi equation or so-called level set equation, while
topology optimization usually uses mathematical programming techniques.
The most significant difference is interpretation of the scalar field. Topology opti-
mization directly interprets the field distribution as the shape and the topology of
the structure, while the level set method extracts the iso-contour, usually at a zero
value (i.e., the zero-level set) as the surface of the structure thus regarding one side
of the contour as the material domain. Theoretically, this means that there exists no
intermediate material state in between air and solid (unlike topology optimization),
and thus, the level set method is free of the so-called gray scale problem.
Figure 4.9 illustrates the overall idea of the level set method. The level set function,
which is shown in the colored height map in the figure, is updated as time progresses.
There is a black line on the level set function which indicates the iso-contour at a
certain value such as zero. At the bottom, the iso-contour is projected on the plane
and inside of the contour is filled to indicate that material domain which is defined
by the level set function. As the level set function is updated, the material domain
also is updated. Sometimes, this update process involves topological changes, such
as splitting or merging of the domains or dissipation of holes.
Algorithm of the Level Set Method
Hereafter, a brief introduction of a level set-based topology optimization method
described in a paper by Yamasaki et al. [47] is given.
4.1 Topology Optimization 53
Material domain
Time
Fig. 4.9 Time-dependent update of level set function and represented shape
Let ρd represent the fixed design domain, which entirely envelopes the structure
to be designed. The designed structure is filled with a material, so that it is equivalent
to the material domain, ρm , and the complementary domain is the void domain.
The interface of the two domains is defined by the level set contour of the level set
function, φ(x)
⎧
⎪
⎨φ(x) > 0 for ∀x ∈ ρm \ σρm ,
φ(x) = 0 for ∀x ∈ σρm , (4.14)
⎪
⎩
φ(x) < 0 for ∀x ∈ ρd \ ρm ,
where x stands for a position in ρd , and σρm is the boundary of ρm . By using the
above level set function, an arbitrary shape and topology of the material domain,
ρm , in the fixed design domain, ρd , can be represented.
The level set function is then projected to material density distribution by the
following Heaviside function,
0 (φ < 0)
H(φ) = (4.15)
1 (φ ≥ 0).
Minimize Fo := l(u)dρd (4.16a)
φ ρd
Subject to g1 := H̃(φ)dρd − vu ≤ 0 (4.16b)
ρd
g2 := α(φ(x))dρm − P̄ < 0, (4.16c)
ρm
where l is the compliance of the structure, u is the physical state variable, vu is upper
bound of material volume, and P̄ is upper bound for the perimeter of the structure
for regularization of the complexity of the shape (i.e., so-called perimeter control).
With the classic formulation, the level set function is updated with the following
level set equation
σφ (x, t)
+ v N (x, t) |∇φ(x, t)| = 0, (4.17)
σt
where t is fictitious time and v N (x, t) is the normal velocity provided by sensitivity
analysis. Solving Eq. (4.17) is not a trivial task since it may contain singular
points. Therefore, Yamasaki et al. [46] proposed a technique called geometrical
re-initialization, which extracts the level set contour and calculates the signed dis-
tance from each node to the level set contour geometrically. With this method,
|∇φ(x, t)| = 1 is always maintained, and thus, Eq. (4.17) is reduced to a simpler
form
σφ (x, t)
+ v N (x, t) = 0. (4.18)
σt
Then, this expression is discretized with respect to t using forward finite difference
φ(x, t + ηt) − φ(x, t)
+ v N (x, t) = 0, (4.19)
ηt
and finally, this yields the following recursive update equation
Since the level set method itself does not have a framework to handle constraints,
the Lagrangian is built with the objective function and constraints
L(φ) = Fo (φ) + λ̂i gi (φ), (4.21)
i
4.1 Topology Optimization 55
10 cm
where λ̂i is the Lagrange multiplier for the constraint, gi . A natural extension of
this method is to use the augmented Lagrangian method for update of the Lagrange
multipliers.
The normal velocity is calculated using design sensitivity information of the
Lagrangian
σL (φ)
vN = (4.22)
σφ
where ηx is the maximum distance of adjacent mesh nodes; refer to [47] for addi-
tional details.
Level Set Method – Structural Example
Consider a plane-stress short cantilever benchmark problem to minimize compliance
to the boundary condition depicted in Fig. 4.10. The BC, σρu , is fixed and a 10
N load is vertically applied downward on σρt . The elastic modulus of the solid
portion of the material domain is E o = 2.1 × 1011 N/m2 and the Poisson’s ratio is
ν = 0.3. The allowed volume fraction of the material domain is 0.5. The level set
method cannot be started with a completely flat structure. So, a trivial initial design
is given, as shown in Figs. 4.11 and 4.12 shows the optimal configuration obtained
with the level set method [47]. As shown in the figure, the method yields a reasonable
black and white result.
In the above section of this chapter, an emphasis was placed on the topology optimiza-
tion of distributed parameter (or continuum) systems, which is made mathematically
feasible through the breakdown of the continuum into many discrete finite elements
56 4 Optimization Methods for Electromechanical Systems
Fig. 4.11 Provided structure (for the level set method algorithm). With kind permission from
Springer Science+Business Media: Structural and Multidisciplinary Optimization, A level set-based
topology optimization method using the discretized signed distance function as the design variables,
Yamasaki et al. [47, Fig. 4], © Springer-Verlag 2009
Fig. 4.12 Optimal structure of the example for the 2-D minimum compliance problem. Perimeter
constraint is imposed, initial structure is Fig. 4.11 and number of elements = 100 × 50. With kind
permission from Springer Science+Business Media: Structural and Multidisciplinary Optimization,
A level set-based topology optimization method using the discretized signed distance function as
the design variables, Yamasaki et al. [47, Fig. 7b], © Springer-Verlag 2009
Ls Ls Ls
xs ts
u
s =3 s =2 s =1 xs
Fig. 4.13 A cantilevered beam structure consisting of three segments each having a hollow square
cross section. The beam consists of m segments each of length L s . Each segment has a thickness,
ts , with a side length of xs for segment s = 1, . . . , m for the design variables, as shown in the beam
end view image on the right
main focus in this book is on the application of automated parametric sweeps to the
informed design and structural optimization of 3-D multiphysics systems involving
multiple discrete variables (related to the physical size of components) within a finite
element framework. While such tools do not necessarily employ gradient-based opti-
mizers to arrive at a globally optimal structural configuration, they allow for rapid
visualization and performance comparison of a range of possible designs that are
selectable based on one or more other determining factors such as global manufac-
turability, cost, or packaging constraints. These tools are shown to complement the
structural topology optimization techniques described earlier in this chapter. Two
representative multiphysics parametric sizing analysis/optimization studies are pre-
sented in Chap. 5 in relation to the composite laminate layer thickness determination
of a thermal-structural substrate assembly for an electronics packaging application
and the selection of a specific channel aspect ratio for an electronics cold plate
assembly.
References
1. Andreasen CS, Gersborg AR, Sigmund O (2008) Topology optimization for microfluidic mix-
ers. Int J Numer Meth Fl 61:498–513. doi:10.1002/fld.1964
2. Aage N, Poulsen TH, Gersborg-Hansen A, Sigmund O (2008) Topology optimization
of large scale Stokes flow problems. Struct Multidiscip O 35:175–180. doi:10.1007/s00158-
007-0128-0
3. Bendsøe MP, Kikuchi N (1988) Generating optimal topologies in structural design using
a homogenization method. Comput Method Appl M 71:197–224. doi:10.1016/0045-
7825(88)90086-2
4. Bendsøe MP, Sigmund O (1999) Material interpolation schemes in topology optimization.
Arch Appl Mech 69:635–654. doi:10.1007/s004190050248
5. Bendsøe MP, Sigmund O (2003) Topology optimization: theory, methods, and applications,
2nd edn. Springer, Berlin
6. Borvall T, Petersson J (2003) Topology optimization of fluids in Stokes flow. Int J Numer Meth
Fl 41:77–107. doi:10.1002/fld.426
7. Chen Y, Zhou S, Li Q (2010) Multiobjective topology optimization for finite periodic structures.
Comput Struct 88:806–811. doi:10.1016/j.compstruc.2009.10.003
58 4 Optimization Methods for Electromechanical Systems
30. Nishiwaki S, Frecker MI, Min S, Kikuchi N (1998) Topology optimization of compliant mech-
anisms using the homogenization method. Int J Numer Meth Eng 42:535–559. doi:10.1002/
(SICI)1097-0207(19980615)42:3<535:AID-NME372>3.0.CO;2-J
31. Olesen LH, Okkels F, Bruus H (2006) A high-level programming-language implementation
of topology optimization applied to steady-state Navier-Stokes flow. Int J Numer Meth Eng
65:975–1001. doi:10.1002/nme.1468
32. Osher S, Sethian JA (1988) Fronts propagating with curvature-dependent speed: algo-
rithms based on Hamilton-Jacobi formulations. J Comput Phys 79:12–49. doi:10.1016/0021-
9991(88)90002-2
33. Qian X, Sigmund O (2013) Topological design of electromechanical actuators with robustness
toward over- and under-etching. Comput Method Appl M 253:237–251. doi:10.1016/j.cma.
2012.08.020
34. Reddy JN, Gartling DK (2000) The finite element method in heat transfer and fluid dynamics,
2nd edn. CRC Press, Boca Raton
35. Rozvany GIN (2009) A critical review of established methods of structural topology optimiza-
tion. Struct Multidiscip O 37:217–237. doi:10.1007/s00158-007-0217-0
36. Sigmund O, Torquato S (1996) Composites with extremal thermal expansion coefficients. Appl
Phys Lett 69:3203. doi:10.1063/1.117961
37. Sigmund O, Torquato S (1999) Design of smart composite materials using topology optimiza-
tion. Smart Mater Struct 8. doi:10.1088/0964-1726/8/3/308
38. Sigmund O (2001) Design of multiphysics actuators using topology optimization - Part I:
one-material structures. Comput Method Appl M 190:6577–6604. doi:10.1016/S0045-
7825(01)00251-1
39. Sigmund O (2001) Design of multiphysics actuators using topology optimization - Part II:
two-material structures. Comput Method Appl M 190:6605–6627. doi:10.1016/S0045-
7825(01)00252-3
40. Sigmund O, Jensen JS (2003) Systematic design of phononic band-gap materials and structures
by topology optimization. Philos T Roy Soc A 361:1001–1019. doi:10.1098/rsta.2003.1177
41. Sigmund O, Maute K (2013) Topology optimization approaches. Struct Multidiscip O 48:1031–
1055. doi:10.1007/s00158-013-0978-6
42. Spillers WR, MacBain KM (2009) Structural optimization. Springer, Dordrecht
43. Svanberg K (1987) The method of moving asymptotes—a new method for structural optimiza-
tion. Int J Numer Meth Eng 24:359–373. doi:10.1002/nme.1620240207
44. Svanberg K, Svärd H (2013) Density filters for topology optimization based on the Pythagorean
means. Struct Multidiscip O 48:859–875. doi:10.1007/s00158-013-0938-1
45. Suzuki K, Kikuchi N (1991) A homogenization method for shape and topology optimization.
Comput Method Appl M 93:291–318. doi:10.1016/0045-7825(91)90245-2
46. Yamasaki S, Nishiwaki S, Yamada T, Izui K, Yoshimura M (2010) A structural optimization
method based on the level set method using a new geometry-based re-initialization scheme.
Int J Numer Meth Eng 83:1580–1624. doi:10.1002/nme.2874
47. Yamasaki S, Nomura T, Kawamoto A, Sato K, Izui K, Nishiwaki S (2010) A level set based
topology optimization method using the discretized signed distance function as the design
variables. Struct Multidiscip O 41:685–698. doi:10.1007/s00158-009-0453-6
48. Yin L, Ananthasuresh GK (2002) A novel topology design scheme for the multi-physics prob-
lems of electro-thermally actuated compliant micromechanisms. Sensor Actuat A-Phys 97–
98:599–609. doi:10.1016/S0924-4247(01)00853-6
49. Yoon GH, Jensen JS, Sigmund O (2007) Topology optimization of acoustic-structure interaction
problems using a mixed finite element formulation. Int J Numer Meth Eng 70:1049–1075.
doi:10.1002/nme.1900
50. Yoon GH (2010) Topological design of heat dissipating structure with forced convective heat
transfer. J Mech Sci Technol 24:1225–1233. doi:10.1007/s12206-010-0328-1
51. Yoon GH (2012) Topological layout design of electro-fluid-thermal-compliant actuator. Com-
put Method Appl M 209–212:28–44. doi:10.1016/j.cma.2011.11.005
Chapter 5
Electromechanical System Simulation
and Optimization Studies
Thermal management of electronic systems is crucial for their reliable operation and
efficient design tools for heat sinks, cold plates, bus bars, and packaging materials
are increasingly becoming important as assembly size decreases and power density
increases. Accordingly, five design optimization examples are covered below with the
first topic centered on exploring the trade-offs between minimizing electrical resis-
tance versus thermal resistance in the multiobjective design optimization of electrical
conductors. Thermomechanical systems are then addressed in the next section, where
a direct-bonded copper (DBC) substrate composite is designed for increased relia-
bility by minimizing the maximum stress present at the bonded interface between the
substrate and an attached electronics device during high temperature operation. From
there, the focus shifts to traditional conjugate heat transfer thermal-fluid systems and
applications to cold plate design, and building on this, a more complex three-physics
fully coupled system is subsequently proposed, which involves the design layout
of a heat sink for a localized heat source considering thermomagnetic convective
heat transfer. In the final part of this section, heat flow control in structures that
exhibit anisotropic material physical properties (specifically, thermal conductivity)
is introduced to illustrate how structural optimization methods may be applied at
multiple physical length scales. From a broad perspective, the different simulation
and optimization techniques presented in this section represent a suite of numerical
tools for enhanced electro-thermomagnetic system concept development and design.
5.1 Electronic System Component Analysis and Design 63
Logically, many of these tools may be synthesized and/or adapted to other related
engineering applications.
λT
γC p = ◦ · k (σ ) ◦T − h(σ )(T − To ) + Q, (5.1)
λt
◦ · J = Q j, (5.2)
J = ρ (σ )E + Je , (5.3)
E = −◦ν, (5.4)
where the state variables include the conductor temperature, T , current density, J,
and electric scalar potential, ν. In Eqs. (5.1) and (5.3), the conductor-specific heat
capacity, thermal conductivity, and electrical conductivity are given as C p , k(σ ), and
64 5 Electromechanical System Simulation and Optimization Studies
ρ (σ ), where the latter two material physical parameters are defined as a function of
the design density variable, σ . The heat source term, Q, in Eq. (5.1) arises from the
electric current determined via Eqs. (5.2)–(5.4), where Je and E are the externally
applied current and electric field strength, respectively. For steady-state analysis (as
considered below), the term on the left-hand side of Eq. (5.1) disappears.
In the bus bar design optimization problem, both the thermal conductivity and
electrical conductivity material physical parameters are interpolated using a SIMP
approach to determine the design of the conductor as follows:
where ko and ρo are the thermal conductivity and electrical conductivity of the
assumed solid conductor material, and p is a standard penalization parameter set
equal to three following [6].
Since bus bars are typically not completely covered by a thermal/electrical insu-
lation layer, convection to the environment may occur through the external surface
of the conductor, and this effect plays a role in heat transfer during operation. In
the design optimization example that follows, this convective heat transfer effect is
captured by enforcing a surface convection heat transfer boundary condition, refer
to Sect. 3.2, at the solid-to-void interface of the conductor. Here, an approach similar
to [44] is adopted, where the design-dependent heat transfer coefficient, h(σ ), is
defined via a smoothed hat-function as
h × 10−2
if σ < σl − Ωb
o 1 σ −σ
∂(σ −σl )
h o 2 + 2Ωb + 2∂ sin if σl − Ωb < σ < σl + Ωb
1
l
Ωb
h(σ ) = h o
if σl + Ωb < σ < σu − Ωb (5.7)
h o 1 − σ −σu − 1 sin ∂(σ −σu )
if σu − Ωb < σ < σu + Ωb
2 2Ωb 2∂ Ωb
0 if σu + Ωb < σ
Fo = w1 A + w2 B, where (5.8)
⎧ ⎪ ⎨ 2 ⎩
1 λ T λ2T
A= k(σ ) + dςd , and (5.9)
2 λx2 λ y2
ςd
⎧ ⎪ ⎨ 2 ⎩
1 λ ν λ 2ν
B= ρ (σ ) + dςd . (5.10)
2 λx2 λ y2
ςd
A general energy formulation is used to obtain both terms in Eq. (5.8), where A
is related to the thermal resistance of the system, while B is associated with the elec-
trical resistance of the conductor. Similar to the thermal-fluid design optimization
in Sect. 5.1.3 that follows later in this chapter, w1 and w2 in Eq. (5.8) are weight-
ing values that scale the respective thermal and electrical portions of the objective
function.
Thus, the full electrothermal optimization problem may be stated as
Find σ
Minimize Eq. (5.8)
Subject to Eqs. (5.1) − (5.4)
σ = H̃(φ̃)
−R 2f ◦ 2 φ̃ + φ̃ = φ
−1
∗φ ∗1
σ dςd − vu ∗ 0
ςd
Given Eqs. (5.5) − (5.7).
x Terminal 3
(Ground)
The design domain was meshed with approximately 8,790 triangular finite ele-
ments having √35,600 DOF, and the optimization problem was solved using com-
mercial finite element solvers [15] coupled with a MMA optimizer [90]. Values of
σl = 0.45, σu = 0.75, and Ωb = 5 × 10−3 were used in Eq. (5.7). A Helmholtz
filtering routine with a filter radius of R f = 5 × 10−3 was also implemented to
enforce a minimum design length scale, as described in Sect. 4.1.
As an initial trial, weighting values of w1 = 0 and w2 = 1 in Eqs. (5.8)–(5.10)
were assumed along with a 20 % conductor volume fraction constraint (vu = 0.2 m2 )
to solve the pure electrical resistance minimization problem. The optimal solution
was found in 50 iterations and took less than 5 minutes of computational time on a
two-core laptop with 2.50 GHz processors and 3 GB of RAM. The bus bar design
solution is shown on the left in Fig. 5.2 including normalized temperature contours
and normalized total current density vectors. The surface convection heat transfer
coefficient distribution is shown on the right in Fig. 5.2. Observe that a bus bar with
the minimum length possible is constructed between the current source (Terminal 1)
and the nearest ground (Terminal 2) to minimize the electrical resistance of the
structure. Furthermore, Fig. 5.2 highlights the manner in which the surface heat
transfer coefficient follows the conductor surface boundary with a zero value inside
the solid material and a two-order of magnitude lower value outside of the bus bar,
as specified per Eq. (5.7).
Extending the above results, two additional objective function weighting value
pairs were examined to visualize the trade-off between the electrical and thermal
terms in Fo , Eqs. (5.8)–(5.10). The bus bar design results for w1 and w2 pair values
corresponding to Case I = [0, 1], Case II = [0.35, 0.65], and Case III = [0.5, 0.5]
are shown in Fig. 5.3 on the left, center, and right, respectively.
The calculated bus bar electrical resistance, R, in terms of the voltage drop, V ,
and electrical current, i, is
5.1 Electronic System Component Analysis and Design 67
Fig. 5.2 Optimal bus bar design for 20 % conductor volume fraction and minimum electrical
resistance, [w1 = 0, w2 = 1], with normalized temperature contours and total current density
vectors (red arrows) shown on the left. Note the red colored temperature contours indicate the
highest temperature; dark regions = solid material; light regions = void. The corresponding final
heat transfer coefficient distribution on the domain is shown on the right in units of W/(m2 K)
Fig. 5.3 Optimal bus bar designs for 20 % conductor volume fraction and weighting value (i.e., w1
and w2 ) pairs of [0, 1], [0.35, 0.65], and [0.5, 0.5] are shown, respectively, on the left, center, and
right with normalized temperature contours and total current density vectors (red arrows). Note the
red colored temperature contours indicate the highest temperature; dark regions = solid material;
light regions = void
V
R= , (5.11)
i
while convective thermal resistance, Rth(cnv) per [45], is expressed in terms of the
applied heat transfer coefficient, h, and conductor surface area, As , as
1
Rth(cnv) = . (5.12)
h As
These two metrics are provided in Table 5.1 for each design shown in Fig. 5.3. Con-
vective thermal resistance is used as a metric here since the sole heat loss mechanism
68 5 Electromechanical System Simulation and Optimization Studies
Table 5.1 Electrical and convective thermal resistances for Cases I–III shown in Fig. 5.3 (normal-
ized by the respective Case I values)
Case I Case II Case III
common in today’s hybrid vehicles, they are also found in a variety of other sustain-
able engineering applications that involve energy generation and storage including
wind turbines, solar energy installations, and future wireless charging systems. In
all of these cases, multiphysics topology optimization has the potential to play a key
role in performance enhancement.
Substrate damage is a typical end result of high-temperature-induced stress found
in power electronics packaging. Direct bond copper substrates have been used in
recent power module designs due to superior thermal performance. A DBC substrate
consists of an insulating ceramic layer sandwiched between two copper layers. Often,
additional materials may be deposited on top of the copper layers. The different
materials have distinct thermal properties including different CTEs. Copper has a
CTE that is 3–5 times larger than those of ceramic materials found in power module
substrates. The CTE mismatch results in significant thermally induced stress at the
interfaces between the electronics device and the substrate. Therefore, for reliable
power modules, it is critical to minimize the thermally induced stress in the die and
the corresponding bonding layer. Various efforts related to fabrication and material
design have been devoted to solving such thermal stress problems. Stress reduction
may be achieved by parametrically varying the thickness of the ceramic or metal
layers. However, beyond simple material changes or thickness modifications, a new
approach to thermal stress reduction and delamination prevention is to modify the
overall substrate dimensions (e.g., size) and introducing a step layer at the edges of
the metal layer; see for example [72]. This interesting approach suggests the use of
numerical optimization methods in the resolution of thermal stress and delamination
issues in power module substrates.
In this section, a parametric study of a sample electronics package is presented,
where the thickness of the top and bottom metal layers of the DBC substrate is varied
in order to minimize the stress in the substrate to the device bonding layer. While
bonding layer stress reduction may be achieved via asymmetric DBC metal layer
thickness, this reduction comes at the cost of increased complexity in terms of the
base fabrication of the DBC structure. Instead, level set-based topology optimization
is applied to the design of the top copper layer of a DBC substrate, per [76], having
symmetric metal layer thickness in order to better minimize the thermally induced
stress that occurs at the bonding layer between the electronics device and the top
metallic layer of the DBC. Here, we adopt a 2-D pattern design approach instead
of a cross-sectional configuration design strategy to avoid additional changes to the
fabrication process. Thus, following the method introduced in Sect. 4.1.1, a 2-D
level set field is first prepared for the design variables, and then the level set contour
is projected into a 3-D space to construct a model for finite element analysis. The
optimization code was implemented following the level set topology optimization
approach with geometric reinitialization scheme from [93]. Finally, a 2-D design
is directly obtained, which is utilized as a photo mask pattern for fabrication of a
prototype DBC.
70 5 Electromechanical System Simulation and Optimization Studies
Uniform DBC
Bottom Metal Layer
Symmetry Plane Symmetry Plane (0.3 mm to 0.4 mm)
Fixed Displacement
Convective Heat Flux
Fig. 5.4 One-quarter symmetry solid model geometry of electronics package used for parametric
sizing study
The deformation and resulting thermal stress state of a laminated composite such
as the electronics package described above is highly dependent on the CTE of the
individual layers of the composite structure as well as the thickness of each lamina.
Here, this electronics package is studied initially via a parametric sizing analysis to
explore this stress-to-layer-thickness relationship.
The model geometry utilized in this study is shown in Fig. 5.4, where a one-
quarter symmetry model is adopted. Standard material properties were assumed for
the various material layers, which moving from bottom to top include copper (in-
plane √17 mm square, thickness √0.3 mm), aluminum nitride (in-plane √21 mm
square, thickness √0.6 mm), copper (in-plane √17 mm square, thickness √0.3 mm),
solder (in-plane √9 mm × 12 mm, thickness √0.1 mm), and silicon (in-plane
√9 mm ×12 mm, thickness √0.3 mm). A heat transfer coefficient of h = 2,000
W/(m2 K) at a reference temperature of To = 293 K was applied to the bottom of
the package. A device power of 50 W was applied to the silicon layer.
In this investigation, the top metal layer plus the bottom metal layer of the DBC
are each varied in size from 0.3 to 0.4 mm in 0.01 mm increments to determine the
optimal configuration that minimizes the maximum Von Mises stress in the package
solder bond layer. Thus, a total of 121 simulations were run in parametric fashion
over 9 h using an eight-core workstation with 2 GHz processors and 64 GB of
RAM to determine the response of the composite structure for all possible layer
thickness combinations. The computational mesh settings adopted for the parametric
study involved an approximate total number of DOF ranging from 3.2 ×105 to 4.2
×105 depending on the structural configuration (i.e., layer thickness geometry) under
consideration.
5.1 Electronic System Component Analysis and Design 71
The deformed state of the package is shown in Fig. 5.5, where the DBC top and
bottom metal layers have thicknesses of 0.3 and 0.4 mm, respectively. The structural
response for the 121 different possible layer thickness combinations was computed,
and configuration #11 shown in Fig. 5.5 produces the lowest maximum Von Mises
stress, 52.1 MPa, in the solder bond layer with a maximum device temperature of
118.8 ⊥ C. The response surface from the full parametric study is shown in Fig. 5.6,
where we see that asymmetric package designs comprising a DBC with a top metal
layer that is thinner than the bottom metal layer consistently perform better. Nonethe-
less, such asymmetric DBC designs may pose potential problems in terms of defor-
mation due to residual stresses during high temperature substrate fabrication, and a
package design strategy that retains the symmetry of the DBC metal layer thickness
is desirable.
Optimization Model and Results
From Sect. 3.3, the governing equations for the steady-state thermal stress analysis
of a structure may be written as
− ◦ · k (σ ) ◦T = Q, (5.13)
− ◦ · η = Fα , (5.15)
where T is the temperature state variable, and η is the matrix of thermally induced
stresses. The thermal conductivity, k(σ ), stiffness tensor, C(σ ), and coefficient of
thermal expansion, α(σ ), are all design-dependent material physical parameters.
Note that To in the coupling expression, Eq. (5.14), is a reference temperature required
for the calculation of the resultant thermal strains.
For the thermal stress design optimization problem, the design-dependent material
physical parameters are interpolated as follows:
In all three of the above expressions, the subscript, o, indicates a reference value,
which in this case, corresponds to void versus the subscript, s, which corresponds to
a solid material physical parameter. Note that the stiffness tensor, C(σ ) in Eq. 5.14,
depends on the design-dependent elastic modulus, E(σ ), in Eq. 5.16. Additionally,
H̃(σ ) is the relaxed (differentiable) Heaviside function, as described in Sect. 4.1.
A metal layer pattern for the DBC structure is designed by minimizing the objec-
tive function,
72 5 Electromechanical System Simulation and Optimization Studies
[mm]
[°C]
Fig. 5.5 Representative deformed package geometry with superimposed total displacement
contours (upper image) and temperature contours (lower image) for parametric case #11 with
the DBC top metal layer thickness equal to 0.3 mm and the DBC bottom metal layer thickness set
to 0.4 mm
5.1 Electronic System Component Analysis and Design 73
Maximum
[MPa]
Minimum
Bottom Metal Layer Top Metal Layer
Thickness [mm] Thickness [mm]
Fig. 5.6 Response surface from the parametric size analysis for the DBC with top and bottom metal
layers of varying thickness ranging from 0.3 to 0.4 mm in 0.01 mm increments for a total of 121
possible configurations
⎧Ω
1
Fo = η : dΩ; Ω = ◦u + (◦u)T , (5.19)
2
0
which is defined on the bonding layer between the device and DBC. In Eq. (5.19),
η is the Cauchy stress tensor, Ω is a vector of the resulting strains, and u is the
displacement vector in the domain of interest.
A one-quarter symmetry model of an assumed multilayer electronics package
including the device at top and a laminated DBC structure on bottom with a solder
bond attachment layer positioned in between is shown in Fig. 5.7. Here, the top
metal layer of the DBC structure is specifically designed to minimize the strain
energy density, via Eq. (5.19), in the mechanical bonding layer between the DBC
and device.
The numerical model for the optimization is similar to the one presented above;
that is, the same device power dissipation, heat transfer convective flux, and boundary
conditions were assumed. Additionally, the same material physical parameters were
used for each layer in the numerical model. In this case, the top metal layer was
meshed in 2-D using 1,764 quadrilateral elements, while the subsequent full 3-D
structure was meshed using √ 45 × 103 hexahedral elements. The optimal solution
was found using commercial finite element software [15] in 100 iterations and took
74 5 Electromechanical System Simulation and Optimization Studies
Top view
DBC metal layer Quarter view
(top side, Symmetry
design domain) boundaries
Device
(heat source)
Bottom view
Metal layer
(bottom side,
heat flux boundary)
Bond layer
(extremely thin)
Other surfaces
(adiabatic)
Fig. 5.7 One-quarter symmetry model of DBC structure including silicon device at the top of the
stack, an extremely thin middle solder bond layer, and a laminated DBC structure
Fig. 5.8 Package geometry assumed for initialization of the level set algorithm (left) and final
structure including optimized top metal layer of DBC (right)
Photo resistor
Cu
Insulation
Cu
Fig. 5.9 Prototype DBC structures batch fabricated using the level set optimization result with a
standard DBC etching process, as shown schematically at the top of the figure
[MPa]
Fig. 5.10 Overlaid response surfaces from the parameter studies for the DBC with and without
optimal topology top metal layer
Discussion
As a follow-up to the earlier parametric sizing study, a similar study was performed
using the optimized DBC top metal layer topology, and the response surfaces from
the two studies are plotted together in Fig. 5.10. For the same range of DBC top
76 5 Electromechanical System Simulation and Optimization Studies
and bottom metal layer thicknesses, optimizing the topology of the top metal layer
consistently produces lower overall stress levels. In fact, in all cases where the DBC
has symmetric top and bottom metal layer thickness, the DBC with optimized top
metal layer topology demonstrates enhanced performance. For example, the upper
image in Fig. 5.11 shows the deformed state of the package with optimized DBC top
metal layer topology with superimposed total displacement contours for parametric
case #1. Here, the DBC top metal layer thickness is 0.3 mm, and the DBC bottom
metal layer thickness is 0.3 mm. In this case, the maximum total displacement, 32 μm,
occurs at the corners of the structure with plate bending in two dimensions; note that
this bending mode is consistent across all 121 parametric models. The specific layer
thickness combination for the structure shown in Fig. 5.11 results in a Von Mises
stress of 54.7 MPa in the DBC to device solder bond layer, which is within 5 %
of the best case scenario using a DBC with asymmetric metal layer thickness. The
corresponding thermal contours for this case are also shown in the lower image of
Fig. 5.11, for reference, and the maximum device temperature is 120.2 ⊥ C.
A physical explanation for the enhanced performance of the optimized DBC is that
the finger-like structure introduces greater compliance and establishes an equal stress
state in the top metal layer of the DBC, which relieves stress concentrations in the
solder bond layer. This idea is supported by the larger maximum total displacement
(i.e., 32 vs. 26 μm) of the package plus the overall package Von Mises stress contours,
as shown in Fig. 5.12. Additionally, Fig. 5.13 shows the volumetric strain contours at
the solder bond to device interface surface and highlights lower strain at the corners of
this critical interface in the optimized package (upper image) versus the well-known
stress concentration that typically occurs in corners of regularly shaped rectangular
or square DBC packages (lower image). Note that the color contours in Fig. 5.13 are
plotted using the same range of volumetric strain values (with higher strain levels
shown in red and lower strain levels in blue) for direct comparison of the images.
[mm]
[°C]
Fig. 5.11 Representative deformed package geometry with superimposed total displacement con-
tours (upper image) and temperature contours (lower image) for parametric case #1 with the DBC
top metal layer thickness equal to 0.3 mm and the DBC bottom metal layer thickness set to 0.3 mm
78 5 Electromechanical System Simulation and Optimization Studies
function. The first term is related to the domain average temperature, while the second
term is associated with the fluid power dissipated (or flow resistance) in the domain.
Following the material distribution approach outlined in Chap. 4, the effective inverse
permeability and thermal conductivity of an assumed porous medium are specified
as a function of the material design variable, σ ; see [8, 95] for additional details.
Related methods have previously been used in the design optimization of ducts [77],
heat transfer surfaces [20, 44], and microfluidic mixers [4].
Following [19, 24, 77, 83], and similar to the single physics fluid flow example
from Sect. 4.1, the governing equations for steady-state Navier-Stokes flow in an
idealized porous medium are given as
◦ · v = 0, and (5.20)
γ (v · ◦v) = −◦ P + ◦ · η ◦v + (◦v) T − α̃ (σ ) v. (5.21)
Equation (5.20) represents the fluid incompressibility constraint, while Eq. (5.21)
describes laminar fluid flow. In these expressions, γ and η are again the fluid density
and dynamic viscosity, respectively. The inverse permeability of the porous medium,
α̃(σ ), is assumed to approximately represent an actual porous medium, per [77]. The
state variables once more include the fluid pressure, P, and velocity field terms in
the vector, v.
In addition to fluid flow, the governing equation for steady-state convection-
diffusion heat transfer is
γC p (v · ◦T ) = ◦ · k (σ ) ◦T + Q, (5.22)
where C p represents the heat capacity, k(σ ) is the thermal conductivity of the fluid,
and Q is the volumetric power density.
Following the approach in Sect. 5.1.1, a multi-term objective function, Fo , was
implemented to optimize for both heat transfer and fluid flow. Specifically, the objec-
tive was specified to minimize the mean temperature and total fluid power dissipated
in the system,
Fo = w1 A + w2 B, where (5.23)
⎧ ⎡ ⎣
A= k (σ ) (◦T )2 + γC p [T (v · ◦T )] dςd , and (5.24)
ςd
⎧
1 ⎦ ⎨ λvi ⎩
λv j 2 ⎦
B= ⎤ η + + α̃ (σ ) vi2 dςd . (5.25)
2 λx j λ xi
ςd i, j i
As described in [19, 24], The terms w1 and w2 in Eq. (5.23) are user selected
weighting values that scale the respective thermal and fluid portions of the objective
function. ’Tuning’ of these values assists in convergence and modifies the result-
ing optimal topology by affecting the dominance of one physical process relative
80 5 Electromechanical System Simulation and Optimization Studies
to another. For simplicity, these weighting values were selected manually in the
numerical example that follows, although an adaptive scaling strategy for automatic
determination of the weighting of the individual objective function terms could alter-
natively be adopted for this problem; more information is provided later in Sect. 5.4.1.
Focusing on the optimal steady-state fluid flow and channel layout, the thermal
conductivity and inverse permeability of the porous medium were, respectively, inter-
polated using a penalty method and convex interpolation scheme; refer to Chap. 4.
These effective properties, k and α̃, were interpolated via the main design parameter,
σ , which varied from 0 (minimally porous, lower thermally conductive solid) to 1
(high thermally conductive fluid), respectively, as follows:
q(1 − σ )
α̃ (σ ) = α̃min + (α̃max − α̃min ) . (5.27)
q +σ
where the penalization power, p, is set to 3 and the tuning parameter, q, is set equal
to 0.1; see Sect. 4.1 for additional details and [24] for more information regarding
the specific problem setup and assumptions.
Thus, the full optimization problem may be formulated as
Find σ
Minimize Eq. (5.23)
Subject to (5.20)−(5.22)
Eqs.
σ dςd − vu ∗ 0
ςd
0∗σ ∗1
Given Eqs. (5.26) and (5.27).
Fig. 5.14 Assumed 2-D model with loads and boundary conditions (the model represents a thin
square plate with a center fluid inlet and fluid outlets along all four sides. The design domain is
subject to uniform heat generation). Reprinted from [24, Fig. 5a], Copyright (2012), with permission
from American Society of Mechanical Engineers (ASME)
Fig. 5.15 Optimal cooling channel topologies for the 2-D model from Fig. 5.14 obtained using
coarse (left) and fine (right) computational meshes. The color contours indicate the normalized
fluid velocity magnitude with larger velocity at higher elevations
both fine and coarse meshes having approximately 2,900 and 25,000 DOF, respec-
tively, to examine the effect of mesh refinement on the final topology since a filtering
routine was not implemented. The refined mesh result required approximately 1.4 h
computational time on a dual-core laptop with 2.5 GHz processors and 3 GB of RAM.
The optimal channel topologies for the model described above are shown in
Fig. 5.15, where the coarse mesh result is on the left and the fine mesh result is
on the right [24]. A 40 % fluid volume fraction was used in both cases. Observe that
a branching channel structure is obtained, and greater mesh refinement produces a
higher level of branching complexity.
Figure 5.16 shows the normalized temperature contour results for both mesh con-
ditions. Note that the ratio of the weightings for the thermal to fluid portions of
the objective function is approximately 30:1; refer to [24] for further details. These
82 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.16 Normalized temperature contours for the optimal cooling channel topologies from
Fig. 5.15 using coarse (left) and fine (right) computational meshes. Reprinted from [24, Figs. 9a,
b], Copyright (2012), with permission from American Society of Mechanical Engineers (ASME)
results illustrate the manner in which the branching channel structure decreases the
temperature of the design domain through effective fluid delivery to all portions of
the design domain. Consequently, the insulative porous solid material is forced out
toward the edges of the domain regardless of the magnitude of the inverse permeabil-
ity. The main diagonal branches of both structures then work to reduce the maximum
temperature, which tends to occur in the isolated corners of the domain.
Extensions to 3-D Design
The above topology optimization formulation for thermal-fluid systems may be
extended to handle 3-D laminar flow conjugate heat transfer problems with increased
computational cost. In the following brief example, a similar square thin plate (width
to thickness ratio = 10:1) optimization domain is considered with a center, fixed tem-
perature and velocity (laminar flow) fluid inlet and four, zero pressure side fluid
outlets. A heat flux was applied to the bottom wall of the domain, and a material
interpolation strategy similar to [95] was utilized. Here, solid regions are assigned a
greater thermal conductivity value than fluid regions (e.g., copper and water, respec-
tively) to better capture the 3-D fin effect of the solid portion of the heat exchange
structure.
A Pareto front for the design optimization problem is shown in Fig. 5.17, where
the normalized fluid temperature contours are superimposed on an isosurface of the
optimal cooling channel flow path with a slice through the solid domain (i.e., dark
regions). In this figure, design results with greater weighting of the thermal portion
of the objective function are found higher up along the ordinate axis. In contrast, the
design results with greater weighting of the fluid flow resistance term of the objective
function in Eq. (5.23) are shown further to the right along the abscissa. As a point
5.1 Electronic System Component Analysis and Design 83
0.8
Thermal Term Weight
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
Fluid Term Weight
Fig. 5.17 Pareto front for 3-D topology optimization of a thin plate heated from below and having
a center, fixed temperature and velocity fluid inlet with four, zero pressure side outlets. Normalized
temperature contours are superimposed on an isosurface of the optimal cooling channel flow path
with a slice through the solid domain (i.e., dark regions)
of reference, each solution shown in Fig. 5.17 was again obtained using commer-
cial software [15] with a computational mesh having √480,000 DOF and required
approximately 4 hr computational time on a six-core workstation with 3.46 GHz
processors and 48 GB of RAM.
Observe that a simple cross flow structure that minimizes the distance from the
fluid inlet to outlets is obtained when priority is placed on minimizing the fluid power
dissipated in the domain; refer to the lower right solution in Fig. 5.17. In contrast,
increasing the priority of the thermal term in the objective function (i.e., minimizing
the mean temperature of the domain, as well) produces branching cooling channel
topologies, as seen in the middle two images shown in Fig. 5.17. These cooling
84 5 Electromechanical System Simulation and Optimization Studies
channel solutions, which are very similar to those studied in the literature [89],
effectively manage the domain corner hot spots, as explained in detail in [24]. Finally,
by placing the most emphasis on the thermal term in the objective function, while
decreasing the weighting of the fluid flow resistance term, diagonally extending fluid
channels are obtained, as illustrated in the top left image in Fig. 5.17. Interestingly,
these results show similar characteristics to the 2-D results described in [24] and
presented above, which indicates that, for certain classes of problems, a 2-D modeling
approach may be sufficient in providing effective, manufacturable design solutions
in less time due to reduced computational cost. However, the fluid expansion (i.e.,
a 3-D effect) seen toward the edges of the design domain in the images in Fig. 5.17
indicates that for design domains of larger thickness, a fully 3-D method is likely the
best approach to capturing the relevant structural topology.
In this section, the derived 2-D cooling channel topology for a square design domain
presented above is used to demonstrate the utility of parametric size analysis tech-
niques for heat sink design; the study is based on the work in [21]. Here, the effect of
channel aspect ratio on the performance of the branching channel cold plate design
is examined, and the resulting trade-off between heat transfer and pumping power is
explained.
The coarse mesh topology optimization result (on the left in Fig. 5.15) was used
in the development of the CAD geometry for a cold plate structure, as illustrated
in Fig. 5.18, where a one-quarter symmetry model is assumed for computational
efficiency. The model geometry consists of a jet plate (shown transparent) with a
center nozzle that directs the coolant downwards into the channel structure, through
the radial branching channels, and then out through the sides of the cold plate.
5.1 Electronic System Component Analysis and Design 85
The coarse mesh topology optimization result from Sect. 5.1.3 was employed
due to the relatively small overall size of the cold plate and to reduce the overall
complexity of the microscale channel system. The size of the full plate is 17.2 mm
square with a 1.26 mm initial thickness. The initial channel height parameter, h c ,
value was set to 0.5 mm and was swept from this initial value up to 2.0 mm in
0.5 mm increments using the ’parametric sweep’ feature in a commercially available
software package. The cold plate channel width, wc , ranged from approximately 0.66
to 2.42 mm.
For comparison purposes, a straight microchannel cold plate was also studied.
A symmetry model of a single microchannel representing a 17.2 mm square cold
plate with an initial thickness 1.26 mm is shown in Fig. 5.19. The channel width
was set equal to the minimum channel width of the branching cold plate design, i.e.,
0.66 mm, while the channel height was again swept from 0.5 to 2 mm in 0.5 mm
increments.
In both cases, the channel aspect ratio was assumed as the channel height divided
by the minimum channel width, ω = h c /wc_min . Therefore, four separate values of
the channel aspect ratio (√0.75, 1.5, 2.25, and 3) were considered for each of the
described parametric computational models.
Parametric Model and Results
For each of the finned cold plate channel structures, shown in Figs. 5.18 and 5.19, a
thermally conductive heat sink material was assumed having a thermal conductivity,
specific heat capacity, and density of 160 W/(m K), 900 J/(kg K), and 2700 kg/m3 ,
respectively. The coolant fluid thermal conductivity, specific heat capacity, density,
and dynamic viscosity were, respectively, set to 0.654 W/(m K), 4182 J/(kg K),
982 kg/m3 , and 4.4E-4 Pa s. An insulative material was also assumed for the jet
plate and straight microchannel cap with a thermal conductivity of 0.26 W/(m K), a
specific heat capacity of 1700 J/(kg K), and a density of 1150 kg/m3 .
Appropriate symmetry boundary conditions were applied to each of the models
described in Figs. 5.18 and 5.19. Additionally, a uniform heat flux, q ∅∅ = 100 W/cm2 ,
was applied to the bottom side of each channel structure. The fluid inlet was set to an
86 5 Electromechanical System Simulation and Optimization Studies
383
375
373
365
363
355
353
345
353
338 338
Tmax
Fig. 5.20 Temperature contour results for branching microchannel cold plate with 0.15 L/min flow
rate and 0.5 mm channel height (shown transparent for clarity on left). Temperature contours for
straight microchannel cold plate with 0.15 L/min flow rate and 0.5 mm channel height (shown
transparent for clarity on right)
elevated temperature above 323 K. A zero pressure, convective flux condition was
set at each fluid outlet boundary.
An automated parameter study was set up using commercial finite element soft-
ware [16] to sweep the height parameter, h c , in each model from 0.5 to 2 mm, as
previously described. The computational models were solved at each channel height
parameter step over a range of fluid inlet volumetric flow rates spanning 0.025–
0.15 L/min. For each case, appropriate mesh refinement was selected with greater
element density in the fluid channels and at the channel walls.
The temperature contour results for the branching cold plate with a 0.5 mm channel
height at maximum inlet fluid flow rate are shown in Fig. 5.20 on the left. Likewise,
the results for the straight channel structure are shown in Fig. 5.20 on the right.
Observe that the maximum temperature, Tmax , for the branching channel design, 383
K, occurs at the outside corners of the square plate. For the straight channel design
the maximum temperature, Tmax = 391 K, occurs at the end of the channel after the
thermal boundary layer has developed.
Two performance metrics are used to evaluate the effectiveness of each of the
cold plate designs over the various height parameters and flow rate ranges. The unit
thermal resistance, Rth∅∅ , and pressure drop, P, are calculated for each cold plate as
∅∅ Tmax − Tin
Rth = , and (5.28)
q ∅∅
P = Pin − Pout . (5.29)
In Eq. (5.28), Tmax is the maximum cold plate temperature, Tin is the fluid inlet
temperature, and q ∅∅ is the applied heat flux. The pressure drop in Eq. (5.29) is
calculated as the difference between the cold plate inlet pressure, Pin , and outlet
pressure, Pout .
Figure 5.21 (left side image) provides the unit thermal resistance as a function
of flow rate for both cold plates over all channel height parameter values. The unit
thermal resistance of the branching microchannel cold plate for h c = 0.5 mm is less
5.1 Electronic System Component Analysis and Design 87
Fig. 5.21 Comparison of unit thermal resistance (left) and pressure drop (right) for branching (BC)
and straight microchannel (MC) cold plates
than that for the straight microchannel system. However, as the height and channel
aspect ratio is increased, the unit thermal resistance of the straight microchannel
system drops to values that are approximately one-half of those for the branching
microchannel system. Additionally, the branching channel cold plate unit thermal
resistance increases as channel height is increased from 0.5 to 1 mm. However, as
the channel height is further increased to 1.5 and 2 mm, the unit thermal resistance
then progressively decreases.
The pressure drop as a function of flow rate for both cold plate designs is also
shown in Fig. 5.21 (right side image). Observe that at the 0.5 mm channel height
the branching cold plate design exhibits a maximum pressure drop value of 0.24
kPa compared with 2.1 kPa for the straight microchannel design. This trend in mag-
nitude is consistent across the various channel aspect ratios where the branching
microchannel cold plate generally exhibits one to two orders of magnitude lower
pressure drop.
While the straight microchannel system outperforms the branching channel struc-
ture in terms of heat transfer at larger channel aspect ratios, the considerably larger
P may be prohibitive in applications where pumping power is limited. This char-
acteristic of microchannels is well known, and it is a primary obstacle to widespread
application [37].
The combined use of topology optimization with parametric design studies may
lead to interesting solutions to well-established engineering problems in terms of the
development of novel prototypes and products [25, 26, 28]. As a first example, the
heat transfer performance of this branching channel cold plate may be increased by
88 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.22 Layers 1 and 2 of a multi-pass branching microchannel cold plate. Layer 1 is shown
on the bottom row with the corresponding flow channel topology optimization result on the right.
Layer 2 is shown on the top row with the corresponding flow channel topology optimization result
on the right. Arrows indicate the fluid flow direction; color contours indicate the normalized fluid
velocity value with larger values at higher elevations
Outlet Inlet
No visible
bond lines
Fig. 5.23 Top view of fabricated multi-pass branching microchannel heat sink with cross-section
view shown below. The cross-section view is shown with the fluid flow path (blue arrows) and heat
input region (red arrows) overlaid. The multilayer cold plate design is fabricated out of aluminum
material using a diffusion bonding process to create a continuous thermal path through the thickness
of the assembly. The cross-section image highlights the fact that there are no visible bond lines post-
assembly
Thus, the integrated optimization, analysis, and design procedure may be exploited
in conjunction with the fluid manifold design procedure outlined Chap. 4 in the devel-
opment of a variety of efficient real-world heat transfer and fluid flow structures. It
should be noted that while conjugate heat transfer considering forced convection
and laminar flow was the main focus of this section, the formulation may be readily
extended to address additional effects such as natural convection, where the Boussi-
nesq approximation is applied to couple the fluid temperature and velocity fields [2].
Furthermore, the topic of thermal-fluid topology optimization for higher Re number
turbulent flows is also a developing research area [54] that will further increase the
potential of using such methods for realistic cold plate, heat exchanger, and cooling
jacket design.
90 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.24 Manifold hierarchical microchannel cold plate topology optimization results (left) and
prototype structure (right). Arrows on the left side image indicate the fluid flow direction; color
contours indicate the normalized fluid velocity value with larger values at higher elevations
Fig. 5.25 A schematic of a classic closed-loop magnetically controlled convective heat transfer sys-
tem illustrating the thermomagnetic siphoning effect. Reprinted with permission from [59, Fig. 1],
Copyright (2012), AIP Publishing LLC
While cold plate microchannel design for electronics applications was considered
in Sect. 5.1.3, a related concept is presented in this section, where the control of
fluid flow through the use of thermomagnetic field effects for heat transfer purposes
is considered. The use of stationary or time-varying magnetic fields to control the
motion of magnetically susceptible fluids is a well-known phenomenon [10, 36, 69,
80, 84, 85], and the use of this phenomenon represents an interesting concept for
the thermal management of electronic systems. Here, an optimal magnetic field dis-
tribution is shown to enhance thermomagnetic siphoning and consequently provide
increased heat transfer for localized heat sources. While implementation and assur-
ing correct fluid properties for maximum thermomagnetic siphoning effect [50] is
still challenging, exploiting such physical couplings for the passive pumping of fluid
may offer advantages in future heat transfer devices for electromechanical systems.
A schematic of a classic magnetically controlled convective heat transfer system
is shown in Fig. 5.25. In this system, a magnetic fluid container is subjected to a
5.1 Electronic System Component Analysis and Design 91
Step 1:
To – Initial temperature Final Solution
Br – Permanent magnet strength
Yes
No
Converged?
Update Br design
variable via MMA
Step 2: Step 6:
μr_mf – Magnetic permeability T – Temperature
Step 3: Step 5:
Step 4:
A – Magnetic vector potential u – Fluid velocity
fF – fluid body force
B – Magnetic flux density p – Pressure
Fig. 5.26 Flowchart of computations for the multiphysics analysis of a magnetically controlled
convective heat transfer system. Note the dashed portions indicate optimization steps. Reprinted
from [60, Fig. 2], Copyright (2012), with permission from American Society of Mechanical Engi-
neers (ASME)
temperature gradient, where the top side of the system is fixed at a cold temperature,
Tcold , and the bottom side is held at a hot temperature, Thot . Additionally, a magnetic
field gradient, −dH /dy, is assumed along the y-axis, as shown. As discussed in [84],
the magnetic susceptibility, ψ , of the fluid is inversely proportional to the temperature.
Thus, cold fluid at the top of the container is more strongly magnetized and drawn to
the region of higher magnetic field strength thus displacing the hot fluid at the bottom
of the container. This unique system involves three physics, including magnetic fields
generated either by permanent magnets (PMs) or electromagnets in air and external
to the magnetic fluid container, fluid flow emanating from the associated magnetic
fluid body force, and heat transfer related to the aforementioned fluid flow.
The flowchart for the multiphysics steady-state computational analysis of the
system just described is shown in Fig. 5.26. In this work, it is assumed that the
magnetic field gradient is generated by a PM positioned in an air domain outside
of the magnetic fluid container. Thus, the non-linear analysis procedure begins with
Step 1, where the temperature gradient (or heat source plus heat sink) is specified
along with the magnet residual magnetic flux density, Br , and an initial condition
system temperature, To . These inputs determine the magnetic fluid permeability,
μr _m f , via the temperature-dependent fluid susceptibility. After Step 2, a magnetic
field analysis is performed to compute the magnetic vector potential, A, and flux
density, B, (i.e., Step 3) which allows for the calculation of the magnetic body force,
92 5 Electromechanical System Simulation and Optimization Studies
f F , acting on the fluid (i.e., Step 4). A fluid flow analysis follows in Step 5, where
the fluid velocity vector, v, and pressure, P, state variables are determined. Finally,
the system temperature, T , distribution is determined in Step 6. Note that the dashed
portions of the flowchart in Fig. 5.26 refer to the optimization procedure that is
explained following the description of the governing equations below.
Since the magnetic field is generated by a PM, there is no external current in the
system, and Maxwell’s equations reduce to2
⎨ ⎩ ⎨ ⎩
1 1
◦× ◦ ×A =◦ × Br ; (B = ◦ × A) , (5.30)
μr μr
where μr is the relative magnetic permeability (including the surrounding air, perma-
nent magnet, and magnetic fluid regions of the analysis domain), A is the magnetic
vector potential, and the magnetic flux density and residual magnetic flux density of
the PM are B and Br , respectively.
Following [9], the permeability of the magnetic fluid, μr _m f , is dependent on the
fluid temperature, T , as follows:
⎨ ⎩
C ·γ
μr _m f = μo (1 + ψ ) , with ψ = 4∂ , (5.31)
T · Mw
where ψ is the fluid magnetic susceptibility, C is the Curie constant of the magnetic
fluid, γ is the fluid density, and Mw is the molecular weight.
Using Eq. (5.31), the magnetic fluid body force is then calculated as
1 1
f F = μo ψ ◦ H , with H =
2
◦ × A , (5.32)
2 μr _m f
where μo is the permeability of free space, and H is the magnitude of the applied
magnetic field. This body force term is then utilized in the fluid flow analysis for which
the Navier–Stokes equation and standard incompressibility constraint are assumed.
γ (v · ◦v) = −◦ P + ◦ · η ◦v + (◦v) T + f F , and (5.33)
◦ · v = 0. (5.34)
In the above expressions, the magnetic fluid dynamic viscosity is given as η, while
v is the fluid velocity vector, and P is the pressure state variable.
In terms of heat transfer within the system, the standard convection-diffusion
equation is assumed inside the fluid domain
γC p (v · ◦T ) = ◦ · (k◦T ) + Q, (5.35)
2 Reprinted, with permission, from [59], Copyright (2012), AIP Publishing LLC.
5.1 Electronic System Component Analysis and Design 93
where k and C p are the thermal conductivity and specific heat capacity of the magnetic
fluid, respectively, T is the temperature, and Q is the applied volumetric power
density. For any solid domain, Eq. (5.35) reduces to
0 = ◦ · (ks ◦T ) + Q s , (5.36)
where the subscript, s, indicates material physical parameters and thermal loads
associated with the solid regions of the system [83].
The multiphysics analysis of the system is performed by solving Eqs. (5.30)–
(5.36) using the finite element method. Using the finite element formulation,
Eqs. (5.30)–(5.36) can be represented in matrix form as
K M (T) · A = F M , (5.37)
K F (v) · (v P) = F F (A),
T
(5.38)
KT (v) · T = FT , (5.39)
where K M , K F , and KT are the global system stiffness matrices, and F M , F F , and
FT are the corresponding global force vectors for the magnetic, fluid, and thermal
equations, respectively. The coupled non-linear equations, Eqs. (5.37)–(5.39), are
then solved numerically using commercial finite element software [15].
As described in [59], two design variables, σ1 and σ2 , are assigned to each finite
element in the design domain for the optimization problem. The x and y direction
components of the residual flux density, Br , of each element are defined and inter-
polated as
Br,x = Br _P M · σ13 · sin{720⊥ · (σ2 − 0.5)}, (5.40)
In Eqs. (5.40)–(5.41), the first design variable, σ1 , controls the strength of the
permanent magnet, Br _P M . When σ1 = 0 the corresponding finite element represents
a void (i.e., air), and when σ1 is equal to unity, the element represents a PM. The
second design variable, σ2 , then controls the PM magnetization direction, where the
design variable σ2 is constrained between 0 and 1 (similar to σ1 ) as is usual in topology
optimization. When σ2 = 0.5 in Eqs. (5.40)–(5.41), the angle of the magnetization
direction relative to the y-axis is zero, and thus the magnetization direction points
along the positive y-direction.
The optimal PM locations and magnetization directions that minimize the maxi-
mum system temperature, Tmax , of the system,
Fo = Tmax , (5.42)
94 5 Electromechanical System Simulation and Optimization Studies
are then found as per the procedure described in [59]. Thus, the full thermomagnetic
convective optimization problem may be stated as
Find σ1 and σ2
Minimize Eq. (5.42)
Subject to Eqs. (5.37) − (5.39)
0 ∗ σ1 ∗ 1
0 ∗ σ2 ∗ 1
Given Eqs. (5.40) − (5.41).
Tcold
Thot
PM
Fig. 5.27 Assumed 2-D simulation domain for a classic closed-loop magnetically controlled
convective heat transfer system; refer to the corresponding conceptual schematic shown in Fig. 5.25
Fig. 5.28 Simulation results for the example system from Fig. 5.27: a magnetic vector potential
field lines; b normalized magnitude of the magnetic fluid body force distribution (note lighter regions
indicate larger magnitude); c magnetic fluid velocity contours with streamlines; and d temperature
contours
from the top side down towards the bottom side of the container. The corresponding
temperature contours are provided in Fig. 5.28d to illustrate the related convection
effect.
96 5 Electromechanical System Simulation and Optimization Studies
Heat source
Magnetic fluid container
PM design domain (2X)
Fig. 5.29 Assumed 2-D simulation domain for magnetically controlled convective heat transfer
system. Circular heat source is positioned in the center of a square magnetic fluid container. Two PM
design domains are shown (i.e., one each on the left and right sides) separated from the magnetic fluid
container by a thin aluminum layer. Convective cooling is assumed on the top and bottom boundaries
of the fluid container in addition to a surrounding air environment (not shown). Reprinted from
[60, Fig. 3], Copyright (2012), with permission from American Society of Mechanical Engineers
(ASME)
Note that the thermomagnetic siphoning described here is similar to the exten-
sively studied Rayleigh-Bénard convection transport phenomenon [84]. However,
in this case the buoyancy force due to thermal expansion of the fluid is replaced
by a similar magnetic fluid body force that is related to the temperature-dependent
magnetic fluid susceptibility. This convective instability is exploited in the topology
optimization example that follows.
fluid container. Without fluid motion, the maximum temperature of the system is
449.6 K, for reference.
To mesh the analysis domain, approximately 6,240 quadratic elements with 82,900
DOF were used including a surrounding air region (not shown in Fig. 5.29) for the
magnetic field analysis. It was found that higher element orders produced accurate
solutions for the more complex fluid flow patterns being considered. Thus, the fluid
velocity, v, and pressure, P, terms were interpolated using third and second order
finite elements, respectively. For the temperature, T , and magnetic vector potential,
A, second order finite elements were utilized. The solution of the full optimization
problem required 450 iterations or approximately 12.5 h on a quad-core workstation
with 3.0 GHz processors and 8 GB of RAM.
The optimal set of design variables, σ1 and σ2 , is obtained by solving the previously
described topology optimization problem, and the results are visualized in Fig. 5.30a.
The darker regions in the design domain indicate PM material, while the light colored
regions indicate air. Again, a volume constraint was not used to restrict the amount of
the PM material, and the solution is almost entirely PM with the exception of a small
air-gap magnetic flux barrier built in the middle of each design domain. The arrows
in Fig. 5.30a also show the magnetization direction of the PM material. As can be
seen in Fig. 5.30b, the optimal PM arrangement generates a strong magnetic field on
the fluid side of the system, while minimizing the magnetic field in the surrounding
air region on the opposite side; the resulting normalized magnetic fluid body force
distribution is shown in Fig. 5.30c, for reference. This special arrangement of the PM
magnetization direction, sometimes referred to as a Halbach array [39], is obtained
as the optimal PM layout for the computed magnetic field.
The fluid velocity contours with streamlines are shown in Fig. 5.30d, where the
average and maximum fluid velocities are, respectively, 3.37 × 10−3 m/s and 1.58 ×
10−2 m/s. These values are more than four times larger in magnitude when compared
with the best performer from non-optimized test cases. Additionally, a very clear four
eddy flow pattern is evident. Observe that the PM in each quadrant of the system acts
as a combined fluid ’pump’ plus ’conveyor’ to, respectively, draw the fluid in toward
the PM and then move the fluid along the length of the PM and vertical sides of the
container thus creating symmetric convective flow patterns that efficiently transport
the heat from source to sink. This result is further visualized in Fig. 5.30e, where the
maximum temperature of the system is 402 K (129 ⊥ C). Here, the heat is siphoned
off of the cylindrical heat source and distributed along the width of each heat sink at
the top and bottom of the magnetic fluid container.
For comparison purposes, several PM layouts were manually tested to check the
effectiveness of the thermomagnetic convection phenomenon without optimization.
The PM domains were split up into four blocks, and each quadrant domain was
filled with a 0.1 T PM with various assumed magnetization directions. The num-
ber of the possible combinations is large even if the magnetization directions are
constrained to be parallel to the Cartesian axes. So, for simplicity, only the four best
performing designs obtained through a manual study are presented here. The selected
magnetization directions are shown in Table 5.2, and the corresponding numerical
results are shown in Table 5.3. The latter table provides the computed values for
98 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.30 Design optimization results for system from Fig 5.29: a PM distribution and magneti-
zation direction (note dark regions = PM; light regions = air; arrows = magnetization direction);
b magnetic vector potential field lines; c normalized magnitude of the magnetic fluid body force
distribution (note lighter regions indicate larger magnitudes); d magnetic fluid velocity contours
with streamlines; and e temperature contours. Reprinted from [60, Fig. 4], Copyright (2012), with
permission from American Society of Mechanical Engineers (ASME)
5.1 Electronic System Component Analysis and Design 99
Table 5.2 Four PM magnetization direction layouts (indicated by arrows) assigned for comparison
with optimization results in Fig. 5.30
Note ‘h’ = horizontal; ‘v’ = vertical; ‘par’ = parallel; ‘sym’ = symmetric; ‘ant’ = anti-symmetric
Table 5.3 Fluid velocity contours plus magnetic vector potential lines for the cases from Table 5.2
Note ‘h’ = horizontal; ‘v’ = vertical; ‘par’ = parallel; ‘sym’ = symmetric; ‘ant’ = anti-symmetric;
‘rot’ = rotating; ‘K’ = temperature units; ‘m/s’ = velocity units
the maximum temperature, Tmax , average velocity magnitude of the fluid, vavg , and
maximum velocity magnitude of the fluid, vmax . The figures included in Table 5.3
show the corresponding fluid velocity contours with magnetic vector potential field
lines and arrows denoting the PM magnetization directions.
The results in Table 5.3 indicate that an increased number of recirculation zones
is generally not a good predictor of better heat transfer. For example, while the best
performer is Case I, which has six recirculation zones, with a 421.8 K maximum
temperature and 7.20 × 10−4 m/s average velocity, the second best performer (Case
II), has four recirculation zones, with a 24 % lower average velocity and a maximum
temperature that is only 0.9 K higher than Case I. Additionally, the average velocity
for Case II is smaller in magnitude than the average velocity for Case III which also
has six recirculation zones, yet Case II still results in a lower maximum temperature.
The maximum velocity is even less of a predictor of good heat transfer performance.
Furthermore, since the lowest maximum temperature achieved (i.e., Case I) is √20 K
higher than the optimized system shown in Fig. 5.30, intuition fails to provide an
optimum PM layout, and in all likelihood, there exists no analytical method to solve
this design problem.
100 5 Electromechanical System Simulation and Optimization Studies
Inclusion Matrix
x
z
Fig. 5.31 Anisotropic thermal-composite material with zoomed view of unit cell showing a ther-
mally conductive non-spherical inclusion embedded in a thermally conductive matrix. With kind
permission from [30, Fig. 1], © Springer-Verlag Berlin Heidelberg 2013
where the thermal conductivity, k, is equal in all directions and the volumetric heat
generated inside the domain, Q, flows in a direction perpendicular to the isothermal
surface passing through the point in space under consideration [79].
In contrast, two-phase composite materials comprising nonspherical inclusions
embedded in a surrounding matrix, such as the arrangement shown in Fig. 5.31, are
expected to exhibit anisotropic material thermal conductivity, where heat conduction
in the Cartesian coordinate system is governed by
⎨ ⎩
λ λT
Q=− ki j . (5.44)
λ xi λx j
For an anisotropic solid, the heat does not necessarily flow perpendicular to a
given isothermal surface, and the thermal conductivity involves nine conductivity
coefficients, ki j , which are components of a symmetric second-order tensor,4 k, [79]
as follows:
k11 k12 k13
ki j = ⎤ k21 k22 k23 , where (5.45)
k31 k32 k33
ki j = k ji i, j = 1, 2, 3. (5.46)
In the case of a thin square plate structure with a temperature gradient, −ϒx ,
along the x-axis shown in Fig. 5.32, the 3-D formulation above may be reduced to
describe heat conduction in two dimensions (2-D). Here, the x-y plane dependence
of the thermal conductivity on the inclusion angle, υ , may be approximated via the
4 Note that k and K throughout this section, i.e., Sect. 5.1.5 refer, respectively, to the second-
order thermal conductivity tensor and thermal conductivity matrix instead of local/global stiffness
matrices.
102 5 Electromechanical System Simulation and Optimization Studies
Inclusion
angle
φ
Heat flow
θ direction
-τx x
t
w
Fig. 5.32 A thin (i.e., t << w) square composite plate with anisotropic thermal conductivity
comprising elliptic cylinder inclusions embedded in a surrounding matrix and subjected to a tem-
perature gradient along the x-axis at an angle to the inclusion direction. With kind permission from
[30, Fig. 2], © Springer-Verlag Berlin Heidelberg 2013
where k1 and k2 are, respectively, the unit cell principal thermal conductivities parallel
and normal to the inclusion direction.
In contrast to [22] where k12 was assumed to be zero, shear terms are now repre-
sented by Eq. (5.49). Furthermore, elliptic cylinder type inclusions were arbitrarily
selected for the two-phase composite material microstructure; see Fig. 5.32. Thus,
instead of using a simplified composite slab model [22], the unit cell principal k1 and
k2 values may be determined for any shape and aspect ratio inclusion, per [42], via
f i (Ki − Km ) Km
Kuc = Km + . (5.50)
(1 − f i ) (Ki − Km ) S + Km
In this expression, Ki is the thermal conductivity matrix for the isotropic inclusion,
Km is the thermal conductivity matrix of the surrounding isotropic media, f i is the
inclusion volume fraction, and S is a second-order tensor analogous to the Eshelby
5.1 Electronic System Component Analysis and Design 103
tensor [34, 42, 43]. For an elliptic cylinder type inclusion, expressions for the com-
ponents of S in terms of the ellipse major, a1 , and minor, a2 , radii may be found
in [42], and the inclusion aspect ratio may be defined as
a1
αi = . (5.51)
a2
This detailed microstructure formulation allows for the investigation of the effect
of the inclusion aspect ratio on the guiding of heat flux within a composite. Also,
note that while a two-phase composite with elliptic cylinder type inclusions was
selected, [42] provides a microstructure formulation for three-phase systems as well
as the S tensor components for a variety of other standard inclusion geometries.
As indicated in Fig. 5.32, the heat flow within the composite does not necessarily
coincide with the direction of the temperature gradient, per [41], and the angle, ϕ,
between the inclusion direction and temperature gradient may be found using
This expression may be used to accurately compute the heat flow inside the composite
material.
Formulation of the Optimization Problem
The formulation of the material thermal conductivity design optimization problem
follows the detailed numerical procedure provided in [22]. A standard finite element
solver [15] is utilized for the steady-state conduction heat transfer problem, and a
gradient-based Method of Moving Asymptotes optimizer [90] is used to design the
material microstructure by determining the optimal inclusion angle, υ , that minimizes
a particular objective function.
To find the optimal inclusion angle distribution for the anisotropic material thermal
conductivity, the optimization problem is formulated as
Find σ
Minimize Fo
Subject to Eq. (5.44)
−1 ∗ σ ∗ 1
Here, Fo is a single or multiterm objective function, and the material design variable,
σ , is constrained on the interval [−1, 1]. Since the scalar function σ can take any
point-wise value, the raw function may produce oscillating designs and a Helmholtz
PDE-based filter is implemented in the optimization routine, per [51]. This filter acts
to produce a smooth scalar function, φ̃, that is then transformed into υ in Eqs. (5.47)–
(5.49) to cover the range of inclusion design angles covering −180⊥ to +180⊥ ; i.e.,
υ√ = φ̃∂ . Thus, the entire optimization problem is formulated as follows:
104 5 Electromechanical System Simulation and Optimization Studies
Adiabatic
Find σ
Minimize Fo
Subject to Eq. (5.44)
Eqs. (5.47) − (5.49)
υ = φ̃∂
−R 2f ◦ 2 φ̃ + φ̃ = σ
−1 ∗ σ ∗ 1
λ2T λ2T λT λT
Fo_g = k11 + k22 2 + 2k12 , (5.53)
λx 2 λy λx λy
where the terms k11 , k22 , and k12 in this expression are defined in Eqs. (5.47)–(5.49).
In contrast to [22] where an orthotropic material was assumed with k12 = 0, shear
terms are now represented by the third term on the right-hand side (RHS) in Eq. (5.53).
Heat Flux Shielding
Heat flux shielding may be realized by minimizing the magnitude of the temperature
gradients and heat flow within a selected region of a system. It can be achieved by
minimizing thermal compliance, Fo_g , in Eq. (5.53). Thus, the objective function for
shielding, Fo_s , may be adopted as
Note that in the case of the center domain with isotropic thermal conductivity, ko ,
shown in Fig. 5.33, k11 = k22 = ko and k12 = 0 resulting in the simplified objective
function ⎨ 2 ⎩
λ T λ2T
Fo_s_i = ko + . (5.55)
λx2 λ y2
and in case of the center domain with isotropic thermal conductivity in Fig. 5.33, a
simplified objective function is likewise obtained
q = −k◦T, (5.58)
where q is the heat flux vector at a spacial position. Expanding the RHS of Eq. (5.58)
gives
k11 λλTx + k12 λλTy
q=− . (5.59)
k21 λλTx + k22 λλTy
Assuming the heat flux vector, q, is aligned with a unit vector, ê, that lies along
a target angle direction in the global coordinate system; Eq. (5.59) suggests that the
direction of heat flow through a given spatial position may be controlled by defining
an objective function that realigns the heat flux vector with a target unit vector, −ê,
pointed in the opposite direction. Thus, the following multi-term objective function
is proposed:
Fo_r = w1 T1 + w2 T2 , (5.60)
where w1 and w2 are assumed weighting values and T1 and T2 are given by
T1 = ê · q, and (5.61)
2
T2 = ê × q . (5.62)
Physically, heat flux directivity is achieved by taking the scalar product of ê and
q via Eq. (5.61). By minimizing T1 , the component of q that is anti-parallel to ê is
maximized, However, this criteria does not strictly make q anti-parallel to ê since
it does not minimize the normal component of q directly. In order to fully align
the angle of q to ê, the normal component should be minimized. By including the
norm of the cross
product via Eq. (5.62), the angular error between ê and q is further
reduced (i.e., ê ≈q≈ sin (ϑ) n̂ ∀ 0 as ϑ ∀ 0 or ∂ ).
Numerical Studies
In this section, the validity of the objective functions from the previous section for
heat flux shielding, focusing, and reversal is examined via numerical experiments
involving the benchmark problem from Fig. 5.33. The goal is to control the flow of
heat within the isotropic center domain of Fig. 5.33 by designing the structural layout
of the anisotropic composite ring. Comparisons are made with the results reported
in [70].
Uniform hot and cold side temperature boundary conditions of 314 K and 273 K,
respectively, were applied to the structure shown in Fig. 5.33. As described in [70],
for the composite ring to blend into the background and not disturb the external
5.1 Electronic System Component Analysis and Design 107
180 º
0º
314 K
273 K
α= 1 α = 10 α = 100
Fig. 5.34 Composite ring heat flux shielding design results. The results in left, middle, and right
columns correspond to inclusion aspect ratios of αi = 1, 10, and 100, respectively. The images in
the top row show the absolute value of the composite structure design variable field. The images in
the bottom row show the temperature field with heat flux streamlines. With kind permission from
[30, Fig. 4], © Springer-Verlag Berlin Heidelberg 2013
temperature field, the thermal resistance of the surrounding medium (including cen-
ter region) should be close to the reduced average of the materials comprising the
two-phase ring, i.e., ki km √ ko2 . Thus, similar values of 2.6 W/(m K), 0.13 W/(m K),
and 0.56 W/(m K) were chosen for the thermal conductivity of the inclusion, matrix,
and surrounding media, respectively. An inclusion volume fraction of 0.5 with initial
vertical orientation was assumed, and various inclusion aspect ratios were exam-
ined in the numerical examples that follow. For all of the numerical results, a filter
radius, R f = 0.002, was used, and the optimization objective function was defined
exclusively on the center domain enclosed by the composite design domain.
Design Results: Heat Flux Shielding
The images in the left, center, and right columns in Fig. 5.34 show the heat flux
shielding design results for the composite ring using inclusion aspect ratios of αi = 1,
10, and 100, respectively. The top row in Fig. 5.34 shows the absolute value of the
structural design result, where dark colored regions indicate an inclusion angle of
zero degrees and light colored regions indicate an inclusion angle of 180⊥ . The
108 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.35 Composite ring structural design result for heat flux shielding with αi = 100. The
computed streamlines show the inclusion angle orientation and indicate that the inclusions are
aligned to form concentric rings. With kind permission from [30, Fig. 6], © Springer-Verlag Berlin
Heidelberg 2013
temperature field with heat flux streamlines, computed using Eq. (5.52), is shown in
the bottom row in Fig. 5.34.
Observe that the results in the leftmost column of Fig. 5.34 were obtained after a
single iteration and heat flux within the composite ring is not altered; a temperature
gradient (left to right from point A to B in Fig. 5.33) of 7.24 K exists in the center
domain. This result is due to an inclusion aspect ratio of αi = 1 that defines a circular
particle, and thus the material within the ring behaves in an isotropic fashion with
heat flux unmodified within the ring, as expected. As the particle aspect ratio is
increased to 10 and then 100, the solution shows increased heat flux shielding (as
evidenced by the streamlines in Fig. 5.34) with reduced temperature gradients in the
center domain of 3.33 K and 2.07 K, respectively.
In the limit, as αi ∀ ∈, the k1 and k2 principal thermal conductivity values
defined by the composite microstructure model, Eq. (5.50), approach values obtained
using a layered slab approach; refer to [22]. Thus, for higher inclusion aspect ratios
the solution for the composite ring quickly converges toward a structure comprising
approximately concentric laminations, as shown in Fig. 5.35. This structure matches
the composite ring design originally proposed in [70], where the temperature field
outside of the ring is unaffected and all of the heat flux manipulation occurs inside
the composite.
Design Results: Heat Flux Focusing
Several inclusion aspect ratios were examined for the case of heat flux focusing, where
the optimal structure was found using Eq. (5.56) after approximately 50 iterations.
The temperature gradient (left to right from point A to B in Fig. 5.33) in the center
domain for αi = 10 and 100 was 12.4 K, and 14.6 K, respectively, with over 100 %
increase in the heat flux for the higher aspect ratio. For reference, the design results
for αi = 100 are shown in Fig. 5.36, where again the temperature contours outside
of the ring are unaffected and all of the bending occurs within the composite ring.
Higher aspect ratio inclusions lead to a greater amount of heat flux focusing, and
5.1 Electronic System Component Analysis and Design 109
314 K 180 º
273 K 0º
Fig. 5.36 Composite ring heat flux focusing design results with inclusion aspect ratio of αi = 100.
The image on the left shows the temperature field with heat flux streamlines. The image on the right
shows the absolute value of the composite structure design variable field. With kind permission
from [30, Fig. 7], © Springer-Verlag Berlin Heidelberg 2013
Fig. 5.37 Composite ring structural design result for heat flux focusing with αi = 100. The
computed streamlines show the inclusion angle orientation and indicate that the inclusions are
aligned to form radial spokes. With kind permission from [30, Fig. 8], © Springer-Verlag Berlin
Heidelberg 2013
the resultant optimal structure in Fig. 5.37 resembles a composite ring composed of
radially extending laminated spokes, which is again similar to the assumed structure
for heat flux concentration in [70].
Design Results: Heat Flux Reversal
The structural design results for heat flux reversal are shown in Fig. 5.38. For the
case of heat flux reversal, it was found that an order-of-magnitude larger inclusion
thermal conductivity value, 50 W/(m K) instead of 2.6 W/(m K), was required to
trigger a twisting effect in the structural design result. Also, higher inclusion aspect
ratios were found to produce a greater heat flux reversal effect, and specifying the
weighting values w1 and w2 along with the target angle (i.e., ê-vector) in Eqs. (5.60)–
(5.62) allowed for specific modification of the twist angle.
Each set of results was obtained within 100 numerical iterations using αi = 104 .
The images in the columns of Fig. 5.38 (from left to right) are for different weighting
110 5 Electromechanical System Simulation and Optimization Studies
180 º
0º
314 K
273 K
w1 = 1 w1 = 1 w1 = 1 w1 = 1
w2 = 6.5 x 10-3 w2 = 6.5 x 10-3 w2 = 1 x 10-2 w2 = 1 x 10-2
Target Angle = 90° Target Angle = 120° Target Angle = 135° Target Angle = 180°
Fig. 5.38 Composite ring heat flux reversal design results with inclusion aspect ratio of αi = 104 .
The results in columns correspond to different w1 , w2 , and target angle (i.e., ê-vector) parameters
in Eqs. (5.60)–(5.62). The computed streamlines in the images in the top row show the inclusion
angle orientation and indicate that the inclusions are aligned to form approximate spiral structures.
The images in the middle row show the absolute value of the composite structure design variable
field. The images in the bottom row show the temperature field with heat flux streamlines. With
kind permission from [30, Fig. 9], © Springer-Verlag Berlin Heidelberg 2013
values and target angle combinations that produce a 90⊥ , 120⊥ , 135⊥ , and 180⊥ turn
in the heat flow. Complementary views of the composite structural design results are
shown in the top and middle rows of Fig. 5.38, while the temperature field with heat
flux streamlines is shown in the bottom row.
5.1 Electronic System Component Analysis and Design 111
Inspection of the lower left temperature field image in Fig. 5.38 reveals that
prioritizing heat flow in the vertical direction (w1 = 1; w2 = 6.5 × 10−3 ; ê = [0 1])
produces an exact 90⊥ turn in the heat flux. The corresponding structural design
results (left middle and top images) show that a spiral shape laminated structure is
formed to guide the heat flux across the center domain. However, the spiral is not
fully wrapped, and the upper left and lower right corners of the ring are designed such
that the inclusions follow a discontinuous path that guides the heat flux vertically
inward.
Moving to the second column from the left in Fig. 5.38, by specifying w1 = 1,
w2 = 6.5 × 10−3 , and ê = [0.5 0.866], the composite ring is designed to produce
a 120⊥ turn in the direction of the heat flux. Here, the spiral shape structure is less
discontinuous in the upper left and lower right regions of the ring, which creates a
smoother heat flow path around the annulus. This smoothing of the spiral continues
as seen in the third column from the left in Fig. 5.38, where the parameters w1 = 1,
w2 = 1 × 10−2 , and ê = [0.707 0.707]) generate a 135⊥ turn in the direction of the
heat flux.
Finally, a complete 180⊥ heat flux reversal across the center domain is obtained
by setting w1 = 1, w2 = 1 × 10−2 , and ê = [1 0]) in Eqs. (5.60)–(5.61); refer to
the lower right temperature field image in Fig. 5.38 where a temperature gradient
(left to right from point A to B in Fig. 5.33) of −4.3 K exists. Observe in the corre-
sponding structural design results (right middle and top images in Fig. 5.38) that the
arrangement of the inclusions is nearly continuous and approximates a logarithmic
spiral. Despite some bending of the temperature field outside of the composite ring
due to the larger assumed inclusion thermal conductivity value, this structural design
result matches the layout of the ’fan-like’ composite from [70] which was obtained
via an assumed coordinate transformation approach.
While the fundamental experimental results in [70] provide a satisfactory level of val-
idation for the presented optimization method for heat flux shielding (i.e., cloaking),
focusing, and reversal, the demonstration of such devices in ultra-thin composite
structural embodiments is a critical step to implementation in practical electronics
systems. The details of such representative structures may be found in [29], and are
summarized in this section. Specifically, two-material composite structures fabri-
cated using standard PCB manufacturing techniques are described. These example
composites were designed for heat flux cloaking, focusing, and reversal by etching
high thermal conductivity copper trace patterns in an engineered fashion on low
thermal conductivity PCB substrates. A schematic (including dimensions) of such a
PCB structure for heat flux cloaking is provided in Fig. 5.39a, where the left and right
ends of the structure are designed as thermal bus bars for the application and removal
of heat, respectively, to the structure. The center region of the structure consists of
concentric copper ring traces, as suggested by the design results shown in Fig. 5.35,
while the surrounding region comprises a copper mesh designed to cloak the cen-
112 5 Electromechanical System Simulation and Optimization Studies
Region A wm Region B
Region A
(a) z
y
Hot
p1 p2
x Cold
W
l1 T
wuc
Region C
L l1
Fig. 5.39 a A representative schematic of an ultrathin anisotropic composite device for heat flux
cloaking. Dimensions device length, L = 115 mm; device width, W = 50 mm; device thick-
ness, T =578 μm; mesh width, wm =200 μm; unit cell width, wuc = 2.5 mm; and thermal bus
bar length, l1 = 37.5 mm. Zoomed view image of the center region of each fabricated compos-
ite: b baseline device without center heat flow control region, c heat flux cloaking device (ring
OD √1.85 cm, ID √1 cm), d heat flux focusing device (ring OD √1.85 cm, ID √1 cm), and e heat
flux reversal device (ring OD 2.5 cm, ID √1 cm). Reprinted, with permission, from [29, Fig. 1],
Copyright (2013), AIP Publishing LLC
ter heat flow manipulation region; refer to [29] for details. In Fig 5.39b, a baseline
fabricated composite structure without a center heat flow control region is shown,
while Figs. 5.39c, d, and e, show PCB structures for heat flux cloaking, focusing,
and reversal, respectively. The focusing and cloaking composites, respectively, uti-
lize radial spoke traces and spiral traces in the center heat flux manipulation region
consistent with the design results shown in Figs. 5.37 and 5.38.
As explained in [29], a fixed temperature gradient of 35 K was applied across the
center region of each device using a copper block plus cartridge heater mounted to
the hot side thermal bus bar and a thermoelectric cooler connected to the cold side
thermal bus bar. The temperature field of the center region of each device was then
imaged using a calibrated infrared camera. Additional heat transfer simulations were
performed to verify the response of each PCB composite considering both conduction
and convection effects. The numerical and experimental results are shown in Fig. 5.40
with computational results shown on the top and experimental results shown on the
bottom; refer to [29] for further details.
Observe in Fig. 5.40 that the measured response of each thermal composite is
consistent with the heat flux cloaking, focusing, and reversal effects anticipated
in the prior optimization study. These results indicate that such anisotropic heat
5.1 Electronic System Component Analysis and Design 113
318 K
(a) (b) (c) (d)
p1 p2 p1 p2 p1 p2 p1 p2
283 K
318 K
(e) (f) (g) (h)
p1 p2 p1 p2 p1 p2 p1 p2
283 K
Fig. 5.40 Thermal contours of each device center region: a baseline—numerical, b heat
flux cloaking—numerical, c heat flux focusing—numerical, d heat flux reversal—numerical, e
baseline—experimental, f heat flux cloaking—experimental, g heat flux focusing—experimental,
and h heat flux reversal—experimental. Note points p1 and p2 indicate locations for temperature
gradient measurement; refer to Fig. 5.39. Reprinted, with permission from [29, Fig. 3], Copyright
(2013), AIP Publishing LLC
conduction models and associated optimization tools may be used with good accuracy
in the development of advanced PCB thermal composites for electronics thermal
management applications. Logical ongoing extensions of this work include the co-
design of multilayer PCBs for electrothermal functionality as well as laminated
composites for a variety of electromechanical applications in consumer electronics,
automotive, and aerospace fields.
where Nc is the number of coil winding turns. In terms of an electrical circuit, the
magnetic flux linkage is again defined as
where i is the coil electrical current. By equating Eqs. (5.63) and (5.64), the induc-
tance is then derived as ⎧
Nc
L= Bda. (5.65)
i
It should be noted that the magnetic flux density, B, inside the winding area is utilized
in the calculation of Eq. (5.65).
Method 2—Magnetic Energy In the second approach, the inductance is calculated
using a magnetic energy formulation. The magnetic energy, Wmag , in a linear material
may be assumed to be ⎧
1
Wmag = B 2 dv, (5.66)
2μ
5.2 Magnetic Component Analysis and Design 115
Fig. 5.41 The two-loss phenomena of ferromagnetic material in an inductor core: a eddy current
losses are due to induced eddy currents in the conductive material which generates Joule heating
losses; b hysteresis losses are represented by the total shaded area encompassed within the B-H
curve
where μ is the magnetic permeability and v is the volume of the analysis domain.
The magnetic energy of the electric circuit is defined as
1 2
Wmag = Li . (5.67)
2
By equating the two energy formulations, Eqs. (5.66) and (5.67), the inductance is
derived as ⎧
1
L= 2 B 2 dv (5.68)
μi
Here, the magnetic flux density, B, in the entire analysis domain is used for the
calculation.
Heat Loss Calculation
The heat losses in magnetic components can be classified into coil losses (i.e., load
losses) and core losses (i.e., no-load losses). Coil losses are related to Joule heating,
which results from the flow of the electric current through the coil, that is, the elec-
trically conductive material. The coil power loss, Pcoil , in Watts can be represented
as
Pcoil = Ri 2 , (5.69)
where R is the electric resistance of the coil. As the operating frequency of a magnetic
component increases, the electrical resistance increases due to skin and proximity
effects, which occur due to the magnetic field passing through the coil winding area.
As shown in Fig. 5.41, core loss is related to two separate effects that arise when
dealing with any ferromagnetic core material: (1) eddy current effects, and (2) hys-
teresis effects. Eddy current effects develop when electric currents are induced within
an electrically conductive material due to a change in the magnetic field. Associated
eddy current power loss, Peddy , is related to Joule heating as follows:
116 5 Electromechanical System Simulation and Optimization Studies
Air
Core
4mm
10mm
Coil Coil
2 3
mm mm
7.5mm 15mm
Fig. 5.42 Magnetic component analysis example: a 2-D inductor model comprising two E-shaped
core components and coils
⎧ ⎦
Peddy = ke f p2 B 2p dv, (5.70)
p
where kh is a coefficient determined using the hysteresis curve of the core material, n
is also a material coefficient usually between 1.5 and 2.5, and v is the domain volume
of the core material. The total core loss is then calculated as the sum of eddy current
losses, Eq. (5.70), and hysteresis losses, Eq. (5.71), which can be calculated using
the amplitude of p-th sinusoidal variation of magnetic flux density, B p , inside the
core material. Finally, the summation of the coil loss and total core losses produces
the total loss, Ploss_total , of a given magnetic component
(1GHz) Air-core
High
Small Small
Nano magnetic material
(Granular Film - CoZrO)
Future Development
Frequency
Operating
Ferrite
current
Size
(NiZn, MnZn)
- Amorphous material
- Permalloy
- Kool Mu
- Iron Powder Large
Low Large
Silicon Iron
(1kHz)
Low Loss High
Fig. 5.43 Frequency versus core loss for existing inductor core materials
the air-gap between two identically sized E-shaped core pieces is set to zero. The
thickness of the inductor core is defined as 10 mm. The coil is wound using 60 turns
of AWG 15 wire (diameter = 1.450 mm). The coil current is set as 10 A, and the
frequency of the current is set to 20 kHz.
Figure 5.43 shows an assessment of existing core materials for inductor design.
The inductor core loss is basically proportional to the device operating frequency
since eddy current and hysteresis losses are a function of the operating frequency;
refer to Eqs. (5.70) and (5.71). As the frequency increases, the inductor size can be
miniaturized because the reactance is proportional to the frequency. However, the
operating current and power is then limited due to increased core loss. Thus, an air-
core may be used for an inductor that operates at a very high frequency (e.g., GHz
regime).
Table 5.4 shows the assumed material physical parameters of the inductor com-
ponents, which are needed for the magnetic and thermal analysis in this example. A
magnetic analysis with harmonic excitation is performed by solving the governing
equations explained in Sect. 3.5. In the analysis, an air-region is included since the
permeability of the air and core material is of the same order, and an insulating condi-
tion is applied at the air boundary. From the analysis, the magnetic field distribution
is obtained. Then, the inductance and heat losses are calculated using the method
explained in Sect. 5.2.1. The thermal analysis is performed by solving the governing
118 5 Electromechanical System Simulation and Optimization Studies
(T)
(T)
Fig. 5.44 Two-dimensional magnetic field analysis result for the E-shaped core inductor of
Example 1. The color contours show the strength of the magnetic flux density, while the lines
indicate magnetic equipotential lines
heat diffusion equation introduced in Sect. 3.2. For simplicity in the thermal analy-
sis considered here, the air-region is excluded and a fixed temperature condition is
applied at the boundary of the core material (assumed here to be 293 K). However,
a convection heat transfer boundary condition could also be applied to the edge of
the core material, or the model could be further coupled to a third physical process
and include natural convection via a conjugate heat transfer analysis that includes
the surrounding air domain.
The numerical model was discretized using 6,762 second-order triangular ele-
ments, and the coupled magnetic and thermal analysis required approximately 5 s to
solve on a quad-core workstation with a 3.4 GHz processor and 16 GB of RAM. The
magnetic field distribution for the E-shaped core inductor is presented in Fig. 5.44.
Due to the eddy currents inside the core material, the magnetic field tends to be con-
centrated near the surface of the core material. The maximum magnetic flux density,
0.327 T, appears around the corner between the core and coil area. The inductance
is calculated using both methods explained above. The inductance is calculated as
0.331 μH using the first method based on the flux linkage formulation, Eq. (5.65),
and the inductance is 0.365 μH as determined using the second method based on the
magnetic energy formulation, Eq. (5.68). The two inductance calculation methods
give similar results. It should be noted again that the disadvantage of the first method
is that the inductance calculation result can vary depending on the domain defined for
5.2 Magnetic Component Analysis and Design 119
(K)
(K)
Fig. 5.45 Two-dimensional thermal analysis result for the E-shaped core inductor of Example 1.
The color contours show the temperature distribution
the flux linkage calculation. Thus, the second method might be preferred in that the
finite element analysis is based on an energy formulation, and the magnetic energy
for the entire domain is used for the calculation (i.e., there is no need to decide which
domain to use).
From the magnetic field distribution, the eddy current and hysteresis losses are,
respectively, calculated using Eqs. (5.70) and (5.71). In this example, the material-
dependent parameters ke and kh , and n in Eqs. (5.70) and (5.71) are, respectively, set
to 3.5, 1.5 and 1.8. The coil loss is simply calculated via Eq. (5.69). These three losses
are considered as the heat source for the ensuing thermal analysis. Then, a steady-state
thermal analysis is performed to calculate the inductor temperature distribution. The
analysis result is presented in Fig. 5.45, where the maximum temperature (302.1 K)
appears at the center of the inductor.
Extending the prior analysis, a two-dimensional inductor is designed using the struc-
tural topology optimization method. For the performance analysis of the inductor,
Maxwell’s equations are modified for a low-frequency time-harmonic analysis and
solved using the finite element method; see Sect. 3.5 for the formulated governing
120 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.46 Two-dimensional inductor design domain: a Example 1—coil only designed, b Example
2—coil and core designed simultaneously
equations. Here, two different examples are presented for the design optimization of
inductors.
In the first example, the inductor coil shape is designed using a topology opti-
mization approach with a fixed core shape, as shown in Fig. 5.46a. The second design
example then aims to find the optimal coil and core shape of an inductor simulta-
neously. The design domain for the second example is shown in Fig. 5.46b. The
design objective, i.e., to maximize the inductance, L, is the same in both numerical
examples, where the inductance is calculated from the magnetic field distribution
using Eq. (5.68). Thus, the objective function may be expressed as
⎧
1
Fo = B 2 dv, (5.73)
μi 2
and consequently, the optimization problems for the two examples are, respectively,
summarized as follows.
For Example 1, where the coil only is designed
Find σ
Minimize Eq. (5.73)
Subject to (3.34) − (3.38)
Eqs.
σ dςd − vu ∗ 0
ςd
0∗σ ∗1
Given ρ = ρcopper σ p ,
and for Example 2, where the coil and core are designed simultaneously
5.2 Magnetic Component Analysis and Design 121
Find σ1 and σ2
Minimize Eq. (5.73)
Subject to (3.34) − (3.38)
Eqs.
σ1 dςd − v1 ∗ 0, ςd σ2 dςd − v2 ∗ 0
ςd
0 ∗ σ1 ∗ 1, 0 ∗ σ2 ∗ 1
p p
Given ρ = ρcopper σ1 , μr = μr _steel σ2 .
In the latter optimization problem statement, σ1 and σ2 are design variables assigned
to each finite element, ρ is the electric conductivity, and μr is the relative magnetic
permeability. Here, the coil volume of the second optimization problem, that is v2 ,
is fixed to 10 % of the whole design domain.
The design results for the first example are presented in Fig. 5.47a– c, where
the full structure is shown by reflecting the results about the model symmetry line.
Here, the relative permeability and electric conductivity of the core material are
set as 150 and 1 S/m, respectively, and the electric conductivity of the copper is
set to 5.96 × 107 S/m. The design optimization is performed for three different
operating frequencies including 1 kHz, 5 kHz, and 10 kHz. The numerical model for
this example, as well as the following example, was discretized using 8,800 linear
quadrilateral elements, and each iteration for the optimization loop required 3 s to
solve on a quad-core workstation with a 3.4 GHz processor and 16 GB of RAM.
As can be seen in the design results in Fig. 5.47, the coil shape becomes thinner as
frequency increases. This result may be attributed to fact that the skin effect is more
prominent as frequency increases. Thus, several thin wires are more effective than
fewer thick wires in increasing the total current flowing through the coil, and this
consequently maximizes the inductance.
Finally, the design results for the second example are presented in Fig. 5.48. In
this case, the operating frequency is fixed at 500 kHz, and the relative permeability
and electric conductivity of the core material are again, respectively, set as 150 and
1 S/m. The optimization process is performed for different core volume constraint
values leading to the design results shown in Fig. 5.48. A circular shape is obtained
as the optimal coil design regardless of the core volume fraction. A rounded shape
is also obtained for the optimal core designs. As the volume of the core material
is increased from 40 to 90 %, the inductance, L, increases from 14.5 mH to 36.3
mH. However, the inductance does not increase much after reaching an 80 % core
volume fraction. This implies that the design engineer might be able to reduce the
core volume without losing the inductance after reaching a specific core volume
percentage, and the additional core might not be necessary for the magnetic field
distribution.
122 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.47 Design Example 1, inductor coil design results (top row) and current density distribution
(bottom row) for each operating frequency: a 1 kHz, b 5 kHz, and c 10 kHz. Note that for the design
results shown in the top row, the dark gray regions indicate the coil conductor material, while light
gray regions indicate the core material. For the current density distribution shown in the bottom
row, red-colored contours indicate a high current density and blue-colored contours indicate a low
current density
Fig. 5.48 Design Example 2, inductor core and coil design results: a 40 % core volume (L =
14.5 mH), b 60 % core volume (L = 27.4 mH), c 80 % core volume (L = 36.0 mH), and d 90 %
core volume (L = 36.3 mH). Note that the light gray regions indicate the coil conductor material,
while dark regions indicate the core material
5.3 RF Component Analysis and Design 123
In this section, we introduce two design studies for RF devices. The first study
pertains to a microstrip line device design, which involves conductive material design
with a frequency domain analysis. The second study then focuses on a dielectric
resonator antenna (DRA) with a time domain analysis plus a multiphysics extension.
where Z s is the surface impedance of the boundary, and H1 and H2 are, respectively,
the magnetic field on one side and the other side of the boundary. The above equation
balances the contribution of continuity across the boundary and conduction on the
surface via Z s . This means, this boundary condition is capable of representing both
conductor and insulator by changing the value of Z s as follows:
$
n̂ × (H1 − H2 ) = 0 (Z s ∀ ∈) Insulator
" # (5.77)
0 = n̂ × n̂ × E (Z s ∀ 0). Conductor
where the appended subscripts s_I and s_M , respectively, indicate values for an insu-
lator and conductor.
Formulation of Optimization Problem
Figure 5.49 shows one example of an analysis model used for designing a microstrip
patch antenna. As shown in this figure, the analysis domain contains a rectangular
box, PCB substrate, and planar square area, D, that is placed on the top side of the
box. The analysis domain is truncated by scattering boundaries, and the substrate
is suspended at the center. The bottom of the substrate is defined as a conductive
surface to represent the ground plane. On the top side of the substrate, a microstrip
line and square, D, is drawn. The input port is on the edge of the substrate.
126 5 Electromechanical System Simulation and Optimization Studies
where |S11 (ωi )| is the return loss at frequency point i, and wi is the weighting value
for each frequency point.
On the other hand, for designs dealing with power distribution, it is natural to use
the Poynting vector as the objective function
⎦
n
" #
Fo = wi n̂i · Re E × H∞ , (5.80)
i=1
Find σ
Minimize Eq. (5.79) or (5.80)
Subject to Maxwell’s equations; refer to Chap. 3, Sect. 3.6
0∗σ ∗1
Given Eq. (5.75) or (5.78).
Fig. 5.50 Projection path from design variable to surface impedance for the microstrip optimization
problem
Find σ
Minimize Eq. (5.79) or (5.80)
Subject to Maxwell’s equations; refer to Chap. 3, Sect. 3.6
σ = H̃(φ̃)
−R 2f ◦ 2 φ̃ + φ̃ = φ
−1 ∗ φ ∗ 1
Given Eq. (5.75) or (5.78),
where R f is the filter radius for Helmholtz filtering [51], as explained in Sect. 4.1
of Chap. 4. The projection and interpolation steps from design variable to surface
impedance are summarized in Fig. 5.50.
As previously described, this optimization problem can be solved by any type
of mathematical programming method. In this study, gradient-based non-linear pro-
gramming methods, i.e., Sequential Linear Programming, SLP [96], and Method of
Moving Asymptotes, MMA [90], are used with design sensitivity information that
is obtained efficiently by the adjoint method.
During the optimization procedure, the design variable may contain intermediate
values, or so-called grayscale. This is important to give slackness to the optimization
problem; however, such values should be eliminated in the final solution. In the
formulation above, it is possible to control the amount of the grayscale by tightening
the transition bandwidth, h t , in the relaxed Heaviside function Eq. (4.9). Here, h t is
128 5 Electromechanical System Simulation and Optimization Studies
where h t_min and h t_max are, respectively, the minimum and maximum values for h t ,
n is the current iteration step, w is the half bandwidth of transition in iteration steps,
and p is the control parameter of the curve profile.
Device Setup
A microstrip coupler for vertical interconnect on a multilayer PCB is designed using
a Heaviside projection formulation combined with MMA [90], as described above.
Figure 5.51 shows a schematic of the analysis model used for the coupler design
problem. The coupler is designed to connect microstrip lines on different layers of
a multilayer PCB. The PCB has three metal layers and two 0.5 mm thick dielectric
substrates between the metal layers. The input microstrip line, Line 1, is on the
bottom metal layer and connects Port 1 to the design domain. The design domain
consists of 10 mm squares on each metal layer. The middle metal layer is used as
the ground plane. On the top metal layer, the output microstrip lines, Lines 2 and 3,
are connected from the design domain and terminated by Ports 2 and 3, respectively.
All microstrip lines are 1.6 mm in width.
The Poynting vector-based formulation, Eq. (5.80), is used for the objective func-
tion. Each Poynting vector is calculated at the cross-section of each microstrip line
at the middle of the line. The weighting value is 0.1, −0.45, and 0.45 for Line 1, 2
and 3, respectively. The analysis frequency is 4.0 GHz.
5.3 RF Component Analysis and Design 129
Top layer
Middle layer
Bottom layer
Fig. 5.52 Density distribution of the initial and final microstrip coupler structure
Optimization Results
The optimization starts from a very trivial initial design, as shown in the left-hand
side column of Fig. 5.52. The value of the initial design variable, φinit , is all zero,
i.e., σ = 0.5. However, in order to avoid electric field concentration along the border
of the design domain, the value of σ on the boundaries is controlled to be a specified
value such that σ = 0 on insulative area boundaries and unity on conductive area
boundaries. This is done by enforcing a Dirichlet boundary condition of the projection
equation (i.e., scalar Helmholtz equation).
Figure 5.52 shows the initial and optimized density distribution on each layer. The
initial design just shows a vague distribution gradually varying from either metal or
insulator on the constrained boundary to intermediate material. On the other hand,
the optimized result shows a clear black and white topology with a very narrow
transition area from metal to insulator. Figure 5.53 shows the computed electric field
distribution at the middle of the lower and upper substrates of the device at 4.0 GHz.
Fabrication and Experiments
A multilayer printed circuit board fabrication technique is used for preparing a mea-
surement sample. Since the design obtained by topology optimization is a scalar field,
it has to be converted into a vectorial drawing that is compatible with CAD-based
fabrication systems. As observed in Fig. 5.52, the obtained design has a relatively
clear boundary between material and void regions, thanks to the transition bandwidth
130 5 Electromechanical System Simulation and Optimization Studies
Initial Optimal
Middle-Top
Bottom-Middle
0 2500[V/m] 0 3500[V/m]
y
Fig. 5.53 Electric field distribution at the middle of the lower and upper substrates at 4.0 GHz
control. Thus, the vectorial data can be extracted using ordinary contour extraction
algorithms. The extracted device CAD geometry is illustrated in Fig. 5.54 and a
representative prototype fabricated using this geometry is shown in Fig. 5.55.
To confirm the performance of the final device design, the response of the coupler
was predicted via a follow-up simulation using the synthesized CAD geometry. In
this model, the pattern is defined geometrically with spline curves. The microstrips
are extended with bends to have Port 1 and Port 2 at different positions on an edge
for attaching SMA connectors for measurements. The ground plane is also added on
the top and bottom side of the multilayer substrate.
The computed response of the microstrip coupler device is shown in Fig. 5.56.
The power level of the top layer is one-half of the power level on the bottom layer, and
power is equally distributed to the right and left on the top layer. The bottom row of
Fig. 5.56 shows a zoomed view of the coupler region; the field distribution is almost
the same as that shown in Fig. 5.53. Additionally, Fig. 5.57 shows the simulated
and experimentally measured S-parameters of the generated coupler device. The RF
microstrip coupler design has a reasonable working bandwidth with equal power
distribution, and a decent match is obtained between the simulation and experiment.
In summary, we proposed a flexible and efficient structural optimization method
for microstrip components in this section. The method is based on topology optimiza-
tion with adjoint analysis. The microstrip configuration is modeled by the conductiv-
ity of the surface. Numerical examples and representative prototypes were described
to demonstrate the practicality of the method.
5.3 RF Component Analysis and Design 131
z
y
x
Fig. 5.54 Extracted vector outlines of the various substrate layers (top row) and real microstrip
coupler device CAD model (middle row: top view, bottom row: perspective view)
132 5 Electromechanical System Simulation and Optimization Studies
Board layout
Top Bottom
Fig. 5.55 Images of a fabricated microstrip coupler device including board layout (top row) and
zoomed views of top and bottom side of the board in a local coupler region
DRAs are a class of antennas made with dielectric materials first proposed by Long
et al. in 1983 [63]. DRAs utilize the radiation phenomena of dielectric resonators
in open space, and they usually consist of a dielectric material, such as ceramic, in
a simple geometrical solid shape, which is fed by a probe or a microstrip. Since
mainly analytical approaches have been applied to evaluate DRA performance, the
proposed DRA shapes have been limited to simple geometries such as cylinders [63],
rectangular parallelepipeds [64], hemispheres [65], and various flat solids such as
rectangular discs [68], circular discs [61], and triangular shapes [62]. To enhance
the operational bandwidth of these antenna designs, several modifications of sim-
ple DRA shapes have been proposed, such as stacking resonators having different
permittivity [52], combining resonators having different resonance frequencies [35],
introducing an insert or an air-gap between the ground plane and the resonator [18],
making a cylindrical resonator annular [67], and providing a rectangular resonator
with notches [82].
Topology optimization is usually solved using the finite element method since it
is mathematically well established. Meanwhile, for engineering applications involv-
ing high frequency electromagnetics, different numerical analysis methods are fre-
quently used. Time domain methods are one group of such methods. Usually, the
finite element method is combined with a frequency domain formulation, as discussed
in Sect. 3.6, whereas the time domain method literary solves Maxwell’s equations
5.3 RF Component Analysis and Design 133
Top
Bottom
Bottom Top
x 0 30000[V/m]
Fig. 5.56 CAD model simulation of electric field distribution at the middle of the lower and upper
substrates at 4.0 GHz
0
-3
-10
S-parameters (dB)
-20
S11
S21
-30 S31
S11 (Sim.)
S21 (Sim.)
S31 (Sim.)
-40
3 3.5 4 4.5 5 5.5 6
Frequency (GHz)
where air and solid are the permittivity in air and the solid dielectric material, respec-
tively. While the SIMP method introduces a penalization parameter, p, to penalize
intermediate design values by modifying the relationship between weight and stiff-
ness for structural mechanics problems, here we use a simple linear interpolation by
setting p to 1.0, since permittivity is considered to change linearly with respect to
the normalized density, σ (x). In the presented numerical examples, we confirm that
improvements in the optimal configurations using different values of p are of mixed
utility and not entirely practical, when compared with the case where p is set to 1.0.
Formulation of Optimization Problem
Here, the reflection coefficient |S11 (ω)|, i.e., the return loss is formulated by the
following:
Pref (ω)
|S11 (ω)| = , (5.83)
Pin (ω)
which is a measure of the radiation efficiency for the electromagnetic waves radiating
from the design domain, where Pin (ω) and Pref (ω) are the input power and reflected
power in the frequency domain, respectively. That is, by minimizing |S11 (ω)|, the
radiated power, Prad (ω), which is a measure of the radiation performance of the
antenna, is maximized, since we are operating under the lossless assumption formu-
lated as
Pin (ω) = Prad (ω) + Pref (ω). (5.84)
Note that, Eq. (5.83) is most commonly used as a performance measure for antenna
development.
To calculate Pref (ω), an electric field in the form of a Gaussian pulse is fed in
at the feeding plane in transverse electric and magnetic (TEM) mode. This electric
field propagates in both upward and downward directions, but the power directed
136 5 Electromechanical System Simulation and Optimization Studies
downward toward the bottom end of the cable is absorbed. The rest of the power
directed upward propagates through the observation plane and the inactivated top
end absorbing boundary condition (ABC), and a portion of it radiates into the air and
is ultimately absorbed by the PML. Some portion of the field is reflected, reenters the
coaxial cable, propagates through the observation plane downward and is absorbed
by the ABC at the bottom end of the coaxial cable. The history of the electric field
at the observation plane is computed as Eobs (x, t) with respect to location, x, and
time, t, and we can obtain Pref (t) as a function of time by the following equation:
⎧
1
Pref (t) = (Eref (x, t))2 dςm , (5.85)
2
ςm
where Eref (x, t) denotes the reflected component of the electric field at point x, ςm
is the area of observation plane and Pref (ω) is obtained by Fourier transformation
as follows:
⎧∈
Pref (ω) = Pref (t)e−iωt dt, (5.86)
−∈
where Ein (x, t) is the electric field at the observation plane excited by the input pulse
alone, without reflection, and is calculated separately.
Since Pin (ω) is constant in this optimization problem, minimizing |S11 (ω)| is
equivalent to minimizing Pref (ω). When the input pulse has a power spectrum with
the same peak frequency as the antenna’s working frequency, and has the required
bandwidth, minimizing Pref (t) subject to the chosen input pulse minimizes Pref (ω)
in the specified bandwidth, with the input pulse weighted to have its peak at the
working frequency. Thus, the objective function is formulated as an integral value
of Pref (t) with respect to time, as follows:
⎧t f
Fo = Pref (t)dt, (5.88)
0
where t f is a fixed final time, and a Gaussian pulse having a center frequency equal to
the target frequency, and a −3 dB specified power bandwidth, is selected for the input
pulse. Therefore, the optimization problem is formulated as described in Sect. 4.1,
where the above objective function is minimized subject to a volume constraint, the
material interpolation function Eq. (5.82), and Maxwell’s equations, as discussed in
Chap. 3.
5.3 RF Component Analysis and Design 137
Adjoint Analysis
Following the related discussion in Sect. 4.1, adjoint analysis can also be carried out
in the time domain. First, the objective function expressed in Eq. (5.88) is rewritten
as follows:
⎧ ⎧t f
Fo = G(Ē, σ )dtdςm , (5.89)
ςm t=0
1
G(Ē, σ ) = |Eobs (t) − Ein (t)|2 (5.90)
2
On the other hand, by discretizing Maxwell’s equations in the time domain, the
following equation can be obtained:
¨ − [B(σ )]{Ē}
[M(σ )]{Ē} ˙ + [K (σ )]{Ē} = {Q}, (5.91)
where dot denotes the derivative with respect to time, t. [M(σ )], [B(σ )] and [K (σ )]
are system matrices and {Q} is the initial condition value for the time transient
equation. Here, the residual, R(E), is introduced as follows:
¨ + [B(σ )]{Ē}
R(E) = [M(σ )]{Ē} ˙ + [K (σ )]{Ē} − {Q}. (5.92)
where λ̂ is the time-dependent adjoint variable vector. Then, taking the first variation
of the above equation yields
⎧ ⎪
λFo∞ λFo∞
δFo∞ = {δ Ē} + {δσ } dςm . (5.94)
λ ĒT λσ T
ςm
where λ̂ can be obtained by solving the above equation. The sensitivity vector can
be calculated by substituting this value into Eq. (5.94). Equation (5.95) is equivalent
˙ term sign is negative, which implies that this is a
to Eq. (5.91), except that the {λ̂}
⎡ ⎣T
terminal value problem, and λG is a terminal value of this time transient equation.
λ Ē
Therefore, Eq. (5.95) can be solved using the same coefficient matrices by using
⎡ ⎣T
the backward time direction scheme, with λG as the input current pulse in the
λ Ē
observation domain of the objective function.
When the FDTD method is used in practical cases, these equations are approxi-
mated using finite difference grids in space and time. The design sensitivity of the
objective function with respect to design variable, σi , is calculated as
⎧ ⎧t f
λFo λR(Ē(t))
≈− Ê(t)T · dtdςm , (5.96)
λσi λσi
ςm_i 0
where Ē(t) is the electric field obtained by the forward FDTD method, Ê(t) is the
adjoint variable with respect to Ē(t), ςm_i is a domain of a Yee cell in which σi
is located, and the bar over E is added to indicate that the values of Ē have been
previously calculated by forward FDTD, and are not variables.
Single Band Design
Here, several design cases of FDTD topology optimization are presented.
Figures 5.59, 5.60, 5.61 and 5.62 show the optimized configurations (on left) and
return loss of the DRAs (on right). Each configuration shows a bottom view of the
antenna. For the return loss plots, the thin lines show the input pulse (green) and
reflected pulse (red), while the thick lines show the return loss of the initial structure
(green) and optimal structure (red). As seen in the return loss plots, the working
frequency band (i.e., the frequency band where the return loss is less than −10 dB)
broadly encompasses each target frequency.
The obtained structure shows clear separation of the void and material regions.
This is achieved by a technique termed ’subcutoff pulse’ [74]; refer to Fig. 5.63.
Without any treatment, wideband optimization is likely to yield solutions exhibit-
ing a large portion of grayscale area; however, the penalization scheme in SIMP
does not work, since the objective function does not monotonically depend on the
amount of material (in contrast with many structural design problems). However,
if the target frequency is unfeasibly low, the objective function shows a monotonic
relationship to the amount of material, because a larger resonator shows a lower
resonant frequency, such that, the optimizer tries to fill all of the design domain with
material. The subcutoff pulse technique utilizes this phenomenon. By mixing low
frequency components into the input pulse, a monotonic problem is obtained. If the
mixed frequency is below the cutoff frequency, it is infeasible. Additionally, since
it is necessary to keep the void areas to reduce the objective at the target frequency
5.3 RF Component Analysis and Design 139
0 1
Input
Return
-5 Initial 0.8
Optimal
-15 0.4
-20 0.2
-25 0
3 4 5 6 7 8 9
Fig. 5.59 DRA design targeted at 4.0 GHz. Bottom view of antenna shown (on left) with return
loss plot (on right), where the horizontal axis indicates frequency in units of GHz
0 1
Input
Return
-5 Initial 0.8
Optimal
Return Loss (dB)
-10 0.6
-15 0.4
-20 0.2
-25 0
3 4 5 6 7 8 9
Fig. 5.60 DRA design targeted at 5.0 GHz. Bottom view of antenna shown (on left) with return
loss plot (on right), where the horizontal axis indicates frequency in units of GHz
0 1
Input
Return
-5 Initial 0.8
Optimal
Return Loss (dB)
-10 0.6
-15 0.4
-20 0.2
-25 0
3 4 5 6 7 8 9
Fig. 5.61 DRA design targeted at 6.0 GHz. Bottom view of antenna shown (on left) with return
loss plot (on right), where the horizontal axis indicates frequency in units of GHz
140 5 Electromechanical System Simulation and Optimization Studies
0 1
Input
Return
-5 Initial 0.8
Optimal
-15 0.4
-20 0.2
-25 0
3 4 5 6 7 8 9
Fig. 5.62 DRA design targeted at 7.0 GHz. Bottom view of antenna shown (on left) with return
loss plot (on right), where the horizontal axis indicates frequency in units of GHz
band, an implicit upper volume constraint effect arises. These two effects help in
separating the void and material domains.
Multiphysics Design
As we see in Figs. 5.59, 5.60, 5.61 and 5.62, one side of each obtained structure is
floating from the ground. This is not a very good design in terms of installation of
the antenna since the gap distance can vary and sometimes antenna performance is
sensitive to such distance. In order to avoid such situations, multiphysics design is an
attractive solution. If the antenna structure is stiff enough in response to mechanical
loads, it is possible to avoid installation issues, while improving the reliability of the
ceramic part. For this purpose, the DRA design problem is solved in conjunction
with a structural mechanics problem. Additionally, by doing so, we can expect that
SIMP penalization becomes effective [73]. This may help to remove any remaining
grayscale elements observed in Fig. 5.60 or 5.62.
The schematic for the structural mechanics setup is shown in Fig. 5.64. Assuming
that the four corners of the bottom of the DRA are fixed, a vertical load is applied at
the top center of the structure. The DRA is designed to have minimum compliance
to show the maximum stiffness against the load.
5.3 RF Component Analysis and Design 141
Find σ (5.99a)
Pr e f l(u)
Minimize Fo = w + (1 − w) (5.99b)
σ Prinit
ef
l(u)init
Subject to 0 ∗ σ ∗ 1 (5.99c)
⎧
v= σ dς ∗ vu (5.99d)
D
Maxwell’s equations; refer to Chap. 3, Sect. 3.6 (5.99e)
Eq. (5.97) (5.99f)
Given Eqs. (5.100)–(5.101), (5.99g)
where w is the weighing coefficient, and the superscript, init , stands for an objective
value with initial design, i.e., σ = 0.5 everywhere. The governing equations are
coupled with the following material interpolation functions
142 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.65 Flowchart of Pareto front calculation for optimization problem, Eq. (5.99)
Here, the upper bound of material volume is constrained as Eq. (5.99d). This is to
utilize the penalization effect for grayscale values via a SIMP scheme. The SIMP
scheme shows an effect when an upper volume constraint is combined with a penalty
parameter value larger than unity. The penalty value is set to 3 as in Eq. (5.101). The
upper volume bound is set to 0.5.
Since this setup uses SIMP for grayscale penalization, the subcutoff pulse tech-
nique is not used. Instead, a multi-band design problem is specified using pulse
mixing. The input pulse is a summation of two sine modulated Gaussian pulses
whose center frequencies are 5.0 GHz and 7.0 GHz.
Figure 5.65 shows a flowchart for the Pareto front calculation. Starting from
w = 0, the multi-objective topology optimization problem is solved using both
FDTD (for the electromagnetics analysis) and FEM (for the structural analysis).
5.3 RF Component Analysis and Design 143
With both numerical methods, state variables are solved, and objective functions for
each physics and constraints are calculated. If the analyzed solution exhibits suffi-
cient performance, the design iteration loop is terminated. Otherwise, a sensitivity
analysis is performed. The electromagnetics problem involves adjoint analysis with
time backward calculation, while the structural problem does not require additional
analysis since it is a self adjoint problem. Using the obtained sensitivity information,
the design variables are then updated using non-linear programming methods; in
this study, SLP is used. When the solution is converged, the weighting coefficient
is increased and the optimizer is restarted from a new blank initial design. When
w > 1, all procedures are terminated.
Figure 5.66 shows the Pareto front obtained by the discussed optimization
problem setup. The vertical axis shows the normalized structural objective value,
l(u)/l(u)init , and the horizontal axis shows the normalized electromagnetics objec-
tive value, Pr e f /Prinit
e f . The markers show the values for w = 0.0, 0.01, 0.02,
0.1, 0.5, 0.7, 0.9, and 1.0. The corresponding antenna design configuration is also
shown in the figure. The red area indicates solid material, while transparent areas
indicate grayscale regions. Dotted lines radiating from coordinate (1, 1) show the
convergence path of each marker.
The Pareto front appears to be convex, considering the selected weighting values.
With a weighting value less than 0.5, the design demonstrates desirable stiffness,
and with w greater than or equal to 0.98, the solution shows good electromagnetics
performance. However, in terms of grayscale penalization, designs with w ≥ 0.90
are not sufficient. With precise analysis, we found that w = 0.03 is the solution that
best balances the multiple objectives and grayscale penalization. Figure 5.67 shows
the obtained antenna configuration with w = 0.03. The left side image shows the
top view and the right side image shows a section view of a cut plane positioned at a
distance of 3.5 mm from the base of the antenna. The obtained structure is a table-like
structure with four legs and no grayscale exists in the solution. Figure 5.68 shows
the frequency characteristics of the obtained design. The antenna has a sufficient,
−10 dB, working bandwidth at both 5.0 GHz and 7.0 GHz, the two target center
144 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.67 Antenna configuration obtained using w = 0.03. Top view (left image) and cut plane
view at 3.5 mm level (right image)
frequencies. For reference, Fig. 5.69 also shows the antenna configurations with top,
sectioned, and bottom views obtained for weighting values of 0, 0.5, 0.9, and 1 from
left to right, respectively.
Fabrication Constraint
In the previous section, a multiband DRA without grayscale was successfully
designed with a weighted sum multi-objective formulation and SIMP. However,
the process for finding the appropriate weighting coefficient is inefficient. Addition-
ally, as shown in Figs. 5.67 and 5.69, the obtained designs are hard to realize since
they contain irregular cavities and holes. In this section, a problem formulation that
provides for a fabrication constraint is explained.
The constraint method is another well-known approach for multi-objective prob-
lems. The constraint method transforms the multi-objective optimization problem
into a single objective problem by converting objectives into constraints with the
exception of the one single main objective. In our case, the electromagnetics objec-
tive is converted into a constraint since the maximum allowable reflection energy is
easier to estimate than the maximum allowable elastic compliance. In this way, it is
not necessary to calculate the Pareto front to find an optimum design that satisfies the
5.3 RF Component Analysis and Design 145
Top view
Cut view
Bottom view 0.00 0.50 0.90 1.00
Fig. 5.69 Antenna configurations for w = 0, 0.5, 0.9, and 1 shown in columns from left to right,
respectively. Top view (top row), cut view (middle row), and bottom view (lower row)
designer’s requirement. Thus, the first term in Eq. (5.99b) is converted as follows:
Pr e f Pu
F¯o = init < init , (5.102)
Pr e f P
where H is the Heaviside step function. When the exclusive area is completely filled
with solid material, there exists no grayscale. This condition can be written as
vex = v, (5.104)
where v is the volume ( D σ dς) commonly used in topology optimization, e.g.
Eq. (5.99d). In the implementation, this expression is rewritten in a relaxed form
⎧ ⎧
ṽex − v = H̃(σ (x))dς − σ dς < Ωin f , (5.105)
D D
146 5 Electromechanical System Simulation and Optimization Studies
Opposite
taper
feature
Mold
(bottom)
where H̃ is the relaxed Heaviside function, Eq. (4.9), introduced in Chap. 4, and
Ωin f is an infinitesimal value. By enforcing the above constraint, grayscale can be
eliminated from the design.
Assume that the antenna is made of high permittivity thermoplastic material, such
as ceramic filled polymer composite, and is fabricated using molding technologies;
refer to Fig. 5.70. The outer shape should reside in a specified area, due to cosmetic
or packaging concerns. Also, the inner or bottom shape should have no negative taper
to assure the product can be removed from the mold; see Fig. 5.71.
For the outer surface, one straightforward approach is to prepare a fixed design
domain which satisfies the specification. However, this is not very flexible for adopt-
ing specification changes. Instead, the following constraints may be used:
where, Q̄(x) is a level set function to represent a given outer shape ( Q̄(x) > 0 is the
domain side), thus Q̄(x) itself is constant throughout the optimization. Note that the
design domain is discretized using a structured cubic grid since FDTD is used for
the electromagnetics analysis, thus finding the adjacent design variable is trivial in
contrast with using a finite element mesh.
5.3 RF Component Analysis and Design 147
Negative taper occurs if a point has more material density than a point above the
point of interest. Therefore, in order to eliminate the negative taper, the following
constraint may be used:
where z is the grid size in z-direction. Note that the condition, Q̄(x) ≥ z is a
necessary and sufficient condition to avoid numerical contradiction, so that it can
be replaced by other necessary conditions, as we show in the following complete
formulation.
Here, we specify Q̄(x) as a spherical shell with center at xc and radius r0 , that is
Q̄(x) = r0 − |x − xc |. By combining all of the constraints above, the optimization
problem may be formulated as follows:
Find σ (5.108a)
l(u)
Minimize Fo = (5.108b)
0<σ <1 l(u)init
Pr e f Pu
Subject to F¯o = init < init , (5.108c)
Pr e f P
vex = v Exclusive volume (5.108d)
R1 = σ (x) < 0 if r0 < |x − xc |
Outer shape control (void domain)
(5.108e)
T1 = σ (x − z) − σ (x) < 0 if |x − xc | < r0 − r
Taper control (5.108f)
Maxwell’s equations; refer to Chap. 3, Sect. 3.6 (5.108g)
Eq. (5.97) (5.108h)
Given Eqs. (5.100)–(5.101). (5.108i)
Here, constraint Eqs. (5.108c) and (5.108d) are non-linear, whereas Eqs. (5.108e)
and (5.108f) are linear constraints, where r ≥ z. SLP is used as the non-linear
programming method, the same as the previous problems. For the non-linear con-
straints, an augmented Lagrangian method is used. The linear constraints in the
formulation are pointwise constraints so that the total number of them is approxi-
mately the same as the number of design variables. These constraints are enforced
in subproblems in SLP.
Figure 5.72 shows the flowchart for the full optimization problem, Eq. (5.108).
Compared with the previous flowchart, this process is more simple with only a single
148 5 Electromechanical System Simulation and Optimization Studies
loop, thus it is more computationally efficient. However, the number of the iterations
using a single loop process may take more time because of the non-linear constraints.
Still, this modified process is beneficial since the user does not have to decide an
appropriate weighting coefficient set to cover the Pareto front.
Figure 5.73 shows the optimized configuration with and without the fabrica-
tion constraints. The case without the fabrication constraints is solved without
Eqs. (5.108d), (5.108e) and (5.108f), but instead, with a normal upper volume con-
straint, i.e., Eq. (5.99d). The upper row shows a view from the top side, the middle
row shows the same top side view with a cross section obtained by a vertical cut
plane, and the lower row shows a view from the bottom side. The left column shows
the result without fabrication constraints, while the right column shows the result
with fabrication constraints. Both methods produce a configuration that has a small
amount of grayscale area. This indicates the constraint-based approach for grayscale
control is a better way to balance two objectives, in this case. The cross-sectional view
reveals that the left side configuration has small cavities or inverse taper, whereas
5.3 RF Component Analysis and Design 149
Fig. 5.73 Antenna config- Without fabrication constraint With fabrication constraints
urations obtained with and
without fabrication constraints
Top view
Cut view
Bottom view
the right hand configuration does not, and maintains a mechanically stable structure
with supports at the four bottom corners.
In this section, the analysis and design optimization of actuators is explained. Two
separate 2-D numerical examples are provided covering common topics including the
magnetostructural design of a solenoid and the design of an actuator with respect to
the optimal arrangement of the permanent magnet, coils, and ferromagnetic material.
The multiphysics nature of the first study elucidates the relationship between the
magnetic force generated by a solenoid and the resulting structural deformation
that occurs within the system. The second study then provides a detailed example
of how to optimally arrange multiple components in an actuator to maximize the
average magnetic force acting on the moving component (i.e., plunger). Here, the
magnetic field is logically connected to the structural response. Both studies have
clear ramifications in the development electromechanical assemblies, where over-
built and bulky structures that utilize excessive amounts of expensive rare-earth
materials are detrimental to the design of highly efficient systems.
150 5 Electromechanical System Simulation and Optimization Studies
In the first example, the solenoid actuator shown in Fig. 5.74 is optimally designed
for the improvement of coupled magnetostructural performance. More specifically,
a structural topology optimization approach is applied for the design of an actuator
armature, and the design goal is chosen to maximize the total magnetic force acting
on the armature and to minimize the mechanical deformation of the armature caused
by the distributed magnetic force [56]. A magnetostatic analysis is performed first,
and the distribution of the magnetic force is obtained. Next, the structural analysis is
performed when the calculated force distribution is applied to the actuator. Finally,
the mechanical compliance is computed. For the design optimization, the sensitivity
analysis is performed using the adjoint method, and the optimization problem is
formulated and solved using the sequential linear programming (SLP) method.
There are two distributed force calculation methods: (a) surface force distribution
and (b) body force distribution. Using the first method, the force distribution on the
surface of the magnetic materials is found, while using the second method, the dis-
tributed body force inside the magnetic materials is calculated using a virtual air-gap
scheme [12, 13]. Both methods utilize general force calculation formulations such
as Maxwell Stress Tensor (MST) or Coulomb Virtual Work (CVW) schemes. There
has been inconclusive debate about whether the surface or body force distribution is
a truer representation of the real situation. The surface and body force distributions
obtained from the same magnetic field distribution are completely different. Thus,
the mechanical deformation generated by two force distribution results may be dif-
ferent, as well. Here, the design optimization results for both surface and body force
distributions are compared. Differences in optimization results are identified, and
it is validated that each result is optimal in its force distribution. Such an example
makes clear the multiphysics nature of actuator design.
5.4 Actuator Analysis and Design 151
Magnetostructural Analysis
The coupled magnetostructural analysis is performed using the finite element method.
The analysis steps can be summarized as follows: (1) The magnetostatic analysis is
carried out using the magnetic vector potential, A; (2) The magnetic force distri-
bution, f M_d , is calculated from the magnetic field distribution; (3) The magnetic
force distribution yields the structural nodal force vector, f S . Then, the full structural
analysis for the magnetic body force gives the mechanical displacement vector, U,
and the compliance, l. The detailed procedure for each step is explained below.
For the first step, the magnetostatic governing equation (see Sect. 3.5) is formu-
lated in matrix form using the finite element method as
KM A = FM , (5.109)
where K M is the magnetic stiffness matrix, A is the nodal vector for magnetic vector
potential, and F M is the magnetic load vector. After solving Eq. (5.109), the magnetic
flux density B is calculated as
B = ◦ × A. (5.110)
Second, from the magnetic field, B, the surface or body force distribution, fb , is
calculated using either the MST or CVW method. In this example, the MST method
is applied for computing both surface and body force distributions. For the calculation
of the surface force distribution, the MST formulation is applied along the surface
of the magnetic material. For the body force distribution, a virtual air-gap is inserted
between the small bodies contacting each other because the MST method is defined
only on an object surrounded by air. The magnetic flux density, B, is decomposed
into normal and tangential components along the surface lines as follows:
B = Bn n̂ + Bt t̂, (5.111)
where n̂ is a normal unit vector and t̂ is a tangential unit vector with respect to the
surface lines. Then, the MST method gives the distributed magnetic force
⎪( ⎪(
1 1
f M_d = (Bn Bt )ds n̂ + (B − Bt )ds t̂.
2 2
(5.112)
μo 2μo n
Third, dividing this expression by the small body volume provides the distributed
body force, fb , which is used to obtain the structural nodal force vector, f S . In the
j
finite element formulation, the nodal force at the j-th node, f S , is represented as
⎧
fS = N j fb dς, (5.113)
152 5 Electromechanical System Simulation and Optimization Studies
where N j is the shape function of j-th node. Then, the finite element equation for
the structural analysis can be written as
KS U = FS , (5.114)
where K S is the assembled structural stiffness matrix, and U is the mechanical dis-
placement vector. The assembled structural force vector, F S , is a function of the
magnetic vector potential, A. Therefore, the structural analysis is coupled with the
magnetostatic analysis.
By solving Eq. (5.114), U is obtained, and finally the mechanical compliance, l,
is calculated as
l = FTS U. (5.115)
where the objective function, Fo , is defined as the scaled sum of the mechanical
compliance and the total magnetic force.
The scaling factors, s1 and s2 , are controlled in each iteration step using the
adaptive scaling strategy proposed in [86] and written as
1/(Ftot + Ωin f )
s1 = s1o , and (5.117)
1/(FTS U + Ωin f ) + 1/(Ftot + Ωin f )
1/(FTS U + Ωin f )
s2 = s2o . (5.118)
1/(FTS U + Ωin f ) + 1/(Ftot + Ωin f )
the other term of the objective function. When the compliance, FTS U, is dominantly
reduced, the scaling factor, s2 , increases to reduce the total force, Ftot , and vice-
versa. The term Ωin f is an additional infinitesimal number, which is utilized to avoid
a ’divided-by-zero’ situation. Then, the optimum solution that satisfies both objective
function terms can be obtained.
In sum, the optimization problem for the coupled magnetostructural problem can
be formulated as
Find σ
Minimize Eq. (5.116)
Subject to Eqs. (5.109) − (5.115)
ςd σ dςd − v1 ∗ 0
0∗σ ∗1
Given μr = μr _steel σ p
E = E steel σ p .
Here, the volume of the design domain is constrained to a specified volume fraction
v1 , μr is the relative magnetic permeability, and E is the elastic modulus. The above
optimization problem is solved by SLP, which requires the design sensitivity.
For this example, the sensitivity of the objective function is analytically derived
using the adjoint variable method. The derivative of the function Fo with respect to
the density design variable, σ , is written as
⎪
dFo d(F S (A(σ ), σ )T U(σ )) d(Ftot (A(σ ), σ ))
= s1 − s2 . (5.119)
dσ dσ dσ
By using the adjoint variable method with two weakly coupled finite element equa-
tions, Eqs. (5.109) and (5.114), the first differentiation term in Eq. (5.119) is derived
as
dFTS U λF S (A(σ ), σ )T λKM K
S
=2 U + λ̂T1 A + λ̂T2 − U . (5.120)
dσ λσ λσ λσ
The adjoint variables, λ̂1 and λ̂2 , are calculated by solving the following adjoint
equations:
λFT
K M λ̂1 = −2 S U, and (5.121)
λA
K S λ̂2 = F S . (5.122)
The second differentiation term of Eq. (5.119) is also derived using the adjoint
variable method as
154 5 Electromechanical System Simulation and Optimization Studies
d Ftot λ Ftot λK M
= + λ̂T3 A. (5.123)
dσ λσ λσ
λ Ftot
T
K M λ̂3 = − . (5.124)
λA
Finally, the analytical sensitivity of the objective function is obtained from
Eqs. (5.119)–(5.124).
Optimization Model and Results
The solenoid actuator shown in Fig. 5.74 is designed using the above analysis and
optimization approach. First, the actuator is simplified into a 2-D model, as shown
in Fig. 5.75. The design domain is defined at the armature, and the displacement of
the armature right boundary is fixed to zero. The design volume fraction is set as 0.6,
and the SLP move limit is set as 0.004.
The two armature design results with equipotential lines are presented in Fig. 5.76.
The numerical model was discretized using approximately 9,200 linear quadrilateral
elements. One iteration for the optimization routine required about 5 s to solve on a
quad-core workstation with a 3.4 GHz processor and 8 GB of RAM. The first design,
AD1 (shown on the left), is obtained when the distributed magnetic force is calculated
as a body force, and the second design, AD2 (shown on the right), is obtained when
it is calculated as a surface force. Both design results, AD1 and AD2, have a similar
general structure configuration. The structure connecting the left upper and bottom
boundaries is designed to maximize the total force, Ftot . Then, most of the magnetic
field flows through the armature, as is visualized through the equipotential lines.
On the other hand, the structure in a right upper region is designed to minimize the
mechanical deformation.
5.4 Actuator Analysis and Design 155
Fig. 5.76 Armature Design AD1 (left) obtained when the distributed magnetic force is calculated
as a body force, and Armature Design AD2 (right) obtained when the distributed magnetic force
is calculated as a surface force. Note that the magnetic equipotential lines are overlaid, for clarity.
© 2010 IEEE. Reprinted with permission from [56, Figs. 5 and 6]
Fig. 5.77 Body force distribution of design AD1 (left) and design AD2 (right) in grayscale. © 2010
IEEE. Reprinted with permission from [56, Fig. 7]
Despite the similar general configuration, the two design results also show the
following (considerable) differences. The design AD1 contains a large hole in the
bottom part to raise the location connecting the armature and fixed boundary, while
the design AD2 has a sharp corner in this same region. In addition, the design AD1
consists of a thin structure and small hole in the right upper region, while the design
AD2 has a thick structure in this location. These differences are due to the minimiza-
tion of the compliance via different force distributions.
To validate that each design result is optimal in its own force distribution, the
body and surface force distribution of both design results are investigated. The body
force distributions of the two design results AD1 and AD2 are presented as grayscale
images in Fig. 5.77. The black colored regions correspond to a strong body force, and
156 5 Electromechanical System Simulation and Optimization Studies
Table 5.5 Comparison of the total magnetic force and mechanical compliance for AD1 and AD2
actuator designs
Design AD1 Design AD2
Body force Ftot,x = −2680.22 N Ftot,y = −2680.22 N Ftot,x = −2674.19 N Ftot,y = −2047.13 N
Compliance = 0.6237 × 10−3 m/N Compliance = 0.1355 × 10−2 m/N
Surface force Ftot,x = −2690.49 N Ftot,y = −2054.20 N Ftot,x = −2685.01 N Ftot,y = −2051.56 N
Compliance = 0.4747 × 10−3 m/N Compliance = 0.2515 × 10−3 m/N
© 2010 IEEE. Reprinted with permission from [56], Table I
the white color regions indicate a weak body force. The surface force distributions of
two design results are presented in Fig. 5.78. Table 5.5 compares the x and y-direction
total force and mechanical compliance of each design.
If the body force shown in Fig. 5.77 is the true force distribution, the design AD1
is much stiffer than the design AD2 as can be seen in the compliance comparison in
Table 5.5. The body force of the design AD1 is evenly distributed around the center
of the structure, and the compliance due to this distribution can be minimized by the
raised connecting location. In contrast, a sharp corner of the design AD2 concentrates
the body force, and consequently deteriorates the stiffness. On the other hand, if the
surface force distribution shown in Fig. 5.78 is true, the design AD2 is much stiffer
than the design AD1. The surface force of two designs is similarly distributed at
the left upper and bottom boundaries. To minimize the compliance of this force
distribution, a thick structure in the right upper part is required as can be seen in the
design AD2. The total forces of all cases are almost identical as shown in Table 5.5.
Therefore, both design AD1 and AD2 can satisfy the first design goal, that is, the
maximization of the total force, whichever of the force distribution is true. However,
each design can satisfy the compliance minimization only when the corresponding
force distribution is true.
5.4 Actuator Analysis and Design 157
KM A = FM , (5.125)
where K M is the magnetic stiffness matrix, A is the nodal vector for magnetic vector
potential, and F M is the magnetic load vector. After solving Eq. (5.125), the magnetic
flux density B is calculated as
B = ◦ × A. (5.126)
The previous example in Sect. 5.4.1 deals with using the distributed force for the
structural analysis. However, this example addresses only the total force acting on
the plunger. Note that the integration of the distributed force along the body is the
total magnetic force acting on the body. The MST method gives the total magnetic
force, Ftot , as
⎪( ⎪(
1 1
Ftot = (B − Bt )ds n̂ +
2 2
(Bn Bt )ds t̂, (5.127)
2μo n μo
158 5 Electromechanical System Simulation and Optimization Studies
where Bn and Bt are, respectively, the normal and tangential components of magnetic
flux density, B, while n̂ and t̂ are, respectively, unit vectors normal and tangential
to the integration path enveloping the body subject to the magnetic force, i.e., the
actuator plunger.
The above magnetic force calculation process is repeated as the plunger is stepped
through 20 equal 1 mm increments, which account for the full 20 mm plunger move-
ment. In such a way, the average magnetic force acting on the plunger is computed
over its entire range of travel.
Design Optimization
Here, the topology design optimization method is presented for the simultaneous
co-design of the actuator PMs, coils, and ferromagnetic material. To apply topology
optimization, the material properties in the design domain are interpolated using a
method similar to that proposed in [14, 81]. Four design variables, σ1 through σ4 ,
in each finite element control the material properties and allow for determination of
the PM, coils or ferromagnetic material portions of the design domain. The relative
permeability, μr , the residual magnetic flux density, Br , and the current density Jz
are defined as
p p
μr (σ1 ) = μr _ f err o σ1 1 + μr _air (1 − σ1 1 ) (5.128)
p
Br = σ2 2 Br _P M (5.129)
Jz = σ3 3 Jz∞ + σ4 3 (−Jz∞ )
p p
(5.130)
where p1 , p2 , and p3 are the penalization parameters (respectively, set to 0.5, 3 and
3), μr _ f err o and μr _air are the relative permeability of the ferromagnetic material
and air (respectively, set to 1500 and 1), Br _P M is the residual flux density of the
PM, and Jz∞ is the loaded external current density at the coil. Here, the design vari-
ables σ1 and σ2 control the relative permeability and the strength of the PM of each
finite element, respectively. Two design variables σ3 and σ4 are, respectively, used
to represent equal volume positive and negative direction current density coils.
m ,
The design goal is to maximize the average axial direction magnetic force, Ftot_x
acting on the plunger. Therefore, a objective function for the optimization problem
can be expressed as
1 ⎦ m
20
Fo = Ftot_x , (5.131)
20
m=1
where the objective function, Fo , is defined as the average of the axial direction
m , which can be calculated using Eq. (5.127) at 20 equally spaced plunger
force, Ftot_x
locations. The optimization problem is solved using the SLP method.
5.4 Actuator Analysis and Design 159
For this example, the sensitivity of the objective function, Fo , with respect to the
design variables, σ1 through σ4 , can be written as
1 ⎦ d Ftot_x
20 m
dFo
= (k = 1, ..., 4). (5.132)
dσk 20 dσk
m=1
d Fm
The sensitivity of the magnetic force at the m-th plunger location, dσtot_x
k
, is analyti-
m , with the
cally derived using the adjoint variable method. The magnetic force, Ftot_x
adjoint term may be written as
T
m
Ftot_x (Am ) + λ̂m [Km
M (σ1 )A − F M (σ1 , σ2 , σ3 , σ4 )],
m m
(5.133)
where λ̂m is the adjoint variable. The differentiation of Eq. (5.133) with respect to
σ1 gives the sensitivity of Ftot_x as
T λK M (σ1 ) m T λF M (σ1 , σ2 , σ3 , σ4 )
m
d Ftot_x m m
= λ̂m A − λ̂m . (5.134)
dσ1 λσ1 λσ1
The adjoint variable, λ̂m , is obtained by solving the adjoint equation derived as
T m
m (Am )T
λ Ftot_x
M λ̂ = −
Km . (5.135)
λAm
Likewise, the sensitivity of the objective function, Fo , with respect to σ2 -σ4 is, respec-
tively, derived as
T λF M (σ1 , σ2 , σ3 , σ4 )
m
d Ftot_x m
= −λ̂m (k = 2, 3, 4). (5.136)
dσk λσk
The adjoint equation needed to calculate λ̂m is derived following the prior procedure
used to arrive at Eq. (5.135); refer to [58] for additional details.
Design Result 1: Linear Actuator Plunger Design
The above analysis and optimization method is applied to the design of the PM
actuator model shown in Fig. 5.79. This actuator model is composed of a yoke with
coils and a plunger that is composed of ferromagnetic material and PMs. The shape
and location of the ferromagnetic material and PMs in the plunger is simultaneously
optimized. In this example, the magnetization direction of the PMs is fixed in a ver-
tical orientation. The material properties of the ferromagnetic material and PMs are
interpolated using two design variables, σ1 and σ2 ; refer to Eqs. (5.128) and (5.129).
The optimization goal is to maximize the average force during the 20 mm plunger
movement. Thus, the optimization problem may be stated as
160 5 Electromechanical System Simulation and Optimization Studies
Coil
Yoke
Plunger Magnetic
Force
Permanent
Magnet
Fig. 5.79 Conceptual diagram of a linear actuator, Design 1, that consists of a stationary yoke with
coil and a moving plunger with permanent magnets
Find σ1 and σ2
Maximize Eq. (5.131)
Subject to Eqs. (5.125) − (5.127)
σ dς − v1 ∗ 0
ςd 1 d
ςd σ2 dςd − v2 ∗ 0
0 ∗ σ1 ∗ 1
0 ∗ σ2 ∗ 1
Given Eqs. (5.128) − (5.129).
Here, the v1 and v2 are volume fraction for the ferromagnetic material and PM,
respectively. For the optimization problem, the SLP move limit is set to 0.015. The
volume constraints for the PMs and ferromagnetic material are, respectively, set to
5 % and 50 % of the design domain. An external current of 3 A flows through the
coil, which is assumed to have 480 turns, and the residual magnetic flux density of
the PM, Br _P M , is set to 0.4 T.
The 2-D design and analysis domain for the PM actuator is shown in Fig. 5.80. The
plunger movement is simulated by moving the plunger from d = 0 to d = 20 mm
in 20 equal, 1 mm, steps. Due to symmetry, one-half of the actuator is analyzed
using appropriate boundary conditions. In this example, PM actuators with two dif-
ferent configurations for the magnetic field source, i.e., the PM and external current
are studied; see Fig. 5.81. The figure shows two combinations of the magnetization
direction of the plunger PM and the direction of the current in the yoke coils. Model
1 has an anti-symmetric magnetic field source configuration, and thus one large mag-
netic field loop passes through all three components. For Model 1, a zero Neumann
condition is applied to the vector potential, A z , at the symmetry edge for zero tan-
gential (i.e., x-direction) magnetic flux density, Bx . On the other hand, Model 2 has
5.4 Actuator Analysis and Design 161
y
x
Yoke
Air Coil
Symmetry
Fig. 5.80 Two-dimensional design and analysis domain for linear actuator Design 1
X X X X X X X X X X X X X X
N N
S S
N S
S N
• • • • • • • X X X X X X X
Model 1 Model 2
Fig. 5.81 Two linear actuator Design 1 models assumed to have different magnetic field source
configurations (initial designs shown)
a symmetric field source configuration, which generates two symmetric and smaller
magnetic field loops. For Model 2, a zero Dirichlet condition is applied to the vector
potential, A z , at the symmetry edge for a zero normal (i.e., y-direction) magnetic
flux density, B y . In both models, the optimal shape and locations of the PMs and
ferromagnetic material in the plunger is found, and the average forces of the optimal
designs are compared.
The plunger design results for the two models are compared in Fig. 5.82. In each
case, the analysis model was discretized using approximately 9,000 linear quadri-
lateral elements, and a single iteration for the optimization loop required about 10 s
to solve on a quad-core workstation with a 3.4 GHz processor and 16 GB of RAM.
Observe in Fig. 5.82 that the optimal shape and location of the PMs (dark gray
structures) are almost identical. Also, a small air-gap region appears to the right of
162 5 Electromechanical System Simulation and Optimization Studies
• • • • • • • • • • • • • •
X X X X X X X X X X X X X X
N N
S S
N S
S N
X X X X X X X • • • • • • •
• • • • • • • X X X X X X X
Model 1 Model 2
Fig. 5.82 Comparison of optimal linear actuator Design 1 results for two different models. Note
that the dark gray regions indicate PM material, light gray regions indicate ferromagnetic material,
white regions indicate air, and magnetic equipotential lines are overlaid, for clarity
the PM regions in both models. The role of this air-gap is to act as a flux barrier
that increases the magnetic force by reducing the magnetic flux leakage that occurs
inside the plunger. Although the PMs and flux barrier optimization results are almost
identical in both models, the plunger ferromagnetic material (light gray structure)
result is different due to the dissimilar magnetic field source configurations. The
result for Model 1 shows a large cavity in the center of the plunger, but this cavity is
not present in the results Model 2. The different optimization results, however, gen-
erate an almost identical magnetic flux distribution around the air-gap, as shown in
Fig. 5.82. The magnetic field surrounding the air-gap mainly influences the magnetic
force magnitude. Therefore, the different plunger designs result in a nearly identical
optimal magnetic flux distribution around the air-gap, which results in a maximized
magnetic force.
The force profiles with respect to the plunger locations are compared in Fig. 5.83.
Table. 5.6 compares the average force from the initial design and optimized design.
Although the average forces of two optimized designs (285.3 and 284.9 N) vary by
less than 0.2 %, the initial design average forces (87.8 and 105.1 N) differ by 19.7 %.
This demonstrates that the magnetic field source configuration affects the magnetic
force in some structures, as expected. Lastly, in Model 1, the average force on the
plunger is significantly increased by 325 %, while in Model 2 the average force is
increased by 271 %.
5.4 Actuator Analysis and Design 163
Force (N)
Force (N)
Fig. 5.83 Comparison of optimized and initial average force versus plunger distance traveled for
linear actuator Design 1
Find σ1 − σ4
Maximize Eq. (5.131)
Subject to Eqs. (5.125) − (5.127)
σ dς − v1 ∗ 0, ςd σ2 dςd − v2 ∗ 0
ςd 1 d
ςd σ3 dςd − v3 ∗ 0, ςd σ4 dςd − v4 ∗ 0
0 ∗ σ1 ∗ 1, 0 ∗ σ2 ∗ 1
0 ∗ σ3 ∗ 1, 0 ∗ σ4 ∗ 1
Given Eqs. (5.128) − (5.130).
Here, the volume constraints for the PM, positive direction coil, negative direction
coil, and ferromagnetic material are set to 0.015, 0.025, 0.025, and 0.25, respectively.
164 5 Electromechanical System Simulation and Optimization Studies
Coil
Yoke Permanent
Magnet
Plunger Magnetic
Force
Fig. 5.84 Conceptual diagram of a linear actuator, Design 2, that consists of a stationary yoke with
coils and PMs plus a moving plunger
y
Yoke – Design Domain
x
- Permanent Magnet Coil
Fig. 5.85 Two-dimensional design domain and finite-element analysis domain with boundary con-
ditions for linear actuator Design 2. © 2011 IEEE. Reprinted with permission from [58, Fig. 3]
In Eq. (5.130), an external current of 3 A flows through the coil, which is assumed to
consist of 200 turns, and the residual magnetic flux density of the permanent magnet,
Br _P M , is similarly set to 0.4 T for Eq. (5.129).
The 2-D design and analysis domain of the Design 2 actuator is shown in Fig. 5.85.
The plunger movement is simulated by moving the plunger from d = 0 to d = 20 mm
in 20 equal, 1 mm, steps. Due to symmetry, one-half of the actuator is analyzed
using appropriate boundary conditions. In this example, a zero Neumann condition
is applied to the vector potential, A z , at the symmetry edge for the zero tangential (i.e.,
x-direction) magnetic flux density, Bx . Using the topology optimization method, the
optimal location and shape of the yoke components are found. In addition, the effect
of the PM magnetization direction and the external current density strength on the
optimization results is also investigated.
The initial and optimized results for linear actuator Design 2 are presented in
Fig. 5.86 on the top left and right, respectively. In this example, the model was
discretized using 9,600 linear quadrilateral elements, where one iteration of the opti-
mization routine took slightly longer, 12 s, to solve on a quad-core workstation with
a 3.4 GHz processor and 16 GB of RAM.
5.4 Actuator Analysis and Design 165
• • • • • • • • •
X X X X X X X X X
S S
X X • •
X X X • • •
N N
S S
X X X • • •
X X • •
N N
X X X X X X X X X
• • • • • • • • •
Fig. 5.86 Comparison of initial (top left) and optimal (top right) linear actuator Design 2 results
at d = 10 mm, and their corresponding force profiles versus plunger location, d, (bottom center).
Note that the dark gray regions indicate PM material, light gray regions indicate ferromagnetic
material, white regions indicate air, coil regions are shown graphically, and magnetic equipotential
lines are overlaid, for clarity. © 2011 IEEE. Reprinted with permission from [58, Fig. 4 and 5]
In the optimal structure in Fig. 5.86, both the PMs (dark regions) and the coils
appear near the air-gap. The optimal structure maximizes the magnetic force by
reducing the leakage magnetic field outside of the air-gap and plunger, which can
be observed by comparing the magnetic flux lines of the initial and optimal designs
in Fig. 5.86. Note that the magnetic field surrounding the air-gap mainly influences
the magnetic force magnitude. The magnetic force profile with the plunger location,
d, and the average force is also compared in Fig. 5.86 (bottom center image). The
average force of the optimal design is expected to be 41.9 % higher than that of the
initial design.
In addition, the influence of the PM magnetization direction on the optimization
result is investigated here. The optimization process was repeated for four differ-
ent magnetization directions (45⊥ , 90⊥ , 135⊥ , and 180⊥ with respect to the y-axis).
Figure 5.87 shows the optimization results with magnetic flux lines for each mag-
netization direction. In every result, the coils appear near the air-gap to minimize
leakage, and the PMs appear such that their magnetization directions are aligned
with the direction of the magnetic field generated by the coil and passing through
166 5 Electromechanical System Simulation and Optimization Studies
Fig. 5.87 Linear actuator Design 2 optimization results with magnetic flux lines (at d=10 mm) for
different PM magnetization directions (first and second rows) and magnetic force profile comparison
(bottom center). Note that the dark gray regions indicate PM material, light gray regions indicate
ferromagnetic material, white regions indicate air, coil regions are shown graphically, and magnetic
equipotential lines are overlaid, for clarity. © 2011 IEEE. Reprinted, with permission, from [58,
Fig. 6]
the ferromagnetic material. The magnetic force profile and average force for each
actuator design is further compared in Fig. 5.87 (bottom center image). When the
magnetization direction is 180⊥ , the PMs are located closest to the air-gap and con-
sequently, the smallest leakage occurs along with the highest average force on the
plunger.
5.4 Actuator Analysis and Design 167
Fig. 5.88 Linear actuator Design 2 optimization results with magnetic flux lines (at d=10 mm) for
a fixed PM strength and varying external current density strengths. Note that the dark gray regions
indicate PM material, light gray regions indicate ferromagnetic material, white regions indicate air,
coil regions are shown graphically, and magnetic equipotential lines are overlaid, for clarity. © 2011
IEEE. Reprinted with permission from [58, Fig. 7]
Finally, the effect of the external current strength on the optimization result is
investigated. It is hypothesized that the ratio of the amount of PM to the external
current density strength affects the optimization result. Thus, in order to investigate
this effect, the optimization process is performed using four different current density
strengths (i.e., 1 A, 3 A, 5 A, and 7 A) with a fixed PM strength of 0.4 T. Here, the
magnetization direction is held fixed at 90⊥ with respect to the y-axis. Figure 5.88
shows the various optimization results with magnetic flux lines. As the external coil
current strength increases, the optimization process tends to generate a structure with
smaller magnetic reluctance in order to maximize the magnetic field generated by the
coils. Typically, a PM has high magnetic reluctance since its permeability is almost
the same as that of air. Compared to the optimization result with 1 A, the structure
with 3 A contains a thin PM structure, which has smaller magnetic reluctance despite
higher leakage of the PM magnetic field. When the external current increases to 5
A, a thick ferromagnetic structure is formed to accommodate the magnetic field
generated by the coil, and a thin additional structure is created for the PM magnetic
168 5 Electromechanical System Simulation and Optimization Studies
Rotor
field. Even for the highest current case (7 A), the PM is a minor magnetic field source,
and consequently, the ferromagnetic material is used only to minimize the magnetic
reluctance associated with the coil source.
In this section, two different types of electric motors are considered including:
(1) switched reluctance motors and (2) interior permanent magnet motors. Here,
switched reluctance motors (SRMs) are first considered as an alternative to interior
permanent magnet motors due to their advantage of not requiring expensive rare-earth
materials in their construction. However, one well-known challenge with SRMs is
the torque ripple associated with their operation, and the minimization of this torque
ripple for an SRM is considered as a design objective in Sect. 5.5.1. In contrast, while
scarcity and supply chain issues drives the cost of rare-earth materials up, an interior
permanent magnet motor exhibits a significantly higher torque to motor volume ratio
(or torque density) due to the enhanced magnetic properties of the associated per-
manent magnets [31]; such motors are found in many applications where increased
performance in a compact package is required. Thus, for completeness, the multi-
physics analysis of a more traditional interior permanent magnet motor is presented
in Sect. 5.5.2.
The following SRM design optimization study utilizes a 2-D model, as shown in
Fig. 5.89 [55]. The design goal is to minimize the previous described torque ripple
effect, which again is the main disadvantage of a switched reluctance motor. To
5.5 Electric Motor Analysis and Design 169
achieve the design goal, not only the rotor/stator geometry of the motor but the electric
variables (voltage on/off angles) are optimally designed. In addition, the copper loss
is restrained by introducing a constraint on the root-mean-square (RMS) value of the
phase current. The optimization problem is formulated for the above objective and
constraint. The optimization problem is solved by using an SLP method, and then
the optimal material distribution of the rotor/stator and the optimal voltage on–off
angles can be found. The 2-D model of a 6/4 SRM (6 stator poles and 4 rotor poles)
is chosen in the following design example. The analysis model takes the magnetic
saturation effect into consideration, which means that the magnetic permeability is
dependent on the magnetic flux density.
Performance Analysis
Motor performance is explicitly represented using a mathematical approximation
method. This analytical representation of the performance is derived for the analyt-
ical sensitivity, which is required for gradient-based optimization. The performance
analysis is conducted using the following procedure. Step 1: the flux linkage, λlink ,
is approximated with respect to current, i, and rotor angle, υr , by using Fourier series
expansions and piecewise quadratic polynomials. Step 2: the phase current curve is
obtained by solving the voltage equation of an equivalent circuit. Step 3: the torque
profile is determined using the global virtual work method.
Step 1—Flux Linkage Model: The flux linkage model represents the relation
between the flux linkage, λlink , the current of the stator coil, i, and the rotor posi-
tion, υr . To calculate the flux linkage, the magnetic energy, Wmag , at discrete current
values and rotor positions is calculated by solving the governing equations for low-
frequency electromagnetics with magnetic saturation effect, which is explained in
Sect. 3.5. This non-linear equation is solved using the finite element method with
Newton–Raphson iteration. From the equation solution, the magnetic energy, Wmag ,
can be obtained as ⎧ ⎨⎧ B ⎩
Wmag = H · dB dv, (5.137)
allspace 0
where H is magnetic field intensity and B is magnetic flux density. This process
is performed at various phase current values and rotor positions, and the magnetic
energy, (Wmag ) jk , at a discrete j-th phase current and k-th rotor position is calculated.
The next step is to find the flux linkage curves with respect to the phase current,
i, axis at fixed rotor position, υr . From the approximation using piecewise quadratic
polynomials, the flux linkage, (λlink )km , at a fixed k-th rotor position and m-th piece-
wise interval can be written as
where C1km –C3km are the polynomial coefficients. The reason why piecewise
quadratic polynomials are chosen is as follows. First, any piecewise function may
reasonably represent both a line and a curve. The second-order polynomial gives the
analytical solution of the voltage equation. Finally, the derivative of the flux linkage,
170 5 Electromechanical System Simulation and Optimization Studies
λlink , is not required for the performance analysis. Therefore, the C 1 continuity that
can be considered in a cubic polynomial does not need to be satisfied. Figure 5.90
shows the flux linkage curve with respect to the phase current at a fixed rotor angle.
To find the three coefficients, C1km –C3km , three conditions are required at each k-th
rotor position and m-th piecewise interval. Two conditions are obtained from the rela-
tion between the flux linkage, (λlink )km , and the magnetic energy, (Wmag ) jk , which
can be represented as
⎧
ij
(Wmag ) jk = i j · (λlink )km i j − (λlink )km idi (5.139)
0
In the above equation, the discrete current position, j, is set as 2m − 1, and 2m. The
last condition is the continuous condition at the current position, j = 2m + 1. From
these three conditions, the system of equations to find the coefficients C1 jm –C3 jm
is derived as
2 3
3 i 2m + 3 i 2m−1 2 i 2m + 2 i 2m−1 i 2m−1
1 3 1 2 1 2
C1 jm A1
⎤ 2i3 × ⎤ C2 jm = ⎤ A2 , (5.140)
3 2m+1 + 3 i 2m−1 2 i 2m+1 + 2 i 2m−1 i 2m−1
1 3 1 2 1 2
2
i 2m−1 i 2m−1 1 C3 jm A3
where
A1 = (Wmag ) j (2m)
⎦
m−1
1
+ ( C1 jk i 2k+1
3
− i 2k−1
3
3
k=1
1
+ C2 jk i 2k+1
2
− i 2k−1
2
+ C2 jk i 2k+1 − i 2k−1 ), (5.141)
2
A2 = (Wmag ) j (2m−1)
5.5 Electric Motor Analysis and Design 171
⎦
m−1
1
+ ( C1 jk i 2k+1
3
− i 2k−1
3
3
k=1
1
+ C2 jk i 2k+1 2
− i 2k−1
2
+ C2 jk i 2k+1 − i 2k−1 ), (5.142)
2
A3 = C1 j (m−1) i 2m−1
2
+ C2 jm i 2m−1 + C3 jm . (5.143)
Next, a Fourier series expansion is applied to Eq. (5.138) for the approximation with
respect to the rotor position, υr , axis. Then, the model representing the flux linkage,
λlink , with respect to phase current, i, and rotor angle, υr , is represented as
1 ⎦
m
Fqm,0 = (Cq( j+1)m + Cq( j)m )(υr ( j+1) − υr ( j) ) , (5.145)
2∂/Pr
j=1
m ⎪
2 ⎦ cos(n Pr υr ( j+1) ) − cos(n Pr υr ( j) )
Fqm,n = Cq( j+1)m · sin(n Pr υr ( j+1) ) +
n∂ n Pr (υr ( j+1) − υr ( j) )
j=1
cos(n Pr υr ( j+1) ) − cos(n Pr υr ( j) )
− Cq( j)m · sin(n Pr υr ( j) ) + . (5.146)
n Pr (υr ( j+1) − υr ( j) )
Figure 5.91 shows the attained flux linkage model. The filled circular marks represent
the flux linkage at discrete current and rotor angle positions, where the finite element
analysis is carried out. As shown in Fig. 5.91, the flux linkage is smoothly approx-
imated using the proposed method for the analytical representation with respect to
current, i, and rotor angle, υr .
Step 2—Phase Current Curve: After calculating the flux linkage, the current curve
with respect to rotor angle is obtained by solving the voltage equation. The voltage
equation of SRMs is given by
d(λlink (i, υr )
V = Ri + , (5.147)
dt
where V is the source voltage, and R is the stator coil resistance.
172 5 Electromechanical System Simulation and Optimization Studies
(θ r )
(θ r )
In order to solve the voltage equation analytically, the voltage drop due to the stator
coil resistance is neglected. This assumption might cause errors in the performance
analysis result. However, the voltage drop due to the resistance might be much less
than the back electromotive force term by the time variation of the flux linkage.
Thus, the error caused by the assumption is expected to be acceptable. Additionally,
we assume steady-state rotation and zero phase current at the voltage-on angle. The
voltage waveform used in this work is shown in Fig. 5.92. By solving Eq. (5.147)
with the prior assumptions, the explicit expression of phase current curve is obtained
as )
where
⎦
NF
Xi = Fim,n cos(n Pr υr ), (5.149)
n=0
5.5 Electric Motor Analysis and Design 173
(θ r )
υr (0) is the voltage-on angle, λlink(0) is the flux linkage at υr (0) , and ω is the constant
angular velocity [55]. Figure 5.93 shows the phase current profiles of typical 6/4
SRMs assuming non-linearly saturated material. In this example, typical rotor and
stator shapes, as shown in Fig. 5.89, are used.
Step 3—Torque Curve: The torque profile of SRMs is calculated using the global
virtual work method. The global virtual work method is based on the principle of
conservation of energy and virtual displacement. This method uses the co-energy at
a set of closely spaced rotor positions. The co-energy, Wco , is defined as
⎧
i
Wco = λlink d ĩ . (5.150)
0 υr =constant
By taking the derivative of the total co-energy with respect to the rotor angle, υr , the
torque, T , profile is obtained as
λ Wco
T (i, υr ) = . (5.151)
λυr i=constant
Substituting Eqs. (5.144) and (5.150) into Eq. (5.151), the torque profile is explicitly
given by
NF
⎦ 1 1
T (i, υr ) = − F1m,n (i 3 − i 2m−1
3
)+ F2m,n (i 2 − i 2m−1
2
)
3 2
n=0
+ F3m,n (i − i 2m−1 ) n Pr sin(n Pr υr )
⎦⎦
m−1 NF
1 1
− F1k,n (i 2k+1
3
− i 2k−1
3
) + F2k,n (i 2k+1
2
− i 2k−1
2
)
3 2
k=1 n=0
+ F3k,n (i 2k+1 − i 2k−1 ) n Pr sin(n Pr υr ). (5.152)
174 5 Electromechanical System Simulation and Optimization Studies
(θ r )
It is noted that the function for the current, i, is obtained from Eq. (5.148). From
Eq. (5.152), we can obtain the single-phase torque profile of SRMs (i.e., the torque
profile when the voltage of a single pair of the stator poles is switched on and off).
Since one phase is overlapped with the previous or the next phase in SRMs, the two
consecutive single-phase torque profiles are summed in the overlapped range. After
summing the torque profiles in all phases, the total torque profile is determined.
The overlapped range depends on the number of rotor poles, Pr , and stator poles,
Ps . The rotor pole repeats its position after a rotor period of 2∂/Pr . During this
period, the rotor pole passes Ps /2 stator poles. Therefore, the single-pulse torque
profile is overlapped every Ptorque calculated as
2∂/Pr 4∂
Ptorque = = . (5.153)
Ps /2 Pr Ps
In a steady rotation, the torque profiles of all phases are identical. Therefore, the total
torque, Ttot , is given by
⎦r
2∂/P
Ttot (υr ) = T (υr + (m − 1)Ptorque ). (5.154)
m=1
Figure 5.94 shows three consecutive single-phase torque profiles and the total torque
profile of the SRM model shown in Fig. 5.89. In this analysis, the voltage on-off
angles are set to 0⊥ and 50⊥ , respectively. Single-phase torque exists between these
two angles, and repeats every 30⊥ (i.e., Ptorque is 30⊥ ). As illustrated in Fig. 5.94,
SRMs have inherently large torque ripples.
Optimization Model and Results
As we have discussed throughout this book, structural topology optimization allows
one to find the optimal arrangement of a structure by setting material densities as
5.5 Electric Motor Analysis and Design 175
design variables. The material densities are relaxed as continuous variables in order
to use a mathematical programming method such as SLP. The material becomes air
when the density is zero, and it becomes steel when the density is one. The density
distribution of a design domain is optimized to achieve specific design goals. Here,
the design goals are to minimize torque ripple and the mass of the rotor
n ⎨
⎦ ⎩2
Tk
Fo = − 1 + w1 Fc (σ ) + w2 v(σ ), (5.155)
T∞
k=1
while the following constraint is applied to the RMS value of the current curve:
⎛ , ⎞
-
-1 ⎦n
ir ms ⎝≈ . i k2 dt ⎠ ∗ ir∞ms , (5.156)
T
k=1
where
⎦
ND
σu (1 − σu )
Fc (σ ) = , (5.157)
ND
u=1
1
σu A u
v(σ ) = 1u , (5.158)
u Au
Tk is the torque value at k-th rotor angle υr (k) (k = 1, . . . , n), n is the number of
discrete rotor angles, T ∞ is the target average torque, ir ms is the RMS value of phase
current, Au is the area of u-th element, w1 and w2 are the weighting values for the
density convergence function, Fc (σ ), and the normalized rotor mass function, v(σ ),
respectively, and N D is the number of design variables.
To minimize the torque ripple, the torque, Tk , at discrete rotor angles, υr (k) , is
handled to become the constant target average torque, T ∞ . The density convergence
function, Fc (σ ), is added to the objective function to enforce the convergence of
intermediate densities to zero or one, and the normalized rotor mass function, v(σ ),
is added to minimize the mass of the rotor as an additional objective function. In
the constraints, the RMS value of current curve is restricted to ir∞ms to confine the
copper loss, which is one of main heat sources, i.e., (Pcopper = Rir2ms ). Thus, the
optimization problem can be stated as
Find σu
Minimize Eq. (5.155) with Eqs. (5.157) − (5.158)
Subject to Eq. (5.156)
ςd σu dςd − v ∗ 0
0 ∗ σu ∗ 1.
176 5 Electromechanical System Simulation and Optimization Studies
Stator
Periodic Domain
Rotor
Periodic Domain
The structural design domain is shown in Fig. 5.95 [57]. Since a 6/4 SRM is
chosen, a 1/12 part of the stator and a 1/8 part of the rotor are considered as the geo-
metric design domain. Then, the material densities in the geometric design domain
are reflected and copied to the corresponding symmetric and periodic parts of the
rotor and stator. Using this approach, the symmetries in the rotor and stator of the
considered SRM are effectively exploited.
The formulated optimization problem is solved using the SLP method requir-
ing the sensitivity of objective and constraint functions. The proposed performance
analysis model gives the analytical representation of the current curve and torque
profile. Therefore, the analytical sensitivity of objective and constraint functions with
respect to design variables, σu , is available.
The sensitivity is analytically derived using direct differentiation and the chain
rule. The sensitivity of the objective function, Fo , in Eq. (5.155) with respect to the
design variable, σu , is given as
⎦
n ⎨ ⎩
dFo Tk λTk 1 − 2σu
= 2 − 1.0 + w1 + w2 Au . (5.159)
dσu T∞ λσu ND
j=1
A 2-D 6/4 SRM is designed using the proposed optimization method. The target
average torque, T ∞ , is set as the average torque of the typical design, which is
shown in Fig. 5.89. The constraint value, ir∞ms , on the RMS value of phase current is
also set as the RMS current of the typical design. The move limits of the SLP method
are 0.004 for the density design variables and 0.0004 for the angle design variables.
The grayscale structure obtained as the optimization result is presented in Fig. 5.96.
The dark colored regions represent steel, and the light colored regions represent air.
The 2-D SRM computational model was discretized using 20,672 linear quadrilateral
elements, and a single iteration of the optimization routine required roughly 9 s when
solving the problem on a quad-core workstation with a 3.4 GHz processor and 8 GB
of RAM.
5.5 Electric Motor Analysis and Design 177
Fig. 5.96 Optimal topology of a 6/4 SRM rotor and stator. Note that dark colored regions indicate
steel, while light colored regions indicate air; refer to Fig. 5.89 for further descriptions. © 2010 Lee
J. Reprinted, with permission, from [55, Fig. 6.4.2c]
As shown in Fig. 5.96, the density design variables have clearly converged to
zero (void) or one (steel) by applying the density convergence function, Fc (σ ), in
Eq. (5.157). The rotor and stator poles are optimized to have a notched shape near
the air-gap. The notched shape of rotor and stator poles manipulates the inductance
profile, and consequently, minimizes the SRM torque ripples.
Table 5.7 shows the design parameters of the typical and optimized SRM and also
compares the average torques and the RMS values of the current curve. As shown in
Table 5.7, the arc lengths of both the rotor and the stator poles are increased in order
to minimize torque ripples. The range of the voltage on–off angles is also widened.
The average torque of the optimized motor successfully satisfies the target average
torque (e.g. 6 % lower than target average torque). The RMS value of the current
curve in the designed motor is well constrained.
In Fig. 5.97, the current curve and total torque profile of the optimized design
(solid lines) are compared to those of the typical design (dashed lines). As shown
in Fig. 5.97a, the operating ranges of phase current increase with a widened voltage
waveform. The narrow current curve of the typical model gets spread in the optimized
model while satisfying the current constraint. As shown in Fig. 5.97b, the total torque
profile of the optimized design is nearly matched with the target torque, and the torque
ripple is remarkably reduced.
178 5 Electromechanical System Simulation and Optimization Studies
Table 5.7 Design parameters and performance analysis result for SRM motors
Typical shape Optimized shape
Stator pole arc 27⊥ Stator pole arc 39⊥
Rotor pole arc 42⊥ Rotor pole arc 45⊥
Voltage-on angle 0⊥ Voltage-on angle −4.43⊥
Voltage-off angle 50⊥ Voltage-off angle 63.50⊥
Average torque 184.7 N m Average torque 172.9 N m
RMS of current 57.2 A RMS of current 56.4 A
( )
tot)
(
( )
Fig. 5.97 Performance comparison of the typical 6/4 SRM shown in Fig. 5.89 and optimized 6/4
SRM shown in Fig. 5.96: a Current curve, b total torque profile. © 2010 Lee J. Reprinted with
permission from [55, Fig. 6.4.3a, b]
The interior permanent magnet synchronous motor (IPMSM) has been the most
promising candidate for the powertrain of hybrid and electric vehicles due its high
power density and efficiency. In this final simulation example of the text, the mul-
tiphysics analysis of an IPMSM is presented. The geometry and material properties
of the IPMSM used in this study is chosen based on the IPMSM found in a typi-
5.5 Electric Motor Analysis and Design 179
Permanent Magnet
cal commercial hybrid vehicle [32, 78]. Here, the modeling of the IPMSM is first
explained, and then a coupled electromagnetic/heat transfer analysis is performed to
predict the torque and thermal performance of the IPMSM.
Modeling of IPMSM
In order to carry out the multiphysics simulation of an IPMSM, an analysis model
was prepared using the commercial software COMSOL®5 [17]. The geometric para-
meters and material properties of the IPMSM components are chosen based on the
data presented in [32, 78]. Figure 5.98 shows the assumed 2-D geometry of the
IPMSM. The motor consists of the 3-phase distributed windings, a 48-pole stator,
and an 8-pole rotor with 16 NdFeB PM pieces. The axis in the middle of the rotor PM
poles is called the direct-axis, or simply d-axis, and the intermediate axis between
d-axes is called the quadrature axis, or simply q-axis. The q-axis is leading the d-axis
by 90⊥ electrical degrees.
The assumed geometric parameters of the IPMSM model are summarized in
Table 5.8. The configuration of the three-phase winding is presented in Fig. 5.99.
The phase A winding is located at the outermost position, and the phase C winding is
(a) A+
(b) (c)
B- C-
B+ C+
A- A-
B+ C+
B- C-
A+ A+
B- C-
B+ C+
A- A-
B+ C+
A+ B- C-
located at the innermost position closest to the rotor. As mentioned above, there are
48 winding slots, and each slot area is divided into three for a three-phase distributed
winding. The winding for each phase is excited by an AC current with different phase
angles. The current for each phase is modeled as
where I peak is the amplitude of the AC current, and f is the frequency of the AC
current, which is synchronous with the rotating speed of the rotor. When applying
the above currents to each phase, the eight poles of the stator electromagnet rotates
proportionally to the frequency of the AC current. The wire size is set as American
Wire Gauge (AWG) 20, which has a diameter of 0.812 mm. The winding area for
one phase in one slot is approximately 50 mm2 . In this example, the winding is
modeled as 12 parallel bundled wires, and 3 bundled wires are positioned inside the
area designated for one phase in one slot. When applying this wire configuration, the
RMS current density is calculated as 10.2 A/m2 when the phase current is 100 RMS
Ampere.
The material physical parameters for the PM, iron, and coils are summarized
in Table 5.9. Due to the magnetic saturation effect of the iron material, the relative
permeability, μr , of the iron material gradually decreases as the magnetic flux density,
B, increases. For the electromagnetic analysis, a quasistatic assumption is adopted for
simplicity, and thus electrical conductivity, ρ , and permittivity, , material physical
parameters are not used in this example. For the thermal analysis, only steady-state
heat conduction in solid materials is considered, and thus, the density, γ, and heat
capacity, C p , material physical parameters are also not used.
5.5 Electric Motor Analysis and Design 181
Table 5.9 Material physical parameters for the various components in the assumed IPMSM
Material Relative permeability (μr ) Thermal conductivity (k)
Iron in stator and rotor 1200 (when B is less than 1 T) 20 W/(m K)
PM in rotor 1 (Constant) 5 W/(m K)
Copper in winding 1 (Constant) 400 W/(m K)
Multiphysics Analysis
A coupled electromagnetic and thermal analysis is performed using the aforemen-
tioned commercial FEA software [17]. The 2-D IPMSM model was discretized using
75,658 second order triangular elements. Here, the single computational analysis
required for the non-linear magnetostatic-thermal problem required approximately
60 seconds to solve on a quad-core workstation with a 3.4 GHz processor and 8 GB
of RAM.
In the solution process, a low-frequency electromagnetic field analysis is first
performed by solving the magnetostatic governing equations; refer to Sect. 3.5. A
synchronous motor rotates by maintaining the electrical degrees between rotor poles
and stator poles. Figure 5.100a, b allow for a comparison of the magnetic field
distribution when the electrical degrees between the rotor PM pole and stator pole
are 0⊥ and 90⊥ , respectively. As can be seen in Fig. 5.100, the electrical degrees
are 0⊥ when the quadrature axis (q-axis) current is a maximum and the direct axis
(d-axis) current is zero, and they are 90⊥ when the q-axis current is zero and the
d-axis current is a maximum.
From the magnetic field distribution, the torque applied to the rotor is calculated
using the Maxwell Stress Tensor method. Figure 5.101 presents the torque curves for
the motor with respect to electrical degrees. The total torque of an IPMSM motor is
composed of two components, which include the reluctance torque and PM torque.
Reluctance torque is generated due to the interaction between the stator electromagnet
and rotor iron. The magnetic reluctance of the rotor varies in different rotor positions
because of the uneven shape of the rotor iron, which generates the torque. As shown
in Fig. 5.101, the reluctance torque is zero when the stator current is matched with
either the q-axis (i.e., electrical degrees = 0⊥ ) or d-axis (i.e., electrical degrees = 90⊥ ).
PM torque is generated due to the interaction between the stator electromagnet and
the rotor PM. PM torque is maximized when the stator current along the q-axis is
maximized (i.e., electrical degrees = 0⊥ ), and it becomes zero when the stator current
is aligned along the d-axis (i.e., electrical degrees = 90⊥ ). Because the periods of the
torque components are different, the electrical angle maximizing the total torque is
located between 0⊥ and 90⊥ depending on the amplitude of each torque term.
Next, the thermal analysis of the IPMSM is performed, where heat is generated
by the coil loss in the winding plus core loss in the stator iron. In the steady-state
operation of a synchronous motor, the magnetic field in the rotor winding does not
vary with respect to time, and thus the core loss in the rotor iron is ignored in this
example. The coil loss, Pcoil , in the winding can be calculated as Joule loss
182 5 Electromechanical System Simulation and Optimization Studies
(a) (b)
Fig. 5.100 Two-dimensional low-frequency electromagnetic field analysis results for an IPMSM
with 48 stator poles and 8 rotor poles. Note that lighter colored contours indicate a higher current
density, and the magnetic equipotential lines are overlaid, for clarity
Fig. 5.101 Torque curves for
an IPMSM with respect to
electrical angles
⎧
Pcoil = γe J 2 dv, (5.163)
Physteresis = kh B 2 f, (5.165)
where ρ is the electrical conductivity of the iron (1×107 S/m), d is the thickness of the
laminated iron (0.305 mm), kh is the hysteresis coefficient (set as 143 Ws/T2 /m3 ), and
f is the fundamental frequency (266.7 Hz corresponding to 2000 RPM). A convective
cooling condition is applied at the motor outer boundary, where the heat transfer
coefficient is set as 50 W/(m2 K) with a reference temperature of 293.2 K. The
distribution of the temperature is shown in Fig. 5.102, and the maximum temperature,
363.8 K, is expected at the end of stator pole near the air-gap.
References
1. Aage N, Mortensen NA, Sigmund O (2010) Topology optimization of metallic devices for
microwave applications. Int J Numer Meth Eng 83:228–248. doi:10.1002/nme.2837
2. Alexandersen J, Andreasen CS, Aage N, Lazarov BS, Sigmund O (2013) Topology optimisation
for coupled convection problems. 10th world Congress on structural and multidisciplinary
optimization, Orlando, 19–24 May 2013
3. Andkjær J, Nishiwaki S, Nomura T, Sigmund O (2010) Topology optimization of grating
couplers for the efficient excitation of surface plasmons. J Opt Soc Am B 27:1828–1832.
doi:10.1364/JOSAB.27.001828
184 5 Electromechanical System Simulation and Optimization Studies
4. Andreasen CS, Gersborg AR, Sigmund O (2008) Topology optimization for microfluidic mix-
ers. Int J Numer Meth Fl 61:498–513. doi:10.1002/fld.1964
5. Bendsøe MP (1989) Optimal shape design as a material distribution problem. Struct Multidiscip
O 1:193–202. doi:10.1007/BF01650949
6. Bendsøe MP, Sigmund O (2003) Topology optimization: theory, methods, and applications,
2nd edn. Springer, Berlin
7. Berenger J (1994) A perfectly matched layer for the absorption of electromagnetic waves. J
Comput Phys 114:185–200. doi:10.1006/jcph.1994.1159
8. Borvall T, Petersson J (2003) Topology optimization of fluids in stokes flow. Int J Numer Meth
Fl 41:77–107. doi:10.1002/fld.426
9. Boulware JC, Jensen S (2010) Thermomagnetic siphoning on a bundle of current-carrying
wires. In: Proceedings of the COMSOL conference 2010, Boston, 7–9 Oct 2010
10. Braithwaite D, Beaugnon E, Tournier R (1991) Magnetically controlled convection in para-
magnetic fluid. Nature 354:134–136. doi:10.1038/354134a0
11. Bruggi M, Cinquini C (2011) Topology optimization for thermal insulation: an application to
building engineering. Eng Optimiz 43:1223–1242. doi:10.1080/0305215X.2010.550284
12. Choi HS, Park IH, Lee SH (2006) Electromagnetic body force calculation based on virtual air
gap. J Appl Phys 99:08H903. doi:10.1063/1.2173207
13. Choi HS, Park IH, Lee SH (2006) Concept of virtual air gap and its applications for force
calculation. IEEE T Magn 42:663–666. doi:10.1109/TMAG.2006.871594
14. Choi JS, Yoo J (2009) Simultaneous structural topology optimization of electromagnetic
sources and ferromagnetic materials. Comput Method Appl M 198:2111–2121. doi:10.1016/
j.cma.2009.02.015
15. COMSOL AB (2008) COMSOL Multiphysics, Ver. 3.5a. Stockholm
16. COMSOL AB (2010) COMSOL Multiphysics, Ver. 4.0a. Stockholm
17. COMSOL AB (2011) COMSOL Multiphysics, Ver. 4.2. Stockholm
18. Cooper M, Petosa A, Wight JS, Ittipiboon A (1996) Investigation of dielectric resonator anten-
nas for L-band communications. Paper presented at the antenna technology and applied elec-
tromagnetics symposium, ANTEM ‘96
19. Dede EM (2009) Multiphysics topology optimization of heat transfer and fluid flow systems.
In: Proceedings of the COMSOL conference 2009, Boston, 8–10 Oct 2009
20. Dede EM (2010) Multiphysics optimization, synthesis, and application of jet impingement
target surfaces. In: Proceedings of the 12th IEEE intersociety conference on thermal and ther-
momechanical phenomena in electronic systems, Las Vegas, 2–5 June 2010. doi:10.1109/
ITHERM.2010.5501408
21. Dede EM (2010) The influence of channel aspect ratio on performance of optimized thermal-
fluid structures. In: Proceedings of the COMSOL conference 2010, Boston, 7–9 Oct 2010
22. Dede EM (2010) Simulation and optimization of heat flow via anisotropic material thermal
conductivity. Comp Mater Sci 50:510–515. doi:10.1016/j.commatsci.2010.09.012
23. Dede EM (2011) Experimental investigation of the thermal performance of a manifold hier-
archical microchannel cold plate. In: Proceedings of the ASME 2011 Pacific Rim technical
conference and exhibition on packaging and integration of electronic and photonic systems,
MEMS and NEMS, vol 2. Portland, 6–8 July 2011. doi:10.1115/IPACK2011-52023
24. Dede EM (2012) Optimization and design of a multipass branching microchannel heat sink for
electronics cooling. J Electron Packaging 134:041001. doi:10.1115/1.4007159
25. Dede EM (2012) Jet impingement heat exchanger apparatuses and power electronics modules.
US Patent, 8,199,505 B2
26. Dede EM (2013) Power electronics modules and power electronics module assemblies. US
Patent, 8,391,008 B2
27. Dede EM, Liu Y (2013) Experimental and numerical investigation of a multi-pass branching
microchannel heat sink. Appl Therm Eng 55:51–60. doi:10.1016/j.applthermaleng.2013.02.
038
28. Dede EM, Liu Y (2013) Cold plate assemblies and power electronics modules. US Patent,
8,427,832 B2
References 185
29. Dede EM, Nomura T, Schmalenberg P, Lee JS (2013) Heat flux cloaking, focusing, and
reversal in ultra-thin composites considering conduction-convection effects. Appl Phys Lett
103:063501. doi:10.1063/1.4816775
30. Dede EM, Nomura T, Lee J (2014) Thermal-composite design optimization for heat flux shield-
ing, focusing, and reversal. Struct Multidiscip O 49:59–68. doi:10.1007/s00158-013-0963-0
31. Dent PC (2012) Rare earth elements and permanent magnets (invited). J Appl Phys 111:07A721.
doi:10.1063/1.3676616
32. Dorrell DG, Knight AM, Evans L, Popescu M (2012) Analysis and design techniques applied to
hybrid vehicle drive machines—assessment of alternative ipm and induction motor topologies.
IEEE T Ind Electron 59:3690–3699. doi:10.1109/TIE.2011.2165460
33. Erentok A, Sigmund O (2008) Three-dimensional topology optimized electrically-small con-
formal antenna. Paper presented at the IEEE antennas and propagation society international
symposium, San Diego, 5–12 July 2008. doi:10.1109/APS.2008.4619200
34. Eshelby JD (1957) The determination of the elastic field of an ellipsoidal inclusion, and related
problems. P R Soc A 241:376–396. doi:10.1098/rspa.1957.0133
35. Fan Z, Antar Y, Ittipiboon A, Petosa A (1996) Parasitic coplanar three-element dielectric
resonator antenna subarray. Electron Lett 32:789–790. doi:10.1049/el:19960514
36. Finlaysan BA (1970) Convective instability of ferromagnetic fluids. J Fluid Mech 40:753–767.
doi:10.1017/S0022112070000423
37. Garimella SV, Singhal V (2004) Single-phase flow and heat transport and pumping
considerations in microchannel heat sinks. Heat Transfer Eng 25:15–25. doi:10.1080/
01457630490248241
38. Gelet J-L, Gerlaud A (2013) Control of Joule heating extends performance and device life.
IEEE Spectrum 5:S-24–25
39. Halbach K (1980) Design of permanent multipole magnets with oriented rare earth cobalt
material. Nucl Instrum Methods 169:1–10. doi:10.1016/0029-554X(80)90094-4
40. Harpole GM, Eninger JE (1991) Micro-channel heat exchanger optimization. In: Proceedings
of the 7th IEEE SEMI-THERM symposium, Phoenix, 12–14 February 1991. doi:10.1109/
STHERM.1991.152913
41. Hasselman DPH, Bhatt H, Donaldson KY, Thomas JR (1992) Effect of fiber orientation and
sample geometry on the effective thermal conductivity of a uniaxial carbon fiber reinforced
glass matrix composite. J Compos Mater 26:2278–2288. doi:10.1177/002199839202601506
42. Hatta H, Minoru T (1986) Equivalent inclusion method for steady state conduction in compos-
ites. Int J Eng Sci 24:1159–1172. doi:10.1016/0020-7225(86)90011-X
43. Hull D, Clyne TW (1996) An introduction to composite materials, 2nd edn. Cambridge Uni-
versity Press, Cambridge
44. Iga A, Nishiwaki S, Izui K, Yoshimura M (2009) Topology optimization for thermal conductors
considering design-dependent effects, including heat conduction and convection. Int J Heat
Mass Tran 52:2721–2732. doi:10.1016/j.ijheatmasstransfer.2008.12.013
45. Incropera FP (1999) Liquid cooling of electronic devices by single-phase convection. Wiley-
Interscience, New York
46. Incropera FP, DeWitt DP, Bergman TL, Lavine AS (2007) Introduction to heat transfer, 5th
edn. Wiley, Hoboken
47. Jensen JS, Sigmund O (2004) Systematic design of photonic crystal structures using topology
optimization: low-loss waveguide bends. Appl Phys Lett 84:2022. doi:10.1063/1.1688450
48. Jin N, Rahmat-Samii Y (2005) Parallel particle swarm optimization and finite-difference time-
domain (pso/fdtd) algorithm for multiband and wide-band patch antenna designs. IEEE T
Antenn Propag 53:3459–3468. doi:10.1109/TAP.2005.858842
49. Johnson JM, Rahmat-Samii Y (1999) Genetic algorithms and method of moments (ga/mom)
for the design of integrated antennas. IEEE T Antenn Propag 47:1606–1614. doi:10.1109/8.
805906
50. Karimi-Moghaddam G, Gould R, Bhattacharya S (2012) Numerical investigation of electronic
liquid cooling based on the thermomagnetic effect. In: Proceedings of the ASME 2012 inter-
national mechanical engineering Congress & Exposition, Houston, 9–15 Nov 2012. doi:10.
1115/IMECE2012-87764
186 5 Electromechanical System Simulation and Optimization Studies
74. Nomura T, Sato K, Taguchi K, Kashiwa T, Nishiwaki S (2007) Structural topology optimization
for the design of broadband dielectric resonator antennas using the finite difference time domain
technique. Int J Numer Meth Eng 71:1261–1296. doi:10.1002/nme.1974
75. Nomura T, Ohkado M, Schmalenberg P, Lee J, Ahmed O, Bakr M (2013) Topology optimization
method for microstrips using boundary condition representation and adjoint analysis. European
microwave conference (EuMC), Nuremberg, 6–10 Oct 2013
76. Nomura T, Yoon SW, Lee J, Dede EM (2013) Level set based topology optimization of directly
bonded copper substrates targeting thermal stress minimization on die-substrate bonding line.
10th World Congress on structural and multidisciplinary optimization, Orlando, 19–24 May
2013
77. Olesen LH, Okkels F, Bruus H (2006) A high-level programming-language implementation
of topology optimization applied to steady-state Navier-Stokes flow. Int J Numer Meth Eng
65:975–1001. doi:10.1002/nme.1468
78. Olszewski M (2011) Evaluation of the 2010 Toyota Prius hybrid synergy drive system.
ORNL/TM-2010/253, Oak Ridge National Laboratory, Oak Ridge
79. Özişik MN (1993) Heat conduction, 2nd edn. Wiley, New York
80. Ozoe H (2005) Magnetic convection. Imperial College Press, London
81. Park S, Yoo J, Choi JS (2009) Simultaneous optimal design of the yoke and the coil in the
perpendicular magnetic recording head. IEEE T Magn 45:3668–3671. doi:10.1109/TMAG.
2009.2023878
82. Petosa A, Ittipiboon A, Antar Y, Roscoe D, Cuhaci M (1998) Recent advances in dielectric-
resonator antenna technology. IEEE Antenn Propag M 40:35–48. doi:10.1109/74.706069
83. Reddy JN, Gartling DK (2000) The finite element method in heat transfer and fluid dynamics,
2nd edn. CRC Press, Boca Raton
84. Rosensweig RE (1985) Ferrohydrodynamics. Cambridge University Press, New York
85. Schwab L, Hildebrandt U, Stierstadt K (1983) Magnetic Bénard convection. J Magn Magn
Mater 39:113–114. doi:10.1016/0304-8853(83)90412-2
86. Seo JH (2009) Optimal design of material microstructure for convective heat transfer in a
solid-fluid mixture. Dissertation, University of Michigan, Ann Arbor
87. Sigmund O (2001) A 99 line topology optimization code written in Matlab. Struct Multidiscip
O 21:120–127. doi:10.1007/s001580050176
88. Sobhan CB, Garimella SV (2001) A comparative analysis of studies on heat transfer and fluid
flow in microchannels. Microscale Therm Eng 5:293–311. doi:10.1080/10893950152646759
89. Sullivan PF, Ramadhyani S, Incropera FP (1992) Extended surfaces to enhance impingement
cooling with single circular round jets. Joint ASME/JSME conference on electronic packaging,
Milpitas, 9–12 Apr 1992
90. Svanberg K (1987) The method of moving asymptotes—a new method for structural optimiza-
tion. Int J Numer Meth Eng 24:359–373. doi:10.1002/nme.1620240207
91. Tsuji Y, Hirayama K, Nomura T, Sato K, Nishiwaki S (2006) Design of optical circuit devices
based on topology optimization. IEEE Photonic Tech L 18:850–852. doi:10.1109/LPT.2006.
871686
92. Tuckerman DB, Pease RFW (1981) High-performance heat sinking for VLSI. IEEE Electr
Devices L 2:126–129. doi:10.1109/EDL.1981.25367
93. Yamasaki S, Nishiwaki S, Yamada T, Izui K, Yoshimura M (2010) A structural optimization
method based on the level set method using a new geometry-based re-initialization scheme.
Int J Numer Meth Eng 83:1580–1624. doi:10.1002/nme.2874
94. Yee K (1996) Numerical solution of initial boundary value problems involving Maxwell’s
equations in isotropic media. IEEE T Antenn Propag 14:302–307. doi:10.1109/TAP.1966.
1138693
95. Yoon GH (2010) Topological design of heat dissipating structure with forced convective heat
transfer. J Mech Sci Technol 24:1225–1233. doi:10.1007/s12206-010-0328-1
96. Zienkiewicz OC, Campbell JS (1973) Shape optimization and sequential linear programming.
In: Gallagher RH, Zienkiewicz OC (eds) Optimum structural design. Wiley, New York, pp
109–126
Chapter 6
Extensions to New Topics
A 2-D modeling approach was assumed for many of the multiphysics computa-
tional models and the majority of the optimization studies presented herein, in order
to balance solution accuracy with overall computational cost. The validity of this
assumption is critical to relevance of the numerical results, and it is clear that certain
classes of problems would certainly benefit from the enriched design information that
a 3-D modeling approach can provide. For example, in the case of conjugate heat
transfer, performing an optimization study in two dimensions does little to resolve
the 3-D effects of fin geometry on heat transfer enhancement.
In the case of structural optimization for multiphysics systems, the major obstacle
in adopting a fully 3-D modeling strategy continues to be the excessive computational
cost associated with obtaining design solutions. From a computational mechanics
perspective, the co-design for many physical processes in a greater number of spatial
dimensions equates to more degrees of freedom and larger systems of equations that
must be solved. Thus, while it is expected that increased computing power (e.g., HPC
systems) will assist in making this problem tractable, there exists an even greater need
for better numerical techniques that enable the solution of large-scale multiphysics
optimization problems in three dimensions. Some possible solutions to these large
scale problems may include parallel processing for topology optimization [1, 6], the
use of novel optimization schemes involving polyhedral meshes [8], or isogeometric
analysis for topology optimization [5, 21, 24].
While moving from two to three dimensions in the spatial domain for multi-
physics simulation and optimization opens new research areas, the handling of solu-
tions beyond steady-state analysis must also be considered. The issue of dynamic
response is of growing importance for higher efficiency electromechanical systems,
see for example [13, 16], and the expansion of multiphysics design optimization tools
into the temporal domain for transient analysis or frequency domain for eigenmode
analysis is an area ripe for further exploration.
Additionally, most of the design examples presented in this book have dealt with
linear material properties. However, one example in which this limiting assumption
may be drawn into question is in relation to the evaluation of thermally induced
stresses in electronics solder bond materials. It is well known that, when thermal
cycled due to device power dissipation, solder bond attachment materials found in
electronics can behave in a viscoplastic manner [29], where life cycle prediction is a
function of cumulative damage. Therefore, the incorporation of more advanced non-
linear or temperature-dependent materials and appropriate modeling strategies [4]
into multiphysics design optimization studies is critical to greater understanding of
performance and reliability.
Building off of the discussion of linear versus non-linear materials, we have shown
in Sect. 5.1.5 that the design of anisotropic materials and structures has the potential
to open up new design spaces for certain classes of multiphysics problems. In the
example in Sect. 5.1.5, a second-order tensor (i.e., the anisotropic material thermal
conductivity tensor) was designed to manipulate the flow of heat in a composite
structure. Logical extensions of this work to the design of structural systems involving
higher order tensors is of relevance to both aerospace and automotive industries for
composite design, and the co-design of structural topology and material layout, as
proposed in [12, 17], represents another interesting area for future investigation.
Finally, another topic for greater development in the scaling-up of systems is
the design optimization of multi-body electromechanical assemblies. Several of the
actuator and motor examples described in Chap. 5 represent a starting point, and
supplementing multi-body design approaches with the optimization of additional
physical processes, richer 3-D models, and faster parallel solution techniques has
great potential.
A graphical summary of the various above ideas related to the scaling-up of
multiphysics simulation and design optimization tools and techniques is provided in
the chart in Fig. 6.1.
6.2 Treatment of Surfaces and Interfaces 191
Scaling-Up of Multiphysics
Simulation & Optimization
Fig. 6.1 Concepts related to the scaling-up of multiphysics simulation and design optimization
tools and techniques
Stator
Transition from + +++
Air gap +
laminar to turbulent + +
+ +
Rotor +
+ Reduced
++ + +
Current + current
concentration flow inside
Airfoil near surface conductor
Solid-to-fluid interface Motor air gap Skin effect in high frequency conductors
(flow field transition shown) (zoomed view shown) (conductor cross-section shown)
Fig. 6.2 Example problems for which the treatment of surfaces and interfaces is important in
simulation and optimization
Opposite Specified
taper appearance
feature
Mold Mold
(bottom) (bottom)
Fig. 6.3 A representative example of an undesired negative taper structure (left) versus a constrained
mold fabrication process resulting in a desired external appearance (right). While shown with a ‘box-
like’ external appearance, the outer surface of the structure illustrated on the right may take most
any desired shape, as long it is feasible in terms of manufacturing
approach outlined in Chap. 4, and recently reviewed in [27], has good potential for
enabling advancements in this area. An early study related to boundary-tracking
with an Arbitrary Lagrangian-Eulerian (ALE) mesh for imposing specific bound-
ary conditions in metallic waveguides is found in [28]. This method was also used
with good results in the design optimization and performance validation of metallic
infrared filters, where the transmission spectrum strongly depends on the propaga-
tion of surface plasmon polaritons across the metal-to-dielectric interface [18, 19].
Other numerical approaches, such as the ‘universal mesh’ scheme described in [23]
and the explicit interface representation explained in [3], are being explored. Logical
target applications of this technology include boundary mesh generation for topol-
ogy optimization in fluidics and the development of robust methods for topological
optimization using 3-D conformal meshes.
The scaling-up of systems paired with the treatment of important surfaces and inter-
faces is expected to provide a richer model definition to better capture and predict
the behavior of electromechanical devices. Another feature that should be considered
in conjunction with these enrichment strategies is the implementation of additional
design constraints. Complex electromechanical systems are subject to performance
constraints beyond volume or mass reduction [25], and as described in Sect. 5.3.2,
fabrication constraints are often important to arrive at designs that are feasible from
a manufacturing point of view; see Fig. 6.3.
An approach to this problem is to enforce a minimum length scale constraint
using a projection method, as demonstrated in [9] for 2-D truss structural design. A
further advanced example of this technique for smaller scale devices may be found
in [20], where the authors mimic a MEMS fabrication process (see Fig. 6.4) with
mathematical projection in the design of a process photomask for a structural rib;
sample fabricated devices versus various optimal rib topologies are shown in Fig. 6.5
6.3 Free Versus Constrained Systems-Toward Manufacturability 193
Fig. 6.4 Topology optimization flow using the micromachining process like mapping. Reprinted
from [20, Fig. 3], Copyright (2013), with permission from Elsevier
on the left and right, respectively. Note that a similar strategy was used in Sect. 5.1.2,
where a 2-D etching process photomask pattern was optimized and projected into
3-D for the design of a larger scale DBC substrate that minimized the thermally
induced stresses present at the device bonding layer in an electronics package.
194 6 Extensions to New Topics
Fig. 6.5 SEM images of fabricated top electrodes: a initial dot-array pattern, and b optimized
topologies. c Variation in photomask topology with number of optimization iterations starting from
dot-array pattern. Reprinted from [20, Figs. 6 and 7], 6 pages, Copyright (2013), with permission
from Elsevier
Fig. 6.6 Representative process for synthesizing model geometry—from optimization to rapid
prototype. The figure shows a pin fin heat sink optimized in 3-D for air jet impingement cooling
References
1. Aage N, Lazarov BS (2013) Parallel framework for topology optimization using the method
of moving asymptotes. Struct Multidiscip O 47:493–505. doi:10.1007/s00158-012-0869-2
2. Azernikov S, Fischer A (2004) Efficient surface reconstruction method for distributed CAD.
Comput Aided Des 36:799–808. doi:10.1016/j.cad.2003.09.006
3. Christiansen AN, Nobel-Jørgensen M, Aage N, Sigmund O, Baerentzen JA, (2014) Topol-
ogy optimization using an explicit interface representation. Struct Multidiscip O 49:387–399.
doi:10.1007/s00158-013-0983-9
4. Darveaux R (2000) Effect of simulation methodology on solder joint crack growth correlation.
In: Proceedings of the 50th electronic components and technology conference, Las Vegas,
21–24 May 2000. doi:10.1109/ECTC.2000.853299
5. Dedè L, Borden MJ, Hughes TJR (2012) Isogeometric analysis for topology optimization with
a phase field model. Arch Comput Method E 19:427–465. doi:10.1007/s11831-012-9075-z
6. Evgrafov A, Rupp CJ, Maute K, Dunn ML (2008) Large-scale parallel topology optimiza-
tion using a dual-primal substructuring solver. Struct Multidiscip O 36:329–345. doi:10.1007/
s00158-007-0190-7
7. Fischer A (2011) Engineering-oriented geometry methods for modeling and analyzing scanned
data. J Comput Inf Sci Eng 11:021002. doi:10.1115/1.3593415
8. Gain A (2013) Polytope-based topology optimization using a mimetic-inspired method. Dis-
sertation, University of Illinois at Urbana-Champaign
9. Guest JK, Prévost JH, Belytschko T (2004) Achieving minimum length scale in topology
optimization using nodal design variables and projection functions. Int J Numer Meth Eng
61:238–254. doi:10.1002/nme.1064
10. Haslinger J, Kočvara M, Leugering G, Stingl M (2010) Multidisciplinary free material opti-
mization. SIAM J Appl Math 70:2709–2728. doi:10.1137/090774446
11. Jansen M, Lazarov BS, Schevenels M, Sigmund O (2013) On the similarities between
micro/nano lithography and topology optimization projection methods. In: 10th world congress
on structural and multidisciplinary optimization, Orlando, 19–24 May 2013
12. Lee J, Nomura T, Dede EM (2012) Heat flow control in thermo-magnetic convective systems
using engineered magnetic fields. Appl Phys Lett 101:123507. doi:10.1063/1.4754119
13. Lee J, Yoon SW (2012) Topology design optimization of electromagnetic vibration energy
harvester to maximize power output. J Magn 18:283–288. doi:10.4283/JMAG.2013.18.3.283
14. Lipson H, Kurman M (2013) Fabricated: the new world of 3D printing. Wiley, Indianapolis
15. Malone E, Lipson H (2008) Multi-material freeform fabrication of active systems. In: Proceed-
ings of the ASME 2008 9th biennial conference on engineering systems design and analysis,
Haifa, 7–9 July 2008. doi:10.1115/ESDA2008-59313
16. Noh JY, Yoon GH (2012) Topology optimization of piezoelectric energy harvesting devices
considering static and harmonic dynamic loads. Adv Eng Softw 53:45–60. doi:10.1016/j.
advengsoft.2012.07.008
17. Nomura T, Dede EM, Lee J, Yamasaki S, Matsumori T, Kawamoto A, Kikuchi N (2014)
General topology optimization method with continuous and discrete orientation design using
isoparametric projection (submitted)
18. Ohkado M, Nomura T, Matsumori T, Kawamoto A, Fujikawa H, Sato K, Yamasaki S, Nishi-
waki S (2011) Structural optimization of SPPs color filter using a level set-based topology
optimization method incorporating the ALE method. In: 9th world congress on structural and
multidisciplinary optimization, Shizuoka, 13–17 June 2011
19. Ohkado M, Nomura T, Miura A, Fujikawa H, Ikeda N, Sugimoto Y, Nishiwaki S (2013)
Structural optimization of metallic infrared filters based on extraordinary optical transmission.
T MRS Jap 38:167–170. doi:10.14723/tmrsj.38.167
20. Ozaki T, Nomura T, Fujitsuka N, Shimaoka K, Akashi T (2013) Topology optimization using
multistep mapping from 2D photomask to 3D structure for designing reinforcing rib. Sensors
Actuat A-Phys (in press). doi:10.1016/j.sna.2013.08.033
References 197
21. Qian X (2013) Topology optimization in B-spline space. Comput Method Appl M 265:15–35.
doi:10.1016/j.cma.2013.06.001
22. Qian X, Sigmund O (2013) Topological design of electromechanical actuators with robustness
toward over- and under-etching. Comput Method Appl M 253:237–251. doi:10.1016/j.cma.
2012.08.020
23. Rangarajan R, Lew AJ (2012) Universal meshes: a new paradigm for computing with noncon-
forming triangulations. http://arxiv.org/abs/1201.4903v1. Accessed 19 March 2014
24. Seo YD, Kim HJ, Youn SK (2010) Isogeometric topology optimization using trimmed spline
surfaces. Comput Method Appl M 199:3270–3296. doi:10.1016/j.cma.2010.06.033
25. Sigmund O (2001) Design of multiphysics actuators using topology optimization-Part
I: one-material structures. Comput Method Appl M 190:6577–6604. doi:10.1016/S0045-
7825(01)00251-1
26. Sigmund O (2009) Manufacturing tolerant topology optimization. Acta Mech Sinica 25:227–
239. doi:10.1007/s10409-009-0240-z
27. van Dijk NP, Maute K, Langelaar M, van Keulen F (2013) Level-set methods for structural
topology optimization: a review. Struct Multidiscip O 48:437–472. doi:10.1007/s00158-013-
0912-y
28. Yamasaki S, Nomura T, Kawamoto A, Sato K, Nishiwaki S (2011) A level set-based topology
optimization method targeting metallic waveguide design problems. Int J Numer Meth Eng
87:844–868. doi:10.1002/nme.3135
29. Ye H, Lin M, Basaran C (2002) Failure modes and FEM analysis of power electronic packaging.
Finite Elem Anal Des 38:601–612. doi:10.1016/S0168-874X(01)00094-4
Chapter 7
Appendix: Sample Multiphysics
Optimization Code
In this appendix, example numerical code is provided for a 2-D electrothermal design
optimization problem similar to the example presented in Sect. 5.1.1. Here, Joule
heating considering side/surface convection (as described in [3]) of the conductor is
once more assumed. However, in this case, the solid/void surface convection bound-
ary condition is interpolated using the first derivative of the logistic function, for
simplicity, instead of the smoothed hat-function, Eq. (5.7).
The finite element implementation of the problem follows the description in [1]
for the finite element analysis (FEA) of a thermoelectric device, with the excep-
tion that Seebeck and Peltier effects are neglected. The code is written in a custom
MATLAB®1 script that exploits a sequential solver for the FEA portion of the prob-
lem, where the electrical analysis is first performed and the generated resistive heating
thermal loads are passed to the following heat transfer analysis. The code incorpo-
rates the OC optimizer and filtering routine provided in [4] for minimum length scale
control. The code is constructed as a continuous script instead of a series of functions
for clarity, and numerous comments have been added throughout the script to explain
the various steps implemented in the program.
A 2-D schematic of the problem set up is provided in Fig. 7.1, where a square
design domain is specified with a current source assumed at Terminal 1 (top left)
plus electrical ground and a fixed temperature condition at Terminal 2 (bottom right).
Material properties of copper were assumed for both the material electrical conduc-
tivity and thermal conductivity. Further interpretation of the sample code is left to the
reader, and additional FEA implementation details may be found in [2]. The reader
is also referred to Sect. 5.1.1 for a slightly modified, yet related, example problem
that was implemented instead in a commercial software environment.
Terminal 1
(Current Source)
Design Dependent
Convective Heat
Transfer, h( )
Terminal 2
(Electrical Ground,
y
Fixed Temperature)
x All Remaining
Boundaries Adiabatic
Fig. 7.1 Assumed 2-D simulation domain for electrothermal optimization example written in
MATLAB® programming language. A current source is assumed at Terminal 1 with electrical
ground and thermal ground at Terminal 2. Design dependent convective heat transfer, h(γ ), is
assumed inside the domain with all remaining external boundaries considered adiabatic
1 %---------------------------------------------------------------
2 % Multiphysics Topology Optimization Code
3 % Written by: E.M. Dede, T. Nomura, J. Lee
4 % Last Modified: April 1st, 2014
5 % Description: This code optimizes a 2-D electrothermal system
6 % with current input to the top left terminal of the square
7 % domain and electrical/thermal ground at the bottom right
8 % terminal. The electrical power is used as the internal heat
9 % generation for the thermal portion. The electrical
10 % conductivity, thermal conductivity, and surface convection
11 % coefficient are interpolated on the domain. A sequential FEA
12 % solution process is assumed. An OC optimizer is used
13 %---------------------------------------------------------------
14 clc % Clear the workspace
15 clear all % Clear the variables
16
17 maxiter = 100; % Max # of iterations
18
19 nelx = 50.0; % # of elements in x-dir
20 nely = 50.0; % # of elements in y-dir
21 volfrac = 0.2; % Define volume fraction
7.1 MATLAB® Example Program for Multiphysics Topology 201
66 Ue = sparse((nely+1)*(nelx+1),1);
67
68 % Loop to form the global stiffness matrix
69 for ely = 1:nely % Loop over # el y-dir.
70 for elx = 1:nelx % Loop over # el x-dir.
71 % Upper left el node number in GCS
72 n1 = (nely+1)*(elx-1)+ely;
73 % Upper right el node number in GCS
74 n2 = (nely+1)* elx +ely;
75 % Element dof
76 edof = [n1; n2; n2+1; n1+1];
77 % Populate electrical system K matrix
78 Ke(edof,edof) = Ke(edof,edof) + ...
79 (0.001+0.999*x(ely,elx)ˆpenal)*KEe;
80 end
81 end
82
109 Qe=sigma*(px.ˆ2+py.ˆ2);
110
111 % Thermal FEA
112 %-----------------------------------------------------------
113 k = 400; % Thermal conductivity
114 KEt = k*[ 2/3 -1/6 -1/3 -1/6% Conductivity matrix
115 -1/6 2/3 -1/6 -1/3
116 -1/3 -1/6 2/3 -1/6
117 -1/6 -1/3 -1/6 2/3];
118
119 KEh = ho*[ 2/3 1/6 0 1/6 % Convection matrix
120 1/6 2/3 1/6 0
121 0 1/6 2/3 1/6
122 1/6 0 1/6 2/3];
123 FEh = ho*Tamb; % Convective loads
124
125 % Initialize thermal system K matrix plus global thermal
126 % load, convective load, and temperature displacement
127 % vectors
128 Kt = sparse((nelx+1)*(nely+1), (nelx+1)*(nely+1));
129 Ft = sparse((nely+1)*(nelx+1),1);
130 Fh = sparse((nely+1)*(nelx+1),1);
131 Ut = sparse((nely+1)*(nelx+1),1);
132
133 % Loop to form the global stiffness matrix
134 for ely = 1:nely % Loop over # el y-dir.
135 for elx = 1:nelx % Loop over # el x-dir.
136 % Upper left el node number in GCS
137 n1 = (nely+1)*(elx-1)+ely;
138 % Upper right el node number in GCS
139 n2 = (nely+1)* elx +ely;
140 % Element dof
141 edof = [n1; n2; n2+1; n1+1];
142 % Logistic function for ’h’ interpolation
143 lfi = (1/(1+exp(-lfp*(x(ely,elx)-lfoff))));
144 % Populate thermal K matrix
145 Kt(edof,edof) = Kt(edof,edof) + ...
146 (0.001+0.999*x(ely,elx)ˆpenal)*KEt + ...
147 (lfi*(1-lfi))*KEh;
148 end
149 end
150
208
209 % Filtering of sensitivities
210 %-----------------------------------------------------------
211 % Initialize modified sensitivity matrix
212 dcn=zeros(nely,nelx);
213 for i=1:nelx % Loop # el x-dir.
214 for j=1:nely % Loop # el y-dir.
215 sum=0.0;
216 for k = max(i-round(rmin),1):...
217 min(i+round(rmin),nelx)
218 for l = max(j-round(rmin),1):...
219 min(j+round(rmin),nely)
220 % Compute conv. factor
221 fac = rmin-sqrt((i-k)ˆ2+(j-l)ˆ2);
222 sum = sum+max(0,fac);
223 % Modified sensitivity
224 dcn(j,i) = dcn(j,i) + ...
225 max(0,fac)*x(l,k)*dc(l,k);
226 end
227 end
228 % Finish compute modifying sensitivities
229 dcn(j,i)=dcn(j,i)/(x(j,i)*sum);
230 end
231 end
232 dc=dcn;
233
234 % Design update by the optimality criteria method
235 %-----------------------------------------------------------
236 % Bounds for lagrange multi. & move-limit
237 l1=0; l2=100000; move=0.2;
238 % Lagrange multiplier convergence criteria
206 7 Appendix: Sample Multiphysics Optimization Code
259
260 % Print topology results
261 %-----------------------------------------------------------
262 change = max(max(abs(x-xold)));
263 disp([’ It.: ’ sprintf(’%4i’,loop) ...
264 ’ Obj.: ’ sprintf(’%10.4f’,c) ...
265 ’ Vol.: ’ sprintf(’%6.3f’,xsum/(nelx*nely))...
266 ’ ch.: ’ sprintf(’%6.3f’,change) ])
267
268 %Plot the density
269 colormap(gray); imagesc(-x); axis equal; ...
270 axis tight; axis off; pause(1e-8);
271 title(’Density [Dark Regions = Solid; Light Regions = Void]’)
272
273 end % Loop for iterations
274 toc
275
276 % Post-processing
277 %---------------------------------------------------------------
278 % Plot the electric potential
279 Ueplot = zeros(nely+1,nelx+1);
280 for i = 1:nelx+1
281 for j = 1:nely+1
7.1 MATLAB® Example Program for Multiphysics Topology 207
287
288 % Plot the E-field
289 figure; quiver(xx,nely-yy,-px,py,3); ...
290 axis equal; axis tight; axis off
291 title(’Electric Field Vector Plot’)
292
293 % Plot the electrical power
294 figure; imagesc(Qe); axis equal; axis tight; axis off; colorbar
295 title(’Power [SI Units: W]’)
296
297 % Plot the temperature field
298 Utplot = zeros(nely+1,nelx+1);
299 for i = 1:nelx+1
300 for j = 1:nely+1
301 Utplot(j,i) = Ut((i-1)*(nely+1)+j);
302 end
303 end
304 figure; imagesc(Utplot); axis equal; axis tight; axis off; colorbar
305 title(’Temperature [SI Units: K]’)
306
307 % Plot the ’h’ field
308 Fhplot = zeros(nely+1,nelx+1);
309 for i = 1:nelx+1
310 for j = 1:nely+1
311 Fhplot(j,i) = Fh((i-1)*(nely+1)+j)/Tamb;
312 end
313 end
314 figure; imagesc(Fhplot); axis equal; axis tight; axis off; colorbar
315 title(’Surface Convection Coefficient [SI Units: W/(mˆ2*K)]’)
208 7 Appendix: Sample Multiphysics Optimization Code
References
1. Antonova EE, Looman DC (2005) Finite elements for thermoelectric device analysis in
ANSYS. Paper presented at the 24th international conference on thermoelectrics, Clemson
University, Clemson, 19–23. doi:10.1109/ICT.2005.1519922
2. Kwon YW, Bang H (2000) The finite element method using MATLAB, 2nd edn. CRC Press,
Boca Raton
3. Seo JH (2009) Optimal design of material microstructure for convective heat transfer in a solid-
fluid mixture. Dissertation, University of Michigan, Ann Arbor
4. Sigmund O (2001) A 99 line topology optimization code written in Matlab. Struct Multidiscip
O 21:120–127. doi:10.1007/s001580050176
Index
A transition, 126
Actuator, 2, 4, 7, 11, 12, 17–19, 26, 42, 61, velocity, 30, 50, 80
62, 150, 151, 155, 157–168, 190 Boussinesq approximation, 89
Additive manufacturing, 195 Bus bar, 62–64, 66, 67, 112, 113
Adjoint method, 43–46, 124, 125, 128, 131,
138, 139, 144, 151, 154, 155, 158,
160 C
Altair, 8 Co-energy, 174
Ampère’s law, 30, 34, 35 Coefficient of thermal expansion (CTE), 9,
Anisotropic materials, 61, 62, 100, 102–107, 13, 28, 69, 70
113, 190 Cold plate, 13, 57, 61, 62, 84–90
ANSYS, 5, 8 Commercial software, 3, 5, 6, 8, 23, 43, 45,
Arbitrary Lagrangian-Eulerian (ALE) mesh, 46, 66, 73, 80, 83, 85, 86, 93, 180,
192 182
Compliance, 76, 157
mechanical, 151–153, 157
B minimization, 7, 55, 141, 142, 145, 153,
Boundary conditions, 5, 21, 24–29, 32, 34, 154, 156, 157
46, 50, 55, 81, 125, 126, 161, 165 thermal, 106
absorbing, 137 Compliant mechanism, 7, 42
adiabatic, 46 Composites, 27, 42, 57, 62, 70, 100, 102–
convection, 26, 64, 65, 96, 119, 184 105, 107–114, 147, 158, 190
current density, 26, 33, 65 Computer numerically controlled (CNC),
displacement, 23, 55 195
electric field potential, 26, 31 Computer-aided design (CAD), 6, 23, 84, 85,
heat flux, 26, 28, 30, 65, 73, 80, 86 130–132, 134, 194
impedance, 124, 126 Computer-aided engineering (CAE), 2, 6
magnetic field potential, 31 COMSOL, 5, 8, 43, 45, 180
perfect electric conductor (PEC), 37, Conservation of energy, 25, 174
124, 135 Conservation of linear momentum, 29
radiation, 26 Conservation of mass, 28
Robin, 37 Converter, 18, 68
Sommerfeld radiation, 38 Coulomb forces, 17
surface stress, 30, 50, 80, 86 Coulomb Virtual Work (CVW), 151
symmetry, 85 Coupled-field analysis
temperature, 26, 28, 30, 46, 80, 91, 107, direct, 5
119 segregated solver, 5
F
Fabrication constraints, 62, 134, 145, 149, I
150, 192 Inductor, 11, 14, 15, 31, 62, 114, 116–121,
Faraday’s effect, 17 123
Faraday’s law, 30, 34, 35 Infrared, 15, 16, 113, 124
Filter, 47–49, 66, 68, 104, 105, 108, 128, 199 Inverter, 4, 18, 68
Finite difference-time domain (FDTD), 34, Isight, 5
35, 124, 134, 139, 143, 147 Isogeometric analysis, 190
Finite element analysis (FEA), 1–3, 5, 21, Isoparametric element, 37
23, 24, 56, 68, 182, 199 Isotropic materials, 7, 21, 22, 27, 41, 42, 100,
post-processing, 6, 24, 43 103, 106, 107, 125
Index 211
J O
Jet impingement, 84, 85, 88, 195 Ohm’s law, 25
Joule heating, 13–15, 18, 21, 25, 26, 30, 63, Optical tweezers, 17
116, 199 Optimality criteria (OC), 68, 199
Optimization
multi-objective, 8, 42, 64, 79, 81, 82, 84,
K 104, 106, 107, 153, 199
Kerr effect, 16 shape, 1, 6, 8, 9, 41, 126
size, 1, 6, 8, 41, 70, 84–87
topology, 1, 4–8, 41–44, 46–49, 51, 52,
L 63, 68, 76, 90, 100, 120, 125, 133, 145,
Lagrange multiplier, 55 153, 158, 169, 190, 191, 195, 200
Laminar flow, 28, 49, 79, 80, 82, 89
Level set method, 41, 42, 52–56, 69, 74, 75,
147, 191 P
Linear-elastic materials, 21, 22, 27 Parallel processing, 190
Logistic function, 199 Pareto front, 8, 82, 83, 143–145, 149
Peltier effect, 199
Perimeter control, 54
M
Permanent magnet (PM), 17, 91–94, 96–
Magnetic components, 11, 14, 15, 18, 31, 61,
100, 158–161, 163, 165–168,
62, 114–117
180, 182
Magnetic energy, 14, 18, 114, 115, 119, 170
Permittivity, 16, 17, 31, 35, 124, 133, 135,
Magnetic flux linkage, 115
136, 147, 181
Magnetic permeability, 14, 18, 31, 35, 92,
Phononic materials, 7, 42
116, 122, 154, 170
Photolithography, 194
Magnetostatic analysis, 33, 34, 151–153,
Photomask, 192
158, 182
Plane-strain, 24
Material non-linearity, 17, 27, 174, 190
Plane-stress, 55
MATLAB, 8, 46, 68, 199, 200
Pockel’s effect, 16
Maxwell Stress Tensor (MST) method, 151,
Point cloud data, 194
152, 158, 182
Poisson’s ratio, 24, 55
Maxwell’s equations, 30, 31, 34, 92, 120,
Power electronics, 69
133, 137, 138
Powertrain, see drive-train
Mechanical stress, 1, 17, 22, 56
Poynting vector, 127, 129
Mesh refinement, 24, 46, 48, 81, 191
Printed circuit board (PCB), 13, 100, 112–
Metamaterial, 100, 158
114, 125, 126, 129, 130
Method of moving asymptotes (MMA), 43,
Proximity effects, 116
46, 66, 94, 104
Microelectromechanical systems (MEMS),
11, 17, 26
Microwave heating, 16 R
Motor, 2, 4, 11, 12, 17–19, 61, 62, 68, 169, Radiation efficiency, 136
170, 178–180, 182, 184, 190, 191 Rayleigh-Bénard convection, 96
Multi-axis machining, 195 Reactive impedance, 15
Multi-body systems, 190 Resistive heating, see Joule heating
Reynolds number, 49, 89
RF components, 11, 12, 15, 124
N Robotics, 17
Nédélec conformal element, 36 Robust design, 48, 195
NASTRAN, 1
National Aeronautics and Space Administra-
tion (NASA), 1 S
Neumann conditions, 26, 30, 161, 165 Seebeck effect, 199
212 Index