Improved non-overlapping domain decomposition algorithms for the eddy current problem

0.64 MB·

by  Y. Boubendir

6 Improved non-overlapping domain decomposition algorithms for
1
0 the eddy current problem
2
n Y. Boubendira, H. Haddarb M. K. Riahia
a aDepartment of mathematical science, New Jersey Institute of Technology, New Jersey USA.
J
bINRIA, Ecole Polytechnique (CMAP) and Universit´e Paris Saclay, Route de Saclay, Palaiseau 91128 Cedex FRANCE.
1
3 Received*****;acceptedafterrevision+++++
PresentedbyAˆ£Aˆ£Aˆ£Aˆ£Aˆ£
]
A
N
h. Abstract
t
a A domain decomposition method is proposed based on carefully chosen impedance transmission operators for a
m hybrid formulation of the eddy current problem. Preliminary analysis and numerical results are provided in the
[ sphericalcaseshowingthepotentialoftheseconditionsinacceleratingtheconvergencerate. To cite this article:
Y. Boubendir, H. Haddar and M.K. Riahi, C. R. Acad. Sci. Paris, Ser. I 340 (2015).
1
v R´esum´e
4
Nousproposonsunem´ethodeded´ecompositiondedomainebas´eesurunchoixparticulierdel’´ecrituredescondi-
9
tions de transmission pour une formulation hybride du mod`ele courant de Foucault. Nous donnons des r´esultats
2
0 analytiques et num´eriques pr´eliminaires dans le cas sph´erique qui montrent le potentiel de ces conditions dans
0 l’acc´el´eration de la vitesse de convergence d’une r´esolution it´erative du probl`eme. Pour citer cet article : Y.
. Boubendir, H. Haddar et M.K. Riahi, C. R. Acad. Sci. Paris, Ser. I 340 (2015).
2
0
6
1
:
v
i
X Version franc¸aise abr´eg´ee
r
a
Le mod`ele de courant de Foucault s’´ecrit pour le champ ´electrique E et le champ magn´etique H

curlH−σE =0 dans Ω

curlE−iωµH =0 dans Ω (1)
n×H =n×J sur ∂Ω.
e
Email addresses: [email protected](Y.Boubendir),[email protected](H.Haddar),[email protected]
(M.K.Riahi).
URL: https://web.njit.edu/ riahi/(M.K.Riahi).
PreprintsubmittedtotheAcad´emiedessciences 2f´evrier2016
ou` la perm´eabilit´e magn´etique µ est une fonction `a valeurs r´eelles qui peut d´ependre de l’espace, et J
e
repr´esente un terme source. La conductivit´e´electrique σ est nulle en dehors de la partie conductrice Ω .
C
Ces equations sont compl´et´ees par les conditions d’interface traduisant la continuit´e des composantes
tangentielles de E et H `a travers Γ, l’interface entre Ω et le vide Ω :=Ω\Ω .
C I C
Lechampmagn´etiqueHest`arotationnelnuldansΩ .Ceciimplique,lorsquecedernierestsimplement
I
connexe, que H=∇p (voir par exemple [3]), et le mod`ele s’´ecrit sous la forme
(cid:40)
curlcurlE−κ2E=0 dans Ω
C (2)
∇·(µ∇p)=0 dans Ω
I
avec κ2 :=iωµσ. Les conditions de continuit´e `a l’interface Γ s’´ecrivent
−−→ ∂p 1
curlE×n=iωµcurl p = curl E sur Γ (3)
Γ ∂n iωµ Γ
−−→
ou` n est un vecteur normal unitaire sur Γ qui pointe vers l’ext´erieur de Ω et ou` curl p:=∇ p×n et
C Γ Γ
curl E:=div (E×n)avec∇ etdiv d´esignantrespectivementlegradientetladivergencesurfaciques.
Γ Γ Γ Γ
Nousremarquonsquelecasnonsimplementconnexepourrait´egalementˆetretrait´eaumoyenderajouter
`aplescontributionsd’unnombrefinidefonctionsharmoniques([3]).Danscetravailnousproposonsune
m´ethode de d´ecomposition de domaine bas´e sur la d´ecomposition de Ω en partie conductrice et le vide
mais en remplac¸ant les conditions d’interface naturelle par les conditions de continuit´e suivantes
∂p −−→ 1
+β curl curl p= curl (E+β curlE×n) (4)
∂n I Γ Γ iωµ Γ I
(cid:18) (cid:19)
−−→ −−→ ∂p
curlE×n+β curl (curlE·n)=iωµcurl p+β (5)
C Γ Γ C∂n
ou` β and β sont des nombres complex appropri´ees. On montre que ces conditions sont ´equivalentes
I C
aux conditions originales sous certaines conditions sur β et β (Lemme 1.1). Notre algorithme it´eratif
C I
est d´ecrit par (14)-(15). Dans la cas sph´erique, on montre que asymptotiquement, l’op´erateur d’it´eration
(16) est une contraction plus rapide que celle correspondant au choix β = β = 0 (Proposition 2.2).
I C
Nous concluons cette note par un exemple num´erique 3D montrant une meilleure vitesse de convergence
pour des valeurs de β et β diff´erentes de 0.
C I
1. Introduction
The eddy current approximation of the Maxwell equations, for the electric field E and the magnetic
field H reads

curlH−σE =0 in Ω

curlE−iωµH =0 in Ω (6)
n×H =n×J on ∂Ω.
e
where the magnetic permeability µ is a real valued function that may depend on space and J stands
e
for the source excitations. The electric conductivity σ has support only in the conductive material Ω .
C
These equations are supplemented with the continuity of tangential components of E and H across Γ,
the interface between conductive and non conductive regions.
The field H is curl free in the insulating region Ω := Ω\Ω . This implies in the case of a simply
I C
connected topology that the magnetic field is a gradient of a harmonic scalar field i.e. H = ∇p (see [3]
and references therein), which leads to the problem
2
(cid:40)
curlcurlE−κ2E=0 in Ω
C (7)
∇·(µ∇p)=0 in Ω .
I
where κ2 :=iωµσ. The interface continuity conditions across Γ lead to
−−→ ∂p 1
curlE×n=iωµcurl p = curl E on Γ (8)
Γ ∂n iωµ Γ
−−→
where n is a unitary normal vector on Γ pointing toward the exterior of the conductor. We use curl p:=
Γ
∇ p×n and curl E:=div (E×n) with ∇ and div respectively being the surface gradient and the
Γ Γ Γ Γ Γ
surface divergence. To simplify the presentation we also assume that ∂Ω⊂∂Ω and close the problem by
I
imposingaboundaryconditionp=f on ∂Ω.Letusalreadyremarkthatthecaseofnonsimplyconnected
domaincanalsobetreatedinasimilarwaybutwithadditionaltechnicalitiesrelatedtotheincorporation
of (finite number of) divergence and curl free functions.
Manymethods,suchaspotentialandhybridformulations[4,2,9],havebeendevelopedtodealwiththis
type of problem. An extensive overview of these methods can be found in [3]. Generally speaking, these
methods fall in the class of direct formulations gathering the problem in Ω and Ω into the same linear
C I
system to solve. A nice discussion of the advantages/disadvantages of such methods and others can be
alsofoundin[2],wheretheauthorsproposeanewnumericaldiscretizationtoovercomethesedifficulties,
mainly related to the size of the linear systems. Although the approach in [2] represents an important
contribution in solving these difficulties, it remains strongly dependent on the use of preconditioners
because of the ill-conditioning which is further worsen by the high contrast created by the conductivity
in Ω .
C
Domain decomposition methods [8,10] are well suited for this problem since they allow the problems
in Ω and Ω to be solved separately with appropriate approaches. Up to our knowledge, the only
C I
domaindecompositionalgorithmdevelopedforformulation(7)isthegivenin[3]andexploitstransmission
conditions (8). We here propose to improve this method by modifying these conditions under the form
∂p −−→ 1
+β curl curl p= curl (E+β curlE×n) (9)
∂n I Γ Γ iωµ Γ I
(cid:18) (cid:19)
−−→ −−→ ∂p
curlE×n+β curl (curlE·n)=iωµcurl p+β (10)
C Γ Γ C∂n
where β and β denote appropriate (complex) numbers. For similar ideas in different contexts we refer
I C
the reader to [5,7,6].
We establish in this part a consistency property for the above impedance boundary condition.
Lemma 1.1 The following conditions
Re{−β β }≥0, or Im{β β }=(cid:54) 0 (11)
C I C I
ensure consistency between the original conditions (8) and the new ones (9)-(10).
Proof. Let us define the following quantities on Γ
−−→ ∂p 1
Ξ :=curlE×n−iωµcurl p, Ξ := − curl E, (12)
C Γ I ∂n iωµ Γ
which are zero if interface conditions (8) hold. Interface conditions (10) and (9) can be written as
−−→ 1
Ξ +iωµβ curl Ξ =0, Ξ + β curl Ξ =0. (13)
C C Γ I I iωµ I Γ C
−−→
Therefore, Ξ − β β curl curl Ξ = 0. with implies, using a variational form, that (cid:107)Ξ (cid:107)2 −
C C I Γ Γ C C L2(Γ)
3
β β (cid:107)curl Ξ (cid:107)2 = 0. If (11) is satisfied then Ξ = 0 which also implies Ξ = 0 and one recov-
C I Γ C L2(Γ) C I
ers interface conditions (8). (cid:3)
The problem at hand is then reduced to the following iterative algorithm where the two problems are
solved separately (with appropriate variational formulations and discretization)

∆p(k+1) =0 in ΩI
∂p(k+1) −−→ 1 (14)
 ∂n +βIcurlΓcurlΓp(k+1) = iωµcurlΓ(E(k)+βIcurlE(k)×n) onΓ

curlcurlE(k+1)−κ2E(k+1) =0 in ΩC
−−→ −−→ ∂p(k) (15)
curlE(k+1)×n+βCcurlΓ(curlE(k+1)·n) =iωµcurlΓ(p(k)+βC ∂n ) onΓ.
It order to ensure well posed problems for p(k+1) and E(k+1), a variational study of problems (14) and
(15) show that sufficient conditions are respectively Reβ ≤0 and Reβ ≥0, Imβ ≤0.
I C C
2. Iteration operator in the case of concentric spheres
Thissectionisdedicatedtothecomputationoftheeigenvaluesoftheiterationoperator,denotedbyT,
inthecaseofconcentricspheres.Themaingoalistostudythedependenceoftheseeigenvaluesonβ and
I
β and show that non zero values of β and β improve the convergence of the iterative procedure. In
C I C
the case f =0 (meaning the solution is 0) better behavior correspond with eigenvalues of modulus closer
to 0 in the case of a successive iterative procedure. Consider the sphere Ω=B ⊂R3 with radius R>1.
R
The insulating and conducting regions are respectively given by B \B and B , where B is the unit
R 1 1 1
sphere (the case of a sphere of radius r <R can be easily deduced using an appropriate scaling). Assume
that g (resp. g ) represents the right side on Γ in (15) (resp. (14)), and let us define g = (g ,g )T.
C I C I
Performing one iteration consists in computing
  
0 TC g(Cn)
g(n+1)=Tg(n)=   (16)
TI 0 g(n)
I
with
(cid:16) (cid:17) 1 (cid:16) (cid:17) (cid:16) (cid:17) −−→ (cid:16) (cid:17)
T g(n) := curl E(n+1)+β curlE(n+1)×n , T g(n) :=iωµcurl p(n+1)+β ∂ p(n+1)
I C iωµ Γ I C I Γ C n
(17)
1
Let Ynm, n = 0,1,..., −n ≤ m ≤ n denote the spherical harmonics and set Umn := (cid:112) ∇S2Ynm
n(n+1)
and Vm(xˆ):=xˆ×Um(xˆ), xˆ ∈S2. We remark from the expression of g that this field belongs to span
n n C
{Vm,n = 1,2...,−n ≤ m ≤ n}. We also denote by j the spherical Bessel function of the first kind of
n n
order n.
Proposition 2.1 Assume that
∞ n ∞ n
(cid:88) (cid:88) (cid:88) (cid:88)
g = gn,mVm and g = gn,mYm(xˆ).
C C n I I n
n=1m=−n n=1m=−n
Then we have
∞ n ∞ n
(cid:88) (cid:88) (cid:88) (cid:88)
T (g )= (T )mgn,mVm and T (g )= (T )mgn,mYm(xˆ)
C I C n I n I C I n C n
n=1m=−n n=1m=−n
4
with
(cid:32) (cid:33)
1
(T )m :=−iωµ B +β n(n+1)A /(A +β B ) (18)
C n (cid:112) I C I I I I
n(n+1)
(cid:32) (cid:33)
(T )m := 1 −1 B +β (cid:112)n(n+1)A (cid:46)(A +β B ) (19)
I n iωµ (cid:112)n(n+1) C I C C C C
and A := −n(cid:0)1+R−2n(cid:1) , B := n(n+1)(cid:0)1−R−2n(cid:1), A := −(cid:112)n(n+1)(j (κ)+κj(cid:48)(κ)) and B :=
I I C n n C
−(n(n+1))23jn(κ).
Proof.Thesolutionpofproblem(14)canbeexpressedasp(x)=(cid:80)∞ (cid:80)n (cm|x|n+dm|x|−n)Ym(xˆ).
n=1 m=−n n n n
From the homogeneous exterior boundary conditions we get cm =−dmR−2n. Using the boundary condi-
n n
tion (9) we obtain
(A +β B )dm =gn,m.
I I I n I
The expression of (T )m then follows from evaluating the expression of T (g ).
C n C I
i
ThesolutionEofproblem(15)canbeexpressedasE(x)=(cid:80)∞ (cid:80)n amMm(x)− bmcurlMm(x)
n=1 m=−n n n κ n n
with Mm(x):=curl(cid:0)xj (κ|x|)Ym(xˆ)(cid:1). Since
n n n
(cid:18) (cid:19)
Mm(x)=−(cid:112)n(n+1)j (κ|x|)Vm(xˆ) and curlMm(x)=−(cid:112)n(n+1) 1 j (κ|x|)+κj(cid:48)(κ|x|) Um(xˆ)
n n n n |x| n n n
we get at |x|=1,
∞ n
curlE×n=(cid:88) (cid:88) −am(cid:112)n(n+1)(cid:16)j (κ)+κj(cid:48)(κ)(cid:17)Vm(xˆ)+iκbm(cid:112)n(n+1)j (κ)Um(xˆ).
n n n n n n n
n=1m=−n
In addition, because curl Um =0, we obtain at |x|=1,
Γ n
∞ n
−c−u→rlΓcurlΓE= (cid:88) (cid:88) −amn (n(n+1))23 jn(κ)Vmn(xˆ). (20)
n=1m=−n
Therefore, combination of (20),(20) and the boundary condition (10) leads to
∞ n
−−→ (cid:88) (cid:88)
curlE×n+β curl curl E= (A +β B )Vm(xˆ). (21)
C Γ Γ C C C n
n=1m=−n
This implies in particular that
am =gn,m/(A +β B ), bm =0. (22)
n C C C C n
The expression (T)m then directly follows from evaluating T (g ). (cid:3)
n I C
Using the structure of T one observes that the operator T2 is diagonal with eigenvalues
(T)m :=(T )m·(T )m. (23)
n I n C n
Proposition 2.2 The leading term in the asymptotic expansion of the amplification coefficient for large
n and any β ,β is given by
I C
(cid:12) (cid:12)
|(T)m|∼(cid:12)(cid:12)1−nβC(cid:12)(cid:12) as n→∞.
n (cid:12)1+nβ (cid:12)
C
5
Proof. The asymptotic expansion for the Bessel functions of the first kind with complex argument [1,
Formula 9.3.1] for a fixed complex argument κ and a large integer n is given by jn(κ) ∼ √21πn(cid:0)2enκ(cid:1)n
and j(cid:48)(κ) ∼ nj (κ).We then have the following asymptotic expansions of the coefficients appearing in
n κ n
Proposition 2.1 A ∼ −n,B ∼ n2, A ∼ −n2j (κ) and B ∼ −n3j (κ). Plugging these asymptotics in
I I C n C n
formula (18) and (19) respectively give for large n
1
(T )m ∼iωµ(1−nβ )/(1−nβ ) and (T )m ∼ (1−β n)/(−1−β n).
C n C I I n iωµ I C
This gives the announced result. (cid:3)
We then conclude that as long as Reβ >0 the coefficient (T)m has a modulus strictly smaller than 1
C n
for large n. It is surprising that the asymptotic behavior is independent from the values of β which may
I
suggestthatthisparameterhaslessimportantinfluenceinacceleratingtheiterations.Moreimportantly,
the asymptotic formula shows that |(T)m(Reβ >0)|<|(T)m(Reβ =0)| for n sufficiently large which
n C n C
would indicate for the cases Reβ >0 better convergence properties than for the natural choice β =0
C I
andβ =0.Wealsoremarkthatpurelyimaginaryvaluesofβ donotimprovetheconvergenceproperties
C C
for large modes.
We plot in Figure 1 the amplification coefficient (23) of the iterative procedure for each mode n,m.
We observe that the above conclusions also hold for all modes and that better behavior is observed for
Reβ >0. The asymptotic behavior of |(T)m| is also confirmed by Figure 1 right.
C n
Figure1.Plotsofn(cid:55)→|(T)mn|forthechoiceof(βI,βC)={(-1.e-2,1.e-2),(-1.e-2,1.e-1)(left)anddifferentchoicesofβI =−βC
(right)}.
3. 3D Finite Elements preliminary experiments
Preliminary results of the new formulation are presented in this section. The insulating region is given
by B \B where R = 2 and r = 1, and the conduction region by B . The electromagnetic coefficients
R r r
are chosen to be; σ = 1 for the conductivity and µ = 1 for the permeability. The frequency is set to
ω =π/4. The coefficients β and β similar to the ones used in Figure 1 (e.g. β ,β )={(−1.e−2,1.e−
I C I C
2),(−1.e−2,1.e−1)), and the source term is given by f = sin(x+iy). For problem (15) (resp. (14)),
the N´ed´elec (resp. P -Lagrange) finite elements are used. Figure 2 exhibits convergence of the residual
1
for the iterative procedures. Clearly the case (β ,β )=(0,0) is less performant.
I C
Moreadvancednumericalanalysisoftheproposedschemewillbethesubjectofforthcomingwork.We
shall in particular numerically and theoretically discuss optimal choices for the impedance parameters
6
(a) (b) (c)
Figure2.Computedsolutionforthepotentialpintheinsulator(a)andforelectricfieldintheconductor(b).Theresidual
oftheDDMiterations(c).
β and β and explore the possibility of using parameters that are operators of an appropriate negative
C I
orderthatwouldprovideanasymptoticlimitof|(T)m|strictlysmallerthan1(thisisclearlynotthecase
n
of constant parameters as indicated by Proposition 2.2).
Acknowledgements
Y. Boubendir gratefully acknowledges support from NSF through grant No. DMS-1319720.
References
[1] MiltonAbramowitzandIreneAStegun.Handbookofmathematicalfunctions:withformulas,graphs,andmathematical
tables. Number55.CourierCorporation,1964.
[2] Ana Alonso Rodr´ıguez, Enrico Bertolazzi, Riccardo Ghiloni, and Alberto Valli. Finite element simulation of eddy
currentproblemsusingmagneticscalarpotentials. J. Comput. Phys.,294:503–523,2015.
[3] Ana Alonso Rodr´ıguez and Alberto Valli. Eddy current approximation of Maxwell equations, volume 4 of MS&A.
Modeling, Simulation and Applications. Springer-VerlagItalia,Milan,2010. Theory,algorithmsandapplications.
[4] AnaAlonsoRodr´ıguezandAlbertoValli. Finiteelementpotentials. Appl. Numer. Math.,95:2–14,2015.
[5] YassineBoubendir.Ananalysisofthebem-femnon-overlappingdomaindecompositionmethodforascatteringproblem.
Journal of computational and applied mathematics,204(2):282–291,2007.
[6] YassineBoubendir,XavierAntoine,andChristopheGeuzaine. Aquasi-optimalnon-overlappingdomaindecomposition
algorithmforthehelmholtzequation. Journal of Computational Physics,231(2):262–280,2012.
[7] MartinJGander,Fr´ed´ericMagoules,andFr´ed´ericNataf.Optimizedschwarzmethodswithoutoverlapforthehelmholtz
equation. SIAM Journal on Scientific Computing,24(1):38–60,2002.
[8] AlfioQuarteroniandAlbertoValli. Domaindecompositionmethodsforpartialdifferentialequations. NumberCMCS-
BOOK-2009-019.OxfordUniversityPress,1999.
[9] AnaAlonsoRodrA˜guez,RalfHiptmair,andAlbertoValli. Ahybridformulationofeddycurrentproblems. Numerical
Methods for Partial Differential Equations,21(4):742–763,2005.
[10]AndreaToselliandOlofWidlund. Domain decomposition methods: algorithms and theory,volume3. Springer,2005.
7