Finite Element Methods for Maxwell’s Equations Peter Monk Department of Mathematical Sciences University of Delaware Research funded by AFOSR and NSF. 75 years of Math. Comp. Peter Monk (UD) FEM for Maxwell MC-75 1 / 36
Table of Contents 1 Introduction A model scattering problem Variational formulation Continuous FEM Edge elements 2 Solar cells – Photonics Derivation of the model Discretization Implementation and Results Peter Monk (UD) FEM for Maxwell MC-75 2 / 36
Table of Contents 1 Introduction A model scattering problem Variational formulation Continuous FEM Edge elements 2 Solar cells – Photonics Derivation of the model Discretization Implementation and Results Peter Monk (UD) FEM for Maxwell MC-75 3 / 36
Time Harmonic Maxwell’s Equations E : Electric field and H : Magnetic field (both complex vector valued functions of position) The linear time harmonic Maxwell system at angular frequency ω > 0 is: − i ωǫ E + σ E − ∇ × H = − J , − i ωµ H + ∇ × E = 0 , where ǫ is the electric permittivity, µ the magnetic permeability, σ the conductivity and J is the applied current density. We assume µ = µ 0 > 0 and solve for E : � ǫ � + i σ ω 2 µ 0 ǫ 0 E − ∇ × ∇ × E = − i ωµ 0 J ǫ 0 ǫ 0 ω Define the wave number κ = ω √ ǫ 0 µ 0 and complex relative permitivitty � � ǫ 0 + i σ ǫ ǫ r = . In applications ǫ r := ǫ r ( x , ω ). ǫ 0 ω Peter Monk (UD) FEM for Maxwell MC-75 4 / 36
Scattering Problem E i : Known incident field E s : Scattered electric field E : Total electric field κ > 0 : Wave-number ǫ r = 1: Outside D is air Data (plane wave): E i = p exp( i κ d · x ) ( d ⊥ p , | d | = 1) Known incident field: Equations (no source, ǫ r = 1): ∇ × ∇ × E − κ 2 E = 0 in R 3 \ D Maxwell’s Equations: E = E i + E s in R 3 \ D Total field: Boundary Conditions: ( ν is unit outward normal) Perfect Electric Conductor (PEC): ν × E = 0 on Γ lim r →∞ (( ∇ × E s ) × x − i κ r E s ) = 0. Silver-M¨ uller Radiation Condition: For a bounded Lipschitz domain D with connected complement, this problem is well posed. Peter Monk (UD) FEM for Maxwell MC-75 5 / 36
Reduction to a Bounded Domain Introduce a closed surface Σ containing the scatterer (e.g. a sphere of radius R large enough) : Let Ω be the region inside Σ and outside D Need an appropriate boundary condition on the artificial bound- ary Σ. Peter Monk (UD) FEM for Maxwell MC-75 6 / 36
Absorbing Boundary Condition (simplest case) The simplest method is to ap- ply the Silver-M¨ uller Radiation Condition on Σ. Note ν is the outward normal on Σ. Abusing notation, E now de- notes the approximate field on the truncated domain. Incident Field: E i ( x ) = p exp( i κ x · d ), p · d = 0, and | d | = 1 ∇ × ( ∇ × E ) − κ 2 E = 0 in Ω ∇ × E i � × ν − i κ E i � ( ∇ × E ) × ν − i κ E T = T on Σ E × ν = 0 on Γ Here E T = ( ν × E ) × ν is the tangential trace of E . Peter Monk (UD) FEM for Maxwell MC-75 7 / 36
Function Spaces u ∈ ( L 2 (Ω)) 3 | ∇ × u ∈ ( L 2 (Ω)) 3 � � H ( curl ; Ω) = u ∈ H ( curl ; Ω) | ν × u | Σ ∈ ( L 2 (Σ)) 3 , ν × u = 0 on Γ � � = X with norms � � u � 2 ( L 2 (Ω)) 3 + �∇ × u � 2 � u � H ( curl ;Ω) = ( L 2 (Ω)) 3 � � u � 2 H ( curl ;Ω) + � ν × u � 2 � u � X = ( L 2 (Σ)) 3 Let � � ( u , v ) = u · v dV , � u , v � Σ = u · v dA . Ω Σ To avoid complications, from now on we assume Ω has two connected boundaries, Σ, Γ, and D has connected complement. Peter Monk (UD) FEM for Maxwell MC-75 8 / 36
Galerkin Method Multiply by the conjugate of a test function φ such that φ × ν = 0 on Γ, and integrate over Ω: � ∇ × ( ∇ × E ) − κ 2 E � � 0 = · φ dV Ω � � ∇ × E · ∇ × φ − κ 2 E · φ dV + = ν × ∇ × E · φ dA . Ω ∂ Ω The boundary terms are replaced using boundary data or the vanishing trace on Γ: � � ν × ∇ × E · φ dA = − i κ E T · φ dA ∂ Ω Σ � ν × ∇ × E i − i κ E i � � + · φ dA . T Σ We arrive at the variational problem of finding E ∈ X such that ( ∇ × E , ∇ × φ ) − κ 2 ( E , φ ) − i κ � E T , φ T � Σ = � F , φ T � Σ for all φ ∈ X where F = ( ∇ × E i ) × ν + i κ E i T . Peter Monk (UD) FEM for Maxwell MC-75 9 / 36
Existence and Approximation Problem: curl has a large null space. We use the Helmholtz decomposition: p ∈ H 1 (Ω) | p = 0 on Γ , p constant on Σ � � Define: S = then ∇ S ⊂ X . Choosing φ = ∇ ξ , ξ ∈ S as a test function: − κ 2 ( E , ∇ ξ ) = 0 so E is divergence free. Next we prove uniqueness. Then, using the subspace of divergence free functions ˜ X ⊂ X , the compact embedding of X in L 2 and the Fredholm alternative we have: ˜ Theorem Under the previous assumptions on D, there is a unique solution to the variational problem for any κ > 0 . If Σ is a sphere of radius R , and B is a fixed domain inside Σ and outside D then for R large enough � E true − E truncated � H ( curl ; B ) ≤ C R 2 Peter Monk (UD) FEM for Maxwell MC-75 10 / 36
Standard continuous elements [1980-90’s] Suppose that Ω has been covered by a regular tetrahedral mesh denoted by T h (tetrahedra having a maximum diameter h ). An obvious choice: use three copies of standard continuous piecewise linear finite elements. If we construct a finite element subspace X h ⊂ X using these continuous elements (note X h ⊂ ( H 1 (Ω)) 3 ), we find that the previously defined variational formulation gives incorrect answers due to lack of control of the divergence. Peter Monk (UD) FEM for Maxwell MC-75 11 / 36
Standard continuous elements [1980-90’s] Suppose that Ω has been covered by a regular tetrahedral mesh denoted by T h (tetrahedra having a maximum diameter h ). An obvious choice: use three copies of standard continuous piecewise linear finite elements. If we construct a finite element subspace X h ⊂ X using these continuous elements (note X h ⊂ ( H 1 (Ω)) 3 ), we find that the previously defined variational formulation gives incorrect answers due to lack of control of the divergence. We can try to add a penalty term. Choose γ > 0 sufficiently large and seek E h ∈ X h such that ( ∇ × E h , ∇ × φ h ) + γ ( ∇ · E h , ∇ · φ h ) − κ 2 ( E h , φ h ) − i κ � E h , T , φ h , T � Σ = � F , φ h , T � for all φ h ∈ X h Peter Monk (UD) FEM for Maxwell MC-75 11 / 36
A numerical analyst’s nightmare However, if Ω has reentrant corners, we may compute solutions that converge as h → 0 but to the the wrong answer! 1 The correct space for the problem is X N = X ∩ H ( div , Ω) but H N = H 1 (Ω) 3 ∩ X is a closed subspace of X N in the curl+div norm. If you want to use continuous elements, consult Costabel & Dauge ( γ needs to be position dependent) or more recently Bonito’s papers 2 . To handle this problem and discontinuous fields due to jumps in ǫ r , we can elec 3 (see also Whitney). use vector finite elements in H ( curl ) due to N´ ed´ 1M. Costabel, M. Dauge, Numer. Math. 93 (2002) 239-277. 2A. Bonito et al., Math. Model. Numer. Anal., 50, 1457–1489, 2016. 3J.C. N´ ed´ elec, Numer. Math. 35 (1980) 315-341. Peter Monk (UD) FEM for Maxwell MC-75 12 / 36
Finite Elements in H ( curl ) [N´ ed´ elec 1980, 1986] The lowest order edge finite element space is X h = { u h ∈ H ( curl ; Ω) | u h | K = a K + b K × x , a K , b K ∈ C 3 , � ∀ K ∈ T h . The degrees of freedom (unknowns) for this � element are e u h · τ h ds for each edge e of each tetrahedron where τ e is an appropriately oriented tangent vector. Note: N´ ed´ elec describes elements of all or- ders and in a later paper a second family of elements. 4 Engineering codes often use 2nd or higher order elements. 4P. Monk, Finite Element Methods for Maxwell’s Equations , Oxford University Press, 2003. Peter Monk (UD) FEM for Maxwell MC-75 13 / 36
Edge Element Method Let X h be the discrete space consisting of edge finite elements. We now seek E h ∈ X h such that ( ∇ × E h , ∇ × φ h ) − κ 2 ( E h , φ h ) − i κ � E h , T , φ h , T � Σ = � F , φ h , T � for all φ h ∈ X h . Peter Monk (UD) FEM for Maxwell MC-75 14 / 36
Discrete Divergence Free Functions Recall that if S = { p ∈ H 1 (Ω) p = 0 on Γ , p = constant on Σ } then ∇ S ⊂ X and this property enabled control of the divergence. Peter Monk (UD) FEM for Maxwell MC-75 15 / 36
Discrete Divergence Free Functions Recall that if S = { p ∈ H 1 (Ω) p = 0 on Γ , p = constant on Σ } then ∇ S ⊂ X and this property enabled control of the divergence. An important property of N´ ed´ elec’s elements is that they contain many gradients. In the lowest order case, if S h = { p h ∈ S | p h | K ∈ P 1 , ∀ K ∈ T h } , then ∇ S h ⊂ X h . Peter Monk (UD) FEM for Maxwell MC-75 15 / 36
Discrete divergence free We write a discrete Helmholtz decomposition X h = X 0 , h ⊕ ∇ S h . Functions in X 0 , h are said to be discrete divergence free . X 0 , h = { u h ∈ X h | ( u h , ∇ ξ h ) = 0 , for all ξ h ∈ S h } . Note X 0 , h �⊂ X 0 . Peter Monk (UD) FEM for Maxwell MC-75 16 / 36
Recommend
More recommend
