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

код для вставкиСкачатьEM MODELING AND SIMULATION OF MICROWAVE ELECTRONIC COMPONENTS AND DEVICES WITH MULTI-SCALE AND MULTI-PHYSICS EFFECTS Dissertation Presented in Partial Fulﬁllment of the Requirements for the Degree Doctor of Philosophy in the Graduate School of The Ohio State University By Jue Wang, B.Eng Graduate Program in Electrical and Computer Engineering The Ohio State University 2015 Dissertation Committee: Prof. Jin-Fa Lee, Advisor Prof. Kubilay Sertel Prof. Fernando L. Teixeira ProQuest Number: 10001908 All rights reserved INFORMATION TO ALL USERS 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 2015 Abstract This work investigates various numerical methods for modeling and analyzing microwave components and electronic devices with multi-scale structures and multiphysics eﬀects. In the ﬁrst 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 eﬀects in high-frequency circuit systems is explored. To be speciﬁc, 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 ii 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. iii This work is dedicated to the ones that I love. iv Acknowledgments 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 services. I am particularly grateful to Dr. Zhen Peng. His wide knowledge has given me signiﬁcant 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 Mahaﬀey, Wan-Chun Liao, Cagdas Gunes, Dr. Dongwei Li, Jiaqing Lu, Shaoxin Peng, Xuezhe Tian and many more others. v 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 ﬁrst-class company industrial atmosphere, I received a lot of unselﬁsh help and friendship from diﬀerent 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 ﬁnd 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. vi Vita 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. Publications Research Publications D. Ansari, J. Wang, Z. Peng, J.-F Lee “A universal array approach for ﬁnite 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 vii Table of Contents Page Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Vita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1. Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 1.2 1.3 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . Dissertation Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 1 4 8 Literature Review on Numerical Methods in Multi-scale and Multi-physics Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2. 2.1 2.2 2.3 2.4 2.5 3. Overview . . . . . . . . . . . . . FDTD and FETD methods . . . DGTD methods . . . . . . . . . Homogenization of Heterogeneous Multi-physics coupling . . . . . . . . . . . 9 10 11 12 13 A Universal Array Approach for Finite Elements with Continuously Material Properties and Curved Boundaries . . . . . . . . . . . . . . . . . . 14 3.1 15 Introduction . . . . . . . . . . . . . . . . . . . . . Composites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii 3.2 3.3 3.4 3.5 3.6 3.7 4. 5. . . . . . . 17 18 19 30 31 39 Ferromagnetic Nano-wires: Homogenization and Applications . . . . . . 42 4.1 4.2 4.3 4.4 42 44 52 64 . . . . . . . . . . . . . . . . . . . . . . . . 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 5.1 5.2 6. Methodology . . . . . . . . . . Notation . . . . . . . . . . . . . Universal Arrays . . . . . . . . CIMPT Approach for the PML Numerical Results . . . . . . . Conclusion . . . . . . . . . . . DGTD integrated with SPICE solver . . . . . . . . . . . . . . . . . 69 DGTD integrated with IBIS model . . . . . . . . . . . . . . . . . . 101 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 ix List of Tables Table Page 3.1 Conformal-DDM Related Notation . . . . . . . . . . . . . . . . . . . 27 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 ﬁeld 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 ﬁeld vs. Normalized Magnetization . . . . . . 56 5.1 Geometrical description of the pulse generator . . . . . . . . . . . . . 97 5.2 Step-recovery diode description . . . . . . . . . . . . . . . . . . . . . 99 x List of Figures Figure Page 3.1 (a) The 3D reference element. (b) An actual curved element. . . . . 20 3.2 (a) Edge representation of the curved element. (b) Face representation of the curved element. Only two faces are plotted here. . . . . . . . 28 3.3 The cubic Lagrangian interpolation of the curved tetrahedron. . . . 28 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. . . . 29 3.5 (a) Problem geometry. (b) |S11 | versus frequency. . . . . . . . . . . . 33 3.6 (a)Electrical ﬁeld intensity. (b) Conformal-DDM’s Partition of the Mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.7 Schematic diagram showing the conﬁguration of the quarter-Luneburg lens scattering problem. Curved elements were applied to the surface of the lens only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.8 Schematic diagram showing the conﬁguration of the R = 6λ0 Luneburg lens antenna problem. Curved elements were applied to the surface of the lens only. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.9 A three-slice ﬁeld plot of the E · ŷ ﬁeld 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. . . . . . . . . . 37 3.10 The radiation pattern (dB) of Luneburg lens antennas of various radius R as function of θ at ϕ = 0. . . . . . . . . . . . . . . . . . . . . . . . 37 xi 3.11 Schematic diagram showing the conﬁguration of the PML-lens antenna problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.12 (a) Snapshots of the magnitude of the electric ﬁeld distribution on yzplane. (b) Snapshots of the magnitude of the electric ﬁeld distribution on zx-plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.13 The radiation pattern (dB) of the Luneburg lens antenna with R = λ0 at ϕ = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.1 Sketch of FMNW substrate structure. . . . . . . . . . . . . . . . . . . 45 4.2 Schematic hysteresis loop. . . . . . . . . . . . . . . . . . . . . . . . . 46 4.3 Flow-chart of FMNW based metamaterial homogenization process and application enumeration. . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.4 Comparison of eﬀective relative permeability tensor of CoFe nanowire membrane obtained by diﬀerent methods. . . . . . . . . . . . . . . . 54 4.5 Comparison of eﬀective relative permeability tensor of NiFe nanowire membrane obtained by diﬀerent methods. . . . . . . . . . . . . . . . 54 4.6 Microstrip line sketch. . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.7 Major hysteresis loop (solid blue) and minor hysteresis curves (dashed red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.8 Left: contour plot of the S21 parameter from simulated result. Right (a) and (b): reference experiment and model results. . . . . . . . . . 58 4.9 Electric ﬁeld distributions along the microstrip line with FMNW subext strate and biased with the same external bias ﬁeld, Hdc = 0.5KOe, at two diﬀerent resonance frequencies. . . . . . . . . . . . . . . . . . . . 59 4.10 FMNW isolator structure. . . . . . . . . . . . . . . . . . . . . . . . . 60 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. 61 xii 4.12 S12 parameter comparison. . . . . . . . . . . . . . . . . . . . . . . . . 62 4.13 S21 parameter comparison. . . . . . . . . . . . . . . . . . . . . . . . . 63 4.14 Electric ﬁeld distributions along the microstrip line of the isolator at two diﬀerent working frequencies. . . . . . . . . . . . . . . . . . . . . 64 4.15 Detailed description of a dual-band self-biased circulator. . . . . . . . 65 4.16 Full-wave simulation result of the conceptual microstrip circulator–S parameter plotting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.17 Full-wave simulation result of the conceptual microstrip circulator– Electric ﬁeld plotting. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.1 Sketch of EM-Circuit System . . . . . . . . . . . . . . . . . . . . . . 78 5.2 (a) EM coupled to SPICE (b) SPICE coupled to EM . . . . . . . . . 79 5.3 Geometrical illustration of the impedance surface. . . . . . . . . . . . 81 5.4 SPICE netlist example. . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.5 Flowchart of DGTD integrated with SPICE . . . . . . . . . . . . . . 87 5.6 Integration by Gaussian quadrature rule . . . . . . . . . . . . . . . . 89 5.7 Microstrip line with circuit network. . . . . . . . . . . . . . . . . . . . 91 5.8 CST multi-line port experiment. . . . . . . . . . . . . . . . . . . . . . 92 5.9 Result comparison at near end of microstrip. . . . . . . . . . . . . . . 93 5.10 Result comparison at far end of microstrip. . . . . . . . . . . . . . . . 94 5.11 Microstrip line with nonlinear circuit network. . . . . . . . . . . . . . 94 5.12 Comparisons of numerical results with diﬀerent schemes. . . . . . . . 95 5.13 Schematic structure of SRD pulse generator . . . . . . . . . . . . . . . 97 xiii 5.14 Circuit structure of SRD pulse generator. . . . . . . . . . . . . . . . . 98 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 buﬀer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.18 Output buﬀer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 xiv Chapter 1: Background and Motivation 1.1 Background Computational electromagnetics (CEM) is the scientiﬁc subject of modeling the interaction of electromagnetic ﬁelds 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 reﬁne 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. 1 How to accurately and eﬃciently model electromagnetic ﬁelds 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 diﬀerential equation (PDE) based and integral equation (IE) based methods. Those methods are further distinguished as frequency and time domain methods. Among numerous methods, ﬁnite element methods (FEM) [2] have been developed intensely and have achieved huge success in many diﬀerent areas. They solve the Maxwell’s equations starting by discretizing the computational domain into ﬁnite elements, then deﬁne a set of linear basis functions on each element, and approximate the solution in a function space with ﬁnite dimension. FEM formulates the problem from governing PDE, which gives great ﬂexibility in dealing with complex materials. What’s more, the grids can be irregular, it provides a way of more accurate and ﬂexible 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 ﬁeld 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, reﬁning mesh is one way, using higher order basis functions is another way. In [8], the construction of higher-order tangential vector ﬁnite 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 ﬁeld. 2 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 diﬀerence 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 suﬀer from ineﬃciency 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 diﬃcult to maintain, or in other situations the accuracy is sacriﬁced. 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 3 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 eﬃciency 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 diﬀerent groups of mesh with diﬀerent densities can be chosen diﬀerently, instead of being bounded by the smallest mesh size globally [28]. Furthermore, DGTD method can be hybridized with many other solvers with minor modiﬁcation, 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 systems. 1.2 Problem Statement Although numerous CEM methods have achieved great performance in the past, it is still diﬃcult to achieve both certain level of accuracy and low memory cost at 4 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 eﬀective 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 diﬀerent phenomena. That will aggravate the diversity of the system eigen spectrum, causing an ill-conditioned system. Another diﬃculty 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 ﬂowing through conductors. Therefore, it is essential to account for both electrical and thermal eﬀects 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 ﬂuctuations are usually modeled using piecewise 5 constant functions within each element. However, using an appropriate polynomial functions to replace the piecewise representation yields obvious beneﬁts. 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-proﬁle microwave devices with dual-band magnetic eﬀect. In previous work [29] [30], the authors present a Maxwell-Garnett like formulation to analytically extract the eﬀective permeability tensor of FMNWs, while this research work uses a 3D full-wave numerical process more generally to extract the macroscopic eﬀective 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 computations. 6 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 buﬀer information speciﬁcation) 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 eﬃcient since it only involves local modiﬁcation of the EM solver. 7 1.3 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 directions. 8 Chapter 2: Literature Review on Numerical Methods in Multi-scale and Multi-physics Applications 2.1 Overview 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 diﬀerent physical scales and the physics changes sigiﬁcantly over these length-scales. Indeed, solving electromagnetics problem in this enrivonment is a challenging task [31] [32]. To tackle this type of problems eﬃciently, it is essential to distinguish the features of very detailed electrical structure from the rest of the computational domain, which is on a diﬀerent frequency scale. Generally speaking, two broad directions can be taken. One is to locally reﬁne the mesh in areas where ﬁne features occur and ﬁeld 9 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 ﬁelds near the ﬁne feature and then coupled to the global solution elsewhere. In this situation, the diﬃculties come from both obtaining local solutions and coupling process. 2.2 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 ﬁne grid with a small step size only around ﬁne features. An interpolation in space and time is used to calculate the ﬁeld 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 ﬂexibility in multi-scale modeling. Although FDTD method and its variations have achieved huge success during the early attempting, they suﬀer from the inﬂexibility 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 ﬁnite elements, however, compared to FDTD, it provides substantial ﬂexibility 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 10 (PML) is presented. A time-domain ﬁnite element method for modeling complex broad-band antennas is presented in [37]. 2.3 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 ﬂuxes and leap-frog marching in time. A DGTD method with upwind ﬂuxes using a low storage high-order Runge-Kutta marching scheme is developed in [25]. Compared to central ﬂux scheme, upwind ﬂux, 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 diﬀerential equations (ODEs) is obtained. When using an explicit time marching style, the system is conditionally stable; the time step needs to be conﬁned 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 diﬀerent time step sizes can be applied to diﬀerent groups of elements with diﬀerent mesh densities. In this fashion, the computations can be largely speed 11 up. With all the beneﬁts, it is noted that DGTD methods, with the introduction of diﬀerent sets of unknowns on the element interface, suﬀer from a larger number of degrees of freedom compared to FETD methods. 2.4 Homogenization of Heterogeneous Composites When studying the behavior of electromagnetic ﬁeld in a material presenting heterogeneous microstructures, direct computations can often be very diﬃcult due to the fast oscillation of ﬁelds on the very ﬁne 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 ﬁne features are obvious: ﬁrst, 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 reﬂection and transmission coeﬃcients on ﬁnite lengths of the electromagnetic metamaterials. In [47], the equivalent material properties is extracted based on averaging the local ﬁelds. 12 2.5 Multi-physics coupling Apart from the multi-scale features that widely exist in modern EM RF problems, multi-physics eﬀects 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 diﬀerent physics (that matter) couple to each other, based on a common data model is of practical impact [48]. To list only a few eﬀorts already made here: in [49], the authors give an overview of reconﬁgurable 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. 13 Chapter 3: A Universal Array Approach for Finite Elements with Continuously Material Properties and Curved Boundaries 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 ﬂexibility and eﬃciency 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 beneﬁts: a) more accurate physical models for problems with CIMPTs and curved features; 14 and, b) better utilization of computer resources when higher-order curved elements replace lower-order elements of smaller physical dimensions. 3.1 Introduction 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 ﬂexible 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 ﬁsh-eye lens, semiconductor devices and media where material properties are aﬀected or governed by diﬀusion 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 15 physics. In cases as such, material properties follow continuous spatial ﬂuctuations 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 simpliﬁes 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 brieﬂy discussed followed by numerical results. 16 3.2 Methodology Systematically, the evaluation process of FE stiﬀness and mass matrices can be cast into the following basic steps: 1. Assume a reference element with a ﬁxed 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 deﬁned between the reference element and a presumed physical element. 4. Express the required integral-diﬀerential 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 deﬁned (depend) on the reference element. 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. 17 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 veriﬁcation, is used for the solution of a few example problems including an ungraded Luneburg lens. 3.3 Notation 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 T 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 deﬁned 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 deﬁned 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 suﬃciently smooth real function ϱ deﬁned on either the reference or the physical element the column vector ∂ϱ ∂ϱ form of the gradient is denoted as [∇]x ϱ , [ ∂x ] or [∇]ξ ϱ , [ ∂ξ ]. Furthermore, various i i physical quantities and functional spaces are used in this work for which a brief listing can be found in Table. (3.1). 18 3.4 Universal Arrays This section presents a step-by-step realization of the concepts introduced. 3.4.1 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 analytically. 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 deﬁned 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 deﬁnition. [ ] [ ] [ ] [ ] ∂xi ∂ xi◦z(ξ) ∂ξi ∂ ξi◦z−1 (x) J= = ↔ K= = ∂ξj ∂ξj ∂xj ∂xj (3.1) Straightforward application of the chain rule reveals that the following relations hold between column vector representation of the gradients deﬁned over Kr and Kp . [∇]x ϱ = KT × [∇ξ ]ϱ (3.2) [∇]ξ ϱ = JT × [∇x ]ϱ (3.3) 19 Figure 3.1: (a) The 3D reference element. (b) An actual curved element. 3.4.2 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ξ = ∑ ξ βu,j ∇ξ ξj (3.4) j where the u subscript signiﬁes 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 diﬃcult 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) 20 satisfy the necessary tangential continuity conditions for a H(curl, Kr ) conforming basis on Kp [64]. β⃗ux = ∑ ξ βu,j ∇x ξj (3.5) j From Eq. (3.5), one writes: x ]i [βu,i = ∑ [ ∂ξj ] j ∂xi [ ] [ x] ξ ξ βu,j → βu,i = KT × βu,j (3.6) i 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 ]. 3.4.3 The Mass Matrix Using the electrical ﬁeld formulation, the entries of the FE mass matrix T for an electromagnetic radiation/propagation problem can be formulated as ˆ Tuv = ˆ Kp β⃗ux · (ε̄¯ · β⃗vx )dKp = Kp x T ¯ x [βu,i ] × ε̄ ×[βv,j ]dKp (3.7) x 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 as: ˆ Tuv = Kr ξ T ξ [βu,i ] ×K× ε̄¯×KT [βv,j ] det J dKr (3.8) Next, if we deﬁne [ςij ] as Eq. (3.9), T is obtained as Eq. (3.10). [ςij ] , (det J) K× ε̄¯×KT ∑ˆ ξ ξ Tuv = [βu,i ]i [ςij ]ij [βv,j ]j dKr i,j Kr 21 (3.9) (3.10) 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 (ξ)}. Namely, [ςij ](x) ≈ ∑ [ςij ]|x=xs ℓs (ξ ◦ z−1 (x)) (3.11) s 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 deﬁned 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 deﬁne ςijs to represent as [ςij ]|x=xs and have: ( ) ςijs , [ςij ]ij |x=xs = (det J) K× ε̄¯×KT x=xs (3.12) Hence, deﬁning κ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 ]j Tuv = ∑ ∑ (3.13) ˆ ςijs 0≤s<S i,j 22 Kr κui,vj ℓs (ξ) dKr (3.14) 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 = ˆ (ςijs±ςjis) i≤j,0≤s<S ◦ ±uivjs T ˆ , Kr κui,vj ±κuj,vi ℓs dKr 1 + δij Kr (3.15) κui,vj ±κuj,vi ℓs dKr 2 + 2δij (3.16) C± ijs , ςijs ±ςjis (3.17) ∑ 1 (Tuv ±Tvu ) = C± ijs 2 i≤j,0≤s<S ◦ ±uivjs T Implementation of the above formulation leads to total (3.18) 9×(6+3) S 2 multiply-add operations per matrix entry as it seeks to ﬁnd 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. 3.4.4 The Stiffness Matrix Similar to 3.4.3, the electrical ﬁeld formulation for electromagnetic radiation/propagation problem yields the following expression for the FE stiﬀness matrix S: ˆ Suv = Kp ¯−1 ∇x ∧ β⃗vx )dKp ∇x ∧ β⃗ux · (µ̄ 23 (3.19) ¯ 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 suﬃciently smooth function ϕ: ∑ ξ ξ ∇x ∧ β⃗vx= ∇x [βv,m ]m ∧∇x ξm +[βv,m ]m ∇x ∧∇x ξm = ∑ m ξ ]m ∧ ∇x ξm ∇ [βv,m x (3.20) m Knowing that ∇x f = ∑ ∂f x m ∂ξm ∇ ξm ∇x ∧ β⃗vx = we get: ξ ∑ ∂[βv,n ]n m,n ∂ξm ∇x ξm ∧∇x ξn Since ∇x ξm ∧∇x ξn = −∇x ξn ∧∇x ξm , it further develops that: ( ) ξ ξ ∑ ∂[βv,n ∂[β ] ] n m v,m ∇x ∧ β⃗vx = − ∇x ξm∧∇x ξn ∂ξ ∂ξ m n m<n (3.21) (3.22) In Eq. (3.22) the summation over m, n is replaced with m < n which is due to the ¯−1 anti-symmetric nature of the exterior product in ∇x ξm ∧∇x ξn . Next, given the µ̄ ∑ ¯−1 = p,q x̂p µ̄ ¯−1 ¯−1 x ⃗ x tensor as µ̄ p,q x̂q we expand µ̄ ·∇ ∧ βv as: ¯−1 ·(∇x ∧ β⃗vx ) = µ̄ ) ( ξ ξ ∑ ∑ ∂[βv,n ]n ]n ∂[βv,n −1 ¯pq − (x̂q ·∇x ξm∧∇x ξn) x̂p µ̄ ∂ξm ∂ξm p,q m<n (3.23) ¯−1 ·∇x ∧ β⃗vx , κuij,vmn and Now following Eq. (3.19), in order to evaluate ∇x ∧ β⃗ux · µ̄ ςijp,mnq are deﬁned 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 ]m κuij,vmn, − − ∂ξi ∂ξj ∂ξm ∂ξn 24 (3.24) (3.25) Hence, it is important to observe the following symmetries for κuij,vmn and ςijp,mnq : κuij,vmn = κvmn,uij (3.26) ςijp,mnq = ςmnq,ijp (3.27) ςijp,mnq κuij,vmn = ςmnq,ijp κvmn,uij (3.28) 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 ¯−1 µ̄ pq ςijp,mnq (κuij,vmn ±κvij,umn ) det J dKr (3.29) Kr ¯−1 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 reﬂected 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 signiﬁes evaluation at x = xs . In other words we have: ¯−1 ¯−1 µ̄ pqs , µ̄pq (x)|x=xs (3.30) ςijp,mnq,s , ςijp,mnq (x)|x=xs (3.31) det Js , det J|x=xs (3.32) Hence, Eq. (3.29) is reformulated as: Suv ±Svu = ∑ i < j, p m < n, q 0≤s<S 25 ˆ ¯−1 µ̄ pqs ςijp,mnq,s det Js ℓs (ξ)(κuij,vmn ±κvij,umn )dKr (3.33) Kr Now, taking advantage of the (ij) ↔ (mn) symmetry, Eq. (3.33) is written as Eq. (3.34). ∑ ˆ κuij,vmn ±κvij,umn 1 + δij,mn Kr i < j, p Suv ±Svu = m < n, q ij ≤ mn 0≤s<S ¯−1 µ̄ pqs (ςijp,mnq,s ±ςmnp,ijq,s )ℓs (ξ) det Js dKr (3.34) ¯−1 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 Kr i<j m<n ij ≤ mn p≤q 0≤s<S ¯−1 ¯−1 (µ̄ pqs ± µ̄qps )(ςijp,mnq,s ±ςmnp,ijq,s ) det Js ℓ(ξ)dKr 1 + δpq (3.35) In Eq. (3.35), the bound ij < mn reﬂects 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 coeﬃcients C± ijmns as formulated in Eq. (3.37) followed by the multiply-add operations encoded in Eq. (3.38). ˆ ◦± κuij,vmn ± κvij,umn Suij,vmn,s , ℓs (ξ)dKr 2 + 2δij,mn Kr −1 −1 ∑ (µ̄ ¯ ¯ ± µ̄ )(ς ± ς ) ijp,mnq mnp,ijq pq qp ± Cij,mn,s, det J 1 + δpq p≤q 1 (Suv ±Svu )= 2 ∑ ◦± C± ij,mn,s Suij,vmn,s (3.36) (3.37) x=xs (3.38) i<j m<n ij ≤ mn 0≤s<S To evaluate the actual number of multiply-add operations required for evaluation of the stiﬀness matrix, one needs to observe that C− ij,mn,s = 0 for cases where ij = mn. 26 ¯−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. Notation Deﬁnition ΓD part of the boundary with essential BCs Ωi ith sub domain of Ω ¯ri µ̄ relative magnetic permeability tensor in Ωi Ei electrical ﬁeld in sub domain i ei tangent component of Ei , i.e. πτ (Ei ) n̂i 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 2 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 3.4.5 Curvature The UA approach presented here is intrinsically capable of handling elements with curved geometry. Also, our experiments show that the method can be successfully 27 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. 28 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 transﬁnite interpolation [66] of the faces as depicted in Fig. (3.2). The resulting face expressions are then used to build a tetrahedral transﬁnite interpolation of the physical volume Kp . Hence, the ﬁnal transﬁnite mapping conforms exactly to the surface of the sphere. As depicted in Fig. (3.3), the resulting coordinate mapping is then ﬁtted 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 . 29 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. 3.5 CIMPT Approach for the PML Accurate realization of conformal PML layer conditions leads to artiﬁcial MPTs that are continuous functions of global coordinates and often involve curved element geometries. As a direct consequence, the CIMPT approach will be very proﬁtable 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 ﬁeld 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 brieﬂy 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 Λ̄ µ̄ (3.39) ¯ ε̄¯ = ε0 Λ̄ (3.40) 30 ¯ can be written as Eq. (3.41) in In a 3D spherical coordinates system (r, θ, ϕ), Λ̄ which sr (r) is deﬁned as in Eq. (3.42). The detailed derivation is conducted in [71]. ( )2 r̃ 1 r̂ + ϕ̂sr (r)ϕ̂ + θ̂sr (r)θ̂ r sr (r) σ(r) sr (r) = α(r) + ȷ ω ¯ Λ̄ r,θ,ϕ = r̂ r = Rin + ζ3 (3.41) (3.42) (3.43) In Eq. (3.41), r̃ is a variable obtained by means of complex coordinate stretching as r̃ = Rin + ζ̃3 where ζ̃3 is deﬁned 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 diﬀerence 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 = ζ3 ( ˆ s(κ)dκ = 0 0 σ(κ) α(κ) + ȷ ω ) dκ (3.44) ¯ 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 eﬀectively handled using the CIMPT approach. Some numerical results concerning the proposed CIMPT PML are presented in section 3.6. 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 31 R 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 TM R Xeon⃝ processors. The Krylov subspace solver’s iterations were truncated at a relative tolerance of 10−4 . 3.6.1 Waveguide 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 deﬁnition 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. Method CIMPT 3-Layer DoF 4370 2570 Elements 864Tet 504Tet CPU Time 3sec 1sec Field Basis 2 2 Mat. Basis 2 2 5-layer 3290 648Tet 2sec 2 2 7-layer 5434 1067Tet 3sec 2 2 Ilic et al. 205 3Hex N.A. (4, 2, 4), (4, 2, 7) 1 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 ﬁeld 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 32 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 ﬁeld solution will be more complex compared to the piece-wise constant material case. Thus, ﬁeld 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 ﬁeld 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 diﬀerent 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 ﬁeld basis along (x̂, ŷ, ẑ) directions. 33 Figure 3.6: (a)Electrical ﬁeld intensity. (b) Conformal-DDM’s Partition of the Mesh. 3.6.2 Luneburg Lens The Luneburg lens problem discussed here is a dielectric sphere with a relative permittivity tensor of ε̄¯ = (2 − r2 ¯ )Ī 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/reﬂector devices operated in microwave and millimeter wave frequencies [72] [73]. Hence, the Luneburg lens is a good example for the purpose of verifying the ﬁdelity 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 eﬀects 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). 34 Figure 3.7: Schematic diagram showing the conﬁguration 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 ﬁeld intensity. The plane-wave focusing eﬀect 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 ﬁnal solution was obtained in 88 iterations (13 seconds). 35 The Lens Antenna Our next experiment involves exciting lens antennas of various radius and calculating their far ﬁeld 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 z 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 ﬁeld solution for the R = 6λ0 case is plotted in Fig. (3.9). One can observe the ﬂattening of the phase front as the outgoing waves propagate away from the aperture of the waveguide. Figure 3.8: Schematic diagram showing the conﬁguration of the R = 6λ0 Luneburg lens antenna problem. Curved elements were applied to the surface of the lens only. 36 Figure 3.9: A three-slice ﬁeld plot of the E · ŷ ﬁeld 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 37 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 ﬁnish 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 }. 3.6.3 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 ﬁeld 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 38 Figure 3.11: Schematic diagram showing the conﬁguration of the PML-lens antenna problem. Figure 3.12: (a) Snapshots of the magnitude of the electric ﬁeld distribution on yzplane. (b) Snapshots of the magnitude of the electric ﬁeld distribution on zx-plane. strong attenuation of the radiated ﬁeld which validates its proper implementation. Fig. (3.13) plots the resulting radiation pattern in the ϕ=0 plane. 3.7 Conclusion An eﬃcient UA approach capable of handling curved geometries and continuously varying MPTs is introduced in this chapter. Compared to the conventional approach 39 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 ﬂexible 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 veriﬁcation 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 ﬂuctuations; b) in doped semiconductors 40 and plasma antennas/reﬂectors, dielectric properties may follow a continuous proﬁle following various physical phenomena. Hence, the introduction of the UA method will play an inﬂuential role in future developments associated to the CIMPT approach. Finally, it worth mentioning that complete utilization of the beneﬁts of the proposed CIMPT and UA methods can be achieved by means of implementations with higher orders of ﬁeld basis functions. 41 Chapter 4: Ferromagnetic Nano-wires: Homogenization and Applications 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 eﬀective 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 eﬀects of FMNWs employed in RF devices. 4.1 Introduction Metamaterials can be broadly deﬁned as artiﬁcial eﬀective 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 artiﬁcially structural materials is at the stage where interdisciplinary interest is taking place, such as multi-scale 42 metamaterials, reconﬁgurable metamaterials and the combination with nanotechnology. Compared with the earlier generation of metamaterials, the new generation of metamaterials incorporate ﬁner structural scales, for instance, the nanometer or even atomic scale, which suﬀer less from scattering loss as well as fabrication diﬃculties. While the potentials of metamaterials and their possible applications are appealing, it is often diﬃcult to model them accurately due to inherently the very ﬁne 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 reﬂection and transmission coeﬃcients on ﬁnite lengths of the electromagnetic metamaterials; in [47], the authors have determined the eﬀective material constants based on averaging the local ﬁelds; and in [78], a homogenized description is derived of periodic metamaterials, which are made of magnetodielectric inclusions, and subsequently, closed-form expressions for the eﬀective 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 ﬁelds are deﬁned. 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), 43 anisotropy control and self-bias. Compared against traditionally bulky ferromagnets, the utilization of FMNW metamaterials can signiﬁcantly 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 eﬀective 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 eﬀective permeability tensor, by ﬁrst computing the static magnetic ﬁeld 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 deﬁnition 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. 4.2 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 eﬀective dynamic permeability of the array in a complex tensorial form. 44 Figure 4.1: Sketch of FMNW substrate structure. The proposed homogenization procedure can be divided into the following major steps. 4.2.1 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 ﬁeld 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 ﬁeld Hext dc . Fig. (4.2) shows a schematic hysteresis curve. 4.2.2 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 (4.1) ρs = Mdc · n̂, (4.2) and 45 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: 1 Φ(⃗r) = 4π ˆ ρv 1 ′ ′ dv + r − ⃗r | 4π V |⃗ ˆ ρs ′ ′ ds . r − ⃗r | S |⃗ (4.3) Since in our case Mdc is constant, it results in ρv = 0; thus Eq. (4.3) can be simpliﬁed as, 1 Φ(⃗r) = 4π ˆ ρs ′ ′ ds . r − ⃗r | S |⃗ (4.4) 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]. 46 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). (4.5) 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 |], 2 (4.6) 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, respectively. Φ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 4π |⃗r − ⃗r′ | ′ (4.7) ⃗ 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 ′ 1 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 signiﬁcantly using 2D-FFT for the summation 47 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 ﬁeld Hm is then readily available for every ferromagnetic nanowire from the magnetic potential: Hm (⃗r) = −∇Φ(⃗r). 4.2.3 (4.8) Computation of the small signal solution of the LLG equation Assuming that the corresponding ac magnetization ﬁeld h is much smaller than the external DC bias, a decomposition of the total eﬀective magnetization ﬁeld Hef f into the DC and the small ac signal hence proves to be beneﬁcial. Namely, Hef f = Hdc + h, Hdc = Hext dc + Hm . (4.9) Note that other contributions to the eﬀective magnetic ﬁeld Hef f , such as the exchange interaction, anisotropic ﬁeld and the thermal ﬁeld, are neglected since they are much less signiﬁcant compared against the demagnetization ﬁeld. Similarly, the total magnetization moment M is separated into DC and ac components as well: M = Mdc + m. 48 (4.10) 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: αG ∂M ∂M = −µ0 |γg |M × (Hef f + ), ∂t µ0 |γg |Ms ∂t (4.11) 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 ẑ × m, Ms (4.12) with M being the magnitude of M. Eq. (4.12) can be written into: M my Ms M mx = −µ0 |γg | (M hx + mz Hx − mx Hz ) − jωαG Ms jωmx = −µ0 |γg | (−M hy + my Hz − mz Hy ) + jωαG jωmy jωmz = −µ0 |γg | (mx Hy − my Hx ) The 2 × 2 matrix form of x and y components is: 49 (4.13) [ Axx = Axx Axy Ayx Ayy ][ mx my [ ] hx = µ0 |γg | M hy ] −jHx Hy Ms |γg |2 µ0 2 +Hz Ms |γg |µ0 ω−jM αG ω 2 ωMs 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 ωMs (4.14) 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 deﬁned similarly. 4.2.4 Extraction of the equivalent/effective permeability tensor With the small signal magnetization moment, m, computed in response to given external applied ac magnetic ﬁeld, h (through the excitation of three diﬀerent polarizations), it follows that the 3×3 magnetic susceptibility tensor χω can be determined. Consequently, the relative equivalent/eﬀective homogenized permeability tensor µω of FMNWs is also available, and we note that the extracted permeability tensor is nonreciprocal and dispersive. m(ω) = χω h(ω) 50 χω = µ0 |γg |M Axx Axy Ayx Ayy 0 0 0 (4.15) and, µxx µxy 0 µω = χω + I = µyx µyy 0 . 0 0 1 (4.16) 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 bulk. 4.2.5 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. 51 Figure 4.3: Flow-chart of FMNW based metamaterial homogenization process and application enumeration. 4.3 4.3.1 Numerical Results Validation of extracted permeability tensor We ﬁrst validate the proposed numerical homogenization procedure by comparing the eﬀective material property computed against an analytical model proposed in [82] for a nanowired substrate of diﬀerent magnetic alloys. Speciﬁcally, a nonreciprocal microstrip phase shifter on a substrate with two diﬀerent 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 52 diﬀerent heights of the membrane thickness, with hCoF e = 0.68 and hN iF e = 0.84, respectively. The analytical homogenized model proposed in [82] for the nanowires has the form: µ −jκ 0 µ = jκ µ 0 , 0 0 1 (4.17) with 2πMs γP hfF M R (m2 + 1) , fF2 M R − f 2 4πMs γP hf m κ = , fF2 M R − f 2 µ = 1+ (4.18) (4.19) 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 ﬁeld 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 ﬁnally P is the porosity. Figs. (4.4) and (4.5) below show the comparisons of the eﬀective 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 ﬁeld, after saturation), so it is deﬁnitely still valid for applying a small signal model analysis. From these ﬁgures, 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. 53 Figure 4.4: Comparison of eﬀective relative permeability tensor of CoFe nanowire membrane obtained by diﬀerent methods. Figure 4.5: Comparison of eﬀective relative permeability tensor of NiFe nanowire membrane obtained by diﬀerent methods. 54 4.3.2 Double ferromagnetic resonance phenomenon FMR has been studied extensively due to their applications into microwave devices. Speciﬁcally, 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 ﬁeld 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 ﬁelds 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 ﬁgure, 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 ﬁelds. Figure 4.6: Microstrip line sketch. 55 We took the hysteresis loop ﬁgure 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 ﬁgure, the pairs of the applied external DC magnetic ext ﬁeld 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 ﬁeld from −5.0KOe to 5.0KOe with a step of 500Oe, following the speciﬁed branch indicated in Fig. (4.7). Table ext (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 red). Table ext Hdc (KOe) Mn ext Hdc (KOe) Mn 4.1: External DC Magnetic ﬁeld vs. -2.5 -2.0 -1.5 -1.0 -0.5 -1.0 -0.96 -0.9 -0.82 -0.68 1.0 1.5 2.0 2.5 3.0 -0.002 0.24 0.48 0.69 0.85 56 Normalized Magnetization 0.0 0.5 -0.48 -0.25 3.5 4.0 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 eﬀective permeability tensors of the FMNW from 20GHz to 40GHz. Subsequently, the extracted permeability tensors are employed to compute the S21 coeﬃcients through full-wave electromagnetic computations. However, it is worth mentioning that in order to account for the fabrication imperfections, fringeﬁeld eﬀects, and the measurement uncertainties, we adopted a similar adjustment as suggested in [29] [83]. Namely, µc = I + q(µh − I), (4.20) where µc is the eﬀective 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 ﬁlling 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 eﬀective permeability given by: µc = 1 + q(µh − 1), (4.21) 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 eﬀects of the FMNWs. Although, with the adjustment introduced, the non-reciprocal eﬀects induced will not be signiﬁcant. For example, we have the eﬀective permeability, µc ext at 33GHz and the external applied magnetic ﬁeld Hdc = −5KOe, as: 1.068 − j0.2 −0.2 − j0.055 0. µc = 0.2 + j0.055 1.068 − j0.2 0. . 0. 0. 1. 57 (4.22) 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 ﬁgure, 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 ﬁgure, 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 ﬁgures in Fig. (4.8), that with the external applied ext DC magnetic ﬁeld 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 ﬁeld. However, it is in the region where the applied ﬁeld is below saturation, 58 ext 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 ﬁeld distributions along the microstrip line, with a single external bias ﬁeld set to be 0.5KOe. As can be seen from the ﬁgure, there are two diﬀerent frequencies, 22GHz and 30GHz, at which the electromagnetic energy gets absorbed almost completely by the FMNW substrate. Figure 4.9: Electric ﬁeld distributions along the microstrip line with FMNW subext strate and biased with the same external bias ﬁeld, Hdc = 0.5KOe, at two diﬀerent resonance frequencies. 4.3.3 Dual-band integrated edge-mode isolator In this section we ﬁrst 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 59 the substrate as well as the geometrical description of the isolator are included in Fig. (4.10). 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 ﬁeld 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. 60 (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 61 the isolator obtained by diﬀerent 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 process. 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 ﬁeld distributions along the microstrip line of the isolator is plotted in Fig. (4.14) for these two center frequencies. 4.3.4 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 62 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 eﬀective permeability tensor of FMNW substrate obtained through the proposed homogenization procedure, and the relative permittivity of porous alumina membrane with and without magnetic nanowire ﬁllings 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 coeﬃcients, with Port 1 being excited, over the frequency range, from 20GHz to 40GHz. As evidenced from the ﬁgure, the circulator has two working frequency bands: 25GHz to 29GHz and 31.5GHz to 33GHz. In the 63 Figure 4.14: Electric ﬁeld distributions along the microstrip line of the isolator at two diﬀerent working frequencies. ﬁrst 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 ﬁeld 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. 4.4 Conclusion A numerical homogenization procedure is proposed in this chapter to study the unique properties of FMNW metamaterials. Numerical results validate the proposed 64 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 ﬁeld computations. 65 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 ﬁeld plotting. 66 Chapter 5: Self-Consistent Integration of Discontinuous Galerkin Time Domain Method and Circuit Simulator for Microwave Components and Devices with Circuit Network Analysis 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 eﬃcient 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 ﬁeld and circuit device simulations. This demands sophisticated and eﬃcient 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. 67 Furthermore, allowing for discontinuity of the discrete trial and test spaces, DGTD can provide a substantial amount of ﬂexibility 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 eﬃciently. 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 buﬀer information speciﬁcation (IBIS) standard. IBIS is a speciﬁcation that deﬁnes the format on the information about the input/output buﬀers 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 eﬃcient solution for wide-ranged applications. Be it SPICE or IBIS, the interface between the EM and circuit parts is similar. Generally speaking, the interface of the two parts is deﬁned as a planar surface port, on which the voltage across is calculated in EM simulation at each time step, and 68 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 co-simulation. 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 ﬁrst. It is then followed by the illustration of its integration to DGTD solver. Again some interesting numerical examples will be presented. 5.1 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 diﬀerent port establishing strategies are investigated for the EM-circuit couplings. Several numerical examples are shown to validate the accuracy and application of the proposed method. 5.1.1 Introduction 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 ﬁeld and circuit simulations [85] [86]. On one hand, EM ﬁelds at high frequencies would cause wave phenomena such as EM radiation, signal 69 delay and distortion, skin and proximity eﬀects 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 eﬀorts have been made to incorporate EM eﬀects into IC simulation. These works can be roughly classiﬁed into two categories. One is a hybrid method of bringing the circuit eﬀects 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 eﬀects 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 eﬀects as distributed circuit sources, and simulate them together with the devices within only circuit simulator. While worthwhile, the latter method suﬀers from accuracy problems, especially when at high frequencies. In contrast, the ﬁrst 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] oﬀers 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 70 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 eﬃcient 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 modiﬁed 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 ﬁeld-eﬀect 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 ﬁeld-eﬀect transistors (MOSFETs), has to be redeﬁned 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 diﬀerent characterization models and parameters (such as the gate oxide thickness, temperature dependence coeﬃcient, etc.), and the number keeps growing. These 71 model levels should be indispensable to circuit analysis of real-life devices, but reimplementing these parameters and deﬁnitions 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 deﬁnitions 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 diﬀerence (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 72 DGTD and the circuit part is simulated in SPICE simultaneously with respect to self-consistency in an iterative manner. The developed method oﬀers transient and wide frequency-band computations and the inclusion of arbitrary linear/nonlinear passive/active devices easily, for complex conﬁgured integrated circuit simulation and analysis. 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 diﬀerent viewpoints. Finally numerical results to validate the developed method are presented. 5.1.2 Formulation 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 ∂t ∂E ∇ × H = ϵ̄¯ ∂t ¯ ∇ × E = −µ̄ (5.1) ¯(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 73 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 deﬁne 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 (5.2) πτ (ui ) = n̂i × ui × n̂i |∂Ki (5.3) 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 ﬁrst order Maxwell’s equations in time domain. By deﬁning the ﬁnite-dimensional discrete space as Vhp = {v ∈ [L2 (Ω)]3 : v|K ∈ [Pp (K)]3 , ∀K ∈ Th }, with L2 (Ω) the linear space of square-integrable vector ﬁelds 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 µ̄ dΩ w· ∇×E+ ∂τ Ω ˆ 74 ( ) ϵ̄¯ · ∂E − v· ∇×H+ dΩ ∂τ Ω ˆ ˆ + {{v}} · [[H]]γ dS − {{w}} · [[E]]γ dS Fh Fh ˆ ˆ [[v]]π · [[E]]π dS + f [[w]]π · [[H]]π dS = 0, +e ˆ Fh Fh ∀ (w, v) ∈ Vhp × Vhp . (5.4) where the notations are deﬁned as follows, ( ) {{u}} = πτ (ui ) + πτ (uj ) /2 [[u]]γ = γτ (ui ) + γτ (uj ) on FhI [[u]]π = πτ (ui ) − πτ (uj ) {{u}} = πτ (u) B [[u]]γ = γτ (u) on Fh [[u]]π = πτ (ui ) (5.5) (5.6) √ 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 ﬂux. 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 ∑ (i) (i) the local electric and magnetic ﬁelds 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 coeﬃcients for the basis functions. By separating the w, v testing in Eq. 75 (5.4), we are able to build a local semi-discrete system within element Ki as, Mϵ ∂ei = Se hi − Fiie hi + ePiie ei ∂t ij −Fij e hj + ePe ej Mµ (5.7a) ∂hi = −Sh ei + Fiih ei + f Piih hi ∂t ij +Fij h ej + f Ph hj (5.7b) where ei (t) and hi (t) are the time dependent coeﬃcient vectors, respectively. j is the indices of the neighboring elements of element i. The local matrices are expressed as ˆ i (Mϵ )nm = vin · ϵ̄¯i vim dV Ki ˆ i (Se )nm = vin · ∇ × wim dV Ki ˆ 1 ii (Fe )nm = πτ (vin ) · γτ (wim )dS 2 FhI ˆ 1 ij (Fe )nm = πτ (vin ) · γτ (wjm )dS 2 FhI ˆ ii (Pe )nm = − πτ (vin ) · πτ (vim )dS FhI ˆ ij (Fe )nm = πτ (vin ) · πτ (vjm )dS FhI and ˆ (Miµ )nm ¯i wim dV win · µ̄ = ˆ Ki (Sih )nm win · ∇ × vim dV = ˆ 1 πτ (win ) · γτ (vim )dS 2 FhI ˆ 1 πτ (win ) · γτ (vjm )dS 2 FhI ˆ − πτ (win ) · πτ (wim )dS FhI ˆ πτ (win ) · πτ (wjm )dS Ki (Fiih )nm = (Fij h )nm = (Piih )nm = (Fij h )nm = FhI 76 The ﬁrst-order time derivatives in the resulting diﬀerential Eqs. (5.7) can be discretized using central diﬀerence approximation with second order accuracy. The electric ﬁeld unknowns are evaluated at tn = nδt and the magnetic ﬁeld unknowns are evaluated at tn+1/2 = (n + 1/2)δt. For the two penalty terms arising from upwind ﬂux n+ 1 n+ 1 2 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 i n+ 21 −δtFij e hj n+ 32 Miµ hi n + eδtPij e ej n+ 12 = (Miµ + f δtPiih )hi (5.8a) + δt(−Sih + Fiih )en+1 i n+ 21 n+1 +δtFij + f δtPij h ej h hj (5.8b) 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]. 5.1.3 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 ﬁrst decompose the whole system into two subsystems, one including the 77 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 ′ ′ AB AB VEM = VCKT ′ ′ AB AB IEM = ICKT (5.9a) (5.9b) AB through a planar surface port as shown in Fig. (5.1). In Eq. (5.9), VEM represents AB the voltage across the surface port, and IEM represents the current ﬂowing on the ′ ′ AB planar surface in the EM structure. VCKT denotes the voltage across the circuits, ′ ′ AB 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 78 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 ) (5.10a) VEM = g(IEM ) (5.10b) By applying one of the self-consistent conditions Eq. (5.9b), we can obtain VEM = g(IEM ) = g(ICKT ) = g(f (VCKT )) (5.11) 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. f g Information exchange between DGTD and SPICE goes in a VCKT − →I− → VEM fashion. Note that Eq. (5.9a) requires VCKT and VEM be the same. If only Eq. (5.9b) 79 is ensured, the diﬀerence of VCKT and VEM could be noticeably large after such a information ﬂowing 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 diﬀerence of VCKT and VEM suﬃciently 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 ﬁnd 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 ϵ. f g k k k k • Start from VCKT − → IEM (ICKT )− → VEM . k k • Compare VCKT with VEM ; If their diﬀerence is bigger than ϵ, denote k+1 k VCKT = VCKT + ∆V (5.12) Otherwise, stop looping. k+1 k+1 • To ﬁnd VCKT = g(f (VCKT )), apply the Taylor’s expansion of the function k VEM = g(f (VCKT )) at VCKT , ignore the second and higher order terms, and we get k+1 k+1 k )) = g(f (VCKT )) + = g(f (VCKT VCKT ∂gf | k ∆V ∂V VCKT (5.13) k k • Since VEM = g(f (VCKT )), associating with Eq. (5.12), we have ∆V = k k VCKT − VEM ∂gf | k −1 ∂V VCKT 80 (5.14) 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 modiﬁcation 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. 81 Starting from the boundary conditions of the electric and magnetic ﬁelds on the port, we can obtain the relationships on surface ΓCKT as follows, n̂i × Hi + n̂j × Hj = JCKT = ICKT w n̂i × Ei + n̂j × Ej = 0 (5.15) (5.16) The magnetic ﬁeld relationship Eq. (5.15) is used to represent the currents from the circuits. Meanwhile, the electric ﬁelds 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, Residuals Function Space (1) RΓCKT = [[H]]γ − JCKT (2) RΓCKT = [[E]]γ (1) ∈ H(divτ , ΓCKT ) (5.17) ∈ H(divτ , ΓCKT ) (5.18) (2) 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 (2) 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 + dΩ ∂τ Ω ( ) ˆ ϵ̄¯ · ∂E − v· ∇×H+ dΩ ∂τ Ω 82 ˆ + Fh {{v}} · ([[H]]γ − JCKT ) dS ˆ − Fh {{w}} · [[E]]γ dS = 0 (5.19) Separating the w and v testing, we can obtain local the semi-discrete equations for element Ki as ∂ei = Se hi − Fiie hi + BJe JCKT − Fij e hj ∂t ∂hi Mµ = −Sh ei + Fiih ei + Fij h ej ∂t Mϵ where (BJe )n = 1 2 ´ FhI (5.20a) (5.20b) πτ (vin )· l̂ and l̂ denotes the unit vector along the length direction of the port. Likewise as Eq. (5.8), we apply the central diﬀerence approximation to 1 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 i n+ 1 n+ 12 +BJe JCKT2 − δtFij e hj n+ 32 Miµ hi n+ 12 = Miµ hi (5.21a) + δt(−Sih + Fiih )en+1 i n+1 +δtFij h ej (5.21b) 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. 83 To clarify the implementation process, we assume the following responses of EM n n , vcap , and circuit systems are known at time tn = nδt. These responses include VCKT n+1 n 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 speciﬁcally with all the input/output variables as, n+ 1 n+1 n+1 n+1 n n , iind ) = fSP ICE (VCKT , VCKT , vcap , inind ) (ICKT2 , vcap (5.22) The relationship in Eq. (5.22) indicates the SPICE simulation from tn to tn+1 . n+1 n The voltage pair VCKT , VCKT are used to construct a piecewise linear voltage pulse n+1 n in SPICE with VCKT at tn and VCKT at tn+1 , in order to drive the whole circuit. n vcap 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 ﬁle is shown in Fig. (5.4), n+ 1 As for the output shown in Eq. (5.22), the computed port current ICKT2 is ﬂowing n+1 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 steps: 84 Figure 5.4: SPICE netlist example. n n • Start from time tn , and assume VEM = VCKT . • Set the iteration number k = 0, and make a guess of the initial voltage value n+1,k n−1 n n across the circuit as VCKT = VCKT + (VCKT − VCKT ) • Reset the SPICE conﬁguration commands by updating the voltage source and the initial condition of the circuit, including all historical circuit status parameters such as node voltages, current ﬂows through all inductors and all DC values of the sources at tn . • SPICE takes voltage as a piece-wise linear voltage source deﬁned from tn to n+1,k n tn + δt, with the range from VCKT to VCKT . 85 • 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 n+1,k Eq. (5.8) and Eq. (5.21). Then compute VEM from the line integral of the electric ﬁeld through the circuit port. n+1,k+1 • Calculate VCKT by using the relations of Eq. (5.14) and Eq. (5.12). • If n+1,k+1 n+1,k |VCKT − VCKT | n+1,k |VCKT | <ϵ (5.23) n+1 n+1 the nth to (n + 1)th time step marching is completed. We save VCKT (VEM ) and the states of the circuit at tn+1 , and begin the co-simulation at next time step. If Eq. (5.23) is not satisﬁed, increase k and go back to continue looping until convergence. Fig. (5.5) shows the ﬂowchart 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 ﬁrst, like 86 Figure 5.5: Flowchart of DGTD integrated with SPICE on vector ﬁtting basis [109,110], and then write the netlists for the corresponding S/Yparameter or admittance matrices. More speciﬁcally, 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]. 5.1.4 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 ﬁnite 87 width to a circuit of zero width, and to approximate the problem of reality in a most physical way. Based on diﬀerent 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 ﬁrst 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 deﬁne the voltage across the surface port in an average sense, ´W ´L ´W E(x, y) · l̂dydx V (x)dx 0 0 VEM = = 0 , ´W W dx 0 (5.24) 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), VEM n 1 ∑ = w(xi )V (xi ) W i=1 (5.25) 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)) 88 (5.26) 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: ˆ ˆ W I= J(x)dx = 0 W F (VEM (x)) dx (5.27) 0 Similarly, we express the integration into a weighted summation by Gaussian integration rule, and obtain: ˆ W I = 0 = W 2 = W 2 F (VEM (x)) dx ( )) ˆ 1 ( W W F VEM x+ dx 2 2 −1 ( ( )) n ∑ W W xi + wi F VEM 2 2 i=1 89 = n ∑ w(xi )I(xi ) (5.28) i=1 As a result, the total current on the port can be seen as a weighted summation of the currents at diﬀerent locations along the width direction of the port. To implement this idea, the voltages are calculated at diﬀerent 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 diﬀerent locations. Therefore, the second way takes longer time compared to the average voltage method. Comparing these two schemes, the ﬁrst method evaluates the voltage as a weighted summation of the voltages at diﬀerent locations of the port, while the latter method expresses the current as a weighted summation of the currents at diﬀerent 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. 90 5.1.5 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. Microstrip We ﬁrst 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. 91 As for the reference, the discrete port deﬁnition 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 deﬁnitions in tetrahedron meshes. However, we can modify the line port of CST a bit to let it ﬁt in a surface port deﬁnition. 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 deﬁned 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 92 proposed DGTD-SPICE co-simulation method agree very well with the results from CST. Figure 5.9: Result comparison at near end of microstrip. To investigate the diﬀerence 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 93 Figure 5.10: Result comparison at far end of microstrip. Figure 5.11: Microstrip line with nonlinear circuit network. 94 Figure 5.12: Comparisons of numerical results with diﬀerent 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 oﬀers 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). 95 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 reﬂect 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 eﬀects. 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 diﬀerent 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 diﬀerence actually should be determined by the width of the surface port and the voltage diﬀerence from one side to the center. In real-life PCB, the width of the circuit elements would be smaller and the voltage diﬀerence 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, 96 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 conﬁguration. 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 s 0.50 mm h 0.40 mm d 0.20 mm 97 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 ﬁrst 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 ﬁlter and reaches at the SRD package at Port 5. The SRD is used as a charge controlled switch for waveform sharpening. In the ﬁrst 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, ﬁltered by the DC blocking 1.5-pF capacitor, and produces a sub-nanosecond pulse 98 Table 5.2: Step-recovery diode description Name Description Parameter IS Saturation current per unit 0.5 pA area RS Ohmic series resistance 0.13 Ω N Emission coeﬃcient 1.3 VJ Contact potential at area 0.5 V junction M Grading coeﬃcient at area 0.235 junction BV Reverse breakdown voltage 60 V IBV Current at breakdown volt- 10 µA age TT Transit time 30 ns FC Coeﬃcient 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) conﬁrm 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. 5.1.6 Conclusion 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 99 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 generator. 100 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 diﬀerent 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. 5.2 5.2.1 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 buﬀers 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]). 101 Figure 5.17: Input buﬀer. In general, the data in an IBIS ﬁle is used to construct a buﬀer 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 buﬀer I-V (current versus voltage) and switching (output voltage versus time) characteristics. 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 buﬀer’s capacitance) and the package parasitic parameters of lumped RLC combinations. For the input buﬀer, the required information is: 1. The buﬀer’s I-V characteristics 2. The buﬀer’s capacitance 3. Package parameters 102 Figure 5.18: Output buﬀer. Fig. (5.18) shows the structure of output buﬀer. 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 buﬀer’s capacitance) and the package parasitics. An output or I/O buﬀer is characterized using the following information: 1. The buﬀer’s output I-V characteristics when the output is in the logic low state 2. The buﬀer’s output I-V characteristics when the output is in the logic high state 3. The buﬀer’s output I-V characteristics when the output is forced below ground and above the power supply rail 4. The time it takes a buﬀer’s output to switch logic states from low to high and high to low 5. The buﬀer’s capacitance 6. Package parameters As far as the IBIS ﬁle format is concerned, it is a readable ACSII ﬁle that contains the ﬁrst part of the header general information about the ﬁle, the device and the 103 company, then the second part of the device name, pinout and pin-to-buﬀer 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 included. 5.2.2 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 buﬀers are detailed below in 5.2.3 and 5.2.4, respectively. 5.2.3 Transient representation of IBIS input buffer model Compared to output buﬀer model, the input buﬀer 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 buﬀer when the GND clamp and the power clamp diodes are turned on, respectively. The GND clamp is active when the buﬀer 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. 104 5.2.4 Transient representation of IBIS input buffer model As previously pointed out in Fig. (5.18), the output buﬀer 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 ﬁxed resistance direc whose value is speﬁciﬁed in the ﬁle. 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 brieﬂy given, the details of this method is referred in [115]. The current ﬂowing through pullup transistor is determined by the voltage across it referenced to Vcc , noted as Ipu (V ), while the current ﬂowing 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-oﬀ state and 1 turn-on state. The change of multiplier from 0 to 1 represents the transistor status change from turn-oﬀ 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). 105 Figure 5.19: Illustration of the introduced multiplier switching behavors. Figure 5.20: Equivalent analog behavior model. With the introduced multiplier coeﬃcients, 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 buﬀer output signal. The current source Ipc (V ) represents the Vcc -clamping eﬀect and the current source Igc (V ) represents the Gnd-clamping eﬀect. 106 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 ﬁxed 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 (5.30) 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 coeﬃcients, 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 107 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 ﬁle are supplied over the range of voltages the output could possibly generate or experience. However, if a buﬀer 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]. 5.2.5 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 buﬀer 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 beﬀer gate. 108 Figure 5.21: Microstrip line with IBIS models. Figure 5.22: Package top view. 109 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 simpliﬁed 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 ﬁrst 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 ﬁles mpc52001.ibs and t37zi t.ibs can be downloaded from the homepages of Freescale and 110 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. 111 Figure 5.26: Result comparison of voltages at the microcontroller (port1). Micron, respectively. Here in this example the ouput buﬀer of the microcontroller and the input buﬀer 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 112 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. 5.2.6 Conclusion 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 modiﬁcation of DGTD. 113 Chapter 6: Conclusion The complexity of modern electromagnetic applications demand eﬀective and efﬁcient numerical methods. Various numerical methods are out there with strength emphasizing diﬀerent aspects. Therefore, to apply the appropriate solver for a speciﬁc problem is of vital importance. Usually diﬀerent solvers can be applied to solve the same problem, but the convenience for modeling, the eﬃciency 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 eﬃcient 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 ﬁeld 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 114 Chapter 5, full-wave transient EM circuit phenomena with linear/non-linear eﬀects are investigated through DGTD method coupled with SPICE and IBIS circuit models in a self-consistent integration scheme. Numerical examples all validate the proposed methods. 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 suﬃcient anymore. It is important to adaptively choose an appropriate step size at each time stepping of the process, so that the non-linear eﬀect 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 diﬀerent 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 eﬀects from the minor factors that are responsible for the whole system is therefore important. It can greatly enhance the eﬃciency 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 115 can somewhat automatically use the optimal solver for the speciﬁc type of problme is some path we endeavor for. 116 Bibliography [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, 2002. [3] J.-f. Lee, D.-K. Sun, and Z. Cendes, “Full-wave analysis of dielectric waveguides using tangential vector ﬁnite 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 ﬁnite 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 ﬁnite 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 ﬁnite 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 ﬁnite elements for electromagnetic ﬁeld 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 ﬁnite 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. 117 [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 ﬁnite 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 ﬁnite 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 ﬁnite-diﬀerence approximation of the time-domain electromagnetic-ﬁeld equations,” Electromagnetic Compatibility, IEEE Transactions on, no. 4, pp. 377–382, 1981. [17] S. Dey and R. Mittra, “A locally conformal ﬁnite-diﬀerence time-domain (FDTD) algorithm for modeling three-dimensional perfectly conducting objects,” Microwave and Guided Wave Letters, IEEE, vol. 7, no. 9, pp. 273–275, 1997. [18] T. G. Jurgens and A. Taﬂove, “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 ﬁnite-element methods,” IEEE transactions on antennas and propagation, vol. 45, no. 3, pp. 430–442, 1997. [20] S. D. Gedney and U. Navsariwala, “An unconditionally stable ﬁnite 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 ﬁnite-element method for electromagnetic simulations,” Antennas and Propagation, IEEE Transactions on, vol. 50, no. 11, pp. 1624–1632, 2002. 118 [22] D. A. White, “Orthogonal vector basis functions for time domain ﬁnite 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 ﬁnite-diﬀerence 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-diﬀusion 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 ﬁnite 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, “Eﬀective 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. 119 [33] S. S. Zivanovic, K. S. Yee, and K. K. Mei, “A subgridding method for the time-domain ﬁnite-diﬀerence method to solve Maxwell’s equations,” Microwave Theory and Techniques, IEEE Transactions on, vol. 39, no. 3, pp. 471–479, 1991. [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 ﬁnite-diﬀerence and ﬁnite-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) ﬁnite-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 ﬁnite 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. 120 [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 ﬁnite 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 eﬀective permittivity and permeability of metamaterials from reﬂection and transmission coeﬃcients,” Phys. Rev. B, vol. 65, no. 19, 2002. [47] D. R. Smith and J. P. Pendry, “Homogenization of metamaterials by ﬁeld 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 reconﬁgurable 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 ﬁnite-element time-domain method,” Microwave Theory and Techniques, IEEE Transactions on, vol. 60, pp. 2684–2692, Sept. 2012. [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. 121 [54] J. S. Wang and N. Ida, “Curvilinear and higher order edge ﬁnite elements in electromagnetic ﬁeld 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 ﬁnite 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 ﬁnite elements in nonlinear magnetic ﬁeld problems,” IEEE Trans. Magn., vol. 33, Sept. 1997. [57] M. Dufresne and P. P. Silvester, “Universal matrices for the N-dimensional ﬁnite 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 ﬁnite 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 ﬁnite element method (NEFEM),” International Journal for Numerical Methods in Engineering, vol. 76, pp. 56–83, Feb. 2008. 122 [66] A. Perronnet, “Interpolation transﬁnie 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, “Diﬀerential forms, metrics, and the reﬂectionless 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, 2000. [69] F. Teixeira, K.-P. Hwang, W. Chew, and J. Jin, “Conformal PML-FDTD schemes for electromagnetic ﬁeld simulations: A dynamic stability study,” Antennas and Propagation, IEEE Transactions on, vol. 49, no. 6, pp. 902–907, 2001. [70] B. Donderici and F. L. Teixeira, “Conformal perfectly matched layer for the mixed ﬁnite 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 ﬁsh-eye lens antennas,” IEEE Trans. Ant. Prop., vol. 56, pp. 3058–3062, Sept. 2008. [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. 123 [77] J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic ﬁelds,” Science, vol. 312, 2006. [78] A. Alu, “First-principle homogenization theory for periodic metamaterials,” Phys. Rev. B, vol. 84, 2011. [79] I. Tsukerman, “Eﬀective 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, Taﬂove, 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. 1994. [88] V. A. Thomas, M. E. Jones, M. Piket-May, A. Taﬂove, 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. 124 [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 ﬁnite-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 ﬁnite 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 ﬁnite element method for the time-dependent ﬁrst order Maxwell’s equations,” Antennas and Propagation, IEEE Transactions on, vol. 58, pp. 4085–4090, Dec. 2010. [96] P. Li and L. J. Jiang, “A hybrid electromagnetics-circuit simulation method exploiting discontinuous Galerkin ﬁnite 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 ﬁnite 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 indeﬁnite time-harmonic Maxwell’s equations,” Numer. Math., vol. 100, p. 485518, Mar. 2005. 125 [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 modiﬁed 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. . http://www.synopsys.com, 2008. [106] Ngspice Users Manual. http://ngspice.sourceforge.net, 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 ﬁtting,” 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 buﬀer 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 ﬁrst approach to IBIS models: What they are and how they are generated,” AN-715Application Note, Analog Device, 2004. 126 [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. 127

1/--страниц