Broadband Finite-Element Impedance Computation for Parasitic Extraction
EElectrical Engineering manuscript No. (will be inserted by the editor)
Broadband Finite-Element Impedance Computation for ParasiticExtraction
J. Stysch · A. Klaedtke · H. De Gersem
Received: date / Accepted: date
Abstract
Parasitic extraction is a powerful tool in the de-sign process of electromechanical devices, specifically aspart of workflows that ensure electromagnetic compatibil-ity. A novel scheme to extract impedances from CAD de-vice models, suitable for a finite element implementation,is derived from Maxwell’s equations in differential form. Itprovides a foundation for parasitic extraction across a broadfrequency range and is able to handle inhomogeneous per-mittivities and permeabilities, making it more flexible thanexisting integral equation approaches. The approach allowsfor the automatic treatment of multi-port models of arbi-trary conductor geometry without requiring any significantmanual user interference. This is achieved by computing aspecial source current density from its given divergence onchosen terminal surfaces, subsequently using this currentdensity to compute the electric field, and finally calculat-ing the impedance via a scalar potential. A mandatory low-frequency stabilization scheme is outlined, facilitating thefinite element implementation. Two useful quasistatic ap-proximations and the advantageous special case of perfectelectric conductors are treated theoretically. The magneto-quasistatic approximation is validated against an analyticalmodel in a numerical experiment. Moreover, the intrinsic ca-pability of the method to treat inhomogeneous permittivitiesand permeabilities is demonstrated with a simple capacitor-coil model including dielectric and magnetic core materials.
Keywords
Parasitic effects · Finite element method · Quasistatics
J. Stysch · A. KlaedtkeRobert Bosch GmbH71272 Renningen, GermanyE-mail: [email protected]. De GersemInstitute for Accelerator Science and Electromagnetic FieldsTechnical University of Darmstadt
Increasing switching frequencies in power electronics, minia-turization and stricter electromagnetic compatibility (EMC)regulations pose a big challenge for the development of elec-tromechanical devices. Numerical simulations used in de-sign and optimization workflows are commonplace. Simula-tions do not only allow to assess a design before prototypingbut also to analyze the electromagnetic behavior of modelsthat are difficult to access in measurements due to, e.g., verycompact dimensions. An important task in EMC analysis isto quantify the influence of the non-functional elements likeinterconnects or chassis of a device under test (DUT) on itsintended functionality. This is typically done by extracting parasitic components from a model of the DUT in order tocomplement circuit models that contain the functional ele-ments. Parasitic components can either be simply frequency-independent resistances, inductances, and capacitances [1]or frequency-dependent (reduced order) models of transferfunctions like the impedance, admittance and scattering ma-tricies Z , Y , and S , respectively [2].The partial element electric circuit (PEEC) method [1]and further developments thereof [2, 3] are classical solutiontechniques in this context. This class of methods is based onMaxwell’s equations in integral form, using Green’s func-tions in the solution process, which forbids any straightfor-ward treatment of spatially inhomogeneous permittivities ε or permeabilities µ . In their standard form the methods ofthis class require a spatial discretization into cuboid volumeelements, which can be impractical to represent the complexgeometries encountered in industrial applications. General-izing these methods with less restrictive volume elementsmay be cumbersome [4].These disadvantages are avoided with an approach basedon the finite element (FE) method, which allows for a highlyflexible spatial discretization (commonly using tetrahedral a r X i v : . [ c s . C E ] S e p J. Stysch, A. Klaedtke, H. De Gersem elements), and an inherent treatment of inhomogeneous ma-terial parameters. In [5] an FE-based method for imped-ance computation was introduced in a rudimentary fashion,lacking an appropriate theoretical treatment, and consider-ing only the special case of homogeneous material param-eters and lossless conductors in the Darwin approximation.Nevertheless, a successful application of this method servedas the foundation of a geometrical sensitivity analysis in [6]and confirmed its potential. The aim of this paper is to pro-vide a solid theoretical foundation for a more general imped-ance computation method based on the approach outlinedin [5]. The such computed impedance matrix Z may subse-quently be used to extract parasitic lumped elements or tocompute the Y and S transfer functions.Section 2 describes how the impedance of a conductorsegment can be calculated from field quantities by introduc-ing a reduced voltage, which can ultimately be calculatedfrom the electric scalar potential φ . The main section 3 pro-vides the derivation of the field-theoretical model determin-ing φ via the E-field formulation of Maxwell’s equations.A suitable excitation current density is introduced, and anyunwanted inductive influence of this excitation current is de-embedded from the result for φ by a compensation term.Two quasistatic approximations are provided, and the spe-cial case of lossless conductors for a high-frequency approx-imation of the inductance is discussed. Finally, a necessarylow-frequency stabilization scheme is outlined. In section 4numeric results of the proposed method are compared to an-alytical impedance and inductance values of a wire, and thecapability of the method to handle inhomogeneous permit-tivities ε and permeabilities µ is demonstrated with a modelfeaturing a dielectric and magnetic core. A brief summaryconcludes this paper in section 5. To reconcile the path-dependent voltage concept of elec-tromagnetic field theory with the path-independent voltageconcept of electrical circuits can pose a challenge. In thissection, we consider the example of a one dimensional wireas discussed in [7, Chapter 14.16] and heuristically general-ize the result to three-dimensional conductors.Fig. 1 displays a path c along a wire, which ends in theterminals T a and T b . A return path r connecting the two ter-minals forms a closed loop with c , that encloses the surface S . Integrating Faraday’s law in frequency domain over S and applying Stokes’ theorem yields (cid:73) ∂ S EEE · dlll = (cid:90) r EEE · dlll + (cid:90) c EEE · dlll = − j ω (cid:90) S BBB · dSSS , (1)with EEE and
BBB denoting the electric field strength and themagnetic flux density, respectively. The voltage V that has cT a rT b Fig. 1: Wire c with terminals T a and T b and return path r .to be applied at the terminals in order to move charges from T a to T b against the self-induced electric field EEE must be thenegative of the integral over r , which yields V = − (cid:90) r EEE · dlll = (cid:90) c EEE · dlll + j ω (cid:90) S BBB · dSSS . (2)This is the voltage an ideal voltmeter would measure be-tween the terminals. The first integral on the right-hand sidecaptures the voltage drop due to the internal impedance ofthe wire (i.e., for a one-dimensional wire, the ohmic resis-tance) while the second integral is associated with the ex-ternal reactance of the loop [7]. For frequencies sufficientlybelow the first resonance of the system, the surface integralover BBB divided by the current I causing the magnetic field isthe inductance of the loop, L loop = I (cid:90) S BBB · dSSS . (3)However, for the purpose of extracting parasitics to be usedin circuit simulation, we are not interested in the inductiveresponse of the closed loop but only in the part of the re-sponse due to the wire. In an experiment, the return currentpath r closing the loop would correspond to either a currentsource or a voltmeter connected to the terminals. We there-fore want to exclude the part of the inductive response due tothe path r between the terminals and, moreover, avoid hav-ing to specify the path closing the loop altogether.This is achieved by employing the concept of partial in-ductance (see, e.g., [8]): The partial inductances L c and L r of the respective segments of the loop can be calculated byexpressing the magnetic field with the vector potential AAA as BBB = curl AAA and applying Stokes’ theorem: L loop = I (cid:73) ∂ S AAA · dlll = I (cid:90) c AAA · dlll + I (cid:90) r AAA · dlll = : L c + L r . (4)The values of these partial inductances are depend on thegauge condition chosen for the magnetic vector potential AAA and electric scalar potential φ .Subtracting from the voltage of (2) the inductive con-tribution j ω IL r of the path between the terminals allows to roadband Finite-Element Impedance Computation for Parasitic Extraction 3 define a reduced “conductor voltage” V c that contains onlythe resistive and inductive response of the wire c , V c : = V − I j ω L r = (cid:90) c EEE · dlll + j ω (cid:90) c AAA · dlll . (5)Expressing the electric field in (5) as EEE = − grad φ − j ω AAA (6)yields the simplification V c = − (cid:90) c grad φ · dlll = φ ( T b ) − φ ( T a ) . (7)Thus, the reduced voltage is in fact just the potential dif-ference of the scalar potential φ between the two terminalsand thereby path independent. As it incorporates a partialinductance, V c formally depends on the gauge condition ofpotentials φ and AAA .The path independence of V c facilitates a heuristic gen-eralization to three-dimensional conductors, where the ter-minals are surfaces instead of points. The reduced voltagecan be calculated analogously to the one-dimensional caseof (7) by simply averaging the potential over the respectivesurfaces, V c : = A ( T b ) (cid:90) T b φ dS − A ( T a ) (cid:90) T a φ dS . (8)Here, A ( T a ) and A ( T b ) denote the surface areas of the twoterminal surfaces.In the general case, there may be several conductors thatmay each have multiple terminals. To calculate an N × N im-pedance matrix Z , a topology of N branches connecting theterminals must be provided. A series of N numerical experi-ments can then be conducted, in each of which the current I flows through one of the branches j . The element Z i j of theinductance matrix is given through the voltage V i j betweenthe two terminals of branch i (which is calculated using (8)), Z i j = V i j I . (9) Z iscalculated from the electric scalar potential φ . For a systemexcited by a given source current density JJJ s , two possiblestrategies to calculate φ are available:The first option is to solve Maxwell’s equations directly in apotential formulation for the electric scalar potential φ andthe magnetic vector potential AAA . The second option is to first calculate
EEE by solving the ‘E-field formulation’ (see, e.g,[9]), and subsequently calculate φ in a second step using (6)and a gauge condition. It is more advantageous to use the lat-ter “ EEE approach” for several reasons: The two fields
EEE and φ can be calculated in sequence (except for the Darwin ap-proximation case discussed in section 3.4), thereby avoidinga computationally more expensive coupled boundary valueproblem (BVP), which occurs in the “ AAA - φ approach”. Fur-thermore, the EEE approach allows for an easy treatment ofconductors modeled as perfect electric conductors (PECs),which is an important special case discussed in section 3.6.Finally, the use of the E-field formulation allows for an el-egant stabilization of the low-frequency instability inherentto all FE solutions of Maxwell’s equations in frequency do-main, which is discussed in section 3.7.The E-field formulation can be derived combining Am-p`ere’s law and Faraday’s law. To provide a complete bound-ary value problem (BVP) that forms the basis of a FE so-lution, the boundary ∂ Ω of the computational domain Ω isassumed to be the union of an electric and a magnetic bound-ary, Γ el and Γ mag , on which electric (Dirichlet) and magnetic(Neumann) boundary conditions, respectively, are to be ap-plied: ∂ Ω = Γ el ∪ Γ mag . (10)Commonly either Γ el or Γ mag is empty, such that ∂ Ω is en-tirely electric or magnetic. The BVP of the E-field formula-tion thereby reads,curl ν r curl EEE + j ω µ σ EEE − ω c ε r EEE = − j ω µ JJJ s in Ω , (11a) nnn × EEE = Γ el , (11b) nnn · ε r EEE = Γ mag . (11c)Here, σ denotes the electric conductivity, µ permeabilityof the vacuum, and c the speed of light, while ε r and ν r are the spatially variable relative permittivity and relativereluctivity, respectively, and nnn is the normal vector on ∂ Ω .The magnetic boundary condition for EEE (11c) can be derivedfrom the magnetic boundary condition for the magnetic fieldstrength
HHH , nnn × HHH =
0, using Amp`ere’s law and demanding nnn · ( σ EEE + JJJ s ) = Γ mag .The partial differential equation (PDE) to determine φ isfound by using (6) to eliminate the vector potential AAA froma gauge condition. Here, the Lorenz gauge [10],div ε r AAA + j ω c φ = , (12)is chosen since it enables the calculation of the inductivecompensation term introduced in section 3.3. Together (6)and (12) yield the PDEdiv ε r grad φ + ω c φ = − div ε r EEE (13)
J. Stysch, A. Klaedtke, H. De Gersem
The full BVP to compute the electric scalar potential is givenin section 3.3, after modeling the source current density
JJJ s and introducing the compensation term.3.2 Modeling the Source Current DensityThe remaining task is to identify a suitable source currentdensity JJJ s modeling a current source connected to two of theterminals of the DUT. The source current density JJJ s shouldinject a constant current I at the terminal surface T b and ex-tract the same current again at T a , such that it flows alongthe integration path c in Fig. 1. Inside the conductor, thiscurrent I is transported between the terminals by a conduc-tion current density JJJ c , which unlike JJJ s directly couples to EEE via Ohm’s law,
JJJ c = σ EEE , and thereby captures the resis-tive response to the enforced current flow. For the currentdensities
JJJ a and JJJ b on the terminal surfaces homogeneousdistributions are chosen, JJJ a = ˆ nnn a I A ( T a ) and JJJ b = − ˆ nnn b I A ( T b ) , (14)with ˆ nnn i denoting the unit normal vector pointing out of theconductor, and A ( T i ) again the area of the respective termi-nals. The associated divergence of JJJ s must hence be givenbydiv JJJ s = | JJJ a | δ ( ) on T a , −| JJJ b | δ ( ) on T b , , (15)with the delta distribution δ . This expression describes thesituation that at the terminals the source current density J s takes over the task of transporting the current I from theconduction current density JJJ c .Our method to calculate JJJ s has to be so general that itcan easily produce a suitable JJJ s independently of how theterminal surfaces are positioned in relation to e.g. the outerboundary ∂ Ω or any of the conducting areas of the DUT.Generally, an expression for the divergence of the vectorfield JJJ s is not sufficient to determine JJJ s . However, choos-ing a gradient field ansatz, JJJ s = grad ξ , (16)enables a computation of JJJ s from its divergence given in(15) without further specifying the path of the source cur-rent. The BVP determining the underlying potential ξ (andtherefore JJJ s ) is then given bydiv grad ξ = div JJJ s in Ω , (17a) ξ = const . on Γ el , (17b) nnn · grad ξ = Γ mag . (17c)The boundary conditions (17b) and (17c) are chosen in thisway for consistency with the BVP (11) determining EEE . 3.3 Compensating the Inductive Influence of the SourceCurrent DensityThe gradient-field source current density proposed in theprevious subsection does not model an ideal current sourcein one respect: Ideal current sources must not influence theDUT inductively. In section 2, a reduced voltage V c was de-fined in which the part of the inductive response related tothe return path of the current was eliminated. In addition, theelectromagnetic fields causing V c should not capture any in-ductive influence of the source current JJJ s . The electric fieldcalculated with (11) includes this unwanted influence of thesource current, since any source current that has a compo-nent parallel to the DUT at a finite distance must be ex-pected to have a direct inductive influence on the fields in theDUT (such parallel components are unavoidable for three-dimensional conductors of arbitrary shape).It is, however, possible to quantify and eliminate (andthereby de-embed) the contribution of this unwanted directinfluence of the source current density in the calculation ofthe scalar potential with (13). To this end, the “total” elec-tric field EEE of (11), incorporating inductive effects relatedto both the conductor currents and the source current, is ex-pressed as the difference of the compensated field
EEE c cap-turing only the influence of the conductors and the field EEE s capturing the counteractive inductive influence of the sourcecurrent density, EEE = EEE c − EEE s . (18)The scalar potential to be calculated with (13) must be thecompensated potential φ c excluding the unwanted inductiveinfluence of JJJ s ,div ε r grad φ c + ω c φ c = − div ε r EEE c = − div ε r EEE − div ε r EEE s . (19)The compensation term div ε r EEE s can be determined employ-ing the potential formulation of Maxwell’s equations in Lorenzgauge,curl ν r curl AAA − ε r grad div ε r AAA − ω c ε r AAA = µ JJJ , (20a)div ε r grad φ + ω c φ = − ε ρ . (20b)This formulation of Maxwell’s equations enables the inde-pendent calculation of the vector potential AAA from the cur-rent density
JJJ and the scalar potential φ from the chargedensity ρ in the case that all currents in a given system aresource currents, JJJ = JJJ s , that do not couple to the electricfield via Ohm’s law.To calculate div ε r EEE s , the non-physical (since not chargeconserving) situation JJJ = JJJ s and ρ = roadband Finite-Element Impedance Computation for Parasitic Extraction 5 the conduction current density JJJ c is disregarded. The cor-rection term calculated from these sources therefore cap-tures the isolated effects of the source current density JJJ s .By only considering source currents and no source charges,the “source” scalar potential φ s associated with EEE s vanishesand the correction term only depends on the “source” vectorpotential AAA s ,div ε r EEE s = − j ω div ε r AAA s (21)A boundary value problem to calculate the scalar field div ε r AAA s directly is found by applying the divergence operator to (20a)div ε r grad div ε r AAA s + ω c div ε r AAA s = − µ div JJJ s , (22)and supplementing the same boundary conditions on the outerboundary ∂ Ω as for the computation of ξ in (16). The termdiv JJJ s is given with (15). To simplify the notation in the fol-lowing, the scalar field g : = − µ div ε r AAA s (23)is defined. Expressing (22) with g and supplying boundaryconditions yields the BVPdiv ε r grad g + ω c g = div JJJ s in Ω , (24a) g = const . on Γ el , (24b) nnn · ε r g = Γ mag . (24c)Using the definition of g in (19) and supplementing bound-ary conditions yields the BVP determining the compensatedscalar potential φ c ,div ε r grad φ c + ω c φ c = − div ε r EEE − j ω µ g in Ω , (25a) φ c = const . on Γ el , (25b) nnn · ε r grad φ c = Γ mag . (25c)This concludes the derivation of the field theoretical model.Thus, there are in total three steps to compute the compen-sated scalar potential φ c needed for the impedance calcula-tion:1. Compute ξ and g from the div JJJ expression of (15) using(17) and (24), respectively.2. Compute
EEE from
JJJ s = grad ξ using (11).3. Compute φ c from EEE and g using (25).This procedure is illustrated in Fig 2.The such computed fields are full-wave solutions con-sidering all effects of Maxwell’s equations, including re-tardation. Considering wave effects in a FE context gener-ally requires a strategy preventing reflections from the outer input:div JJJ s compute ξ compute gJJJ s = grad ξ compute EEE compute φ c Fig. 2: Steps for the computation of φ c .boundary of the finite computational domain, usually em-ploying perfectly matched layers (see, e.g., [11]). However,if the structures of a DUT are small compared to the wave-lengths associated with frequencies relevant to its EMC anal-ysis, wave effects can generally be disregarded. To not un-necessarily complicate the FEM implementation and to avoida possible greater numerical cost due to a finer space dis-cretization needed for the layers, wave solutions can in suchcases already be eliminated on the level of the underlyingPDEs by employing the quasistatic Darwin approximationintroduced in the following section 3.4.For applications that only require to consider inductiveeffects and ohmic losses the more restrictive magnetoqua-sistatic approximation introduced in section 3.5 is applica-ble, further reducing the complexity and numerical cost ofthe implementation.3.4 Darwin ApproximationDarwin’s approximation [12], sometimes also referred to as“full quasistatics”, neglects wave solutions while still fullycapturing resistive, inductive, and capacitive effects. In thissense it behaves as the “natural” field-theoretical equivalentto the circuit model, which also precludes wave effects. Theapproximation arises from eliminating the magnetic vectorpotential contribution in both the displacement current andin Gauss’ law, and is hence dependent on the gauge con-dition of the potentials. In the Lorenz gauge, it amounts toneglecting the ω -terms in (20). Hence, also the ω -termin (24a) disappears in Darwin’s approximation, yielding the J. Stysch, A. Klaedtke, H. De Gersem frequency-independent BVPdiv ε r grad g = div JJJ s in Ω , (26a) g = const . on Γ el , (26b) nnn · ε r g = Γ mag . (26c)A comparison with (17) determining the potential ξ of thesource current JJJ s shows that in Darwin’s approximation g can be used to express the source current, JJJ s = grad ξ = ε r grad g . (27)This is a great advantage if a model has to be evaluated atseveral frequency points, since in the Darwin case only oneBVP must be solved to calculate both JJJ s and g , in contrast to1 + N f BPVs in the non-approximated case, with N f beingthe number of frequency points to be evaluated.As the magnetic vector potential contribution to the dis-placement current is neglected in the Darwin approximation,the E-field formulation changes tocurl ν r curl EEE + j ω µ σ EEE + ω c ε grad φ = − j ω µ JJJ s . (28)Hence, it cannot be solved independently but only in a cou-pled BVP together with (25), which readscurl ν r curl EEE + j ω µ σ EEE + ω c ε r grad φ c = − j ω µ ε r grad g in Ω , (29a)div ε r EEE + div ε r grad φ c + ω c φ c = − j ω µ g in Ω , (29b) nnn × EEE = φ c = const . on Γ el , (29c) nnn · ε r EEE = nnn · ε r grad φ c = Γ mag . (29d)3.5 Magnetoquasistatic ApproximationThe magnetoquasistatic (MQS) approximation neglects thedisplacement current altogether and thereby also precludescapacitive effects in addition to wave effects. In contrast toDarwin’s approximation, the ω -terms of both (29a) and(29b) are neglected in the MQS formulation proposed here;while the first term is part of the displacement current, thesecond term results from the Lorenz gauge condition. Ne-glecting the latter is equivalent to choosing the Coulombgauge condition for the calculation of φ . This step is taken inthe MQS approximation to establish a consistency at higher-frequencies between the two PDEs, which in numerical ex-periments proves to be necessary to obtain plausible results at arbitrarily high frequencies. The BVP to determine EEE inthe MQS approximation is hence given bycurl ν r curl EEE + j ω µ σ EEE = − j ω µ ε r grad g in Ω , (30a) nnn × EEE = Γ el (30b) nnn · ε r EEE = Γ mag . (30c)The compensated scalar potential φ c is subsequently calcu-lated with the BVPdiv ε r grad φ c = − div ε r EEE − j ω µ g in Ω , (31a) φ c = const . on Γ el , (31b) nnn · ε r grad φ c = Γ mag . (31c)As in Darwin’s approximation, for (30) and (31) the scalarfield g is determined with the frequency-independent BVP(26).3.6 PEC CaseDue to the skin and proximity effects, the parasitics rep-resented by the impedances Z extracted using the general-case equations of the previous subsections are dependent onthe frequency in a general way. For some applications ofEMC analysis, however, it may suffice to provide the high-frequency limit which corresponds to a particular frequencydependence and to a decoupled pair of a frequency-indepen-dent inductance matrix L and a frequency-independent ca-pacitance matrix C . Approximating the parasitic effects withsimple frequency-independent lumped elements enables anespecially straightforward combination with the functionalelements in a joint circuit and a convenient and fast simula-tion thereof.It is possible to extract the high-frequency inductanceat a moderate numerical cost by modeling the conductorsof the DUT as perfect electric conductors (PECs). This en-forces a fully developed skin effect in the conductors suchthat the electric field EEE vanishes in the conducting domains.Then, at any frequency (below the first resonance) the sameinductance matrix, i.e. its high-frequency limit, is obtained(Fig. 4a).The computational cost of this approach is much lowerthan for the general case for two reasons: First, since nofrequency dependent behavior occurs, the BVPs have to besolved at only a single frequency point. Second, it is muchcheaper to solve the E-field formulation (11) if the conduc-tors are PECs since in this case (11a) is only enforced in thenon-conducting region Ω = Ω \ Ω c , supplemented by elec-tric boundary conditions for EEE on the boundary of the con-ducting region Ω c . Furthermore, the ohmic loss term j ω µ σ EEE disappears from the E-field form, such that the whole systemof equations only needs to be solved for the imaginary part roadband Finite-Element Impedance Computation for Parasitic Extraction 7 of the fields
EEE and φ c , since Re ( EEE ) = ( φ c ) = JJJ s is chosen real. This leads to purely real operator matricesafter FE discretization. In the PEC case, Darwin’s approximation captures induc-tive and capacitive effects while precluding ohmic lossesand wave effects. After solving (26) for g in the full domain Ω , the compensated potential φ c is computed with the BVPcurl ν r curl EEE + ω c ε r grad φ c = − j ω µ ε r grad g in Ω , (32a)div ε r EEE + div ε r grad φ c + ω c φ c = − j ω µ g in Ω , (32b) nnn × EEE = ∂ Ω c ∪ Γ el , (32c) φ c = const . on Γ el , (32d) nnn · ε r EEE = nnn · ε r grad φ c = Γ mag . (32e)Note that while the vectorial equation (32a) is only enforcedin the non-conducting subdomain Ω , the scalar equation(32b) is enforced in the full domain Ω . Modeling the conductors as PECs in the MQS approxima-tion yields a system that captures only inductive effects whileprecluding capacitive effects, ohmic losses, and wave ef-fects. The corresponding BVP to determine
EEE is given bycurl ν r curl EEE = − j ω µ ε r grad g in Ω , (33a) nnn × EEE = ∂ Ω c ∪ Γ el , (33b) nnn · ε r EEE = Γ mag , (33c)It is enforced only in the non-conducting subdomian Ω .The BVPs determining g and φ c , (26) and (31), respectively,remain unchanged and are enforced in the full domain Ω .Capturing only inductive effects, the reactance calcu-lated with (8) and (9) from this φ c has a linear frequencydependence and its corresponding inductance is frequencyindependent (Fig. 4a). A real, frequency-independent set ofequations to determine the constant inductance of the MQSPEC case directly can therefore be obtained by dividing theBVPs (33) and (31), and (8) and (9) by j ω (introducing thescaled fields EEE (cid:48) = EEE / j ω and φ (cid:48) c = φ c / j ω ). 3.7 Finite Element DiscretizationWe discretize the BVPs of the previous subsections by ex-pressing the scalar fields ξ , g , and φ c with the H -conform-ing and the vector field EEE with the H ( curl ) -conforming ba-sis functions given in [13]. To this end, the domain Ω ismeshed with tetrahedral elements. Testing in a Galerkin ap-proach the scalar and vectorial PDEs with the scalar and vec-torial basis functions, respectively, and integrating over thedomain Ω discretizes the BVPs into sparse linear systemsof equations.FE discretizations of the frequency-domain Maxwell equa-tions are notorious for resulting in singular operator ma-tricies at lower frequencies [14, 15, 16]. Generally, a low-frequency stabilization scheme is necessary to ensure thatthe associated FE method linear system is solvable at all fre-quencies. In [16] Eller et al. derived a stable weak formula-tion based on the E-field formulation (11). It is very generalin its scope, such that it can be applied to the Darwin andMQS approximations, and to the PEC case in a straightfor-ward way.The basic approach of this method to split the Sobolevspace H ( curl , Ω ) of the trial and test functions of EEE intothree parts H ( curl , Ω ) = V ⊕ W ⊕ U (34)with V : = { v ∈ H ( curl , Ω ) : curl v (cid:54) = } , (35) W : = { w ∈ H ( curl , Ω ) : curl w = , w (cid:54) = Ω c } , (36) U : = { u ∈ H ( curl , Ω ) : curl u = , u = Ω c } . (37)The electric field is thereby decomposed into one part sim-ilar to a vector potential, EEE V ∈ V , and two parts that are es-sentially gradient fields produced by scalar potentials, EEE W ∈ W and EEE U ∈ U , the latter of which vanishes in the conduct-ing subdomain Ω c , EEE = j ω EEE V + ( j ω ) / EEE W + EEE U . (38)For edge-based FE basis functions of the lowest order p = V are found using atree-cotree split [17]. The higher orders of the hierarchicalbasis functions of [13] are already split into a part with non-vanishing curl and a gradient field part, and thereby give thehigher-order functions of V explicitly. J. Stysch, A. Klaedtke, H. De Gersem the academic example of a straight wire of circular crosssection. An analytic expression for the internal partial im-pedance Z int of a wire of length l and radius r is given in[18], Z int = R + j ω L int = jl π r (cid:114) ω µσ (cid:18) Ber ( q ) + j Bei ( q ) Ber (cid:48) ( q ) + j Bei (cid:48) ( q ) (cid:19) (39)with q = r √ ω µσ , and Ber, Bei, Ber (cid:48) and Bei (cid:48) denoting thereal and imaginary Kelvin functions and their derivatives,respectively. This expression captures the ohmic resistance R = Re ( Z ) and the contribution of internal partial inductance L int to the reactance X = Im ( Z ) . It can be combined with theexpression for the external partial inductance L ext of a wiregiven in [8], L ext = µ l π (cid:32) arsinh (cid:18) lr (cid:19) − (cid:114) + (cid:16) rl (cid:17) + rl (cid:33) , (40)to produce an analytic approximation Z ana that describes theohmic resistance and the inductive contribution to the im-pedance of the wire, Z ana : = Z int + j ω L ext = R + j ω ( L int + L ext ) . (41)Neglecting both wave propagation and capacitive effects, Z ana is an MQS approximation.The example of a wire of length l =
50 mm and radius r = ∂ Ω = Γ el , and values producedwith a magnetic boundary, ∂ Ω = Γ mag , as suggested by thestrategic dual image technique [19], and the finite compu-tational domain Ω was chosen very large compared to themodel size.Fig. 3a compares the resistance over frequency calcu-lated both analytically with (41) and numerically with theMQS approximation of section 3.5 (the R values of the Dar-win approximation and the non-approximated system are thesame as the MQS values). The analytical and numerical val-ues are virtually identical until f = f >
10 MHz,however, the numerical values quickly approach a constantupper limit due to the FE space discretization becoming un-able to resolve the current density skin depth. This demon-strates the general infeasibility of numerical methods requir-ing a volumetric discretization for high-frequency resistancecomputations. In the context of EMC analysis, however, theresistive effects are generally only of minor importance.Fig. 3b compares the modulus of the reactance X of theanalytical values and both quasistatic approximations in ahigher frequency interval. While the analytical and numer-ical MQS values show the identical unchanging linear re-sponse, the Darwin approximation is capable of capturing resonant behavior resulting form the interplay of inductiveand capacitive effects. This illustrates on the one hand thatthe MQS approximation is suitable for precisely those ap-plications in which only the inductive effects are of interest,and on the other hand that the Darwin approximation maybe employed for a resonance analysis (which has been in-vestigated in [20]).In the MQS approximation, the inductance L is simply L = X / ω . Fig. 4 compares the inductance values of the nu-merical MQS approximation over frequency to the analyt-ical values calculated from (41). The values are in goodagreement; unlike for the resistance R in Fig. 3a the finiteminimal skin depth due to the space discretization does notlead to a major qualitative difference between the analyticaland numerical values at higher frequencies. The finite dis-cretization of the conductor in the numerical method man-ifests itself merely in a sightly altered curvature in the in-terval between 100 kHz and 10 MHz. The values of the PECcase equation system (33) displayed in Fig. 4a show howthis modification of the general case method produces thehigh-frequency limit of L directly at any frequency. Fig. 4bdemonstrates the effect of the inductive compensation termderived in section 3.3 by comparing the values of the regularmethod with compensation to values calculated without thecompensation term. The compensation leads to a frequency-independent shift of the inductance of ≈ .
85 nH. The un-compensated values are therefore on average ≈
13 % lowerthan both the analytical and compensated numerical values.4.2 Capacitor Coil ModelThe capability of the proposed FE-based method for imped-ance computation to handle inhomogeneous permittivities ε and permeabilities µ is demonstrated with a model featuringa parallel-plate capacitor with a dielectric, and a coil with aclosed magnetic core, as displayed in Fig. 5. Similar modelswithout the dielectric and magnetic core have been consid-ered by other authors, e.g. in [14].Modeling the conductors of the model as PECs as de-scribed in section 3.6 allows to extract a frequency-indepen-dent capacitance C for the capacitor and inductance L for thecoil. The capacitance C is calculated with the Darwin systemof section 3.6.1 by choosing the surfaces of the capacitorplates as terminal surfaces and computing the reactance X FE between the two plates at a suitably low test frequency f .The capacitance is then given by C = − / ( π f X FE ( f )) .As an example for such a capacitance extraction, themodel of Fig. 5 is considered with the choice ε r =
100 inthe dielectric and µ r = f =
100 Hz the capacitance is computed as roadband Finite-Element Impedance Computation for Parasitic Extraction 910 − − − − f [Hz] R [ Ω ] MQS approximationanalytical values(a) 10 f [Hz] | X | [ Ω ] MQS approximationanalytical valuesDarwin approximation(b)
Fig. 3: Frequency dependent resistance R (left) and modulus of the reactance X (right). The FEM space discretization severelylimits the applicability of the numerical method for high-frequency resistance computation. While the Darwin approximationis capable of capturing resonant behavior in the reactance, the response of the MQS approximation is only inductive andhence linear. f [Hz] L [ n H ] MQS approximationMQS PEC caseanalytical values(a) 10 f [Hz] L [ n H ] MQS w/ compensationMQS w/o compensationanalytical values(b)
Fig. 4: Frequency dependence of the inductance L . The numerical values of the MQS approximation are in good agreementwith the analytical values. The PEC case gives the high-frequency limit directly (left plot). The effect of the compensationterm derived in section 3.3 is demonstrated in the right plot.Fig. 5: Model featuring a capacitor with a dielectric and acoil with a closed magnetic core. The conducting material isdisplayed in yellow. The total length of the model is 40 mm . C ≈ .
184 pF. The relative error e X : = (cid:12)(cid:12)(cid:12)(cid:12) X FM − X C X FE (cid:12)(cid:12)(cid:12)(cid:12) (42)between the reactance X C : = − / ( π fC ) associated withthe extracted capacitance and the FE solution X FE is dis-played in Fig. 6. The error e X stays below 10 − for fre-quencies f < X C perfectly capturesthe behavior of the FE solution X FE and that the frequenciesof this range are suitable test frequencies for the extractionof the capacitance C . Above f = X FE becomes measurable. − − − − f [Hz] e rr o r e X Fig. 6: Frequency dependence of the relative error e X be-tween the reactance X C of the extracted capacitance C andthe FE solution X FE , as defined in (42).However, the relative error e X stays well below 10 − for fre-quencies as high as f = C of thecapacitor normalized by the relative permittivity ε r of thedielectric for a parameter sweep of ε r (with the permeabilityof the magnetic core fixed at µ r = ε r the normalized capacitance approaches a constant value C / ε r ≈ .
908 pF. This demonstrates the proportionality of C to relative permittivity of the dielectric that takes effect whenthe field lines of the electric field EEE begin to be localized inthe dielectric at values ε r > L ofthe coil normalized with the relative permeability µ r of thecore (the terminal surfaces being the right capacitor plateand the right end surface of the coil, with the permittivityof the dielectric fixed at ε r = L is computed directly with the MQS approx-imation as described in section 3.6.2. The behavior of thenormalized inductance L / µ r is qualitatively similar to thatof C / ε r in Fig. 7a: For large values of µ r the normalized in-ductance approaches a constant value L / µ r ≈ . L to µ r for values µ r (cid:29)
1, atwhich the field lines of the magnetic flux
BBB are concentratedin the core.
In his work, the theoretical foundation of an FE-based nu-merical impedance computation method was developed, pro-viding a basis for parasitic extraction for models of com-plex geometry. After clarifying the relationship of the im-pedance matrix Z and the scalar potential φ , a the field equa-tions for the calculation of φ were formulated and a suit-able source current density was modeled. A de-embedding approach was pursued to eliminate any unwanted inductiveinfluence of the current source on the result. The FE dis-cretization including a necessary low-frequency stabiliza-tion scheme was briefly discussed. While the derivation didnot rely on any approximation to Maxwell’s equations, twoquasistatic approximations especially relevant for applica-tion were also provided, as well as the important special caseof perfect electric conductors giving a high-frequency induc-tance approximation. The method was validated by compar-ing the quasistatic approximations and the special PEC caseto analytic values for a wire model. The results demonstratethat the finite minimal skin depth of the FE discretizationdoes not negatively influence the quality of high-frequencyinductance results in a substantial way. The capability of themethod to handle inhomogeneous permittivities and perme-abilities was demonstrated with a capacitor coil model fea-turing a dielectric and a magnetic core. References