OPTIMIZATION OF ELECTROMAGNETIC DEVICES AND MATERIALS
By
Kazuko Fuchi

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Mechanical Engineering
2012

ABSTRACT
OPTIMIZATION OF ELECTROMAGNETIC DEVICES AND MATERIALS
By
Kazuko Fuchi
Topology optimization is a computational design methodology that is used to find topologies of optimal design, that match design goals. Application of topology optimization in the
area of electromagnetics is extremely interesting, since for many electromagnetic design problems the topology of inclusion structures in materials or components in devices can have a
dramatic effect on the interaction of the designed structure with the electromagnetic waves.
The success of the use of topology optimization for design of metamaterials, antennas, waveguides, etc. reported in literature suggests the importance of research into effective uses of
topology optimization in electromagnetic design problems.
The purpose of this dissertation is to investigate efficient electromagnetic analysis methods that can be combined with the method of topology optimization, as well as to develop a
new way of using topology optimization to design electromagnetic devices with transformable
complex surface geometry, using the concept of origami.
The first part of the dissertation focuses on the investigation of a topology optimization
framework for design of periodic structures for electromagnetic applications, using a rigorous
and efficient finite element analysis method. The second part introduces a topology optimization based design method for origami folding patterns. To illustrate its use, the origami
design method is employed to design frequency selective surfaces for electromagnetic applications, that can transform in their surface geometry to alter their working frequency through
folding and unfolding.

To my parents, Hiroyuki Fuchi and Kyoko Makino

iii

ACKNOWLEDGMENTS

I have been supported and encouraged by many people throughout my time at Michigan
State University. I would like to acknowledge those who have directly helped me through
my Ph.D. program.
Foremost, I would like to express my deepest gratitude to my academic advisor, Prof.
Alejandro R. Diaz. He has given me all that I could asked for, and far beyond. He has taken
me around the world for professional meetings; he has taught me everything about research
and being an academic that I know now; and he has been an exceptional friend.
The best thing about my research project was its multidisciplinary nature. I had the
honor to work with Prof. Edward Rothwell from electrical engineering, who has kept me
ecstatic with his numerous insightful and creative suggestions. I am very appreciative of the
members of the Metamaterials Group, Dr. Raoul Ouedraogo and Junyan Tang for the shared
enthusiasm towards our collaborative work, and Khoa Nguyen and Benjamin Crowgey for
their help on parts of the project. Extended thanks to my fellow students, John Tam, Isil
Berkun, Xiankai Song, Tingli Cai, and Smruti Panigrahi for their continuous encouragement
and invaluable friendship. Special thanks to Prof. Alan Haddow, who was my mentor
for my first TA assignment, Prof. Jongeun Choi and Prof. Seungik Baek, who agreed to
serve on my Ph.D. committee, and the graduate coordinator, Prof. Andre Benard, for their
generous support and helpful discussions. My sincere gratitude goes to my collaborators,
Prof. Shinji Nishiwaki and Prof. Takayuki Yamada, who have introduced me to a new
topology optimization method and stimulated me.
I am grateful to Prof. Martin Berz for offering me a tremendous amount of help getting
adopted to the life in the US, my mother, Prof. Kyoko Makino, for being an extraordinary
iv

role model and an inspiration, and my father, Dr. Hiroyuki Fuchi, for always understanding
me and motivating me towards a meaningful life. Finally, my boyfriend, Eric Wolf deserves
my loving thanks for being on my side and leading me to becoming a happier and better
person, and the Wolf family for making me feel at home.

v

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . .

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

ix

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

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

xi

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

Finite Element Formulation for Time Harmonic Electromagnetic Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Constitutive Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Analysis of Periodic Media . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Two-dimensional domain, 1-D periodicity . . . . . . . . . . . . . . . .
2.2.2 Three-dimensional domain, 2-D periodicity . . . . . . . . . . . . . . .
Finite Element Mesh Truncation Techniques . . . . . . . . . . . . . . . . . .
2.3.1 Absorbing boundary condition . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Perfectly matched layer . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.3 Plane wave expansion . . . . . . . . . . . . . . . . . . . . . . . . . . .
Finite Element Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
9
9
10
11
12
13
14
17

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

18
19
20
20
21
24
25
26
29
29
29
30
31
32
33
34

List of Figures
Chapter 1
Chapter 2
2.1
2.2

2.3

2.4

Chapter 3 Topology Optimization
3.1 Problem Setup . . . . . . . . .
3.2 Level Set Approach . . . . . . .
3.2.1 Material distribution . .
3.2.2 Optimization problem .
3.2.3 Numerical examples . .
3.2.3.1 Example 1 . . .
3.2.3.2 Example 2 . . .
3.2.4 Conclusion . . . . . . . .
3.3 Density Approach . . . . . . . .
3.3.1 Material distribution . .
3.3.2 Optimization problem .
3.3.3 Optimality condition . .
3.3.4 Sensitivity analysis . . .
3.3.5 Influence of parameter p
3.3.6 Numerical examples . .

of 2D
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
vi

Periodic Dielectric Structures
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

35
38
40
44
48

Chapter 4 Topology Optimization of 3D Periodic Structures
4.1 Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Material distribution . . . . . . . . . . . . . . . . . . . . . . .
4.3 Optimization problem . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Sensitivity analysis . . . . . . . . . . . . . . . . . . . .
4.3.2 Projection of the effective density . . . . . . . . . . . .
4.3.3 Numerical example . . . . . . . . . . . . . . . . . . . .
4.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

50
50
52
53
54
56
56
61

Applications
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .

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

62
63
63
65
67
67
69
70
71
71
74
77
77
78
80
84
92
94

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

98
100
102
104
105
106
106
107
110
112

3.3.7
3.3.8

3.3.6.1 Example 1 . . . . .
3.3.6.2 Example 2 . . . . .
3.3.6.3 Example 3 . . . . .
Designs with parallel bands
Conclusion . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

Chapter 5 Origami Tunable Surfaces for Electromagnetic
5.1 Origami Tunable Frequency Selective Surfaces . . . . . . .
5.1.1 Miura-ori . . . . . . . . . . . . . . . . . . . . . . .
5.1.2 Performance evaluation . . . . . . . . . . . . . . . .
5.1.3 Choice of conducting element types . . . . . . . . .
5.1.3.1 Concentric double ring copper prints . . .
5.1.3.2 Cross-shaped copper prints . . . . . . . .
5.1.3.3 Summary . . . . . . . . . . . . . . . . . .
5.1.4 Angle dependency . . . . . . . . . . . . . . . . . . .
5.1.5 Polarization dependency . . . . . . . . . . . . . . .
5.1.6 Experiment . . . . . . . . . . . . . . . . . . . . . .
5.1.7 Conclusions . . . . . . . . . . . . . . . . . . . . . .
5.2 Origami Tunable Metamaterials . . . . . . . . . . . . . . .
5.2.1 Design of an origami tunable metamaterial . . . . .
5.2.2 Demonstration of tuning of a metamaterial . . . . .
5.2.3 Principle of operation . . . . . . . . . . . . . . . . .
5.2.4 Experiment . . . . . . . . . . . . . . . . . . . . . .
5.2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . .
Chapter 6 Topology Optimization for Origami
6.1 Construction of the Ground Structure . . . . .
6.2 Crease Type Assignment and a Folded State .
6.3 Foldability Conditions . . . . . . . . . . . . .
6.4 Optimization of Origami Design . . . . . . . .
6.4.1 Geometric properties . . . . . . . . . .
6.4.2 Objective function . . . . . . . . . . .
6.4.3 Optimization problem . . . . . . . . .
6.4.4 Sensitivity analysis . . . . . . . . . . .
6.5 Numerical examples . . . . . . . . . . . . . . .
vii

Design
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.

6.6
6.7

6.5.1 Example 1 .
6.5.2 Example 2 .
Application: design
Conclusions . . . .

. . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . .
of polarization insensitive tunable
. . . . . . . . . . . . . . . . . . .

. . .
. . .
FSS
. . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

113
115
117
120

Chapter 7

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122

Chapter A

Finite Element Equations for 2D Periodic Structures . . . . . . 126

Chapter B

Finite Element Equations for 3D Periodic Structures . . . . . . 132

Bibliography . . . . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . . . 141

LIST OF TABLES

Table 3.1

Parameters for example 1 . . . . . . . . . . . . . . . . . . . . . . . .

25

Table 3.2

Parameters for example 2 . . . . . . . . . . . . . . . . . . . . . . . .

26

Table 3.3

Parameters used in Example 1 . . . . . . . . . . . . . . . . . . . . .

35

Table 3.4

Parameters used in Example 2 . . . . . . . . . . . . . . . . . . . . .

38

Table 3.5

Parameters used in Example 3 . . . . . . . . . . . . . . . . . . . . .

41

Table 3.6

Parameters used for the angle dependency

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

47

Table 4.1

Parameters used in the example . . . . . . . . . . . . . . . . . . . .

58

Table 5.1

Parameters used in sample F-band design . . . . . . . . . . . . . . .

81

Table 5.2

Resonance frequencies found using HFSS . . . . . . . . . . . . . . .

83

ix

LIST OF FIGURES

Figure 1.1

Design of an SRR . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

Figure 2.1

A two-dimensional representative cell . . . . . . . . . . . . . . . . .

10

Figure 2.2

A three-dimensional representative cell . . . . . . . . . . . . . . . .

11

Figure 2.3

PML setup for 1-D and 2-D periodic structures . . . . . . . . . . . .

13

Figure 3.1

Design Domain

19

Figure 3.2

Optimal structure and the field distribution for example 1. For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this thesis (or dissertation). 26

Figure 3.3

Optimal structure and the field distribution for example 2. . . . . .

27

Figure 3.4

Frequency sweep for the optimum structure for example 2. . . . . .

27

Figure 3.5

A known band-gap structure and the frequency sweep. . . . . . . . .

28

Figure 3.6

Example 1. Optimum solution and iteration history . . . . . . . . .

36

Figure 3.7

Example 1. Optimum solution with vacuum background material . .

37

Figure 3.8

Example 2. Solution and its wavelength sweep . . . . . . . . . . . .

39

Figure 3.9

Example 3. Initial design and optimum solution after postprocessing

41

Figure 3.10

Example 3. Wavelength and angle sweeps (Nm = 10) . . . . . . . .

42

Figure 3.11

Example 3. Angle sweep computed using one mode only (Nm = 0) .

43

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

x

Figure 3.12

Example 3. Optimum structure starting from a homogeneous material distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Figure 3.13

Wave propagation through vertical and inclined bands . . . . . . . .

45

Figure 3.14

Optimization using slanted bands for the initial guess . . . . . . . .

46

Figure 3.15

Angle dependency of the path length dL . . . . . . . . . . . . . . .

47

Figure 4.1

Design Domain

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

51

Figure 4.2

Projection function . . . . . . . . . . . . . . . . . . . . . . . . . . .

57

Figure 4.3

Optimal solution and iteration history . . . . . . . . . . . . . . . . .

59

Figure 4.4

Processed solutions . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 5.1

Miura-ori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

Figure 5.2

Unit cell of Miura-ori . . . . . . . . . . . . . . . . . . . . . . . . . .

64

Figure 5.3

Analysis setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Figure 5.4

Concentric ring conducting element . . . . . . . . . . . . . . . . . .

67

Figure 5.5

Simulated transmission coefficient |S21 | for concentric rings at normal
incidence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

Figure 5.6

Skewed cross-shaped conducting element . . . . . . . . . . . . . . .

69

Figure 5.7

Simulated (Sim.) and measured (Exp.) transmission coefficient |S21 |
at normal incidence . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

Simulated and measured transmission coefficient |S21 | at oblique incidence at θ=30° . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

Figure 5.9

Simulated and measured resonance frequencies, fres in GHz . . . .

73

Figure 5.10

Fabricated origami tunable FSS . . . . . . . . . . . . . . . . . . . .

74

Figure 5.8

xi

Figure 5.11

Experiment setup for normal incidence . . . . . . . . . . . . . . . .

75

Figure 5.12

Experiment setup for oblique incidence at θ = 30° . . . . . . . . . .

76

Figure 5.13

Corrugated sheet with SRRs, α = 0° . . . . . . . . . . . . . . . . . .

78

Figure 5.14

Corrugated sheet with SRRs, α = 10° . . . . . . . . . . . . . . . . .

79

Figure 5.15

Origami metamaterial unit cell design . . . . . . . . . . . . . . . . .

80

Figure 5.16

Transmission coefficient, |S21 | . . . . . . . . . . . . . . . . . . . . .

81

Figure 5.17

Reflection coefficient, |S11 | . . . . . . . . . . . . . . . . . . . . . . .

82

Figure 5.18

Equivalent circuit of an SRR . . . . . . . . . . . . . . . . . . . . . .

85

Figure 5.19

Electric field magnitude in one unit cell for α = 0° . . . . . . . . . .

86

Figure 5.20

Electric field magnitude in one unit cell for α = 2° . . . . . . . . . .

86

Figure 5.21

Electric field magnitude in one unit cell for α = 8° . . . . . . . . . .

87

Figure 5.22

Electric field magnitude in one unit cell for α = 10° . . . . . . . . .

87

Figure 5.23

Geometry of offset parallel conducting strips . . . . . . . . . . . . .

89

Figure 5.24

Geometry of unfolded origami metamaterial unit cell . . . . . . . . .

89

Figure 5.25

Relative change in resonance frequency (dashed curves are best-fit
lines) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90

Figure 5.26

Derivative of relative change in resonance frequency . . . . . . . . .

91

Figure 5.27

Two halves of fabricated origami metamaterial unit cell with the unit
cell width w = 6.85 mm . . . . . . . . . . . . . . . . . . . . . . . . .

93

Origami metamaterial unit cell placed into an F-band waveguide sample holder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

Figure 5.28

xii

Figure 5.29

Figure 5.30

Figure 5.31

Measured transmission coefficients for an origami metamaterial unit
cell placed into an F-band waveguide . . . . . . . . . . . . . . . . .

95

Resonance frequency of an origami metamaterial unit cell placed into
an F-band waveguide . . . . . . . . . . . . . . . . . . . . . . . . . .

96

Transmission coefficients for an origami metamaterial unit cell placed
into an F-band waveguide from HFSS simulations . . . . . . . . . .

97

Figure 6.1

Origami grid systems . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Figure 6.2

Variations of ground structures . . . . . . . . . . . . . . . . . . . . . 101

Figure 6.3

Flat and folded states of a single-vertex crease . . . . . . . . . . . . 103

Figure 6.4

Single-vertex crease . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 6.5

Objective function f . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Figure 6.6

Flowchart of foldline elimination algorithm . . . . . . . . . . . . . . 111

Figure 6.7

Labeled vertices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

Figure 6.8

Candidates of fixed folding angles . . . . . . . . . . . . . . . . . . . 114

Figure 6.9

Design 1, Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 6.10

Design 2, Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure 6.11

Design 1, Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure 6.12

Design 2, example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

Figure 6.13

Conducting decoration for FSS . . . . . . . . . . . . . . . . . . . . . 118

Figure 6.14

Transmission coefficient . . . . . . . . . . . . . . . . . . . . . . . . . 119

Figure 6.15

Transmission coefficient for the two-layer design . . . . . . . . . . . 120

Figure B.1

Brick element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
xiii

Chapter 1
Introduction
Topology optimization is a computational design methodology that is used to find topologies
of optimal designs that match design goals. It was first introduced in the context of structural
design by Bendse and Kikuchi in [1]. Optimal topologies of structures are found through
indicator functions, used to represent the effective “density” of a material in each element
in the finite element formulation. In [1] effective material properties of each element are
computed using a homogenization method and used in the constitutive equations. Within
the last two decades, the use of topology optimization has been extended to material design
[2] with target properties involving many different physics [3], including thermal [4] and fluid
mechanics [5] [6] [7], acoustics [8] [9], photonics and electromagnetics [10]-[26].
Application of topology optimization in the area of electromagnetics is extremely interesting, since for many electromagnetic design problems the topology can have a dramatic effect
on the interaction of the designed structure with the electromagnetic waves. An example of
such instance is design of metamaterials. Metamaterials are engineered materials composed
of metallic or dielectric inclusions within a unit cell of a periodic array, and they exhibit
1

Figure 1.1: Design of an SRR

extraordinary electromagnetic properties. A particular type of metamaterial is known to exhibit negative index of refraction, a property that is not possessed by any natural occurring
material. The negative index of refraction material is theorized by Pendry in [27], and a composite design in the radio frequency (RF) range was proposed by Smith [28] and fabricated
and tested in [29] where the possibility of synthesizing such materials was demonstrated.
In [28] and [29], a combination of metallic elements; a wire for negative electric permittivity and a split ring resonator (SRR) for negative magnetic permeability, is used to achieve a
negative index of refraction. One example of SRRs is the concentric rings with gaps, shown
in Fig.1.1. The topology of such an element is of crucial importance, and a small variation
to its topology could alter the design’s electromagnetic behavior dramatically. For instance,
filling the gaps of the SRRs and closing the rings would create a short, and the negative
magnetic permeability will no longer be observed. The ability to find the topology of one
design with such a unique character is of significant value, as the design can be modified
through e.g., parametric or shape optimization, to be incorporated into various applications.
In case of the negative index refraction metamaterial in [28] and [29], the topology of the
original design was found through physical insights. Topology optimization can be an effi2

cient alternative to finding the starting point of a design process. The work done by Diaz and
Sigmund in [19] revealed that metamaterials with negative permeability can be designed via
topology optimization, and not all the negative permeability elements have to have a similar
geometry as the SRRs.
A popular choice of a design optimization algorithm for electromagnetic problems is a
genetic algorithm (GA). Its main advantage is simplicity in implementation, especially for
design of advanced electromagnetic devices that requires complex analysis. An introduction
to GA for electromagnetic problems and details on the implementation of sample problems
are provided in [30]. In [31] Ouedraogo et al. use a GA-based topology optimization method
to design a miniaturized antenna. A rectangular section placed above a loop antenna is
treated as a design domain, which is divided into rectangular subdomains. Each subdomain is assigned a 1 or 0, which indicates whether copper should be placed(1) or not(0). A
combination of 1’s and 0’s constitutes a chromosome, and a best fit chromosome is selected
by means of natural selection through mating, mutation and selection of best performing
chromosomes at each iteration. Similar approaches have been used to design frequency selective surfaces, metamaterials and antennas in [32]-[34]. A notable disadvantage of designing
electromagnetic devices using GA-based topology optimization is that many iterations are
necessary to obtain an acceptable design. Moreover, the stopping criteria are not based on
optimality conditions, and the solution is not guaranteed to be optimal.
Rigorous topology optimization methods using gradient-based algorithms have also been
used to design electromagnetic devices in [10]-[26]. Gradient-based topology optimization
methods can lead to faster convergence, and solutions satisfy optimality conditions. The
challenge in using gradient-based methods is the computation of sensitivity of the perfor3

mance with respect to the design. Special care needs to be taken in choosing the design
control and performance measure, and analysis methods that allow for efficient computation
of sensitivity need to be selected. Finite element methods (FEM) and finite difference time
domain (FDTD) methods are the two common analysis methods that have been used for
analysis in the frequency domain and time domain, respectively, when using topology optimization as the design method. Earlier works employ the finite element method for the
availability of an efficient sensitivity analysis method based on an adjoint variable problem
formulation. A typical choice of the design variable is an indicator function used to describe
the constituent dielectric or conducting material distributed within the design domain. Examples of the use of topology optimization in conjunction with the finite element analysis
include the design of patch antenna substrates by Kiziltas et al. [11], photonic waveguides
by Jensen and Sigmund [12], [14], and periodic structures for transmission power control by
Nomura et al. [20]. In those works the density function describing the dielectric material
distribution is defined over two-dimensional domains and used to control the design through
an interpolation scheme. Aage et al. [25] and Erentok and Sigmund [26] used a similar formulation but with the density function defined over a three dimensional domain describing
the distribution of conducting material for design of antennas. Development of the adjoint
variable method for the sensitivity analysis for FDTD in [35]-[37] let us to choose FDTD as
an analysis method when appropriate. In [16] Nomura et al. designed dielectric antennas
via topology optimization using a density function as the design variable that describes the
dielectric distribution within a 3D design domain, and the FDTD method to find the electric
and magnetic fields, which are used in the performance measure.

Increasing demand for an improved performance in advanced electromagnetic devices
4

and materials indicates the importance of the development of topology optimization methods
suited for electromagnetic design problems. To work with state of the art systems, an efficient
use of the available computational resources is essential. The purpose of this dissertation
is to investigate efficient electromagnetic analysis methods that can be combined with the
method of topology optimization, as well as the effective use of topology optimization in
solving electromagnetic design problems e.g., a good choice of design variables.
This dissertation consists of two parts. Chapters 2-4 focus on the investigation of a
topology optimization framework for design of periodic structures for electromagnetic applications, using a rigorous and efficient finite element analysis method. Chapters 5-6 discuss a
new design method that incorporates topology optimization and the concept of origami for
design of electromagnetic devices involving a complex 3D geometry. A short description of
each chapter is given below.
In chapter 2, the finite element framework for the electromagnetic analysis is discussed.
The focus is on the analysis of the electromagnetic wave propagation through periodic media.
The discussion extends to the implementation of a rigorous and efficient mesh truncation
technique based on plane wave expansions.
In chapter 3, topology optimization problems for design of two dimensional periodic dielectric structures are discussed. Two variations of the topology optimization problem solved
in [20] are solved using alternative formulations. The goal of the problem is to find periodic dielectric structures that exhibit desired electromagnetic wave transmission/reflection
characteristics in the frequency domain. The first variation uses a level set function for the
representation of the material distribution, and the second variation uses a different type of
mesh truncation method in the finite element formulation.
5

In chapter 4, topology optimization problems for design of three dimensional periodic
dielectric structures are discussed.
In chapter 5, a new approach to designing tunable electromagnetic devices is introduced,
which makes use of a transformation of complex surface geometry through folding. Sample
designs of frequency selective surfaces (FSSs) that can be tuned in the working frequency
through a folding and unfolding motion are provided.
In chapter 6, a method to design origami folding patterns based on topology optimization
is introduced. The mathematics of origami is used to carry out the analysis of the transformation of a foldable sheet from its flat configuration to the folded configuration. The
design method follows the “ground structure” approach of structural topology optimization,
where optimal folding patterns are found by assigning presence and type of folds to a set of
lines drawn on a two-dimensional domain. The origami design method is then applied to the
design of origami-based tunable FSSs. Folding patterns with desired geometric properties for
the design of tunable FSSs are obtained using the origami design method, and conducting
elements are added as a decoration to the foldable surface to achieve FSSs that can be folded
and unfolded to tune the working frequency.

6

Chapter 2
Finite Element Formulation for Time
Harmonic Electromagnetic Waves
Finite element methods can be used to analyze many complex systems and to efficiently
compute the sensitivity of the performance of the system with respect to the design, which
is essential in topology optimization. In this chapter a finite element formulation for the
analysis of electromagnetic wave propagation in the frequency domain is discussed, with an
emphasis on the mesh truncation method based on a plane wave expansion for periodic problems. There are books available on finite element analysis of time harmonic electromagnetic
waves with details of the formulation and implementation, for example, refer to [38]-[40].

2.1

Constitutive Equations

The electromagnetic wave behavior is governed by Maxwell’s equations:
7

∂B
∂t

(2.1)

∂D
+J
∂t

(2.2)

∇•D=ρ

(2.3)

∇•B=0

(2.4)

∂ρ
∂t

(2.5)

∇×E=−
∇×H=

∇•J=−

where E and H are the electric and magnetic fields, D and B are the electric and magnetic
flux density, J is the electric current density, and ρ is the electric charge density.
When the time dependency is assumed to take the form exp (jωt), Maxwell’s equations
reduce to vector wave equations

∇ × E + jωµH = 0
∇ × H − jωǫE = 0
where j =

√

(2.6)

−1 and ω is the angular frequency. Relevant material properties are described

by µ and ǫ, which denote magnetic permeability and electric permittivity, respectively. These
equations can be decoupled and solved as a problem of finding either the electric field or
magnetic field. If the electric field is taken as unknown, the governing equation becomes

∇×

1
jσ
2
∇ × E − k0 ǫr −
µr
ωǫ0

E=0

(2.7)

where k0 , σ, ǫr and µr are the wavenumber in vacuum, electrical conductivity, relative
8

permittivity and relative permeability such that ǫ = ǫr ǫ0 and µ = µr µ0 with

ǫ0 ≈ 8.854 × 10−12 F/m
µ0 ≈ 4π × 10−7 A · m

(2.8)

For a polarized wave, Maxwell’s equations can be simplified to a scalar Helmholtz equation. For a T M z (transverse magnetic to z) polarization, the z−component of the electric
field Ez can be taken as the unknown variable and found by solving the Helmholtz equation

2
∇2 Ez + k0 ǫr Ez = 0

(2.9)

Other electric field components can be expressed in terms of Ez and found by postprocessing.

2.2

Analysis of Periodic Media

Numerical analysis of electromagnetic wave propagation through periodic media has a great
importance as its applications are found in various areas such as modeling of radar, communication and sensing systems as well as design of engineered materials such as metamaterials.
Analysis setup of structures with 1-D and 2-D periodicity are described here.

2.2.1

Two-dimensional domain, 1-D periodicity

A dielectric medium with a periodicity in one direction (1-D periodicity) can be described by
a distribution of electric permittivity ǫ in a two-dimensional representative domain Ω shown
9

Ω

d

Ezin

x
y
Figure 2.1: A two-dimensional representative cell
in Fig.2.1. The material is assumed periodic in the x-direction with a tiling vector e =(d, 0).
A transverse magnetic polarization can be assumed with the electric field component Ez
taken as the unknown variable. The constitutive equation is the scalar Helmholtz equation
in Eq.2.9. To describe a field distribution in a periodic medium with a tiling vector e =(d, 0),
Ez is constrained to follow a Bloch-Floquet condition:

Ez (x + md, y) = Ez (x, y)e−jmα0 d ,

m∈Z

(2.10)

where α0 = k0 sinθ is the x-component of the propagation vector k of the incident field
Ezin = exp(−jα0 x).

2.2.2

Three-dimensional domain, 2-D periodicity

A medium with a periodicity in two directions (2-D periodicity) can be described by a
distribution of constituent materials in a three-dimensional representative domain Ω shown in
Fig.2.2. The structure is periodic in the x- and y-directions with tiling vectors ex =(T x, 0, 0)
and ey =(0, T y, 0). Constituent materials are described by values of electric permittivity ǫ,
magnetic permeability µ and conductivity σ.
For the 2-D periodic structure, the full-wave analysis is necessary. The electric field
10

Tx

Ty

Ω

z y
x
Ein
Figure 2.2: A three-dimensional representative cell
formulation of Maxwell’s equations can be adopted, and the vector wave equation in Eq.2.7
ˆ
ˆ
ˆ
can be solved for the electric field E = Ex x + Ey y + Ez z. The Bloch-Floquet condition for
the electric field in this case is

E x + mTx , y + nTy , z = E (x, y, z) exp −jmα0 Tx − jnβ0 Ty ,

m, n ∈ Z

(2.11)

where α0 and β0 are x− and y−components of the propagation vector k of the incident field
Ein = exp[−j(α0 x + β0 y)].

2.3

Finite Element Mesh Truncation Techniques

In the numerical formulation discussed above, an incident electric field is applied at the input
boundary and exits at the output boundary, where no physical boundary exists. To avoid
11

artificial reflections at the input and output boundaries, a mesh truncation technique with
appropriate boundary conditions must be used. Some of the frequently used mesh truncation
techniques are discussed in [39].

2.3.1

Absorbing boundary condition

A relatively simple method is to use an absorbing boundary condition, in which boundary
conditions ensure that any plane wave normally incident to the boundaries is absorbed.
Absorbing boundary condition is formulated as an approximation to Sommerfeld radiation condition as
n × (∇ × E) + jk0 n × (ˆ × E) ≈ 0
ˆ
ˆ
n

(2.12)

where n is the unit vector normal to the truncation boundary.
ˆ
When choosing absorbing boundary condition as the mesh truncation technique, the
corresponding (input and output) boundaries must be placed far enough from scatterer or
sources inside of the analysis domain. A rule of thumb is to place the input or output
boundary at least a half wavelength from any scatterer or sources.
For the scalar formulation with the Ez as the unknown variable, Eq.2.12 can be written
explicitly as
n • Ez + jkEz = jk (1 + cos θ) Ezin
ˆ

(2.13)

at the input boundary, with the incident electric field Ezin and

n • Ez + jkEz = 0
ˆ
at the output boundary.
12

(2.14)

PML
PML

PML
Ω

Ω
Ezin

PML
Ein

Figure 2.3: PML setup for 1-D and 2-D periodic structures

2.3.2

Perfectly matched layer

Another commonly used mesh truncation technique is to place layers of theoretical materials
called perfectly matched layers (PMLs) at boundaries. Examples of PML setups for 1-D and
2-D periodic structures are illustrated in Fig.2.3, where PMLs are placed at the input and
output boundaries, shown in gray.
A PML is a theoretical material designed so that reflection is reduced for a wide range of
frequencies, polarizations and angles of incidence. Reduction of absolute value of reflection
coefficient

|R (θ)| = exp −2kx

L
0

s (z)dzcosθ

(2.15)

L
can be achieved by increasing the integral 0 s (z)dz. For analysis with a wide range of
frequencies, choosing

s (z) =
13

σ (z)
ωǫ

(2.16)

will eliminate the frequency dependency of the PML. PMLs can be combined with a perfect
electric conductor boundary condition

n×E=0
ˆ

(2.17)

at the outer ends of the PML regions to truncate the mesh. In the case of a TM polarized
formulation with Ez as the variable, Eq.2.17 can be written simply as:

Ez = 0

2.3.3

(2.18)

Plane wave expansion

A rigorous mesh truncation technique is to use a plane wave expansion at the truncation
boundaries, a method discussed in [41] and used to solve an inverse problem.
In this approach, the electric field E at the input boundary is expressed as a combination
of the incident and reflected field as

E x, y, zin = Ein + Eref

(2.19)

ˆ
ˆ
z
where the incident field is written in terms of the propagation vector k0 = α0 x + β0 y + γ00ˆ
ˆ
and a constant unit vector p along the polarization direction of the electric field as

Ein = exp −j α0 x + β0 y + γ0 z − zin
The reflected field is expressed as a Fourier series as
14

p
ˆ

(2.20)

+∞
Eref =

+∞

m=−∞ n=−∞

rmn exp jγmn z − zin

exp [−j (αm x + βn y)],

n, m ∈ Z
(2.21)

At the output boundary, there is only the transmitted field

E (x, y, zout ) = Etr

(2.22)

which is also expressed as a Fourier series as

+∞

+∞

Etr =
m=−∞ n=−∞

tmn exp[−jγmn (z − zout )] exp[−j (αm x + βn y)],

n, m ∈ Z
(2.23)

In the three-dimensional formulation, the coefficients rmn and tmn are constant vectors,
and they become scalars for the two-dimensional formulation.
The propagation constant for each mode γmn can be found by




where

2
2
2
2
2
2
k0 − αm − βn ; αm + βn ≤ k0
γmn =

 −j α2 + β 2 − k 2 ; α2 + β 2 > k 2
m
n
m
n
0
0

2πm
αm = α0 −
Tx
2πn
βn = β0 −
Ty
15

(2.24)

(2.25)

It can be seen that only a finite number of modes are propagating modes, and higher modes
will be evanescent. This mesh truncation technique is rigorous since higher modes as well as
the dominant (m = n = 0) mode are included in the analysis. In particular, it is important
to include higher propagating modes, if there are any, as those modes can potentially have
a significant effect. All the propagating modes and some of the evanescent modes should be
included, but in practice any higher modes may be eliminated from the computation of the
sums in Eq.2.21 and 2.23, as they will have a negligible effect on the analysis accuracy.
In the two-dimensional formulation, a similar expansion can be used to truncate the
mesh. At the input boundary Ez can be expressed with reflection coefficients rm as

Ez (x, y) = exp[−jχ0 (y − yin )] exp(−jα0 x)
+

∞

m=−∞

rm exp[jχm (y − yin )] exp(−jαm x)

(2.26)

Similarly, at the output boundary, Ez can be expressed with transmission coefficients tm as

Ez (x, y) =

∞
m=−∞

tm exp[−jχm (y − yout )] exp(−jαm x)

(2.27)

Coordinates yin and yout are the y-coordinates of the incidence and exit boundaries, and
αm and χm are defined as

αm = α0 +
and
16

2πm
d

χm =

2.4





;


 −j

2
2
αm − k0

2
(nL k0 )2 − αm ≥ 0

;

2
(nL k0 )2 − αm

2
2
k0 − αm < 0,

(2.28)

m∈Z

Finite Element Formulation

The finite element equations can be obtained, following [42], by deriving the weak form
associated with Eq.2.7 or 2.9 and using an approximation

Ee ≈

Ne
i=1

Ne E e
i i

(2.29)

within each element e. The details of the derivation of the weak form and the finite element
equations are given in appendices A and B.
In the following chapters, finite element formulations will be used to analyze electromagnetic wave propagation in periodic structures.

17

Chapter 3
Topology Optimization of 2D Periodic
Dielectric Structures
In this chapter, topology optimization methods are used to design structures with periodicity in one direction. The goal is to find 1-D periodic structures that exhibit desirable
transmission and reflection characteristics by distributing two dielectric materials of distinct
electric permittivity values over a rectangular representative cell. The wave propagation is
along the orthogonal direction of the direction of periodicity, and the material is assumed
to be uniform in the third direction. For a polarized time-harmonic electromagnetic wave,
the analysis reduces to finding a two-dimensional distribution of one of the electric field
components by solving a scalar Helmholtz equation. This is a problem solved in [20] for
minimum power transmission using a density function to describe the design and the finite
element method as the analysis method. The problem can be solved as a minimization or
maximization problem of transmission (or reflection) to design structures that can be used
for many potential applications, e.g., frequency selective surfaces and radomes.
18

y=yin

Γb

y=yout

Γin

Γout
ΩD

d

Ezin

x
y

Γa

Ω

Figure 3.1: Design Domain
In this chapter, variations of the problem are solved using different design description
and analysis methods to investigate an efficient way to solve this relatively simple problem. In section 1, an approach based on a level set function to describe the distribution of
the constituent materials is used, combined with the finite element method with absorbing
boundary condition for the mesh truncation. The analysis method employed here is the
same as in [20]. In section 2, a standard density approach is used to describe the material
distribution, combined with the finite element method using the plane wave expansion at
the mesh truncation boundaries.

3.1

Problem Setup

The optimization problem is set up as follows. The objective of the problem is to find a 1-D
periodic structure, described by the distribution of a material with relative electric permittivity ǫa in a background material of relative electric permittivity ǫb , that minimizes/maximizes
the electromagnetic power flow. The problem is modeled on a 2D rectangular domain Ω
shown in Fig.3.1. An incident field Ezin enters the domain from the left at Γin , and the
power flow through the structure is measured at the right boundary Γout .
19

3.2

Level Set Approach

First, a level set function is used to represent the material distribution. Allaire, et al. [43]
and Wang, et al. [44] are among the first authors that proposed a formulation of topology
optimization problems based on level set methods. Level set methods use a so-called level
set function to assign a material at each location in the problem domain. While a “classical”
density approach permits the use of a mixture of two constituent materials in regions of
intermediate density (“gray” regions), a level set function is combined with a step function
to express material distribution, permitting no regions with intermediate material properties.
However, in computations, even in level set methods there will be mixtures of the two phases
near material interfaces, as a smoothed approximation of a step function is used. To facilitate
convergence to binary solutions while avoiding complex designs, the work here follows the
phase-field method introduced by Bourdin and Chambolle in [45], using a density approach,
combined with a level-set representation of the design.

3.2.1

Material distribution

A level set function Φ(x, y) : Ω → ℜ is used to describe the distribution of dielectric material
ǫr within domain Ω as

ǫr (Φ) = ǫb + H(Φ) ǫa − ǫb

(3.1)

Here H is the Heaviside function defined as

H(Φ) =



 1 if


 0 if
20

Φ≥0
Φ<0

(3.2)

In summary, we have ǫr = ǫa wherever the value of the function Φ is positive and ǫr = ǫb
wherever the value of is negative. The zero-level set of Φ defines the interfaces between ǫa
and ǫb . In computation, a smooth approximation of the Heaviside function is used to avoid
numerical difficulties.

3.2.2

Optimization problem

The goal of the optimization problem is to find a structure that minimizes power flow through
it. The problem is stated as follows: Find Φ(x, y) : Ω → ℜ that
minimizes
1
E E ∗ dΓ
J (Ez ) = 2 Γ
out z z

(3.3)

subject to
A (Φ) = Ω H (Φ) dΩ = A0
∗
where Ez is the complex conjugate of Ez , A (Φ) is the total area covered by material ǫa and
A0 is the prescribed area specifying the amount of material where ǫa is used.
A penalty term is added to the objective function to facilitate convergence to binary
solutions while avoiding rapid variations in Φ that may result in complex shapes and poor
¯
convergence. The augmented objective function J is defined as:

1
π
¯
J =J+
τ |∇Φ|2 + τ √ H (Φ) (1 − H (Φ)) dΩ
2 Ω
4 2

(3.4)

The penalty parameter τ is adjusted according to the relative importance of the penalty
term. A larger value of τ will result in a solution with simpler features or a slower variation
in a material distribution. The Lagrangian associated with the optimization problem using
21

¯
J as the objective is:

¯
L=J +λ

Ω

H (Φ) dΩ − A0

(3.5)

where a Lagrange multiplier λ associated with the are constraint is introduced. A gradientbased method is used to evolve the level set function Φ as

∂Φ
δL
∂ǫr
π
= −K
= −K β + λ
− τ ∇2 Φ + τ √ (δ(Φ) − 2δ(Φ)H(Φ))
∂t
δΦ
∂Φ
4 2

(3.6)

where K > 0 is a proportionality constant that may be used to control the step size, and β
is

2 ∂ǫr E W = k 2 (ǫ − ǫ ) δ(Φ)E W
β = k0
z
a
z
0 a
∂Φ

(3.7)

Here W is the solution of the adjoint problem

2
−∇2 W − k0 ǫr W = 0

(3.8)

n • ∇W + jk0 W = 0
ˆ

(3.9)

n • ∇W + jk0 W = Ez
ˆ

(3.10)

with boundary conditions:

on the input boundary Γin and

on the output boundary Γout , and δ is the Dirac-delta function.
22

Time evolution corresponds to an iterative scheme if the derivative of Φ is discretized as

1
∂Φ
≈
Φt+∆t − Φt
∂t
∆t

(3.11)

with time increment ∆t. The discretized version of Eq.3.6 is then:

1
Φt+∆t − Φt =
∆t
π
− K β + λ∆ǫδ Φt − τ ∇2 Φt+∆t + τ √ δ Φt H Φt
4 2

(3.12)

Collecting terms with Φt+∆t on the left-hand side, we obtain an equation to update Φt+∆t :

− α∇2 Φt+∆t + Φt+∆t =
απ
−K∆tβ − √ δ Φt − 2δ Φt H Φt
4 2

+ Φt

− λ K∆t∆ǫδ Φt

(3.13)

with the Dirichlet boundary condition:

Φt+∆t = c0

(3.14)

on the left and right boundaries of the design domain, where c0 is a constant, fixing the
material on the left and right boundaries of ΩD and the periodicity condition

∇Φt+∆t (x, y = 0) = Φt+∆t (x, y = d)
and α = Kτ ∆t controls the step size.
23

(3.15)

To solve Eq.3.13 for Φt+∆t , first f0 and f1 are defined as follows.

απ
f0 = −K∆tβ − √ δ Φt − 2δ Φt H Φt
4 2
f1 = K∆t∆ǫδ Φt

+ Φt
(3.16)

The solution of Eq.3.13 can be expressed as

Φt+∆t = Φ0 − λΦ1

(3.17)

where Φ0 and Φ1 are the solutions to:

− α∇2 Φt+∆t + Φt+∆t = f0
0
0
− α∇2 Φt+∆t + Φt+∆t = f1
1
1

(3.18)

The Lagrange multiplier λ is determined by finding λ that satisfies

A Φt+∆t = A0

3.2.3

(3.19)

Numerical examples

The governing equation Eq.2.9 is solved for the electric field component Ez . The problem
is modeled and solved using the finite element method in a commercial software COMSOL
[46] using the absorbing boundary condition, as described in Sec2.3.1. In all of the numer24

Table 3.1: Parameters for example 1
Parameter
d
D ǫa θ A0
Value
0.3λ0 4d 5.0 0° 40%
ical examples investigated here, vacuum is used as the background material i.e. ǫb = 1.
Homogeneous regions with ǫr = ǫb and depth 0.5d are placed to the right and left of the
design domain ΩD , so that any waves scattered off of a dielectric inclusion in ΩD may be
approximated as normally incident waves at Γin and Γout . The parameter τ is adjusted so
that the relative magnitude of the phase-field term and the electromagnetic power flow term
stays between 1/50 and 1/100. The parameter K∆t is adjusted according to the mesh size.

3.2.3.1

Example 1

This example uses parameters shown in Tab.3.2.3.1, where λ0 is the wavelength of the
incident electromagnetic wave in vacuum, corresponding to the prescribed incident frequency
of 100GHz. A structure with uniform permittivity is used as an initial guess. Normal
incidence θ = 0° is used, and 40% of the area is allowed to have the inclusion material
(ǫr = ǫa ).
The optimal structure and the electric field distribution through it are shown in Fig.3.2. In
Fig.3.2(b) the red and blue regions correspond to high and low electric field strength. The
power transmitted J is computed as in Eq.3.3, and compared to the value of J through
vacuum, J0 . The ratio

RJ =

J
J0

(3.20)

is used to report the performance of the optimum structure. The structure shown in Fig.3.2
achieves RJ = 2.24 × 10−3 .
25

(a) Optimum structure

(b) Electric field magnitude through the optimum structure
Figure 3.2: Optimal structure and the field distribution for example 1. For interpretation
of the references to color in this and all other figures, the reader is referred to the electronic
version of this thesis (or dissertation).
Table 3.2: Parameters for example 2
Parameter
d
D
ǫa
θ
A0
Value
0.3λ0 4d 11.56 0° 19.6%
3.2.3.2

Example 2

Parameters shown in Tab.3.2.3.2 are used for this example. Normal incidence is used, and
this time only 19% of the area is allowed to have the inclusion material (ǫr = ǫa ), but a
higher dielectric constant is used. A structure with uniform permittivity is used as an initial
guess again.
The optimal structure is shown in Fig.3.3, through which RJ = 1.04 × 10−4 is achieved.
The objective function J is plotted over a range of incident frequencies in Fig.3.4. The
plot shows that the solution has the characteristics of a band-gap structure, with a gap of low
transmission in a frequency range between 70-120GHz. The structure shown in Fig.3.3(a)
indeed resembles a known photonic band-gap, discussed by Sigmund and Hougaard in [47].
See Fig.3.5(a) and (b) for the photonic band-gap structure and its performance, computed
26

(a) Optimum structure

(b) Electric field magnitude through the optimum structure
Figure 3.3: Optimal structure and the field distribution for example 2.

0.5
0.4

J

0.3
0.2
0.1
0
50

100
150
Frqeuency (GHz)

200

Figure 3.4: Frequency sweep for the optimum structure for example 2.

27

(a) Band-gap structure
0.5
0.4

J

0.3
0.2
0.1
0
50

100
150
Frqeuency (GHz)

200

(b) Frequency sweep for the band-gap structure
Figure 3.5: A known band-gap structure and the frequency sweep.

using the parameters displayed in Tab.3.2.3.2.
Comparing the performance of the two structures, we observe that the solution obtained
through optimization performs better than the band-gap structure at the prescribed target
frequency of 100GHz (RJ = 2.44 × 10−3 for the band-gap structure). In this example,
the algorithm found a structure that resembles a photonic band-gap structure with high
performance. Looking at the structure in Fig.3.3(a), we also observe that the structure is
not symmetric across the middle of the design domain, even though the incident angle is 0°,
and the domain is periodic in the vertical direction. We suspect that this lack of symmetry
may be caused by small numerical discrepancies in the periodic boundary conditions applied
28

during the finite element analysis.

3.2.4

Conclusion

In computations, smoothing of the level set function is necessary when using a gradient-based
optimization algorithm, and in the end solutions still suffer from having mixture (gray)
regions at interfaces. It was also observed that several parameters affecting the analysis
accuracy and optimization algorithm need to be adjusted in an ad hoc manner but with
great care to avoid numerical instability.

3.3

Density Approach

Next, an interpolation scheme similar to the SIMP approach [48], [49], [50] is used to characterize the material distribution. This is a “classical” approach used in topology optimization
problems. To facilitate convergence a filter is used.

3.3.1

Material distribution

The relative permittivity at location (x, y) in Ω is expressed as

ǫr (ρ) = ǫb + (ǫa − ǫb )ρp

(3.21)

where ρ = ρ(x, y) ∈ [0, 1] is the effective density of material ǫa and p ≥ 1 is a coefficient
used to facilitate convergence. The effect of the parameter p is explained later. Upon
discretization, ρ is piece-wise constant, taking a constant value ρe within each element e.
The computational domain Ω may include a subdomain Ω0 where ρ is fixed, i.e., where the
29

material is not designed.

3.3.2

Optimization problem

The analysis of the electric field distribution is done using the plane wave expansion at the
input and output boundaries (Γin and Γout in Fig.3.1), as described in Sec.2.3.3. Details of
the finite element formulation are given in Appendix A. Using this formulation, one obtains
the modal transmission and reflection coefficients as a part of the solution. That is, the
solution to the finite element equation (Eq.3.22)








 KEE KEr KEt   Ez  

 


 

 KrE
dI
0  r  = 

 


 

KtE
0
dI
t



fE 


fr 


0

(3.22)

includes not only the field distribution Ez but the reflection r and transmission t coefficients
associated with all the modes included in the analysis. The objective of the optimization
problem can be defined by extracting the relevant reflection/transmission coefficients. The
optimization problem is stated as follows: Find ρ that

minimizes
J = J(z(ρ), z∗ (ρ)) = z∗T Az
subject to
0≤ρ≤1
where A is a diagonal matrix with constant coefficients
30

(3.23)

Akk = { 0, · · · , 0, a−Nm , · · · , a0 , · · · , aNm , b−Nm , · · · , b0 , · · · , bNm }
Nn terms
2Nm +1 terms
2Nm +1 terms

(3.24)

Matrix A is used to accommodate different design goals by selecting different values for the
coefficients ak and bk . For instance, setting ak = 1 and bk = 0 sets the objective function J
to be the sum of absolute values of reflection coefficients rm . If reflection coefficients of only
the propagating modes are included in the sum, the value of the objective function becomes
the reflection coefficient at the incident boundary. Similarly, setting ak = 0 and bk = 1
corresponds to minimization of transmission at the exit boundary.

3.3.3

Optimality condition

The Lagrangian L associated with (4.6) with non negative multipliers µ1 and µ2 is

L = J(ρ) + µ1 (ρ − 1) + µ2 (−ρ)

(3.25)

If ρ∗ is a local optimum, the KKT conditions require that for each design variable ρ∗
e

∂L
∂J
∗ = ∂ρ∗ + µ1 − µ2 = 0
∂ρe
e

(3.26)

If ρ∗ = 0, µ1 = 0 and µ2 > 0, then from (3.26) ∂J = µ2 > 0. If ρ∗ = 1, µ1 > 0 and
e
e
∂ρ∗
e
µ2 = 0, then ∂J = −µ1 < 0. Otherwise, µ1 = µ2 = 0 and ∂J = 0. In summary, if ρ∗ is
∂ρ∗
∂ρ∗
e
e
a local optimum, then for all e = 1, · · · , Ne
31

dJ > 0 if
ρe = 0
dρe
dJ < 0 if
ρe = 1
dρe
dJ = 0 if 0 < ρ < 1
e
dρe

(3.27)

where Ne is the number of elements.

3.3.4

Sensitivity analysis

Gradients of the objective function J are computed using a standard adjoint variable approach, starting from

dJ
∂J ∂z∗
∂z T ∂J T
∗
= ∗
+
( ) = λT Pe + λ∗T Pe
e
e
dρe
∂z ∂ρe ∂ρe ∂z

(3.28)

for all e = 1, · · · , Ne . The adjoint variable λe is obtained by solving the adjoint problem

KT λ = (

∂J T
) = AT z∗
∂z

(3.29)

Vector Pe is computed using

Pe = −

∂Ke
ze
∂ρe

(3.30)

Note that the portion of the stiffness matrix that corresponds to the area integral, KEE , is
the only component of K that depends on the design variable ρ and therefore
32




∂KEEe
0 0 
 ∂ρe


∂Ke 

=
0
0 0 


∂ρe


0
0 0
with
ij
∂KEEe
∂ρe

3.3.5

=

ij
∂KEEe ∂ǫr
∂ǫr

∂ρe

= k0 2 pρe p−1 (ǫa − ǫb )

Ωe

˜
j
φe i φe dxdy

Influence of parameter p

The factor p in Eq.4.5 does not have the same effect as in a typical topology optimization
problem with a volume constraint. Here the factor is used only to speed up convergence
to a binary solution. In numerical experiments, it is observed that convergence in regions
of elements with high values of ρ (near ρ = 1) is faster than convergence of elements with
low ρ values (near ρ = 0). What is also notable is that a strictly binary solution has better
performance for a minimum transmission problem than a solution of similar topology where
a thin layer of intermediate permittivity (0 < ρ < 1) is present at the interface between the
two phases. This can be explained using an analogy to camera lens design. Often, a camera
lens is coated with a matching layer of material with permittivity value between air and
the lens material to improve transmission of light into the lens. Without the intermediate
layer, the large mismatch between the indices of refraction would cause a strong reflection
at the interface. For a minimum transmission problem the opposite is true, a high contrast
is desirable and thus removing “gray” elements at the interface improves performance.
Based on those observations, the following scheme is implemented to speed up conver33

gence to a binary solution. The value p = 1 is used until the optimality conditions (Eq.3.27)
are satisfied. If the solution is not essentially binary (less than a small percentage of elements
is “gray”), p is increased e.g. up to p = 3. Iterations continue until a maximum number
of iterations is reached or until convergence to an essentially binary solution is achieved. In
most problems, setting p = 1 is sufficient. Occasionally faster convergence is achieved by
using p = 3. No instance was observed where the final solution depends on the value of p.

3.3.6

Numerical examples

The use of the approach is illustrated using simple examples that focus on different features
of the problem. A material with relative electric permittivity ǫa = 12 (Si) is distributed in a
background material with ǫb = 4.5 (SiO2 ). The permittivity of a vacuum is assumed for the
material outside of the analysis domain Ω. Without loss of generality, the magnitude of the
excitation |Ezin | is set to unity for convenience. Other parameters such as the wavelength
of the incident wave λ, angle of incidence θ, the dimension of the cell in the propagation
direction D, and the highest mode order Nm used in modal expansions are chosen in each
case to highlight different features of the problem or its solution. All examples are solved
using the MMA (Method of Moving Asymptotes) [51]. Except as noted, the initial design
has uniform permittivity with ρ = 0.2 throughout the design domain. A filter is used to
avoid complex designs and facilitate convergence to binary solutions. A “cone” filter [52],
[53] with rmin = 1.5 is used. Iterations start with p = 1 and may increase to p = 3.
34

Table 3.3: Parameters used in Example 1
D λ/d
θ
Nm
4d 3.7 40 °
9

3.3.6.1

Example 1

In this example we explore the solution of the problem solved in [20]. Slabs of non-designable
material of width d0 = 0.5d and fixed permittivity ǫb are placed at the incidence and exit
boundaries. The objective is to find a configuration that minimizes transmission, which
corresponds to setting ak = 0 in Eq.3.24. Coefficients bk are set as follows:

bk =



 1 if


 0 if

χk is real
χk is imaginary

Mode numbers k take integer values in [−Nm , Nm ]. χk can be computed for each mode k
as in Eq.A.4, based only on design independent parameters d, k0 , and θ. Note that a real
χk corresponds to a propagation mode and an imaginary χk corresponds to an evanescent
mode.
An electromagnetic wave with wavelength λ = 3.7d (k0 = 2π in Eq.2.9) is incident at the
λ
boundary Γin at θ = 40°. The highest mode order included in the wave expansion Nm is set
to 9. This is the maximum value allowed under the Nyquist condition so that the wavelength
of mode Nm is larger than the length of two elements. Values of other parameters used in
this example are shown in Table 3.3.6.1.
The optimum material layout and the iteration history are shown in Fig.3.6(a) and (b).
The solution is feasible and satisfies optimality conditions (Eq.3.27) within a prescribed
35

(a) Optimum structure

1
0.8

J0 = 8.78×10-1

J

0.6

J50 = 3.69×10-3

0.4
0.2
0
0

10

20

30

iteration

40

50

(b) Iteration history
Figure 3.6: Example 1. Optimum solution and iteration history

36

Figure 3.7: Example 1. Optimum solution with vacuum background material
tolerance, i.e., replacing Eq.3.27 by

dJ > 0 if
ρe ≤ δ
dρe
dJ < 0 if ρ ≥ 1 − δ
e
dρe
| dJ | ≤ δ
otherwise
dρe
with a threshold δ = 0.001. In this problem the same solution is obtained if only one term
is kept in the plane wave expansion, i.e., if the problem is solved with Nm = 0. This is
not surprising since higher modes become less significant when homogeneous slabs of length
0.5d are placed at the incidence and exit boundaries. Moreover, in the design in Fig.3.6(a)
all the material interfaces are parallel to the tiling vector. As a result, only the dominant
mode is present, i.e. there is no coupling into higher modes. This confirms the validity
of the single mode assumption used in [20] and [21]. However, the authors in [20] report
an instability in the optimization procedure when transmission at the boundary is used as
the objective function. Because of this instability the authors are forced to use a surrogate
objective function based on an area integral of the field. In the present work, however,
no such instability is observed, as verified by the iteration history in Fig.3.6(b). It should
also be noted that the material layout found here is qualitatively different from the solution
obtained in [20].
To investigate the influence of the host medium, the problem is now solved setting the
background material to vacuum. Figure 3.7 shows the optimum solution. The two have
37

Table 3.4: Parameters used in Example 2
D
λ/d
θ
Nm
2d 2.5 − 3.3 40 ° 10
qualitatively similar features (straight bands parallel to the tiling vector), but a much lower
objective function value of J = 4.98×10−5 (compared to J = 3.69×10−3 ) is achieved when
vacuum is used. The improvement in performance is due to the higher contrast in indices of
refraction, which causes stronger reflection.

3.3.6.2

Example 2

In this example we seek a configuration that minimizes transmission for wavelengths in the
range λ/d = 2.5 − 3.3. Material can be designed everywhere in Ω. The function J (i) =
(i)
z∗T A(i) z(i) is computed for each (λ/d)(i) in the interval [2.5,3.3] at 7 equally spaced

sample points. Terms are then added up to obtain the total objective J = 1 (J (1) + J (2) +
7
· · · + J (7) ). Coefficients ak and bk are set as in Example 1 for each (λ/d)(i) . Values for
other parameters used in this example are shown in Tab.3.3.6.2.
After 200 iterations, a solution that satisfies the optimality conditions with δ = 0.001 is
obtained. The solution has a small number of gray elements near the interfaces between the
two materials, which are removed by applying a simple threshold at ρ = 0.5 with negligible effect on performance. Figure 3.8(a) shows the solution after postprocessing. For this
solution the total transmission coefficient is calculated as

T =
˜
k

|t˜ |2
k

(3.31)

˜
where the sum is over all modal transmission coefficients of propagating modes k. T is
38

(a) Optimum solution after postprocessing
1
0.8

T

0.6
0.4
0.2
0
2.5

3.3

target

4
λ/d
λ /d = 2.5−3.3

5

(b) Wavelength sweep
Figure 3.8: Example 2. Solution and its wavelength sweep

39

plotted in Fig.3.8(b) for a range of wavelengths. Note that T = 1 corresponds to total
transmission while T = 0 corresponds to total reflection of the incident wave. A drop in
transmission is apparent in the range λ/d = 2.5 − 3.3, verifying that the solution has the
propagation characteristics set as the target of the optimization problem.

3.3.6.3

Example 3

In the preceding examples solutions achieve optimum performance through the use of material arranged in straight bands of varying widths and spacing parallel to the tiling vector.
One can change the number, width and spacing of the bands of material to achieve the
desired transmission characteristics. The preceding examples show that the optimization
algorithm essentially uses this strategy to find optimum combinations of the band parameters. Coincidentally, the configurations found have homogeneous slabs at the incidence and
exit boundaries, as well as a simple overall material distribution inside the computational
domain. As a result, a single mode assumption turns out to be sufficient in the analysis of
these problems.
In this example, an instance where the single mode assumption is not valid, is shown. This
solution is obtained by biasing the optimization algorithm with a non uniform initial guess.
A number of arbitrary starting layouts were tested, and a layout that yielded a qualitatively
different optimum layout was used. This starting design converges to a more complicated
material layout that requires a higher number of modes be included in the analysis.
In this example, the goal is to achieve a low transmission configuration at any angle of
incidence between 0°and 90°. A weighted sum of the objective function J evaluated over a
finite number of angles is used in computations. The objective function J (i) is computed for
40

Table 3.5: Parameters used in Example 3
D λ/d
θ
Nm
2d 3.0 0°, 20°, 40°, 45°, 60°, 80° 10

(a) Initial design

(b) Optimum solution after postprocessing
Figure 3.9: Example 3. Initial design and optimum solution after postprocessing
1
each angle θ(i) , and the total objective is calculated as J = N (J (1) + J (2) + · · · + J (Na ) ),
a
where Na is the number of angles included. It should be noted that the angles should be
selected with some care. For instance, the set 0°, 20°, 40°, 60°and 80° yields an undesirable
peak in transmission at 45°(the remedy is straight forward: simply add 45° to the set). The
coefficients ak and bk are set as in Example 2 when computing each J (i) . Values for other
parameters used in this example are shown in Tab.3.3.6.3.
The initial guess and the optimum solution after postprocessing are shown in Fig.3.9(a)
and (b). Again, the solution (before postprocessing) satisfies the optimality conditions with
a threshold δ = 0.001. Figure 3.10(a) shows the total transmission coefficient T averaged
over all angles of incidence, plotted for wavelengths in a range around the target wavelength
41

0.8
0.7
0.6

T

0.5
0.4
0.3
0.2

solution after threshold
optimum solution

0.1
0
2

3

4

5

λ/d
target λ/d = 3.0
(a) Wavelength sweep
0.03

solution after threshold
optimum solution

0.025

T

0.02
0.015
0.01
0.005
0
0

20

40

60

80

angle (degrees)
(b) Angle sweep
Figure 3.10: Example 3. Wavelength and angle sweeps (Nm = 10)

42

0.7
0.6
0.5

T

0.4
0.3
0.2
0.1
0
0

20

40
60
angle (degrees)

80

Figure 3.11: Example 3. Angle sweep computed using one mode only (Nm = 0)

λ = 3.0d. The solid and dashed lines show the wavelength response curves of the solution
before and after postprocessing, respectively. There is no significant difference between the
two curves, and in both cases transmission is suppressed at the target wavelength λ = 3.0d,
as desired. Objective function values are J = 1.32 × 10−2 and J = 1.28 × 10−2 , respectively
(The solution after postprocessing is better. Note that a filter was used and this prevents
the algorithm from reaching a binary solution). Figure 3.10(b) shows the response of the
optimum configuration for various angles. T is plotted in the range 0°-90° with λ = 3.0d for
the solution before (dashed) and after (solid) postprocessing. Again, no significant difference
is observed between the two designs. In both cases, even though transmission is not uniform
for all angles, it is suppressed over the entire range of interest 0°-90°.
With the more complex material distribution in Fig.3.9(b), one must be careful to account for higher modes in the analysis. Complexity of the material distribution, especially
inhomogeneity near the incident and exit boundaries, may result in significant deviations
43

Figure 3.12: Example 3. Optimum structure starting from a homogeneous material distribution
from the single mode assumption. This is illustrated by Fig.3.11, where the response for
different incidence angles is shown for the design in Fig.3.9(b), computed now using only
one mode (Nm = 0). The figure clearly shows that using only a single mode analysis the
predicted response is quite different from that obtained when Nm = 10 is used (Fig.3.10(b)).
Using a single mode expansion in this case could lead to convergence to erroneous results.
To conclude this example the problem is solved again, starting now from a homogeneous
initial design with ρ = 0.2. Figure 3.12 shows an optimum solution starting from such
design. This solution has an objective value of J = 1.08 × 10−2 , compared to the value
J = 1.28 × 10−2 obtained using the non uniform initial design. The new configuration in
Fig.3.12 differs from the one in Fig.3.9 in that it allows for analysis using a single mode
assumption, simply because of the homogeneous material placed at the incidence and exit
boundaries.

3.3.7

Designs with parallel bands

In most instances, layouts with bands parallel to the tiling vector were found as optimum
solutions. As the incoming wave arrives at the first “black” band, the magnitude of the
transmitted wave is reduced due to reflection (see Fig.3.13 (a)). At the next interface the
44

reflection
interference

Einc

(a) Bands parallel to the tiling vector

scattering

Einc

(b) Bands perpendicular to the incidence
Figure 3.13: Wave propagation through vertical and inclined bands

magnitude of the wave is reduced further due to reflection as well as destructive interference
with a wave that is reflected off of an adjacent interface. Going through multiple bands in
this manner, the transmitted wave at the exit boundary will be of reduced magnitude. This
explanation applies to any angle of incidence.
In contrast, a layout with slanted bands has sharp corners at the incidence boundary
(see Fig.3.13 (b)). As the incoming wave enters the domain, reflection at each band reduces
the magnitude of the wave as before. However, as the incident wave scatters at a corner,
the overall wave propagation behavior is significantly altered. In particular, electromagnetic
power flow into the direction parallel to the bands is allowed in the form of surface wave
modes, thus the wave may have high transmission at the exit boundary.
To investigate the effectiveness of a layout with bands parallel to the tiling vector, a
45

(a) Initial guess

(b) After 10 iterations

(c) Optimum structure
Figure 3.14: Optimization using slanted bands for the initial guess

46

4

path length in "black"
path length in "white"
θ=60°
θ=30°

L

d /d

3

2

θ=0°

1

0
0

20

40

60

80

angle (degrees)
Figure 3.15: Angle dependency of the path length dL
Table 3.6: Parameters used for the angle dependency
D λ/d
θ
Nm
1.9d 1.4 variable 11
minimum transmission problem is solved using as the initial guess a slanted grating perpendicular to the incident wave at angle of incidence θ = 45°(Fig.3.14 (a)). A layout with bands
parallel to the tiling vector is obtained (Fig. 3.14 (b) and (c)) as predicted from the earlier
discussion.
Even though qualitatively very similar designs are effective for any angle of incidence θ,
optimum solutions do depend on θ. The number of bands and the spacing between them are
adjusted depending on the angle of incidence. To illustrate this point, a problem is solved
with the parameters shown in Tab.3.3.7 for angles of incidence in the interval [0°,80°] with
a 5° increment. Lower electric permittivity values ǫb = 1 and ǫa = 1.5 are used to amplify
the angle dependency.
The optimum structures for θ = 0°, θ = 30°, θ = 60°are shown in Fig.3.15. The length
47

of the path that the refracted wave travels inside each band in the optimum layout, dL , is
plotted over angles of incidence in Fig.3.15. For this simple geometry, dL can be computed
analytically from the material properties and Snell’s law. The solid and dotted lines show
the path length inside “black” (ǫa = 1.5) and “white” (ǫb = 1) bands respectively. As the
angle of incidence increases, the spacing between the bands is adjusted and the path length
increases. The curves are continuous even though the number of bands in the optimum
solution is not. The number of bands decreases as the spacing between bands becomes too
large to fit within the domain (the number of bands depends on D).

3.3.8

Conclusion

For minimum transmission problems the solution strategy proposed is stable and robust,
even when only boundary terms are included in the objective. The formulation presented
can accommodate complex material distributions near boundaries by considering a variable
number of modes kept in the wave expansion at the boundaries. The problem has potentially
many local optima. Straight bands of material were found as optimum configurations whenever uniform initial configurations were used. In those cases a single mode approximation is
sufficient for accurate analysis: the same result would be obtained whether higher modes are
included or not. More complex material distributions appeared in optimum configurations
when iterations started from non uniform designs. In such cases it is crucial to include higher
modes in the analysis, as the single mode assumption could result in inaccurate evaluations
of performance. An advantage of the method introduced here is that it does not rely on an a
priori knowledge of how many modes should be kept to maintain accuracy. The optimization
method works equally well and with only marginally different effort whether one or more
48

modes are used in the analysis.

49

Chapter 4
Topology Optimization of 3D Periodic
Structures
In this chapter, a topology optimization method is used to design structures with periodicity
in two directions. The problem is setup in a similar manner as in chapter 3, for the design
of 1-D periodic structures. The goal of the design problem is to find periodic structures
that exhibit desirable transmission and reflection characteristics. Dielectric materials are
distributed in a 3D design domain by assigning an effective density at each location within
the domain.

4.1

Problem Setup

The optimization problem is set up as follows. The objective of the problem is to find a
2-D periodic structure, described by a distribution of a dielectric material within the design
domain shown ΩD in Fig.4.1, that minimizes/maximizes the electromagnetic power flow.
The electromagnetic propagation is analyzed within a 3D domain Ω, which is a representative
50

Tx

Ty

Γout

z=zout

ΩD

z y
z=zin

Γin

x
Ω

Ein
Figure 4.1: Design Domain

cell of a periodic structure with periodicity in the x− and y−directions with tiling vectors
ex =(T x, 0, 0) and ey =(0, T y, 0). An incident field Ein enters the domain from the bottom
at Γin , and the power flow through the structure is measured at the top boundary Γout .
The finite element method is used to carry out the analysis, with a plane wave expansion
at the truncation boundaries (Γin and Γout in Fig.4.1), as described in Sec.2.3.3. Details of
the finite element formulation are give in appendix B.
The electric field is assumed to take the form

E (x, y, z) = Ein + Eref
p
ˆ
= exp −j α0 x + β0 y + γ0 z − zin
+∞
+∞
+
rmn exp jγmn z − zin ψmn
m=−∞ n=−∞
51

(4.1)

at the input boundary Γin and

E (x, y, z) = Etr
+∞

+∞

=
m=−∞ n=−∞

tmn exp[−jγmn (z − zout )]ψmn

(4.2)

at the output boundary Γout . In Eqs.4.1 and 4.1 ψmn = exp [−j (αm x + βn y)], αm and βn
are x− and y−components of the propagation vector, and γmn is the propagation constant.
Infinite sums in the expansion are truncated in computation, as in the 2D structure case.
The modal transmission and reflection coefficients (tmn and rmn ) can be computed as
1
∗
Eψmn dΓ
Tx Ty Γout

(4.3)

1
∗
(E − Ein )ψmn dΓ
Tx Ty Γ
in

(4.4)

tmn =

and

rmn =

where tmn and rmn are complex constant vectors.

4.2

Material distribution

The relative permittivity at location (x, y, z) in ΩD is expressed as

ǫr (ρ) = ǫb + (ǫa − ǫb )ρ
52

(4.5)

in terms of relative permittivities of constituents, ǫa and ǫb . Here ρ = ρ(x, y, z) ∈ [0, 1] is the
effective density of the inclusion material ǫa . Upon discretization, ρ is piece-wise constant,
taking a constant value ρe within each element e.
The effective density ρ may be controlled through a piece-wise constant mapping X. One
choice is to consider a material distribution that only varies in the x− and y−directions and
is uniform in the “thickness” (z−) direction, so that ρ(x, y, z) = X(x, y). A more complex
case would be a variable material distribution in all three directions. Only the former
case is considered here, and the latter case is not investigated for its excess complexity in
computation and fabrication.

4.3

Optimization problem

The objective function of the optimization problem can be defined in terms of reflection or
transmission coefficients. The optimization problem is stated as follows: Find X(x, y) that

minimizes
J = J(S(ρ(X)), S∗ (ρ(X))) = S∗T AS

(4.6)

subject to
0 ≤ ρ(X) ≤ 1
where S refers to either tmn or rmn , and A is a constant diagonal matrix of size (2Nm +
1)(2Nn + 1) such that
53








A=





a1

0

0

0

0

a2

0

0

0

0

...

0

0

0

0

a(2N +1)(2N +1)
m
n













(4.7)

where Nm and Nn denote the highest modes included in the sums in Eqs.4.1 and 4.2. In
essence, A is used to collect only the modal reflection or transmission coefficients corresponding to the propagating (energy carrying) modes by setting

ak =



 1 ; if k corresponds to a propagating mode

(4.8)


 0 ; otherwise

in Eq.4.7. Note that for a mode with indices m and n, a real γmn indicates a propagation
mode and an imaginary γmn indicates an evanescent mode.

4.3.1

Sensitivity analysis

Gradients of the objective function J with respect to the element material effective density
ρe are computed using an adjoint variable approach as

∂J ∂S∗
∂S T ∂J T
dJ
∗
=
+
( ) = λT Pe + λ∗T Pe
e
e
∗ ∂ρe
dρe
∂S
∂ρe ∂S

(4.9)

dJ
= λT Pe + (λT Pe )∗
e
e
dρe

(4.10)

Adjoint variable λ can be obtained by solving the adjoint problem
54

∂J ∂S∗
∂J T
+
=
∂E
∂S∗ ∂E

KT λ =

∂S T
∂E

∂J T
∂S

T
(4.11)

and

∂J T
= AS
∂S

∂J
∗T
∗ = S A,
∂S

(4.12)

where

1
∂S
∗
=
φψmn dΓ
∂E
Tx Ty Γ
in

(4.13)

1
∂S
φψ ∗ dΓ
=
∂E
Tx Ty Γout mn

(4.14)

if S = rmn and

if S = tmn .
Vector Pe is computed using

Pe = −

∂Ke
2 p−1
˜
φe φe dΩEe
Ee = k0 pρe (ǫa − ǫb )
∂ρe
Ωe

(4.15)

Gradients with respect to the design variable X are computed using the chain rule as

dJ ∂ρ
dJ
=
dX
dρ ∂X
where
55

(4.16)

∂ρe
∂Xi

4.3.2



 1 ; if ρe is controlled by X
i
=

 0 ; otherwise

(4.17)

Projection of the effective density

To facilitate convergence to a binary solution, a projection technique is used [54], [55]. The
effective density ρ is projected to ρ using a smooth Heaviside function Hβ :
¯

ρ = Hβ (ρ) =
¯

tanh(βη) + tanh(β(ρ − η))
tanh(βη) + tanh(β(1 − η))

(4.18)

which approaches a step function as β −→ ∞. Parameter β is used to control smoothness
and, and η is used to control the threshold value. The projection function is plotted for
various values of β and η in Fig.4.2(a) and (b), respectively, to illustrate the effect of β and
η.
The projected density ρ is used to describe the material distribution in the analysis.
¯

4.3.3

Numerical example

The use of the approach is illustrated using a simple example. A material with relative electric permittivity ǫa = 2.25 is distributed in a vacuum background (ǫb = 1), within a design
domain ΩD that occupies the entire computational domain Ω. Domain Ω is described by
tiling vectors ex = (T x, 0, 0) = (1.5µm, 0, 0) and ey = (0, T y, 0) = (1.5µm, 0, 0) and thickness 1.0µm, in the z−direction, and discretized by 20x20x20 brick elements. The material
below the input boundary Γin in Fig.4.1 is vacuum, and the space above the output boundary Γout is assumed to be filled with the inclusion material with ǫa . Normal incidence with
the vacuum wavelength λ = 0.6µm = 0.4Tx and the magnitude |Ein | = 1 is applied. Values
56

1

0.8

ρ

0.6
β=1
β=5
β=10

0.4

0.2

0
0

0.2

0.4

0.6

0.8

1

ρ

(a) η=0.5
1

0.8

ρ

0.6

0.4
η=0.2
η=0.5
η=0.8

0.2

0
0

0.2

0.4

0.6

0.8

ρ

(b) β=5
Figure 4.2: Projection function

57

1

Table 4.1: Parameters used in the example
Parameter
Tx
Ty
λ
ǫa ǫb
Value
1.5µm 1.5µm 0.6µm 2.25 1
of the parameters used in the example are summarized in Table 4.1. The highest mode order
included in the wave expansion Nm is set to 2.
The initial design has uniform permittivity with ρ = 0.5 ∀(x, y) ∈ ΩD . The projection
function in Eq.4.18 is used to facilitate convergence to a binary solution. Parameter η is set
to 0.5, and β is set to 1 and is increased gradually during the optimization iterations.
The objective is to find a configuration that minimizes transmission, which corresponds
to setting S = tmn , and an optimal solution is found using the MMA [51].
The optimal solution, the 3x3 tile associated with the optimal solution and the iteration
history are shown in Fig.4.3(a), (b) and (c), respectively, where black color corresponds to
material with ǫa and white is vacuum in Fig.4.3(a) and (b). The peaks at 51, 101, 151, etc.
iterations seen in the iteration history in Fig.4.3(b) are due to the increased value of parameter β used in the projection function in Eq.4.18. As β increases, the projection function
becomes closer to a step function, allowing a smaller number of effective densities to take
intermediate values. The rate at which β changes should be adjusted carefully; increasing β
too slowly will require many iterations until convergence is achieved, and increasing β too
quickly will push intermediate effective densities (ρ) towards 0 or 1 too aggressively, leading
to poor convergence. Since Hβ for a large β has small gradients near ρ = 0 and ρ = 1,
once ρ is pushed to a value close to either 0 or 1, it will get “stuck” at that value. Parameter β is increased by multiplying the previously set value by 1.355 whenever the maximum
change in design variables is smaller than 0.01 and the change in the objective function is
smaller than 10 × 10−4 . The transmission coefficient of the binary solution obtained by
58

(a) Optimal solution

(b) 3x3 tile of the solution

0.7

Transmission coefficient

0.6
0.5
0.4
0.3
0.2
0.1
0
0

200

400
Iteration

600

800

(c) Iteration history
Figure 4.3: Optimal solution and iteration history

59

T = 0.128

T = 0.056

(a) Optimal solution

(b) 12 pixels altered

T = 0.077

T = 0.081

(d) 1 pixel altered, case1

(e) 1 pixel altered, case 2

T = 0.093

(c) 3 pixels altered
T = 0.076

(f) 1 pixel altered, case 3

Figure 4.4: Processed solutions
applying a threshold at ρ = 0.5 is shown for every 50 iterations, using an asterisk (∗) in
Fig.4.3(c). It can be observed that the objective value of the thresholded design approaches
that of the unprocessed design towards the end of the optimization process. After 800 iterations, β = 4.08 × 1011 , and we obtain an optimal solution with no “gray” regions, with the
transmission coefficient T = 0.056.
In Fig.4.3(a), we observe that several pixels are surrounded by pixels of a different color.
To evaluate the performance of simpler designs in the neighborhood of the optimal solution,
some designs are obtained by altering some of the pixels that are surrounded by a different
color, indicated by red squares, and the transmission coefficient is measured. The optimal
60

solution with its transmission coefficient is shown in Fig.4.4(a) again. The alternate designs
and their transmission coefficients are shown in Fig.4.4(b)-(d). It can be seen that as a
larger number of pixels is altered, performance becomes worse. To obtain simpler designs,
one needs to use a filtering technique during the optimization process instead of relying on
image processing after obtaining an optimal solution.

4.3.4

Conclusions

A topology optimization method was used to design a 3D structure with periodicity in the
x− and y−directions. The material distribution is controlled by a mapping that varies in
the xy-plane and uniform in the z−direction. A structure that minimizes transmission for
an incidence of a prescribed frequency was found. The effective density ρ was projected,
using a smooth Heaviside function, to ρ to be used in the analysis. Use of the parameter β
¯
to adjust the sharpness of the projection function needs great care, in order to achieve a fast
and stable convergence.

61

Chapter 5
Origami Tunable Surfaces for
Electromagnetic Applications
The use of the concept of origami for engineering applications is gaining interest as an integral
part of innovative design processes. Origami is a Japanese word for “art of paper folding”, a
tradition that appears in many cultures [56]. In origami design, the construction of a 2D or
3D geometry from a flat, usually square sheet is achieved through a combination of simple
folding steps.
Today, use of origami in engineering includes collapsible structures, structures that can be
folded to take up only a small fraction of the space taken by the final, unfolded configuration
of the structure. Examples from the aerospace industry include foldable satellite structures
such as solar panels [57][58].
In this chapter, the use of the concept of origami in design of electromagnetic devices
in [59] and [60] is discussed. In particular, shifting of resonance frequencies of frequency
selective surfaces and tunable metamaterials is discussed.
62

5.1

Origami Tunable Frequency Selective Surfaces

Frequency selective surfaces (FSSs) are used as band-pass or band-stop surfaces to filter
electromagnetic signals to enhance the operation of various electromagnetic systems such
as radar, communication and sensing systems. The capability of adjusting the working frequency of FSSs while in operation expands the utility of devices in advanced applications.
Tuning is commonly achieved using lumped components such as varactors [61], [62]. Another
method is to use special types of substrates such as liquid crystals [63]-[65] or ferrites [66],
wherein a 15-20% shift of a resonance frequency can be achieved by altering the properties
of the substrate through the application of an external excitation. Concepts based on changing the geometric configuration of metallic inclusions by means of switches have also been
proposed [67]. A drawback of such designs of tunable FSSs is that they require continuous
excitation from an external source, which increases the energy consumption of the system.
A new approach to designing tunable FSSs based on the concept of origami is introduced
in [59]. A layer of a conventional FSS, composed of a periodic array of conducting or
dielectric elements printed on a flat substrate, is folded into a periodic pattern in an origamilike fashion. The transformation of the surface geometry of the layer through a process of
folding and unfolding is applied to invoke tuning of the working frequency of FSSs. Such an
FSS requires an external energy input only at instances when the working frequency is to be
adjusted.

5.1.1

Miura-ori

An example of a folding pattern that can be used to design a tunable FSS is the well known
Miura-ori, shown in Fig.5.1(a). This pattern has one control parameter that determines a
63

Figure 5.1: Miura-ori

x
lx

y
z

α
β

b

ly

a
Figure 5.2: Unit cell of Miura-ori
folding state that ranges from flat to highly folded and remains periodic throughout the
motion. Varying this parameter changes the relative location of the prints of conducting
material, and thus the frequency of resonance. Computer simulations and experiments are
conducted to test the performance of the origami tunable FSS.
The unit cell of the chevron origami structure of Fig.5.1(a) consists of four parallelogramshaped facets, as shown in Fig.5.1(b). The cell is repeated periodically along two tiling
directions x and y with lengths lx and ly . The angle β determines the folding state. The
ˆ
ˆ
other geometric parameters, a, b and α, describe the parallelogram that forms the unit cell.
64

When β = 0° the sheet is flat. As β is increased, lengths lx and ly change, as do the
orientations of elements printed on the sheet. The length parameters are expressed in terms
of β as:

lx = 2a cos γ

(5.1)

ly = 2b(sin γ cos α + cos γ cos β sin α)
where γ = tan−1 [1/ (tan α cos β)] is an intermediate value defined for convenience. Each
facet (a parallelogram plane in case of Miura-ori) within the unit cell may be decorated, for
example, with an element made of a thin conducting print to make the folded surface into a
working FSS in the radio frequency (RF) range.

5.1.2

Performance evaluation

Performance of the designed FSSs is done by analyzing the electromagnetic wave propagation within the representative cell of a periodic structure and computing the reflection and
transmission coefficients (|S11 | and |S21 |) using a full-wave solver HFSS [68], based on finite
element analysis. An example of a setup is shown in Fig.5.3.
For Miura-ori, the projection of the representative cell follows a chevron shape. The
periodic boundary conditions are used at all six sides. At the input (top) and output (bottom) boundaries, PMLs described in 2.3.2 are used. A thin dielectric sheet as well as the
conducting prints are modeled using the impedance boundary condition.
65

Ein
PML

Conducting
print

Dielectric
sheet

z
y

PML
x

Figure 5.3: Analysis setup

66

w
Ro
b
Ri

a
Figure 5.4: Concentric ring conducting element

5.1.3

Choice of conducting element types

Variations of commonly used conducting prints are evaluated for their performance in terms
of the amount of shift in the resonance frequency and the strength of resonance.

5.1.3.1

Concentric double ring copper prints

First, concentric rings, shown in Fig.5.4 are considered.
Consider a design with a unit cell described by a = 25mm, b = 20mm and α = 60°.
Copper elements with Ro =5.45mm, Ri =4.8mm and w=0.4mm are placed on facets of the
folded sheet. The sheet is assumed to be very thin, and the dielectric constant of unity is
used in simulations. The FSS is illuminated by a normally-incident plane wave with the
electric field polarized along the x-axis and the propagation vector along the z-axis. The
transmission responses of an FSS in folded states described by β = 0°, 20°, 40° and 60° are
plotted in Fig.5.5. Computation of the transmission coefficient S21 was undertaken using
HFSS.
67

0
β = 0 deg
β = 20 deg
β = 40 deg
β = 60 deg

|S 21 (dB)|

-5

-10

-15
8.4

8.6

8.8
9
Frequency (GHz)

9.2

9.4

Figure 5.5: Simulated transmission coefficient |S21 | for concentric rings at normal incidence

68

Ca
Wb

b

Cb
Wa
a

Figure 5.6: Skewed cross-shaped conducting element
The amount of frequency shift that occurs as β is changed between angles β1 and β2 is
computed as
fβ1 − fβ2
∆f
=
fmid
0.5 fβ1 + fβ2

(5.2)

where fmid is the midpoint between frequencies fβ1 and fβ2 . Using this formula, the shift
in resonance observed in Fig.5.5 is 1.4%, as β is changed from 0° to 60°. The value of
transmission coefficient is kept under -10dB, except at β =40°.

5.1.3.2

Cross-shaped copper prints

Skewed cross prints, shown in Fig.5.6 are considered next, where the arms of the crosses are
aligned with the sides of the parallelogram.
Consider a design with a unit cell described by a = 19.6mm, b = 14.7mm and α = 45°.
Copper elements with Wa = Wb = 7.3mm and Ca = Cb = 2.2mm are placed onto the
facets of the cells. A sheet of paper is chosen as the substrate, and a dielectric constant of 3
and thickness 0.1mm are used in simulations. The FSS is illuminated again by a normallyincident plane wave with the electric field polarized along the x-axis and the propagation
vector along the z-axis. The transmission responses of an FSS in folded states described by
69

0

|S 21 (dB)|

-5

Sim.
Sim.
Sim.
Exp.
Exp.
Exp.

-10
-15

β=15 deg
β=30 deg
β=45 deg
β=15 deg
β=30 deg
β=45 deg

-20
-25

8

9

10
Frequency (GHz)

11

12

Figure 5.7: Simulated (Sim.) and measured (Exp.) transmission coefficient |S21 | at normal
incidence
β = 15°, 30° and 45° are plotted in Fig.5.7. As β increases, the resonance shifts continuously
to a higher frequency, with a shift of 19% occurring as β changes from 0° to 60°. The value
of transmission coefficient is kept well under -10dB at all values of β examined. Both in
terms of the amount of shift in resoance frequency and the strength of resonance, skewed
cross-shaped elements are more desirable, compared to concentric rings.

5.1.3.3

Summary

The frequency shift observed is due both to the changes in interactions between the incident field and the metallic structures caused by the change in unit cell size, and to mutual
interactions between the structures. The interactions and thus the resonance behavior are
70

highly dependent on the type of folding pattern, the shape of the conducting elements, and
the properties of the substrate. It has been shown that a sheet folded into Miura-ori and
decorated with concentric rings does not show much tunability. This is because with circular
rings, the location of the resonance depends primarily on the circumference of the rings, but
the spacing between the printed elements is not as important as with cross-shaped prints.

5.1.4

Angle dependency

To investigate the resonance frequency dependency on angle of incidence, the performance of
the folded FSS is evaluated for oblique incidence. The transmission responses of an oblique
incidence case are plotted in Fig.5.8. Here the original z-directed propagation vector is
rotated to lie in the y-z plane plane at an angle θ = 30° to the z-axis, with the electric
field still polarized along the x-axis. A similar trend as in the normal incidence case is
observed; the resonance frequency increases as the folding angle β increases. The three
curves shown in Fig. 5.9 summarize the dependency of resonance frequency on β for three
different angles of incidence θ= 0°, 30° and 60°. One can observe the same trend of increasing
resonance frequency with β in all three curves. It is worth mentioning that there is a band of
frequencies shared by all three curves. Within this frequency band, there is a β that brings
the resonance frequency to a desired location for any (of the three) angles of incidence. This
range of frequencies is indicated as a gray band in Fig.5.9.

5.1.5

Polarization dependency

Next, polarization dependency of a folded FSS is investigated. An incidence with the electric
field polarized along the y-axis is used in simulation. In this case, much weaker resonances
71

0
-5
Sim.
Sim.
Sim.
Exp.
Exp.
Exp.

|S 21 (dB)|

-10
-15
-20
-25

β=15 deg
β=30 deg
β=45 deg
β=15 deg
β=30 deg
β=45 deg

-30
-35

9

10

11
Frequency (GHz)

12

13

Figure 5.8: Simulated and measured transmission coefficient |S21 | at oblique incidence at
θ=30°

72

13.5
Sim.
Sim.
Sim.
Exp.
Exp.

Resonance frequency (GHz)

13
12.5
12

θ = 0 deg
θ = 30 deg
θ = 60 deg
θ = 0 deg
θ = 30 deg

10

20
30
40
Folding angle β in degrees

11.5
11
10.5
10
9.5

0

50

60

Figure 5.9: Simulated and measured resonance frequencies, fres in GHz

73

Figure 5.10: Fabricated origami tunable FSS
are observed, and therefore the usability of the folded FSS is significantly degraded for the
orthogonal polarization. A different design strategy using symmetric designs may resolve
this problem.

5.1.6

Experiment

A prototype FSS with an origami geometry identical to that used in simulations was constructed and its performance was measured to validate the origami tuning concept. The
folding crease pattern of a 40 by 40 element array was printed on a sheet of paper, and the
sheet was folded by hand, producing the chevron structure of Fig.5.1. When flat (β=0°),
the array measures 55.3 cm by 58.7 cm. Cross-shaped elements were chemically etched on
strips of copper tape, and one element was attached to each facet of the folded pattern.
The fabricated foldable FSS is shown in Fig.5.10. A more sophisticated fabrication method
relying on automated procedures is envisioned for eventual manufacture of foldable FSSs.
Measurement of the transmission properties of the prototype FSS was undertaken as
follows. The surface was folded to a desired angle β and attached to a foam sheet using
removable pins. For small values of β, clamps were used to hold the sheet in place, as
74

Figure 5.11: Experiment setup for normal incidence

shown in Figs.5.11 and 5.12. The sheet was placed on an adjustable stand in the center
of an anechoic chamber of size 12 by 12 by 24 feet (3.66 by 3.66 by 7.32 m). American
Electronics Laboratory model H-1498TEM-horn antennas were placed against opposite walls
and attached through rigid coaxial cable to an HP 8510C vector network analyzer, which
was used to measure the transmission coefficient S21 through the surface. Nylon dielectric
lenses were placed in front of the antennas to create a focused beam with a roughly Gaussian
profile at the position of the FSS; at 10 GHz the spot diameter of the beam is about 20 cm,
and the phase across the spot is approximately uniform. Focusing of the beam is important
to eliminate termination effects, such as edge diffraction, and interaction of the beam with
the supporting structures. Calibration was accomplished using a through-measurement with
the surface absent, first subtracting the background response measured with a metal plate
blocking the beam.
75

Figure 5.12: Experiment setup for oblique incidence at θ = 30°

Transmission responses measured using normal incidence (experimental setup shown in
Fig.5.11) and an oblique incidence angle of θ=30° (experimental setup shown in Fig.5.12) are
plotted in Fig.5.7 and Fig.5.8, respectively. The measured resonance frequencies for various
folding angles β are very nearly those predicted by simulation, as summarized in Fig.5.9.
Note that the measured resonance curves are not as sharp and pronounced as predicted
by simulations. This could be because the focused incident field used in the experiment
illuminates only a portion of the FSS array, and simulations assume the array is infinite and
plane-wave illuminated. Also, the loss in the paper used to construct the surface and in the
glue used to hold the crosses in place may lead to resonances with a lower quality factor than
predicted in simulation. Measurement at an incidence angle of θ=60° was not undertaken
because the FSS profile presented to the incident beam is too narrow at this angle to ignore
edge effects. Also, the orthogonal polarization state was not measured since it was shown in
76

simulation that in this case resonances are not strongly excited.

5.1.7

Conclusions

The feasibility of tuning an FSS by simple mechanical means is demonstrated using a Miuraori structure with cross-shaped metallic prints in simulations and experiments. By changing
the folding angle from 0° to 60°, the resonance frequency may be shifted by 19%. Tunability
performance was observed to be dependent on the type of conducting prints.

5.2

Origami Tunable Metamaterials

A well-known realization of metamaterials is achieved using SRRs, proposed by Pendry [27].
For instance, edge-coupled SRRs are formed by two concentric, coplanar split rings made of
a thin layer of conducting material printed on a dielectric substrate. It is possible to arrange
SRRs in a three-dimensional array to make a metamaterial that exhibits resonance behavior
when excited by an electromagnetic wave with a magnetic field oriented normally to the
plane of the rings. This resonance is enabled by a capacitance dominated by the proximity
of the inner and outer edges of the rings, and by a loop inductance. The frequency at which
the resonance occurs can be determined based on the geometry of the rings, the material
properties of the conducting prints and the substrate, and on the thickness of the substrate.
One disadvantage of metamaterials based on SRRs is a narrow bandwidth. To expand
the utility of these types of metamaterials, researchers have investigated various ways to tune
the resonance frequency of metamaterials in response to changes in environment or operating
mission. Similar to the design of tunable FSSs, incorporation of lumped components such
as varactors, diodes and potentiometers, in an array of resonating structures is a common
77

Figure 5.13: Corrugated sheet with SRRs, α = 0°
approach to achieve tunability [69]-[74]. For example, a varactor loaded split ring resonator
has been used to achieve a 30% shift in resonance in an RF range [70]. A wider range of
tuning (65%) has been achieved by using a combination of multiple switches and a varactor
placed in a spiral resonator [75]. Another approach is to use a host medium whose dielectric
properties can be altered [76]-[79], again, in a similar manner as in the design of FSSs.
A new approach to designing tunable metamaterials based on the concept of origami
([60]) is introduced here, where SRRs are printed on a surface folded into a periodic pattern
that can be un-folded by controlling one free folding parameter. The inner and outer rings are
printed on different surfaces so that when the folding parameter is varied, the gap between
the rings, and thus the capacitance of the resonators, is altered.

5.2.1

Design of an origami tunable metamaterial

A simple example of an arrangement that can be used for the design of tunable metamaterials
is illustrated in Fig.5.13. The folding pattern is that of a sheet folded in a corrugated
arrangement, which may be unfolded to form the surface shown in Fig.5.14.
78

Figure 5.14: Corrugated sheet with SRRs, α = 10°

A thin dielectric sheet is decorated with SRRs, with each ring in an SRR pair placed on
opposite faces of the corrugation. When the structure is completely folded, the rings lie on
the same plane and form an SRR centrally embedded within a dielectric sheet of thickness
2d. As the sheet unfolds, an air gap opens between the adjacent dielectric surfaces and the
rings separate. A unit cell containing a single ring pair is shown in Fig.5.15, and this cell
is repeated in the direction perpendicular to the corrugation. Thus the folding pattern is
periodic in one direction and uniform in the other, perpendicular direction. The pattern
has one free parameter, represented here by the angle α, that can be varied to change the
capacitance of the rings. When α = 0◦ , the array of SRRs is oriented vertically to the plane of
the corrugation. When the edges of the folded sheet are pulled apart along the corrugation
direction, the angle α increases (Fig.5.14) and the separation between the two split rings
increases. Excitation is by a plane wave incident from above such that the magnetic field is
aligned along the axis of the fully folded rings (α = 0◦ ).
79

Ro
G
Wg

Wr

w

Ri

w
Figure 5.15: Origami metamaterial unit cell design

5.2.2

Demonstration of tuning of a metamaterial

The full-wave solver HFSS was used to analyze the electromagnetic behavior of the origami
metamaterial, using periodic boundary conditions and PMLs as in Sec.5.1.2. An incident
plane wave with propagation direction perpendicular to the sheet was used to excite the
rings, with its magnetic field oriented normal to the plane of the rings in the α = 0◦
configuration. The transmission coefficients (S21 ) were computed using the field solution
found in simulations.
A sample design was analyzed where the resonance frequency of the fully folded (α = 0°)
structure is placed in the lower portion of the F-band, at 5.29 GHz. The unit cell parameters
from Fig.5.15 are given in Tab.5.1, while the spacing in Fig.5.13 is set to s = 4mm. Note
that Wg is the separation of the ring edges when α = 0◦ . Both rings were assigned the
conductivity of copper, and the substrate was chosen to be 0.5 mm thick with the properties
of Rogers Duroid 5870 (dielectric constant ǫr = 2.33 and loss tangent tan δ = 0.0012.)
Transmission and reflection coefficients are plotted in Figs.5.16 and 5.17, respectively, for
α=0°-10° in 2◦ increments. The transmission and reflection curves show that the resonance
80

Table 5.1: Parameters used in sample F-band design
Parameter Value (mm)
w
Ro
Ri
Wr
Wg
G

7
2.2
1.82
0.2
0.18
0.25

0

α=10°

−10

α=8°
α=6°

S

21

in dB

−5

−15

α=4°

−20

−25
4

α=2°

α=0°

5

6
7
Frequency in GHz

8

Figure 5.16: Transmission coefficient, |S21 |

81

9

0
−5
−10

−20

α=10°
α=8°
α=6°
α=4°

S

11

in dB

−15

−25
−30

α=2°

−35
−40
4

α=0°
5

6
7
Frequency in GHz

8

Figure 5.17: Reflection coefficient, |S11 |

82

9

Table 5.2: Resonance frequencies found using HFSS
α (deg) fr (GHz) fr (GHz) fr (GHz)
ǫr = 10.2 ǫr = 2.33
ǫr = 1
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

2.56
3.79
4.12
4.33
4.50
4.62
4.73
4.82
4.90
4.97
5.02
5.07
5.12
5.16
5.21
5.24

5.29
6.32
6.78
7.13
7.46
7.68
7.87
8.07
8.09
8.26
8.37
8.45
8.51
8.59
8.67
8.69

7.97
8.22
8.53
8.82
9.11
9.38
9.64
9.87
10.03
10.19
10.18
10.29
10.44
10.47
10.59
10.62

frequencies shift upwards as α is increased. Resonance frequencies for α between 0◦ and 10◦
are tabulated in Tab.5.2.
The shift in resonance frequency resulting from unfolding the sheet from its fully folded
state (α = 0◦ ) to an angle α can be quantified using the formula

δ(α) =

fr (α)
fr (α) − fr (α = 0)
=
−1
fr (α = 0)
fr (0)

(5.3)

A 2◦ change in α produces a resonance shift of 28%, a 5◦ change produces a shift of 45%,
and a 10◦ change produces a shift of 58%. Thus, the upward shift in resonance frequency is
quite rapid as the sheet begins to unfold, with the magnitude of the shift tapering off as α
increases. The result is a metamaterial that can be tuned mechanically over a wide band of
frequencies with only slight physical movement. The principle behind this effect is explored
83

in the next section.

5.2.3

Principle of operation

The proposed origami-based tunable metamaterial can be broadly tuned with small movements because the resonance frequency of the split-rings is highly sensitive to folding angle
when the angle is small. To understand the origin of this effect, recall that the resonance
frequency of a set of concentric split rings may be described using a simple series RLC circuit
model [80]. Using the circuit shown in Fig.5.18, the resonance frequency is computed using

fr =

1
2π Leq C

(5.4)

For a coplanar SRR, the inductance, Leq , is the inductance of a loop with a radius that
is the average of the two rings, r0 = Ri + Wg /2, while the capacitance, C, is due in part to
the capacitance between the ring edges and in part to the capacitance of the gaps. Usually,
the capacitance of the gaps is neglected, but the presence of the gaps causes the ring edges
to act like two capacitors in series, each with the capacitance of a half ring:

C=

Ceq Ceq
Ceq
=
Ceq + Ceq
2

(5.5)

where

Ceq = πr0 CP U L

(5.6)

Here CP U L is the per-unit-length capacitance of the ring pair, which depends on the properties of the dielectric substrate as well as the geometry of the rings. When the rings are used
84

Leq

Req
Ceq

Ceq

Figure 5.18: Equivalent circuit of an SRR
in the construction of an array of unit cells, Eq.5.4 can be used to compute the resonance
frequency of the resulting structure, if the effects of the mutual interaction between the array
elements are included.
When the sheet is completely folded (α = 0°), each unit cell consists of a set of coplanar
rings centrally embedded in a dielectric sheet. Because the capacitance is dominated by the
field lines extending between the adjacent edges of the rings, the majority of the electric flux
is within the dielectric, and the capacitance of a single set of rings may be approximated
by assuming that the rings are embedded in an unbounded dielectric. As the surface is
unfolded and the adjacent dielectric sheets move apart, the capacitance changes for two
distinct reasons. First, as the coplanar rings separate the capacitance decreases due to both
the change in alignment of the adjacent ring surfaces and to the increase in distance between
the edges. Second, as the dielectric surfaces unfold, an air gap appears between the surfaces,
and the electric flux begins to concentrate in this air gap region. This effect can be seen
in Figs.5.19-5.22, which are three-dimensional plots of the electric field strength computed
using HFSS for the parameters of Tab.5.1 and a folding angles of α = 0°-8°. Strong field
is shown in red color. The field is clearly strongest in the air region immediately adjacent
to rings, at a point 90◦ from the gaps in the rings. Because the dielectric constant of
85

Figure 5.19: Electric field magnitude in one unit cell for α = 0°

Figure 5.20: Electric field magnitude in one unit cell for α = 2°
air is significantly smaller than that of the substrate on which the rings are mounted, the
capacitance decreases rapidly as the extent of this region increases. An accurate computation
of how this effect compares with the reduction in capacitance due to caused by the increase
in distance between the rings requires a full wave simulation of the type already described
in Sec.5.2.2. However, it is possible to assess the significance of the contribution of the air
gap effect to the tunability of the structure using some simple analysis.
Of predominant importance is the change in resonance frequency with folding angle,
described by equation Eq.5.3. Assuming that the inductance is not seriously affected by
small changes in α, the change in fr may be described entirely in terms of the change in
capacitance by employing Eq.5.4:
86

Figure 5.21: Electric field magnitude in one unit cell for α = 8°

Figure 5.22: Electric field magnitude in one unit cell for α = 10°

87

δ(α) =

Ceq (α = 0)
−1
Ceq (α)

(5.7)

For small folding angles the array effect shouldn’t change significantly as α changes, so
the change in capacitance can be studied by observing a unit cell. To determine the effect
of the ring separation on capacitance, consider a set of concentric rings immersed in an
unbounded dielectric. As the rings unfold, an estimate of the capacitance of the inclined
rings can be made using the per-unit-length capacitance, CP U L , of the simple canonical
structure shown in Fig.5.23. Here two infinite parallel conducting strips of width Wr are
offset in both the x and y directions, and CP U L can be computed as a function of ∆ and
h using simple numerical techniques [81]. Ignoring the slight angle between the strips, the
capacitance of each half of the inclined rings shown in Fig.5.24 can be found by assuming
that ∆ = Wr + Wg and that CP U L only depends on h:

Ceq (α) =

π
0

CP U L (h)r0 dφ

(5.8)

where
h(φ, α) =

w
+ (2Ro − 2Wr − Wg ) cos φ cos α.
2

(5.9)

Fig.5.25 shows plots of δ(α) found using Eq.5.8 in Eq.5.7 for the geometrical parameters
in Tab. 5.1. It is seen that as α is increased from zero, δ(α) increases at a fairly constant rate
due to decreasing capacitance. In comparison, when δ(α) is plotted from Eq.5.3 by using the
resonance frequencies from the HFSS simulations described in Section 5.2.2, a much more
rapid increase is seen for small α. This suggests that at small folding angles, the decrease in
capacitance is dominated by the appearance of the air gap rather than the separation of the
88

y

Wr

h
Wr

x
∆

Figure 5.23: Geometry of offset parallel conducting strips

Figure 5.24: Geometry of unfolded origami metamaterial unit cell

89

1

Full-wave

εr=10.2

δ

0.8

Full-wave

εr=2.33

0.6

Full-wave
no dielectric

0.4
0.2

Capacitance model
0
0

5

10

15

α
Figure 5.25: Relative change in resonance frequency (dashed curves are best-fit lines)

rings. However, for larger α most of the electric flux resides in the air gap, and thus as α
is increased from larger values, the rate of increase in resonance frequency should approach
that of two rings separating in free space. Observing the two curves it is seen that the slopes
begin to resemble each other as α approaches about 8 degrees. To verify this, the derivative
of δ(α) with respect to α was computed for both the capacitance model and the full-wave
analysis by finding the slope of the δ(α) curves. Because slight variations in the HFSS results
preclude using a simple finite-differences derivative, the data was first fit using the empirical
model

4
δ(α) =

An α1/n

(5.10)

n=1
and then this expression was analytically differentiated. The resulting derivatives are plotted
in Fig.5.26. Clearly the relative changes in resonance frequency begin to coalesce above
90

0.15
Full-wave
εr=10.2

δ‘

0.1

Full-wave
0.05 εr=2.33

Capacitance model

0
0

2

4

6

8

10

12

14

α
Figure 5.26: Derivative of relative change in resonance frequency
α = 8◦ , and thus above this angle the change in resonance frequency is due primarily to the
separation of the rings and not to the increasing air gap.
To study the capacitance effect further, the HFSS simulations were repeated for the
same geometrical configuration as shown in Tab.5.1 but with the dielectric parameters of
the substrate set first to those of free space (ǫr = 1, tan δ = 0) and then to those of
Rogers Duroid 6010 (ǫr = 10.2, tan δ = 0.0023.) The resulting resonance frequencies are
shown in Tab.5.2, and δ(α) found using Eq.5.3 is plotted in Fig.5.25. When the dielectric
constant is set to that of free space, the effect of the opening air gap is not present, and the
increase in resonance frequency is due entirely to the decrease in capacitance caused by the
separation of the rings. The resulting change in resonance frequency, shown as the crosses
in Fig.5.25, follows closely that predicted by the capacitance model for rings opening in free
space, even though the resonance frequencies found from HFSS with ǫr = 1 are significantly
91

higher than those found by HFSS with ǫr = 2.33. This result also validates the assumption
that the resonance behavior of the unfolding origami metamaterial array can be deduced by
observing a single unit cell, making it easier to devise the experiment described in Sec.5.2.4.
Finally, when the dielectric constant is increased to ǫr = 10.2 from ǫr = 2.33, the contrast
between the permittivity of the substrate and that of the air gap that opens as the rings are
unfolded increases, and the impact of the opening air gap on the capacitance of the system
is magnified. This results in a more rapidly rising value of δ with increasing α at low values
of angle, as is clearly seen in Fig.5.25, and a greater overall fractional change in resonance
frequency at any value of α. Even so, as with the ǫr = 2.33 substrate, the rate of change
levels out and becomes the same as two rings opening in free space, as shown in Fig.5.26.

5.2.4

Experiment

An experiment was conducted to test the validity of the design concept. Two copper rings
were etched on separate 0.38 mm-thick sheets of Rogers Duroid 5870 (dielectric constant
ǫr = 2.33 and loss tangent tan δ = 0.0012) as in Fig.5.27 to form an origami metamaterial
unit cell when placed together. The unit cell has dimensions w = 6.85, and the etched rings
have the geometry values shown in Tab.5.1. The dielectric sheets were glued along an edge
to form a hinge, and the sheets were opened to a set angle α by placing a small spacer
between the sheets at the open edge. The origami metamaterial unit cell was then placed
onto a piece of Styrofoam and inserted into a section of WR-159 F-band (4.90-7.05 GHz)
waveguide (Fig.5.28.) The hinge was oriented vertically in the guide so that the horizontal
magnetic field of the dominant T E10 mode is along the axes of the rings, thus exciting
the ring resonance. Finally, the transmission parameters |S21 | of the waveguide section
92

Figure 5.27: Two halves of fabricated origami metamaterial unit cell with the unit cell width
w = 6.85 mm

Figure 5.28: Origami metamaterial unit cell placed into an F-band waveguide sample holder
were measured to determine the frequency of resonance. Because the unit cell behaves
as a negative permeability metamaterial near resonance, the propagation constant of the
dominant mode becomes imaginary and the wave becomes evanescent. Resonance is thus
indicated by a stop band with a dip in |S21 |.
Fig.5.29 shows the measured values of |S21 | for several values of folding angle α. It is
clearly seen that as the unit cell is unfolded, the resonance frequency quickly increases. A
plot of resonance frequency versus α is shown in Fig.5.30, indicating that the resonance
93

frequency may be shifted across nearly all of F-band with a slight increase in α from 0◦
to 2◦ . As a comparison, the unit cell placed into a waveguide with perfectly conducting
walls was simulated using HFSS, and values of |S21 | are shown for several values of α in
Fig.5.31. A similar trend upward in resonance frequency with increasing α may be seen
in this figure. The resonance frequencies from the simulations are plotted in Fig.5.30 and
compared to the measured resonance frequencies. Interestingly, it appears that the measured
resonance frequencies increase with α at a rate even greater than predicted by simulation.
This is probably due to a slight unintended gap near the hinge which exists even for very
small α. As the sheets are unfolded, this small gap becomes less important compared to
the gap opening between the sheets, and the curves agree more closely. The case of α = 0°
was measured last, with the two sheets glued togther to ensure no gap. At this angle the
resonance frequencies from experiment and simulation are nearly identical.

5.2.5

Conclusions

The reflection and transmission characteristics of a corrugated sheet decorated with SRRs
were investigated using full-wave simulations. When the folding angle α is increased, resonance is shown to shift rapidly to higher frequencies. A simple analysis reveals that the
dominant effect on resonance at small values of α is the decrease in capacitance due to the
air gap that opens between the rings as α is increased from zero. Experimental results using a unit cell placed into a waveguide section verify this effect. These results suggest the
possibility of tuning folded metamaterial structures using slight mechanical manipulations.

94

1
0

|S 21 | (dB)

−1
−2

α=0°
α=1.84°

−3

α=0.16°

−4

α=0.49°

−5
−6
5

5.5

6
6.5
frequency (GHz)

α=0.94°

7

Figure 5.29: Measured transmission coefficients for an origami metamaterial unit cell placed
into an F-band waveguide

95

7.5

F-band waveguide measurement

resonance freq (GHz)

7

6.5

F-band waveguide simulation
6

5.5

5
0

1

2
α (deg)

3

4

Figure 5.30: Resonance frequency of an origami metamaterial unit cell placed into an F-band
waveguide

96

0

−1

|S 21 | (dB)

−2

−3

α=0°
α=0.25°
α=0.5°
α=1°
α=1.5°

−4

−5

−6
5

5.5

6
6.5
frequency (GHz)

α=2°
α=2.5°
α=3°
α=3.5°
α=4°
7

7.5

Figure 5.31: Transmission coefficients for an origami metamaterial unit cell placed into an
F-band waveguide from HFSS simulations

97

Chapter 6
Topology Optimization for Origami
Design
The use of the folding and unfolding motion of origami for resonance frequency tuning in
electromagnetic applications was introduced in chapter 5. In the designs used in chapter
5 and most of other applications of origami in engineering design, such as self-assembly of
microdevices such as biomedical devices [82], electronics [83], and microfluidics [84], simple
folding patterns with well-known geometric properties are used.
The development of design algorithms and the mathematics of origami in the last two
decades allows for systematic design of quite intricate folding patterns [56], [85]-[91]. For
example, Lang ( [86] [87]) developed a computer program where the origami design is decomposed into two processes: design of a base and design of the complete model. An origami
designer draws a stick figure in the program indicating, e.g., the number and lengths of legs.
The algorithm then finds a crease pattern that provides the desired base, and the designer
completes the design using artistic skills. This is called the tree method. Such methods
98

are commonly used for art work that may involve very complex geometry. Using a different
approach, Tachi [89] [90] developed an origami design method aimed at industrial applications based on tuck-folding. The method systematically finds a folding pattern that follows
a prescribed polyhedral surface by tucking some part of the sheet under another to adjust
the surface geometry.

The work presented in this chapter introduces an optimization-based method for origami
design. The method used here is similar to an origami design technique that involves the use
of a sheet of paper with lines drawn on it. The lines act as a guideline for origami designers
to decide where to create folds; some of the lines become folds and some remain flat. The
idea is that, if one starts with a sheet with many lines, many different origami designs
can be constructed by creating folds along some of those lines. This is analogous to the
ground-structure approach to topology optimization of truss structures [92] [93] . The idea
there is to start from a ground structure, a potential structure that includes a sufficiently
rich set of truss elements, and eliminate unnecessary elements based on an optimization
algorithm. In this work we present an origami design counterpart of a ground structure for
truss design, constructed by drawing a set of lines on a 2D domain, and an optimization
strategy to find patterns with folds along these lines that result in a folded geometry with
desirable, target geometric properties. Examples of target properties include a prescribed
distance between two points on the sheet and a prescribed angle between two planes, in the
folded configuration. Origami designs considered here are restricted to those that may be
constructed using only simple folds, namely, those made by folding a sheet along a line to
make a mountain or a valley shape. More advanced folding techniques, such as sliding a flap
into a pocket, are not considered.
99

The proposed origami design method can be incorporated into the existing folding design
problems such as the design of deployable structures, packaging materials, self-assembled microdevices and dynamically alterable electromagnetic devices to improve their performance.
Using sheet materials with customized properties, which may be achieved through the inclusion of nano-particles, the use of the origami deign method can be extended to material
design, e.g., for the design of thermoelectric materials, membranes, peristaltic porous media,
battery or fuel-cell microarchitectures.

6.1

Construction of the Ground Structure

There are several types of grids traditionally used as a guideline for origami design. One
simple example is a square or triangular grid system, used in [94] for design of origami
tessellations, formed by a repetition of folded patterns that make up a flat or curved surface.
For instance, an origami waterbomb, an example of origami tessellation that follows a curved
surface, can be constructed using the triangular grid system shown in Fig.6.1(a). A grid
system that results in a more intricate design is the angular grid system studied in [91]. This
grid is constructed by drawing lines through reference points, defined as either the corners
of the sheet or intersections of lines already drawn. At each reference point, 2n lines with
a separation angle 90°/n are drawn, creating new intersection points. Fig.6.1(b) shows an
example of angular grid systems for n = 4, where lines at each reference point are drawn
22.5° apart from each other. Maekawa constructed very complex and popular origami designs
such as a peacock, a dinosaur and a devil (shown in [95] with instructions), using this type
of grid system.
In this work, ground structures similar to those used in structural design are used. A
100

(a) Triangular grid system

(b) 22.5° grid system

Figure 6.1: Origami grid systems

(a)

(b)

(c)

(d)

Figure 6.2: Variations of ground structures

101

ground structure is constructed by placing points along the boundary of the design domain
and drawing connections between all point pairs, creating new vertices at intersections. This
is a simple method that provides full control of where vertices appear at the boundary. It
is advantageous to have such control for origami design problems where constraints (e.g.,
periodic boundary conditions for tessellations) are applied at the boundary. Using this technique, one can construct symmetric ground structures such as the one shown in Fig.6.2(a).
This structure is appropriate for design of origami tessellations, since it can be tiled in the
horizontal and vertical directions with vertices and lines consistent with the periodicity assumption at the boundary. Asymmetric ground structures can be constructed as well, for
example, by adding a point to the grid shown in Fig.6.2(a), to construct a structure shown
in Fig.6.2(b). The ground structure can be refined to accommodate more complex designs
by adding points as in Fig.6.2(c). Ground structures on a non square domain can also be
created, as in Fig.6.2(d). Miura-ori, an origami tessellation used for engineering applications
in [57][59][96] and several other engineering designs found in literature, can be constructed
from the ground structure shown in Fig.6.2(a).

6.2

Crease Type Assignment and a Folded State

Here we consider a rigid origami. A rigid origami is a developable surface made of origami
facets and foldlines that can be replaced by rigid panels and hinges, respectively, and can
be flattened without distorting the sheet. The model in [88] is used to analyze the geometry
of folding patterns in rigid origami. The process of creating a fold in this model can be
considered as a rotation of a facet by an angle with respect to its adjacent facet. This
angle of rotation is called a folding angle. Fig.6.3(a) shows a circular section of a ground
102

ρ<0

ρ>0
ρ<0
ρ=0
ρ<0
(a) Flat

(b) Folded

Figure 6.3: Flat and folded states of a single-vertex crease

structure around a vertex, called a single-vertex crease, with n = 5 lines extending from
it. An example of a folded geometry is shown in Fig.6.3(b) for a given set of n = 5 folding
angles, denoted ρ. The type of fold at each segment is determined by the value of folding
angle: zero, negative and positive folding angles correspond to flat, mountain and valley
fold, respectively, as shown in Fig.6.3. Assignment of fold type on the entire domain is done
by assigning folding angles to all segments in the ground structure.

The folding angles will be used as design variables to find a folding pattern that meets
the target geometric properties. A ground structure with folding angles assigned to all
segments represents a folded state: the configuration of the folded geometry of a sheet.
One can process information encoded in the folding angles and connectivity of the ground
structure to find the geometry of the folded surface. Geometric features relevant to the
specific design problem can then be extracted and used within an optimization framework
to search for folding patterns with the target properties. The use of continuous variables as
design variables allows the use of gradient based optimization algorithms.
103

6.3

Foldability Conditions

An important constraint in origami design is that the designed pattern is in fact foldable.
Necessary conditions for foldable origami design were investigated by belcastro and Hull in
[85] and require that the sheet does not stretch or rip during the folding process, and that
each face remain flat. The condition to avoid a sheet intersecting itself and penetrating to the
other side (self-intersection) is not addressed in [85] (this is still an open problem in origami
design). The necessary conditions for foldability in [85] were used to simulate a motion of
rigid origami in [88]. A similar approach is used in this work. The foldability condition of a
crease around the k th vertex, having n crease lines, is expressed as

k
k
k
Fk ρk , θk = Rρk Rθ1 Rρk Rθ2 · · · Rρk Rθn = I
n
1
2

(6.1)

Matrices Rθk ’s and Rρk ’s are rotation matrices

k
Rθi

Rρ k
i

and

k k
k
θ1 , θ2 , · · · θn


k − sin θk 0

cos θi
i




 sin θk cos θk 0
=

i
i




0
0
1


0
0 
1




= 0 cos ρk − sin ρk 

i
i




k cos ρk
0 sin ρi
i


(6.2)

(6.3)

are angles of rotation about the axis perpendicular to the (flat) sheet,

measured between adjacent crease lines, while ρk , ρk , · · · ρk
n
1 2
104

are folding angles, as shown

ρ3
ρ2
l2

l3
θ3

θ2
θ1
θ4

l4

l1

ρ4
ρ1
Figure 6.4: Single-vertex crease

in Fig.6.4. For a multi-vertex crease Eq.6.1 has to be satisfied at all vertices k = 1, 2, · · · M ,
where M is the number of vertices. Note that the notation ρk refers to the ith crease around
i
vertex k i.e. indexing used in Eq.6.1 is local. It has to be mapped to the global indexing
going from 1 to N , where N is the total number of design variables. For design of origami
tessellations, foldability conditions for vertices along the edge need to account for crease lines
connected to their periodic pairs in neighboring tiles.

6.4

Optimization of Origami Design

The goal of the optimization problem is to find a combination of folding angles that results
in a folded state with target geometric properties. Examples of target properties include
a prescribed distance between two points on the sheet and a prescribed angle between two
planes, in the folded configuration.
105

6.4.1

Geometric properties

Each geometric property to be controlled is expressed in terms of the coordinates X of
vertices in a folded state as

J = J (X)

(6.4)

The coordinates X in a folded state are computed for a given combination of folding angles.
For instance, the orientation of the crease line l2 in a single-vertex crease shown in Fig.6.4
can be found through a rotation of l1 about the axis of the circle by θ1 , followed by a rotation
about the axis along l1 by ρ1 . One can work around one vertex at a time, computing the
coordinates of crease ends relative to the vertex. This process can be repeated for all vertices
to find X.

6.4.2

Objective function

A simpler design has fewer foldlines and a larger number of “off” foldlines that remain flat
(ρ = 0) throughout the folding process. It is our interest to design a folding pattern that
achieves the target geometric properties using only a small number of “on” foldlines. This
way, any device that is designed based on the folding pattern requires a small number of
components and less effort during the fabrication process.
To find folding patterns with a small number of “on” foldlines, one may use in the
optimization problem an objective function that favors designs with larger number of “off”
foldlines. One such function is f , defined as
106

ρ2w(ρ)
1
a=1.0
0.5

a=0.5

0

−π

π

ρ

Figure 6.5: Objective function f

f=
i

ρ2 w ρi
i

(6.5)

Here w ρi is a weight function that takes the form
2
w ρi = Ce− ρi /a

(6.6)

with a constant C = e/a2 and a parameter a ∈ (0, ∞). Curves in Fig.6.5 show ρ2 w(ρ) for
a = 0.5 and a = 1. As maxima occur at ρ = ±a, their location can be adjusted by choosing
a. The objective function penalizes folding angles near ρ = ±a, forcing ρ away from ±a,
towards zero or towards larger values.

6.4.3

Optimization problem

If a folding process is considered as a time evolution of geometry of a sheet, an optimal
solution should describe the geometry of a folded sheet at one instant during the folding
process. To specify the instant at which the geometric properties of a folded sheet are
examined, one of the folding angles ρr is prescribed to take a fixed value ρ0 , r ∈ {1, 2, · · · N }.
r
107

After ρr is fixed, the optimization problem is
r
Optimization Problem P0 :
r
ρi for i ∈ U0 that

Find

Minimizes f (ρ)
Subject to gk = 0 for k = 1, 2, · · · M

(6.7)

hi = 0 for i = 1, 2, · · · Meq
−π ≤ ρi ≤ π

r
for ∀i ∈ U0

r
r
where U0 = {1, 2, · · · N } \ r is the set of indices of the free folding angles in P0 . In Problem
r
P0 ,

• Constraints gk express the foldability condition. gk is a measure of the deviation of
a given design from a foldable design, computed by collecting the three independent
components of a 3x3 matrix Gk as

T
g k = Gk , Gk , Gk
23 31 12

(6.8)

where

Gk =
ij

2
1
Fk ρk , θk − Iij ,
ij
2

i, j = 1, 2, 3

(6.9)

• Constraints hi are used to express how well the design meets a target property Ji .
They are defined by setting
108

hi =

1
∗
J − Ji 2
2 i

(6.10)

∗
for each geometric specification Ji and corresponding target value Ji .

r
Clearly, there is no guarantee that a feasible solution to P0 exists for a given target
r
property and fixed ground structure. If no feasible solution to P0 is found for any fixed
angle (r and ρ0 ), the problem setting must be reconsidered, possibly by modifying the
r
ground structure.
r
If there is a feasible solution to P0 , let ρ∗ be the solution. Note that ρ∗ may still have
many small but nonzero folding angles. A sequence of optimization problems is used to set
these angles to zero while retaining foldability.
To find a feasible solution with a reduced number of “on” foldlines near ρ∗ , the foldline
with the smallest folding angle

ρq = mini∈U r ρ∗
0 i

(6.11)

r
r
is eliminated. A new set of free folding angles U1 = U0 \ q is defined and ρq is set to zero.
A new optimization problem is solved, using now the reduced set of design variables. The
process is repeated until no feasible solution is found. The pth step in this process solves
109

r
Optimization Problem Pp :

Find

r
ρi for i ∈ Up that

Minimizes f (ρ)
Subject to gk = 0 for k = 1, 2, · · · M

(6.12)

hi = 0 for i = 1, 2, · · · Meq
r
−π ≤ ρi ≤ π for ∀i ∈ Up
∗p−1
r
ρi − ρi
≤ δ ∀i ∈ Up
r
r
r
where Up = Up−1 \ qp−1 is the set of indices of the free folding angles in Pp . Note that an
∗p−1
additional constraint ρi − ρi
≤ δ has been added to restrict the search within a box
of width δ centered at the previous solution ρ∗p−1 .
r
If there is no feasible solution to Problem Pp , ρ∗p−1 is a solution to the problem for the
given choice of fixed fold angle ρr . The algorithm is summarized in the flowchart in Fig.6.6.

6.4.4

Sensitivity analysis

To facilitate numerical implementation, gradients based on analytical expressions may be
computed and provided to the optimizer. The foldability constraint depends explicitly on
folding angle ρ, and the gradients of Gk can be computed as
ij

k
k
dGk
k − I dF = Fk − I Rρk Rθk Rρk Rθk · · · dRρi Rθk · · · Rρk Rθk
= F
n n
1 1 2 2
i
dρk
dρk
dρk
i
i
i
110

(6.13)

Start
r

Solve optimization problem P 0
Is there a
solution?

N

End
No solution

Y
p← 1
Identify smallest foldline q
Fix ρq=0
r

Solve optimization problem P p
Is there a
solution?

N

End

ρ∗p-1 is the solution

Y
p ← p+1
Figure 6.6: Flowchart of foldline elimination algorithm

111

if the ith crease is connected to the k th vertex and set to 0 otherwise. The constraints h
associated with the geometric properties are implicit functions of X, which depends on ρ.
First gradients of vertex coordinates with respect to ρ are computed, then the chain rule is
applied to find the gradients of the constraints as

∂h dX
dh
=
dρi
∂X dρi

(6.14)

The gradients of f are



 0

; i ∈ Up
/ r
df
=
2
dρi 
 2Ca−2 x a2 − ρ2 e− ρi /a
r
; i ∈ Up
i

6.5

(6.15)

Numerical examples

Two origami tessellations are designed to test the method. The ground structure shown in
Fig.6.2(a) with N = 88 design variables is used. Note that the periodicity condition requires
that folding angles corresponding to foldlines on the right boundary are the same as those
on the left boundary, and those on the top boundary are the same as those on the bottom
boundary. The number of design variables in this particular ground structure equals the total
number of segments minus four. The optimization problem is solved using the interior-point
approach described in [97]-[99], implemented in the optimization tool fmincon in MATLAB
[100].
112

Xb2

X n4

Xa1

Xn3

Xa2

Xn1

Xb1

X n2

Figure 6.7: Labeled vertices

6.5.1

Example 1

The first example seeks a folding pattern that changes the size of the tessellation tile as it
folds, while maintaining the aspect ratio. This property is expressed using

J=

T
T
1
1
Xb − Xb
X a1 − X a2
Xb − Xb
X a1 − X a2 −
1
2
1
2
2
2

(6.16)

with target value J ∗ = 0. The locations of midpoints Xa1 , Xa2 , Xb and Xb in the flat
1
2
configuration are shown in Fig. 6.7.
A random initial guess ρ0 with a small magnitude (< π/3000) is used. The problem is
solved multiple times, each time using a different fixed folding angle ρr = ρ0 , choosing r
r
from the segments in Fig.6.8. Due to symmetry, only the 13 segments highlighted in Fig.6.8
need to be investigated.
113

Figure 6.8: Candidates of fixed folding angles

(a) Crease Pattern

(b) Tessellation

Figure 6.9: Design 1, Example 1

(a) Crease Pattern

(b) Tessellation

Figure 6.10: Design 2, Example 1
114

Figures 6.9 and 6.10 show crease patterns and folded tessellations corresponding to two
of the 13 solutions found, which satisfy constraints g and h within the specified tolerance
ǫ = 0.05. Parameters a = π/4, δ = π/12 and ρ0 = π/4 are used. The dotted and solid lines
r
in the crease patterns represent mountain and valley folds, respectively.

6.5.2

Example 2

For the second example, additional constraints are used to ensure that the projected shape
of the folded sheet follows a flat, square surface. The first geometric property J1 is the same
as before (as in Eq.6.16). To keep the tile shape square,

J2 =

1
2

X a1 − X a2

T

2
(6.17)

Xb − Xb
1
2

∗
should be maintained at J2 = 0. To design a tessellation that follows a flat surface (non
conforming) as it is folded and unfolded, it is also necessary that

1
J3 =
2

Xb1 − Xa1 × Xb2 − Xa1
Xb1 − Xa1 × Xb2 − Xa1

(Xa2 − Xa1 ) × Xb2 − Xa1
−
(Xa2 − Xa1 ) × Xb2 − Xa1

2
(6.18)

J4 =

2
1
Xa2 − Xa1 − Xn2 + Xn1
2

(6.19)

J5 =

2
1
Xa2 − Xa1 − Xn3 + Xn4
2

(6.20)

115

(a) Crease Pattern

(b) Tessellation

Figure 6.11: Design 1, Example 2

(a) Crease Pattern

(b) Tessellation

Figure 6.12: Design 2, example 2

J6 =

2
1
Xb2 − Xb1 − Xn4 + Xn1
2

(6.21)

J7 =

2
1
Xb2 − Xb1 − Xn3 + Xn2
2

(6.22)

∗
∗
∗
be maintained at J3 = J4 = · · · J7 = 0.
Figures 6.11 and 6.12 show crease patterns and folded tessellations corresponding to two
of the solutions found, which satisfy the constraints within the specified tolerance ǫ = 0.05
as before.
116

6.6

Application: design of polarization insensitive tunable FSS

In this section, the usability of the origami design method is demonstrated in an electromagnetic problem. Frequency selective surfaces (FSSs) that can be tuned in frequency based
on the folding and unfolding motion are designed. The origami design method is used to
design folding patterns than can be used to construct polarization insensitive tunable FSSs.
Designs considered are thin substrates decorated with conducting elements and folded into
tessellations. The representative unit of periodicity is designed, using appropriate boundary
conditions.
The goal of the optimization problem is to find a folding pattern that changes the size of
the tessellation tile as it folds, while maintaining the aspect ratio, which correspond to the
geometric constraints used in Sec.6.5.2. Two sets of rectangular strips are placed orthogonal
to each other on a square sheet of substrate to make a unit cell of an FSS that is symmetric to
two orthogonal polarizations as shown in Fig.6.13(a). The the mechanism of resonance tuning
of the design here is similar to the one used in [59], in which cross-shaped conducting elements
are printed on each facet of a substrate folded into Miura-ori, and the change in the length of
periodicity as a result of folding shifts the resonance frequency. The main difference is that
the design used here is perfectly symmetric for the x- and y-directions in its flat configuration,
while in [59] the layout of decoration is asymmetric even in its flat configuration, and thus
the FSS’s performance dramatically changes for different polarizations. By requiring that
the substrate folds in a way that changes the lengths of periodicity in the x- and y-directions
by the same amount, the resonance behavior is expected to be nearly symmetric for incidence
117

l

l

w

w
y

y

x

x
(a) Two orthogonal sets of coupled strips

(b) One set of coupled strips

Figure 6.13: Conducting decoration for FSS

with the electric field polarized along both x- and y-directions.
The full evaluation is done by computing the transmission and reflection coefficients
using HFSS. In simulations, the tile size in the flat configuration is set to 200mm, and
the conducting elements illustrated in Fig.6.13 with w = 20mm and l = 120mm and the
conductivity of copper are used. Since the strips are larger than the areas of the triangles
and quadrilaterals enclosed by lines in the ground structure, the strips will have crease lines
through them in the FSS’s folded states. The substrate is assumed to be very thin, and
the properties of air are used in simulations. A plane wave incidence with the propagation
directed normally to the plane of the sheet is used.
Figure 6.14 shows the transmission coefficients |S21 | for the polarization-insensitive FSS
that satisfies the geometric constraints. The solid and dotted curves show |S21 | for the
incidence with the electric field polarized along the x-direction in their flat and folded configurations, and circles and squares show the data for the orthogonal polarization where the
electric field is in the y-direction. One can observe that the performance is nearly the same
for the two orthogonal polarizations. In both polarizations, the resonance shift of 8% - 9%
118

0

|S21| (dB)

−5
−10
−15

Flat, E along x
Folded, E along x
Flat, E along y
Folded, E along y

−20
−25
1.2

1.4
1.6
Frequency (GHz)

1.8

Figure 6.14: Transmission coefficient
is observed, covering nearly the same range of frequencies.
The decoration is designed such that current is induced on the coupled strips along the
direction of the electric field with the first resonance occurring near λ/2 = l, where λ is
the vacuum wavelength of the incidence. Strip width w, separation of the elements and the
folding pattern will determine the exact location of the resonance. As the sheet is folded,
size of the periodic cell decreases by 9%, changing the coupling of elements, resulting in a
shift in resonance frequency. A large shift in resonance can be achieved when strips are small
compared to the size of the cell, because the change in the cell size is more significant relative
to the size of the strips. However, use of very small strips will degrade the strength of the
resonance. The decoration was designed in an ad hoc manner based on those considerations,
such that a strong resonance with a large shift is achieved.
An alternative design is to use two layers of sheets decorated with one set of strips shown
in Fig.6.13(b), rotated 90° from each other. The transmission coefficients are computed for
such a design and shown in Fig.6.15. It is seen that the performance is again nearly the
119

0

|S21| (dB)

−5
−10
−15
−20
−25
1.2

Flat, E along x
Folded, E along x
Flat, E along y
Folded, E along y
1.4
1.6
Frequency (GHz)

1.8

Figure 6.15: Transmission coefficient for the two-layer design
same for the two orthogonal polarizations.

6.7

Conclusions

A method to design origami patterns based on topology optimization is introduced. Folding
patterns with desired geometric properties are found by assigning presence and fold types
to crease lines in a “ground structure”, using a topology optimization method with folding
angles as design variables. Usability of the proposed origami design method is demonstrated
in design of tunable FSSs that are insensitive to polarization of the incident electromagnetic
wave, which can be tuned in frequency by 8%.
One of the challenges in implementing the origami design method was to eliminate small
folds. The technique used in the proposed method involves an objective function designed to
aid the selection of a small fold to eliminate. It is required that the optimization problem is
solved many times, once after one fold is eliminated. The selection of a proper combination
120

of unnecessary folds may be made more efficiently through a formulation involving integer
programming.
Another potential issue is that for a given set of geometric constraints, an inappropriate
choice of the ground structure or the value of the fixed fold may lead to a problem with no
feasible solution. In the future work, an investigation of the relation between the existence
of a solution for a certain set of geometric constraints and the ground structure, as well as
the value of the fixed folding angle is needed.
The usability of the origami design method is demonstrated in an electromagnetic problem, where tunable FSSs based on the folding and unfolding motion are designed. The
origami design method is used to design folding patterns than can be used to construct
polarization insensitive tunable FSSs.

121

Chapter 7
Conclusions
In the first part of this dissertation, methods of topology optimization were used to design
periodic structures, composed of dielectric materials, with prescribed transmission/reflection
characteristics. Two dimensional and three dimensional design problems were solved.
A formulation using the level set function to describe the material distribution was developed for design of 2D structures, in an attempt to develop a formulation that results in
solutions with no mixture (gray) regions. In computation, smoothing of the level set function was necessary for a gradient-based optimization algorithm, and in the end solutions
still suffered from having mixture regions at interfaces. It was also observed that several
parameters affecting the analysis accuracy and optimization algorithm need to be adjusted
in an ad hoc manner but with great care to avoid numerical instability. Another formulation
based on the density approach was investigated. A rigorous mesh truncation method using
a plane wave expansion was applied at the input and output boundaries in this formulation.
The solution strategy proposed was stable and robust, and the formulation presented could
accommodate complex material distributions near boundaries by including higher modes
122

into consideration.
A similar formulation based on the density approach was used to design a 3D structure
with periodicity in the x− and y−directions. A periodic structure with the material varying
in the xy-plane and uniform in the z−direction was considered. A structure that minimizes
transmission for an incidence of a prescribed frequency was found. The effective density was
projected, using a smooth Heaviside function, to a near binary density in the analysis. It
was found that use of the parameter, that adjusts the sharpness of the projection function,
needs great care, in order to achieve a fast and stable convergence.
In the second part of this dissertation, a new design method based on topology optimization and the concept of origami was introduced for design of electromagnetic devices
involving transformations of complex 3D geometries for tuning of their working frequencies.
First the feasibility of tuning the working frequency of electromagnetic devices via a
folding and unfolding motion was demonstrated using two sample designs. An FSS based
on a Miura-ori structure with cross-shaped metallic prints was shown to shift its resonance
frequency as the folding configuration was altered, mainly due to the change in the periodic
cell size. A metamaterial based on a corrugated dielectric sheet decorated with SRRs was
also shown to achieve tunability in resonance frequency as an unfolding motion was applied.
The two split-rings are arranged vertically on the same plane in the initial configuration, and
as an unfolding motion is applied, the planes hosting the rings become separated, shifting
the resonance frequency rapidly. A simple analysis revealed that the main cause of the rapid
shift was the decrease in capacitance due to the air gap that opened between the split-rings.
A method to design origami patterns based on topology optimization was then introduced. Folding patterns with desired geometric properties were found by assigning presence
123

and fold types to crease lines in a “ground structure”, using a topology optimization method
with folding angles as design variables. Usability of the proposed origami design method was
demonstrated in design of tunable FSSs that are insensitive to polarization of the incident
electromagnetic wave. An appropriate folding pattern was designed by finding a pattern
that changes the size of the tessellation tile as the sheet is folded, while maintaining the aspect ratio. Copper strips were placed symmetrically on the sheet to construct a polarization
insensitive tunable FSSs. In the optimization formulation, finding desirable folding patterns
without having many small folds was found to be challenging. The proposed method requires that the optimization problem is solved many times in order to avoid many small
folds. Development of a more efficient approach for removing unnecessary small folds is
needed. Another potential issue with the proposed formulation is that for a given set of
geometric constraints, an inappropriate choice of the ground structure or the value of the
fixed fold may lead to a problem with no feasible solution. In the future work, investigation
of the relation between the existence of a solution for a certain set of geometric constraints
and the ground structure, as well as the value of the fixed folding angle is needed.

124

APPENDICES

125

Appendix A
Finite Element Equations for 2D
Periodic Structures
The governing equation is the scalar Helmholtz equation in Eq.2.9:

2
∇2 Ez + k0 ǫr Ez = 0

(A.1)

The field is assumed to take the form

+∞
Ez (x, y) = exp[jχ0 (y − yin )] exp(−jα0 x) +

m=−∞

rm exp[−jχm (y − yin )]ψm

(A.2)

at the input boundary Γin and
+∞
Ez (x, y) =
m=−∞

tm exp[−jχm (y − yin )]ψm

(A.3)

at the output boundary Γout , where rm and tm are modal reflection and transmission
126

coefficients, respectively, and ψm = exp(−jαm x). The propagation vector components αm
and χm are computed as

αm = α0 +

2πm
d

and





2
2
k0 − αm
χm =

 −j α2 − k 2
m
0

;
;

2
2
k0 − αm ≥ 0
2
2
k0 − αm < 0,

(A.4)
m∈Z

Consider fields Ez and w that satisfy Bloch-Floquet boundary conditions such that

Ez ∈ S1 := {Ez ∈ H 1 (Ω) : Ez (x + nd, y) = Ez (x, y)e−jnα0 d },

n∈Z

w ∈ S2 := {w ∈ H 1 (Ω) : w(x + nd, y) = w(x, y)e+jnα0 d },

n∈Z

(A.5)

The weak form can be obtained by multiplying both sides by a test function w and integrating
as

Ω
=⇒

Ω

2
w∇2 Ez + k0 ǫr wEz dΩ = 0
2
−∇w∇Ez + k0 ǫr wEz dΩ +

The boundary terms are
127

Γ

∇Ez • nwdΓ = 0 ∀w ∈ S2
ˆ

(A.6)

Γ

∇Ez • nwdΓ =
ˆ
+

Γin
Γa

∇Ez • (−ˆ)wdΓ +
y

∇Ez • (−ˆ)wdΓ +
x

Γout
Γb

∇Ez • y wdΓ
ˆ

∇Ez • xwdΓ
ˆ

(A.7)

and the terms from the periodic boundaries (Γa and Γb ) vanish as

−

Γa

∇Ez • xwdΓ +
ˆ

=−

Γa

∇Ez • xwdΓ +
ˆ

Γb
Γa

∇Ez • xwdΓ
ˆ
exp(+jα0 d) exp(−jα0 d)∇Ez • xwdΓ = 0
ˆ

(A.8)

since w|Γ = exp(+jα0 d) w|Γa and ∇Ez |Γ = exp(−jα0 d) ∇Ez |Γa .
b
b
The term from the input boundary (Γin ) is

Γin

∇Ez • (−ˆ)wdΓ = −
y





+∞

∂Ez
jχ ψ − j
χm rm ψm  wdΓ
wdΓ =
0 0
∂dy
Γin
Γin
m=−∞
(A.9)

using Eq.A.2. The term from the output boundary (Γout ) is

Γout

∇Ez • y wdΓ =
ˆ



+∞



∂Ez
−j
χm tm ψm  wdΓ
wdΓ =
∂dy
Γout
Γout
m=−∞

using Eq.A.3. Substituting Eq.A.8-A.10 into Eq.A.6, we obtain
128

(A.10)

Ω

2
(−∇Ez (x, y) • ∇w(x, y) + k0 ǫr Ez (x, y)w(x, y))dΩ
+∞

−j

+∞
χ m rm

m=−∞

= −jχ0

Γin

Γin

ψm (x)w(x, yin )dΓ − j

χm t m
m=−∞

Γout

ψm (x)w(x, yout )dΓ

ψ0 (x)w(x, yin )dΓ

(A.11)

Equations for the reflection and transmission coefficients rm and tm are obtained by
∗
multiplying both sides of Eq.A.2 and A.3 by ψm and integrating as
˜

Γin

∗
Ez ψm dΓ =
˜

Γin

+∞

∗
ψ0 ψm dΓ +
˜

rm

Γin

∗
ψm ψm dΓ and
˜

m=−∞
+∞
∗
∗ dΓ =
ψm ψm dΓ
Ez ψm
tm
˜
˜
Γout
Γout
m=−∞

(A.12)

Using the orthogonality relationship

Γ

∗
ψm ψm dΓ =
˜

we obtain

rm d −

Γin



 0 ; m=m
˜

(A.13)


 d ; m=m
˜

∗
Ez (x, yin )ψm (x)dΓ = −δ0m d and

tm d −

Γout

∗
Ez (x, yout )ψm (x)dΓ = 0

129

(A.14)

Equations A.11 and A.14 will be discretized, following a finite element formulation with
linear rectangular elements, to find an approximate electric field component Ez and reflection
and transmission coefficients (rm and tm ). Using linear shape functions { φm } Nn , where
m=1
Nn is the number of degrees of freedom, Eqs.A.11 and A.14 are approximated by finite
element equations

Kz = F

(A.15)

where z = [Ez , r, t]T , and F = fE , fr , 0 T and




 KEE KEr KEt 




K= K
dI
0 

 rE


KtE
0
dI

(KEE )sn =

Ω

2
˜
˜
(−(∇φn ) · (∇φs ) + k0 ǫr φn φs )dΩ

(KEr )nm = −jχm
(KEt )nm = −jχm

Γout

(KrE )mn = −
(KtE )mn = −

Γin

˜
ψm (x)φn (x, yout )dΓ

Γin

Γout

(fE )n = −jχ0

˜
ψm (x)φn (x, yin )dΓ

∗
φn (x, yin )ψm (x)dΓ

∗
φn (x, yout )ψm (x)dΓ

Γin

˜
ψ0 (x)φn (x, yin )dΓ
(fr )m = −δ0m d

130

where s, n ∈ {−Nn , −Nn +1, · · · , 0, · · · , Nn −1, Nn } and m ∈ {−Nm , −Nm +1, · · · , 0, · · · , Nm −
˜
1, Nm }. Here φ and φ are bilinear shape functions associated with Ez and w, and the corresponding DOFs on the periodic boundaries at x = 0 and x = d are enforced to satisfy
conditions Ez ∈ S1 and w ∈ S2 .

131

Appendix B
Finite Element Equations for 3D
Periodic Structures
The governing equation is the vector wave equation in Eq.2.7:

∇×

jσ
1
2
∇ × E − k0 ǫr −
µr
ωǫ0

E=0

(B.1)

The electric field is assumed to take the form

E (x, y, z) = Ein + Eref
p
ˆ
= exp −j α0 x + β0 y + γ0 z − zin
+∞
+∞
+
rmn exp jγmn z − zin ψmn
m=−∞ n=−∞
at the input boundary Γin and
132

(B.2)

E (x, y, z) = Etr
+∞

+∞

=
m=−∞ n=−∞

tmn exp[−jγmn (z − zout )]ψmn

(B.3)

at the output boundary Γout , where ψmn = exp [−j (αm x + βn y)] and rmn and tmn are
modal reflection and transmission coefficients that can be computed as:

rmn =

1
∗
∗
E(x, y, z)ψmn − pψ00 ψmn dΓ
ˆ
Tx Ty Γ
in

(B.4)

1
∗
E(x, y, z)ψmn dΓ
Tx Ty Γout

(B.5)

and

tmn =

The propagation constant for each mode γmn can be found by




where

2
2
2
2
2
2
k0 − αm − βn ; αm + βn ≤ k0
γmn =

 −j α2 + β 2 − k 2 ; α2 + β 2 > k 2
m
n
m
n
0
0

2πm
Tx
2πn
βn = β0 −
Ty

(B.6)

αm = α0 −

Consider fields E and V that satisfy Bloch-Floquet boundary conditions such that
133

(B.7)

E ∈ S1 := {E ∈ H 1 (Ω) : E(x + mTx , y + nTy )
= E(x, y, z) exp(−jmα0 Tx − jnβ0 Ty )},

n, m ∈ Z

V ∈ S2 := {V ∈ H 1 (Ω) : V(x + mTx , y + nTy )
= V(x, y, z) exp(+jmα0 Tx + jnβ0 Ty )},

n, m ∈ Z

(B.8)

The weak form is obtained through multiplying both sides by a test function V and
integrating as

Ω

V • (∇ ×

1
jσ
2
∇ × E)dΩ −
k0 ǫr −
µr
ωǫ0
Ω

V • EdΩ = 0 ∀V ∈ S2

(B.9)

This can be organized as

jσ
1
2
V • E dΩ
(∇ × V) • (∇ × E) − k0 ǫr −
ωǫ0
Ω µr
1
V • (ˆ × ∇ × E) dΓ = 0 ∀V ∈ S2
n
+
Γ µr

(B.10)

using the first vector Green’s theorem (see pp. 711 [39]):

Ω

a • (∇ × u∇ × b)dΩ =

Ω

u(∇ × a) • (∇ × b)dΩ −
134

Γ

uˆ • (a × ∇ × b)dΓ
n

(B.11)

and the vector identity
a • (b × c) = −c • (b × a)

(B.12)

The boundary terms in Eq.B.10 are

1
n • (V × ∇ × E) dΓ
ˆ
Γ µr
1
1
=
(−ˆ) • (V × ∇ × E) dΓ +
x
x • (V × ∇ × E) dΓ
ˆ
Γlef t µr
Γright µr
1
1
+
(−ˆ) • (V × ∇ × E) dΓ +
y
y • (V × ∇ × E) dΓ
ˆ
Γf ront µr
Γback µr
1
1
(−ˆ) • (V × ∇ × E) dΓ +
z
z • (V × ∇ × E) dΓ
ˆ
+
Γout µr
Γin µr

(B.13)

and the terms from the periodic boundaries (Γa and Γb ) vanish, similarly to the 1-D periodicity case. The weak form will then become

jσ
1
2
V • E dΩ
(∇ × V) • (∇ × E) − k0 ǫr −
ωǫ0
Ω µr
1
1
V • (ˆ × ∇ × E) dΓ −
z
V • (ˆ × ∇ × E) dΓ = 0 ∀V ∈ S2
z
+
Γin µr
Γout µr
135

(B.14)

The term from the input boundary (Γin ) is
1
V • (−ˆ × ∇ × E) dΓ
z
Γin µr


+∞
+∞
1
=−
p
rmn ψmn  dΓ
V • z × ∇ × exp[−j (α0 x + β0 y + γ0 )]ˆ +
ˆ
Γin µr
m=−∞ n=−∞
+∞
+∞
1
1
∗
Vψmn dΓ •
z×∇×
ˆ
Eψmn dΓ exp[+jγmn (z − zin )]
=
µr
Tx Ty
Γin
m=−∞ n=−∞ Γin
+
−

Γin

1
V • z × ∇ × exp[−jγ00 (z − zin )]ψ00 p dΓ
ˆ
ˆ
µr

Γin

1
V • z × ∇ × exp[+jγ00 (z − zin )]ψ00 p dΓ
ˆ
ˆ
µr

(B.15)

using Eq.B.2 and B.4.
The term from the output boundary (Γout ) is
1
V • (ˆ × ∇ × E) dΓ
z
Γout µr


+∞
+∞
1
=
ˆ
V • z × ∇ ×
tmn ψmn  dΓ
Γout µr
m=−∞ n=−∞
+∞
+∞
1
1
∗
Eψmn dΓ exp[−jγmn (z − zout )]
Vψmn dΓ •
z×∇×
ˆ
=
µr
Tx Ty
Γout
m=−∞ n=−∞ Γout
(B.16)

using Eq.B.3 and B.5.
Equation B.14 will be discretized, following a finite element formulation using brick edge
elements as in [39], where the electric field components within each element e can be expressed
as
136

4

6

8

3

12

11
10

9

z

2

5

y

7

1

x
Figure B.1: Brick element

e
Ex =

4
i=1

e
φe E i ,
i

e
Ey =

8
i=5

e
φe E i ,
i

e
Ez =

12
i=9

e
φe E i
i

(B.17)

and φi is a shape functions that corresponds to the tangential field component along edge i,
as numbered in Fig.B.1.
137

1
φe = e e
1 l l
yz
1
φe = e e
2 l l
yz
1
φe = e e
3 l l
yz
1
φe = e e
4 l l
yz
1
φe = e e
5 l l
zx
1
φe = e e
6 l l
zx
1
φe = e e
7 l l
zx
1
φe = e e
8 l l
zx
1
φe = e e
9 l l
xy
1
φe = e e
10 l l
xy
1
φe = e e
11 l l
xy
1
φe = e e
12 l l
xy

e
e
e l
e + ly − y
zc + z − z
yc
2
2
e
e
e l
e + ly + y
zc + z − z
−yc
2
2
e
e
e l
e + ly − y
−zc + z + z
yc
2
2
e
e
e + ly + y
e l
−yc
−zc + z + z
2
2
e
le
e l
zc + z − z
xe + x − x
c
2
2
e
le
e l
−zc + z + z
xe + x − x
c
2
2
e
le
e l
zc + z − z
−xe + x + x
c
2
2
e
le
e l
−zc + z + z
−xe + x + x
c
2
2
e
le
e ly
xe + x − x
yc + − y
c
2
2
le
−xe + x + x
c
2

e
yc +

e
ly

−y
2
e
le
e ly
xe + x − x
−yc + + y
c
2
2
e
le
e ly
−xe + x + x
−yc + + y
c
2
2
(B.18)

The finite element equations corresponding to Eq.B.14 are

KE = F
138

(B.19)

in which the element matrices are

e
Kij =

jσ ˜
1
2
˜
φi φj
∇ × φi • ∇ × φj − k0 ǫr −
ωǫ0
Ωe µr
+∞
+∞
1 ˜
1
−
φ ψ dΓ •
φ ψ ∗ dΓ
z×∇×
ˆ
e µe i mn
e j mn
Tx Ty
Γin
m=−∞ n=−∞ Γin r
+∞
+∞
1
1 ˜
∗
φj ψmn dΓ
φ ψ dΓ •
z×∇×
ˆ
+
e
e µe i mn
Tx Ty
Γout
m=−∞ n=−∞ Γout r

(B.20)

and

Fie =

1 ˜e
φ • z × ∇ × exp[−jγ00(z−z ) ]ψ00 p dΓ
ˆ
ˆ
e µe i
in
Γin r
1 ˜e
−
φ • z × ∇ × exp[+jγ00(z−z ) ]ψ00 p dΓ
ˆ
ˆ
e µe i
in
Γin r

(B.21)

˜
where φ and φ are shape functions associated with the electric field E and the test function
V, respectively, and the corresponding DOFs on the periodic boundaries are enforced to
satisfy the conditions E ∈ S1 and V ∈ S2 .

139

BIBLIOGRAPHY

140

BIBLIOGRAPHY
[1] Bendse M. P. and Kikuchi N. Generating optimal topologies in structural design using a
homogenization method. Computer Methods In Applied Mechanics And Engineering, Vol.
71, pp. 197-224, 1988.
[2] Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. International Journal of Solids and Structures, Vol. 31 No. 17, pp. 2313-2339,
1994.
[3] Bendse M. P. and Sigmund O. Topology optimization : Theory, Methods, and Applications. Springer: Berlin, Heidelberg, 2003.
[4] Sigmund O. and Torquato S. Design of materials with extreme thermal expansion using
a three-phase topology optimization method. Journal of the Mechanics and Physics of
Solids, Vol. 45, No. 6, pp. 1037-1067, 1997.
[5] Olesen L.H., Okkels F. and Bruus H. Topology optimization of Navier-Stokes flow in
microfluidics. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2004), Neittaanmaki P., Rossi T., Korotov S., Onate E., Periaux
J. and Knorzer D. (eds.) Jyvaskyla, 24-28 July 2004.
[6] Gersborg-Hansen A., Sigmund O. and Haber R.B. Topology optimization of channel flow
problems. Structural and Multidisciplinary Optimization, Vol. 30, No. 3, pp. 181-192, 2005.
[7] Bruns T.E. Topology optimization of convection-dominated, steady-state heat transfer
problems. International Journal of Heat and Mass Transfer, Vol. 50, No. 15-16, pp. 28592873, 2007.
[8] Wadbro E., Berggren M. Topology optimization of an acoustic horn. Computer Methods
in Applied Mechanics and Engineering, Vol. 196, No. 1-3, pp. 420-436, 2006.
[9] Dhring M.B., Jensen J.S., Sigmund O. Acoustic design by topology optimization. Journal
of Sound and Vibration, Vol. 317, No.3-5, pp. 557-575, 2008.
141

[10] Chung Y.S. and Cheon C. Optimal shape design of dielectric structure using FDTD and
topology optimization. Microwave Symposium Digest, 2001 IEEE MTT-S International,
Vol. 3, pp. 2063 - 2066, 2001.
[11] Kiziltas G., Psychoudakis D., Volakis J. L. and Kikuchi N. Topology design optimization
of dielectric substrates for bandwidth improvement of a patch antenna. IEEE Transactions
On Antennas And Propagation, Vol. 51, No. 10, pp. 2732-2743, 2003.
[12] Jensen J.S. and Sigmund O. Systematic design of photonic crystal structures using
topology optimization: Low-loss waveguide bends. Applied Physics Letters, Vol. 84, No.
12, pp. 2022-2024, 2004.
[13] Borel P.I., Harpoth A., Frandsen L.H., Kristensen M., Shi P., Jensen J.S. and Sigmund
O. Topology optimization and fabrication of photonic crystal structures. Optics Express,
Vol. 12, No. 9, pp. 1996-2001, 2004.
[14] Jensen J.S. and Sigmund O. Topology optimization of photonic crystal structures: a
high-bandwidth low-loss T-junction waveguide. Journal of the Optical Society of America
B, Vol. 22, No. 6, pp. 1191-1198, 2005.
[15] Tsuji Y., Hirayama K., Nomura T., Sato K. and Nishiwaki S. Design of optical circuit
devices based on topology optimization. IEEE Photonics Technology Letters, Vol. 18, No.
5-8, pp. 850-852, 2006.
[16] Nomura T., Sato K., Taguchi K., Kashiwa T. and Nishiwaki S. Structural topology
optimization for the design of broadband dielectric resonator antennas using the finite
difference time domain technique. International Journal for Numerical Methods in Engineering, Vol. 71, pp. 1261-1296, 2007.
[17] Sigmund O. and Hougaard K. Geometric properties of optimal photonic crystals. Physical Review Letters, Vol. 100, No. 15, pp. 153904, 2008.
[18] Riishede J. and Sigmund O. Inverse design of dispersion compensating optical fiber
using topology optimization. Journal of the Optical Society of America B, Vol. 25, No. 1,
pp. 88-97, 2008.
[19] Diaz A.R. and Sigmund O. A topology optimization method for design of negative
permeability metamaterials. Structural and Multidisciplinary Optimization, Vol. 41, No.2,
pp. 163-177, 2009.
142

[20] Nomura T., Nishiwaki S., Sato K. and Hirayama K. Topology optimization for the
design of periodic microstructures composed of electromagnetic materials. Finite Elements
in Analysis and Design, Vol. 45, No. 3, pp. 210-226, 2009.
[21] Fuchi K., Diaz A.R., Yamada T. and Nishiwaki S. A level set-based topology optimization method for power flow control in electromagnetics. 8th Congress on Structural and
Multidisciplinary Optimization, Lisbon, Portugal, 2009.
[22] Dhring M.B. and Sigmund O. and Feurer T. Design of photonic bandgap fibers by
topology optimization. Journal of the Optical Society of America B, Vol. 27, No. 1, pp.
51-58, 2010.
[23] Matzen R., Jensen J.S. and Sigmund O. Topology optimization for transient response
of photonic crystal structures. Journal of the Optical Society of America B, Vol. 27, No.
10, pp. 2040-2050, 2010.
[24] Andkjr J., Nishiwaki S., Nomura T. and Sigmund O. Topology optimization of grating
couplers for the efficient excitation of surface plasmons. Journal of the Optical Society of
America B, Vol. 27, No. 9, pp.1828-1832, 2010.
[25] Aage N., Mortensen N.A. and Sigmund O. Topology optimization of metallic devices
for microwave applications. International Journal for Numerical Methods in Engineering,
Vol. 83, No. 2, pp. 228-248, 2010.
[26] Erentok A. and Sigmund O. Topology optimization of sub-wavelength antennas. IEEE
Transactions on Antennas and Propagation, Vol. 59, No. 1, pp. 58-69, 2011.
[27] Pendry J.B., Holden A.J., Robbins D.J. and Stewart W.J. Magnetism from conductors and enhanced nonlinear phenomena. IEEE Transactions on Microwave Theory and
Techniques, Vol. 47, No. 11, pp. 2075-2084, 1999.
[28] Smith D.R., Padilla W.J., Vier D.C.,Nemat-Nasser S.C. and Schultz S. Composite
medium with simultaneously negative permeability and permittivity. Physical Review Letters, Vol. 84, No. 18, pp. 4184-4187, 2000.
[29] Shelby R.A., Smith D.R and Schultz S. Experimental Verification of a Negative Index
of Refraction. Science, Vol. 292, No. 5514, pp. 77-79, 2001.
[30] Haupt R. L. An introduction to genetic algorithms for electromagnetics. IEEE Antennas
Propagation Magazine, Vol. 37, pp. 7-15, 1995.
143

[31] Ouedraogo R. O., Rothwell E. J. , Diaz A. R., Chen S. Y., Temme A. and Fuchi K. In
situ optimization of metamaterial-inspired loop antennas. IEEE Antennas and Wireless
Propagation Letters, Vol. 9, pp. 75-78, 2010.
[32] Kern D.J., Werner D.H., Monorchio A., Lanuzza L. and Wilhelm M.J. The design
synthesis of multiband artificial magnetic conductors using high impedance frequency
selective surfaces. IEEE Transactions on Antennas and Propagation, Vol. 53, No. 1, pp.
8-17, 2005.
[33] Chen C.H., Chen P.Y, Wang H., Tsai J.H. and Ni W.X. Synthesis design of artificial
magnetic metamaterials using a genetic algorithm. Optics Express, Vol. 16, No. 17, pp.
12806-12818, 2008.
[34] Bayraktar Z. Gregory M.D., Wang X. and Werner, D.H. Matched impedance thin planar
composite magneto-dielectric metasurfaces. IEEE Transactions on Antennas and Propagation, Vol. 60, No. 4, pp. 1910-1920, 2012.
[35] Chung Y.S., Ryu J., Cheon C., Park I.H. and Hahn S.Y. Optimal design method for
microwave device using time domain method and design sensitivity analysis - part I: FETD
case. IEEE Transactions on Magnetics, Vol. 37, No. 5, pp. 3289-3293, 2001.
[36] Chung Y.S., Cheon C., Park I.H. and Hahn S.Y. Optimal design method for microwave
device using time domain method and design sensitivity analysis - part II: FDTD case.
IEEE Transactions on Magnetics, Vol. 37, No. 5, pp. 3255-3259, 2001.
[37] Nikolova N.K., Tam H.W. and Bakr M.H. Sensitivity analysis with the FDTD method
on structured grids. IEEE Transactions on Microwave Theory and Techniques, Vol. 52,
No. 4, pp. 1207-1216, 2004.
[38] Silvester P.P. and Ferrari R.L. Finite elements for electrical engineers, Third edition.
Cambridge University Press, Cambridge, New York, 1996.
[39] Jin J. The finite element method in electromagnetics, Second edition. IEEE Press, John
Wiley & Sons, Inc., Cambridge, New York, 2002.
[40] Jin J.M. and Riley D.J. Finite element analysis of antennas and arrays. IEEE Press,
John Wiley & Sons, Inc., Hoboken, 2009.
[41] Bao G. Finite-element approximation of time-harmonic waves in periodic structures.
SIAM Journal on Numerical Analysis, Vol. 32, No. 4, pp. 1155-1169, 1995.
144

[42] Dossou K., Byrne M.A. and Botten L.C. Finite element computation of grating scattering matrices and application to photonic crystal band calculations. Journal of Computational Physics, Vol. 219, No. 1, pp. 120-143, 2006.
[43] Allaire G., Jouve F. and Toader A.M.. Structural optimization using sensitivity analysis
and a level-set method. Journal of Computational Physics, Vol. 194, No. 1, pp. 363-393,
2004.
[44] Wang M.Y., Wang X. and Guo D. A level set method for structural topology optimization. Computer Methods in Applied Mechanics and Engineering, Vol. 192, No. 1-2, pp.
227-246, 2003.
[45] Bourdin B. and Chambolle A. The phase-field method in optimal design. IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials: Status
and perspectives, M. P. Bendse, N. Olhoff, and O. Sigmund (Eds.), pp. 207-215, Springer,
2006.
[46] Comsol:www.comsol.com.
[47] Sigmund O. and Hougaard K. Geometric properties of optimal photonic crystals. Physical Review Letters, Vol. 100, No. 15, pp. 3904, 2008.
[48] Bendse M.P. Optimal shape design as a material distribution problem. Structural Optimization, Vol. 1, No. 4, pp. 193-202, 1990.
[49] Rozvany G.I.N., Zhou M. and Birker T. Generalized shape optimization without homogenization. Structural Optimization, Vol. 4, No. 3-4, pp. 250-252, 1992.
[50] Mlejnik H.P. and Schirrmacher R. An engineers approach optimal material distribution
and shape finding. Computer Methods in Applied Mechanics and Engineering, Vol. 106,
No. 1-2, pp. 1-26, 1993.
[51] Svanberg K. The method of moving asymptotes - a new method for structural optimization. International Journal for Numerical Methods in Engineering, Vol. 24, No. 2, pp.
359-373, 1987.
[52] Bruns T.E. and Tortorelli D.A. Topology optimization of non-linear elastic structures
and compliant mechanisms. Computer Methods in Applied Mechanics and Engineering,
Vol. 190, No. 26-27, pp. 3443-3459, 2001.
145

[53] Bourdin B. Filters in topology optimization. International Journal for Numerical Methods in Engineering, Vol. 50, No. 9, pp. 2143-2158, 2001.
[54] Sigmund O. Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization, Vol. 33, No. 4-5, pp. 401-424, 2007.
[55] Wang F., Lazarov B.S. and Sigmund O. On projection methods, convergence and robust
formulations in topology optimization. Structural and Multidisciplinary Optimization, Vol.
43, No. 6, pp. 767-784, 2011.
[56] Demaine ED and O’Rourke J. Geometric folding algorithms: linkages, origami, polyhedra. Cambridge University Press: New York, 2007.
[57] Miura K. Method of Packaging and Deployment of Large Membranes in Space. The
institute of space and astronautical science report, Vol.618, pp. 1-9, 1985.
[58] Freeland R.E., Bilyeu G.D., Veal G.R. and Mikulas M.M. Inflatable deployable space
structures technology summary. IAF 49th Congress, Melbourne, Australia, Sep. 1998.
[59] Fuchi K., Tang J., Crowgey B., Diaz A.R., Rothwell E.J. and Ouedraogo R.O. Origami
tunable frequency selective surfaces. IEEE Antennas and Wireless Propagation Letters,
Vol. 11, pp. 473-475, 2012.
[60] Fuchi K., Diaz A.R., Rothwell E.J., Ouedraogo R.O. and Tang J. An origami tunable
metamaterial. Journal of Applied Physics, Vol. 111, No. 084905, 2012.
[61] Mias C. Varactor-tunable frequency selective surface with resistive-lumped-element biasing grids. IEEE Microwave and Wireless Components Letters, Vol.15, pp. 570572, 2005.
[62] Schoenlinner B., Abbaspour-Tamijani A., Kempel L.C. and Rebeiz G.M. Switchable
low-loss RF MEMS Ka-band frequency-selective surface. IEEE Transactions on Microwave
Theory and Techniques, Vol. 52, No. 11, pp. 2474-2481, 2004.
[63] Li M., Yu B. and Behdad N. Liquid Tunable Frequency Selective Surfaces. IEEE Microwave and Wireless Component Letters, Vol. 20, No. 8, pp. 423-425, 2010.
[64] Lima A.C.D., Parker E.A. and Langley R.J. Tunable frequency selective surface using
liquid substrates. Electronics Letters, Vol. 30, No. 4, pp. 281-282, 1994.
[65] Bossard J.A., Liang X., Li L., Yun S., Werner D.H., Weiner B., Mayer T.S., Cristman P.E., Diaz A. and Khoo I.C. Tunable frequency selective surfaces and negative-zero146

positive index metamaterials based on liquid crystals. IEEE Transactions on Antennas
and Propagation, Vol. 56, No. 5, pp. 1308-1320, 2008.
[66] Chang T.K., Langley R.J. and Parker E.A. Frequency selective surfaces on biased ferrite
substrates. Electronics Letters, Vol. 30, No. 15, pp.1193-1194, 1994.
[67] Bossard J.A., Werner D.H. and Mayer T.S. A novel design methodology for reconfigurable frequency selective surfaces using genetic algorithms. IEEE Transactions on Antennas and Propagation, Vol. 53, No.4, pp.1390-1400, 2005.
[68] HFSS: www.ansys.com.
[69] Gil I., Garcia-Garcia J., Bonache J., Martin F., Sorolla M. and Marques R. Varactorloaded split ring resonators for tunable notch filters at microwave frequencies. Electronics
Letters, Vol. 40, No. 21, pp. 1347-1348, 2004.
[70] Gil M., Bonache J., Selga J., Garcia-Garcia J. and Mart´ F. Broadband resonant-type
ın
metamaterial transmission lines. IEEE Microwave and Wireless Components Letters, Vol.
17, No. 2, pp. 97-99, 2007.
[71] Aydina K. and Ozbay E. Capacitor loaded split ring resonators as tunable metamaterial
components. Journal of Applied Physics, Vol. 101, No. 2, pp. 024911-1-024911-5, 2007.
[72] Li H., Zhang Y. and He L. Tunable filters based on the varactor-loaded split-ring resonant structure coupled to the microstrip line. Microwave and Millimeter Wave Technology
International Conference, Vol. 3, pp. 1580-1582, 2008.
[73] Hand T.H. and Cummer S.A. Frequency tunable electromagnetic metamaterial using
ferroelectric loaded split rings. Journal of Applied Physics, Vol. 103, No. 6, pp. 066105 066105-3, 2008.
[74] Hand T.H. and Cummer S. Controllable magnetic metamaterial using digitally addressable split-ring resonators. IEEE Antennas and Wireless Propagation Letters, Vol. 8,
pp.262-265, 2009.
[75] Ouedraogo R.O. Topology optimization of metamaterials and their applications to RF
component design, (Doctoral dissertation). Michigan State University, East Lansing, MI,
2011.
147

[76] Khoo I.C., Werner D.H., Liang X., Diaz A. and Weiner B. Nanosphere dispersed liquid
crystals for tunable negative-zero-positive index of refraction in the optical and terahertz
regimes. Optics Letters, Vol. 31, pp. 25922594, 2006.
[77] Werner D.H., Kwon D.H., Khoo I.C., Kildeshev A.V. and Shalaev V.M. Liquid crystalclad near-infrared metamaterials with tunable negative -zero-positive refractive indices.
Optics Express, Vol. 15, pp. 33423347, 2007.
[78] Zhao Q., Kang L., Du B., Li B. and Zhou J. Electrically tunable negative permeability
metamaterials based on nematic liquid crystals. Applied Physics Letters, Vol. 90, pp.
011112-1 - 011112-3, 2007.
[79] Khoo I.C., Diaz A., Liou J., Stinger M.V., Huang J. and Ma Y. Liquid crystals tunable
optical metamaterials. IEEE Journal of Selected Topics in Quantum Electronics, Vol. 16,
No. 2, pp. 410-417, 2010.
[80] Marqu´s R., Mart´ F. and Sorolla M. Metamaterials with negative parameters: theory,
e
ın
design, and microwave applications, Wiley & Sons, Inc. Hoboken, New Jersey, 2007.
[81] Harrington R.F. Field Computation by Moment Methods, IEEE Press, Piscataway, NJ,
2000.
[82] Randall C.L., Gultepe E. and Gracias D. Self-folding devices and materials for biomedical applications. Trends in Biotechnology, Vol. 30,No. 3, pp. 138-146, 2011.
[83] Gracias D.H., Tien J., Breen T.L., Hsu C. and Whitesides G.M. Forming electrical
networks in three dimensions by self-assembly. Science, Vol. 289, No. 5482, pp. 1170-1172,
2000.
[84] Jamal M., Bassik N., Cho J.H., Randall C.L. and Gracias D.H. Directed growth of
fibroblasts into three dimensional micropatterned geometries via self-assembling scaffolds.
Biomaterials, Vol. 31, No. 7, pp. 1683-1690, 2010.
[85] belcastro S. and Hull T. Modelling the folding of paper into three dimensions using
affine transformations. Linear Algebra and its Applications, Vol. 348, pp. 273-282, 2002.
[86] Lang R.J. The tree method of origami design. In The Second International Meeting of
Origami Science and Scientific Origami, Miura K (ed). Seian University of Art of Design:
Otsu, Japan, pp. 72-82, 1994.
148

[87] Lang R.J. Threemaker 4.0: A Program for Origami Design, 1998.
http://www.langorigami.com/science/computational/treemaker/TreeMkr40.pdf.
[88] Tachi T. Simulation of Rigid Origami. In The Fourth International Conference on
Origami in Science, Mathematics, and Education, Lang R (ed). Pasadena, USA, pp. 175187, 2009.
[89] Tachi T. 3D origami design based on tucking molecule. In The Fourth International
Conference on Origami in Science, Mathematics, and Education, Lang R (ed). Pasadena,
USA, pp. 259-272, 2009.
[90] Tachi T. Origamizing polyhedral surfaces. IEEE Transactions on Visualization and
Computer Graphics, Vol. 16, No. 2, pp. 298-311, 2010.
[91] Tachi T. and Demaine E.D. Degenerative Coordinates in 22.5° Grid System, In The Fifth
International Meeting of Origami Science, Mathematics and Education, Wang-Iverson P,
Lang RJ and Yim M (eds). Singapore Management University: Singapore, pp. 489-497,
2010.
[92] Dorn W.S., Gomory R.E. and Greenberg H.J. Automatic design of optimal structures.
Journal de Mecanique, Vol 3, No 1, pp. 25-52, 1964.
[93] Bendse M.P., Bental A. and Zowe J. Optimization methods for truss geometry and
topology design. Structural and Multidisciplinary Optimization, Vol. 7, No. 3, pp. 141159, 1994.
[94] Gjerde E. Origami Tessellations: Awe-inspiring geometric designs. A K Peters, Massachusetts, 2008.
[95] Maekawa J. Genuine origami: 43 mathematically-based models, from simple to complex.
Japan Publications Trading, Tokyo, 2008.
[96] Elsayed E.A. and Basily B. A continuous folding process for sheet materials. International Journal of Materials and Product Technology, Vol. 21, No. 1/2/3, pp. 217-238,
2004.
[97] Byrd R.H., Hribar M.E. and Nocedal J. An interior point algorithm for large-scale
nonlinear programming. SIAM Journal on Optimization, Vol. 9, No. 4, pp. 877-900, 1999.
149

[98] Byrd R.H., Gilbert J.C. and Nocedal J. A trust region method based on interior point
techniques for nonlinear programming. Mathematical Programming, Vol. 89, No. 1, pp.
149-185, 2000.
[99] Waltz R.A., Morales J.L., Nocedal J. and Orban D. An interior algorithm for nonlinear
optimization that combines line search and trust region steps. Mathematical Programming,
Vol. 107, No. 3, pp. 391-408, 2006.
[100] MATLAB:www.mathworks.com.

150