Análisis de un método BEM–FEM para la resolución numérica de un problema de magnetostática en R3

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) An´alisis de un m´etodo BEM–FEM para la resoluci´on num´erica de un problema de magnetost´atica en R3 P. Salgado 1and V. Selgas2 1Dpto. de Matem´atica Aplicada, Universidad de Santiago de Compostela. E-mail: [email protected]. 2Dpto. de Matem´aticas, Universidad de La Coru˜na. E-mail: [email protected]. Palabras clave: magnetost´atica, elementos finitos, elementos de contorno, potencial escalar reducido, potencial total. Resumen En este trabajo analizamos una formulaci´on BEM–FEM sim´etrica para resolver el problema tridimensional de magnetost´atica utilizando potenciales escalares magn´eticos. 1. Introducci´on El problema de magnetost´atica consiste en calcular el campo magn´etico creado por una corriente dada que no depende del tiempo. Concretamente, suponiendo conocida la densidad de corriente J, el problema consiste en encontrar el campo magn´etico Hdefinido en R3y cumpliendo rot H =J,(1) div (µH)=0,(2) H(x) = O(|x|−1) cuando |x| → ∞,(3) donde µes la permeabilidad magn´etica. Para resolver este modelo, en ingenier´ıa el´ectrica se utilizan distintos m´etodos num´ericos, cuya diferencia fundamental son las inc´ognitas principales del problema; v´ease por ejemplo [5]. Los resultados num´ericos de la literatura indican que las formulaciones en t´erminos de potenciales escalares son las m´as eficientes desde un punto de vista computacional. En particular, la combinaci´on de los denominados potencial escalar reducido ypotencial total parece ser la m´as eficiente. Esta estrategia, que combina el potencial total en los materiales magn´eticos sin corriente y el potential reducido en el aire y en los materiales no magn´eticos que transportan corriente, fue introducida en 1979 por Simkin y Trowbridge [9] para dominios bidimensionales y extendida 1 P. Salgado and V. Selgas posteriormente a dominios tridimensionales. Sin embargo, el an´alisis matem´atico de esta formulaci´on y de su resoluci´on num´erica con un m´etodo de elementos finitos no se ha realizado hasta fechas muy recientes en [2]. Es importante se˜nalar que el an´alisis desarrollado en [2] es para dominios tridimensionales acotados y por tanto requiere a˜nadir al problema (1–2) condiciones de contorno aproximadas, dado que el dominio natural del problema es todo el espacio. Ahora bien, las ecuaciones son homog´eneas con coeficientes constantes en el exterior de una regi´on acotada, propiedad que permite combinar un m´etodo de elementos finitos (FEM) con un m´etodo de elementos de contorno (BEM). As´ı, en este trabajo analizaremos una formulaci´on en t´erminos de los potenciales total y reducido planteada en R3y propondremos un m´etodo BEM–FEM para su resoluci´on num´erica. Cabe se˜nalar que estudiaremos el problema considerando un dominio magn´etico que puede ser m´ultiplemente conexo, lo cual conduce a trabajar con un potencial total multivaluado. Seguiremos las ideas de [2] para el tratamiento del potencial multivaluado, evitando de este modo, las costosas aproximaciones num´ericas de las funciones de base del espacio de campos arm´onicos de Neumann, que se construyen por ejemplo en [6] para un problema cuasi–estacionario. Describiremos la discretizaci´on del problema y mostraremos la convergencia del m´etodo propuesto. 2. El problema en t´erminos de potenciales escalares Sea Ω la regi´on ocupada por los materiales magn´eticos; i.e. Ω :={x∈R3;µ(x)6=µ0 }, donde µ0denota la permeabilidad magn´etica del vac´ıo. El dominio Ω se supone abierto, acotado y conexo, con frontera Γ := ∂Ω conexa y Lipschitz continua, pero no necesariamente simplemente conexo; v´ease la Figura 1. Supongamos que la densidad de corriente J∈L2(R3) cumple J|Ω=0y divJ= 0 en R3. Entonces el campo vectorial Tdefinido a partir de la Ley de Biot–Savart, T(x) := 1 4πZR3 J(y)×(x−y) |x−y|3dy∀x∈R3, pertenece al espacio H1(R3) y satisface rot T =Jen R3,(4) div T= 0 en R3,(5) kTk1,R3≤CkJk0,R3,(6) siendo C > 0 una constante independiente de J; cf. [7, Lemma 2.2(b)]. En particular, rot H =rot T en R3, luego existe un ´unico potencial ϕR∈W1(R3) := {ψdistribuci´on en R3;ψ √1+|x|2∈L2(R3),∇ψ∈L2(R3)}, tal que H=T−∇ϕRen R3.(7) El escalar ϕRse conoce como potencial reducido y permite representar Hseg´un (7) en todo el espacio R3. Sin embargo, para evitar errores de aproximaci´on en los materiales magn´eticos (consecuencia de que Ty∇ϕRson grandes y de magnitud similar en ellos) introduciremos otro potencial escalar definido en Ω. 2 Resoluci´on num´erica de un problema de magnetost´atica En general, Ω no es simplemente conexo, pero suponemos que existe un n´umero finito de superficies Σj(j= 1, . . . , J) tales que Σj⊂Ω, ∂Σj⊂∂Ω, Σj∩Σk=∅(para j6=k) y de modo que el conjunto e Ω := Ω \(∪J j=1Σj) es simplemente conexo y pseudo–Lipschitz. Tambi´en suponemos que para cada j= 1, . . . , J existe una curva γj⊂Γ que corta a Σj una ´unica vez y que es la frontera de una superficie abierta Sj⊂R3\Ω; cf. [1] y v´ease la Figura 1. SoportedeJ Figura 1: Esquema del dominio. Para cada j= 1, . . . , J, fijamos un vector njunitario normal a Σj. Denotamos las dos caras de Σjpor Σ+ jy Σ− j, de modo que njtenga la orientaci´on normal exterior a e Ω sobre Σ+ j. Dada una funci´on e ψ∈H1(e Ω), sea e ∇e ψ∈L2(Ω) la extensi´on de ∇e ψ∈L2(e Ω) a todo Ω, y [ e ψ]Σj:= e ψ|Σ− j−e ψ|Σ+ jel salto de e ψa trav´es de Σja lo largo de la direcci´on de nj. Introducimos el espacio V:= {e ψ∈H1(e Ω) ; [ e ψ]Σj∈R∀j= 1, . . . , J}, que permite caracterizar las funciones con rotacional nulo; concretamente, si q∈L2(Ω) con rotq =0en Ω, entonces existe e ψ∈Vtal que q=e ∇e ψ; v´ease por ejemplo [1]. En particular, dado que rotH =0en Ω, existe un ´unico eϕ∈V/Rtal que H=−e ∇eϕen Ω,(8) yeϕse conoce como potencial total. N´otese que, aplicando el Teorema de Stokes se deduce que [eϕ]Σj=RSjJ·nSjdS =: Ij (j= 1, . . . , J), que son datos del problema. En consecuencia, podemos reescribir eϕcomo eϕ=ϕ+ J X j=1 Ije φj,(9) donde e φjes el ´unico elemento de H1(Ω \Σj) que cumple [e φj]Σj= 1 ; ZΩ\Σj µe ∇e φj·∇ψ= 0 ∀ψ∈H1(Ω), yϕ:= Peϕ∈H1(Ω), con P:V→H1(Ω) el operador de proyecci´on asociado a la suma directa V=H1(Ω) ⊕h{e φj}J j=1i. 3 P. Salgado and V. Selgas An´alogamente, como el campo T∈L2(Ω) tambi´en tiene rotacional nulo en Ω, existe un ´unico potencial eϕT∈V/Rtal que T=−e ∇eϕTen Ω.(10) En virtud de (7), (8) y (10), sabemos que −∇ϕR|Ω=H|Ω−T|Ω=−e ∇(eϕ−eϕT), luego los potenciales eϕ−eϕTyϕR|Ωcoinciden salvo constante aditiva. En particular, eϕ−eϕT∈H1(Ω) y, aplicando la representaci´on (9), deducimos eϕ−Peϕ=eϕT−PeϕT= J X j=1 Ije φjen Ω.(11) A continuaci´on, reescribimos las ecuaciones (1–2) como un problema de transmisi´on usando la siguiente representaci´on del campo magn´etico: H=(−e ∇eϕen Ω, T−∇ϕRen R3\Ω.(12) Para esta representaci´on, la continuidad de las trazas tangencial y normal de HyµHa trav´es de Γ se traduce en las siguientes condiciones de transmisi´on: −e ∇eϕ×n=T×n−∇ϕR×nsobre Γ,(13) −µe ∇eϕ·n=µ0T·n−µ0∇ϕR·nsobre Γ.(14) Con el fin de estudiar la condici´on (13), para cada potencial regular ψ∈H1(Ω) introducimos rotΓψ:= ∇ψ×n, que se denomina rotacional vectorial de ψsobre la superficie Γ. El siguiente resultado puede consultarse en [3, Corollary 3.7]. Lema 2.1 El operador diferencial rotΓ:H1(Ω) →H−1/2(Γ), est´a bien definido y es lineal continuo. Adem´as, su n´ucleo es el conjunto de las funciones de H1(Ω) con traza constante sobre Γ. En virtud de (11) y del Lema 2.1, la condici´on de transmisi´on (13) es equivalente a P(eϕ−eϕT)−ϕR∈Rsobre Γ.(15) Teniendo en cuenta los resultados anteriores, proponemos la siguiente formulaci´on del problema de magnetost´atica. Problema 2.2 Hallar los potenciales eϕ∈V/RyϕR∈W1(R3\Ω), con [eϕ]Σj=Ij (j= 1, . . . , J), que cumplan las siguientes ecuaciones en sentido d´ebil: −div(µe ∇eϕ)=0 en Ω, −4ϕR= 0 en R3\Ω, Peϕ+R=PeϕT+ϕR+Rsobre Γ, −µ∂eϕ ∂n=µ0T·n−µ0 ∂ϕR ∂nsobre Γ. 4 Resoluci´on num´erica de un problema de magnetost´atica 3. Una formulaci´on variacional BEM–FEM sim´etrica En primer lugar, estudiamos el Problema 2.2 en la regi´on no acotada R3\Ω. En esta regi´on, el potencial ϕRes arm´onico y, por tanto, admite la representaci´on integral ϕR(x) = ZΓµ∂Φ(x−y) ∂n(y)ϕR(y)−Φ(x−y)∂ϕR ∂n(y)¶dS(y)∀x∈R3\Ω,(16) donde Φ(x) := −1 4π|x|(x∈R3\{0}) denota la soluci´on fundamental del operador laplaciano en R3. A partir de la representaci´on (16) y de las relaciones de salto de los operadores de simple y doble capa, se deducen las ecuaciones integrales (cf. [4], [8, Theorem 3.1.2]) (1 2I −K)ϕR+V∂ϕR ∂n= 0 y (1 2I+K∗)∂ϕR ∂n+NϕR= 0 sobre Γ,(17) donde VyKdenotan los operadores de simple y doble capa: Vη(x) := ZΓ Φ(x−y)η(y)dS(y) y Kη(x) := ZΓ ∂Φ(x−y) ∂n(y)η(y)dS(y). Adem´as, en (17), K∗eIrepresentan el operador adjunto de Ky la identidad, respectivamente, y Nes el operador hipersingular: Nη(x) := −∂ ∂n(x)µZΓ ∂Φ(x−y) ∂n(y)η(y)dS(y)¶. A continuaci´on se˜nalamos las propiedades fundamentales de VyK; cf. [6, Lemmas 3.3–3.4]. Lema 3.1 Dado s∈[−1/2,1/2], los operadores integrales V:Hs−1/2(Γ) →Hs+1/2(Γ) yK:Hs+1/2(Γ) →Hs+1/2(Γ), son lineales y continuos. Adem´as, el operador Ves H−1/2(Γ)–el´ıptico. Con el fin de estudiar las propiedades del operador N, introducimos el espacio L2 τ(Γ) := {q∈L2(Γ) ; q·n= 0 sobre Γ}, y la traza tangencial πτ:H1/2(Γ) →L2 τ(Γ) dada por πτq:= n×(q|Γ×n)∀q∈ C∞(Ω)3. Denotamos su imagen por Vπ:= πτ(H1/2(Γ)), que es densa en L2 τ(Γ) y en consecuencia podemos definir V0 πel espacio dual de Vπcon espacio pivote L2 τ(Γ). Consideramos πτ ∗: V0 π→H−1/2(Γ) el operador adjunto de πτ, y sea g rotΓ:V→V0 πel operador diferencial g rotΓe ψ:= e ∇e ψ×n∀e ψ∈V. Puede comprobarse que la imagen de una funci´on ψ∈H1(Ω) a trav´es del operador g rotΓdepende ´unicamente de la traza ψ|Γ∈H1/2(Γ). Por lo tanto, podemos introducir el operador diferencial de superficie rotΓ:H1/2(Γ) →V0 πde modo que rotΓ(ψ|Γ) = g rotΓψ:= ∇ψ×n∀ψ∈H1(Ω). Representamos por V:H−1/2(Γ) →H1/2(Γ) el operador de simple capa Vactuando por componentes y definimos b V:= V◦πτ ∗:V0 π→H1/2(Γ). El siguiente resultado queda probado en [6, Lemma 3.3]. 5 P. Salgado and V. Selgas Lema 3.2 El operador N:H1/2(Γ) →H−1/2 0(Γ) es lineal y continuo. Adem´as, para cualesquiera ψ1, ψ2∈H1/2(Γ), se cumple la identidad hNψ1, ψ2i1/2,Γ=hrotΓψ2,πτb VrotΓψ1iV0 π×Vπ. En particular, el operador Ncancela las constantes y es H1/2(Γ)/R–el´ıptico. Introducimos la inc´ognita auxiliar λ:= ∂ϕR ∂n∈H−1/2 0(Γ). Entonces, teniendo en cuenta la condici´on de transmisi´on (15) junto con la propiedad K1 = 1/2 y el Lema 3.2, reescribimos las ecuaciones integrales (17) como hη, (1 2I −K)Peϕi1/2,Γ+hη, Vλi1/2,Γ=hη, (1 2I −K)PeϕTi1/2,Γ∀η∈H−1/2(Γ),(18) hλ, ψi1/2,Γ=hλ, (1 2I −K)ψi1/2,Γ−hNPeϕ, ψi1/2,Γ+hNPeϕT, ψi1/2,Γ∀ψ∈H1/2(Γ).(19) En segundo lugar, hacemos una formulaci´on variacional est´andar del Problema 2.2 en el dominio acotado Ω, aplicando las propiedades (14) y (19): (µe ∇eϕ, ∇ψ)0,Ω−µ0hλ, (1 2I −K)ψi1/2,Γ+µ0hNPeϕ, ψi1/2,Γ= =µ0hNPeϕT, ψi1/2,Γ−µ0hT·n, ψi1/2,Γ∀ψ∈H1(Ω). (20) Consideramos las formas bilineales a1:V×V→R,a2:H−1/2 0(Γ) ×H−1/2 0(Γ) →R yb:H−1/2 0(Γ) ×V→R, definidas como a1(e ψ1,e ψ2) := (µe ∇e ψ1,e ∇e ψ2)0,Ω+µ0hNP e ψ1,Pe ψ2i1/2,Γ, a2(η1, η2) := µ0hη2,Vη1i1/2,Γ;b(η, e ψ) := µ0hη, (1 2I −K)Pe ψi1/2,Γ. Introducimos adem´as las aplicaciones lineales l1:V→Ryl2:H−1/2 0(Γ) →R, dadas por l1(e ψ) := µ0hNPeϕT,Pe ψi1/2,Γ−µ0hT·n,Pe ψi1/2,Γ;l2(η) := b(η, eϕT). Con esta notaci´on y teniendo en cuenta las ecuaciones (11), (18) y (20), proponemos la siguiente formulaci´on variacional BEM–FEM sim´etrica del Problema 2.2. Problema 3.3 Hallar eϕ∈V/Ryλ∈H−1/2 0(Γ), con [eϕ]Σj=Ij(j= 1, . . . , J), tales que a1(eϕ, e ψ)−b(λ, e ψ) = l1(e ψ), a2(λ, η) + b(η, eϕ) = l2(η), para cualesquiera e ψ∈V/Ryη∈H−1/2 0(Γ), con [e ψ]Σj= 0 (j= 1, . . . , J). 6 Resoluci´on num´erica de un problema de magnetost´atica Lema 3.4 El Problema 3.3 tiene una ´unica soluci´on. Adem´as, si eϕ∈V/R,λ∈H−1/2 0(Γ) es la soluci´on del Problema 3.3 y definimos el potencial ϕR ∗(x) := ZΓµ∂Φ(x−y) ∂n(y)(Peϕ(y)−PeϕT(y)) −Φ(x−y)λ(y)¶dS(y)∀x∈R3\Ω, entonces el siguiente campo satisface las ecuaciones (1–2): H:= (−e ∇eϕen Ω, T−∇ϕR ∗en R3\Ω.(21) N´otese que a partir del campo Hse obtiene la inducci´on magn´etica B=µHen R3. 4. Resoluci´on num´erica mediante un m´etodo BEM–FEM En todo lo que sigue supondremos que Ω es un poliedro de Lipschitz. Consideramos una familia {Th}h>0de mallas tetra´edricas del dominio Ω; asumimos tambi´en que cada superficie de corte Σjes poli´edrica y formada por la uni´on de caras de tetraedros, por lo que Thes tambi´en una malla de e Ω. Adem´as, para cada h > 0 denotamos por TΓ hla malla triangular inducida por Thsobre la superficie Γ. Entonces aproximamos las funciones de los espacios VyH−1/2 0(Γ) utilizando los subespacios de dimensi´on finita Vh:= {e ψ∈ C0(e Ω) ; e ψ|K∈P1(K)∀K∈ Th,[e ψ]Σj∈R∀j= 1, . . . , J}, Λh:= {η∈L2(Γ) ; η|T∈R∀T∈ TΓ h,RΓη dS = 0}. Proponemos as´ı la siguiente versi´on discreta del Problema 3.3. Problema 4.1 Hallar eϕh∈Vh/Ryλh∈Λh, con [eϕh]Σj=Ij(j= 1, . . . , J), tales que a1(eϕh,e ψ)−b(λh,e ψ) = l1(e ψ), a2(λh, η) + b(η, eϕh) = l2(η), para cualesquiera e ψ∈Vh/Ryη∈Λh, con [e ψ]Σj= 0 (j= 1, . . . , J). N´otese que la definici´on de binvolucra el c´alculo expl´ıcito de los potenciales e φj. Para evitar esta complicaci´on pr´actica, aprovecharemos que las funciones test presentan la regularidad adicional Λh⊂L2 0(Γ) ⊂H−1/2 0(Γ). Concretamente, esta regularidad da sentido al producto (η, (1 2I − K)e ψ)0,Γpara cualesquiera η∈Λhye ψ∈V; v´ease el Lema 3.1. En consecuencia, podemos introducir la forma bilineal ˆ b:L2 0(Γ) ×V→Ry el operador ˆ l2:L2 0(Γ) →Rcomo las siguientes versiones de byl2que no involucran la proyecci´on P: ˆ b(η, e ψ) := µ0(η, (1 2I −K)e ψ)0,Γ;ˆ l2(η) := ˆ b(η, eϕT). As´ımismo, definimos la forma bilineal ˆa1:V×V→Ry el operador ˆ l1:V→Rcomo ˆa1(e ψ1,e ψ2) := (µe ∇e ψ1,e ∇e ψ2)0,Ω+µ0hg rotΓe ψ2,πτb Vg rotΓe ψ1iV0 π×Vπ, ˆ l1(e ψ) := µ0hg rotΓe ψ, πτb Vg rotΓeϕTiV0 π×Vπ−µ0hT·n,e ψi0,Γ. 7 P. Salgado and V. Selgas Finalmente, seguimos la estrategia propuesta en [2] para aproximar las trazas tangencial y normal del campo T. De este modo, consideramos TND hyTRT h, las interpoladas de Ten los espacios de elementos finitos de N´edelec y de Raviart–Thomas de primer orden, respectivamente; entonces tenemos las aproximaciones T×n∼ =TND h×nyT·n∼ =TRT h·n. Adem´as, si eϕT hdenota la interpolada de eϕTen Vh, entonces TND h=−e ∇eϕT hy, en particular, TND h×n=−g rotΓeϕT hsobre Γ, propiedad que nos permite determinar eϕT hreescribiendo el campo TND hen t´erminos de la base dada en [6, Proposition 4.2]. Por lo tanto, aproximamos ˆ l1(ψ) y ˆ l2(η) por ˆ l(h) 1(ψ) := −µ0hTRT h·n, ψi1/2,Γ+µ0hrotΓ ψ, πτb V(TND h×n)iV0 π ×Vπ;ˆ l(h) 2(η) := ˆ b(η, eϕT h). Proponemos el siguiente esquema para resolver num´ericamente el Problema 3.3. Problema 4.2 Hallar eϕh∈Vh/Ryλh∈Λh, con [eϕh]Σj=Ij(j= 1, . . . , J), tales que ˆa1(eϕh,e ψ)−ˆ b(λh,e ψ) = ˆ l(h) 1(e ψ), a2(λh, η) + ˆ b(η, eϕh) = ˆ l(h) 2(η), para cualesquiera e ψ∈Vh/Ryη∈Λh, con [e ψ]Σj= 0 (j= 1, . . . , J). Proposici´on 4.3 El Problema 4.2 tiene una ´unica soluci´on. Adem´as, si la soluci´on (eϕ, λ) del Problema 3.3 es tal que el campo Hdefinido en (21) presenta la regularidad H|Ω∈ Hr(Ω) yH|R3\Ω∈Wr(R3\Ω)3, con r∈(0,1], entonces keϕ−eϕhk1,Ω+kλ−λhk−1/2,Γ≤C hr³kHkr,Ω+kHkWr(R3\Ω)3+kJk0,R3´. Agradecimientos Los autores agradecen las valiosas sugerencias de los profesores A. Berm´udez, S. Meddahi y R. Rodr´ıguez. Referencias [1] C. Amrouche, C. Bernardi, M. Dauge, and V. Girault, Vector potentials in three-dimensional nonsmooth domains, Math. Meth. Appl. Sci., 21 (1998), pp. 823–864. [2] A. Berm´udez, R. Rodr´ıguez, P. Salgado, A finite element method for the magnetostatic problem in terms of scalar potentials, Preprint DIM 2006-28, Universidad de Concepci´on, Concepci´on, 2006. [3] A. Buffa, M. Costabel, D. Sheen, On traces for H(curl,Ω) in Lipschitz domains, J. Math. Anal. Appl., 276 (2002), pp. 845–867. [4] M. Costabel, Symmetric methods for the coupling of finite elements and boundary elements, The Mathematics of Finite Elements and Applications IV, Academic Press, London 1988. [5] Ch. Magele, H. St¨ogner K. Preis, Comparison of different finite element formulations for 3D magnetostatic problems, IEEE Transaction on Magnetics, 24 (1) (1988), pp. 31–34. [6] S. Meddahi, V. Selgas, A Mixed-FEM and BEM coupling for a three–dimensional eddy current problem, M2AN Math. Model. Numer. Anal., 37 (2) (2003), pp. 291–318. [7] A. J. Meir, P. G. Schmidt, Variational Methods for stationary MHD flow under natural interface conditions, Nonlinear Anal., 26 (1996), pp. 659–689. [8] J.-C. N´ed´elec, Acoustic and electromagnetic equations: Integral representations for harmonic problems, Springer–Verlag, New York, 2001. [9] J. Simkin, C.W. Trowbridge, On the use of the total scalar potential in the numerical solution of field problems in electromagnetics, Int. J. Meth. Eng., 14 (1979), pp. 423–440. 8