вход по аккаунту


EM modeling and simulation of microwave electronic components and devices with multi-scale and multi-physics effects

код для вставкиСкачать
Presented in Partial Fulfillment of the Requirements for
the Degree Doctor of Philosophy in
the Graduate School of The Ohio State University
Jue Wang, B.Eng
Graduate Program in
Electrical and Computer Engineering
The Ohio State University
Dissertation Committee:
Prof. Jin-Fa Lee, Advisor
Prof. Kubilay Sertel
Prof. Fernando L. Teixeira
ProQuest Number: 10001908
All rights reserved
The quality of this reproduction is dependent upon the quality of the copy submitted.
In the unlikely event that the author did not send a complete manuscript
and there are missing pages, these will be noted. Also, if material had to be removed,
a note will indicate the deletion.
ProQuest 10001908
Published by ProQuest LLC (2016). Copyright of the Dissertation is held by the Author.
All rights reserved.
This work is protected against unauthorized copying under Title 17, United States Code
Microform Edition © ProQuest LLC.
ProQuest LLC.
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346
c Copyright by
Jue Wang
This work investigates various numerical methods for modeling and analyzing
microwave components and electronic devices with multi-scale structures and multiphysics effects. In the first part of the work, a universal matrice approach is developed that can handle continuous changing material properties across the element and
curved boundary. This method is implemented in an Interior Penalty based domain
decomposition solver. Several interesting numerical examples including the modeling
of the Luneburg lens and conformal perfectly matched layer validate this method.
In the second part, a full-wave homogenization process of extracting the equivalent
permeability tensor of ferromagnetic nanowire array is presented. By solving for the
small signal model of the Landau-Lifshits equation, the macroscopic material property
can be obtained for the heterogeneous structure. After that, the utilization of ferromagnetic nanowire array into microwave components and devices is studied, showing
interesting properties such as double-band working frequencies and self-bias capability. Coming to the third part of this dissertation, time domain method for transient
linear/non-linear effects in high-frequency circuit systems is explored. To be specific,
a discontinuous Galerkin time-domain method integrated with SPICE circuit solver
and IBIS models, respectively, is developed. The key technique is on the interface
the coupling of EM solver and circuit solver. This work provides a self-consistent
integration so to physically model the whole system. Since the coupling process is
local for DG method, only modest resources are needed. Overall, this dissertation
evaluates and develops several numerical methods with the determination to provide
a platform for multi-scale and multi-physics simulations for EM analysis of modern
radio frequency systems with linear/nonlinear and passive/active factors.
This work is dedicated to the ones that I love.
Over the past years I’m blessed to have received numerous support and encouragement from a great number of individuals.
The deepest gratitude goes to my advisor, Professor Jin-Fa Lee, whose guidance,
support and patience have made this work possible. His rigorous and passionate
attitude towards academia shapes mine. His gracious share of the philosphy in life as
a friend will constantly inspire me throughout the remaining of my life.
I also would like to thank Prof. Fernando Teixeira and Prof. Kubilay Sertel
for being my examination committee members, for their constructive comments and
I am particularly grateful to Dr. Zhen Peng. His wide knowledge has given me
significant help during every critical stage of my research work. What is more, his
unwavering encouragement and trust has made me believe in myself.
During my study, I have also received tremendous help from my colleagues and
friends in the Electroscience Lab, many thanks for their companionship during one
of the most important and unforgetable stages in my life: Dr. Xiaochuan Wang, Dr.
Yang Shao, Dr. Davood Ansari, Dr. Jiangong Wei, Dr. Stylianos Dosopoulos, Dr.
Matt Stephanson, Feiran Lei, Yuanhong Zhao, Caleb Waltz, Kheng-Hwee Lim, Joshua
Mahaffey, Wan-Chun Liao, Cagdas Gunes, Dr. Dongwei Li, Jiaqing Lu, Shaoxin Peng,
Xuezhe Tian and many more others.
In the summer of 2010, thanks to my advisor, I got the opportunity of taking an
internship with ANSYS, Inc. It was another fruitful journey. First of all, I want to
thank Dr. Steve Bardi for being my mentor. Immersed into the first-class company
industrial atmosphere, I received a lot of unselfish help and friendship from different
people, to name a few here: Dr. Kezhong Zhao, Dr. Guanghua Peng and Dr. SeungCheol Lee.
Lastly, I want to thank my family for their unconditional love and trust in me
always. Even across the ocean, thinking about them, I always find inner peace.
Especially my mom, the greatest woman I have ever known and a constantly inspiring
role model all through my life.
The names of the people who have contributed to my progress all those years is
impossible to be completely enumerated here, I sincerely appologize for those who
haven been dropped in this acknowlegement.
July, 2009 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.Eng. Electrical Engineering,
University of Science and Technology
of China, Hefei, China
September, 2009-present . . . . . . . . . . . . . . . . . . . . Graduate Research Associate,
ElectroScience Laboratory,
The Ohio State University.
Research Publications
D. Ansari, J. Wang, Z. Peng, J.-F Lee “A universal array approach for finite elements
with continuously inhomogeneous material properties and curved boundaries”. IEEE
Transactions on Antennas and Propagation, vol. 60, no. 10, pp. 4745-4756, 2012.
Z. Peng, J. Wang, F. R. Lei, J.-F Lee “New computational strategies for electromagnetic modeling of multi-scale heterogeneous composites”. Antennas and Propagation
(EUCAP), Proceedings of the 5th European Conference, pp. 3226-3229, 2011.
J. Wang, Z. Peng, J.-F Lee “ Ferromagnetic nano-wires: homogenization and applications ”. Radio Science Bulletin, no. 349 , 2014.
Fields of Study
Major Field: Electrical and Computer Engineering
Major Field: Electrical Engineering
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Vita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . .
Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . .
Literature Review on Numerical Methods in Multi-scale and Multi-physics
Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Overview . . . . . . . . . . . . .
FDTD and FETD methods . . .
DGTD methods . . . . . . . . .
Homogenization of Heterogeneous
Multi-physics coupling . . . . . .
A Universal Array Approach for Finite Elements with Continuously Material Properties and Curved Boundaries . . . . . . . . . . . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Ferromagnetic Nano-wires: Homogenization and Applications . . . . . .
Introduction . . . . . . . . . . . . . .
The Proposed Homogenization Method
Numerical Results . . . . . . . . . . .
Conclusion . . . . . . . . . . . . . . .
Self-Consistent Integration of Discontinuous Galerkin Time Domain Method
and Circuit Simulator for Microwave Components and Devices with Circuit Network Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Methodology . . . . . . . . . .
Notation . . . . . . . . . . . . .
Universal Arrays . . . . . . . .
CIMPT Approach for the PML
Numerical Results . . . . . . .
Conclusion . . . . . . . . . . .
DGTD integrated with SPICE solver . . . . . . . . . . . . . . . . . 69
DGTD integrated with IBIS model . . . . . . . . . . . . . . . . . . 101
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
List of Tables
3.1 Conformal-DDM Related Notation . . . . . . . . . . . . . . . . . . .
3.2 Solution statistics for the waveguide problem. The “2” in the “Field
Basis“ row indicates the use of a second order curl conforming basis, i.e.
the curl of the basis is complete to the 2nd order. In the reference result
in [55], (m, n, p) represents the orders of the field basis along (x̂, ŷ, ẑ)
directions in hexahedral elements (see Fig. (3.5)). By material (Mat.)
basis we refer to the scalar basis used for the interpolation of the MPT. 32
4.1 External DC Magnetic field vs. Normalized Magnetization . . . . . .
5.1 Geometrical description of the pulse generator . . . . . . . . . . . . .
5.2 Step-recovery diode description . . . . . . . . . . . . . . . . . . . . .
List of Figures
3.1 (a) The 3D reference element. (b) An actual curved element.
. . . .
3.2 (a) Edge representation of the curved element. (b) Face representation
of the curved element. Only two faces are plotted here. . . . . . . .
3.3 The cubic Lagrangian interpolation of the curved tetrahedron.
. . .
3.4 Visualization of a subset of the curved tetrahedral elements used in a
mesh for the Luneburg lens. The surface of the lens is plotted with
some transparency such that inside-sphere edges are also visible. . . .
3.5 (a) Problem geometry. (b) |S11 | versus frequency. . . . . . . . . . . .
3.6 (a)Electrical field intensity. (b) Conformal-DDM’s Partition of the
Mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 Schematic diagram showing the configuration of the quarter-Luneburg
lens scattering problem. Curved elements were applied to the surface
of the lens only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8 Schematic diagram showing the configuration of the R = 6λ0 Luneburg
lens antenna problem. Curved elements were applied to the surface of
the lens only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 A three-slice field plot of the E · ŷ field component in the continuous
index Luneburg lens antenna with R = 6λ0 and f = 10 GHz. The
shaded part of the plot belongs to the dielectric lens. . . . . . . . . .
3.10 The radiation pattern (dB) of Luneburg lens antennas of various radius
R as function of θ at ϕ = 0. . . . . . . . . . . . . . . . . . . . . . . .
3.11 Schematic diagram showing the configuration of the PML-lens antenna
problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 (a) Snapshots of the magnitude of the electric field distribution on yzplane. (b) Snapshots of the magnitude of the electric field distribution
on zx-plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.13 The radiation pattern (dB) of the Luneburg lens antenna with R = λ0
at ϕ = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 Sketch of FMNW substrate structure. . . . . . . . . . . . . . . . . . .
4.2 Schematic hysteresis loop. . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Flow-chart of FMNW based metamaterial homogenization process and
application enumeration. . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Comparison of effective relative permeability tensor of CoFe nanowire
membrane obtained by different methods. . . . . . . . . . . . . . . .
4.5 Comparison of effective relative permeability tensor of NiFe nanowire
membrane obtained by different methods. . . . . . . . . . . . . . . .
4.6 Microstrip line sketch. . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7 Major hysteresis loop (solid blue) and minor hysteresis curves (dashed
red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.8 Left: contour plot of the S21 parameter from simulated result. Right
(a) and (b): reference experiment and model results. . . . . . . . . .
4.9 Electric field distributions along the microstrip line with FMNW subext
strate and biased with the same external bias field, Hdc
= 0.5KOe, at
two different resonance frequencies. . . . . . . . . . . . . . . . . . . .
4.10 FMNW isolator structure. . . . . . . . . . . . . . . . . . . . . . . . .
4.11 Comparison of relative characteristic permeability. Solid lines stand for
analytical model in the reference while marks of triangles and squares
denote the results of the proposed numerical homogenization method.
4.12 S12 parameter comparison. . . . . . . . . . . . . . . . . . . . . . . . .
4.13 S21 parameter comparison. . . . . . . . . . . . . . . . . . . . . . . . .
4.14 Electric field distributions along the microstrip line of the isolator at
two different working frequencies. . . . . . . . . . . . . . . . . . . . .
4.15 Detailed description of a dual-band self-biased circulator. . . . . . . .
4.16 Full-wave simulation result of the conceptual microstrip circulator–S
parameter plotting. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.17 Full-wave simulation result of the conceptual microstrip circulator–
Electric field plotting. . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Sketch of EM-Circuit System . . . . . . . . . . . . . . . . . . . . . .
5.2 (a) EM coupled to SPICE (b) SPICE coupled to EM . . . . . . . . .
5.3 Geometrical illustration of the impedance surface. . . . . . . . . . . .
5.4 SPICE netlist example. . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Flowchart of DGTD integrated with SPICE . . . . . . . . . . . . . .
5.6 Integration by Gaussian quadrature rule . . . . . . . . . . . . . . . .
5.7 Microstrip line with circuit network. . . . . . . . . . . . . . . . . . . .
5.8 CST multi-line port experiment. . . . . . . . . . . . . . . . . . . . . .
5.9 Result comparison at near end of microstrip. . . . . . . . . . . . . . .
5.10 Result comparison at far end of microstrip. . . . . . . . . . . . . . . .
5.11 Microstrip line with nonlinear circuit network. . . . . . . . . . . . . .
5.12 Comparisons of numerical results with different schemes. . . . . . . .
5.13 Schematic structure of SRD pulse generator . . . . . . . . . . . . . . .
5.14 Circuit structure of SRD pulse generator. . . . . . . . . . . . . . . . .
5.15 Result comparison of currents at the input side of the SRD pulse generator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.16 Result comparison of voltages at the output side of the SRD pulse
generator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.17 Input buffer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
5.18 Output buffer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.19 Illustration of the introduced multiplier switching behavors. . . . . . 106
5.20 Equivalent analog behavior model. . . . . . . . . . . . . . . . . . . . 106
5.21 Microstrip line with IBIS models. . . . . . . . . . . . . . . . . . . . . 109
5.22 Package top view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.23 Result comparison of voltages at the driver side of the link. . . . . . . 110
5.24 Result comparison of voltages at the receiver side of the link. . . . . . 111
5.25 Geometry Description of the simple board. (a): top view, (b): side view.111
5.26 Result comparison of voltages at the microcontroller (port1). . . . . . 112
5.27 Result comparison of voltages at the memory chip (port2). . . . . . . 113
Chapter 1: Background and Motivation
Computational electromagnetics (CEM) is the scientific subject of modeling the
interaction of electromagnetic fields with physical objects and the enviroment using
computational approximations to Maxwell’s equations. Today’s technologies, particularly in high-frequency electronics, communication, biomedical and geophysical areas
among others, are developing so fast and changing people’s lives every day. They all
depend on advanced electromagnetic equipment. Most of those modern real-world
electromagnetic problems, due to the irregular or complicated geometries involved,
do not posess closed-form analytical solutions. With computational tools, it allows
engineers to refine and validate designs at a stage where the cost of making changes is
minimal. Thererfore, virtually every industry now incorporates computer-based engineering simulation early in the development process. Indeed, engineering simulation
software has revolutionized electronics design, and as a consequence with the aid of
the fast development in high-performance computational platform, large-scale CEM
problems that require supercomputers and/or clusters have become possible, solving
problems with increasing complexity.
How to accurately and efficiently model electromagnetic fields in complex designs
and systems has gained a lot of importance during the past decades. Various branches
of CEM have emerged since then. From the big picture, CEM methods [1] can be
categorized into partial differential equation (PDE) based and integral equation (IE)
based methods. Those methods are further distinguished as frequency and time
domain methods.
Among numerous methods, finite element methods (FEM) [2] have been developed intensely and have achieved huge success in many different areas. They solve
the Maxwell’s equations starting by discretizing the computational domain into finite
elements, then define a set of linear basis functions on each element, and approximate the solution in a function space with finite dimension. FEM formulates the
problem from governing PDE, which gives great flexibility in dealing with complex
materials. What’s more, the grids can be irregular, it provides a way of more accurate
and flexible geometry modeling. Since the introduction of tangential vector FEMs
(TVFEM) [3] [4] [5] [6], FEM has become a standard method in a variety of applications. Compared to conventional nodal basis functions for discretization [7], TVFEM
eliminates spurious modes by enforcing only the tangential components of the vector
field to be continuous across the element boundaries. However, edge-elements only
provide the lowest-order of accurary in numerical computations; in order to achieve
higher accurary, refining mesh is one way, using higher order basis functions is another way. In [8], the construction of higher-order tangential vector finite elements is
presented. The author in [9] proposes a hierarchal vector basis functions of arbiturary
order with sepate representation of the gradient and rotational parts of the vector
Recent development of domain decomposition methods (DDM) based on FEM
has proved DDM to be the most common paradigm for large-scale simulations. In
DDMs ( [10] [11]), a large problem is reduced to multiple smaller problems, each
of which is easier to solve computationally and most or all of which can be solved
independenly and concurrently. Iterative coupling occurs among the boundaries of
each smaller domains. On the topics of how to reduce the number of interations
and make sure the result converge to a physical solution lies much of its theoretical
interest. With the concept of divide and conquer, many large scale real-life problems
with complicated features can be solved with modest resources [12] [13] [14].
The FEM that solves the time-harmonic Maxwell’s equations is the frequencydomain FEM. In recent years, numerical methods to solve Maxwell’s equations directly in the time domain have also emerged with great promises. Although the timemarching process is usually time-consuming, compared to frequency-domain methods,
time-domain methods are well suited for the broadband phenomena. In a single simulation they provide a solution in a wide frequency range. Furthermore, time-domain
methods are capable to model transient and nonlinear circuits and devices. Finite
difference time-domain (FDTD) method [15] [16] has gained the popularity due to
its simplity. However, it is at the same time well-known that FDTD methods suffer
from inefficiency in modeling complex geometry since they require structured grids
as well as encounter dispersion problems. Although the issues are addressed in later
development [17] [18], the simplicity is very difficult to maintain, or in other situations
the accuracy is sacrificed. Therefore, more sophisticated and versatile time-domain
techniques are needed. Time-domain FEMs (TDFEM) are encouraged under this
scenario. Early development of TDFEM can be found in [19] [20]. The main research
interest in TDFEM since then addresses stability analysis [21], diagonalization of
the system matrices using orthogonal basis functions [22], hybrid method [23], among
many others. From theories to applications, many contributions have been conducted
in TDFEM techniques.
Recently, the discontinuous Galerkin time-domain (DGTD) method has gained
wide attention. DG methods are a class of FEMs that utilize piecewise continuous
basis and testing functions on the discretized elements, and the continuity condition on the boundary of adjacent elements are weakly enforced [24] [25]. they are
well-suited on unstructured meshes for high order approximation with the freedom
of choosing the shape of the grids and the order of basis functions locally. When
applied in time domain, DG results in a global block diagonal mass matrix, therefore
high efficiency in parallelization can be easily achieved [26] [27]. Since information
exchange occurs on adjacent elements, when using an explicit time-marching scheme,
it allows for local time-stepping at a reduced cost; in this scenario, the time stepping
size for different groups of mesh with different densities can be chosen differently, instead of being bounded by the smallest mesh size globally [28]. Furthermore, DGTD
method can be hybridized with many other solvers with minor modification, such as
the integration of circuit solvers. With all those salient features, DGTD method is
a robust method capable of solving complex problems of full model radio frequency
Problem Statement
Although numerous CEM methods have achieved great performance in the past,
it is still difficult to achieve both certain level of accuracy and low memory cost at
fast calculation speed for conventional CEM methods when handling certain types of
real-life problems. One of the challenges comes from the heterogeneity of the geometry involved in many real-life situations. Without attempting to be exhaustive the
following applications are mentioned: the three-dimensional (3D) lithography technology, the effective properties of composite materials increasingly used in enigeering,
advances in interconnect technologies in today’s integrated circuit and package products. This kind of problem contains multi-scale features. Under such disparate scales
Maxwell’s equations show drastically different phenomena. That will aggravate the
diversity of the system eigen spectrum, causing an ill-conditioned system. Another
difficulty comes from the interactions of multi-physical building blocks in one problem.
Almost all practical engineering applications are multi-physics in nature, and various
physical phenomena usually interact with each other. For instance, consider high
power integrated circuits, packages and printed circuit boards, the resistivity of conducting metals increases with the increase of the surrounding temperature resulting
from Joule heating when electrical currents flowing through conductors. Therefore,
it is essential to account for both electrical and thermal effects and their couplings in
order to characterize accurately the performance. This requires the CEM engine to
be versatile enough so to integrate with other computational solvers easily.
While a lot of previous literature have addressed the multi-scale and multi-physics
problems, the scope of the research is delicated to the following main directions:
1. Development of a continuous material property tensor approach for higher order FEM analysis of problems with continuously varying material properties.
Many real-life applications require the modeling of continuous spatial changes
in material properties. Spatial fluctuations are usually modeled using piecewise
constant functions within each element. However, using an appropriate polynomial functions to replace the piecewise representation yields obvious benefits.
First of all, it results in a more accurate physical model for problems where
material properties follow continuous changes. Secondly, it allows for better
utilization of computational resources when combined with higher order curved
elements. This work introduces a universal array approach for handling isoparametric FE matrices for which material properties are allowed to vary as smooth
polynomials across individual elements. With this capability, for multi-physics
problems where dielectric properties change due to the change of other physical
phenomena, the re-meshing process of the problem domain can be avoided.
2. Numerical extraction of the equivalent permeability tensor of ferromagnetic
nanowire (FMNW) array and its applications in microwave/RF components
and devices. FMNW metamaterials exibit novel properties such as double ferromagnetic resonance (FMR), anisotropy control and self-bias. Those properties make FMNW array very appealing in low-profile microwave devices with
dual-band magnetic effect. In previous work [29] [30], the authors present a
Maxwell-Garnett like formulation to analytically extract the effective permeability tensor of FMNWs, while this research work uses a 3D full-wave numerical
process more generally to extract the macroscopic effective permeability tensor
of non-saturated/saturated arrays of axially magnetized FMNWs. After the homogenization process, the utilization of FMNWs based metamaterial structures
in microwave components/devices are studied through full-wave electromagnetic
3. Integration of discontinuous Galerkin time-domain method and circuit solver
for complex circuit and device simulation. Nowadays, with the ever increasing
operating speeds and integration densities of modern integrated circuits (IC),
in order to successfully predict the performance of the practical printed circuit
board (PCB) designs, accurate mixed electromagnetic and circuit co-simulation
have become imperative. This part of the research work uses a DGTD timedomain method integrated with circuit solver to aid the co-design process of the
components, such as microwave structures, and active/passive linear/nonlinear
circuits for PCB products. This method is appealing in several aspects which
will be addressed in more details later. In terms of the circuit solvers chosen to
be integrated with EM solver, there are two directions. One scenario is when
the circuit network structure is known, we use SPICE, the simulation program
with integrated circuit emphasis. SPICE is a well-known general-purpose, open
source analog electronic circuit simulator. It is a powerful program that is
used in integrated circuit and board-level design to check the integrity of circuit designs and to predict circuit behavior. The other scenario, which is more
common in real-life situations, is when the circuit structure is hidden for the
sake of proprietary information protection, IBIS models are used to integrate
with DGTD solver instead. IBIS (Input/Output buffer information specification) is an emerging standard for the behavioral modeling of IC input/output
characteristics. Either using SPICE or IBIS, the coupling between the EM and
circuit parts are addressed in a self-consistent manner with convergence looping
process. This coupling between DGTD and circuit part is very efficient since it
only involves local modification of the EM solver.
Dissertation Outline
The remainder of this dissertation is organized as following. Chapter 2 contains a
brief literature survey of existing approaches for various multi-scale and multi-physics
applications. Chapter 3 is devoted to the universal array approach for handling
geometries with continuously varying material property with curved boundary and
the applications. In Chapter 4, a full-wave homogenization process to extract the
equivalent permeability tensor of FMNW arrays is presented. With the macroscopic
description of the material, microwave devices that utilize FMNW arrays achiving
novel properties are simulated and numerical results are presented. Up till this point,
the scope of the work remains in frequency-domain methods, next in Chapter 5, timedomain method is introduced. Particularly, a self-consistent integration scheme of
DGTD method with circuit solvers for complex circuit and device transient simulation
is proposed. This part of work can further be divided into the integration of SPICE
solver and IBIS models, which will be detailed respectively, followed by numerical
examples. The last chapter closes the work by giving a summary and possible future
Chapter 2: Literature Review on Numerical Methods in
Multi-scale and Multi-physics Applications
Multi-scale and multi-physics modeling play a major role in many important problems today. This is particularly true in electromagnetic compatibility (EMC) problems. EMC is the science which studies the unintentional generation, propagation and
reception of electromagnetic energy with reference to the unwanted electromagnetic
interference (EMI) that such energy may induce. It widely exists in digital computing and wireless communication systems. Due to the ever-increasing complexity in
modern electronic systems, EMC has become one major issue in ICs design. How
to carry out the state of art simulation and modeling on IC, electronic package and
printed circuit board, thus is of great importance. In this kind of problems, it is very
common to have very different physical scales and the physics changes sigificantly over
these length-scales. Indeed, solving electromagnetics problem in this enrivonment is
a challenging task [31] [32].
To tackle this type of problems efficiently, it is essential to distinguish the features
of very detailed electrical structure from the rest of the computational domain, which
is on a different frequency scale. Generally speaking, two broad directions can be
taken. One is to locally refine the mesh in areas where fine features occur and field
varies fast. The other direction is to keep the mesh unchanged at a global resolution
taken into consideration. A local solution is calculated for the fields near the fine
feature and then coupled to the global solution elsewhere. In this situation, the
difficulties come from both obtaining local solutions and coupling process.
FDTD and FETD methods
Here a few representative works will be listed covering the topic. In [33], a subgridding scheme is depicted in which method the entire computational volume is divided
into a coarse grid with a large step size and a fine grid with a small step size only
around fine features. An interpolation in space and time is used to calculate the field
on the boundary of the two sets of grids. This method requires the boundary to
be placed in a homogeneous region, while in [34], the authors present a subgridding
method which allows the grid interface to be the boundary in inhomogeneous regions.
This further has provided more flexibility in multi-scale modeling. Although FDTD
method and its variations have achieved huge success during the early attempting,
they suffer from the inflexibility of modeling geometries since most of them require
regular grids. Herein, FE methods come to rescue. Being naturally constructed
for unstructured grids, FE is suited for numerical solution of PDEs with complex
geometry [35]. FEM shares some common features to FDTD. It also demands volume computational domain and decomposes the domain into finite elements, however,
compared to FDTD, it provides substantial flexibility in supporting various basis functions on unstructured element cell. In [36], an unconditionally stable FETD method
for active microwave circuit analysis with the truncation of perfectly matched layers
(PML) is presented. A time-domain finite element method for modeling complex
broad-band antennas is presented in [37].
DGTD methods
A more recent and viable technique is the DGTD method [24]. This method is able
to model complex geometries with complex materials (anisotropic and/or dispersive
materials) [38]. Since it is applied to time domain, either linear or non-linear materials can be treated. Futhermore, DG method is essentially a domain-decomposition
method [26], therefore, it is a highly-parallel algorithm and can be hybridized with
many other solvers. Since its introduction to solve for the time dependent Maxwell’s
equations, intensive research has been conducted. To list a few here: in [39], the authors have formed an energy conservative DGTD method based on central fluxes and
leap-frog marching in time. A DGTD method with upwind fluxes using a low storage
high-order Runge-Kutta marching scheme is developed in [25]. Compared to central
flux scheme, upwind flux, which closely resembles a Robin type transmission condition in the frequency domain decomposition method [40], results in a slight energy
lossy system, but with the optimal convergence rate. After the space discretization,
a semi-discrete system of ordinary differential equations (ODEs) is obtained. When
using an explicit time marching style, the system is conditionally stable; the time
step needs to be confined by the element size. In many real-life situations where
small features are present, the time step can be very small and it causes a very long
CPU time. To tackle this kind of issue, a local time stepping algorithm is studied
in [41], in which different time step sizes can be applied to different groups of elements
with different mesh densities. In this fashion, the computations can be largely speed
up. With all the benefits, it is noted that DGTD methods, with the introduction of
different sets of unknowns on the element interface, suffer from a larger number of
degrees of freedom compared to FETD methods.
Homogenization of Heterogeneous Composites
When studying the behavior of electromagnetic field in a material presenting heterogeneous microstructures, direct computations can often be very difficult due to
the fast oscillation of fields on the very fine and complex structures, and should be
avoided. Homogenization is a process in which the composite material is replaced
with an equivalent macrospoic, homogeneous material property. The advantages of
using a homogeneous system by averageing out the fine features are obvious: first, it
reduces considerably the cost of numerical simulation; secondly it provides the understanding of the macro dynamics of the considered problem. Analytical treatments of
such problems have been studied for many years and still are an active area [42] [43].
Numerical methods for homogenization have arised recently. One highlight is the heterogeneous multi-scale method (HMM) [44] [45]. A few other methods are listed here:
the authors in [46] present a homogenization method by calculating the reflection
and transmission coefficients on finite lengths of the electromagnetic metamaterials.
In [47], the equivalent material properties is extracted based on averaging the local
Multi-physics coupling
Apart from the multi-scale features that widely exist in modern EM RF problems,
multi-physics effects are often present as well. As a matter of fact, almost every problem is multi-physics in essence, since the world we live in inherently includes multiple
physical phenomena. With that pointed out, how to develop a multi-physics simulation platform where the different physics (that matter) couple to each other, based
on a common data model is of practical impact [48]. To list only a few efforts already
made here: in [49], the authors give an overview of reconfigurable antenna design and
optimizatin. A thermal-aware DC IR-drop co-analysis using non-conformal DDM is
performed in [50]. In [51], a SPICE-based multi-physics simulation technique for integrated MEMS is presented. [52] formulates a SPICE lumped circuit model within
the DGTD method, to solve for transient linear/nonlinear circuit systems, which will
be revisited and improved in later chapter of this dissertation.
Chapter 3: A Universal Array Approach for Finite Elements
with Continuously Material Properties and Curved
Many real life applications involve spatial changes in material property tensors
(MPTs). Clearly, integration of continuously inhomogeneous MPTs (CIMPT) inside
individual element adds to the flexibility and efficiency of FEMs. Curved FEs have
been extensively used to mitigate geometrical non-conformities associated to rectilinear approximation of curvilinear features while CIMPT are conventionally handled
by element-wise constant approximation. FE matrices are traditionally evaluated
through 1) numerical cubature or 2) universal matrices (UMs). In essence, both
methods rely on polynomial integration. Furthermore, complications associated with
evaluation of FE matrices on elements with curved boundaries or CIMPTs are identical in nature, i.e. integrals with non-constant Jacobian or MPT terms. Alternatively,
in this work, a generalized UM approach is proposed, which reduces the cost associated with evaluation of FE matrices with curvilinear features and CIMPTs. Motivated
by simulation of a non-graded Luneburg lens, the conventional element-wise constant
MPTs are replaced with polynomial representations yielding the following benefits:
a) more accurate physical models for problems with CIMPTs and curved features;
and, b) better utilization of computer resources when higher-order curved elements
replace lower-order elements of smaller physical dimensions.
Curved FEs have for long been used to mitigate geometrical-model non-conformities
encountered in computer aided modeling and simulation of real life problems. Isoparametric FEs provide a flexible tool for discretization of curved geometries by means
of polynomial or rational approximation of element coordinate-transform dependent
terms [53] [54]. In general, both numerical cubature and UM methods have been
deployed for evaluation of FE matrix entries. Nevertheless, both numerical cubature and UMs are theoretically based on polynomial approximation of the integrand,
which suggests that they should be equally applicable to various kinds of FE matrix
evaluations. Moreover, in conventional FEM codes, material properties, e.g. permittivity and permeability tensors in electromagnetics (EM), are treated as element-wise
constant functions while some recent works have examined the idea of having nonconstant material properties within individual elements [55] verifying that the techniques can be advantageous. Also in [56], Webb uses a polynomial representation of
the magnetic MPT for evaluation of FE matrices arising from a nonlinear magnetic
problem. While element-wise constant material properties can still be used in a wide
variety of problems, there are cases where material properties do change continuously
in space. The Luneburg lens, Maxwell’s fish-eye lens, semiconductor devices and media where material properties are affected or governed by diffusion of heat, moisture,
dopant etc. are practical examples of problems where the conventional element-wise
constant material approach may fail to provide an accurate model of the underlying
physics. In cases as such, material properties follow continuous spatial fluctuations
that can be modeled more accurately using polynomial approximation, e.g. interpolation.
In this work, we present a generalized UM, i.e. universal array (UA), approach for
evaluation of FE matrices for elements in which material properties are allowed to vary
continuously across individual elements and where element boundaries are allowed to
conform to given curvature. The method can be equally applicable to conventional
rectilinear elements with constant material properties, in which it simplifies to the
well known UM approach [57] with the caveat that our symmetric evaluation of matrix
entries further reduces the computational complexity.
As a test to the presented approach, we have implemented a conformal FEs domain decomposition method (FE-DDM) [40] solution of an ungraded Luneburg lens
operated in the X-Band. The lens problem is an exterior radiation/propagation problem with continuously changing material properties and curved boundaries associated
with the surface of the spherical Luneburg lens.
Accurate expression for conformal perfectly matched layer’s (PML) [58] permittivity and permeability tensors involves continuous functions of the spatial coordinates.
Hence, a PML with continuously varying material proprieties is examined for the
exterior truncation of the wave problem [59]. Obviously, the resulting FE matrices
are, again, evaluated using the proposed CIMPT method. The PML formulation is
briefly discussed followed by numerical results.
Systematically, the evaluation process of FE stiffness and mass matrices can be
cast into the following basic steps:
1. Assume a reference element with a fixed geometry.
2. Develop the required polynomial/vector-polynomial basis on the reference element.
3. Develop the required transformation rules for the basis, its derivations and
the Jacobian defined between the reference element and a presumed physical
4. Express the required integral-differential form on the physical element.
5. Pull back the integration operation onto the reference element, i.e. re-express
the integral over the reference element.
6. Identify the factors that are solely and entirely determined by the geometry
of the physical element from those solely defined (depend) on the reference
7. Being independent from the reference element coordinates, the metric factors
are pulled out of the integral (over the reference element) while the remaining
parts become independent of the metric properties of the physical element.
Hence, the integrals can be pre-calculated and stored in the so called UAs and
the evaluation of FE matrices turns into a sequence of multiply-add operations.
The following section explains how the abovementioned steps are adopted for
evaluation of FE matrices associated with time harmonic Maxwell’s equations. The
FE matrices are then used as an integral part of a conformal-DDM code which, for the
purpose of verification, is used for the solution of a few example problems including
an ungraded Luneburg lens.
Throughout this work, whenever summation bounds of are not explicitly given, it
is implied that they fall in the 0 ≤ index < 3 bounds. Furthermore the superscript
is used to denote matrix transposition. At the same time, column vectors and
matrices are denoted by [ai ] and [aij ] while their individual entries are denoted as
[aij ]ij and [ai ]i , respectively. The J and K symbols are used to denote forward and
backward Jacobian matrices between the reference and the actual physical elements;
and, z denotes a coordinate mapping defined between them. As visualized in Fig.
(3.1), the reference and the physical elements are respectively denoted by Kr and Kp
while ξ and x superscripts stand for the coordinate systems defined on Kr and Kp .
The × sign represents matrix multiplication while, in order to avoid confusion with ×,
∧ is used to denote exterior products in R3 . ∇x and ∇ξ denote the gradient operators
in their respective coordinates, i.e. x and ξ. Hence, for any sufficiently smooth real
function ϱ defined on either the reference or the physical element the column vector
form of the gradient is denoted as [∇]x ϱ , [ ∂x
] or [∇]ξ ϱ , [ ∂ξ
]. Furthermore, various
physical quantities and functional spaces are used in this work for which a brief listing
can be found in Table. (3.1).
Universal Arrays
This section presents a step-by-step realization of the concepts introduced.
The Reference Element Kr
As depicted in Fig. (3.1), the reference tetrahedron is chosen as Kr = {(ξ0 , ξ1 , ξ2 , ξ3 )|ξ0 +
ξ1 + ξ2 + ξ3 = 1, 0 ≤ ξ0 , ξ1 , ξ2 , ξ3 ≤ 1}. Note that the four coordinate variables on Kr
are dependent through ξ0 + ξ1 + ξ2 + ξ3 = 1. Hence, ξ3 is often omitted and replaced
by 1 − (ξ0 + ξ1 + ξ2 ). Monomial and polynomial integrations on Kr can be handled
Given the coordinate mapping z : (ξ0 , ξ1 , ξ2 ) → (x0 , x1 , x2 ) between Kr and a
physical element of interest Kp , the Jacobian matrices J and K can be defined as in
Eq. (3.1) in which ◦ denotes function composition. It is thus obvious that J and K
are the inverse of each other and that det J ̸= 0 and det K ̸= 0 are guaranteed by
] [
] [
∂ xi◦z(ξ)
∂ ξi◦z−1 (x)
↔ K=
Straightforward application of the chain rule reveals that the following relations
hold between column vector representation of the gradients defined over Kr and Kp .
[∇]x ϱ = KT × [∇ξ ]ϱ
[∇]ξ ϱ = JT × [∇x ]ϱ
Figure 3.1: (a) The 3D reference element. (b) An actual curved element.
The H(curl, Kr ) Basis and Transformations
Construction of the H(curl, Kr ) conforming basis (in either of its hierarchical or
spectral forms) has been extensively discussed in the literature [60] [61] [62] [63].
Regardless of the choice of the actual basis, any vector basis function β⃗uξ over the
reference element can be expressed as the sum of three scalar polynomials {βuξ |0 ≤
i < 3} each multiplied by a direction factor {∇ξ ξi |0 ≤ i < 3} [61]. Expressed
explicitly, we have
β⃗uξ =
∇ξ ξj
where the u subscript signifies a particular basis function.
From Eq. (3.4), it is obvious that a basis function β⃗uξ on Kr can be represented
by a column vector representation [βu,i
] with respect to the {∇ξ ξi } direction basis. It
is not difficult to observe that replacing the ∇ξ term in Eq. (3.4) with ∇x yields the
H(curl, Kp ) conforming basis on the physical element Kp . In other words, since ξ can
be written as a function of x (and vice versa), basis functions derived as Eq. (3.5)
satisfy the necessary tangential continuity conditions for a H(curl, Kr ) conforming
basis on Kp [64].
β⃗ux =
∇x ξj
From Eq. (3.5), one writes:
∑ [ ∂ξj ]
[ x]
→ βu,i
= KT × βu,j
Hence, with respect to the {∇x xi } direction basis, the column vector representa[ x]
tion of β⃗ux , i.e. βu,i
, is obtained as KT × [βu,j
The Mass Matrix
Using the electrical field formulation, the entries of the FE mass matrix T for an
electromagnetic radiation/propagation problem can be formulated as
Tuv =
· (ε̄¯ ·
β⃗vx )dKp
x T ¯
] × ε̄ ×[βv,j
where [βu,i
] and ε̄¯ respectively denote the column vector representation of β⃗ux (dis-
cussed in 3.4.2) and the permittivity tensor.
Using Eq. (3.6) and pulling the integration back onto Kr , Eq. (3.7) can be written
Tuv =
ξ T
] ×K× ε̄¯×KT [βv,j
] det J dKr
Next, if we define [ςij ] as Eq. (3.9), T is obtained as Eq. (3.10).
[ςij ] , (det J) K× ε̄¯×KT
Tuv =
]i [ςij ]ij [βv,j
]j dKr
Now, if [ςij ](x) is constant across an element, it can be pulled out of the integral.
Otherwise, [ςij ](x) can be interpolated using a Lagrangian polynomial basis {ℓs (ξ)}.
[ςij ](x) ≈
[ςij ]|x=xs ℓs (ξ ◦ z−1 (x))
In Eq. (3.11), {ℓs } denotes the appropriate Lagrangian basis functions corresponding to an interpolation nodal set {ξ s }. It is clear that the Lagrangian polynomials
are originally defined over Kr and transformed onto Kp . Hence, given the interpolation nodal set {ξ s } on Kr , the interpolation can be built from interpolant values
at the corresponding physical node set {xs |xs = z(ξ s )}. Therefore, for the sake of
compactness we define ςijs to represent as [ςij ]|x=xs and have:
ςijs , [ςij ]ij |x=xs = (det J) K× ε̄¯×KT x=xs
Hence, defining κui,vj as Eq. (3.13) we end up with Eq. (3.14) in which S stands
for the number of interpolation points or dim{ℓs } equivalently. In this work, a 2nd
and a 5th -order equidistant interpolation nodal sets, i.e S = 10 and S = 56, have
been used for rectilinear and curvilinear elements, respectively. However, it is clear
that the order of the polynomial basis {ℓs } can be increased/decreased based upon
necessity. In other words, UAs with various polynomial orders for the Lagrangian
{ℓs } basis can be precalculated, stored and used based on the geometrical/material
complexity of individual elements.
κui,vj , [βu,i
]i [βv,j
Tuv =
∑ ∑
0≤s<S i,j
κui,vj ℓs (ξ) dKr
Next, our attention is turned to the exploitation of the algebraic symmetries of
Eq. (3.14). First it is observed that if the integral part of Eq. (3.14) is precalculated,
evaluation of each Tuv requires 9 × 9S multiply-add operations which is due to 9
operations for each ςijs entry times 9 operations for the summation in the ij indices
as formulated in Eq. (3.14). However, closer examination of Eq. (3.13) reveals that
κui,vj=κvj,ui . This suggests that evaluation of (Tuv+Tvu )/2 and (Tuv−Tvu )/2 instead of
Tuv and Tvu can be achieved with reduced computational complexity. Hence, using
some algebraic manipulation the more symmetric formulation of Eq. (3.15) through
Eq. (3.18) is obtained:
Tuv±Tvu =
◦ ±uivjs
κui,vj ±κuj,vi
ℓs dKr
1 + δij
κui,vj ±κuj,vi
ℓs dKr
2 + 2δij
ijs , ςijs ±ςjis
(Tuv ±Tvu ) =
◦ ±uivjs
Implementation of the above formulation leads to total
operations per matrix entry as it seeks to find two T entries. i.e. Tuv + Tvu and
Tuv −Tvu , at once and because C−
ijs is identical to zero for cases where i = j.
The Stiffness Matrix
Similar to 3.4.3, the electrical field formulation for electromagnetic radiation/propagation
problem yields the following expression for the FE stiffness matrix S:
Suv =
¯−1 ∇x ∧ β⃗vx )dKp
∇x ∧ β⃗ux · (µ̄
¯ represents the magnetic permeability tensor. Let’s take β⃗vx as one of
In Eq. (3.19), µ̄
the vector polynomial basis over Kp and expand it using Eq. (3.5). We would like to
express the curl of β⃗ux in a handy form. First, we take advantage of the vector identity
∇ ∧ (f F⃗ ) = ∇f ∧ F⃗ + f ∇ ∧ F⃗ and the fact that ∇ ∧ (∇ϕ) = 0 for any sufficiently
smooth function ϕ:
∇x ∧ β⃗vx=
∇x [βv,m
]m ∧∇x ξm +[βv,m
]m ∇x ∧∇x ξm
]m ∧ ∇x ξm
∇ [βv,m
Knowing that ∇x f =
m ∂ξm ∇ ξm
∇x ∧ β⃗vx =
we get:
∑ ∂[βv,n
∇x ξm ∧∇x ξn
Since ∇x ξm ∧∇x ξn = −∇x ξn ∧∇x ξm , it further develops that:
∑ ∂[βv,n
∇x ∧ β⃗vx =
∇x ξm∧∇x ξn
In Eq. (3.22) the summation over m, n is replaced with m < n which is due to the
anti-symmetric nature of the exterior product in ∇x ξm ∧∇x ξn . Next, given the µ̄
¯−1 = p,q x̂p µ̄
¯−1 x ⃗ x
tensor as µ̄
p,q x̂q we expand µ̄ ·∇ ∧ βv as:
¯−1 ·(∇x ∧ β⃗vx ) =
∑ ∂[βv,n
]n ∂[βv,n
(x̂q ·∇x ξm∧∇x ξn)
x̂p µ̄
¯−1 ·∇x ∧ β⃗vx , κuij,vmn and
Now following Eq. (3.19), in order to evaluate ∇x ∧ β⃗ux · µ̄
ςijp,mnq are defined in the following equations:
ςijp,mnq , (∇x ξi ∧∇x ξj · x̂p )(∇x ξm ∧∇x ξn · x̂q )
( ξ
∂[βu,j ]j ∂[βu,i
]i ∂[βv,n
]n ∂[βv,m
Hence, it is important to observe the following symmetries for κuij,vmn and ςijp,mnq :
κuij,vmn = κvmn,uij
ςijp,mnq = ςmnq,ijp
ςijp,mnq κuij,vmn = ςmnq,ijp κvmn,uij
In order to fully exploit the abovementioned symmetries, one must head for the
evaluation of 21 (Suv ±Svu ) instead of Suv .
Suv ±Svu =
i < j, p
m < n, q
pq ςijp,mnq (κuij,vmn ±κvij,umn ) det J dKr
At this point, we turn our attention to the µ̄
pq ςijp,mnq det J part of Eq. (3.29).
This term bares all the physical information on the geometry and material properties
of the element. Hence, as a function of x, it can be replaced with an appropriate
Lagrangian interpolation similar to what was used in section 3.4.3. Thus, as reflected
in Eq. (3.33), an ‘s’ subscript is added to ςijp,mnq and the summation indices. Namely,
we write ςijp,mnq,s with the ‘s’ subscript signifies evaluation at x = xs . In other words
we have:
pqs , µ̄pq (x)|x=xs
ςijp,mnq,s , ςijp,mnq (x)|x=xs
det Js , det J|x=xs
Hence, Eq. (3.29) is reformulated as:
Suv ±Svu =
i < j, p
m < n, q
pqs ςijp,mnq,s det Js ℓs (ξ)(κuij,vmn ±κvij,umn )dKr
Now, taking advantage of the (ij) ↔ (mn) symmetry, Eq. (3.33) is written as Eq.
∑ ˆ κuij,vmn ±κvij,umn
1 + δij,mn
i < j, p
Suv ±Svu =
m < n, q
ij ≤ mn
pqs (ςijp,mnq,s ±ςmnp,ijq,s )ℓs (ξ) det Js dKr
Except for the µ̄
pqs term, the above representation is symmetric in the p ↔ q sense.
Hence, the symmetry can be further exploited to yield Eq. (3.35).
∑ ˆ κuij,vmn ±κvij,umn
Suv ±Svu =
1 + δij,mn
ij ≤ mn
pqs ± µ̄qps )(ςijp,mnq,s ±ςmnp,ijq,s )
det Js ℓ(ξ)dKr
1 + δpq
In Eq. (3.35), the bound ij < mn reflects the cases where 31 i + 30 j < 31 m + 30 n.
Finally, evaluation of 21 (Suv ±Svu ) reduces to the one time evaluation of Suijvmns and
the element-wise evaluation of the coefficients C±
ijmns as formulated in Eq. (3.37)
followed by the multiply-add operations encoded in Eq. (3.38).
κuij,vmn ± κvij,umn
ℓs (ξ)dKr
2 + 2δij,mn
∑ (µ̄
det J
1 + δpq
(Suv ±Svu )=
ij,mn,s Suij,vmn,s
ij ≤ mn
To evaluate the actual number of multiply-add operations required for evaluation
of the stiffness matrix, one needs to observe that C−
ij,mn,s = 0 for cases where ij = mn.
¯−1 ¯−1
At the same time, for evaluation of C−
ij,mn,s itself, the (µ̄pq − µ̄qp ) term in Eq. (3.37)
equals to zero when p = q. This leads to 6 × 6 × S multiply-add operations for the +
case and 3×3×S multiply-add operations for the − case. Thus, considering that these
operations yield both Suv and Svu , the formulation requires a total of (6×6+3×3)S/2
multiply-add operations per matrix entry.
part of the boundary with essential BCs
ith sub domain of Ω
relative magnetic permeability tensor in Ωi
electrical field in sub domain i
tangent component of Ei , i.e. πτ (Ei )
outward normal to Γij
J imp
imposed electrical current
γτ (•)
trace operator n̂∧•
πτ (•)
tangential trace
(n̂∧•)∧ n̂
∑ operator
del operator i ∂xi x̂i
tangent component of the del operator
curl operator
divergence operator
L0 (Ωi )
{u ∈ L2 (Ωi )|u = 0 on ΓD }
H r (Γij )
Sobolev space of regularity r
H(curl, Ωi ) {u ∈ L2 (Ωi )|∇∧u ∈ L2 (Ωi )}
H0 (curl, Ωi ) {u ∈ H(curl, Ωi )|γτ (u) = 0 on ΓD }
H −1/2 (Γij )
dual space to H 1/2 (Γij )
Table 3.1: Conformal-DDM Related Notation
The UA approach presented here is intrinsically capable of handling elements with
curved geometry. Also, our experiments show that the method can be successfully
Figure 3.2: (a) Edge representation of the curved element. (b) Face representation of
the curved element. Only two faces are plotted here.
Figure 3.3: The cubic Lagrangian interpolation of the curved tetrahedron.
Figure 3.4: Visualization of a subset of the curved tetrahedral elements used in a mesh
for the Luneburg lens. The surface of the lens is plotted with some transparency such
that inside-sphere edges are also visible.
used for evaluation of curved FE matrices such as (iso)-parametric FEs and nonuniform rational basis splines (NURBS) enhanced FEs [65]. The key question, however, is how to build the required bijective coordinate transformation z : Kr → Kp .
Often, as it is the case with (iso)-parametric elements, a Lagrangian interpolation
of an appropriate set of nodal position samples is used for z. For the case of the
spherical Luneburg lens it is easy to identify the tetrahedra that share faces and/or
edges with the curved surface of the lens. For such elements, an analytical expression
for the curved and rectilinear edges can be written in terms barycentric coordinates
as it is shown in Fig. (3.2). The edge expressions are then used to build a triangular
transfinite interpolation [66] of the faces as depicted in Fig. (3.2). The resulting face
expressions are then used to build a tetrahedral transfinite interpolation of the physical volume Kp . Hence, the final transfinite mapping conforms exactly to the surface
of the sphere. As depicted in Fig. (3.3), the resulting coordinate mapping is then
fitted into a Lagrangian interpolation on Kr (3rd -order in this case) and the resulting
polynomial mapping is used as the actual coordinate transformation z : Kr → Kp .
Nevertheless, the method discussed here can be used for arbitrary curvature as far as
it is possible to have analytical expressions for the curved edges (or faces). It is worth
mentioning that such curvature information is usually available in terms of NURBS
in the CAD model and can be made available to the FE engine at the element-matrix
evaluation level.
CIMPT Approach for the PML
Accurate realization of conformal PML layer conditions leads to artificial MPTs
that are continuous functions of global coordinates and often involve curved element
geometries. As a direct consequence, the CIMPT approach will be very profitable for
the realization of such wave absorbing media.
Here we consider a spherical shell PML wrapped around a spherical Luneburg
lens. MPTs for the conformal PML medium are obtained by means of analytical
continuation of the coordinate variable that is normal to the mesh truncation surface
and extends into the complex variable space [67]. This is achieved through field
transformations [68] [58]. Numerical implementations of the conformal PML in both
the FDTD and in the FETD can be found in [69] and [70], respectively. Here, the
expressions for the PML MPTs are briefly provided. The anisotropic MPTs in a PML
can be expressed as Eq. (3.39) and Eq. (3.40) where ε0 and µ0 respectively denote
¯ denotes a relative
the free space permittivity and permeability constants and where Λ̄
permeability (permittivity) tensor.
¯ = µ0 Λ̄
ε̄¯ = ε0 Λ̄
¯ can be written as Eq. (3.41) in
In a 3D spherical coordinates system (r, θ, ϕ), Λ̄
which sr (r) is defined as in Eq. (3.42). The detailed derivation is conducted in [71].
( )2
r̂ + ϕ̂sr (r)ϕ̂ + θ̂sr (r)θ̂
sr (r)
sr (r) = α(r) + ȷ
r,θ,ϕ = r̂
r = Rin + ζ3
In Eq. (3.41), r̃ is a variable obtained by means of complex coordinate stretching
as r̃ = Rin + ζ̃3 where ζ̃3 is defined as in Eq. (3.44). In Eq. (3.42), ω and σ(r)
respectively denote the angular frequency and the conductivity of the PML while in
Eq. (3.43) Rin denotes the inner radius of the spherical PML layer and ζ3 is the
difference between r, and Rin signifying the normal distance form the inner surface
of the PML. Finally, in Eq. (3.44) σ(κ) and α(κ) are analytic functions of κ [38].
ζ˜3 =
ζ3 (
s(κ)dκ =
α(κ) + ȷ
¯ is expressed in spherical coordinates and must be transformed
In Eq. (3.41), Λ̄
into the Cartesian coordinate system for computer implementation. Obviously the
resulting MPTs would be continuously varying functions of Cartesian coordinate
variables. These MPTs can be effectively handled using the CIMPT approach. Some
numerical results concerning the proposed CIMPT PML are presented in section 3.6.
Numerical Results
This section presents the numerical results obtained using the CIMPT enabled
Conformal DDM approach. The tetrahedral FE mesh used in all examples were
generated using the Cubit meshing software. The waveguide problem and all HFSS
results are obtained on a laptop computer with 8 GB of RAM and an Intel⃝
CoreTM i7
CPU. Larger examples, e.g. the large Luneburg lens, are solved on a workstation
computer with 48 GB of RAM and two Quad-Core
processors. The Krylov
subspace solver’s iterations were truncated at a relative tolerance of 10−4 .
To validate the consistency of the results obtained using the UA approach, we
begin with the solution of the WR-15 waveguide problem discussed in [55] which
involves calculating the |S11 | for a waveguide segment loaded by a spatially varying
lossy dielectric. Both the definition and the |S11 | solution of the waveguide problem
are depicted in Fig. (3.5). The problem is solved using two methods: 1) a multilayer
piecewise constant material model 2) a CIMPT model handled by the UA method.
CIMPT 3-Layer
864Tet 504Tet
CPU Time 3sec
Field Basis 2
Mat. Basis 2
Ilic et al.
(4, 2, 4), (4, 2, 7)
Table 3.2: Solution statistics for the waveguide problem. The “2” in the “Field Basis“
row indicates the use of a second order curl conforming basis, i.e. the curl of the basis
is complete to the 2nd order. In the reference result in [55], (m, n, p) represents the
orders of the field basis along (x̂, ŷ, ẑ) directions in hexahedral elements (see Fig.
(3.5)). By material (Mat.) basis we refer to the scalar basis used for the interpolation
of the MPT.
It is evident from Fig. (3.5) that the solution of the multilayer model approaches
to that of the CIMPT model. Using the conventional method, one needs many layers
to obtain results comparable to those of the CIMPT approach. On the other hand
Figure 3.5: (a) Problem geometry. (b) |S11 | versus frequency.
when material properties are allowed to vary, say in this case as 1st -order polynomials,
the actual field solution will be more complex compared to the piece-wise constant
material case. Thus, field basis with higher orders of polynomial interpolation will
be needed before one can fully enjoy the DoF savings associated with the CIMPT
approach. This is more evident by looking at the statistics provided in Table. (3.2)
where similar mesh dimensions (around λ0 /4) and same basis orders are used for all
cases. Also, note that the DoF numbers reported by [55] are obtained by means of
an optimal choice of the basis orders along various directions. In other words, [55]
exploits the fact that the actual field solution for the waveguide problem has little or
no variation along the transverse directions. This, allows for a strongly anisotropic
choice of the basis orders along different directions. Hence, anisotropic orders of
(4, 2, 4) and (4, 2, 7) are used in the air and dielectric elements respectively where
(m, n, p) refers to the orders of the field basis along (x̂, ŷ, ẑ) directions.
Figure 3.6: (a)Electrical field intensity. (b) Conformal-DDM’s Partition of the Mesh.
Luneburg Lens
The Luneburg lens problem discussed here is a dielectric sphere with a relative
permittivity tensor of ε̄¯ = (2 −
r2 ¯
where r is the distance from the center of the
sphere, R is the actual radius and ¯Ī is the identity tensor in R3 . Besides Luneburg’s
original development on the optical lens, there is plenty of work on a variety of
Luneburg lens designs used in antenna/reflector devices operated in microwave and
millimeter wave frequencies [72] [73]. Hence, the Luneburg lens is a good example for
the purpose of verifying the fidelity of the CIMPT approach.
Conceptually speaking, the Luneburg lens focuses any incoming plane wave into
a point on the surface of the lens. This can be exploited to build highly directive
antennas by essentially placing a receiving/transmitting element at the focus of a desired direction on the surface of the lens. The focusing effects can be clearly observed
in what will be presented in the proceeding sections. This is helpful as a means for
qualitative validation of the results, e.g. Fig. (3.6).
Figure 3.7: Schematic diagram showing the configuration of the quarter-Luneburg
lens scattering problem. Curved elements were applied to the surface of the lens only.
Scattering from an R = 6λ0 Lens
The example solved here is a 6λ0 (λ0 is the free space wavelength) radius lens
exposed to an incoming plane wave at 10 GHz. Fig. (3.7) provides a schematic
diagram of the problem geometry. In this case, the problem symmetry is exploited
and only a quarter of the lens is used. The problem mesh and geometry are created
using the Cubit mesher with a maximum discretization size of λ/2.5. Fig. (3.6) shows
the associated DDM partition of the FEM mesh as well as a plot of the resulting
electrical field intensity. The plane-wave focusing effect can be clearly seen in Fig.
(3.6). The 64, 592 tetrahedron (723, 070 DoF) mesh was automatically partitioned
into 100 sub-domains with the aid of Metis graph partitioning library. Sub-domain
matrix assemblage took 28 seconds with the memory usage peaking to 1.8 GB. The
preconditioner was assembled in 33 seconds and demanded 3.3 GB of memory (No
disk caching was used) while the final solution was obtained in 88 iterations (13
The Lens Antenna
Our next experiment involves exciting lens antennas of various radius and calculating their far field patterns. By reciprocity, it is understood that a point source on
the surface of the sphere should excite a plane wave radiating from the lens. various
lenses of radius R ∈ {λ0 , 2λ0 , 3λ0 , 4λ0 , 5λ0 , 6λ0 } are excited by a WR-90 waveguide
at the frequency of 10 GHz in the T E10
. As depicted in Fig. (3.8), the aperture of
the waveguide is located 5mm apart from the outer surface of the R = 6λ0 lens. The
resulting near field solution for the R = 6λ0 case is plotted in Fig. (3.9). One can
observe the flattening of the phase front as the outgoing waves propagate away from
the aperture of the waveguide.
Figure 3.8: Schematic diagram showing the configuration of the R = 6λ0 Luneburg
lens antenna problem. Curved elements were applied to the surface of the lens only.
Figure 3.9: A three-slice field plot of the E · ŷ field component in the continuous index
Luneburg lens antenna with R = 6λ0 and f = 10 GHz. The shaded part of the plot
belongs to the dielectric lens.
Figure 3.10: The radiation pattern (dB) of Luneburg lens antennas of various radius
R as function of θ at ϕ = 0.
Using DDM, relatively modest computational resources are needed for computation of relatively large problems. Here, we simulated Luneburg lens antennas with
lens diameters of up to 12λ0 . For the R = 6λ0 case, a mesh size of λ/2.5 results
in a discretization with 423, 133 tetrahedra and 4, 790, 122 DoFs (no symmetry was
exploited here). Unlike the quarter lens, the absorbing boundary condition (ABC)
surface is placed on the surface of a surrounding cube (Fig. (3.8)). The DDM solution
involves 500 automatically generated partitions . Assembling sub-domain matrices
and excitation vectors took 5 minutes with a peak memory usage of 4.4 GB. Six
minutes were spent on the construction of the preconditioner with a peak memory
usage of 18 GB (No disk caching was used). Note that the preconditioner’s peak
memory usage is directly associated with the size of the largest partition obtained by
Metis. Hence, the number of DDM partitions was increased to keep the sub-domain
dimensions more or less comparable to those of section 3.6.2. It took 80 iterations to
finish the solution process amounting to one and a half minute of wall clock time. Fig.
(3.10) plots the resulting radiation patterns for the mentioned set of lens antennas
with radius of R ∈ {λ0 , 2λ0 , 3λ0 , 4λ0 , 5λ0 , 6λ0 }.
Conformal PML using CIMPT
As a last example, we report the implementation of a conformal PML by means
of the CIMPT approach. Fig. (3.11) includes a schematic drawing of the problem
geometry. Fig. (3.12) shows snapshots of the field solution of a 1λ0 radius Luneburg
lens antenna wrapped in a continuous PML. The 0.25λ0 thick PML is placed 0.25λ0
away from the surface of the lens. The attenuation function of the PML is set to
sr = 1 + ȷσmax ζ3 /2πε0 in which ζ3 is the distance from any point within the PML
region to the inner surface of the PML layer. Here we use σmax=1. The geometry is
meshed by Cubit using a maximum discretization size of λ/2.5. The PML exhibits
Figure 3.11: Schematic diagram showing the configuration of the PML-lens antenna
Figure 3.12: (a) Snapshots of the magnitude of the electric field distribution on yzplane. (b) Snapshots of the magnitude of the electric field distribution on zx-plane.
strong attenuation of the radiated field which validates its proper implementation.
Fig. (3.13) plots the resulting radiation pattern in the ϕ=0 plane.
An efficient UA approach capable of handling curved geometries and continuously
varying MPTs is introduced in this chapter. Compared to the conventional approach
Figure 3.13: The radiation pattern (dB) of the Luneburg lens antenna with R = λ0
at ϕ = 0
for evaluation of FE matrices, the presented UA approach assumes that material
properties and element Jacobian related terms can vary as continuous functions across
the element. With the aid of appropriate interpolation schemes the required crosselement continuity or discontinuity of the desired terms can be guaranteed. This leads
to more flexible and accurate modeling of the underlying physics.
Merged with a conformal DDM-FEM solver, the UA approach is validated on a
few problems including Luneburg lenses and a spherical perfectly matched layer for
which the MPTs were again modeled as continuous functions of spatial coordinates.
Relatively large electrical problems with continuously changing material properties
were solved and presented for the purpose of numerical validation.
Despite the current verification on the Luneburg lens problem or the conformal
PML, the CIMPT technique is believed to be useful in other areas such as: a) multiphysical problems where dielectric properties change due to other physical phenomena
while re-meshing of the problem domain needs to be avoided, e.g. dielectric permittivity may change due to temperature/pressure fluctuations; b) in doped semiconductors
and plasma antennas/reflectors, dielectric properties may follow a continuous profile
following various physical phenomena. Hence, the introduction of the UA method will
play an influential role in future developments associated to the CIMPT approach.
Finally, it worth mentioning that complete utilization of the benefits of the proposed
CIMPT and UA methods can be achieved by means of implementations with higher
orders of field basis functions.
Chapter 4: Ferromagnetic Nano-wires: Homogenization and
A numerical homogenization method to extract the equivalent electromagnetic
parameters, in particular, the dispersive and nonreciprocal permeability tensor, of
ferromagnetic nanowires (FMNW) is presented in this chapter. With the effective
material property obtained, the performance of FMNW microwave devices, such as
integrated double-band isolators and self-biased circulators, etc., are analyzed using a
full-wave numerical simulation approach. Numerical examples validate the proposed
homogenization method and demonstrate the integration of the multi-physics solution
strategy to characterize the effects of FMNWs employed in RF devices.
Metamaterials can be broadly defined as artificial effective structures composed
of small units and achieving exotic responses that are not readily available in natural
or conventional materials. Extensive research conducted on metamaterials has led to
many innovative theories and devices over the years [74] [75] [76] [77]. Coming to the
second decade of the new century, the development of artificially structural materials
is at the stage where interdisciplinary interest is taking place, such as multi-scale
metamaterials, reconfigurable metamaterials and the combination with nanotechnology. Compared with the earlier generation of metamaterials, the new generation of
metamaterials incorporate finer structural scales, for instance, the nanometer or even
atomic scale, which suffer less from scattering loss as well as fabrication difficulties.
While the potentials of metamaterials and their possible applications are appealing, it is often difficult to model them accurately due to inherently the very fine
and complex structures involved. However, by observing that metamaterials usually
consist of positioned inclusions (most often periodically) much smaller than the electromagnetic wavelength of interest, a homogenized description of the structure can
be adopted as an accurate alternative. Among numerous work already existing in
the literature on homogenization, a few highlights are listed here: in [46], the authors
present a homogenization method by calculating the reflection and transmission coefficients on finite lengths of the electromagnetic metamaterials; in [47], the authors
have determined the effective material constants based on averaging the local fields;
and in [78], a homogenized description is derived of periodic metamaterials, which
are made of magnetodielectric inclusions, and subsequently, closed-form expressions
for the effective constitutive parameters are obtained. Finally, in [79] the author
have combined the physical insight and the Whitney-like interpolation to form a
new homogenization theory of metamaterials, through which coarse-grained fields are
Here we are interested in ferromagnetic nanowire (FMNW) array and exploit its
magnetic properties as well as investigate its integrations with microwave/RF components and devices. Due to their tailored electromagnetic responses, FMNW metamaterials exhibit novel properties such as double ferromagnetic resonance (FMR),
anisotropy control and self-bias. Compared against traditionally bulky ferromagnets,
the utilization of FMNW metamaterials can significantly reduce the size in microwave
devices. The major previous work referred by this work are [29] and [30], in which an
analytical Maxwell-Garnett like formulation is presented for extracting the effective
permeability tensor for non-saturated/saturated arrays of axially magnetized FMNWs as well as enumerating several applications of FMNWs in microwave devices.
Instead of adopting the analytical model, this work proposes a three-dimensional
(3D) full-wave numerical process more generally to extract the macroscopic effective
permeability tensor, by first computing the static magnetic field distribution with the
acceralation of two-dimensional Fast Fourier Transformation method (2D-FFT) ( [80])
and then conducting a small signal analysis ( [81]). With that, the performance of microwave components/devices integrated with FMNWs based metamaterial structures
are analysed through full-wave electromagnetic computations.
This chapter begins with an introduction to put in context the definition and
development on metamaterial structures, the demand and challenge on modeling
them. It then proceeds to state the homogenization methodology, and numerical
result section is included afterwards.
The Proposed Homogenization Method
We consider an array of ferromagnetic nanowires of diameter d and length L
(L >> d) embedded in porous membrane with height h and inter-wire distance D,
as shown in Fig. (4.1). Our goal is to obtain the effective dynamic permeability of
the array in a complex tensorial form.
Figure 4.1: Sketch of FMNW substrate structure.
The proposed homogenization procedure can be divided into the following major
Setting up the biasing/operating point
The FMNW array is subjected externally to axially-aligned (set to be z direction
in this scenario) DC magnetic field Hext
dc and a transverse ac signal h. The initial
magnetization Mdc , therefore, can be obtained from the hysteresis loop of the FMNW
metamaterials, which relates the magnetization to the applied static magnetic field
dc . Fig. (4.2) shows a schematic hysteresis curve.
Computation of the demagnetization field
With the magnetization Mdc obtained, the volume/surface magnetic charge density ρv /ρs in/on the tube is calculated as:
ρv = −∇ · Mdc
ρs = Mdc · n̂,
Figure 4.2: Schematic hysteresis loop.
with n̂ the surface normal of the tube pointing in the outward direction. Consequently,
the magnetic potential Φ can be computed via:
Φ(⃗r) =
′ dv +
r − ⃗r |
V |⃗
′ ds .
r − ⃗r |
S |⃗
Since in our case Mdc is constant, it results in ρv = 0; thus Eq. (4.3) can be
simplified as,
Φ(⃗r) =
′ ds .
r − ⃗r |
S |⃗
In general there would be millions or even billions of nanowires involved in the
structure, the direct computation of Eq. (4.4) requires prohibitive computational
resources. Fortunately, as stated previously, these nanowires are usually distributed
in a periodic lattice, we can therefore accelerate the computation of Eq. (4.4) by
employing 2D-FFT method [80].
Observing Eq. (4.4), the magnetic potential Φ due to the magnetic charge can be
decomposed into near and far interaction terms:
Φ(⃗r) = Φn (⃗r) + Φf (⃗r).
In Eq. (4.5), Φn is the potential contributed by disk charge residing on the top
and bottom surface of the tube and evaluated on the axis of the disk. A closed form
expression for the near term can be expressed as:
Φn (z) =
ρs √
[ (z − z0 )2 + R2 − |z − z0 |],
where ρs is the surface magnetic charge density, R is the radius of the disk, and z and
z0 are the potential evaluation location and the center of the magnetic charge disk,
Φf in Eq. (4.5) can be derived approximately from point charge of ρs πR2 , with R
being the radius of the cylindrical tube, and assume that R << ⃗r − ⃗r , where ⃗r and
⃗r are the locations of the charge and the observation point, respectively. We have:
Φf (⃗r) ≈
1 ∑ ρs πR2
|⃗r − ⃗r′ |
It has already been pointed out before, the nano-tubes are distributed in a Cartesian grid with source charges of the same value but opposite sign residing on the
top and bottom surfaces of each tube. Combined with the fact that static Green’s
r − ⃗r′ ), it leads to a
function is translational invariant, namely, g(⃗r, ⃗r ) = |⃗r−⃗
′ = g( ⃗
r |
2-D block Toeplitz matrix structure for Eq. (4.7). Therefore, the computations of the
magnetic potentials can be speeded up significantly using 2D-FFT for the summation
of far-term contributions. Adding the near-term correction, Eq. (4.6), at the end
completes the computations of the magnetic potentials on the Cartesian grids.
The demagnetization field Hm is then readily available for every ferromagnetic
nanowire from the magnetic potential:
Hm (⃗r) = −∇Φ(⃗r).
Computation of the small signal solution of the LLG
Assuming that the corresponding ac magnetization field h is much smaller than
the external DC bias, a decomposition of the total effective magnetization field Hef f
into the DC and the small ac signal hence proves to be beneficial. Namely,
Hef f = Hdc + h,
Hdc = Hext
dc + Hm .
Note that other contributions to the effective magnetic field Hef f , such as the
exchange interaction, anisotropic field and the thermal field, are neglected since they
are much less significant compared against the demagnetization field.
Similarly, the total magnetization moment M is separated into DC and ac components as well:
M = Mdc + m.
The DC magnetization, Mdc , is mainly in the z-direction in the current consideration. Moreover, the ac magnetization moment, m, can be computed from the
equation of magnetization motion, or the Landau-Lifshitz-Gilbert (LLG) equation:
= −µ0 |γg |M × (Hef f +
µ0 |γg |Ms ∂t
where µ0 is the permeability of free space, γg is the gyromagnetic ratio, αG is the
damping factor, and Ms is the magnitude of the saturation magnetization.
Substituting Eqs. (4.9) and (4.10) into Eq. (4.11) and assuming a time harmonic
dependence of ejωt , and applying the small signal approximation model in [81] by
neglecting higher order terms, Eq. (4.11) can be written as:
jωm = −µ0 |γg |(M ẑ × h + m × Hdc ) − jωαG
ẑ × m,
with M being the magnitude of M.
Eq. (4.12) can be written into:
M my
M mx
= −µ0 |γg | (M hx + mz Hx − mx Hz ) − jωαG
jωmx = −µ0 |γg | (−M hy + my Hz − mz Hy ) + jωαG
jωmz = −µ0 |γg | (mx Hy − my Hx )
The 2 × 2 matrix form of x and y components is:
Axx =
Axx Axy
Ayx Ayy
= µ0 |γg | M
−jHx Hy Ms |γg |2 µ0 2 +Hz Ms |γg |µ0 ω−jM αG ω 2
Axy =
Ayx =
Ayy =
j(Hx 2 |γg |2 µ0 2 −ω 2 )
j(−Hy 2 |γg |2 µ0 2 +ω 2 )
jHx Hy Ms |γg |2 µ0 2 +Hz Ms |γg |µ0 ω−jM αG ω 2
where Hx , Hy and Hz are the components of Hdc in the x, y and z directions,
respectively, and mx , my and mz are those of m. The notations hx and hy are also
defined similarly.
Extraction of the equivalent/effective permeability tensor
With the small signal magnetization moment, m, computed in response to given
external applied ac magnetic field, h (through the excitation of three different polarizations), it follows that the 3×3 magnetic susceptibility tensor χω can be determined.
Consequently, the relative equivalent/effective homogenized permeability tensor µω
of FMNWs is also available, and we note that the extracted permeability tensor is
nonreciprocal and dispersive.
m(ω) = χω h(ω)
χω = 
µ0 |γg |M
Axx Axy
Ayx Ayy
µxx µxy 0
µω = χω + I =  µyx µyy 0  .
0 1
It is worth commenting that in the proposed procedure, in cases where the magnetization moment locally is not largely homogeneous, the extracted macroscopic effective material properties can further be made as inhomogeneous within the FMNW
Integration with microwave/RF components/devices –
Performance Analysis
The extracted electromagnetic material parameters will be employed directly for
microwave/RF components utilizing FMNW substrates. With the macroscopic material properties readily available, a full-wave 3-D electromagnetic simulation will be
applied to analyze the microwave/RF components, such as isolators, circulators, etc.
Fig. (4.3) summarizes the methodology and scope of this work, including procedures to perform the homogenization technique of FMNW based metamaterials and
its applications as RF components in microwave devices.
Figure 4.3: Flow-chart of FMNW based metamaterial homogenization process and
application enumeration.
Numerical Results
Validation of extracted permeability tensor
We first validate the proposed numerical homogenization procedure by comparing
the effective material property computed against an analytical model proposed in [82]
for a nanowired substrate of different magnetic alloys. Specifically, a nonreciprocal
microstrip phase shifter on a substrate with two different magnetic nano wires (CoFe
and NiFe), was designed and analyzed. The substrate is fabricated via electrodeposition of CoFe and NiFe nanowires into 100µm thick porous anodic alumina membranes
with pore diameter of 35nm and membrane porosity P = 12%. In order to enhance
the nonreciprocal behavior of the device, CoFe and NiFe nanowires are grown into
different heights of the membrane thickness, with hCoF e = 0.68 and hN iF e = 0.84,
The analytical homogenized model proposed in [82] for the nanowires has the
µ −jκ 0
µ =  jκ µ 0  ,
0 1
2πMs γP hfF M R (m2 + 1)
fF2 M R − f 2
4πMs γP hf m
κ =
fF2 M R − f 2
µ = 1+
where fF M R = γ[HDC + 2πMs (1 − 3P )] + jαf is the ferromagnetic resonance
frequency, m = M/Ms is the normalized magnetization, α is the damping factor,
HDC is the applied static field parallel to the axial direction of the nanowires (for
self-bias case, it will be zero), f is the operation frequency, Ms is the saturation
magnetization and finally P is the porosity.
Figs. (4.4) and (4.5) below show the comparisons of the effective relative permeability tensor obtained from the proposed homogenization procedure and the reference
analytical model in self-bias situation. Note these structures have high remanence
(magnetization at zero field, after saturation), so it is definitely still valid for applying a small signal model analysis.
From these figures, we have observed that the two models predict resonances
at almost the same frequency for both CoFe and NiFe nanowire substrates, with
fres CoF e = 22.3GHz and fres N iF e around 10GHz.
Figure 4.4: Comparison of effective relative permeability tensor of CoFe nanowire
membrane obtained by different methods.
Figure 4.5: Comparison of effective relative permeability tensor of NiFe nanowire
membrane obtained by different methods.
Double ferromagnetic resonance phenomenon
FMR has been studied extensively due to their applications into microwave devices. Specifically, in FMNW arrays, the FMR phenomenon draws great attention
mainly due to the following reasons: First of all, the FMR frequency is easily tunable by either adjusting the external magnetic field or applying various demagnetizing
cycles so to set the remanent state programmable. Secondly, it is predicted by analytical model and observed through experiments that for applied magnetic fields below
saturation, there exist two sets of FMR frequencies. The existence of double FMR
frequencies can potentially lead to applications with dual working frequency bands.
A schematic drawing of a microstrip line with the substrate of FMNW embedded
in a porous alumina membrane is included in Fig. (4.6). As shown in the figure, the
microstrip line is 0.5mm wide and 16mm long, and the bottom of the membrane is
covered with metal as ground. We shall study the electromagnetic wave propagates
along the microstrip line over a wide frequency band and under a range of external
applied magnetic fields.
Figure 4.6: Microstrip line sketch.
We took the hysteresis loop figure from [83] and replotted it in Fig. (4.7). Particularly, in Fig. (4.7), we highlight the realized magnetization branch within the
hysteresis loop. Through the figure, the pairs of the applied external DC magnetic
field Hdc
and the normalized magnetization Mn (Mn = Mdc /Ms ) are obtained. Note
that the saturation magnetization Ms for the FMNW is 1400kA/m. We conduct the
studies of FMR by sweeping the external DC bias field from −5.0KOe to 5.0KOe
with a step of 500Oe, following the specified branch indicated in Fig. (4.7). Table
(4.1) summarizes the (Hdc
, Mn ) pairs read from the hysteresis loop.
Figure 4.7: Major hysteresis loop (solid blue) and minor hysteresis curves (dashed
4.1: External DC Magnetic field vs.
-2.0 -1.5 -1.0 -0.5
-0.96 -0.9 -0.82 -0.68
2.0 2.5
-0.002 0.24 0.48 0.69 0.85
Normalized Magnetization
-0.48 -0.25
0.97 1.0
With the initial DC magnetization obtained from the hysteresis loop, we perform
the homogenization process described in full previously and outlined in Fig. (4.3)
to extract the relative effective permeability tensors of the FMNW from 20GHz to
40GHz. Subsequently, the extracted permeability tensors are employed to compute
the S21 coefficients through full-wave electromagnetic computations. However, it is
worth mentioning that in order to account for the fabrication imperfections, fringefield effects, and the measurement uncertainties, we adopted a similar adjustment as
suggested in [29] [83]. Namely,
µc = I + q(µh − I),
where µc is the effective permeability tensor that is used in the full-wave simulations,
and µh is the one extracted from the homogenization process. The parameter q is a
geometrical filling factor, which is determined experimentally according to references
[29] [83]. Herein we have q = 0.13. Note that in references [29] [83], the FMNW is
approximated as a reciprocal and isotropic material with the effective permeability
given by:
µc = 1 + q(µh − 1),
where µh is the diagonal component of the permeability tensor of the FMNW substrate.
By extending Eq. (4.21) to Eq. (4.20), we aim to capture the non-reciprocal effects
of the FMNWs. Although, with the adjustment introduced, the non-reciprocal effects
induced will not be significant. For example, we have the effective permeability, µc
at 33GHz and the external applied magnetic field Hdc
= −5KOe, as:
1.068 − j0.2 −0.2 − j0.055 0.
µc =  0.2 + j0.055 1.068 − j0.2 0.  .
In Fig. (4.8), we compare our numerical results against the measurements and
the results from [29]. The numerical results computed using the proposed approach
is shown on the left of Fig. (4.8), whereas the reference experimental and analytical
model results copied from Figure 2, in [83] are included on the right of Fig. (4.8).
As can be seen from the figure, very good agreements can be established between
these three sets of results. As a consequence, through this example, the validity
of the proposed numerical homogenization process is demonstrated. Furthermore,
judging from the figure, it can be argued that the proposed homogenization procedure
produced results that resembled the measurements more than the analytical model
proposed in the refence.
Figure 4.8: Left: contour plot of the S21 parameter from simulated result. Right (a)
and (b): reference experiment and model results.
It is also observed, from the figures in Fig. (4.8), that with the external applied
DC magnetic field above the saturation, |Hdc
| ≥ 2.5KOe, there would be only a
single FMR peak. Additionally, in the region above the saturation, the resonant
frequency increases almost linearly with the magnitude of the external applied DC
magnetic field. However, it is in the region where the applied field is below saturation,
0 < Hdc
< 2.5KOe, that two sets of peaks corresponding to two FMR frequencies
exist. The existence of two FMR frequencies could result in potential applications for
dual-band microwave devices. In Fig. (4.9), we plot the electric field distributions
along the microstrip line, with a single external bias field set to be 0.5KOe. As can
be seen from the figure, there are two different frequencies, 22GHz and 30GHz, at
which the electromagnetic energy gets absorbed almost completely by the FMNW
Figure 4.9: Electric field distributions along the microstrip line with FMNW subext
strate and biased with the same external bias field, Hdc
= 0.5KOe, at two different
resonance frequencies.
Dual-band integrated edge-mode isolator
In this section we first compute the homogenized relative permeability tensor of
a FMNW substrate previously illustrated in [29]. This substrate is utilized in a
dual-band integrated edge-mode isolator in [30], which will be discussed in detail
following the validation of the homogenized material properties. The dimensions of
the substrate as well as the geometrical description of the isolator are included in Fig.
Figure 4.10: FMNW isolator structure.
Comparison of the homogenized relative permeability of the FMNW substrate
The nanowire array is embedded in a porous alumina membrane by selective
deposition. The wires are distributed in a symmetrical network. The average diameter
and length of the nanowire are 40nm and 200µm, respectively, with inter-wire distance
about 110nm. Referring to Fig. (4.1), we have d = 40nm, L = h = 200µm, and
D = 110nm, as documented in Fig. (4.10).
With the external biased magnetic field being −0.5KOe in the axial direction,
subsequently the corresponding average magnetization moment is 350kA/m (Ms =
1400kA/m). Following the same procedure in [29], the relative characteristic permeability µc is obtained from the homogenized relative permeability tensor via µc =
1 + q(µh − 1), µh being the diagonal component of the tensor here. The comparison of
µc at a few frequency points against the results obtained from [29] is plotted in Fig.
(4.11) below. As can be seen, computed results from the proposed method agree well
with the reference.
Figure 4.11: Comparison of relative characteristic permeability. Solid lines stand for
analytical model in the reference while marks of triangles and squares denote the
results of the proposed numerical homogenization method.
A dual-band integrated edge-mode isolator
After validating the proposed homogenization procedure, we are ready to investigate some interesting properties of FMNWs employed in microwave devices. Herein
consider a dual-band integrated edge-mode isolator designed in [30]. Fig. (4.10) elucidates the geometry and the corresponding dimensions of the isolator. The relative
permittivity of the porous alumina membrane (without FMNWs) is 5.3 and that of
the alumina membrane with FMWN inclusions is 6.15 − j0.25. The operation bias
point of the isolator is set to be Hext
dc = −1.3KOe and the corresponding magnetization moment is −140kA/m. Figs. (4.12) and (4.13) plot the S parameters of
the isolator obtained by different methods. As can be seen, compared against the
measurements, all simulation results have a lower insertion loss, which may be due to
the neglect of loss incurred in the physical realization of the device and measurement
Figure 4.12: S12 parameter comparison.
Moreover, as shown in Figs. (4.12) and (4.13), the isolator works at two isolation
frequency bands, centering at 25GHz and 32GHz, respectively. At the two working
frequency bands, the isolation directions are reversed. Subsequently, the electric field
distributions along the microstrip line of the isolator is plotted in Fig. (4.14) for these
two center frequencies.
A dual-band self-biased circulator
Shown in Fig. (4.15) is a planar dual-band self-biased circulator integrated with
FMNW substrate. The FMNW substrate utilized by the circulator is taken directly
from reference [84]. The porous alumina substrate is 220µm thick with the average
Figure 4.13: S21 parameter comparison.
pore diameter and inter-pore distance are 40nm and 100nm, respectively. Magnetic
nanowires are selectively electro-deposited into the membrane. The saturation magnetization of the nanowire array is 1400kA/m and the remanent magnetization is
Mr = 0.48Ms . The back of the template is covered with perfect electric conductor
(PEC) as the ground plane. Furthermore, a PEC circular mask is aligned on top, and
underneath the disk is the nanowire array region. The radius of the disk is 2.5mm
and the three-port access microstrip is 8mm long and 0.5mm wide. With the effective
permeability tensor of FMNW substrate obtained through the proposed homogenization procedure, and the relative permittivity of porous alumina membrane with and
without magnetic nanowire fillings set to be 6.15 − j0.25 and 5.3 ( [84]), respectively,
the full-wave numerical simulation can be launched to evaluate the performance of
the circulator.
Fig. (4.16) plots the transmission coefficients, with Port 1 being excited, over the
frequency range, from 20GHz to 40GHz. As evidenced from the figure, the circulator
has two working frequency bands: 25GHz to 29GHz and 31.5GHz to 33GHz. In the
Figure 4.14: Electric field distributions along the microstrip line of the isolator at two
different working frequencies.
first frequency band, 25GHz to 29GHz, the electromagnetic wave energy circulates
clockwise, whereas in the higher frequency band, 31.5GHz to 33GHz, it circulates
counter-clockwise. The reverse of the circulation direction at these two frequency
bands is clearly demonstrated in Fig. (4.17), where the distributions of the electric
field are plotted for 27GHz and 32GHz, respectively.
With a modest isolation of roughly 10dB, the performance of the circulator is yet
to be enhanced. However, as a proof of concept, the circulator proposed demonstrates
the capability of utilizing FMNW substrate in integrated microwave devices to achieve
novel properties and new applications.
A numerical homogenization procedure is proposed in this chapter to study the
unique properties of FMNW metamaterials. Numerical results validate the proposed
Figure 4.15: Detailed description of a dual-band self-biased circulator.
approach. The double FMR phenomenon of ferromagnetic nanowired array has been
demonstrated. Moreover, two microwave components, a microwave edge-mode isolator and a dual-band self-biased circulator, have been analyzed using the homogenized
material properties in full-wave time harmonic electromagnetic field computations.
Figure 4.16: Full-wave simulation result of the conceptual microstrip circulator–S
parameter plotting.
Figure 4.17: Full-wave simulation result of the conceptual microstrip circulator–
Electric field plotting.
Chapter 5: Self-Consistent Integration of Discontinuous
Galerkin Time Domain Method and Circuit Simulator for
Microwave Components and Devices with Circuit Network
Previous chapters of this dissertation are devoted to solutions to time-harmonic
Maxwell’s equations in freuquency domain and their applications. In this chapter,
a time-domain discontinuosu Galerkin method integrated with circuit simulator for
modeling transient and nonlinear circuits and devices is presented.
The motivation for developing a stable and efficient solver in this area is selfelaborate. Modern integrated circuits have been working at faster operating speeds
with higer integration densities, traditional analysis based on equivalent circuit model
for the whole structure is not accurate enough to represent the complicated strucuture.
It imposes the accuracy requirement for both electromagnetic field and circuit device
simulations. This demands sophisticated and efficient numerical simulation tools.
With the exibility and capability that DGTD provides, it is well suited for the
modeling of such complex system of electronic packaging. DGTD is appealing due
to several well-known reasons. This method can well address elements of various
types and shapes, non-conformal meshes, and non-uniform degrees of approximation.
Furthermore, allowing for discontinuity of the discrete trial and test spaces, DGTD
can provide a substantial amount of flexibility in the choice of basis functions.
As far as the circuit simulator concerned, SPICE, the simulation program with
integrated circuit emphasis, is a powerful program that is embedded with abundant
model libraries for semiconductors, logical gates and integrated circuits. Therefore
the hybridization of DGTD and SPICE has the capability to tackle complicated and
multi-port circuit devices accurately and efficiently.
However, in many other real-life situations, semiconductor vendors always try
to protect proprietary information and are reluctant to provide transistor designs
that underlie the fabrication process. Without the knowledge of the circuit elements and their connections, SPICE cannot be used anymore; fortunately, there has
been an emerging standard accepted by increasingly many electronic design automation (EDA) vendors, semiconductor vendors and system designers for the behavioral
modeling of IC input/output characteristics: the Input/output buffer information
specification (IBIS) standard. IBIS is a specification that defines the format on the
information about the input/output buffers of the integrated circuit product without revealing propretary, internal processes or architehtural information. Although
IBIS model is not as accurate as SPICE model, since IBIS depicts IC input/output
characteristics, it is faster than SPICE models. Therefore, the integration of fullwave electromagnetic solver with IBIS can provide efficient solution for wide-ranged
Be it SPICE or IBIS, the interface between the EM and circuit parts is similar.
Generally speaking, the interface of the two parts is defined as a planar surface port,
on which the voltage across is calculated in EM simulation at each time step, and
coupled to the circuit solver. The circuit solver then calculates the current on the
port and radiates back to the EM simulation. The kep contribution of this work is to
present a self-consistent scheme with convergent looping process for the EM-circuit
The rest of this chapter is organized as below. In 5.1, the summary of DGTD
formulation and its integration to SPICE solver will be illustrated followed by some
vanlidation examples. Later in 5.2, the IBIS model will be introduced in detail first.
It is then followed by the illustration of its integration to DGTD solver. Again some
interesting numerical examples will be presented.
DGTD integrated with SPICE solver
We present herein an approach to integrate the DGTD and the SPICE solver
in this section. The simulation problem of EM structure and circuit devices is decomposed into an EM subsystem and a circuit subsystem. The two subsystems are
coupled through planar surface ports. A self-consistent scheme is proposed to deal
with the coupling between them within a convergent looping form. Two different
port establishing strategies are investigated for the EM-circuit couplings. Several numerical examples are shown to validate the accuracy and application of the proposed
With the ever increasing operating speeds and integration densities of modern
ICs, successful prediction of practical PCB designs imposes an accuracy requirement
for both electromagnetic field and circuit simulations [85] [86]. On one hand, EM
fields at high frequencies would cause wave phenomena such as EM radiation, signal
delay and distortion, skin and proximity effects that can hardly be generated by
pure circuit simulation. On the other hand, semiconductor chips or devices contain
complicated linear/nonlinear components, and play central and crucial roles in the
proper functioning of IC designs. Therefore, mixed electromagnetic and circuit cosimulation have become imperative and should be utilized in a co-design process of
the components, such as microwave structures, and active/passive linear/nonlinear
circuits, for the PCB products.
Over the past decades, considerable efforts have been made to incorporate EM
effects into IC simulation. These works can be roughly classified into two categories.
One is a hybrid method of bringing the circuit effects into full wave EM simulation,
or co-simulation [52, 87–97]. In such methodology, both EM and circuit simulations
are required and a proper interfacing scheme to couple the EM system and circuit
system is expected. The other one is to represent the EM effects in terms of circuit
simulation [32] [98]. That is, to extract the equivalent circuit models of EM structures
by full wave EM simulation over a range of frequencies, or to export the EM effects
as distributed circuit sources, and simulate them together with the devices within
only circuit simulator. While worthwhile, the latter method suffers from accuracy
problems, especially when at high frequencies. In contrast, the first method can well
capture both EM and circuit information through exact simulation of two systems.
Therefore, we will focus on full wave EM simulation in this work.
The DGTD method [25, 39, 52, 94–97, 99–101] offers an appealing alternative for
solving time-domain Maxwell’s equations. This method can well address elements
of various types and shapes, non-conformal meshes, and non-uniform degrees of approximation. In addition, the method can lead to a fully explicit and conditionally
stable time-marching scheme; the mass matrix arising from discretization will be
block diagonal, which enables the local inversion of matrix in each element. Furthermore, information exchange is required only between adjacent elements regardless
of approximation orders and shapes. This feature makes it highly suitable and efficient for parallelization [102]. Therefore, the application of DGTD in EM and circuit
co-simulation would also be promising.
Another important factor in EM-circuit co-simulation is the method for circuit
simulation. In [93], FETD is employed for EM simulation and modified nodal analysis (MNA) simulation is employed to build the circuit solver. Similarly, in both [96]
and [97], DGTD is used for EM simulation, and MNA is used for the circuit simulation. In these works, the equivalent circuit representations of metal-semiconductor
field-effect transistor (MESFET), for example, were built inside the EM-MNA code
and the MNA equations were coupled directly with FETD/DGTD equations to set
up co-simulation. However, it is well known that a general circuit simulator, such
as SPICE, is itself a MNA solver [103–106]. The drawback of the above-mentioned
methods, therefore, is that they tried to reimplement what SPICE does, but not as
robust as SPICE stands. Because all SPICE models, such as diodes, bipolar junction
transistors (BJTs), metal-oxide-semiconductor field-effect transistors (MOSFETs),
has to be redefined in these EM-MNA solvers, as well as the MNA matrix solvers and
numerical nonlinear solvers related to the convergence options in SPICE, probably in
a same way, as required for real industrial application purpose. As an illustration, the
models of MOSFETs in HSPICE has approximate seventy model levels [105], with
different characterization models and parameters (such as the gate oxide thickness,
temperature dependence coefficient, etc.), and the number keeps growing. These
model levels should be indispensable to circuit analysis of real-life devices, but reimplementing these parameters and definitions inside an EM solver seems to require too
much work. Therefore, we believe designing a hybrid EM-circuit algorithm with a
general-purpose circuit simulator is a better way to deal with real-life devices as how
SPICE does. In this way, we can relate the devices directly to physical parameters
and process characterization, and utilize the circuit solver, models and definitions of
existing SPICE software.
Existing work on DGTD and SPICE co-simulation can be found in [52]. The authors have formulated a lumped circuit subcell model to link the DGTD simulation
with SPICE solver. An exponential time difference (ETD) algorithm are adopted with
a fourth-order Runge-Kutta time integration method to generate a time-marching
scheme. Meanwhile, the circuits are solved by calling external SPICE software. Nevertheless, DGTD and SPICE are coupled through an explicit way and not connected
rigorously. Therefore, the accuracy and performance of the method may become impaired in presence of highly nonlinear devices if a time step is not carefully chosen to
be small enough.
In [94], a discontinuous Galerkin time-domain method is shown by attaching
lumped elements directly within the DGTD solver, where the lumped elements are
treated as planar impedance surfaces. In this section, we extend the previous work
of [94] of modeling lumped elements to a co-simulation scheme by integrating DGTD
and SPICE. A hybrid IC structure consisting of electromagnetic components and circuit devices is divided into two subsystems, namely, one containing electromagnetic
structures, one containing circuit devices. Then, planar surface ports are built to
make links between the two subsystems. After that, the EM part is simulated in
DGTD and the circuit part is simulated in SPICE simultaneously with respect to
self-consistency in an iterative manner. The developed method offers transient and
wide frequency-band computations and the inclusion of arbitrary linear/nonlinear
passive/active devices easily, for complex configured integrated circuit simulation and
The rest of the contents begins with a brief review of the interior penalty DGTD
(IP-DGTD) formulation and time discretization. Next, we describe the way to integrate DGTD simulation with a SPICE simulator. A self-consistent looping scheme
is derived with detailed explanation of the coupling mechanism between two simulators. The strategies to establish the coupling ports are carried out in two different
viewpoints. Finally numerical results to validate the developed method are presented.
Original Boundary Value Problem
The problem in this work we consider is described by the time-dependent Maxwell’s
equations in bounded three-dimensional domain Ω as,
∇ × H = ϵ̄¯
∇ × E = −µ̄
¯(r) vary in space
where the electric permittivity ϵ̄¯(r) and the magnetic permeability µ̄
Ω. On PEC and PMC surfaces inside Ω , we apply the boundary conditions as
n̂ × E = 0
on ΓP EC
n̂ × H = 0
on ΓP M C
where n̂ is the unit normal on these surfaces. On the boundary ∂Ω of the computational domain, we adopt the 1st order absorbing boundary condition (i.e. ∇ × E =
−cµ∇ × ∇ × H and ∇ × H = cϵ∇ × ∇ × E).
DGTD Formulation
As previously discussed in [94], in the computational domain of interest Ω, let
Th be the discretization of Ω into tetrahedral elements {Ki }. We denote the set of
all faces as Fh = FhI ∪ FhB , where FhI represents the set of interior faces of adjacent
elements of Th , ∂Ki ∩ ∂Kj , and FBI the set of boundary faces ∂Ki ∩ ∂Ω. Moreover,
we define two trace operators, the twisted tangential trace operator γτ (·) and the
tangential components trace operator πτ (·), which will be employed throughout our
derivations, as
γτ (ui ) = n̂i × ui |∂Ki
πτ (ui ) = n̂i × ui × n̂i |∂Ki
where n̂i is the boundary normal pointing out of the element Ki .
Following the procedure presented in [99, 100], we can derive the interior penalty
DGTD formulation for the first order Maxwell’s equations in time domain. By defining the finite-dimensional discrete space as Vhp = {v ∈ [L2 (Ω)]3 : v|K ∈ [Pp (K)]3 , ∀K ∈
Th }, with L2 (Ω) the linear space of square-integrable vector fields in domain Ω, and
Pp (K) the polynomial space of order p in element K, we can state the DGTD formulation as:
Find (H, E) ∈ Vhp × Vhp such that
¯ · ∂H
w· ∇×E+
ϵ̄¯ · ∂E
− v· ∇×H+
{{v}} · [[H]]γ dS −
{{w}} · [[E]]γ dS
[[v]]π · [[E]]π dS + f
[[w]]π · [[H]]π dS = 0,
∀ (w, v) ∈
Vhp .
where the notations are defined as follows,
 {{u}} = πτ (ui ) + πτ (uj ) /2
[[u]]γ = γτ (ui ) + γτ (uj )
on FhI
[[u]]π = πτ (ui ) − πτ (uj )
 {{u}} = πτ (u)
[[u]]γ = γτ (u)
on Fh
[[u]]π = πτ (ui )
For the value of e and f , we choose e = 2Z1 Γ and f = 2Y1Γ with ZΓ = 21 ( µi /ϵi +
µj /ϵj ) and YΓ = 21 ( ϵi /µi + ϵj /µj ), to obtain an upwind numerical flux. The
choice will lead to a slightly lossy formulation but provide an optimal O(hp+1 ) convergence rate as shown in [25].
DGTD Time Discretization
Assume the number of degrees of freedom in element Ki is di , we can express
the local electric and magnetic fields as Eh (r, t) = dni ein (t)vin (r) and Hh (r, t) =
∑d i
n hin (t)win (r) in terms of basis functions v, w, where ein (t) and hin (t) are the time
dependent coefficients for the basis functions. By separating the w, v testing in Eq.
(5.4), we are able to build a local semi-discrete system within element Ki as,
= Se hi − Fiie hi + ePiie ei
e hj + ePe ej
= −Sh ei + Fiih ei + f Piih hi
h ej + f Ph hj
where ei (t) and hi (t) are the time dependent coefficient vectors, respectively. j is the
indices of the neighboring elements of element i. The local matrices are expressed as
(Mϵ )nm =
vin · ϵ̄¯i vim dV
(Se )nm =
vin · ∇ × wim dV
(Fe )nm =
πτ (vin ) · γτ (wim )dS
2 FhI
(Fe )nm =
πτ (vin ) · γτ (wjm )dS
2 FhI
(Pe )nm = −
πτ (vin ) · πτ (vim )dS
(Fe )nm =
πτ (vin ) · πτ (vjm )dS
(Miµ )nm
¯i wim dV
win · µ̄
(Sih )nm
win · ∇ × vim dV
πτ (win ) · γτ (vim )dS
2 FhI
πτ (win ) · γτ (vjm )dS
2 FhI
πτ (win ) · πτ (wim )dS
πτ (win ) · πτ (wjm )dS
(Fiih )nm
h )nm =
(Piih )nm =
h )nm =
The first-order time derivatives in the resulting differential Eqs. (5.7) can be
discretized using central difference approximation with second order accuracy. The
electric field unknowns are evaluated at tn = nδt and the magnetic field unknowns are
evaluated at tn+1/2 = (n + 1/2)δt. For the two penalty terms arising from upwind flux
n+ 1
n+ 1
formulation, we apply a backward approximation as ei(j)2 ≈ eni(j) and hn+1
i(j) ≈ hi(j) ,
in order to achieve an explicit time marching scheme. As a result, we can express the
fully discretized local system of equations in a leap-frog way as follows,
n+ 12
Miϵ en+1
= (Miϵ + eδtPiie )eni + δt(Sie − Fiie )hi
n+ 21
e hj
n+ 32
Miµ hi
+ eδtPij
e ej
n+ 12
= (Miµ + f δtPiih )hi
+ δt(−Sih + Fiih )en+1
n+ 21
+ f δtPij
h ej
h hj
In addition, local time stepping technique is employed to reduce the computational
time for multi-scale structures and the stability analysis would be the same as mentioned in [39, 94, 101].
DGTD-SPICE Integration
In this section, we describe how we integrate DGTD simulation with a generalpurpose circuit simulator, such as SPICE. These two solvers are coupled through
planar surface ports, with a self-consistent coupling scheme.
Coupling Scheme
Consider a microstrip line structure terminated with a device or circuit network.
We use this system to illustrate how to address EM-circuit computation we are interested. We first decompose the whole system into two subsystems, one including the
electromagnetic structure, the other including the circuits, as shown in Fig. (5.1).
For the microstrip part, we employ full-wave EM simulation; while for the circuit
Figure 5.1: Sketch of EM-Circuit System
network, we use SPICE simulator. Furthermore, we propose a self-consistent scheme
to couple the two decomposed subsystems, that is to ensure
through a planar surface port as shown in Fig. (5.1). In Eq. (5.9), VEM
the voltage across the surface port, and IEM
represents the current flowing on the
planar surface in the EM structure. VCKT
denotes the voltage across the circuits,
and ICKT
denotes the total current in the circuits in SPICE, respectively.
The coupling mechanism comprises of two directions: from EM coupling to circuit,
we obtain the voltage across the surface port and use it as the equivalent driving source
for the SPICE simulator; from circuit coupling to EM, we take the current calculated
from circuit solver and use it as the equivalent current source on the surface port
for full-wave simulation. The process is demonstrated in Fig. (5.2). This coupling
mechanism alone, however, can only ensure Eq. (5.9b), but not Eq. (5.9a). Note
in [52], similarly Eq. (5.9b) is ensured but not Eq. (5.9a). That is the reason and
where we add the self-consistent scheme in.
Figure 5.2: (a) EM coupled to SPICE (b) SPICE coupled to EM
To illustrate our point, we denote the general current and voltage relation for
SPICE circuit network as Eq. (5.10a). With the current radiating back to DGTD
system, the voltage is calculated from full-wave simulation, and we denote this process
as Eq. (5.10b).
ICKT = f (VCKT )
VEM = g(IEM )
By applying one of the self-consistent conditions Eq. (5.9b), we can obtain
VEM = g(IEM ) = g(ICKT ) = g(f (VCKT ))
Here, f is the general function that describes current and voltage relations for arbitrary circuit network or semiconductor devices, thus it can be either linear or nonlinear; function g links current and voltage in full-wave EM system, and it is linear.
Information exchange between DGTD and SPICE goes in a VCKT −
→ VEM fashion. Note that Eq. (5.9a) requires VCKT and VEM be the same. If only Eq. (5.9b)
is ensured, the difference of VCKT and VEM could be noticeably large after such a
information flowing in presence of a highly nonlinear f , which would result in a nonphysical or unstable solution as time goes by. This problem may be avoided by choose
a small enough time step to keep the difference of VCKT and VEM sufficiently small.
Nevertheless, that degree of “small” is not easy to control and it could cost too much
computational resources.
As a result, Eq. (5.9a) must be kept true throughout co-simulation. Generally
speaking, gf = 1 is not necessarily true. However, we can find certain VCKT and VEM
pairs, satisfying both Eq. (5.9a) and Eq. (5.11) within some preset tolerance, from
numerical side’s view. Thus, we have the following self-consistent looping process:
• set an iteration index as k and a tolerance value as ϵ.
• Start from VCKT
• Compare VCKT
with VEM
; If their difference is bigger than ϵ, denote
+ ∆V
Otherwise, stop looping.
• To find VCKT
= g(f (VCKT
)), apply the Taylor’s expansion of the function
VEM = g(f (VCKT )) at VCKT
, ignore the second and higher order terms, and
we get
)) = g(f (VCKT
)) +
= g(f (VCKT
| k ∆V
• Since VEM
= g(f (VCKT
)), associating with Eq. (5.12), we have
∆V =
| k −1
is calculated as Eq. (5.12), and it is served as the input voltage
• Therefore, VCKT
to circuit network for next iteration, until convergence is achieved.
As can be seen from above, the derived looping process is similar to the Newton0
Raphson method [107]. And here the initial guess is set to be VCKT
. Therefore, the
looping method should also share the property of the Newton-Raphson method that
provides a fast convergence rate for the iterations.
SPICE coupling to DGTD
Within full-wave simulation for Eq. (5.10b), thanks to DG method, only minor
modification is needed for elements in the vicinity of the circuit port with the extra
current source. Retrieving the previous work on lumped elements in DGTD method
[94], in general we model them as planar impedance surfaces, as shown in Fig. (5.3).
Compared to wavelength, the port dimension is relatively small, thus we can assume
the current densities are constant on the surface. Without loss of generosity, we can
address two elements Ki , Kj that share a face ∂Ki ∩ ∂Kj on the surface ΓCKT as
shown in Fig. (5.3).
Figure 5.3: Geometrical illustration of the impedance surface.
Starting from the boundary conditions of the electric and magnetic fields on the
port, we can obtain the relationships on surface ΓCKT as follows,
n̂i × Hi + n̂j × Hj = JCKT =
n̂i × Ei + n̂j × Ej = 0
The magnetic field relationship Eq. (5.15) is used to represent the currents from the
circuits. Meanwhile, the electric fields Eq. (5.16) need to be tangentially continuous
across the surface port since no magnetic currents exist on the surface. As a result, the
two relationships need to be weakly enforced through the interior penalty approach
similar to Eq. (5.4). The corresponding residuals then can be written as,
Function Space
RΓCKT = [[H]]γ − JCKT
RΓCKT = [[E]]γ
∈ H(divτ , ΓCKT )
∈ H(divτ , ΓCKT )
As shown in [94], the residuals RΓCKT and RΓCKT represent the surface error electric
and magnetic currents, respectively. According to the dual pairing principle of form(1)
ing reaction integrals, RΓCKT should be tested with a surface E to form the energy
density E · Jerr , thereby we adopt πτ (v) as the testing functions for them, which lies
in the space H(curlτ , ΓCKT ). Similarly, RΓCKT should be tested with surface H and
we choose the testing functions as πτ (w).
With respect to the above analysis, we can replace the interior penalty formulation
of Eq. (5.4) for the elements that share the faces on circuit surface ΓCKT , as follows,
¯ · ∂H
w· ∇ × E +
ϵ̄¯ · ∂E
− v· ∇×H+
{{v}} · ([[H]]γ − JCKT ) dS
{{w}} · [[E]]γ dS = 0
Separating the w and v testing, we can obtain local the semi-discrete equations for
element Ki as
= Se hi − Fiie hi + BJe JCKT − Fij
e hj
= −Sh ei + Fiih ei + Fij
h ej
where (BJe )n =
πτ (vin )· l̂ and l̂ denotes the unit vector along the length direction
of the port. Likewise as Eq. (5.8), we apply the central difference approximation to
time derivatives and evaluate JCKT at tn+ 2 = (n + 21 )δt. The resulting local timedependent update equations can be written as
n+ 12
Miϵ en+1
= Miϵ eni + δt(Sie − Fiie )hi
n+ 1
n+ 12
+BJe JCKT2 − δtFij
e hj
n+ 32
Miµ hi
n+ 12
= Miµ hi
+ δt(−Sih + Fiih )en+1
h ej
DGTD coupling to SPICE
With regard to Eq. (5.10a), we execute SPICE for circuits until a certain time step,
save the continuous quantities of the circuits like voltages/currents/charges/powers,
and let it wait for running command. After DGTD produces a new value of voltage,
we change the parameters in the circuits, apply the stored quantities of circuits as
initial conditions, and continue the circuit simulation from the previous time step.
To clarify the implementation process, we assume the following responses of EM
, vcap
and circuit systems are known at time tn = nδt. These responses include VCKT
inind , and also VCKT
. vcap
denote the voltages across the capacitors and inind denote
the currents in the inductors in the SPICE circuits at tn = nδt. These are the input
to drive the circuits in SPICE, which are saved from the simulation results of DGTD
from tn to tn+1 and the simulation results of SPICE from tn−1 to tn . We can write
Eq. (5.10a) more specifically with all the input/output variables as,
n+ 1
n+1 n+1
, iind ) = fSP ICE (VCKT
, vcap
, inind )
(ICKT2 , vcap
The relationship in Eq. (5.22) indicates the SPICE simulation from tn to tn+1 .
The voltage pair VCKT
are used to construct a piecewise linear voltage pulse
in SPICE with VCKT
at tn and VCKT
at tn+1 , in order to drive the whole circuit.
and inind are used as the values of initial conditions of capacitors and inductors
to re-obtain their states at previous time step. A typical input file is shown in Fig.
n+ 1
As for the output shown in Eq. (5.22), the computed port current ICKT2 is flowing
back on the planar surface in DGTD as addressed in Eq. (5.21). vcap
and in+1
ind are
saved and will be used as the initial conditions for the simulation of next time step,
namely from tn+1 to tn+2 . Note for a multi-port circuit network, multiple voltages
sources should be built, with multiple port currents coupling back to DGTD.
Hybrid Co-simulation
Combining with the previously mentioned formulations and implementations, we
state the proposed self-consistent DGTD-SPICE method with the following major
Figure 5.4: SPICE netlist example.
• Start from time tn , and assume VEM
• Set the iteration number k = 0, and make a guess of the initial voltage value
across the circuit as VCKT
• Reset the SPICE configuration commands by updating the voltage source and
the initial condition of the circuit, including all historical circuit status parameters such as node voltages, current flows through all inductors and all DC values
of the sources at tn .
• SPICE takes voltage as a piece-wise linear voltage source defined from tn to
tn + δt, with the range from VCKT
• With all initial parameters and the voltage source, SPICE performs transient
n+ 1 ,k
analysis from tn to tn + δt via Eq. (5.22), and JCKT2 is captured and returned
to the DGTD routine.
n+ 1 ,k
• In EM system, perform DGTD simulation with the current source JCKT2 via
Eq. (5.8) and Eq. (5.21). Then compute VEM
from the line integral of the
electric field through the circuit port.
• Calculate VCKT
by using the relations of Eq. (5.14) and Eq. (5.12).
• If
the nth to (n + 1)th time step marching is completed. We save VCKT
) and
the states of the circuit at tn+1 , and begin the co-simulation at next time step.
If Eq. (5.23) is not satisfied, increase k and go back to continue looping until
Fig. (5.5) shows the flowchart of the methodology. Note although we illustrate
our method by using only one port, the co-simulation method can be easily extended
to describe multi-port cases as well.
As pointed out by [97], sometimes only network information such as S/Y-parameters
is known for circuits. This situation requires DGTD solves with circuit network representing by admittance matrix. Note the proposed DGTD-SPICE method intrinsically
have this feature as well. An insightful discussion on SPICE-compatible representations of S-parameter matrices (yet not the only way) can be found in [108]. With the
S/Y-parameters provided or tabular S/Y-parameters measured of real multi-port circuit network or devices, we can employ the rational approximation method first, like
Figure 5.5: Flowchart of DGTD integrated with SPICE
on vector fitting basis [109,110], and then write the netlists for the corresponding S/Yparameter or admittance matrices. More specifically, either S/Y-parameter matrices
or admittance matrices can be represented in SPICE with a combination of E-element
(voltage-dependent voltage sources), F-element (current-dependent current sources),
G-element (voltage-dependent current sources), and H-element (current-dependent
voltage sources), along with the “LAPLACE” option in HSPICE [105] or “XSPICE”
option in NGSPICE [106].
Surface Port
As aforementioned, the information exchange between EM and circuit systems are
achieved through planar surface port. To achieve an accurate and physical response,
the surface port in Fig. (5.3) should be built in a way that resembles the real physical
dimension. This gives rise to the problem on how to connect the port of a finite
width to a circuit of zero width, and to approximate the problem of reality in a most
physical way. Based on different perspectives to view the relationship between the
voltage across the surface and the current density on it, we propose two ways to
establish the port for the integration of DGTD and SPICE.
We call the first way an average voltage method. Recalling that the voltage
distribution at the end of a microstrip-like structure is nonuniform, the voltage in the
middle location of the port is larger than it is at the edges. Therefore, we can define
the voltage across the surface port in an average sense,
´W ´L
E(x, y) · l̂dydx
V (x)dx
= 0
where l̂ denotes the unit vector along the length direction of the port. The integration
of voltage in Eq. (5.24) is calculated on the port by adopting the Gaussian quadrature
rule [107], as shown in Fig. (5.6),
1 ∑
w(xi )V (xi )
W i=1
where n is the number of Gaussian quadrature points, w(xi ) is the weight corresponding to the location of the Gaussian point. The value obtained by Eq. (5.25) is passed
to SPICE through the relation of Eq. (5.12), Eq. (5.14), and Eq. (5.22). Then by
applying the voltage to excite the circuits in SPICE, a single value of current ICKT
can be obtained by circuit simulation, which will radiate back to DGTD framework
as the current density in Eq. (5.21).
Another method is named an average current method in short. We start to build
the port by directly express the link between EM voltage and the current density on
the port as:
J(x) = F (VEM (x))
Figure 5.6: Integration by Gaussian quadrature rule
where VEM (x) is the voltage obtained from EM simulation at location x and serves
as the driving source for circuit simulation and J(x) is the corresponding current
density on the circuit port at x. F (·) is the general function that describes the
current density and voltage relations for arbitrary circuit network on the surface, it
can be either linear or nonlinear. From J(x) to form total current, we have:
J(x)dx =
F (VEM (x)) dx
Similarly, we express the integration into a weighted summation by Gaussian integration rule, and obtain:
I =
F (VEM (x)) dx
ˆ 1 (
xi +
wi F VEM
w(xi )I(xi )
As a result, the total current on the port can be seen as a weighted summation of the
currents at different locations along the width direction of the port.
To implement this idea, the voltages are calculated at different locations along the
width direction of the port surface in DGTD, and passed onto SPICE respectively. In
SPICE, we can obtain the current values corresponding to each voltage value. When
returning those values of currents to DGTD, we do the average of the currents with
regard to the integration points and send back the current density as in Eq. (5.21).
Note in this method, however, SPICE needs to be called multiple times since we need
to evaluate I(x) at different locations. Therefore, the second way takes longer time
compared to the average voltage method.
Comparing these two schemes, the first method evaluates the voltage as a weighted
summation of the voltages at different locations of the port, while the latter method
expresses the current as a weighted summation of the currents at different locations
of the port. When the circuit network is linear, it can be easily concluded that
they are essentially the same; but when function F involves nonlinearity, they will
not produce the same results. From the perspective of energy exchange, the average
current method express the energy information V (x)I(x) between the EM system and
circuit system at distributed locations x(i), which should be more reliable than the
average voltage method which handles energy exchange as a whole. We shall make a
comparison between the two schemes in the numerical examples.
Numerical Results
To test the reliability of the developed method, two examples with linear/nonlinear
circuit components attached are presented: one is a microstrip structure; the other is
a step recovery diode (SRD) pulse generator. The developed DGTD-SPICE method
is used to solve these examples. DGTD with 1st order absorbing boundary condition
is used for the EM simulation while NGSPICE [106] is used as the circuit solver.
The models are also built into CST [111], where the simulation results are used as
reference values.
We first investigate a simple example where a microstrip is connected with linear
circuit network, as shown in Fig. (5.7). The geometrical parameters of the EM
structure are w1 = 1 mm, w2 = 2.5 mm, L = 8 mm, h = 2 mm, and the substrate
material is FR-4 with ϵr = 4.726 and tan δ = 0.0255. The input excitation voltage
source is a 600 MHz sinusoidal signal with amplitude of 1 Volts, and the source
resistance is 50Ω. At the termination end, a circuit network similar to [52] is used.
Figure 5.7: Microstrip line with circuit network.
As for the reference, the discrete port definition of CST is based on hexahedron
meshes, which simulates the port with a single-line sampling of voltage. This is not
compatible with our surface port definitions in tetrahedron meshes. However, we can
modify the line port of CST a bit to let it fit in a surface port definition. It is done
by creating multiple line ports in CST, and redistribute the termination network on a
surface by equalizing it into multiple sub-networks in parallel. The values of lumped
elements like resistors, inductors and capacitors in each sub-network is calculated
based on the ratio of the small port surface area to the original port surface area, as
shown in Fig. (5.8). The width, length and position of the surface we “created” in
CST is the same as the one we defined in the tetrahedron meshes for DGTD.
Figure 5.8: CST multi-line port experiment.
The simulation results of DGTD-SPICE and CST are compared in Figs. (5.9)
and (5.10). Since the circuit network is linear, the average voltage method and the
average current method will produce the same results. As shown, the results of the
proposed DGTD-SPICE co-simulation method agree very well with the results from
Figure 5.9: Result comparison at near end of microstrip.
To investigate the difference of our two port methods discussed in Section IV,
we simulate the same microstrip terminated with a nonlinear circuit network that
contains a diode, as shown in Fig. (5.11), and build the ports with both the averaging
voltage and average current methods. The amplitude of the input signal is increased
to 2 Volts. The diode parameters at the termination end are shown in Table (5.2).
For comparison, we also build the same model in CST with the multi-line “surface”
port. We shall explain a bit more on how it is done for the nonlinear circuit network.
For linear lumped elements like resistors, inductors and capacitors, we simply obtain
the redistributed element values exactly according to the surface ratio. For a nonlinear
element such as diode, its parameters are usually not directly in proportion with
port area. However, as SPICE code/model indicates, the current produced by a
Figure 5.10: Result comparison at far end of microstrip.
Figure 5.11: Microstrip line with nonlinear circuit network.
Figure 5.12: Comparisons of numerical results with different schemes.
diode consists of the volume current plus the sidewall current. The volume current is
proportional to the surface area that the diode occupies while the sidewall current is
proportional to the perimeter of the surface area. Fortunately, CST offers a surface
area scaling factor for diode and we use it to scale the diode area as we need, the
process illustrated in Fig. (5.8). In such a way, the diode in the circuit sub-network
can also be correctly represented.
Finally, we have the results of three schemes: the DGTD-SPICE method with an
average voltage port model, the DGTD-SPICE method with an average current port
model, CST with multiple line-port model. The comparison of the numerical results
is demonstrated in Fig. (5.12).
As shown in Fig. (5.12), the average voltage method shows a very smooth curve
while both the average current method and the CST multi-line port models reflect
many small oscillations on top of the smooth curve. Furthermore, the oscillations
from the two curves obtained by the average current method and CST share similar
phenomena. This validates the capability of the average current model to capture
more accurate nonlinear effects.
Overall, from both theory side and numerical experiment side we arrive at the
conclusion that the average current method is more accurate for dealing with nonlinear
network. However, the average current method involves calling for SPICE multiple
times for each time step. This requires more computational resources than the average
voltage method, not to mention combining the iterative process where per solution
of SPICE needs to be recomputed several times until convergence. Therefore, further
acceleration of the method needs to be tackled, and tentatively, the direction of
parallelization of SPICE and DGTD processes can be branched out. However, the
different results produced by the average voltage method and the average current
method that we observe in this example is being exaggerated to some degree, since
this difference actually should be determined by the width of the surface port and
the voltage difference from one side to the center. In real-life PCB, the width of the
circuit elements would be smaller and the voltage difference would be also smaller if
the circuit is not connected to the end of a microstrip.
SRD Pulse Generator
In this example, we consider the numerical simulation of a pulse generator based
on step-recovery diode as shown in Fig. (5.13). The pulse generator consists of
coplanar waveguide layout (substrate material is FR-4) with traces, grounding vias,
an inductor, two capacitors and a step-recovery diode. The detailed geometrical
Figure 5.13: Schematic structure of SRD pulse generator .
parameters of the pulse generator are listed in Table. (5.1). The circuit elements are
represented at the ports in the configuration. Fig. (5.14) draws the circuit layout of
the pulse generator.
The SPICE model of the step-recovery diode that are passed
Table 5.1: Geometrical description of the pulse generator
w1 0.5 mm
w2 1.25 mm
w3 0.20 mm
l1 3.75 mm
l2 0.50 mm
l3 2.50 mm
l4 0.75 mm
0.50 mm
0.40 mm
0.20 mm
Figure 5.14: Circuit structure of SRD pulse generator.
to NGSPICE are listed in Table (5.2). The capacitor in parallel and the inductor in
series with it represent the parasitic elements of the SMD package.
The pulse generator is driven by a sinusoidal voltage source with an amplitude
of 2.36 Volts and frequency of 60 MHz. By examining the structure carefully in
Fig. (5.14), we can get the first understanding of how the pulse generator works. It
is excited at Port 1 with a low-frequency harmonic signal, which passes through a
LC low-pass filter and reaches at the SRD package at Port 5. The SRD is used as
a charge controlled switch for waveform sharpening. In the first half period of the
sinusoidal signal, the SRD is charged by positive voltages and it would behave as
nearly shorted. While in the second half period of the sinusoidal signal, the SRD is
being discharged by negative voltages. Once the charges is removed, the strong nonlinearity of the SRD will cause the impedance to change very sharply from nearly
zero to an extremely large value. In this way, the signal will be blocked by the diode,
filtered by the DC blocking 1.5-pF capacitor, and produces a sub-nanosecond pulse
Table 5.2: Step-recovery diode description
Name Description
Saturation current per unit 0.5 pA
Ohmic series resistance
0.13 Ω
Emission coefficient
Contact potential at area 0.5 V
Grading coefficient at area 0.235
Reverse breakdown voltage 60 V
Current at breakdown volt- 10 µA
Transit time
30 ns
Coefficient for forward-bias 0.8
depletion bottom-wall capacitance formula
at the output side at Port 2. Numerical results in Figs. (5.15) and (5.16) confirm the
theory, where Fig. (5.15) shows the sampled currents at point A of Port 1, and Fig.
(5.15) shows the sampled voltages at point B of Port 2. Moreover, CST reference
results match the numerical results very well.
In addition, in the simulations of both examples, for a relative error tolerance of
ϵ = 1 × 10−6 , the proposed self-consistent method is able to converge at each time
step within an iteration number k of four.
In this section, we have presented a self-consistent approach to solve the problem
of EM-circuit co-simulation by integrating SPICE with DGTD. Initially, the original EM-circuit problem are decomposed into two subsystems. The EM subsystem
Figure 5.15: Result comparison of currents at the input side of the SRD pulse generator.
Figure 5.16: Result comparison of voltages at the output side of the SRD pulse
is solved with DGTD simulation with the currents from the circuit subsystem. The
circuit subsystem is solved with SPICE by using the voltages obtained in EM structure. The two subsystems are coupled through planar surface port. We have shown,
in order to couple the two systems rigorously, a self-consistent condition must be
ensured throughout simulation. Based on this conclusion, we build a way to keep the
two system self-consistently in each time step with an iterative methodology. Moreover, two different ways to establish the coupling ports are discussed and tested along
with their advantages and disadvantages. The application of the method is examined
with two examples including a microstrip and a pulse generator. Good agreement
of the results with references have demonstrated the capability and accuracy of the
developed method to solve mixed electromagnetic and circuit problems.
DGTD integrated with IBIS model
Introduction to IBIS models
When the circuit network structure for semicnoductors, logical gates and integrated circuits cannot be disclosed, which is often the case in real-life situations for
the sake of protection of propreitary information, SPICE cannot be used anymore.
Instead, IBIS is developed as simulation models to meet a variety of customer needs.
Without revealing the underlying circuits structure or process information, IBIS models provide a standardized way of representing the electrical characteristics of a digital
ICs pins (input, output, I/O buffers and the like) behaviorally. It can be used in a
simulation environment to help solve board-level overshoot, undershoot or crosstalk
problems and so on ( [112] [113] [114]).
Figure 5.17: Input buffer.
In general, the data in an IBIS file is used to construct a buffer model useful
for performing signal integrity (SI) simulations and timing analysis of printed circuit boards. The fundamental information needed to perform these simulations is
the buffer I-V (current versus voltage) and switching (output voltage versus time)
Fig. (5.17) shows the structure of an input. It can be viewed as a receiver. It
is composed of two diodes for system-level electrostatic discharge (ESD) protection,
the die capacitance (the buffer’s capacitance) and the package parasitic parameters
of lumped RLC combinations.
For the input buffer, the required information is:
1. The buffer’s I-V characteristics
2. The buffer’s capacitance
3. Package parameters
Figure 5.18: Output buffer.
Fig. (5.18) shows the structure of output buffer. The model can be viewed as a
driver. It consists of a PMOS transistor and a NMOS transistor, two diodes for ESD
protection, the die capacitance (the buffer’s capacitance) and the package parasitics.
An output or I/O buffer is characterized using the following information:
1. The buffer’s output I-V characteristics when the output is in the logic low state
2. The buffer’s output I-V characteristics when the output is in the logic high state
3. The buffer’s output I-V characteristics when the output is forced below ground
and above the power supply rail
4. The time it takes a buffer’s output to switch logic states from low to high and
high to low
5. The buffer’s capacitance
6. Package parameters
As far as the IBIS file format is concerned, it is a readable ACSII file that contains
the first part of the header general information about the file, the device and the
company, then the second part of the device name, pinout and pin-to-buffer mapping,
followed by the third part the I-V and V-T data for each more. If there are more than
one device, the second part and third part are repeated as many times as devices are
DGTD-IBIS co-simulation
Similar to the co-simulation scheme of DGTD and SPICE, the link of electromagnetic components and circuit devices in the IC structure is built at the pin of the
device as planar port. Voltage calculated from EM simulation is sent into IBIS model
then current derived from IBIS model is coupled back to the EM system. Discussions
on the extraction of dynamic information on the pin from IBIS waveform as both
input and output buffers are detailed below in 5.2.3 and 5.2.4, respectively.
Transient representation of IBIS input buffer model
Compared to output buffer model, the input buffer IBIS model integrated into
the co-simulation is relatively easy. Referring to Fig. (5.17), the voltage at the pin
from EM simulation crosses the parasite R pkg, L pkg and C pkg network, the ESD
diodes and C comp, the sillicon die capacitance. At ESD structure, the GND and
power clamp curves of V-I are provided in the IBIS model, representing the electrical
behavior of the buffer when the GND clamp and the power clamp diodes are turned
on, respectively. The GND clamp is active when the buffer is below ground, and the
power clamp is active when it is above VDD. Therefore, to solve for the current on
pin that radiates back to the EM system is fairly straightforward: with the I-V table
to look up for the ESD diodes, the package parameters and die capacitance values
given, a simply circuit solver is realized to obtain the current.
Transient representation of IBIS input buffer model
As previously pointed out in Fig. (5.18), the output buffer of IBIS behavioral
model consists of pullup transistor, pulldown transistor, power clamping diode and
ground clamping diode, ramp time and pad and package parasitics. In the static
domain, current versus voltage tables are used to characterize each of the output
transistors and each of the clamping diodes. Values for rising and falling ramp time
are used to describe dynamic behavior of the output stage. Note that the rising
and falling waveform VT tables are measured under the load condition without the
package parasitic by connecting a fixed resistance direc whose value is speficified in
the file. The VT tables of an IBIS model are purely based on DC condition and
should not be used for transient simulation. How to express the dynamic information
of switching VT tables of IBIS model is critical. Here, a multiplier based method will
be briefly given, the details of this method is referred in [115].
The current flowing through pullup transistor is determined by the voltage across
it referenced to Vcc , noted as Ipu (V ), while the current flowing though pulldown transistor is determined by the voltage across it referenced to Gnd noted as Ipd (V ). Two
time dependent multipliers Ku (t) and Kd (t) are introduced to describe the switching
behavior of pullup and pulldown transistors respectively. The value range of the multipliers is from 0 to 1, where 0 corresponds to transistor turn-off state and 1 turn-on
state. The change of multiplier from 0 to 1 represents the transistor status change
from turn-off to turn-on. The product of K and I is to describe the transistor switching behavior. Thus, for rising edge, Ku (t) changes from 0 to 1 while Kd (t) goes from
1 to 0, and vice versa for the falling edge, as sketched in Fig. (5.19).
Figure 5.19: Illustration of the introduced multiplier switching behavors.
Figure 5.20: Equivalent analog behavior model.
With the introduced multiplier coefficients, the circuit topology of the equivalent
analog behavioral model is illustrated in Fig. (5.20). The pullup characteristic Ipu (V )
and the pulldown characteristic Ipd (V ) are combined with Ku (t) and Kd (t), to describe
the switching behaviror between the High and Low level of the buffer output signal.
The current source Ipc (V ) represents the Vcc -clamping effect and the current source
Igc (V ) represents the Gnd-clamping effect.
Next, we discuss on how to calculate the values of the introduced multipliers.
From Fig. (5.20), the following relation can be derived in Eq. (5.29), where Vdie (t) is
the voltage at the die point at time point t. The values of Ipu (V ), Ipd (V ), Ipc (V ) and
Igc (V ) at that time point can be determined according to IV tables by using piecewise
linear interpolation. Iout (t) can be represented by the load and voltage relation, with
the values of the package parasitic parameters connecting to the given fixed resistance
in the presence of the external voltage coupled from EM system.
Ku (t) × Ipu (Vdie (t)) + Kd (t) × Ipd (Vdie (t)) + Ipc (Vdie (t)) + Igc (Vdie (t)) = Iout (t) (5.29)
At the beginning time point and the end time point of the transition process, the
summation of Ku and Kd is 1. Therefore, assume at every time point it still holds,
we have Eq. (5.30) below.
Ku (t) + Kd (t) = 1
With Eq. (5.29) and Eq. (5.30), Ku and Kd can be solved. Therefore Iout (t) at
each time point can be calculated on the port and then radiate back to EM system.
Finally, as pointed out in [115], the assumption that the summation of Ku and
Kd is 1 at all time steps does not necessarily hold; it is only an approximation. When
only one pair of switching VT tables are given (at GND state the rising edge and at
VCC state the falling edge data), we need the assumption in Eq. (5.30) to solve for
the coefficients, which is some of the IBIS model cases. However when two pairs of
the switching VT tables are available (it includes both in GND state the rising and
falling edge’s and in VCC state the rising and falling edge’s VT relation data), it
means the practical situations where the pullup and pulldown transistors do not start
and end the switching process simultaneously. In that case, the assumption in Eq.
(5.30) does not exist, the independent multiplier Ku and Kd arrays can be solved.
Remarks on IBIS implementation
Here, it is worth noting that, the co-simulation is again conducted in a self-looping
process at each time step. What’s more, since IBIS models is a behavioral model, it
is faster than SPICE circuit solver. However, since it is not a simulator model, I-V
data in IBIS file are supplied over the range of voltages the output could possibly
generate or experience. However, if a buffer is operating in an enrivonment where
its output could be actively driven beyond the limits, the I-V table must be extened
using extropolation. In that scenario, the IBIS model might not be accurate enough
though [112].
Numerical examples
In this section, numerical examples are given to validate the proposed method.
Validation example – High speed links
Here we consider a simple showcase example where the two ends of a microstrip
are connected to the driver (output) and receiver (input) IBIS models, respectively.
The geometry details are listed in Fig. (5.21) below.
At one end of the mircostrip line, it is connected to a driver, the 1Y pin of the
quadruple bus buffer gate (The SN 74AHCT 125, demonstrated in Fig. (5.22). At
the other end it is connected to a receiver, the 1A pin of another quadruple bus beffer
Figure 5.21: Microstrip line with IBIS models.
Figure 5.22: Package top view.
Figure 5.23: Result comparison of voltages at the driver side of the link.
Below in Fig. (5.23) and Fig. (5.24) the voltages calculated at the driver side V 1
and at the receiver side V 2 are plotted respectively. Excellent agreement is found in
CST and DGTD simulations for both 1Y and 1A ends.
A simple board example
Herein a simplified mictrostrip line coupling through vias in a four-layer board with
IBIS models is studied. The geometry of the problem is demonstrated in Fig. (5.25).
From top to bottom, the first layer is the mask with the permittivity constant of 3.0
and the rest three layers are substrates with the permittivity of 4.2, on/through which
the traces and vias are built. The second layer and the bottom layer are grounded
with PEC.
The microstrip is connected with one microcontroller and one memory chip on
board. They are products of Freescale and Micron. The corresponding IBIS files
mpc52001.ibs and t37zi t.ibs can be downloaded from the homepages of Freescale and
Figure 5.24: Result comparison of voltages at the receiver side of the link.
Figure 5.25: Geometry Description of the simple board. (a): top view, (b): side view.
Figure 5.26: Result comparison of voltages at the microcontroller (port1).
Micron, respectively. Here in this example the ouput buffer of the microcontroller
and the input buffer of the memory chip are used.
Since DG method can naturally support non-conformal meshes, it is very suitable
for this type of IC applications with layer-by-layer stackup structure. Each layer
can be meshed indepently. The interfaces between subdomains are allowed to be
geometrically non-conformal. In this way, the mesh generation becomes less intensive.
With the surrounding six air boxes as truncation (not shown here), in total, we end up
to have 10 domains that utilizes 10 MPI process for parallelization. The voltages at
the two ends are recorded. This model is also computed with CST simulation. Again
good agreement is realized. The computational results are shown in Fig. (5.26) and
Fig. (5.27), respectively.
Although this work for the demonstration purpose of the validation and capability
of the DGTD-IBIS co-simulation, simplies the real-life IC applications to a great
Figure 5.27: Result comparison of voltages at the memory chip (port2).
extend, it demonstrates the promise for solving real-life complicated problems with
the versatility and stability of the platform.
In this work, DGTD is proposed for electromagnetic-circuit co-simulation of PCB
problems with the incorporation of IBIS models. The problem is decomposed into
electromagnetic part and circuit part. The electromagnetic part is simulated based
on DGTD. The I-V tables and V-T tables of IBIS models are employed to build the
circuit part that describe the input/output features of IC chips. The coupling between
the EM and circuit parts are addressed in a self-consistent manner with convergent
looping process. This coupling between DGTD and circuit part only involves local
modification of DGTD.
Chapter 6: Conclusion
The complexity of modern electromagnetic applications demand effective and efficient numerical methods. Various numerical methods are out there with strength
emphasizing different aspects. Therefore, to apply the appropriate solver for a specific
problem is of vital importance. Usually different solvers can be applied to solve the
same problem, but the convenience for modeling, the efficiency of the simulation and
even the credibility of the result vary among them.
This work explores multiple solver technologies to solve a wide range of multiscale and multi-physics microwave, RF and high-speed applications. In Chapter 3,
an efficient universal matrice approach to treat continuous varying material property tensor with curved boundary is developed. This eases the modeling process of
real-life problems with spatially changing material properties, as well as provides a
more physical thus more accurate representation of curved FEs. In Chapter 4, a
numerical extraction method of the equivalent permeability tensor of ferromagnetic
nanowire array metamaterial is presented. This method starts with computing the
static magnetic field distribution with the acceralation of 2D-FFT on the wires, then
performing a small signal analysis of the equation of magnetization motion. With the
extracted material property, the utilization of FMNW array in microwave components
and devices is studied with full-wave electromagnetic simulations. Next, coming to
Chapter 5, full-wave transient EM circuit phenomena with linear/non-linear effects
are investigated through DGTD method coupled with SPICE and IBIS circuit models
in a self-consistent integration scheme. Numerical examples all validate the proposed
Research for future development on solving multi-scale and multi-physics problems
may take the following possible paths. One arises in the time-domain simulation for
nonlinear problems. When handling linear problems, the time step size is either
bounded by the element size when using explicit time marching scheme, or by the
Nyquist sampling frequency when using implicit scheme. However, when highly nonlinear phenomena appear, those criterias are not sufficient anymore. It is important
to adaptively choose an appropriate step size at each time stepping of the process, so
that the non-linear effect can be accurately accounted for while it does not cause too
long the CPU time. What is more, spatially at areas where the time step needs to be
chosen very small, the meshes can be grouped for implicit time integration, while the
rest can still proceed in an explicit fashion, thus to develop a hybrid implicit/explicit
time marching method is preferred. Another direction is noticed from the multiphysics simulation. When a real-life system is complicated and invovling different
physics and features, it is an inter-disciplinary subject. Most of the situation, not
all the disciplines count the same weight. How to extract the dominate physics and
make an approximation of the effects from the minor factors that are responsible for
the whole system is therefore important. It can greatly enhance the efficiency of the
computation as well as shed lights for the understanding of the foundation of the
problem itself. In that regards, to develop a versatile computational platform that
can somewhat automatically use the optimal solver for the specific type of problme
is some path we endeavor for.
[1] A. F. Peterson, S. L. Ray, R. Mittra, I. of Electrical, and E. Engineers, Computational methods for electromagnetics, vol. 2. IEEE press New York, 1998.
[2] J.-M. Jin, The finite element method in electromagnetics. Wiley-IEEE Press,
[3] J.-f. Lee, D.-K. Sun, and Z. Cendes, “Full-wave analysis of dielectric waveguides
using tangential vector finite elements,” IEEE Transactions on Microwave Theory and Techniques, vol. 39, no. 8, pp. 1262–1271, 1991.
[4] D.-K. Sun, J. Manges, X. Yuan, and Z. Cendes, “Spurious modes in finite
element methods,” IEEE Antennas and Propagation Magazine, vol. 37, no. 5,
pp. 12–24, 1995.
[5] R. Dyczij-Edlinger, G. Peng, and J.-F. Lee, “A fast vector-potential method
using tangentially continuous vector finite elements,” Microwave Theory and
Techniques, IEEE Transactions on, vol. 46, no. 6, pp. 863–868, 1998.
[6] G. Peng and J.-F. Lee, “Analysis of biaxially anisotropic waveguides using
tangential vector finite elements,” Microwave and Optical Technology Letters,
vol. 9, no. 3, pp. 156–162, 1995.
[7] L. R. Scott and S. Zhang, “Finite element interpolation of nonsmooth functions
satisfying boundary conditions,” Mathematics of Computation, vol. 54, no. 190,
pp. 483–493, 1990.
[8] J.-F. Lee, D. Sun, and Z. Cendes, “Tangential vector finite elements for electromagnetic field computation,” Magnetics, IEEE Transactions on, vol. 27, no. 5,
pp. 4032–4035, 1991.
[9] J. P. Webb, “Hierarchal vector basis functions of arbitrary order for triangular
and tetrahedral finite elements,” Antennas and Propagation, IEEE Transactions on, vol. 47, no. 8, pp. 1244–1253, 1999.
[10] A. Quarteroni and A. Valli, Domain decomposition methods for partial differential equations. No. CMCS-BOOK-2009-019, Oxford University Press, 1999.
[11] B. Smith, P. Bjorstad, and W. Gropp, Domain decomposition: parallel multilevel methods for elliptic partial differential equations. Cambridge university
press, 2004.
[12] K. Zhao, V. Rawat, and J.-F. Lee, “A domain decomposition method with nonconformal meshes for finite periodic and semi-periodic structures,” IEEE Trans.
Ant. Prop., vol. 55, pp. 2559–2570, Sept. 2007.
[13] Y.-J. Li and J.-M. Jin, “A new dual-primal domain decomposition approach for
finite element simulation of 3-D large-scale electromagnetic problems,” IEEE
Trans. Ant. Prop., vol. 55, pp. 2803–2810, Oct. 2007.
[14] Z. Peng, V. Rawat, and J.-F. Lee, “One way domain decomposition method with
second order transmission conditions for solving electromagnetic wave problems,” J. Comput. Phys., vol. 229, pp. 1181–1197, 2010.
[15] K. S. Kunz and R. J. Luebbers, The finite difference time domain method for
electromagnetics. CRC press, 1993.
[16] G. Mur, “Absorbing boundary conditions for the finite-difference approximation
of the time-domain electromagnetic-field equations,” Electromagnetic Compatibility, IEEE Transactions on, no. 4, pp. 377–382, 1981.
[17] S. Dey and R. Mittra, “A locally conformal finite-difference time-domain
(FDTD) algorithm for modeling three-dimensional perfectly conducting objects,” Microwave and Guided Wave Letters, IEEE, vol. 7, no. 9, pp. 273–275,
[18] T. G. Jurgens and A. Taflove, “Three-dimensional contour FDTD modeling of
scattering from single and multiple bodies,” IEEE Transactions on Antennas
and Propagation, vol. 41, no. 12, pp. 1703–1708, 1993.
[19] J.-F. Lee, R. Lee, and A. Cangellaris, “Time-domain finite-element methods,”
IEEE transactions on antennas and propagation, vol. 45, no. 3, pp. 430–442,
[20] S. D. Gedney and U. Navsariwala, “An unconditionally stable finite element
time-domain solution of the vector wave equation,” Microwave and Guided
Wave Letters, IEEE, vol. 5, no. 10, pp. 332–334, 1995.
[21] D. Jiao and J.-M. Jin, “A general approach for the stability analysis of the
time-domain finite-element method for electromagnetic simulations,” Antennas
and Propagation, IEEE Transactions on, vol. 50, no. 11, pp. 1624–1632, 2002.
[22] D. A. White, “Orthogonal vector basis functions for time domain finite element
solution of the vector wave equation,” IEEE Transactions on magnetics, vol. 35,
no. 3, pp. 1458–1461, 1999.
[23] R.-B. Wu and T. Itoh, “Hybrid finite-difference time-domain modeling of curved
surfaces using tetrahedral edge elements,” Antennas and Propagation, IEEE
Transactions on, vol. 45, no. 8, pp. 1302–1309, 1997.
[24] B. Cockburn and C.-W. Shu, “The local discontinuous Galerkin method for
time-dependent convection-diffusion systems,” SIAM Journal on Numerical
Analysis, vol. 35, no. 6, pp. 2440–2463, 1998.
[25] J. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured
grids,” J. Comput. Phys., vol. 181, pp. 186–221, May 2002.
[26] S. D. Gedney, C. Luo, B. Guernsey, J. A. Roden, R. Crawford, J. Miller, et al.,
“The discontinuous Galerkin finite element time domain method (DGFETD):
A high order, globally-explicit method for parallel computation,” in Electromagnetic Compatibility, 2007. EMC 2007. IEEE International Symposium on,
pp. 1–3, IEEE, 2007.
[27] S. Dosopoulos, J. D. Gardiner, and J.-F. Lee, “An MPI/GPU parallelization of
an interior penalty discontinuous Galerkin time domain method for Maxwell’s
equations,” Radio Science, vol. 46, no. 3, 2011.
[28] L. Angulo, J. Alvarez, F. Teixeira, M. Pantoja, and S. Garcia, “Causal-path
local time-stepping in the discontinuous Galerkin method for Maxwell’s equations,” Journal of Computational Physics, vol. 256, pp. 678–695, 2014.
[29] V. Boucher, L.-P. Carignan, T. Kodera, C. Caloz, A. Yelon, and D. Menard,
“Effective permeability tensor and double resonance of interacting bistable ferromagnetic nanowires,” Phys. Rev. B, vol. 80, 2009.
[30] L.-P. Carignan, A. Yelon, D. Menard, and C. Caloz, “Ferromagnetic nanowire
metamaterials: theory and applications,” IEEE Trans. Microw. Theory and
Tech., vol. 59, pp. 2568–2586, 2011.
[31] C. R. Paul, Introduction to electromagnetic compatibility, vol. 184. John Wiley
& Sons, 2006.
[32] E.-P. Li, X.-C. Wei, A. C. Cangellaris, E.-X. Liu, Y.-J. Zhang, M. D’Amore,
J. Kim, and T. Sudo, “Progress review of electromagnetic compatibility analysis technologies for packages, printed circuit boards, and novel interconnects,”
Electromagnetic Compatibility, IEEE Transactions on, vol. 52, no. 2, pp. 248–
265, 2010.
[33] S. S. Zivanovic, K. S. Yee, and K. K. Mei, “A subgridding method for the
time-domain finite-difference method to solve Maxwell’s equations,” Microwave
Theory and Techniques, IEEE Transactions on, vol. 39, no. 3, pp. 471–479,
[34] M. W. Chevalier, R. J. Luebbers, and V. P. Cable, “FDTD local grid with
material traverse,” Antennas and Propagation, IEEE Transactions on, vol. 45,
no. 3, pp. 411–421, 1997.
[35] F. L. Teixeira, “Time-domain finite-difference and finite-element methods for
Maxwell’s equations in complex media,” IEEE Transactions on Antennas and
Propagation, vol. 56, no. 8, pp. 2150–2166, 2008.
[36] H.-P. Tsai, Y. Wang, and T. Itoh, “An unconditionally stable extended (USE)
finite-element time-domain solution of active nonlinear microwave circuits using
perfectly matched layers,” Microwave Theory and Techniques, IEEE Transactions on, vol. 50, no. 10, pp. 2226–2232, 2002.
[37] Z. Lou and J.-M. Jin, “Modeling and simulation of broad-band antennas using the time-domain finite element method,” Antennas and Propagation, IEEE
Transactions on, vol. 53, no. 12, pp. 4099–4110, 2005.
[38] S. Dosopoulos and J.-F. Lee, “Interior penalty discontinuous Galerkin method
for the time-domain Maxwell’s equations,” Magnetics, IEEE Transactions on,
vol. 16, pp. 3512–3515, Aug. 2010.
[39] L. Fezoui, S. Lanteri, S. Lohrengel, and S. Piperno, “Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell’s equations on unstructured meshes,” ESAIM:M2AN, vol. 39,
pp. 1149–1176, Nov. 2005.
[40] V. Rawat, Finite Element Domain Decomposition with Second Order Transmission Conditions for Time-Harmonic Electromagnetic Problems. PhD thesis,
Ohio State University, 2009.
[41] E. Montseny, S. Pernet, X. Ferrières, and G. Cohen, “Dissipative terms and
local time-stepping improvements in a spatial high order discontinuous Galerkin
scheme for the time-domain Maxwell’s equations,” Journal of Computational
Physics, vol. 227, no. 14, pp. 6795–6820, 2008.
[42] D. Cioranescu and P. Donato, “Introduction to homogenization,” 2000.
[43] A. Bensoussan, J.-L. Lions, and G. Papanicolaou, Asymptotic analysis for periodic structures, vol. 374. American Mathematical Soc., 2011.
[44] E. Weinan, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden, “Heterogeneous multiscale methods: a review,” Commun. Comput. Phys, vol. 2, no. 3,
pp. 367–450, 2007.
[45] A. Abdulle, “The finite element heterogeneous multiscale method: a computational strategy for multiscale pdes,” GAKUTO International Series Mathematical Sciences and Applications, vol. 31, no. EPFL-ARTICLE-182121, pp. 135–
184, 2009.
[46] D. R. Smith, S. Schultz, P. Markos, and C. M. Soukoulis, “Determination of
effective permittivity and permeability of metamaterials from reflection and
transmission coefficients,” Phys. Rev. B, vol. 65, no. 19, 2002.
[47] D. R. Smith and J. P. Pendry, “Homogenization of metamaterials by field averaging,” J. Opt. Soc. Am. B, vol. 23, pp. 391–403, 2006.
[48] J. Parry, C. Bailey, and C. Aldham, “Multiphysics modelling for electronics
design,” in Thermal and Thermomechanical Phenomena in Electronic Systems,
2000. ITHERM 2000. The Seventh Intersociety Conference on, vol. 2, pp. 86–
93, IEEE, 2000.
[49] D. Langoni, M. H. Weatherspoon, E. Ogunti, and S. Y. Foo, “An overview of
reconfigurable antennas: Design, simulation, and optimization,” in Wireless and
Microwave Technology Conference, 2009. WAMICON’09. IEEE 10th Annual,
pp. 1–5, IEEE, 2009.
[50] Y. Shao, Z. Peng, and J.-F. Lee, “Thermal-aware DC IR-drop co-analysis using
non-conformal domain decomposition methods,” in Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 468,
pp. 1652–1675, The Royal Society, 2012.
[51] H. Toshiyoshi, “A SPICE-based multi-physics simulation technique for integrated MEMS,” in International Conference on Simulation of Semiconductor
Processes and Devices, 2011.
[52] B. Zhao, J. C. Young, and S. D. Gedney, “SPICE lumped circuit subcell model
for the discontinuous Galerkin finite-element time-domain method,” Microwave
Theory and Techniques, IEEE Transactions on, vol. 60, pp. 2684–2692, Sept.
[53] E. Martini, G. Pelosi, and S. Selleri, “With a new type of curvilinear mapping
for the analysis of microwave passive devices,” IEEE Trans. Ant. Prop., vol. 51,
pp. 1712–1717, June 2003.
[54] J. S. Wang and N. Ida, “Curvilinear and higher order edge finite elements in
electromagnetic field computation,” IEEE Trans. Magn., vol. 29, pp. 1491–1494,
Mar. 1993.
[55] M. M. Ilic, A. Z. Ilic, and B. M. Notaros, “Continuously inhomogeneous higher
order finite elements for 3-d electromagnetic analysis,” IEEE Trans. Ant. Prop.,
vol. 57, pp. 2798–2803, Sept. 2009.
[56] J. P. Webb, “Universal matrices for high order finite elements in nonlinear
magnetic field problems,” IEEE Trans. Magn., vol. 33, Sept. 1997.
[57] M. Dufresne and P. P. Silvester, “Universal matrices for the N-dimensional finite
element,” Computation in Electromagnetics, Third International Conference on
(Conf. Publ. No. 420), pp. 223–228, 1996.
[58] F. Teixeira and W. Chew, “Analytical derivation of a conformal perfectly
matched absorber for electromagnetic waves,” Microw. Opt. Tech. Lett., vol. 17,
pp. 231–236, 1997.
[59] J. Berenger, “Three-Dimensional perfectly matched layer for the absorbtion of
electromagnetic waves,” J. Comput. Phys., vol. 127, no. 2, pp. 363–379, 1996.
[60] J. C. Nedelec, “Mixed finite elements in R3 .,” Numer. Math., vol. 35, pp. 315–
341, 1980.
[61] D.-K. Sun, J.-F. Lee, and Z. Cendes, “Construction of nearly orthogonal Nedelec
bases for rapid convergence with multilevel preconditioned solvers,” SIAM J.
Sci. Comput., vol. 23, no. 4, pp. 1053–1076, 2001.
[62] D. R. Wilton, R. D. Graglia, and A. Peterson, “Higher order interpolatory vector bases for computational electromagnetics,” IEEE Trans. Ant. Prop., vol. 45,
no. 3, pp. 329–342, 1997.
[63] D. Ansari Oghol Beig and M. S. Leong, “Dual-grid based tree/cotree decomposition of higher-order interpolatory H(∇∧, ω) basis,” Int. J. Numer. Meth.
Engng., 2010.
[64] P. Solin and K. Segeth, Higher-Order Finite Element Methods. Chapman &
Hall, 2004.
[65] R. Sevilla, S. Fernández-Méndez, and A. Huerta, “NURBS-enhanced finite element method (NEFEM),” International Journal for Numerical Methods in
Engineering, vol. 76, pp. 56–83, Feb. 2008.
[66] A. Perronnet, “Interpolation transfinie sur le triangle, le tétraèdre et le
pentaèdre. Application à la création de maillages et à la condition de Dirichlet,”
Comptes Rendus de l’Académie des Sciences - Series I - Mathematics, vol. 326,
no. 1, pp. 117–122, 1998.
[67] F. Teixeira and W. Chew, “Differential forms, metrics, and the reflectionless
absorption of electromagnetic waves,” Journal of Electromagnetic Waves and
Applications, vol. 13, no. 5, pp. 665–686, 1999.
[68] F. Teixeira and W. Chew, “Complex space approach to perfectly matched layers: a review and some new developments,” International Journal of Numerical
Modelling: Electronic Networks, Devices and Fields, vol. 13, no. 5, pp. 441–455,
[69] F. Teixeira, K.-P. Hwang, W. Chew, and J. Jin, “Conformal PML-FDTD
schemes for electromagnetic field simulations: A dynamic stability study,” Antennas and Propagation, IEEE Transactions on, vol. 49, no. 6, pp. 902–907,
[70] B. Donderici and F. L. Teixeira, “Conformal perfectly matched layer for the
mixed finite element time-domain method,” Antennas and Propagation, IEEE
Transactions on, vol. 56, no. 4, pp. 1017–1026, 2008.
[71] F. Teixeira and W. Chew, “Systematic derivation of anisotropic PML absorbing
media in cylindrical and spherical coordinates,” Microwave and Guided Wave
Letters, IEEE, vol. 7, no. 11, pp. 371–373, 1997.
[72] B. Fuchs, O. Lafond, S. Palud, L. L. Coq, M. Himdi, M. C. Buck, and
S. Rondineau, “Comparative design and analysis of luneburg and half Maxwell
fish-eye lens antennas,” IEEE Trans. Ant. Prop., vol. 56, pp. 3058–3062, Sept.
[73] Hossein Mosallaei and Y. Rahmat-Samii, “Nonuniform Luneburg and two-shell
lens antennas: radiation characteristics and design optimization,” IEEE Trans.
Ant. Prop., vol. 49, pp. 60–69, Jan. 2001.
[74] V. G. Veselago, “The electrodynamics of substances with simultaneously negative values of ϵ and µ,” Sov. Phys. Usp., vol. 10, pp. 509–514, 1968.
[75] D. R. Smith, W. Padilla, D. C. Vier, S. C. Nemat-Nasser, and S. Schultz, “A
composite medium with simultaneously negative permeability and permittivity,” Phys. Rev. Usp., vol. 84, pp. 4184–4187, 2000.
[76] C. Caloz and T. Itoh, “Application of the transmission line theory of left-handed
(LH) materials to the realization of a microstrip LH line,” IEEE AP-S Int.
Symp., vol. 2, pp. 412–415, 2002.
[77] J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic fields,”
Science, vol. 312, 2006.
[78] A. Alu, “First-principle homogenization theory for periodic metamaterials,”
Phys. Rev. B, vol. 84, 2011.
[79] I. Tsukerman, “Effective parameters of metamaterials: a rigorous homogenization theory via Whitney interpolation,” JOSA B, vol. 28, pp. 577–586, 2011.
[80] S. M. Seo and J.-F. Lee, “A Fast IE-FFT algorithm for solving PEC scattering
problems,” IEEE Trans. Magn., vol. 41, pp. 1476–1479, 2005.
[81] J. Ramprecht and D. Sjoberg, “Biased magnetic materials in RAM applications,” Progress in Electromagnetics Research, vol. 75, pp. 85–117, 2007.
[82] G. Hamoir, J. De La Torre Medina, L. Piraux, and I. Huynen, “Self-biased
nonreciprocal microstrip phase shifter on magnetic nanowired substrate suitable for gyrator applications,” IEEE Trans. Microw. Theory and Tech., vol. 60,
pp. 2152–2157, 2012.
[83] L.-P. Carignan, V. Boucher, T. Kodera, C. Caloz, A. Yelon, and D. Menard,
“Double ferromagnetic resonance in nanowire arrays,” Applied Physics Letters,
vol. 95, no. 6, 2009.
[84] L.-P. Carignan, T. Kodera, A. Yelon, C. Caloz, and D. Menard, “Integrated and
self-biased planar magnetic microwave circuits based on ferromagnetic nanowire
substrates,” IEEE EuMC, pp. 743–746, 2009.
[85] S. H. Hall and H. L. Heck, Advanced Signal Integrity for High-Speed Digital
Design. Hoboken, NJ: Wiley, 2009.
[86] International Technology Roadmap for Semiconductors. http://www.itrs.
net, 2013.
[87] Piket-May, Melinda, Taflove, Allen, Baron, and John, “FD-TD modeling of digital signal propagation in 3-D circuits with passive and active loads,” Microwave
Theory and Techniques, IEEE Transactions on, vol. 42, pp. 1514–1523, Aug.
[88] V. A. Thomas, M. E. Jones, M. Piket-May, A. Taflove, and E. Harrigan, “The
use of SPICE lumped circuits as sub-grid models for FDTD analysis,” IEEE
microwave and guided wave letters, vol. 4, pp. 141–144, May 1994.
[89] C.-N. Kuo, B. Houshmand, and T. Itoh, “Full-wave analysis of packaged microwave circuits with active and nonlinear devices: An FDTD approach,” Microwave Theory and Techniques, IEEE Transactions on, vol. 45, pp. 819–826,
May 1997.
[90] A. Witzig, C. Schuster, P. Regli, and W. Fichtner, “Global modeling of microwave applications by combining the FDTD method and a general semiconductor device and circuit simulator,” Microwave Theory and Techniques, IEEE
Transactions on, vol. 47, pp. 919–928, June 1999.
[91] C. Guo and T. H. Hubing, “Circuit models for power bus structures on printed
circuit boards using a hybrid FEM-SPICE method,” Advanced Packaging, IEEE
Transactions on, vol. 29, pp. 441–447, Aug. 2006.
[92] K. Guillouard, M. F. Wong, V. F. Hanna, and J. Citerne, “A new global timedomain electromagnetic simulator of microwave circuits including lumped elements based on finite-element method,” Microwave Theory and Techniques,
IEEE Transactions on, vol. 47, pp. 2045–2048, Oct. 1999.
[93] R. Wang and J.-M. Jin, “A symmetric electromagnetic-circuit simulator based
on the extended time-domain finite element method,” Microwave Theory and
Techniques, IEEE Transactions on, vol. 56, pp. 2875–2884, Dec. 2008.
[94] S. Dosopoulos and J.-F. Lee, “Interconnect and lumped elements modeling
in interior penalty discontinuous Galerkin time-domain methods,” J. Comput.
Phys., vol. 229, pp. 8521–8536, Aug. 2010.
[95] S. Dosopoulos and J.-F. Lee, “Interior penalty discontinuous Galerkin finite
element method for the time-dependent first order Maxwell’s equations,” Antennas and Propagation, IEEE Transactions on, vol. 58, pp. 4085–4090, Dec.
[96] P. Li and L. J. Jiang, “A hybrid electromagnetics-circuit simulation method exploiting discontinuous Galerkin finite element time domain method,” Microwave
and Wireless Components Letters, IEEE, vol. 23, pp. 113–115, Mar. 2013.
[97] P. Li and L. J. Jiang, “Integration of arbitrary lumped multiport circuit networks into the discontinuous Galerkin time-domain analysis,” Microwave Theory and Techniques, IEEE Transactions on, vol. 61, pp. 2525–2534, July 2013.
[98] X. Gu, R. Rimolo-Donadio, Z. Yu, F. de Paulis, Y. H. Kwark, M. Cocchini,
M. B. Ritter, B. Archambeault, A. Ruehli, J. Fan, and C. Schuster, “Fast
physics-based via and trace models for signal and power integrity co-analysis,”
in Proc. DesignCon Conference, Feb. 2010.
[99] D. N. Arnold, “An interior penalty finite element method with discontinuous
elements,” SIAM J. Numer. Analy., vol. 19, pp. 742–760, Aug. 1982.
[100] P. Houston, I. Perugia, A. Schneebeli, and D. Schotzau, “Interior penalty
method for the indefinite time-harmonic Maxwell’s equations,” Numer. Math.,
vol. 100, p. 485518, Mar. 2005.
[101] E. Montseny, S. Pernet, X. Ferrires, and G. Cohen, “Dissipative terms and
local time-stepping improvements in a spatial high order discontinuous Galerkin
scheme for the time-domain Maxwell’s equations,” J. Comput. Phys., vol. 227,
pp. 6795–6820, July 2008.
[102] S. Dosopoulos and J.-F. Lee, “Non-conformal and parallel discontinuous
Galerkin time domain method for Maxwell’s equations: EM analysis of IC packages,” J. Comput. Phys., vol. 238, pp. 48–70, Dec. 2012.
[103] C.-W. Ho, A. Ruehli, and P. Brennan, “The modified nodal approach to network
analysis,” Circuits and Systems, IEEE Transactions on, vol. 22, pp. 504–509,
June 1975.
[104] T. L. Quarles, The SPICE3 implementation guide. Univ. California at Berkeley,
Berkeley, CA, Tech. Rep., 1989.
[105] Hspice User Guide,Synopsys, Inc. ., 2008.
[106] Ngspice Users Manual., 2014.
[107] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical
Recipes 3rd Edition: The Art of Scientific Computing. New York, NY, USA:
Cambridge University Press, 3 ed., 2007.
[108] S.-J. Moon and A. Cangellaris, “SPICE-compatible representations of SParameter matrices of passive networks with transport delay,” in Proc. Electr.
Perform. Electron. Packag., Oct. 2007.
[109] B. Gustavsen and A. Semlyen, “Rational approximation of frequency domain
responses by vector fitting,” Power Delivery, IEEE Transactions on, vol. 14,
pp. 1052–1061, July 1999.
[110] D. Saraswat, R. Achar, and M. S. Nakhla, “A fast algorithm and practical considerations for passive macromodeling of measured/simulated data,” Advanced
Packaging, IEEE Transactions on, vol. 27, pp. 1052–1061, Feb. 2004.
[111] Computer Simulation Technology, CST STUDIO SUITE. http://www.cst.
net, 2013.
[112] S. Peters, “IBIS forum I/O buffer modeling cookbook,” Revision 2, vol. 2, 1997.
[113] B. Baker, “The IBIS model: A conduit into signal-integrity analysis, part 1,”
Analog Applications Journal Q, vol. 4, p. 2010.
[114] M. Casamayor, “A first approach to IBIS models: What they are and how they
are generated,” AN-715Application Note, Analog Device, 2004.
[115] Y. Wang and H. N. Tan, “The development of analog SPICE behavioral model
based on IBIS model,” in VLSI, 1999. Proceedings. Ninth Great Lakes Symposium on, pp. 101–104, IEEE, 1999.
Без категории
Размер файла
11 173 Кб
Пожаловаться на содержимое документа