UNIVERSIDAD DE OVIEDO Departamento de Matem´aticas

´ ´ ´ ANALISIS MATEMATICO Y NUMERICO DE LAS ECUACIONES DE MAXWELL ´ CUASIESTATICAS

Tesis Doctoral Virginia Selgas Buznego

´ D. Benjam´ın Dugnol Alvarez, director del Departamento de Matem´aticas de la Universidad de Oviedo, autoriza la presentaci´ on de la Tesis Doctoral titulada “An´alisis matem´atico y num´erico de las ecuaciones de Maxwell cuasiest´aticas”, que ha sido realizada por D˜ na. Virginia Selgas Buznego bajo la direcci´on de D. Salim Meddahi.

Oviedo, febrero de 2006

´ Fdo. D. Benjam´ın Dugnol Alvarez

D. Salim Meddahi, Profesor Titular de Matem´atica Aplicada del Departamento de Matem´aticas de la Universidad de Oviedo, CERTIFICA: que la presente memoria, titulada “An´alisis matem´atico y num´erico de las ecuaciones de Maxwell cuasiest´aticas”, ha sido realizada por D˜ na. Virginia Selgas Buznego bajo su direcci´on, dando su conformidad para que sea presentada para optar al Grado de Doctor por la Universidad de Oviedo.

Oviedo, febrero de 2006

Fdo. D. Salim Meddahi

Agradecimientos Este trabajo ha sido subvencionado por el Ministerio de Educaci´on, Cultura y Deporte a trav´es de una Beca de Formaci´ on de Profesorado Universitario con referencia AP2001–2318. Deseo agradecer todo el apoyo recibido en el Departamento de Matem´aticas de la Universidad de Oviedo. Muy especialmente, quiero expresar mi gratitud al director de esta tesis, el Profesor Salim Meddahi, por sus propuestas as´ı como por todo el tiempo que ha dedicado a ayudarme. Adem´as, debo un sincero reconocimiento al Profesor Javier Vald´es, que me anim´o a comenzar esta tesis y ha confiado en m´ı desde un primer momento. Tambi´en deseo agradecer la hospitalidad que disfrut´e en el Departamento de Matem´atica Aplicada de la Universidad de Santiago de Compostela durante la realizaci´on de los cursos de doctorado. En particular, no puedo olvidar las sugerencias y el est´ımulo de mi codirectora del trabajo de investigaci´ on, la Profesora Pilar Salgado. Asimismo, quiero agradecer la c´alida acogida del Departamento de Matem´aticas de la Universidad de Trento durante los meses que estuve all´ı. Concretamente, debo a los Profesores Alberto Valli y Ana Alonso su respaldo e interesantes ideas. Finalmente, me gustar´ıa dedicar este trabajo a mi familia y a mis amigos, por su apoyo incondicional y por su paciencia y su cari˜ no.

vii

Introducci´ on Hasta el siglo XIX, los efectos el´ectricos y magn´eticos se consideraban fen´omenos f´ısicos independientes, ligados a las cargas el´ectricas y a los imanes, respectivamente. Sin embargo, en 1820 Oersted observ´o que las corrientes el´ectricas pod´ıan influir sobre una aguja imantada. A continuaci´ on, las aportaciones de Amp`ere y el descubrimiento de la inducci´on por Faraday establecieron las bases para una teor´ıa unificada del electromagnetismo. Hacia 1860 Maxwell reuni´o y complet´o los resultados anteriores, sintetizando las teor´ıas el´ectrica y magn´etica en un u ´nico sistema de ecuaciones. En particular, estas ecuaciones predicen que los campos el´ectricos variables generan campos magn´eticos y que, rec´ıprocamente, los campos magn´eticos variables inducen corrientes el´ectricas. Desde un punto de vista pr´actico, este fen´omeno de inducci´ on electromagn´etica resulta fundamental para dise˜ nar adecuadamente un dispositivo el´ectrico. En efecto, las corrientes inducidas pueden ser u ´tiles o bien producir p´erdidas, seg´ un sea la aplicaci´on dada a una m´aquina el´ectrica. As´ı por ejemplo, las corrientes inducidas generan las fuerzas electromotrices necesarias para el funcionamiento de un motor de inducci´on. Por el contrario, reducen la eficiencia de un transformador el´ectrico porque disipan energ´ıa por efecto Joule; cf. [28]. Consecuentemente, en la pr´actica interesa modelar la distribuci´on de las corrientes inducidas y simular num´ericamente su comportamiento. Esto involucra el an´alisis del sistema de ecuaciones de Maxwell, que en su forma completa presenta una gran complejidad. En muchas situaciones, el sistema de ecuaciones de Maxwell se aproxima ix

x

Introducci´on

con modelos simplificados cuya resoluci´on num´erica es menos costosa. En particular, cuando se estudian m´aquinas el´ectricas que trabajan a bajas frecuencias, un razonamiento asint´otico permite despreciar las corrientes de desplazamiento y deducir as´ı el denominado modelo de eddy currents; v´eanse [1, 6, 14]. De hecho, este modelo es muy frecuente en ingenier´ıa el´ectrica y por ello est´a suscitando un gran inter´es en el ´ambito del an´alisis matem´atico y num´erico; cf. [2, 3, 4, 6, 8, 9, 10, 11, 12, 13, 14, 32, 33, 39, 41]. En general, la resoluci´on num´erica del problema de eddy currents no puede hacerse con t´ecnicas est´andar del tipo del m´etodo de elementos finitos (FEM) o del m´etodo de elementos de contorno (BEM). En efecto, por un lado, no puede aplicarse directamente un m´etodo BEM porque en general las ecuaciones no son homog´eneas y sus coeficientes no son constantes. Por otro lado, tampoco puede utilizarse directamente un m´etodo FEM, puesto que el campo electromagn´etico no se reduce u ´nicamente a una regi´on acotada sino que se extiende a todo el espacio. En este sentido, una estrategia habitual consiste en la introducci´on de un dominio computacional acotado que sea suficientemente grande. Sin embargo, a menudo esta t´ecnica no resulta u ´til ya que involucra un dominio computacional tan grande que el coste del m´etodo num´erico resulta excesivo. Incluso cuando no es as´ı, el problema debe completarse con condiciones de contorno sobre la frontera del dominio computacional y la elecci´on adecuada de estas condiciones puede resultar complicada; cf. [9]. Es importante se˜ nalar que, en las ecuaciones del problema de eddy currents, las no homogeneidades y las no linealidades se localizan en una zona acotada, fuera de la cual los coeficientes son constantes. Es precisamente esta propiedad la que permite combinar BEM y FEM para resolverlo num´ericamente. De hecho, el primer paso en un m´etodo BEM–FEM es descomponer el dominio (no acotado) del problema original en dos regiones: un dominio acotado que recoge las no homogeneidades y las no linealidades de las ecuaciones; la correspondiente regi´on exterior no acotada, donde las ecuaciones son homog´eneas y lineales, con coeficientes constantes.

Introducci´on

xi

Entonces el problema es equivalente a un problema de transmisi´ on entre ambas regiones y con condiciones de transmisi´on sobre la frontera com´ un, frontera que se dice de acoplamiento. A continuaci´ on se reformula el problema utilizando una representaci´ on integral de la soluci´on en la zona exterior. As´ı se deduce un problema planteado en el dominio acotado y con condiciones de contorno globales (i.e. no diferenciales) sobre la frontera de acoplamiento. Llegados a este punto, el problema puede resolverse num´ericamente aproximando las funciones definidas en el dominio acotado y las asociadas a la frontera de acoplamiento utilizando elementos finitos y elementos de contorno, respectivamente. Con esta estrategia, se obtienen aproximaciones de la soluci´on en la regi´on acotada y de las variables auxiliares definidas sobre la frontera. Adem´as, utilizando estas aproximaciones junto con dicha f´ormula de representaci´ on integral, tambi´en podemos estimar la soluci´on en la regi´on no acotada. Para ver diversas aplicaciones de este m´etodo, pueden consultarse [13, 25, 32, 34, 36, 37, 38, 39, 40, 42, 43, 49, 50, 51, 52]. Tal y como hemos se˜ nalado, las condiciones de contorno sobre la frontera de acoplamiento se obtienen a partir de una f´ormula de representaci´ on integral en la regi´on exterior. Como estas condiciones no son u ´nicas, distintas elecciones particulares dan lugar a diversos tipos de m´etodos BEM–FEM, siendo los m´as habituales el m´etodo de Johnson y N´ed´elec y el m´etodo sim´etrico. En el primero se impone una u ´nica ecuaci´on sobre la frontera, que se obtiene a partir de la f´ormula de representaci´on y de la relaci´on de salto del potencial de doble capa. En este caso, para analizar el error suele necesitarse que el operador asociado al potencial de doble capa sea compacto y, en consecuencia, que la frontera de acoplamiento sea regular; cf. [34, 36, 37]. En el m´etodo sim´etrico se a˜ nade una segunda ecuaci´on sobre la frontera de acoplamiento, que se deduce de la f´ormula de representaci´on utilizando las relaciones de salto de la derivada normal del potencial de simple capa. Entonces, para analizar el error habitualmente basta que la frontera de acoplamiento sea Lipschitz continua; cf. [25]. La estrategia BEM–FEM se ha utilizado para el problema de eddy currents bajo r´egimen arm´onico en [13, 32, 41], entre otros. En [32] se propone y se analiza un m´etodo BEM–FEM que utiliza como

xii

Introducci´on

frontera de acoplamiento la propia del conductor y opta por una formulaci´ on sim´etrica en t´erminos del campo el´ectrico. En [14] se propone un esquema BEM–FEM no sim´etrico y formulado en t´erminos del campo magn´etico. Ahora bien, dicho esquema no est´a escrito en la forma est´andar del m´etodo de Johnson y N´ed´elec porque elimina la derivada normal de la variable asociada a la frontera utilizando el operador de Steklov– Poincar´e. Este esquema se ha implementado en un c´odigo denominado TRIFOU (cf. [10]), aunque no conocemos ning´ un an´alisis de convergencia de este esquema. En todo caso, tal y como se˜ nalamos anteriormente, el an´alisis de los m´etodos BEM–FEM no sim´etricos suele requerir que la frontera de acoplamiento sea regular. En la pr´actica, esta propiedad resulta excesivamente restrictiva sobre la geometr´ıa del conductor y, por ejemplo, no permite estudiar dominios poli´edricos. Teniendo en cuenta estas dificultades, en [41] optamos por una formulaci´ on BEM–FEM sim´etrica en t´erminos del campo magn´etico. En esta memoria seguimos dicho art´ıculo para deducir y estudiar esta formulaci´ on BEM–FEM. Tambi´en proponemos y analizamos una formulaci´ on alternativa para el tratamiento de conductores con topolog´ıas complicadas y la generalizamos al problema de eddy currents en r´egimen de evoluci´on en tiempo. Damos algunos resultados num´ericos obtenidos con ambos esquemas para el problema en r´egimen arm´onico. M´as concretamente, en el Cap´ıtulo 1 detallamos las ecuaciones de Maxwell y el modelo de eddy currents, tanto en r´egimen de evoluci´ on como en el caso arm´onico en tiempo. El objetivo del Cap´ıtulo 2 es introducir las herramientas matem´aticas necesarias para estudiar el problema de eddy currents. Precisamente, se˜ nalamos algunas propiedades fundamentales de espacios de funciones de trazas tangenciales y de operadores diferenciales sobre superficies Lipschitz continuas; cf. [20]. En el Cap´ıtulo 3 estudiamos el problema en r´egimen arm´ onico en tiempo para un conductor simplemente conexo; cf. [41]. Para obtener una formulaci´on variacional BEM–FEM de este problema, seguimos la estrategia introducida por Bossavit en [14]. De este modo, dividimos el espacio R3 en dos regiones: Un dominio acotado (que representa el conduc-

Introducci´on

xiii

tor) y su exterior no acotado. A continuaci´ on, introducimos un potencial escalar magn´etico en la regi´on no acotada. Este potencial es arm´onico y, en consecuencia, nos proporciona condiciones de contorno no locales para el problema en el dominio acotado. Entonces, en el dominio acotado podemos aproximar el campo magn´etico utilizando elementos finitos de N´ed´elec. Como se˜ nalamos anteriormente, optamos por una formulaci´ on BEM–FEM de tipo sim´etrico motivados por el estudio de conductores con fronteras no regulares. De hecho, los resultados dados recientemente en [20] nos permiten probar que el m´etodo num´erico que proponemos converge con orden ´optimo siempre y cuando la frontera del conductor sea Lipschitz continua. En consecuencia, podemos tomar como frontera de acoplamiento directamente la propia del conductor, reduciendo el tama˜ no del dominio computacional y, por tanto, el coste del m´etodo. Tambi´en detallamos c´omo implementar este m´etodo en el ordenador. En efecto, lo reescribimos como un sistema lineal y proponemos f´ormulas de cuadratura para aproximar las integrales singulares involucradas en el ensamblaje de la matriz del sistema. Implementamos este m´etodo en un c´odigo MATLAB y lo utilizamos para estudiar un problema con soluci´on exacta conocida expl´ıcitamente. Comparando la soluci´on exacta con los resultados num´ericos obtenidos, validamos el c´odigo que hemos desarrollado y al mismo tiempo comprobamos las propiedades de convergencia del esquema num´erico. En el Cap´ıtulo 4, estudiamos el problema de eddy currents en r´egimen arm´onico para conductores no simplemente conexos. En esta situaci´on, los m´etodos num´ericos propuestos en [32, 41] necesitan introducir superficies de corte, que pueden resultar excesivamente complicadas y costosas de calcular en la pr´actica; al respecto, v´ease tambi´en [8]. En [5] se muestra que la formulaci´ on en t´erminos del campo magn´etico (para el problema en r´egimen arm´onico y planteado en un dominio acotado) admite una estructura de punto de silla que no requiere la construcci´on de superficies de corte; v´ease [31] para una estrategia similar. B´asicamente, dicha estructura se obtiene introduciendo un multiplicador de Lagrange asociado al rotacional del campo

xiv

Introducci´on

magn´etico en el diel´ectrico que rodea al conductor. Seguimos esta estrategia para el problema planteado en todo el espacio. Concretamente, tomamos un dominio computacional que contiene la regi´on conductora e introducimos un multiplicador de Lagrange para imponer la restricci´on asociada al rotacional del campo magn´etico en el diel´ectrico. A continuaci´ on, introducimos un potencial escalar magn´etico en el exterior del dominio computacional. Nuevamente, este potencial es arm´onico y nos proporciona condiciones de contorno no locales para el problema en el dominio computacional. De esta forma, deducimos una formulaci´on BEM–FEM mixta que est´a bien planteada. A partir de esta formulaci´on, deducimos un esquema de Galerkin utilizando elementos finitos de N´ed´elec y elementos finitos de Raviart–Thomas. Probamos que el correspondiente esquema discreto est´a bien planteado y analizamos su orden de convergencia. Adem´as, detallamos su implementaci´ on y la desarrollamos en un c´odigo MATLAB. Testeamos este c´odigo utilizando un problema con soluci´on exacta conocida expl´ıcitamente; as´ı validamos el c´odigo y comprobamos las propiedades de convergencia del m´etodo num´erico. Es importante se˜ nalar que, incluso cuando la fuente de corriente es de tipo sinusoidal, hay situaciones en las que el campo electromagn´etico no presenta un comportamiento arm´onico en tiempo. Por esta raz´on, en el Cap´ıtulo 5 extendemos nuestro m´etodo num´erico para el problema de eddy currents de evoluci´ on en tiempo, sin restricciones topol´ogicas sobre el conductor. Bajo hip´otesis de regularidad sobre los datos del problema, utilizamos el Teorema de J.L. Lions sobre la existencia y unicidad de soluciones de problemas lineales parab´olicos para mostrar que la formulaci´on BEM–FEM que proponemos est´a bien planteada. Deducimos un esquema semidiscreto en espacio utilizando elementos finitos de N´ed´elec y elementos finitos de Raviart–Thomas. A continuaci´ on, probamos que este esquema semidiscreto est´a bien planteado y analizamos su orden de convergencia. Finalizamos el estudio de este esquema reescribi´endolo en forma matricial.

´Indice general Agradecimientos

VII

Introducci´ on

IX

1. Modelo de eddy currents

1

1.1. Modelo de eddy currents . . . . . . . . . . . . . . . . . . . . . . . .

3

1.2. Modelo de eddy currents en r´egimen arm´onico . . . . . . . . . . . .

5

2. Herramientas matem´ aticas

7

2.1. Espacios de Sobolev cl´asicos . . . . . . . . 2.1.1. ´Indice entero no negativo . . . . . 2.1.2. ´Indice real no negativo . . . . . . . 2.1.3. ´Indice real negativo . . . . . . . .

. . . . . . . . . . . . . .

8

. . . . . . . . . . . . . .

9

. . . . . . . . . . . . . .

9

. . . . . . . . . . . . . . 10

2.2. Espacios de Sobolev sobre superficies y teor´ıa de trazas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1. Espacios de Sobolev sobre superficies . . . . . . . . . . . . . 11 2.2.2. Trazas en espacios de Sobolev . . . . . . . . . . . . . . . . . 12 2.3. Espacios H(div, Ω) y H(rot, Ω) . . . . . . . . . . . . . . . . . . . . 14 2.4. Operadores diferenciales sobre superficies . . . . . . . . . . . . . . 17 2.5. Trazas tangenciales en H(rot, Ω) . . . . . . . . . . . . . . . . . . . 18 2.6. Inversa del operador rotacional de superficie en dominios simplemente conexos . . . . . . . . . . . . . . . . . . . 20 2.7. Potenciales arm´onicos en regiones no acotadas . . . . . . . . . . . . 25 xv

xvi

Introducci´on 2.7.1. Operadores integrales sobre superficies . . . . . . . . . . . . 27

2.8. Espacios b´asicos para problemas de evoluci´ on . . . . . . . . . . . . 34 3. Eddy currents en conductores simplemente conexos

37

3.1. Formulaci´on variacional . . . . . . . . . . . . . . . . . . . . . . . . 37 3.1.1. Un potencial escalar magn´etico . . . . . . . . . . . . . . . . 39 3.1.2. Una formulaci´on variacional BEM–FEM . . . . . . . . . . . 40 3.2. Esquema discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.1. Elementos finitos de N´ed´elec de primer orden . . . . . . . . 44 3.2.2. Descripci´on y an´alisis del problema discreto . . . . . . . . . 47 3.2.3. An´alisis de la convergencia . . . . . . . . . . . . . . . . . . 49 3.3. Forma matricial del esquema discreto . . . . . . . . . . . . . . . . . 50 b h (Ωc ) . . . . . . . . . . . . . . . . . . . 50 3.3.1. Base expl´ıcita de X 3.3.2. Sistema de ecuaciones asociado al problema discreto . . . . 53 3.3.3. Implementaci´on del esquema discreto . . . . . . . . . . . . . 54 3.4. Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4.1. Descripci´on del problema test . . . . . . . . . . . . . . . . . 55 3.4.2. Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . 55 4. Eddy currents en conductores no simplemente conexos

59

4.1. Espacios funcionales auxiliares . . . . . . . . . . . . . . . . . . . . 60 4.2. Formulaci´on variacional . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2.1. Una formulaci´on variacional global . . . . . . . . . . . . . . 66 4.2.2. Una formulaci´on variacional mixta . . . . . . . . . . . . . . 68 4.2.3. Una formulaci´on variacional BEM–FEM . . . . . . . . . . . 70 4.3. Esquema discreto . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.3.1. Descripci´on del problema discreto . . . . . . . . . . . . . . 76 4.3.2. An´alisis del esquema discreto . . . . . . . . . . . . . . . . . 78 4.3.3. An´alisis de la convergencia . . . . . . . . . . . . . . . . . . 82 4.4. Forma matricial del esquema discreto . . . . . . . . . . . . . . . . . 88 4.4.1. Implementaci´on del esquema discreto . . . . . . . . . . . . . 90 4.5. Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . 90

Introducci´on 5. Eddy currents de evoluci´ on 5.1. Formulaci´on variacional . . . . . . . . . . . . . . . . . . 5.1.1. Una formulaci´ on variacional global . . . . . . . . 5.1.2. Una formulaci´ on variacional mixta . . . . . . . . 5.1.3. Una formulaci´ on variacional BEM–FEM . . . . . 5.2. Esquema semidiscreto . . . . . . . . . . . . . . . . . . . 5.2.1. Descripci´on y an´alisis del problema semidiscreto 5.2.2. An´alisis de la convergencia . . . . . . . . . . . . 5.3. Forma matricial del esquema semidiscreto . . . . . . . .

xvii

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

95 96 96 100 103 108 108 110 117

A. Integraci´ on num´ erica

119

Bibliograf´ıa

129

Cap´ıtulo 1

Modelo de eddy currents El comportamiento de un campo electromagn´etico est´ a regido por las ecuaciones de Maxwell, que escribimos a continuaci´ on en forma diferencial:  e,  −∂t D + rot H = J       ∂t B + rot E = 0 ,       

div D = ρ , div B = 0 .

La notaci´on utilizada es la est´andar: D es la inducci´ on el´ectrica, E es el campo el´ectrico, B es la inducci´ on magn´etica, H es el campo magn´etico, e es la densidad de corriente, J ρ es la densidad de carga. 1

(1.0.1)

2

Cap´ıtulo 1

Estas funciones son reales y dependen del tiempo t > 0 y la posici´on x ∈ R3 . La primera ecuaci´on de (1.0.1) es la ley de Amp`ere–Maxwell y coincide con la ley de Amp`ere salvo en el t´ermino adicional ∂t D introducido por Maxwell. Esta ley, reescrita en forma integral, describe c´omo rodean las l´ıneas del campo magn´etico a una superficie que tiene un flujo el´ectrico variable. La segunda ecuaci´on de (1.0.1) se denomina ecuaci´ on de Faraday y describe c´omo rodean las l´ıneas del campo el´ectrico a una superficie por la que pasa un flujo magn´etico variable. La tercera ecuaci´on de (1.0.1) se conoce como ley de Gauss y f´ısicamente significa que el flujo de la inducci´on el´ectrica a trav´es de una superficie cerrada es igual a la carga neta encerrada dentro de dicha superficie. En otras palabras, las l´ıneas del campo el´ectrico empiezan y terminan en cargas el´ectricas. An´alogamente, la cuarta ecuaci´on de (1.0.1) o ley de Gauss del magnetismo afirma que el flujo de la inducci´on magn´etica a trav´es de cualquier superficie cerrada es nulo. En t´erminos de las l´ıneas del campo magn´etico, quiere decir que dichas l´ıneas no convergen ni divergen en ning´ un punto del espacio; es decir, no existen polos magn´eticos aislados. Completamos el sistema (1.0.1) a˜ nadiendo las ecuaciones constitutivas D = εE, B = µH,

(1.0.2)

y la ley de Ohm generalizada e = J + σE, J

(1.0.3)

donde J representa la densidad de la fuente de corriente. Seg´ un estas ecuaciones, el campo electromagn´etico interacciona con la materia a trav´es de los siguientes par´ametros: la permitividad el´ectrica, ε, la permeabilidad magn´etica, µ,

Modelo de eddy currents

3

la conductividad el´ectrica, σ. En efecto, las ecuaciones constitutivas (1.0.2) describen propiedades intr´ınsecas e de la materia. Adem´as, la ley (1.0.3) afirma que la densidad de corriente total, J, consiste en la densidad de la fuente de corriente, J, m´as un t´ermino proporcional al campo el´ectrico E.

1.1.

Modelo de eddy currents

El sistema de ecuaciones (1.0.1–1.0.3) describe un amplio rango de fen´omenos f´ısicos. Restringi´endonos a situaciones particulares, podemos hacer hip´otesis que simplifican las ecuaciones y deducir as´ı modelos que son m´as sencillos de resolver num´ericamente. Entre ellos, uno de los m´as utilizados es el modelo de eddy currents, que surge de las ecuaciones de Maxwell al despreciar las corrientes de desplazamiento, ∂t D, y que es el centro de estudio de este trabajo. M´as concretamente, estudiamos un conductor que ocupa un dominio acotado Ωc ⊂ R3 con frontera Σ := ∂Ωc Lipschitz continua. Suponemos que est´a rodeado por un medio diel´ectrico, por ejemplo aire seco, y denotamos Ωec := R3 \ Ωc . Entonces buscamos el campo electromagn´etico (E, H) que cumpla  rot H = J + σ E en (0, T ) × R3 ,        en (0, T ) × R3 ,  ∂t (µ H) + rot E = 0 (1.1.1) div(ε E) = 0 en (0, T ) × Ωec ,    Z     ε E · n dS = 0 en (0, T ) , ∀j = 1, . . . , NΣ ,  Σj

donde Σj (j = 1, . . . , NΣ ) representan las componentes conexas de Σ. Adem´as, n denota un vector unitario que es normal a Σ en casi todo punto y est´a orientado hacia el interior de Ωec . El campo electromagn´etico tambi´en debe satisfacer las condiciones iniciales H(0, x) = H0 (x) y E(0, x) = E0 (x)

en R3 ,

(1.1.2)

4

Cap´ıtulo 1

as´ı como las condiciones asint´oticas 1 1 H(t, x) = O( ) y E(t, x) = O( ) |x| |x|

cuando |x| → ∞ ,

(1.1.3)

x uniformemente en casi todas las direcciones (t, |x| ), con t ∈ (0, T ) y x ∈ R3 \ {0}.

Nota 1.1.1 En [14] se justifica el problema (1.1.1) como caso l´ımite de la familia uniparam´etrica de problemas  −α ∂t (ε Eα ) + rot Hα = J + σ Eα en R3 ,        ∂t (µ Hα ) + rot Eα = 0 en R3 ,         

Z Σj

en (0, T ) × Ωec ,

div(ε Eα ) = 0 ε Eα · n dS = 0

en (0, T ) , ∀j = 1, . . . , NΣ

En efecto, el problema f´ısico (1.0.1–1.0.3) se corresponde con α = 1 y el problema (1.1.1) se obtiene con un razonamiento asint´ otico para α → 0. En lo que sigue consideramos u ´nicamente materiales sin memoria que sean lineales e is´ otropos. Adem´as, suponemos que los par´ametros ε y µ dependen u ´nicamente de la variable espacial x ∈ R3 . Entonces ε, µ y σ son funciones escalares, reales y acotadas que cumplen las siguientes propiedades: la permitividad el´ectrica es tal que 0 < ε0 ≤ ε(x) ≤ ε1

en c.t.p. x ∈ Ωc ,

ε(x) = ε0

en c.t.p. x ∈ Ωec ;

(1.1.4)

la permeabilidad magn´etica satisface 0 < µ0 ≤ µ(x) ≤ µ1

en c.t.p. x ∈ Ωc ,

µ(x) = µ0

en c.t.p. x ∈ Ωec ;

(1.1.5)

la conductividad el´ectrica cumple 0 < σ0 ≤ σ(t, x) ≤ σ1

en c.t.p. (t, x) ∈ (0, T ) × Ωc ,

σ(t, x) = 0

en c.t.p. (t, x) ∈ (0, T ) × Ωec ;

(1.1.6)

Modelo de eddy currents

5

para algunas constantes ε0 , ε1 , µ0 , µ1 , σ0 , σ1 ∈ R. Bajo estas hip´otesis, deducimos de la primera ecuaci´on de (1.1.1) que la funci´on J debe cumplir las siguientes condiciones de compatibilidad: Z

div J = 0 J · n dS = 0

Σj

1.2.

en (0, T ) × Ωec , en (0, T ) , ∀j = 1, . . . , NΣ .

(1.1.7)

Modelo de eddy currents en r´ egimen arm´ onico

A diferencia de lo que sucede para el problema de evoluci´ on, la modelizaci´on del problema en r´egimen arm´onico emplea funciones complejas. Por esta raz´on, √ umero complejo introducimos la unidad imaginaria ı := −1. Adem´as, dado un n´ α ∈ C, representamos sus partes real e imaginaria por Re[α] y Im[α], respectivamente, mientras que denotamos su conjugado por α y su m´odulo por |α|. Consideramos una bobina por la que circula una corriente alterna a baja frecuencia y que induce un campo electromagn´etico en un conductor pasivo; v´ease la Figura 1.1 a continuaci´on.

J

Ωc Figura 1.1: Geometr´ıa del problema

6

Cap´ıtulo 1 En esta situaci´on, la densidad de la corriente alterna presenta la forma J(t, x) = Re[eıωt j(x)] ,

donde ω denota la frecuencia angular. Adem´as, j es un campo complejo que depende u ´nicamente de x ∈ R3 y se denomina amplitud compleja de la densidad de corriente. Supongamos que la conductividad σ depende u ´nicamente de x ∈ R3 . Entonces los campos el´ectrico y magn´etico tambi´en son arm´onicos en tiempo: E(t, x) = Re[eıωt e(x)]

H(t, x) = Re[eıωt h(x)] .

y

Por lo tanto, el problema (1.1.1) es equivalente a determinar las amplitudes complejas e y h tales que  rot h = j + σ e en R3 ,        en R3 ,  ıωµ h + rot e = 0 (1.2.1) div(ε e) = 0 en Ωec ,    Z     ε e · n dS = 0 ∀j = 1, . . . , NΣ .  Σj

Tambi´en deben satisfacer las condiciones asint´ oticas (1.1.3), que ahora se traducen en 1 1 h(x) = O( ) y e(x) = O( ) cuando |x| → ∞ , (1.2.2) |x| |x| uniformemente en casi todas las direcciones

x |x| ,

con x ∈ R3 \ {0}.

Adem´as, podemos reescribir las condiciones de compatibilidad (1.1.7) como Z

div j = 0 j · n dS = 0

Σj

en Ωec , ∀j = 1, . . . , NΣ .

(1.2.3)

Para cerrar este cap´ıtulo, se˜ nalamos que el problema (1.2.1–1.2.2) es u ´til en la pr´actica para dise˜ nar y estudiar m´aquinas que trabajan a baja frecuencia. Efectivamente, en [6, Teorema 8.1] se prueba que el problema (1.2.1–1.2.2) aproxima al problema f´ısico (1.0.1–1.0.3) con orden 2 respecto a la frecuencia ω.

Cap´ıtulo 2

Herramientas matem´ aticas A lo largo de esta memoria, las letras en negrita representan vectores o funciones vectoriales en C3 . Denotamos la norma eucl´ıdea por |·| y los productos escalar eucl´ıdeo y vectorial por · y ∧, respectivamente. Usamos repetidamente las siguientes identidades vectoriales:

(a)

a ∧ (b ∧ c) = (a · c) b − (a · b) c ,

(b)

a · (b ∧ c) = b · (c ∧ a) = c · (a ∧ b) ,

(c)

rot rot u = ∇(div u) − ∆u ,

(d)

div(v ∧ u) = u · rot v − rot u · v ,

(2.0.1)

para cualesquiera vectores a, b, c ∈ C3 y funciones vectoriales u, v suficientemente regulares. Utilizamos el super´ındice > para designar la transpuesta de una matriz e introducimos una norma matricial k·k que est´e inducida por una norma vectorial. Adem´as, C (con o sin sub´ındices) es una constante gen´erica, no necesariamente la misma en diferentes estimaciones.

7

8

Cap´ıtulo 2

2.1.

Espacios de Sobolev cl´ asicos

Sea O un dominio abierto y acotado de Rn (n = 2 ´o 3). Suponemos que su frontera ∂O es Lipschitz continua, concepto que detallamos a continuaci´ on; v´ease [45, Definici´on 3.1]. Definici´ on 2.1.1. Decimos que ∂O es Lipschitz continua si para cada x ∈ ∂O podemos tomar un sistema de coordenadas ortogonal y = (y1 , . . . , yn−1 , yn ) ≡ (y 0 , yn ) , para el que existen un vector a ∈ Rn con x ∈ U := {y ∈ Rn ; |yj | < aj ∀j = 1, . . . , n}, an una funci´ on f : U 0 → R Lipschitz continua con |f (y 0 )| < para cualquier 2 n o y 0 ∈ U 0 := y 0 ∈ Rn−1 ; |yj0 | < aj ∀j = 1, . . . , n − 1 , cumpliendo que © ª O ∩ U = y ∈ U ; yn < f (y 0 )

y

© ª ∂O ∩ U = y ∈ U ; yn = f (y 0 ) .

Denotamos por D(O) y D(Rn ) los espacios de las funciones indefinidamente diferenciables que tienen soporte compacto contenido en O ´o en Rn , respectivamente. Entonces definimos © ª D(O) := f |O ; f ∈ D(R3 )

y

¡ ¢3 D(O) := D(O) .

Para cada p ∈ [1, ∞), consideramos el espacio Lp (O) de funciones f : O → C que son medibles en el sentido de Lebesgue y con Z p kf kp,O := |f (x)|p dx < +∞ . Ω

Dado un multi´ındice α = (α1 , . . . , αn )> ∈ Nn , escribimos ∂ α :=

∂ αn ∂ α1 α1 · · · ∂x1 ∂xαnn

y

|α| := α1 + · · · + αn .

Herramientas matem´aticas

9

A continuaci´on introducimos los espacios de Sobolev cl´ asicos en el dominio O: En primer lugar, definimos los espacios de Sobolev de ´ındice entero no negativo. Despu´es los generalizamos a ´ındices reales no negativos y, finalmente, a ´ındices reales negativos.

2.1.1.

´Indice entero no negativo

Sea m ∈ N un entero no negativo. Introducimos el espacio de Sobolev © ª H m (O) := f ∈ L2 (O) ; ∂ α f ∈ L2 (O) ∀α ∈ Nn con |α| ≤ m , que es la compleci´on de D(O) con respecto a la norma  kf km,O := 

X Z

1/2 |∂ α f (x)|2 dx

.

|α|≤m O

En particular, H 0 (O) ≡ L2 (O) y utilizamos la siguiente notaci´on est´andar para su producto interior: Z (f, g)0,O := f (x) g(x) dx . O

2.1.2.

´Indice real no negativo

Existen varias definiciones equivalentes de los espacios de Sobolev de ´ındice real no negativo. Aqu´ı los introducimos con el m´etodo de interpolaci´ on real; cf. [29, P´arrafo 1.1]. M´as concretamente, representamos cada n´ umero s ∈ [0, +∞) no entero bajo la forma s = m + sˆ, con m ∈ N y sˆ ∈ (0, 1). Entonces definimos el espacio de Sobolev H s (O) como la compleci´on de D(O) con respecto a la norma ´1/2 ³ X kf ks,O := kf k2m,O + |∂ α f |2sˆ,O , |α|=m

donde

Z Z |g|2sˆ,O :=

O

O

|g(x) − g(y)|2 dx dy . |x − y|n+2ˆs

10

2.1.3.

Cap´ıtulo 2

´Indice real negativo

Introducimos los espacios de Sobolev de ´ındice negativo por dualidad. En efecto, fijado un n´ umero real s no negativo, denotamos por H0s (O) la clausura de D(O) en H s (O). Entonces H0s (O) es denso en L2 (O) y podemos definir H −s (O) como el espacio dual de H0s (O) con respecto al espacio pivote L2 (O). Representamos su producto de dualidad como sigue: H −s (O) × H0s (O) → C (f, g) 7→ hf, gis,O ; en particular, tenemos que hf, gis,O = (f, g)0,O

∀f ∈ L2 (O), g ∈ H0s (O) .

Adem´as, utilizamos la norma inducida en el espacio H −s (O): kf k−s,O :=

sup

g∈H0s (O), g6=0

|hf, gis,O | , kgks,O

as´ı que las siguientes inclusiones son continuas y densas: H0s (O) ,→ L2 (O) ,→ H −s (O) . En lo sucesivo, denotamos con negrita los espacios vectoriales asociados a los espacios de Sobolev anteriores. De este modo, tomamos L2 (O) := (L2 (O))n ; an´alogamente, Hs (O) := (H s (O))n y Hs0 (O) := (H0s (O))n para todo s ∈ R.

2.2.

Espacios de Sobolev sobre superficies y teor´ıa de trazas

Consideramos un dominio acotado Ω ⊆ R3 cuya frontera Γ := ∂Ω sea Lipschitz continua. Sin p´erdida de generalidad, suponemos que Ω, Ωe := R3 \ Ω y Γ son conexos.

Herramientas matem´aticas

11

Ωe Γ Ω

Figura 2.1: Dominio Ω Denotamos por n el vector unitario normal a la superficie Γ y orientado hacia el interior de Ωe . Sea {Ωj ; j = 0, . . . , m} un recubrimiento abierto de Ω tal que m [

Γ⊆

Ωj

y

Γ ∩ Ω0 = ∅ .

j=1

Para cada j = 1, . . . , m, definimos la aplicaci´on φj : Uj0

→ Γj

y 0 7→ φj (y 0 ) := (y 0 , fj (y 0 )) , donde seguimos la notaci´on de la Definici´on 2.1.1 y Γj := Γ ∩ Ωj . As´ı construido, © ª (Γj , φj ) ; j = 1, . . . , m es un atlas de Γ.

2.2.1.

Espacios de Sobolev sobre superficies

Para cada s ∈ [0, 1], definimos el espacio de Sobolev n o H s (Γ) := λ ∈ D0 (Γ) ; λ ◦ φj ∈ H s (Uj0 ∩ φ−1 (Γ )) ∀j = 1, . . . , m j j

(2.2.1)

y lo dotamos de la norma kλks,Γ

 m X :=  kλ ◦ φj k2 j=1

1/2 

s, Uj0 ∩φ−1 j (Γj )

.

(2.2.2)

12

Cap´ıtulo 2

Puede comprobarse que el espacio H s (Γ) es de Hilbert y que su definici´on es independiente de la elecci´on del recubrimiento {Ωj ; j = 0, . . . , m}. Adem´as, H s (Γ) es denso en H 0 (Γ) ≡ L2 (Γ). Esta propiedad nos permite introducir H −s (Γ) como el espacio dual de H s (Γ) con espacio pivote L2 (Γ). Representamos su producto de dualidad como sigue: H −s (Γ) × H s (Γ) → C (λ, η) 7→ hλ, ηis,Γ ; en particular, tenemos que hλ, ηis,Γ = (λ, η)0,Γ

∀λ ∈ L2 (Γ), η ∈ H s (Γ) .

Consideramos la norma inducida en H −s (Γ): kλk−s,Γ :=

sup η∈H s (Γ), η6=0

|hλ, ηis,Γ | , kηks,Γ

de forma que las siguientes inclusiones son continuas y densas: H s (Γ) ,→ L2 (Γ) ,→ H −s (Γ) .

(2.2.3)

Denotamos con negrita los espacios vectoriales asociados; es decir, tomamos := (H s (Γ))3 para cualquier s ∈ [−1, 1]. Tambi´en introducimos

Hs (Γ)

© ª L2τ (Γ) := λ ∈ L2 (Γ) ; λ · n = 0 sobre Γ , ¡ ¢3 que es un subespacio cerrado de L2 (Γ) := L2 (Γ) . As´ı definido, L2τ (Γ) se identifica con el espacio de funciones de L2 (Γ) que pertenecen al fibrado tangente de Γ en casi todo punto; al respecto, v´eanse [20, P´arrafos 1 y 2].

2.2.2.

Trazas en espacios de Sobolev

El siguiente resultado es cl´asico y puede consultarse en [29, Teorema 1.5] ´o en [35, Teorema 3.37].

Herramientas matem´aticas

13

Teorema 2.2.1. Dado s ∈ (0, 1/2], la aplicaci´ on que a cada funci´ on u ∈ D(Ω) le asocia su restricci´ on u|Γ a Γ se extiende de forma u ´nica a un operador lineal continuo γ : H s+1/2 (Ω) → H s (Γ) . Este operador se denomina traza y es suprayectivo. Adem´ as, admite una inversa −1 a derecha que es continua y que denotamos por γ . Fijado s ∈ (0, 1/2], introducimos la traza vectorial γ : Hs+1/2 (Ω) → Hs (Γ) actuando por componentes. Tambi´en definimos el siguiente operador continuo, que se dice traza tangencial : γτ : Hs+1/2 (Ω) → L2τ (Γ) q 7→ γτ q := γq ∧ n .

(2.2.4)

Denotamos su imagen por Hs× (Γ) := γτ ( Hs+1/2 (Ω)), que es un espacio de Hilbert con la norma del grafo n o kλkHs× (Γ) := ´ınf kqks+1/2,Ω ; q ∈ Hs+1/2 (Ω) con γτ q = λ . El siguiente resultado se cumple por la propia definici´on del espacio Hs× (Γ). Proposici´ on 2.2.2. Para todo s ∈ (0, 1/2], el operador lineal γτ : Hs+1/2 (Ω) → Hs× (Γ) es continuo y suprayectivo. Adem´ as, admite una inversa a derecha que es continua −1 y que representamos por γ τ . Para cada s ∈ (0, 1/2], sabemos que H s (Γ) es denso en L2 (Γ) (v´ease (2.2.3)), as´ı que Hs× (Γ) es denso en L2τ (Γ). En consecuencia, podemos introducir el espacio s 2 H−s en podemos × (Γ) dual de H× (Γ) con respecto al espacio pivote Lτ (Γ). Tambi´ −s s extender la siguiente forma sesquilineal al espacio H× (Γ) × H× (Γ): hλ, ηiτ ,Γ := (λ ∧ n, η)0,Γ

∀λ, η ∈ L2τ (Γ) .

14

Cap´ıtulo 2

Entonces

s H−s × (Γ) × H× (Γ) → C

(λ, η) 7→ hλ, ηiτ ,Γ s define un producto de dualidad entre H−s × (Γ) y H× (Γ). Dotamos al espacio H−s × (Γ) de la norma inducida

kλkH−s (Γ) := ×

sup

η∈Hs× (Γ), η6=0

|hλ, ηiτ ,Γ | . kηkHs× (Γ)

De esta forma, las siguientes inclusiones son continuas y densas: Hs× (Γ) ,→ L2τ (Γ) ,→ H−s × (Γ) .

(2.2.5)

En lo sucesivo, denotamos por convenio H0× (Γ) := Lτ2 (Γ). Adem´as, simplificamos la notaci´on escribiendo γτ para la composici´on γτ ◦ γ −1 .

2.3.

Espacios H(div, Ω) y H(rot, Ω)

Consideramos el espacio © ª H(div, Ω) := u ∈ L2 (Ω) ; div u ∈ L2 (Ω) , que es de Hilbert con la norma del grafo ¡ ¢1/2 kukdiv,Ω := kuk20,Ω + kdiv uk20,Ω . Teorema 2.3.1. El espacio D(Ω) es denso en H(div, Ω). Demostraci´ on. La prueba de este resultado puede consultarse en [29, Teorema 2.4] ´o [45, Teorema 3.22]. ¥ Aplicando el Teorema 2.3.1 junto con la f´ormula de Green (u, ∇v)0,Ω + (div u, v)0,Ω = (u|Γ · n, γv)0,Γ

∀u ∈ D(Ω), v ∈ H 1 (Ω) ,

Herramientas matem´aticas

15

deducimos las siguientes propiedades del operador de traza normal; para una exposici´on m´as detallada, v´ease [29, Teorema 2.5 y Corolario 2.6] ´o bien [45, Teorema 3.24]. Proposici´ on 2.3.2. La aplicaci´ on que a cada funci´ on u ∈ D(Ω) le asocia su componente normal u|Γ ·n sobre Γ se extiende a un u ´nico operador lineal continuo γn : H(div, Ω) → H −1/2 (Γ) , que es suprayectivo y se dice traza normal. Adem´ as, se cumple (u, ∇v)0,Ω + (div u, v)0,Ω = hγn u, γvi1/2,Γ ,

(2.3.1)

para cualesquiera funciones u ∈ H(div, Ω) y v ∈ H 1 (Ω). El siguiente resultado es consecuencia de la f´ormula de Green (2.3.1). Corolario 2.3.3. Consideramos Ω1 , Ω2 ⊆ R3 dos dominios (abiertos pero no necesariamente acotados) con fronteras Γ1 := ∂Ω1 , Γ2 := ∂Ω2 Lipschitz continuas y denotamos por Ω el interior de Ω1 ∪ Ω2 . Sea u ∈ L2 (Ω) una funci´ on tal que u|Ω1 ∈ H(div, Ω1 ) y u|Ω2 ∈ H(div, Ω2 ). Entonces u ∈ H(div, Ω) si y s´ olo si γn (u|Ω1 ) = γn (u|Ω2 )

sobre Γ1 ∩ Γ2 .

Definimos el espacio de funciones vectoriales © ª H(rot, Ω) := u ∈ L2 (Ω) ; rot u ∈ L2 (Ω) , que es de Hilbert con la norma k·krot,Ω inducida por el producto escalar (u, v)rot,Ω := (u, v)0,Ω + (rot u, rot v)0,Ω . En general, para cada s ≥ 0 definimos el espacio Hs (rot, Ωc ) := {u ∈ Hs (Ωc ) ; rot u ∈ Hs (Ωc )} , que es de Hilbert con la norma del grafo ¡ ¢1/2 kukHs (rot,Ωc ) := kuk2s,Ωc + krot uk2s,Ωc . El siguiente resultado es cl´asico y puede consultarse en [29, Teorema 2.10] o´ en [45, Teorema 3.26].

16

Cap´ıtulo 2

Teorema 2.3.4. El espacio D(Ω) es denso en H(rot, Ω). Obs´ervese que, seg´ un la cuarta identidad de (2.0.1) y la f´ormula de Gauss, Z Z (u · rot v − rot u · v) dx = (v ∧ u) · n dSx ∀u, v ∈ D(Ω) . Ω

Γ

Aplicando las dos primeras identidades de (2.0.1), tenemos que (u, rot v)0,Ω − (rot u, v)0,Ω = ((u ∧ n) ∧ n, v ∧ n)0,Γ

∀u, v ∈ D(Ω) .

Entonces, con un razomiento por densidad deducimos que (u, rot v)0,Ω − (rot u, v)0,Ω = hγτ u, γτ viτ ,Γ

∀u, v ∈ H1 (Ω) .

(2.3.2)

Utilizando esta relaci´on junto con el Teorema 2.3.4, podemos extender el operador de traza tangencial γτ a todo el espacio H(rot, Ω), tal y como se˜ nalamos a continuaci´on; cf. [29, Teorema 2.11] ´o [45, Teorema 3.29]. Teorema 2.3.5. La traza tangencial definida por (2.2.4) se extiende de forma u ´nica a un operador lineal continuo −1/2

γτ : H(rot, Ω) → H×

(Γ) .

Adem´ as, se cumple (u, rot v)0,Ω − (rot u, v)0,Ω = hγτ u, γτ viτ ,Γ ,

(2.3.3)

para cualesquiera u ∈ H(rot, Ω) y v ∈ H1 (Ω). El siguiente resultado es consecuencia de la f´ormula de Green (2.3.3). Corolario 2.3.6. Consideramos Ω1 , Ω2 ⊆ R3 dos dominios (abiertos pero no necesariamente acotados) con fronteras Γ1 := ∂Ω1 , Γ2 := ∂Ω2 Lipschitz continuas y denotamos por Ω el interior de Ω1 ∪ Ω2 . Sea u ∈ L2 (Ω) una funci´ on tal que u|Ω1 ∈ H(rot, Ω1 ) y u|Ω2 ∈ H(rot, Ω2 ). Entonces u ∈ H(rot, Ω) si y s´ olo si γτ (u|Ω1 ) = γτ (u|Ω2 )

sobre Γ1 ∩ Γ2 .

Herramientas matem´aticas

2.4.

17

Operadores diferenciales sobre superficies

Introducimos el rotacional vectorial de superficie rotΓ := γτ ◦ ∇ ◦ γ −1 , definici´on que queda justificada gracias al siguiente resultado auxiliar. Lema 2.4.1. Sea u ∈ H 1 (Ω) tal que γu = 0. Entonces γτ (∇u) = 0. Demostraci´ on. Fijamos v ∈ H1 (Ω) arbitrario. Por un lado, la f´ormula (2.3.3) para ∇u ∈ H(rot, Ω) y v ∈ H1 (Ω) garantiza que (∇u, rot v)0,Ω = hγτ (∇u), γτ viτ ,Γ . Por otro lado, aplicando la f´ormula (2.3.1) para las funciones rot v ∈ H(div, Ω) y u ∈ H 1 (Ω) junto con la hip´otesis γu = 0, (∇u, rot v)0,Ω = −(u, div(rot v))0,Ω + hγu, γn (rot v)i1/2,Γ = 0 . En consecuencia, hγτ (∇u), γτ viτ ,Γ = 0 ∀v ∈ H1 (Ω) , −1/2

y concluimos que γτ (∇u) = 0 en H×

(Γ).

¥

Proposici´ on 2.4.2. El operador lineal −1/2

rotΓ : H 1/2 (Γ) → H×

(Γ)

est´ a bien definido y es continuo. Demostraci´ on. N´otese en primer lugar que la definici´on de rotΓ es independiente de la elecci´on del levantamiento γ −1 . En efecto, dados u1 , u2 ∈ H 1 (Ω) tales que γu1 = γu2 , tenemos que u1 − u2 ∈ H 1 (Ω) con γ(u1 − u2 ) = 0; luego, seg´ un el Lema 2.4.1, γτ (∇u1 ) = γτ (∇u2 ) . En lo que respecta a la continuidad de rotΓ , basta observar que es composici´on de aplicaciones continuas (v´eanse el Teorema 2.2.1 y la Proposici´on 2.3.5). ¥

18

Cap´ıtulo 2 El operador dual de rotΓ con respecto al producto de dualidad −1/2

H×

1/2

(Γ) × H× (Γ) → C (λ, η) 7→ hλ, ηiτ ,Γ

se conoce como divergencia de superficie y se denota por divΓ : hϕ, divΓ λi1/2,Γ = hrotΓ ϕ, λiτ ,Γ

1/2

∀λ ∈ H× (Γ) , ϕ ∈ H 1/2 (Γ) .

(2.4.1)

An´alogamente, el operador dual de rotΓ con respecto al producto de dualidad −1/2

H×

1/2

(Γ) × H× (Γ) → C (λ, η) 7→ hλ, ηi1/2,Γ

se denomina rotacional escalar de superficie y se designa por rotΓ : hϕ, rotΓ λi1/2,Γ = hrotΓ ϕ, λi1/2,Γ

1/2

∀λ ∈ H× (Γ) , ϕ ∈ H 1/2 (Γ) .

(2.4.2)

Observemos que la Proposici´on 2.4.2 justifica estas definiciones y garantiza que los operadores lineales 1/2

divΓ : H× (Γ) → H −1/2 (Γ)

y

1/2

rotΓ : H× (Γ) → H −1/2 (Γ)

son continuos.

2.5.

Trazas tangenciales en H(rot, Ω) −1/2

Dada una funci´on λ ∈ H× ϕ ∈ H −1/2 (Γ) tal que

(Γ), decimos que divΓ λ ∈ H −1/2 (Γ) si existe

hγu, ϕi1/2,Γ = hγτ (∇u), λiτ ,Γ

∀u ∈ H 2 (Ω) ,

en cuyo caso tomamos divΓ λ := ϕ. Consideramos el espacio n o −1/2 −1/2 H× (divΓ , Γ) := λ ∈ H× (Γ) ; divΓ λ ∈ H −1/2 (Γ) ,

Herramientas matem´aticas

19

que es de Hilbert con la norma del grafo −1/2

H×

(divΓ , Γ) → R µ λ 7→ kλk2

¶1/2 −1/2

H×

(Γ)

+

kdivΓ λk2−1/2,Γ

.

Proposici´ on 2.5.1. Dada una funci´ on u ∈ H(rot, Ω), se cumple en H −1/2 (Γ) .

divΓ (γτ u) = γn (rot u)

Demostraci´ on. Sean u ∈ D(Ω) y v ∈ H 2 (Ω) cualesquiera. Como ∇v ∈ H1 (Ω), podemos aplicar la f´ormula de Green (2.3.3) junto con las definiciones de los operadores rotΓ y divΓ : (∇v, rot u)0,Ω = hrotΓ (γv), γτ uiτ ,Γ = hγv, divΓ (γτ u)i1/2,Γ . Pero tambi´en podemos integrar por partes con la f´ormula (2.3.1), de modo que (rot u, ∇v)0,Ω = hγn (rot u), γvi1/2,Γ . Utilizando las dos u ´ltimas identidades, tenemos hdivΓ (γτ u), γvi1/2,Γ = hγn (rot u), γvi1/2,Γ

∀v ∈ H 2 (Ω) ,

y deducimos que divΓ (γτ u) = γn (rot u) ∈ H −1/2 (Γ) ∀u ∈ D(Ω) .

(2.5.1)

Recordemos que las Proposiciones 2.3.2 y 2.3.5 garantizan que los operadores γn ◦ rot : H(rot, Ω) → H −1/2 (Γ)

−1/2

y

γτ : H(rot, Ω) → H×

(Γ)

son continuos. En consecuencia, el Teorema 2.3.4 nos permite concluir el resultado del enunciado haciendo un razonamiento por densidad en (2.5.1). ¥ Teorema 2.5.2. El operador lineal −1/2

γτ : H(rot, Ω) → H×

(divΓ , Γ)

es continuo y suprayectivo. Adem´ as, admite una inversa a derecha que es continua −1 y que denotamos por γτ .

20

Cap´ıtulo 2

Demostraci´ on. Fijada u ∈ H(rot, Ω) cualquiera, la Proposici´on 2.5.1 garantiza que divΓ (γτ u) = γn (rot u) ∈ H −1/2 (Γ). En virtud de la Proposici´on 2.3.2 sabemos que el operador γn ◦ rot : H(rot, Ω) → H −1/2 (Γ) es continuo, as´ı que kdivΓ (γτ u)k−1/2,Γ = kγn (rot u)k−1/2,Γ ≤ C kukrot,Ω . Utilizando esta propiedad junto con el Teorema 2.3.5, deducimos que el operador −1/2

γτ : H(rot, Ω) → H×

(divΓ , Γ)

est´a bien definido y es continuo. Adem´as, en [20, Teorema 4.1] se demuestra que este operador es suprayectivo. Entonces podemos aplicar el teorema del grafo cerrado y concluimos que el operador γτ admite una inversa a derecha −1/2

γτ−1 : H×

(divΓ , Γ) → H(rot, Ω) ,

que es continua.

2.6.

¥

Inversa del operador rotacional de superficie en dominios simplemente conexos

A lo largo de este p´arrafo, suponemos que el dominio Ω est´a acotado y es simplemente conexo. En esta situaci´on, estudiamos c´omo invertir el operador −1/2

rotΓ : H 1/2 (Γ) → H×

(Γ) .

Con este fin, caracterizamos su n´ ucleo y su imagen como en [20, Corolario 5.3], y a continuaci´on aplicamos el teorema del grafo cerrado. En primer lugar, recordamos un teorema de representaci´ on que es cl´asico y puede consultarse en [29, Teorema 3.4 y Corolario 3.4]; ´o bien en [45, Teoremas 3.37 y 3.38].

Herramientas matem´aticas

21

Teorema 2.6.1. Si u ∈ H(div, Ω) cumple div u = 0 en Ω

y

hγn u, 1i1/2,Γ = 0 ,

entonces existe un potencial vectorial v ∈ H1 (Ω) con div v = 0 y rot v = u. An´ alogamente, si u ∈ H(rot, Ω) satisface rot u = 0 en Ω , entonces existe un potencial escalar v ∈ H 1 (Ω) tal que ∇v = u y este potencial es u ´nico salvo una constante aditiva. A continuaci´on recordamos un resultado auxiliar sobre el operador divergen© ª cia. Para ello, definimos el espacio L20 (Ω) := f ∈ L2 (Ω) ; (f, 1)0,Ω = 0 de las funciones de cuadrado integrable y media nula. Obs´ervese que el operador div : H10 (Ω) → L20 (Ω) es continuo y suprayectivo. Denotamos su n´ ucleo por ker(div) e introducimos © ª ker(div)⊥ := u ∈ H10 (Ω) ; (u, v)0,Ω + (∇u, ∇v)0,Ω = 0 ∀v ∈ ker(div) , que es el subespacio de H10 (Ω) ortogonal a ker(div). Lema 2.6.2. El siguiente operador lineal define un isomorfismo: div : ker(div)⊥ → L20 (Ω) . Demostraci´ on. Puede consultarse en [29, Corolario 2.4].

¥ −1/2

Teorema 2.6.3. El n´ ucleo del operador rotΓ : H 1/2 (Γ) → H×

(Γ) es

ker(rotΓ ) = C . Demostraci´ on. Para comodidad del lector, reproducimos a continuaci´ on la demostraci´on dada en [20, Corolario 3.7]. Por un lado, la inclusi´on ker(rotΓ ) ⊇ C es trivial.

22

Cap´ıtulo 2

Por otro lado, para estudiar la inclusi´on rec´ıproca, consideramos una funci´on p ∈ H 1 (Ω) con rotΓ (γp) = 0. Sea v ∈ H1 (Ω) arbitraria. Aplicando la f´ormula de Green (2.3.3) a las funciones ∇p ∈ H(rot, Ω) y v ∈ H1 (Ω), (∇p, rot v)0,Ω = hrotΓ (γp), γτ viτ ,Γ = 0 . Adem´as, utilizando la f´ormula de integraci´ on por partes (2.3.1) para las funciones 1 rot v ∈ H(div, Ω) y p ∈ H (Ω), (rot v, ∇p)0,Ω = hγn (rot v), γpi1/2,Γ . Teniendo en cuenta estas dos propiedades, resulta que hγn (rot v), γpi1/2,Γ = 0 ∀v ∈ H1 (Ω) . En consecuencia, habr´ıamos concluido el resultado si prob´asemos que © ª −1/2 γn (rot v) ; v ∈ H1 (Ω) = H0 (Γ) , © ª −1/2 donde H0 (Γ) := η ∈ H −1/2 (Γ) ; hη, 1i1/2,Γ = 0 . Mostraremos esta identidad por doble contenido: La inclusi´on (⊆) es consecuencia directa de la f´ormula de Green (2.3.1). En efecto, basta notar que hγn (rot v), 1i1/2,Γ = (rot v, ∇1)0,Ω + (div(rot v), 1)0,Ω = 0 ∀v ∈ H1 (Ω) . −1/2

Para ver la inclusi´on (⊇), fijamos λ ∈ H0 (Γ) arbitrario. Recordemos que, en virtud de la Proposici´on 2.3.2, existe alg´ un u ∈ H(div, Ω) con γn u = λ. Entonces, integrando por partes con la f´ormula (2.3.1), (div u, 1)0,Ω = hγn u, 1i1/2,Γ = hλ, 1i1/2,Γ = 0 . Es decir, div u ∈ L20 (Ω) y, aplicando el Lema 2.6.2, existe alg´ un w ∈ H10 (Ω) tal que div w = div u. N´otese que γn w = 0, luego hγn (u − w), 1i1/2,Γ = hλ, 1i1/2,Γ = 0 . Como adem´as div(u − w) = 0, el Teorema 2.6.1 garantiza la existencia de alg´ un 1 v ∈ H (Ω) con rot v = u − w y por lo tanto λ = γn u = γn (rot v). ¥

Herramientas matem´aticas

23

Lema 2.6.4. Si u ∈ H 1 (Ω), entonces divΓ ◦ rotΓ (γu) = 0 .

(2.6.1) −1/2

Demostraci´ on. Sea u ∈ D(Ω). N´otese que rotΓ (γu) ∈ H× hγv, divΓ ◦ rotΓ (γu)i1/2,Γ = hγτ (∇v), γτ (∇u)iτ ,Γ

(divΓ , Γ) con

∀v ∈ H 1 (Ω) .

Integrando por partes con la f´ormula de Green (2.3.3), deducimos que hγv, divΓ ◦ rotΓ (γu)i1/2,Γ = 0 ∀v ∈ H 1 (Ω) . Puesto que γ(H 1 (Ω)) = H 1/2 (Γ) (v´ease el Teorema 2.2.1), la identidad anterior significa que (2.6.1) se cumple siempre que u ∈ D(Ω). Finalmente, recordemos que D(Ω) es denso en H 1 (Ω) y que el operador −1/2 rotΓ◦γ : H 1 (Ω) → H× (Γ) es continuo (v´eanse el Teorema 2.2.1 y la Proposici´on 2.4.2). Entonces podemos utilizar un razonamiento por densidad y concluir que (2.6.1) tambi´en se cumple para cada funci´on u ∈ H 1 (Ω). ¥ El resultado que presentamos a continuaci´ on puede consultarse en [16, Proposici´on 2.1] para el caso de poliedros de Lipschitz y en [20, Corolario 5.3] para el caso m´as general de dominios de Lipschitz. Teorema 2.6.5. La imagen del operador rotΓ se caracteriza por −1/2

rotΓ (H 1/2 (Γ)) = ker(divΓ ) ∩ H×

(divΓ , Γ) .

Demostraci´ on. La siguiente demostraci´on es an´aloga a la de [20, Teorema 5.1]. −1/2

La inclusi´on rotΓ (H 1/2 (Γ)) ⊆ ker(divΓ ) ∩ H× directa del Lema 2.6.4.

(divΓ , Γ) es consecuencia

−1/2

Para ver la inclusi´on rec´ıproca, fijamos λ ∈ H× (divΓ , Γ) con divΓ λ = 0. En virtud del Teorema 2.5.2, existe alg´ un u ∈ H(rot, Ω) con γτ u = λ .

(2.6.2)

24

Cap´ıtulo 2

Entonces, utilizando la Proposici´on 2.5.1, tenemos que γn (rot u) = divΓ (γτ u) = divΓ λ = 0 . Definimos la funci´on  rot u g := 0

en Ω , en Ωe := R3 \ Ω .

N´otese que g|Ω ∈ H(div, Ω) y g|Ωe ∈ H(div, Ωe ) con div g = 0 en Ω ∪ Ωe ; y adem´as γn (g|Ω ) = γn (g|Ωe ) = 0 sobre Γ. Luego el Corolario 2.3.3 garantiza g ∈ H(div, R3 ) con div g = 0 en R3 . Sea O ⊂ R3 un dominio abierto, acotado y simplemente conexo tal que Ω ⊂ O. Como g|O ∈ H(div, O) con div(g|O ) = 0 en O y γn (g|O ) = 0 sobre ∂O, el Teorema 2.6.1 implica que g|O = rot v para alg´ un v ∈ H1 (O). En particular, tenemos que rot v|O∩Ωe = 0 y, aplicando de nuevo el Teorema 2.6.1, existe alg´ un 1 e ϕ ∈ H (O ∩ Ω ) con ∇ϕ = v|O∩Ωe . N´otese que rotΓ (γϕ) = γτ (∇ϕ) = γτ (v|O∩Ωe ) = γτ (v|Ω ) . Definimos w := v|Ω − ∇ϕ ∈ H(rot, Ω) ,

(2.6.3)

que cumple rot w = rot v|Ω = g|Ω = rot u . Entonces u−w ∈ H(rot, Ω) con rot(u−w) = 0 y, utilizando de nuevo el Teorema 2.6.1, sabemos que existe alg´ un ψ ∈ H 1 (Ω) de forma que ∇ψ = u − w .

Ahora aplicamos sucesivamente las identidades (2.6.2), (2.6.4) y (2.6.3): λ = γτ u = γτ w + γτ (∇ψ) = γτ (v|Ω ) − γτ (∇ϕ) + γτ (∇ψ) .

(2.6.4)

Herramientas matem´aticas

25

Por lo tanto, recordando la propiedad γτ (v|Ω ) = rotΓ (γϕ) junto con la definici´on del operador rotΓ , concluimos que λ = rotΓ (γψ) . ¥ El siguiente resultado es consecuencia de los dos teoremas anteriores y del teorema del grafo cerrado; cf. [41, Proposici´on 2.5]. Corolario 2.6.6. El operador lineal −1/2

rotΓ : H 1/2 (Γ)/C → ker(divΓ ) ∩ H×

(divΓ , Γ)

define un isomorfismo.

2.7.

Potenciales arm´ onicos en regiones no acotadas

A lo largo de este apartado, suponemos que el dominio Ω es acotado y simplemente conexo. Denotamos Ωe := R3 \ Ω, e introducimos el espacio ( ) ϕ 1 e 0 e 2 e 2 e W (Ω ) := ϕ ∈ D (Ω ) ; p ∈L Ã (Ω ) , ∇ϕ ∈ L (Ω ) . 1 + |x|2 En [46, P´arrafo 2.5.4] se prueba que este espacio es de Hilbert con la norma W 1 (Ωe ) → R ϕ 7→ k∇ϕk0,Ωe . En dicha referencia tambi´en se demuestra que © ª ϕ|O ; ϕ ∈ W 1 (Ωe ) = H 1 (O) , para cada subdominio O ⊂ Ωe abierto y acotado; as´ı que el siguiente resultado es consecuencia del Teorema 2.2.1.

26

Cap´ıtulo 2

Lema 2.7.1. La aplicaci´ on que a cada funci´ on u ∈ D(Ωe ) le asocia su restricci´ on u|Γ a Γ se extiende de forma u ´nica a un operador lineal continuo γ : W 1 (Ωe ) → H 1/2 (Γ) . Esta aplicaci´ on es suprayectiva y se dice traza. Lema 2.7.2. Sea u ∈ H(div, Ωe ) ∩ H(rot, Ωe ) con div u = 0

y

rot u = 0

en Ωe .

Entonces existe un u ´nico potencial arm´ onico ϕ ∈ W 1 (Ωe ) tal que ∇ϕ = u en Ωe . Demostraci´ on. Planteamos el problema variacional   hallar ϕ ∈ W 1 (Ωe ) tal que  (∇ϕ, ∇ψ) e = −hγ u, γψi n 0,Ω 1/2,Γ

∀ψ ∈ W 1 (Ωe ) .

(2.7.1)

N´otese que, en virtud de la Proposici´on 2.3.2 y el Lema 2.7.1, la aplicaci´on lineal W 1 (Ωe ) → C ψ 7→ hγn u, γψi1/2,Γ es continua. Luego el Lema de Lax–Milgram garantiza que el problema (2.7.1) tiene soluci´on y que ´esta es u ´nica. Concluimos observando que ϕ ∈ W 1 (Ωe ) es soluci´on del problema (2.7.1) si y s´olo si ∇ϕ = u en Ωe , en cuyo caso tenemos que 4ϕ = div u = 0 en Ωe . ¥ En [29, Teorema 2.9] se debilitan las hip´otesis del Lema 2.7.2 como en el siguiente resultado. Teorema 2.7.3. Consideramos una regi´ on O ⊂ R3 simplemente conexa (no necesariamente acotada) y con frontera Lipschitz continua. Sea u ∈ H(rot, O). Entonces rot u = 0 en O si y s´ olo si existe un potencial ϕ ∈ W 1 (O) con ∇ϕ = u en O, en cuyo caso este potencial es u ´nico.

Herramientas matem´aticas

27

1 , que es la soluci´on fundamental del 4π|x| operador de Laplace en R3 . El siguiente teorema de representaci´ on puede consultarse en [46, Teorema 3.1.1]. Introducimos la funci´on Φ(x) :=

Teorema 2.7.4. Sea ϕ ∈ W 1 (Ωe ) una funci´ on arm´ onica. Entonces Z Z ∂Φ(x − y) ∂ϕ ϕ(x) = γϕ(y) dSy − Φ(x − y) (y) dSy ∀x ∈ Ωe . (2.7.2) ∂n ∂n y Γ Γ

2.7.1.

Operadores integrales sobre superficies

Definimos los potenciales de simple y doble capa sobre la superficie Γ: Z Vλ(x) := Φ(x − y) λ(x) dSy ∀x ∈ Γ , Γ

Z Kλ(x) :=

Γ

∂Φ(x − y) λ(y) dSy ∂ny

∀x ∈ Γ .

Tambi´en introducimos el operador hipersingular sobre Γ: ·Z ¸ ∂ ∂Φ(x − y) N λ(x) := − λ(y) dSy ∂nx Γ ∂ny

∀x ∈ Γ .

Las siguientes propiedades de los operadores de simple y doble capa son cl´asicas y pueden consultarse en [25, Teoremas 1 y 3]. Teorema 2.7.5. Dado s ∈ [−1/2, 1/2], los operadores lineales V : H −1/2+s (Γ) → H 1/2+s (Γ)

y

K : H 1/2+s (Γ) → H 1/2+s (Γ)

son continuos. Adem´ as, existe una constante α1 > 0 de forma que hη, Vηi1/2,Γ ≥ α1 kηk2−1/2,Γ , para todo η ∈ H −1/2 (Γ). Tambi´en se tiene que, si una funci´ on η ∈ H −1/2 (Γ) cumple que Vη ∈ H 1/2+s (Γ) para alg´ un s ∈ [0, 1/2], entonces η ∈ H −1/2+s (Γ) y kηk−1/2+s,Γ ≤ C (kηk−1/2,Γ + kVηk1/2+s,Γ ) .

28

Cap´ıtulo 2

Denotamos por V : H−1/2+s (Γ) → H1/2+s (Γ) el operador de simple capa actuando por componentes (s ∈ [−1/2, 1/2]). Las siguientes propiedades del operador Vτ := γ τ ◦ V quedan probadas en [19, Teorema 4.2]. Teorema 2.7.6. El operador lineal −1/2

Vτ : H×

1/2

(Γ) → H× (Γ)

es continuo y autoadjunto. Adem´ as, existe una constante α2 > 0 tal que hrotΓ η, Vτ (rotΓ η)iτ ,Γ ≥ α2 krotΓ ηk2

−1/2

H×

(Γ)

,

para cada η ∈ H 1/2 (Γ). Dado η ∈ H 1/2 (Γ), definimos el potencial de doble capa Z ∂Φ(x − y) e η(y) dSy ∀x ∈ Ω ∪ Ωe . Kη(x) := ∂n y Γ Para estudiar la regularidad de este operador, introducimos © ª H 1 (Ω)/C × W 1 (Ωe ) := u ∈ D0 (R3 ) ; u|Ω ∈ H 1 (Ω)/C , u|Ωe ∈ W 1 (Ωe ) , que es un espacio de Hilbert dotado de la norma del grafo H 1 (Ω)/C × W 1 (Ωe ) → R u

7→

³ ´1/2 |∇u|20,Ω + |∇u|20,Ωe .

Tambi´en consideramos ½ ¾ ∂u 1 1 e e Z := u ∈ H (Ω)/C × W (Ω ) ; ∆u = 0 en Ω ∪ Ω , [ ]Γ = 0 sobre Γ , ∂n que es un subespacio cerrado de H 1 (Ω)/C × W 1 (Ωe ) y, en particular, de Hilbert. Dada una funci´on u ∈ H 1 (Ω)/C × W 1 (Ωe ), definimos su salto y el salto de su derivada normal a trav´es de Γ como [u]Γ := γ(u|Ω ) − γ(u|Ωe )

y

[

∂u ]Γ := γn (∇u|Ω ) − γn (∇u|Ωe ) , ∂n

Herramientas matem´aticas

29

e respectivamente. Entonces, para cada η ∈ H 1/2 (Γ), el potencial de doble capa Kη 1 1 e pertenece al espacio H (Ω)/C × W (Ω ) y satisface las siguientes propiedades:  e  ∆(Kη) = 0 en D0 (Ω ∪ Ωe ) ,     e Γ = −η [Kη] sobre Γ , (2.7.3)   e  ∂(Kη)   [ ]Γ = 0 sobre Γ ; ∂n cf. [46, P´arrafo 3.2.2]. Lema 2.7.7. El operador lineal e : H 1/2 (Γ)/C → Z K es un isomorfismo. Demostraci´ on. Haciendo una formulaci´ on variacional de (2.7.3), deducimos que e Kη ∈ Z es soluci´on del problema    hallar u ∈ Z tal que (2.7.4) ∂v   (∇u, ∇v)0,Ω + (∇u, ∇v)0,Ωe = − h , ηi ∀v ∈ Z . ∂n 1/2,Γ N´otese que el Lema de Lax–Milgram garantiza que el problema (2.7.4) tiene una u ´nica soluci´on y que ´esta depende con continuidad de η ∈ H 1/2 (Γ). As´ı que el operador e : H 1/2 (Γ) → Z K est´a bien definido y es continuo. Adem´as, este operador es suprayectivo y su n´ ucleo es C, y por lo tanto concluimos el resultado. ¥ e : Z → L2 (R3 ) dado por Definimos el operador ∇  ∇(u| ) en Ω , Ω e := ∇u ∇(u| e ) en Ωe . Ω

30

Cap´ıtulo 2

Adem´as, denotamos por ∗ el producto de convoluci´ on de distribuciones e intro0 3 ducimos la distribuci´on δΓ ∈ D (R ) caracterizada por hδΓ , viD0 (R3 )×D(R3 ) := h1, v|Γ iD0 (Γ)×D(Γ)

∀v ∈ D(R3 ) .

Proposici´ on 2.7.8. Sea η ∈ H 1/2 (Γ). Entonces e Kη) e = − rot (Φ ∗ (rotΓ η δΓ )) ∇(

en D 0 (R3 ) .

Demostraci´ on. En primer lugar, estudiamos c´omo se relacionan las distribue ∈ D 0 (R3 ) y g := ∇( e Kη) e ∈ L2 (R3 ). Con este fin, integramos por ciones ∇(Kη) partes en e viD0 (R3 )×D(R3 ) = −(Kη, e div v)0,Ω − (Kη, e div v)0,Ωe h∇(Kη),

∀v ∈ D(R3 ) ,

y as´ı obtenemos que e viD0 (R3 )×D(R3 ) = (∇(Kη), e v)0,Ω∪Ωe − ([Kη] e Γ n, γv)0,Γ h∇(Kη),

∀v ∈ D(R3 ) .

Luego, utilizando la definici´on de g y la segunda propiedad de (2.7.3), resulta que e = g + η n δΓ ∇(Kη)

en D 0 (R3 ) ,

(2.7.5)

Por otro lado, la primera ecuaci´on de (2.7.3) significa que div g = 0 en Ω ∪ Ωe , mientras que la tercera se traduce en γn (g|Ω ) = γn (g|Ωe ) sobre Γ. Luego, en virtud del Corolario 2.3.3, tenemos que g ∈ H(div, R3 ) con div g = 0 en R3 . As´ı que, usando la tercera identidad de (2.0.1), rot rot g = −∆g en D 0 (R3 ) . Ahora, aplicando el operador rot rot en ambos miembros de (2.7.5), 0 = −∆g + rot rot(η n δΓ ) en D 0 (R3 ) . En consecuencia, g coincide con la siguiente convoluci´ on: g = −Φ ∗ rot rot(η n δΓ ) ,

(2.7.6)

Herramientas matem´aticas

31

que est´a bien definida porque la funci´on Φ es localmente integrable y la distribuci´on rot rot(η n δΓ ) tiene soporte compacto; v´ease [56, P´arrafo VI.3]. En lo que respecta a la distribuci´on rot(η n δΓ ), observamos que hrot(η n δΓ ), viD0 (R3 )×D(R3 ) = hη n δΓ , rot viD0 (R3 )×D(R3 ) = (η, γn (rot v))0,Γ , para cada v ∈ D(R3 ). Entonces, aplicando el Corolario 2.5.1 junto con la definici´on de divΓ , hrot(η n δΓ ), viD0 (R3 )×D(R3 ) = hη, divΓ (γτ v)i1/2,Γ = hrotΓ η, γτ viτ ,Γ . −1/2

Como rotΓ (H 1/2 (Γ)) ⊆ H×

(Γ), tenemos que

hrotΓ η, γτ viτ ,Γ = hrotΓ η, γvi1/2,Γ , y resulta que hrot(η n δΓ ), viD0 (R3 )×D(R3 ) = hrotΓ η, γvi1/2,Γ

∀v ∈ D(R3 ) .

Es decir, rot(η n δΓ ) = rotΓ η δΓ

en D 0 (R3 ) .

Utilizando esta propiedad junto con (2.7.6), concluimos que g = −Φ ∗ rot (rotΓ η δΓ ) = − rot (Φ ∗ (rotΓ η δΓ )) . ¥ A continuaci´on, utilizamos los resultados anteriores para estudiar la regularidad del operador hipersingular N y para dar una expresi´on expl´ıcita de su forma sesquilineal asociada; cf. [41, Lema 3.4]. Teorema 2.7.9. El operador lineal −1/2

N : H 1/2 (Γ)/C → H0

(Γ)

define un isomorfismo, que se denomina isomorfismo de Steklov–Poincar´e. Adem´ as, se relaciona con el operador Vτ a trav´es de la siguiente identidad: hN η, λi1/2,Γ = hrotΓ λ, Vτ rotΓ ηiτ ,Γ , para cualesquiera η, λ ∈ H 1/2 (Γ).

(2.7.7)

32

Cap´ıtulo 2

Demostraci´ on. En [46, Teorema 3.3.2] se prueba este resultado para fronteras regulares. En esta demostraci´on seguimos el razonamiento expuesto en [46] y lo generalizamos a fronteras Lipschitz continuas. Es conocido que la aplicaci´on lineal −1/2

Z → H0 ∂v v 7→ ∂n

(Γ)

define un isomorfismo; v´ease [46, P´arrafo 2.5.7]. Entonces, aplicando el Lema 2.7.7, deducimos que el operador N : H 1/2 (Γ)/C →

Z

−1/2

→ H0

(Γ)

e e e 7→ N η := − ∂(Kη|Ω ) = − ∂(Kη|Ωe ) η 7→ Kη ∂n ∂n es un isomorfismo. Para demostrar la identidad (2.7.7), utilizamos la tercera ecuaci´on de (2.7.3) e integramos por partes como sigue: hN η, λi1/2,Γ = h

e ∂(Kη) e Γ i1/2,Γ = (∇(Kη), e ∇(Kλ)) e 0,Ω + (∇(Kη), e ∇(Kλ)) e 0,Ωe . , [Kλ] ∂n

N´otese que, en virtud del Lema 2.7.8, e Kη) e = − rot(Φ ∗ (rotΓ η δΓ )) en D 0 (R3 ) . ∇( e Ω ∈ H(rot, Ω), podemos integrar por Como (Φ ∗ (rotΓ η δΓ ))|Ω ∈ H1 (Ω) y ∇(Kλ)| partes con la f´ormula (2.3.3): e ∇(Kλ)) e 0,Ω = −hγτ (∇(Kλ)| e Ω ), γτ (Φ ∗ (rotΓ η δΓ )|Ω )i . (∇(Kη), τ ,Γ An´alogamente, tenemos que e ∇(Kλ)) e 0,Ωe = hγτ (∇( e Kλ)| e Ωe ), γτ (Φ ∗ (rotΓ η δΓ )|Ωe )i . (∇(Kη), τ ,Γ

Herramientas matem´aticas Luego

33

h i e Kλ)) e hN η, λi1/2,Γ = −h γτ (∇( , γτ (Φ ∗ (rotΓ η δΓ ))i Γ

τ ,Γ

.

e y rotΓ junto con la tercera ecuaci´on de Pero, utilizando las definiciones de ∇ (2.7.3), resulta que h i h i e Kλ)) e e γτ (∇( = rotΓ Kλ = − rotΓ λ ; Γ

Γ

adem´as, aplicando las definiciones del producto de convoluci´ on y del operador V, Φ ∗ (rotΓ η δΓ )|Γ = V(rotΓ η) en H1/2 (Γ) . Por lo tanto, concluimos que hN η, λi1/2,Γ = hrotΓ λ, Vτ rotΓ ηiτ ,Γ . ¥ Para cada s ∈ [0, 1], introducimos el operador identidad I : H s (Γ) → H s (Γ) y denotamos por K∗ : H −s (Γ) → H −s (Γ) el operador dual de K: hK∗ λ, ηis,Γ = hλ, Kηis,Γ

∀λ ∈ H −s (Γ), η ∈ H s (Γ) .

Combinando la f´ormula de representaci´ on integral (2.7.2) con las relaciones de salto de los operadores de simple y doble capa, deducimos el siguiente resultado; cf. [46, Teorema 3.1.2]. Teorema 2.7.10. Sea ϕ ∈ W 1 (Ωe ) una funci´ on arm´ onica. Entonces 1 ∂ϕ ( I − K)(γϕ) + V( ) = 0 2 ∂n

y

1 ∂ϕ ( I + K∗ )( ) + N (γϕ) = 0 2 ∂n

sobre Γ . (2.7.8)

34

Cap´ıtulo 2

2.8.

Espacios b´ asicos para problemas de evoluci´ on

En este apartado introducimos algunos espacios funcionales auxiliares que utilizaremos para estudiar el problema de evoluci´ on (1.1.1). Para una exposici´on m´as detallada, puede consultarse [56, P´arrafo V.5]. Consideramos funciones definidas en un intervalo acotado (0, T ) y con valores en un espacio V que sea de Hilbert y separable dotado de la norma k·kV . En lo sucesivo, usamos indistintamente las notaciones df = ∂t f dt para representar la derivaci´on con respecto a la variable temporal t ∈ (0, T ). Sea C 0 ([0, T ]; V ) el espacio de las funciones continuas f : [0, T ] → V . Este espacio es de Banach con la norma kf kC 0 ([0,T ];V ) := m´ax kf (t)kV . t∈[0,T ]

En general, para cada m ∈ N, el espacio ½ ¾ dk f m 0 0 C ([0, T ]; V ) := f ∈ C ([0, T ]; V ) ; ∈ C ([0, T ]; V ) ∀1 ≤ k ≤ m dtk es de Banach con la norma kf kC m ([0,T ];V ) := m´ax k 1≤k≤m

dk f k 0 . dtk C ([0,T ];V )

Definimos el espacio L2 ((0, T ); V ) de las funciones f : (0, T ) → V que son medibles en sentido B¨ochner y tales que Z T kf k2L2 ((0,T );V ) := kf (t)k2V dt < +∞ . 0

Tambi´en introducimos el espacio ½ ¾ Z t 1 2 H ((0, T ); V ) := f (t) = f0 + f1 (s) ds ; f0 ∈ V, f1 ∈ L ((0, T ); V ) . 0

El siguiente resultado queda demostrado en [56, Corolario V.5.2].

Herramientas matem´aticas

35

Lema 2.8.1. Sean V y H dos espacios de Hilbert separables. Si A : V → H es una aplicaci´ on lineal continua, entonces Z Z t ¡ t ¢ f (s) ds , Af (s) ds = A 0

0

para cada f ∈ L2 ((0, T ); V ) y t ∈ (0, T ]. En particular, si denotamos por V 0 el espacio dual de V con respecto al producto h·, ·iV 0 ×V , entonces Z 0

t

hf (s), viV

0 ×V

Z t ds = h f (s) ds, viV 0 ×V , 0

para cualesquiera f ∈ L2 ((0, T ); V 0 ), v ∈ V y t ∈ (0, T ].

Cap´ıtulo 3

Eddy currents en conductores simplemente conexos En este cap´ıtulo estudiamos el problema de eddy currents en r´egimen arm´onico (1.2.1–1.2.2) cuando la topolog´ıa del conductor es sencilla. M´as concretamente, suponemos que la regi´on conductora Ωc es conexa y simplemente conexa, y que su frontera Σ := ∂Ωc es Lipschitz continua, conexa y simplemente conexa. Tambi´en asumimos que la funci´on dato j cumple las condiciones de compatibilidad (1.2.3) y que su soporte no toca la superficie Σ. Entonces, razonando por superposici´on, podemos suponer que el soporte de j est´a enteramente contenido en el conductor Ωc ´o en el diel´ectrico Ωec := R3 \ Ωc . De hecho, aqu´ı supondremos que el soporte de j est´a enteramente contenido en Ωc y nos remitimos al art´ıculo [41] para el caso contrario.

3.1.

Formulaci´ on variacional

En primer lugar, buscamos una formulaci´ on variacional BEM–FEM del problema (1.2.1–1.2.2) en t´erminos del campo h. Para justificar los razonamientos de este apartado, recordamos la siguiente propiedad referente al comportamiento 37

38

Cap´ıtulo 3

asint´otico del campo electromagn´etico. Proposici´ on 3.1.1. Si e y h son soluciones del problema (1.2.1–1.2.2), entonces h(x) = O(

1 ) |x|2

y

e(x) = O(

uniformemente en las direcciones

x |x|

1 ) |x|2

cuando |x| → ∞ ,

con x ∈ R3 \ {0}.

Demostraci´ on. V´ease [6, Proposici´on 3.1].

¥

A continuaci´on, seguimos la estrategia de [14] e introducimos un potencial escalar magn´etico en Ωec . Dado que soporte(j) ⊂ Ωc , la primera ecuaci´on de (1.2.1) junto con la hip´otesis (1.1.6) implican que 1 (rot h − j) σ rot h = 0 e =

en Ωc , en Ωec .

(3.1.1)

Adem´as, aplicando el operador div en la segunda ecuaci´on de (1.2.1) y teniendo en cuenta la hip´otesis (1.1.5), div h = 0 en Ωec .

(3.1.2)

Las propiedades anteriores sugieren tomar una funci´on test q : R3 → C3 suficientemente regular, con rot q = 0 y div q = 0 en Ωec . Utilizando la segunda ecuaci´on de (1.2.1) e integrando por partes, resulta que ıω (µh, q)0,R3 + (e, rot q)0,Ωc = 0 . Entonces eliminamos el campo e de esta ecuaci´on variacional usando la primera identidad de (3.1.1): 1 1 ıω (µh, q)0,R3 + ( rot h, rot q)0,Ωc = ( j, rot q)0,Ωc . σ σ

(3.1.3)

Eddy currents en conductores simplemente conexos

3.1.1.

39

Un potencial escalar magn´ etico

A continuaci´on representamos el campo h en el diel´ectrico Ωec mediante un potencial escalar arm´onico. Con este fin, recordamos que la funci´on h|Ωec satisface las propiedades (3.1.1–3.1.2) y, por lo tanto, est´a bajo las hip´otesis del Lema 2.7.2. As´ı que existe un u ´nico potencial arm´onico ψ ∈ W 1 (Ωec ) tal que ∇ψ = h en Ωec . N´otese que, en virtud de la primera ecuaci´on de (1.2.1), tenemos que h ∈ H(rot, R3 ) y por lo tanto su componente tangencial es continua a trav´es de Σ: rotΣ ψ = γτ (h|Ωc ) sobre Σ . Adem´as, tomando el operador div en la segunda ecuaci´on de (1.2.1), resulta que div(µh) = 0 en R3 . En particular, deducimos que µh ∈ H(div, R3 ) y por lo tanto su componente normal es continua a trav´es de Σ: µ0

∂ψ = γn ((µh)|Ωc ) ∂n

sobre Σ .

An´alogamente, para cada funci´on test q : R3 → C3 que sea suficientemente regular y que tenga div q = 0 y rot q = 0 en Ωec , existe una u ´nica funci´on 1 e arm´onica ϕ ∈ W (Ωc ) con ∇ϕ = q en Ωec . Adem´as, esta funci´on cumple que rotΣ ϕ = γτ (q|Ωc )

sobre Σ .

Nota 3.1.1 Es importante subrayar que podemos aplicar el Lema 2.7.2 gracias a las hip´otesis que estamos haciendo sobre la geometr´ıa del conductor Ωc . En la pr´actica, la geometr´ıa del conductor Ωc puede ser complicada y es habitual on de una funci´on que Ωec no sea simplemente conexo. En este caso, la representaci´ e e u ∈ H(rot, Ωc ) con rot u = 0 en Ωc involucra los denominados campos vectoriales de Neumann arm´ onicos. En consecuencia, si generaliz´asemos directamente

40

Cap´ıtulo 3

la formulaci´on variacional que proponemos en este cap´ıtulo para conductores que no sean simplemente conexos, necesitar´ıamos introducir dichos campos de Neumann; cf. [41]. Desde el punto de vista num´erico, la introducci´on de tales campos de Neumann complica la implementaci´ on y aumenta el coste computacional. Es precisamente esta dificultad la que motiva la estrategia que seguiremos en el Cap´ıtulo 4. Ahora, integramos por partes en Ωec y as´ı reescribimos (3.1.3) como ıω(µh, q)0,Ωc + c(h, q) − ıωµ0 h

∂ψ , ϕi1/2,Σ = l(q) , ∂n

(3.1.4)

donde 1 c(q1 , q2 ) := ( rot q1 , rot q2 )0,Ωc σ

1 y l(q) := ( j, rot q)0,Ωc . σ

La ecuaci´on (3.1.4) sugiere introducir la inc´ognita auxiliar λ := Entonces el potencial ψ ∈ W 1 (Ωec ) es la soluci´on d´ebil del problema  en Ωec ,  ∆ψ = 0

∂ψ sobre Σ. ∂n

 ∂ψ = λ sobre Σ . ∂n Adem´as, aplicando los Teoremas 2.7.4 y 2.7.10, tenemos la representaci´ on integral Z Z ∂Φ(x − y) ψ(x) = γψ(y) dSy − Φ(x − y) λ(y) dSy ∀x ∈ Ωec , (3.1.5) ∂n y Σ Σ junto con las ecuaciones integrales 1 ( I − K)γψ + Vλ = 0 2

3.1.2.

1 y ( I + K∗ )λ + N γψ = 0 2

sobre Σ .

(3.1.6)

Una formulaci´ on variacional BEM–FEM

Nuestro siguiente objetivo es combinar la ecuaci´on (3.1.4) con las ecuaciones integrales (3.1.6). Para ello, hacemos una formulaci´ on variacional de la segunda ecuaci´on de (3.1.6): 1 hλ, γϕi1/2,Σ = h( I − K∗ )λ, γϕi1/2,Σ − hN γψ, γϕi1/2,Σ 2

∀ϕ ∈ W 1 (Ωec ) .

Eddy currents en conductores simplemente conexos

41

Entonces, usando la identidad (2.7.7) cuando rotΣ ψ = γτ h y rotΣ ϕ = γτ q, 1 hλ, γϕi1/2,Σ = hλ, ( I − K)γϕi1/2,Σ − hγτ q, Vτ (γτ h)iτ ,Σ . 2 Esta propiedad nos permite reescribir la ecuaci´on (3.1.4) como ıω(µh, q)0,Ωc + c(h, q) + ıωµ0 hγτ q, Vτ (γτ h)iτ ,Σ 1 − ıωµ0 hλ, ( I − K)γϕi1/2,Σ = l(q) . 2

(3.1.7) −1/2

En consecuencia, buscamos dos funciones h ∈ H(rot, Ωc ) y λ ∈ H0 junto con un potencial arm´onico ψ ∈ W 1 (Ωec ), de forma que 1 ( I − K)γψ + Vλ = 0 , 2 rotΣ ψ = γτ h ,

(Σ),

(3.1.8)

y que adem´as cumplan la ecuaci´on variacional (3.1.7) para cualquier q ∈ H(rot, Ωc ) y el correspondiente potencial arm´onico ϕ ∈ W 1 (Ωec ) con rotΣ ϕ = γ τ q. N´otese que, seg´ un la Proposici´on 2.5.1 y el Lema 2.6.4, γn (rot h) = divΣ (γτ h) = divΣ (rotΣ ψ) = 0 ; es decir h ∈ X(Ωc ), donde X(Ωc ) := {u ∈ H(rot, Ωc ) ; γn (rot u) = 0 sobre Σ} . An´alogamente, puede comprobarse que las funciones test q utilizadas en (3.1.8) pertenecen al espacio X(Ωc ). Proposici´ on 3.1.2. El espacio X(Ωc ) es un subespacio cerrado de H(rot, Ωc ). Demostraci´ on. En virtud de la Proposici´on 2.3.2, el operador H(rot, Ωc ) → H(div, Ωc ) → H −1/2 (Σ) u 7→ rot u 7→ γn (rot u) ,

(3.1.9)

es continuo. Luego su n´ ucleo, X(Ωc ), es un subespacio cerrado de H(rot, Ωc ). ¥

42

Cap´ıtulo 3

1 Por otro lado, es conocido que K1 = y entonces el Teorema 2.7.5 garantiza 2 que el siguiente operador est´a bien definido y es continuo: 1 I − K : H 1/2 (Σ)/C → H 1/2 (Σ)/C . 2 Adem´as, en virtud de la Proposici´on 2.5.1 y del Teorema 2.5.2, −1/2

γτ X(Ωc ) = ker(divΣ ) ∩ H×

(divΣ , Σ) .

Entonces, usando el Corolario 2.6.6, deducimos que el operador lineal 1 1/2 B := ( I − K) ◦ rot−1 (Σ)/C Σ ◦γτ : X(Ωc ) → H 2 est´a bien definido y es acotado. En efecto, existe una constante C > 0 tal que kBqk1/2,Σ ≤ C kγτ qkH−1/2 (Σ) ×

∀q ∈ X(Ωc ) .

(3.1.10)

Utilizando este operador, podemos eliminar los potenciales ψ y ϕ de nuestro problema. M´as concretamente, tenemos que 1 1 Bh = ( I − K)γψ y Bq = ( I − K)γϕ , 2 2 en virtud de la segunda ecuaci´on de (3.1.8) y la condici´on rotΣ ϕ = γ τ q sobre la funci´on test. Entonces, el problema (3.1.7–3.1.8) consiste en buscar dos funciones −1/2 h ∈ X(Ωc ) y λ ∈ H0 (Σ) tales que ıω(µh, q)0,Ωc + c(h, q) + ıωµ0 hγτ q, Vτ (γτ h)iτ ,Σ − ıωµ0 hλ, Bqi1/2,Σ = l(q) , ωµ0 hη, Bhi1/2,Σ + ωµ0 hη, Vλi1/2,Σ = 0 , −1/2

para cualesquiera q ∈ X(Ωc ) y η ∈ H0

(Σ).

b c ) := X(Ωc ) × H −1/2 (Γ) y definimos las formas Introducimos el espacio X(Ω 0 sesquilineales a0 (q1 , q2 ) := ıω(µq1 , q2 )0,Ωc + ıωµ0 hγτ q2 , Vτ (γτ q1 i)τ ,Σ , ³ ´ b2 ) := a0 (q1 , q2 ) + ıωµ0 hη2 , Vη1 i1/2,Σ + hη2 , Bq1 i1/2,Σ − hη1 , Bq2 i1/2,Σ , a(b q1 , q b c ). Con esta notaci´on, b1 = (q1 , η1 ), q b2 = (q2 , η2 ) ∈ X(Ω para cualesquiera q reescribimos el problema variacional anterior:

Eddy currents en conductores simplemente conexos

43

b ∈ X(Ω b c ) tal que Problema 3.1.1 Hallar h b q b c) . b) + c(h, q) = l(q) ∀b a(h, q ∈ X(Ω Teorema 3.1.3. El Problema 3.1.1 tiene soluci´ on y ´esta es u ´nica. b c ), aplicando el Teorema 2.7.5, b = (q, η) ∈ X(Ω Demostraci´ on. Para cada q ¸ · µ 1−ı b) = ( q, q)0,Ωc + hγτ q, Vτ (γτ q)iτ ,Σ + hη, Vηi1/2,Σ , a(b q, q Re ωµ0 µ0 y resulta que

donde

·

¸ 1−ı b) + c(q, q)) ≥ α kb Re (a(b q, q qk2X(Ω b c), ωµ0 ½ α := m´ın 1,

1 , α1 ωµ0 σ1

¾ > 0.

Teniendo en cuenta esta propiedad, podemos reescalar el Problema 3.1.1 de modo que est´e bajo las hip´otesis del Lema de Lax–Milgram y as´ı concluimos el resultado. ¥

3.2.

Esquema discreto

En este apartado proponemos una versi´ on discreta del Problema 3.1.1 y demostramos que ´esta tiene una u ´nica soluci´on. Tambi´en hacemos un an´alisis a priori del orden del error cometido al aproximar la soluci´on exacta por la soluci´on del problema discreto. En lo que sigue, suponemos que el dominio Ωc es poli´edrico y denotamos por {Th (Ωc ) ; h > 0} una familia regular de mallas tetra´edricas de Ωc . Adem´as, definimos Th (Σ) como la triangulaci´on que induce Th (Ωc ) sobre la superficie Σ, Th (Σ) := {F cara de alg´ un T ∈ Th (Ωc ) ; F ⊂ Σ} .

44

Cap´ıtulo 3

3.2.1.

Elementos finitos de N´ ed´ elec de primer orden

Aproximamos el espacio H(rot, Ωc ) utilizando elementos finitos de N´ed´elec de primer orden, tambi´en conocidos como elementos de arista. En la literatura, estos elementos aparecen descritos de forma independiente por varios autores, siendo H. Whitney el primero del que tenemos constancia; cf. [55]. En lo que respecta a los elementos de N´edelec de orden mayor que uno, fueron introducidos por J.-C. N´ed´elec en [47] y de ah´ı su denominaci´on. Fijado un tetraedro T ∈ Th (Ωc ) y dado m ∈ N, sea Pm (T ) el espacio de polinomios definidos en T y de grado no mayor que m. Consideramos © ª N D(T ) := a ∧ x + b ; a, b ∈ C3 , que es un subespacio de (P1 (T ))3 con dimensi´on 6. Para cada arista E de T , denotamos por τ TE un vector unitario tangente a E y definimos Z mTE (p) :=

E

p(x) · τ TE dLx ,

siempre que p : T → C3 sea una funci´on suficientemente regular. © ª Introducimos el conjunto GL(T ) := mTE ; E arista de T . Entonces la terna (T, N D(T ), GL(T )) define un elemento finito, que se dice elemento de N´ed´elec de primer orden ´o elemento de arista; cf. [45, Lema 5.36]. Consideramos el espacio de elementos finitos global asociado: N D h (Ωc ) = {vh ∈ H(rot, Ωc ) ; vh |T ∈ N D(T ) ∀T ∈ Th (Ωc )} .

(3.2.1)

Para cada T ∈ Th (Ωc ), denotamos por I T el operador de interpolaci´ on local 3 asociado al espacio N D(T ). Es decir, dada una funci´on u : T → C suficientemente regular, el elemento I T u ∈ N D(T ) se caracteriza por mTE (I T u) = mTE (u) ∀E arista de T . Ahora podemos definir el operador de interpolaci´ on global asociado al espacio N D h (Ωc ). M´as concretamente, dada una funci´on u : Ωc → C3 que sea

Eddy currents en conductores simplemente conexos

45

suficientemente regular, definimos I h u ∈ N D h (Ωc ) como (I h u)|T = I T (u|T ) ∀T ∈ Th (Ωc ) .

(3.2.2)

En efecto, el siguiente resultado garantiza que esta definici´on es adecuada; cf. [2, Lema 5.1] ´o bien [45, Lema 5.38]. Lema 3.2.1. Para cada s > 1/2, el operador lineal I h : Hs (rot, Ωc ) → N D h (Ωc ) est´ a bien definido y es acotado. En el siguiente resultado recogemos algunas de las propiedades de aproximaci´on del espacio N D h (Ωc ) demostradas en [2, Proposici´on 5.6] y en [45, Teorema 5.41]. Teorema 3.2.2. Para cada s ∈ (1/2, 1], existe una constante C > 0 con kq − I h qkrot,Ωc ≤ C hs kqkHs (rot,Ωc ) , para cualquier funci´ on q ∈ Hs (rot, Ωc ) y todo h > 0 suficientemente peque˜ no. Sean T un tetraedro de Th (Ωc ) y F una cara de T . Definimos Z |F | := 1 dSx , F

que se dice medida de F . Denotamos por nF el vector unitario normal a F con orientaci´on hacia el exterior de T . Adem´as, en casi todo punto de ∂F definimos un vector unitario n∂F que est´a sobre el plano que contiene a F y que es normal a ∂F con orientaci´on hacia el exterior de F . V´ease la Figura 3.1 a continuaci´ on. Lema 3.2.3. Sea u una funci´ on de H2 (T ). Entonces Z X γn (rot u) dSx = (nF ∧ n∂F )|E · τ TE mTE (u) . F

E lado de F

(3.2.3)

46

Cap´ıtulo 3

nF F

n∂F

T

Figura 3.1: Vectores nF y n∂F Demostraci´ on. Consideramos una funci´on u ∈ D(T ). Aplicando sucesivamente la Proposici´on 2.5.1 en el dominio T ⊂ R3 y la f´ormula de Stokes sobre la superficie F , resulta que Z Z Z γn (rot u) dSx = div∂T (γτ u) dSx = (γτ u) · n∂F dLx . F

F

∂F

Ahora bien, seg´ un la definici´on (2.2.4) y las identidades vectoriales (2.0.1), (γτ u)|E · n∂F = (nF ∧ n∂F ) · u|E

y

(nF ∧ n∂F )|E = (nF ∧ n∂F )|E · τ TE τ TE ,

para cada lado E de F . En consecuencia, la ecuaci´on (3.2.3) se cumple para cada funci´on u ∈ D(T ). Finalmente, el Teorema de trazas 2.2.1 nos permite concluir con un razonamiento por densidad. ¥ Lema 3.2.4. Sea p una funci´ on deN D(T ) tal que mTE (p) = 0 para toda arista E de F . Entonces γn (rot p) = 0 sobre F . Demostraci´ on. Consideramos una funci´on p ∈ N D(T ) bajo las hip´otesis del enunciado. Entonces, por un lado, el Lema 3.2.3 garantiza que Z X γn (rot p) dSx = (nF ∧ n∂F )|E · τ TE mTE (p) = 0 . F

E lado de F

Eddy currents en conductores simplemente conexos

47

Por otro lado, tenemos que rot p ∈ C3 es constante en T y, en particular, γn (rot p) ∈ C es constante sobre F . As´ı que Z 1 γn (rot p) = γn (rot p) dSx = 0 sobre F . |F | F ¥

3.2.2.

Descripci´ on y an´ alisis del problema discreto

Para resolver num´ericamente el Problema 3.1.1, aproximamos las funciones del espacio X(Ωc ) con las del espacio de dimensi´on finita Xh (Ωc ) := N D h (Ωc ) ∩ X(Ωc ) .

(3.2.4)

Proposici´ on 3.2.5. Sea u ∈ Hs (rot, Ωc ) ∩ X(Ωc ), con s ∈ (1/2, 1]. Entonces I h u ∈ Xh (Ωc ). Demostraci´ on. Por definici´on, I h u ∈ N D h (Ωc ). As´ı que lo que estudiaremos es si γn (rot(I h u)) se anula sobre Σ. Con este fin, introducimos la funci´on auxiliar uh := u − I h u y consideramos una cara F ∈ Th (Σ) de un tetraedro T ∈ Th (Ωc ). N´otese que mTE (uh |T ) = mTE (u|T ) − mTE (I T (u|T )) = 0 ∀E arista de F . Aplicando el Lema 3.2.3, deducimos que Z γn (rot uh ) dSx = 0 . F

Entonces, utilizando la propia definici´on de uh y la hip´otesis γn (rot u) = 0, Z Z γn (rot(I h u)) dSx = γn (rot u) dSx = 0 . F

F

Adem´as, rot(I h u) es constante en T y, en particular, γn (rot(I h u)) es constante sobre F . Por lo tanto, concluimos que Z 1 γn (rot(I h u)) dSx = 0 sobre F . γn (rot(I h u)) = |F | F ¥

48

Cap´ıtulo 3 −1/2

Por otro lado, aproximamos las funciones de H0 (Σ) con funciones que son constantes sobre cada cara de Th (Σ). Es decir, utilizamos el espacio de dimensi´on finita © ª −1/2 Λh (Σ) := η ∈ L2 (Σ) ; η|F ∈ C ∀F ∈ Th (Σ) ∩ H0 (Σ) . (3.2.5) b h (Ωc ) := Xh (Ωc ) × Λh (Σ), proponemos el siguiente m´etodo de Denotando X Galerkin para aproximar la soluci´on del Problema 3.1.1. b h = (hh , λh ) ∈ X b h (Ωc ) tal que Problema 3.2.1 Hallar h bh, q b h (Ωc ) . b) + c(hh , q) = l(q) ∀b a(h q∈X Proposici´ on 3.2.6. El problema discreto 3.2.1 est´ a bien planteado. Demostraci´ on. Basta recordar que la forma sesquilineal b c ) × X(Ω b c) → C X(Ω b2 ) 7→ (b q1 , q

1−ı b2 ) + c(q1 , q2 )) (a(b q1 , q ωµ0

es el´ıptica; v´ease la demostraci´on del Teorema 3.1.3.

¥

Como consecuencia directa del resultado anterior, se cumple una estimaci´on tipo Lema de C´ea: b−h bhk b b−q bkX(Ω kh ≤ C kh b c) X(Ωc )

b h (Ωc ) , ∀b q∈X

(3.2.6)

siendo C > 0 una constante independiente del par´ametro de discretizaci´on h. Nota 3.2.1 Es importante se˜ nalar que la inversa del operador tangencial rotΓ est´a involucrada en la definici´on de la forma sesquilineal a(·, ·) y, en consecuencia, nuestro m´etodo num´erico parece dif´ıcil de implementar y computacionalmente costoso. Sin embargo, esto no sucede en la pr´actica porque dicho operador puede eliminarse en el c´alculo efectivo de la matriz del sistema lineal asociado al Problema 3.2.1. Al respecto, v´ease el P´arrafo 3.3.

Eddy currents en conductores simplemente conexos

3.2.3.

49

An´ alisis de la convergencia

En este apartado analizamos el error cometido al aproximar la soluci´on del Problema 3.1.1 con la del Problema 3.2.1. Con este fin, primero estudiamos los errores de aproximaci´on asociados a Xh (Ωc ) y Λh (Σ) como subespacios de X(Ωc ) y H −1/2 (Σ), respectivamente. Por un lado, fijado s ∈ (1/2, 1], la Proposici´on 3.2.5 permite acotar el error correspondiente al espacio Xh (Ωc ) con la estimaci´on del error de interpolaci´on: ´ınf

ku − qkrot,Ωc ≤ ku − I h ukrot,Ωc

q∈Xh (Ωc )

∀u ∈ Hs (rot, Ωc ) ∩ X(Ωc ) .

Entonces, utilizando el Teorema 3.2.2, ´ınf

q∈Xh (Ωc )

ku − qkrot,Ωc ≤ C hs kukHs (rot,Ωc )

∀u ∈ Hs (rot, Ωc ) ∩ X(Ωc ) . (3.2.7)

Por otro lado, para estudiar el error asociado al espacio Λh (Σ) utilizamos el siguiente resultado cl´asico. Proposici´ on 3.2.7. Para cada s ∈ (1/2, 1], existe una constante C > 0 con ´ınf

kϕ − ηk−1/2,Σ ≤ C hs kϕk−1/2+s,Σ ,

η∈Λh (Σ)

(3.2.8)

para cualquier ϕ ∈ H −1/2+s (Σ) y siempre que h > 0 sea suficientemente peque˜ no. Ahora, aplicando la desigualdad tipo Lema de C´ea (3.2.6) junto con las estimaciones (3.2.7) y (3.2.8), obtenemos la siguiente estimaci´ on del error asociado al esquema discreto. b = (h, λ) del Problema 3.1.1 perTeorema 3.2.8. Supongamos que la soluci´ on h tenece al espacio Hs (rot, Ωc ) × H −1/2+s (Σ) para alg´ un s ∈ (1/2, 1]. Entonces tenemos la siguiente estimaci´ on de error: ¡ ¢ kh − hh krot,Ωc + kλ − λh k−1/2,Σ ≤ C hs khkHs (rot,Ωc ) + kλk−1/2+s,Σ .

50

Cap´ıtulo 3

3.3.

Forma matricial del esquema discreto

En este apartado reescribimos el Problema 3.2.1 como un sistema lineal. b h (Ωc ) y a Con este fin, primero determinamos una base expl´ıcita del espacio X continuaci´on reformulamos el problema en forma matricial.

3.3.1.

b h (Ωc ) Base expl´ıcita de X

Sea Ah (Ωc ) el conjunto de las aristas de la malla Th (Ωc ) y denotamos por el conjunto de las aristas interiores, es decir,

A0h (Ωc )

A0h (Ωc ) := {E ∈ Ah (Ωc ) ; E ∩ Ωc 6= ∅} . Tambi´en introducimos el conjunto Vh (Ωc ) de v´ertices de la malla Th (Ωc ), as´ı como Vh (Σ) := {v ∈ Vh (Ωc ) ; v ∈ Σ} . Para cada E ∈ Ah (Ωc ), sea qE el elemento de N D h (Ωc ) caracterizado por mE 0 (qE ) = δE,E 0

∀E 0 ∈ Ah (Ωc ) .

An´alogamente, para cada v ∈ Vh (Ωc ), sea ϕv la funci´on de C(Ωc ) que cumple ϕv |T ∈ P1 (T ) ∀T ∈ Th (Ωc )

y

ϕv (v 0 ) = δv,v0 ∀v 0 ∈ Vh (Ωc ) .

Proposici´ on 3.3.1. Fijado v0 ∈ Vh (Σ), el conjunto Bh (Ωc ) :=

© ª qE ; E ∈ A0h (Ωc ) ∪ {∇ϕv ; v ∈ Vh (Σ) \ {v0 }}

define una base de Xh (Ωc ). Demostraci´ on. En primer lugar, estudiamos si Bh (Ωc ) ⊂ Xh (Ωc ) . Por un lado, para cada E ∈ A0h (Ωc ) sabemos que qE ∈ N D h (Ωc ) y que mE 0 (qE ) = 0 para toda arista E 0 ⊂ Σ. Entonces el Lema 3.2.4 garantiza que γn (rot qE ) = 0 sobre cada F ∈ Fh (Σ) y por tanto qE ∈ Xh (Ωc ).

Eddy currents en conductores simplemente conexos

51

Por otro lado, para cualquier v ∈ Vh (Σ), la funci´on ∇ϕv ∈ H(rot, Ωc ) es constante en cada T ∈ Th (Ωc ). Luego ∇ϕ ∈ N D h (Ωc ) y, como rot(∇ϕv ) = 0 en Ωc , resulta que ∇ϕv ∈ Xh (Ωc ). En segundo lugar, estudiamos si Bh (Ωc ) es linealmente independiente. Para ello, fijamos aE , av ∈ C (E ∈ A0h (Ωc ), v ∈ Vh (Σ) \ {v0 }) tales que X X aE qE + av ∇ϕv = 0 en Ωc . (3.3.1) E∈A0h (Ωc )

v∈Vh (Σ)\{v0 }

Por la propia definici´on del operador rotΣ , sabemos que γτ (∇ϕv ) = rotΣ ϕv para cada v ∈ Vh (Σ). Adem´as, dada E ∈ A0h (Ωc ), tenemos que mE 0 (qE ) = 0 para cada arista E 0 ⊂ Σ y, en consecuencia, γτ qE = 0. Luego, tomando el operador γτ en ambos miembros de (3.3.1), resulta que X av rotΣ ϕv = 0 sobre Σ . v∈Vh (Σ)\{v0 }

Esta propiedad junto con el Teorema 2.6.3 garantizan que

av ϕv (v0 ) = 0, deducimos que

v∈Vh (Σ)\{v0 }

X

av ϕv = 0

sobre Σ .

v∈Vh (Σ)\{v0 }

En particular, evaluando en cada v ∈ Vh (Σ) \ {v0 }, X av = av0 ϕv0 (v) = 0 ∀v ∈ Vh (Σ) \ {v0 } . v 0 ∈Vh (Σ)\{v0 }

Ahora la identidad (3.3.1) se reescribe como X aE qE = 0 en Ωc . E∈A0h (Ωc )

En particular, tomando el operador mE para cada E ∈ A0h (Ωc ), X aE = aE 0 mE (qE 0 ) = 0 ∀E ∈ A0h (Ωc ) . E 0 ∈A0h (Ωc )

av ϕv es

v∈Vh (Σ)\{v0 }

X

constante sobre Σ. Como

X

52

Cap´ıtulo 3

Finalmente, estudiamos si Bh (Ωc ) es un sistema generador de Xh (Ωc ). Con este objetivo, consideramos q ∈ Xh (Ωc ) arbitrario. Seg´ un la Proposici´on 2.5.1, divΣ (γτ q) = γn (rot q) = 0

sobre Σ .

Entonces el Corolario 2.6.6 garantiza que existe un potencial ψ ∈ H 1/2 (Σ) con rotΣ ψ = γτ q , y este potencial es u ´nico salvo constante aditiva. N´otese que (rotΣ ψ)|T ∩Σ = (γτ q)|T ∩Σ ∈ C3 (constante)

∀T ∈ Th (Ωc ) .

Luego la funci´on ψ ∈ H 1/2 (Σ) es polinomial a trozos, es decir ψ|F ∈ P1 (F ) ∀F ∈ Th (Σ) , y en consecuencia ψ ∈ C(Σ). Introducimos la funci´on ψ0 := ψ − ψ(v0 ) ∈ C(Σ), que es polinomial a trozos y cumple que ψ0 (v0 ) = 0 y rotΣ ψ0 = γτ q. Definimos X ψ0 (v) ϕv en Ωc , ψe0 := v∈Vh (Σ)\{v0 }

que es un levantamiento de ψ0 polinomial a trozos: ψe0 |Σ = ψ0

sobre Σ

y

ψe0 |T ∈ P1 (T ) ∀T ∈ Th (Ωc ) .

Entonces la funci´on q − ∇ψe0 pertenece al espacio N D h (Ωc ) y satisface γτ (q − ∇ψe0 ) = γτ q − rotΣ (ψe0 |Σ ) = γτ q − rotΣ ψ0 = 0 sobre Σ . © ª Como el espacio N D h (Ωc ) ∩ ker γτ est´ a generado por qE ; E ∈ A0h (Ωc ) , concluimos que X X q= ψ0 (v) ∇ϕv + mE (q − ∇ψe0 ) qE . v∈Vh (Σ)\{v0 }

E∈A0h (Ωc )

¥

Eddy currents en conductores simplemente conexos

53

Para cada F ∈ Th (Σ), denotamos por χF su funci´on caracter´ıstica. Entonces, fijada una cara F0 ∈ Th (Σ), el conjunto ½ ¾ 1 1 ρF := χF − χF ; F ∈ Th (Σ) \ {F0 } |F | |F0 | 0 constituye una base de Λh (Σ).

3.3.2.

Sistema de ecuaciones asociado al problema discreto

Expresamos las inc´ognitas del Problema 3.2.1 en t´erminos de las bases de los espacios Xh (Ωc ) y Λh (Σ) propuestas en el apartado anterior: X X X hh = hE qE + hv ∇ϕv y λh = λF ρF . E∈A0h (Ωc )

v∈Vh (Σ)\{v0 }

F ∈Th (Σ)\{F0 }

Definimos las matrices complejas ı µ c 0 CE,E 0 := − c(qE , qE 0 ) , AΩ E,E 0 := ( µ qE , qE )0,Ωc , ωµ0 0 µ µ cΣ AΩ qE , ∇ϕv )0,Ωc , AΣ ∇ϕv , ∇ϕv0 )0,Ωc , v,v 0 := ( E,v := ( µ0 µ0 (3.3.2) 1 Rv,v0 := (rotΣ ϕv0 , Vτ (rotΣ ϕv))0,Σ , BF,v := (ρF , ( I − K)(γϕv))0,Σ , 2 VF,F 0 := (ρF 0 , VρF )0,Σ , para cualesquiera aristas interiores E, E 0 ∈ A0h (Ωc ) ; v´ertices sobre la frontera v, v 0 ∈ Vh (Σ) \ {v0 } y caras de la frontera F, F 0 ∈ Th (Σ) \ {F0 } . Con esta notaci´on, el Problema 3.2.1 es equivalente a determinar los vectores hΩc := (hE )E∈A0 (Ωc ) , h

hΣ := (hv )v∈Vh (Σ)\{v0 }

soluci´on del sistema lineal  AΩc + C   Ωc Σ > ) (A  0

y

λ := (λF )F ∈Th (Σ)\{F0 } ,

   f Ωc hΩc         AΣ + R B > hΣ =  f Σ  ,     0 −λ B V AΩc Σ

0



(3.3.3)

54

Cap´ıtulo 3

donde el t´ermino independiente viene dado por f Ωc := −

3.3.3.

ı (l(qE ))E∈A0 (Ωc ) h ωµ0

y

f Σ := −

ı (l(∇ϕv ))v∈Vh (Σ)\{v0 } . ωµ0

Implementaci´ on del esquema discreto

Para implementar el m´etodo (3.2.1), necesitamos calcular las integrales involucradas en el ensamblaje de la matriz y del segundo miembro del sistema lineal (3.3.3). Por esta raz´on, dise˜ namos f´ ormulas de cuadratura que sean adecuadas para aproximar las integrales que no podamos evaluar de forma exacta. Es importante subrayar que tales f´ormulas de cuadratura deben mantener un compromiso entre su complejidad computacional y el error cometido al aproximar la integral exacta por la integral num´erica. En el Ap´endice A detallamos la estrategia que seguimos. Adem´as, para implementar el m´etodo (3.2.1) tambi´en necesitamos resolver el sistema (3.3.3) de forma eficiente. Al respecto, recordemos que la matriz del sistema (3.3.3) es sim´etrica, no herm´ıtica, y definida positiva. En esta situaci´on, H.A. van der Vorst y J.B.M. Melissen proponen un m´etodo iterativo basado en una relaci´on recursiva de dos t´erminos muy similar a la del gradiente conjugado; cf. [54]. Tal y como mostramos en el siguiente apartado, los resultados num´ericos obtenidos con este m´etodo son satisfactorios; adem´as, nos gustar´ıa se˜ nalar su simplicidad de implementaci´on y bajo coste computacional.

3.4.

Resultados num´ ericos

Hemos implementado el m´etodo num´erico descrito en este cap´ıtulo mediante un c´odigo MATLAB. A continuaci´ on presentamos algunos resultados num´ericos obtenidos con dicho c´odigo para un problema test con soluci´on conocida expl´ıcitamente. As´ı validamos el c´odigo y comprobamos experimentalmente las propiedades de convergencia del m´etodo num´erico.

Eddy currents en conductores simplemente conexos

3.4.1.

55

Descripci´ on del problema test

N´otese que, si fijamos el campo e, entonces los campos h, j asociados al problema (1.2.1–1.2.2) quedan determinados como sigue: ı rot e ωµ j = rot h − σe

en R3 ,

h =

en Ωc , en Ωec .

j = rot h

Para que el problema test sea adecuado, necesitamos que el campo e tenga divergencia nula en el diel´ectrico Ωec y que cumpla las condiciones asint´ oticas (1.2.2). Adem´as, para poder aplicar nuestro m´etodo necesitamos que soporte(j) ⊂ Ωc . En otras palabras, debemos elegir el campo e de modo que 1 soporte(div e) ∪ soporte(rot( rot e)) ⊂ Ωc , µ y tal que e(x) = O(

1 ) y |x|

rot e(x) = O(

uniformemente en casi todas las direcciones es suficiente que e tenga la forma e := rot (f, f, f )>

1 ) |x| x |x| ,

cuando |x| → ∞ ,

con x ∈ R3 \ {0}. En particular,

en R3 ,

donde f : R3 → C es una funci´on suficientemente regular con soporte(f ) ⊂ Ωc .

3.4.2.

Resultados num´ ericos

A continuaci´on, mostramos algunos resultados num´ericos que hemos obtenido para un problema test como el descrito en el P´arrafo 3.4.1. M´as concretamente, suponemos que el conductor ocupa la regi´on Ωc := (−1, 1)3 y que es homog´eneo, con σ(x) = 1

y µ(x) = µ0 := 1

en Ωc .

56

Cap´ıtulo 3

Tomamos ω := 1 y  (x2 − x2 )α (x2 − x2 )α (x2 − x2 )α 0 1 0 2 0 3 f (x) := 0

si x ∈ (−x0 , x0 )3 , si no ,

siendo x0 ∈ (0, 1) y α ∈ (3, +∞) dos par´ametros reales. Para estudiar el orden de convergencia experimental del esquema num´erico, hemos resuelto los problemas discretos correspondientes a diferentes mallas, refinadas sucesivamente. En particular, la Figura 3.2 representa la malla m´as grosera que hemos empleado.

Figura 3.2: Malla m´as grosera

En la Tabla 3.1 presentamos los resultados obtenidos con nuestro m´etodo cuando x0 = 1 y α = 4. M´as concretamente, identificamos cada malla con su

Eddy currents en conductores simplemente conexos

57

di´ametro, que denotamos por h, y desglosamos el error total en las partes asociadas al campo magn´etico h y al potencial λ. Tambi´en presentamos el grado de convergencia experimental para el campo magn´etico, definido por ²=

logkh − hh krot,Ωc − logkh − heh krot,Ωc , log h − log e h

donde h y e h denotan los di´ametros de dos mallas consecutivas.

h

kh − hh krot,Ωc

kλ − λh k0,Σ

²

0.7733

34.8247

0.3136

—–

0.5330

25.0762

0.1634

0.8825

0.2989

14.0539

0.0257

1.0010

0.2337

11.4203

0.0131

0.8433

Tabla 3.1: Resultados num´ericos

Por otra parte, la Figura 3.3 representa el error correspondiente al campo magn´etico h frente al di´ametro h de la malla. Observamos que el error presenta una dependencia lineal con respecto a dicho di´ametro h, lo cual coincide con lo previsto por el Teorema 3.2.8 ya que la soluci´on es regular.

Cap´ıtulo 3

Error absoluto

58

error absoluto y=Ch 1

10

−0.7

10

−0.5

10

−0.3

10

−0.1

10

Diámetro de la malla

Figura 3.3: Error absoluto para hh en H(rot, Ωc ) vs. h (escala logar´ıtmica)

Cap´ıtulo 4

Eddy currents en conductores no simplemente conexos En este cap´ıtulo, proponemos un m´etodo BEM–FEM mixto para el problema de eddy currents en r´egimen arm´onico (1.2.1–1.2.2) sin restricciones topol´ogicas sobre el conductor. Tal y como se˜ nalamos en la Nota 3.1.1, cuando seguimos el razonamiento del Cap´ıtulo 3 encontramos dificultades para representar el campo magn´etico h en el diel´ectrico Ωec . En efecto, esta representaci´ on involucra campos vectoriales de Neumann arm´onicos y, en consecuencia, el correspondiente m´etodo num´erico es dif´ıcil de implementar y computacionalmente costoso; cf. [41]. Para evitar estas dificultades, seguimos la estrategia de [5] y reformulamos el problema de eddy currents como un problema mixto. Al mismo tiempo, esta estrategia nos permite estudiar aquellas situaciones en las que el soporte de la densidad de corriente j atraviesa la frontera Σ. En primer lugar introducimos un dominio computacional Ω conexo y simplemente conexo, tal que Ωc ⊂ Ω y

soporte(j) ⊂ Ω .

(4.0.1)

Adem´as, tomamos Ω con frontera Γ := ∂Ω conexa, simplemente conexa y Lips59

60

Cap´ıtulo 4

chitz continua. Subdividimos el diel´ectrico Ωec en el dominio acotado Ωd := Ω\Ωc y la regi´on no acotada Ωe := R3 \ Ω; v´ease la Figura 4.1 a continuaci´ on.

Ωe Γ j

Ωd Σ

Ωc

Figura 4.1: Geometr´ıa del problema Consideramos los operadores de trazas tangenciales −1/2

γ ext τ : H(rot, Ωd ) → H×

(divΣ , Σ)

y

−1/2

γ int τ : H(rot, Ωc ) → H×

(divΣ , Σ)

tomados desde Ωd y Ωc , respectivamente. Adem´as, denotamos por −1/2

−1 (γ ext τ ) : H×

(divΣ ,Σ) → H(rot,Ωd ) y

−1/2

−1 (γ int τ ) : H×

(divΣ ,Σ) → H(rot,Ωc )

int dos operadores continuos que sean inversas a derecha de γ ext τ y γ τ , respectivamente; v´ease el Teorema 2.5.2.

4.1.

Espacios funcionales auxiliares

Definimos el espacio © ª H(rot, R3 ) := q ∈ L2 (R3 ) ; rot q ∈ L2 (R3 ) ,

Problema arm´onico en conductores no simplemente conexos

61

que es de Hilbert con la norma del grafo ³ ´1/2 kqkrot,R3 := kqk20,R3 + krot qk20,R3 . Introducimos los espacios © ª X(R3 ) := q ∈ H(rot, R3 ) ; div q = 0 y rot q = 0 en Ωe , © ª U(R3 ) := q ∈ L2 (R3 ) ; div q = 0 y rot q = 0 en Ωe , as´ı como sus subespacios © ª q ∈ X(R3 ) ; rot q = 0 en Ωd , © ª U0 (R3 ) := q ∈ U(R3 ) ; rot q = 0 en Ωec . X0 (R3 ) :=

Observemos que X(R3 ) y X0 (R3 ) son subespacios cerrados de H(rot, R3 ) y, por lo tanto, son espacios de Hilbert con la norma k·krot,R3 . An´alogamente, U(R3 ) y U0 (R3 ) son subespacios cerrados de L2 (R3 ) y, en consecuencia, son espacios de Hilbert con la norma k·k0,R3 . Lema 4.1.1. El espacio X0 (R3 ) es denso en U0 (R3 ). Demostraci´ on. Sea q ∈ U0 (R3 ). Entonces q|Ωd ∈ H(rot, Ωd ) y podemos definir −1 ext p := (γ int τ ) (γ τ q|Ωd ) ∈ H(rot, Ωc ) .

Consideramos la descomposici´on q = q1 + q2 , donde  p en Ωc ∈ X0 (R3 ) , q1 := q| e en Ωe Ωc

q2 :=

 q|

Ωc

0

c

− p en Ωc en

Ωec

∈ U0 (R3 ) .

Como D(Ωc ) es denso en L2 (Ωc ), existe una sucesi´on (un )n∈N ⊂ D(Ωc ) tal que e n ∈ D(R3 ) la extensi´on un → q2 |Ωc en L2 (Ωc ). Para cada n ∈ N, denotamos por u e n )n∈N ⊂ X0 (R3 ) de un por cero a todo el espacio. Entonces la sucesi´on (q1 + u converge a q en L2 (R3 ) y concluimos el resultado. ¥

62

Cap´ıtulo 4

El resultado anterior nos permite introducir X00 (R3 ) como el espacio dual de X0 (R3 ) con espacio pivote U0 (R3 ). De esta forma, las siguientes inclusiones son continuas y densas: X0 (R3 ) ,→ U0 (R3 ) ,→ X00 (R3 ) . Nota 4.1.1 Recordemos que, seg´ un la primera condici´on de compatibilidad de e (1.2.3), div j = 0 en Ωc . En particular, j|Ωec ∈ H(div, Ωec ) y por lo tanto el Teorema 2.3.3 garantiza que γn (jd ) = γn (j|Ωe ) = 0 sobre Γ , donde denotamos jd := j|Ωd . La Nota anterior y las condiciones de compatibilidad (1.2.3) sugieren definir el siguiente subespacio cerrado de H(div, Ωd ): © M(Ωd ) := m ∈ L2 (Ωd ) ; div m = 0 en Ωd , ª γn m = 0 sobre Γ , hγn m, 1iΣj = 0 ∀j = 1, . . . , NΣ . Tambi´en definimos el espacio © ª H(div, R3 ) := q ∈ L2 (R3 ) ; div q ∈ L2 (R3 ) , que es de Hilbert con la norma del grafo ³ ´1/2 kqkdiv,R3 := kqk20,R3 + kdiv qk20,R3 . Consideramos el operador de extensi´on   ∇ρ(x)    Em(x) = m(x)     0

E : M(Ωd ) → H(div, R3 ) dado por si x ∈ Ωc , si x ∈ Ωd , si x ∈ Ωe ,

donde ρ ∈ H 1 (Ωc ) ∩ L20 (Ωc ) es arm´onica y se caracteriza de forma u ´nica por la condici´on de contorno ∂ρ = γnext m sobre Σ . ∂n

Problema arm´onico en conductores no simplemente conexos

63

Lema 4.1.2. El operador lineal E : M(Ωd ) → H(div, R3 ) est´ a acotado y satisface div(Em) = 0

en R3 ,

para cada m ∈ M(Ωd ). Demostraci´ on. Por definici´on, la funci´on Em ∈ L2 (R3 ) cumple div(Em) = 0 en Ωe ∪ Ωd ∪ Ωc , y adem´as γn (Em|Ωc ) =

∂ρ = γnext m = γn (Em|Ωd ) sobre Σ , ∂n

γn (Em|Ωd ) = γn m = 0 = γn (Em|Ωe )

sobre Γ .

Entonces el Corolario 2.3.3 garantiza que Em ∈ H(div, R3 ) y que div(Em) = 0 en R3 . Por otro lado, utilizando el lema de Lax–Milgram, sabemos que el problema ( hallar ϕ ∈ H 1 (Ωc ) ∩ L20 (Ωc ) t.q. (∇ϕ, ∇ψ)0,Ωc = hη, γψi1/2,Σ

∀ψ ∈ H 1 (Ωc ) ∩ L20 (Ωc ) ,

est´a bien planteado para cada η ∈ H −1/2 (Σ). N´otese que ρ ∈ H 1 (Ωc ) ∩ L20 (Ωc ) es ∂ρ ∈ H −1/2 (Σ). Luego la soluci´on de este problema para η := ∂n kρk1,Ωc ≤ C1 k

∂ρ k = C1 kγnext mk−1/2,Σ . ∂n −1/2,Σ

Como el operador lineal γnext : H(div, Ωd ) → H −1/2 (Σ) est´a acotado y adem´as div m = 0 en Ωd , deducimos que kρk1,Ωc ≤ C2 kmkH(div,Ωd ) = C2 kmk0,Ωd . As´ı concluimos que Em ∈ H(div, R3 ) tiene divergencia nula en R3 y cumple kEmk2div,R3 = k∇ρk20,Ωc + kmk20,Ωd ≤ (C22 + 1) kmk20,Ωd .

(4.1.1) ¥

64

Cap´ıtulo 4

Lema 4.1.3. El operador lineal M : M(Ωd ) → X(R3 ) ∩ H1 (R3 ), definido por ³ Z Em(y) ´ 1 Mm(x) := rot dy , 4π R3 |x − y| es continuo y adem´ as rot(Mm) = m en Ωd para cada m ∈ M(Ωd ). 1 es 4π |x| localmente integrable en R3 y que el soporte de Em ∈ L2 (R3 ) es compacto. Estas propiedades nos permiten introducir la distribuci´on Mm := rot u, donde u denota la convoluci´on Z Em(y) 1 u(x) := (Φ ∗ Em)(x) = dy . 4π R3 |x − y| Demostraci´ on. En primer lugar, se˜ nalamos que la funci´on Φ(x) :=

As´ı definida, la distribuci´on u satisface −∆u = Em en D 0 (R3 ) ;

(4.1.2)

al respecto, puede consultarse [56, P´arrafo VI.3]. Sea O un dominio acotado tal que Ω ⊂ O. Entonces u|O ∈ L2 (O) con Z ³ ´ kuk0,O ≤ sup Φ(x − y) dy kEmk0,R3 < ∞ . (4.1.3) x∈O

O

En efecto, utilizando la desigualdad de Cauchy–Swartz, ³Z ´ ³Z ´ 2 |u(x)| ≤ Φ(x − y) dy Φ(x − y) |Em(y)|2 dy en c.t.p. x ∈ R3 ; Ω

Ω

luego, por el criterio de mayoraci´on, u|O ∈ L2 (O) con Z ³ ´ Z ³Z ´ 2 kuk0,O ≤ sup Φ(x − y) dy Φ(x − y) dx |Em(y)|2 dy x∈O

O

Ω

O

y obtenemos la propiedad (4.1.3). Ahora, en virtud de las propiedades (4.1.2) y (4.1.3), podemos aplicar [35, Teorema 6.4] y deducir as´ı que ϕu ∈ H2 (O) para cada ϕ ∈ D(O). En particular, u|Ω ∈ H2 (Ω) con Z ³ ´ kuk2,Ω ≤ C (kuk0,O + kEmk0,O ) ≤ C sup Φ(x − y) dy + 1 kEmk0,R3 . x∈O

O

Problema arm´onico en conductores no simplemente conexos

65

Entonces, utilizando la propiedad (4.1.1), kMmk1,Ω ≤ kuk2,Ω ≤ C1 kmk0,Ωd .

(4.1.4)

Por otro lado, tenemos que div u = Φ∗(div Em) = 0 en R3 ; as´ı que, aplicando la tercera identidad de (2.0.1), rot(Mm) = −∆u + ∇(div u) = Em en R3 .

(4.1.5)

En particular, se cumple la propiedad del enunciado: rot(Mm) = m en Ωd . Tambi´en tenemos, en virtud de (4.1.5), que rot(Mm) = 0 en Ωe . Luego el Teorema 2.7.3 garantiza que (Mm)|Ωe = ∇ψ para un potencial ψ ∈ W 1 (Ωe ). Como div(Mm) = 0 en Ωe , deducimos que la funci´on ψ ∈ W 1 (Ωe ) es arm´onica y es la soluci´on del problema variacional  1 e   hallar ψ ∈ W (Ω ) tal que Z  ∇ψ · ∇ϕ dx = hγn (rot u), γϕi1/2,Γ ∀ϕ ∈ W 1 (Ωe ) .  Ωe

En virtud del Lema de Lax–Milgram, este problema est´a bien planteado y, en consecuencia, kMmk0,Ωe = k∇ψk0,Ωe ≤ C2 kγn (rot u)k−1/2,Γ . Utilizando la Proposici´on 2.3.2 junto con la propiedad (4.1.4), resulta que kMmk0,Ωe ≤ C3 krot uk0,Ω ≤ C4 kmk0,Ωd .

(4.1.6)

Teniendo en cuenta que div(Mm) = 0 en R3 y aplicando los resultados (4.1.4), (4.1.5) y (4.1.6), deducimos que Mm ∈ H(rot, R3 ) ∩ H(div, R3 ) con kMmkrot,R3 + kMmkdiv,R3 ≤ C5 kmk0,Ωd . En [29, Nota 2.6] se demuestra que el espacio H(rot, R3 ) ∩ H(div, R3 ) se incluye con continuidad en H1 (R3 ). Por lo tanto, Mm ∈ H1 (R3 ) con kMmk1,R3 ≤ C6 kmk0,Ωd y concluimos que el operador M : M(Ωd ) → X(R3 ) ∩ H1 (R3 ) est´a acotado.

¥

66

4.2.

Cap´ıtulo 4

Formulaci´ on variacional

En este apartado deducimos una formulaci´ on variacional del problema (1.2.1– 1.2.2). M´as concretamente, primero obtenemos una formulaci´ on variacional en 3 t´erminos del campo magn´etico en todo el espacio R . A continuaci´ on, introducimos un multiplicador de Lagrange para imponer la condici´on asociada al rotacional del campo magn´etico en el diel´ectrico. Finalmente, introducimos una variable auxiliar sobre la frontera Γ y as´ı podemos reformular el problema en el dominio acotado Ω.

4.2.1.

Una formulaci´ on variacional global

Comenzamos estudiando el espacio funcional adecuado para buscar el campo magn´etico h. Con este fin, notamos que, seg´ un la primera ecuaci´on de (1.2.1), rot h = 0 en Ωe . Adem´as, aplicando el operador divergencia en la segunda ecuaci´on de (1.2.1), div h = 0 en Ωe . Partiendo de estas dos ecuaciones junto con la condici´on asint´ otica (1.2.2), y utilizando la tercera identidad de (2.0.1), deducimos que   ∆h = 0 en R3 \ BR ,     div h = 0 en R3 \ BR ,   1    h(x) = O( ) uniformemente cuando |x| → ∞ , |x| donde BR representa una bola centrada en el origen con radio R suficientemente grande. En [6, Proposici´on 3.1] se hace un desarrollo del campo h en polinomios esf´ericos arm´onicos para mostrar que, bajo las condiciones anteriores, h(x) = O(

1 ) uniformemente cuando |x| → ∞ . |x|2

Problema arm´onico en conductores no simplemente conexos

67

En particular, tenemos que h ∈ L2 (R3 ). Esta propiedad, junto con la primera ecuaci´on de (1.2.1) y el Lema 4.1.3, sugiere buscar el campo h en el espacio af´ın Mjd + X0 (R3 ). Para buscar una formulaci´ on variacional del problema (1.2.1–1.2.2), tomamos el producto escalar de la segunda ecuaci´on de (1.2.1) por una funci´on gen´erica q ∈ X0 (R3 ), integramos en el dominio Ω y usamos la f´ormula de Green (2.3.3): ıω (µh, q)0,R3 + (e, rot q)0,Ωc = 0 .

(4.2.1)

La primera ecuaci´on de (1.2.1) y la hip´otesis (1.1.6) garantizan que e=

1 (rot h − jc ) en Ωc , σ

siendo jc := j|Ωc . As´ı que podemos eliminar el campo el´ectrico de (4.2.1) y obtenemos ıω(µh, q)0,R3 + c(h, q) = l(q) ∀q ∈ X0 (R3 ) , donde 1 c(q1 , q2 ) := ( rot q1 , rot q2 )0,Ωc σ

1 y l(q) := ( jc , rot q)0,Ωc , σ

para cualesquiera q1 , q2 , q ∈ X(R3 ). Por lo tanto, proponemos la siguiente formulaci´ on d´ebil del problema (1.2.1) en t´erminos del campo magn´etico. Problema 4.2.1 Buscar h ∈ Mjd + X0 (R3 ) tal que ıω(µh, q)0,R3 + c(h, q) = l(q) ∀q ∈ X0 (R3 ) . Teorema 4.2.1. Si jd ∈ M(Ωd ), entonces el Problema 4.2.1 est´ a bien planteado. Demostraci´ on. En t´erminos de la variable auxiliar h∗ := h − Mjd ∈ X0 (R3 ), el problema 4.2.1 es equivalente a   buscar h∗ ∈ X0 (R3 ) tal que (4.2.2)  ıω(µh∗ , q) 3 + c(h∗ , q) = l∗ (q) ∀q ∈ X (R3 ) . 0 0,R

68

Cap´ıtulo 4

Aqu´ı, el operador l∗ : X0 (R3 ) → C se define como l∗ (q) := l(q) − ıω (µMjd , q)0,R3 − c(Mjd , q) ∀q ∈ X0 (R3 ) . N´otese que la siguiente forma sesquilineal es el´ıptica: X0 (R3 ) × X0 (R3 ) → C (q1 , q2 ) 7→ (1 + ı) (

µ 1−ı q1 , q2 )0,R3 + c(q1 , q2 ) . µ0 ω µ0

En efecto, razonando como en la demostraci´on del Teorema 3.1.3, deducimos que · ¸ ¢ 1−ı¡ b)0,R3 + c(q, q) ≥ α kqk2rot,R3 ∀q ∈ X0 (R3 ) , Re ıω(µb q, q ω µ0 donde

½ α := m´ın 1,

¾ 1 > 0. ωµ0 σ1 En consecuencia, podemos reescalar el problema (4.2.2) de modo que est´e bajo las hip´otesis del Lema de Lax–Milgram. Como este problema es equivalente al Problema 4.2.1, concluimos el resultado. ¥

4.2.2.

Una formulaci´ on variacional mixta

En este apartado imponemos la condici´on rot h|Ωd = jd introduciendo un multiplicador de Lagrange. Con este fin, definimos la forma sesquilineal d(m, q) := (m, rot q)0,Ωd

∀q ∈ X(R3 ), m ∈ M(Ωd ).

Sea Υ : X(R3 ) → M0 (Ωd ) el operador lineal caracterizado por hΥq, miM0 (Ωd )×M(Ωd ) := d(m, q) ∀q ∈ X(R3 ), m ∈ M(Ωd ) , donde M0 (Ωd ) representa el espacio dual de M(Ωd ) con espacio pivote L2 (Ωd ). N´otese que este operador est´a acotado; en efecto, basta observar que kΥqkM0 (Ωd ) =

|d(m, q)| ≤ krot qk0,Ωd m∈M(Ωd ), m6=0 kmk0,Ωd sup

∀q ∈ X(R3 ) .

Denotamos por Υ∗ : M(Ωd ) → X0 (R3 ) el operador dual de Υ: hΥ∗ m, qiX0 (R3 )×X(R3 ) = d(m, q) ∀q ∈ X(R3 ), m ∈ M(Ωd ) .

Problema arm´onico en conductores no simplemente conexos

69

Lema 4.2.2. El operador lineal Υ∗ : M(Ωd ) → X0 (R3 )0 es un isomorfismo, siendo © ª X0 (R3 )0 := ` ∈ X0 (R3 ) ; h`, qiX0 (R3 )×X(R3 ) = 0 ∀q ∈ X0 (R3 ) . Demostraci´ on. En virtud del Lema 4.1.3, sup q∈X(R3 ), q6=0

kmk20,Ωd |d(m, q)| |d(m, Mm)| ≥ = kqkrot,R3 kMmkrot,R3 kMmkrot,R3

∀m ∈ M(Ωd ) .

Entonces, utilizando que el operador M : M(Ωd ) → X(R3 ) ∩ H1 (R3 ) est´a acotado, resulta que sup q∈X(R3 ), q6=0

1 |d(m, q)| ≥ kmk0,Ωd kqkrot,R3 kMk

∀m ∈ M(Ωd ) .

Adem´as, la inclusi´on © ª rot q|Ωd ; q ∈ X(R3 ) ⊂ M(Ωd ) garantiza que X0 (R3 ) es el n´ ucleo del operador Υ : X(R3 ) → M0 (Ωd ). Por lo tanto, el resultado es consecuencia de [29, Lema 4.1].

¥

Problema 4.2.2 Buscar h ∈ X(R3 ) y r ∈ M(Ωd ) tales que ıω(µh, q)0,R3 + c(h, q) + d(r, q) = l(q) d(m, h) = (jd , m)0,Ωd

∀q ∈ X(R3 ) , ∀m ∈ M(Ωd ) .

Teorema 4.2.3. Bajo las hip´ otesis del Teorema 4.2.1, el Problema 4.2.2 tiene 3 una u ´nica soluci´ on h ∈ X(R ), r ∈ M(Ωd ). Adem´ as, h pertenece al espacio af´ın 3 Mjd + X0 (R ) y es la u ´nica soluci´ on del Problema 4.2.1. Demostraci´ on. La segunda ecuaci´on del Problema 4.2.2 significa que rot h = jd en Ωd , es decir, h ∈ Mjd +X0 (R3 ). Luego, testeando la primera ecuaci´on de 4.2.2 con funciones q ∈ X0 (R3 ), deducimos que h coincide con la u ´nica soluci´on del

70

Cap´ıtulo 4

Problema 4.2.1. As´ı que el resultado se reduce a mostrar la existencia y unicidad de un multiplicador r ∈ M(Ωd ) tal que d(r, q) = l(q) − ıω(µh, q)0,R3 − c(h, q) ∀q ∈ X(R3 ) .

(4.2.3)

Para ello, introducimos el funcional G ∈ X0 (R3 ) definido por hG, qiX0 (R3 )×X(R3 ) := l(q) − ıω(µh, q)0,R3 − c(h, q) ∀q ∈ X(R3 ) . Como h es la soluci´on del Problema 4.2.1, sabemos que hG, qiX0 (R3 )×X(R3 ) = 0 para cada q ∈ X0 (R3 ). Es decir G ∈ X0 (R3 )0 y, por tanto, el Lema 4.2.2 garantiza que la ecuaci´on Υ∗ r = G tiene una u ´nica soluci´on r ∈ M(Ωd ). Ahora concluimos notando que esta u ´ltima ecuaci´on es equivalente a (4.2.3). ¥

4.2.3.

Una formulaci´ on variacional BEM–FEM

Dado un vector q ∈ X(R3 ), el Lema 2.7.2 garantiza que podemos representarlo en el dominio exterior como q = ∇ϕ en Ωe ,

(4.2.4)

para un u ´nico potencial arm´onico ϕ ∈ W 1 (Ωe ). En particular, como ϕ es arm´onico en Ωe , admite una representaci´on integral de la forma (2.7.2) y cumple las ecuaciones integrales (2.7.8). Definimos n o X(Ω) := q ∈ H(rot, Ω) ; γn (rot q) = 0 en H −1/2 (Γ) , que es un subespacio cerrado de H(rot, Ω); v´ease la Proposici´on 3.1.2. Adem´as, razonando como en el P´arrafo 3.1.2, tenemos que rotΓ : H 1/2 (Γ)/C → γτ X(Ω)

(4.2.5)

Problema arm´onico en conductores no simplemente conexos

71

es un isomorfismo y que el operador lineal 1 1/2 B := ( I − K) ◦ rot−1 (Γ)/C Γ ◦γτ : X(Ω) → H 2 est´a acotado; tambi´en tenemos la siguiente estimaci´on an´aloga a (3.1.10): kBqk1/2,Γ ≤ C kγτ qkH−1/2 (Γ)

∀q ∈ X(Ω).

(4.2.6)

×

−1/2 b Introducimos el espacio X(Ω) := X(Ω)×H0 (Γ) y las formas sesquilineales

a0 (q1 , q2 ) := ıω(µq1 , q2 )0,Ω + ıωµ0 hγτ q2 , Vτ (γτ q1 )iτ ,Γ , ³ ´ b2 ) := a0 (q1 , q2 ) + ıωµ0 hη2 , Vη1 i1/2,Γ + hη2 , Bq1 i1/2,Γ − hη1 , Bq2 i1/2,Γ , a(b q1 , q b b1 := (q1 , η1 ) y q b2 := (q2 , η2 ) en X(Ω). para cualesquiera q N´otese que, bajo la hip´otesis (1.1.5) sobre el coeficiente µ, los Teoremas 2.7.5 y 2.7.6 y la propiedad (4.2.6) garantizan que ³ b2 )| ≤ C1 kq1 k0,Ω kq2 k0,Ω + kγτ q1 kH−1/2 (Γ) kγτ q2 kH−1/2 (Γ) |a(b q1 , q ×

×

+ kη1 k−1/2,Γ kη2 k−1/2,Γ + kη2 k−1/2,Γ kγτ q1 kH−1/2 (Γ) ×

´ + kη1 k−1/2,Γ kγτ q2 kH−1/2 (Γ) , (4.2.7) ×

b) ∈ R con y que adem´as a(b q, q ³ b) ≥ C2 kqk20,Ω + kγτ qk2 a(b q, q

−1/2

H×

(Γ)

´ + kηk2−1/2,Γ ,

(4.2.8)

b b1 := (q1 , η1 ), q b2 := (q2 , η2 ) y q b = (q, η) en X(Ω). para cualesquiera q Proponemos la siguiente formulaci´ on variacional del problema (1.2.1–1.2.2). b := (h, λ) ∈ X(Ω) b Problema 4.2.3 Buscar h y r ∈ M(Ωd ) tales que b q b) + d(r, q) + c(h, q) = l(q) a(h, d(m, h) = (jd , m)0,Ωd

b ∀b q ∈ X(Ω) , ∀m ∈ M(Ωd ) .

72

Cap´ıtulo 4

A continuaci´on damos un resultado auxiliar que nos permitir´a analizar el Problema 4.2.3 y relacionarlo con el Problema 4.2.2. Proposici´ on 4.2.4. El siguiente operador de restricci´ on define un isomorfismo: R : X(R3 ) → X(Ω) q 7→ Rq := q|Ω . Demostraci´ on. Claramente R : X(R3 ) → X(Ω) es un operador lineal continuo. Para mostrar que R es inyectivo, consideramos q ∈ X(R3 ) con Rq = 0; es decir, q = 0 en Ω. Sea ϕ ∈ W 1 (Ωe ) el potencial arm´onico asociado a q como en (4.2.4). N´otese que, por la propia definici´on del operador rotΓ , rotΓ ϕ = γτ (∇ϕ) = γτ q = 0 . Luego el Teorema 2.6.3 garantiza que la traza de ϕ es constante sobre Γ. Ahora, como la funci´on ϕ ∈ W 1 (Ωe es arm´onica y su traza es constante sobre Γ, deducimos que ϕ = 0 en Ωe . Por lo tanto, q se anula en todo R3 . Para mostrar que R es suprayectivo, tomamos q ∈ X(Ω). Como (4.2.5) define un isomorfismo, existe un potencial η ∈ H 1/2 (Γ) con rotΓ η = γτ q y este potencial es u ´nico salvo constante aditiva. Sea ϕ ∈ W 1 (Ωe ) la u ´nica fune ci´on arm´onica en Ω cuya traza coincide con η sobre Γ salvo constante aditiva. Definimos  q en Ω , p := ∇ϕ en Ωe , que cumple γτint p = γτint q = rotΓ η = γτext p sobre Γ . Entonces el Teorema 2.3.6 garantiza que p ∈ X(R3 ), con Rp = q. Ahora podemos aplicar el teorema del grafo cerrado y as´ı concluimos que R : X(R3 ) → X(Ω) es un isomorfismo. ¥

Problema arm´onico en conductores no simplemente conexos

73

Teorema 4.2.5. Bajo las hip´ otesis del Teorema 4.2.1, el Problema 4.2.3 tiene b ∈ X(Ω), b una u ´nica soluci´ on h r ∈ M(Ωd ). Adem´ as, se relaciona con la soluci´ on ext b = (Rh, γ h) y el multiplicador de Lagrange es el h del Problema 4.2.2 seg´ un h n mismo para ambos problemas. Demostraci´ on. Sea h ∈ X(R3 ), r ∈ M(Ωd ) la soluci´on del Problema 4.2.2. b∗ = (h∗ , λ) := (Rh, γ ext h) ∈ X(Ω), b Entonces h r ∈ M(Ωd ) es una soluci´on del n Problema 4.2.3: b∗ ∈ X(Ω) b Es claro que h satisface la segunda ecuaci´on del Problema 4.2.3, as´ı que estudiamos si tambi´en satisface la primera ecuaci´on. Para ello, consideramos una funci´on gen´erica q ∈ X(R3 ) y denotamos por ψ, ϕ ∈ W 1 (Ωe ) los potenciales arm´onicos asociados a los campos h, q ∈ X(R3 ) v´ıa (4.2.4). Teniendo en cuenta la definici´on de rotΓ y la continuidad de la traza tangencial a trav´es de Γ, γτ h∗ = rotΓ ψ

y

γτ q = rotΓ ϕ .

(4.2.9)

∂ψ Definimos λ := γnext h = . Entonces, utilizando las ecuaciones integrales (2.7.8) ∂n para el potencial arm´onico ψ y recordando la definici´on del operador B, Vλ + Bh∗ = 0

y

1 λ = −N γψ + ( I − K∗ )λ 2

sobre Γ .

(4.2.10)

A continuaci´on aplicamos la f´ormula de Green (2.3.1) en (µh, q)0,R3 = (µh∗, q)0,Ω + µ0 (∇ψ, ∇ϕ)0,Ωe , de modo que (µh, q)0,R3 = (µh∗, q)0,Ω − µ0 hλ, ϕi1/2,Γ . Combinando esta ecuaci´on con la segunda de (4.2.10) y con (2.7.7), 1 (µh, q)0,R3 = (µh∗, q)0,Ω + µ0 hrotΓ ϕ, Vτ (rotΓ ψ)iτ ,Γ − µ0 hλ, ( I − K)γϕi1/2,Γ . 2 Utilizando (4.2.9) junto con la definici´on del operador B, deducimos que (µh, q)0,R3 = (µh∗, q)0,Ω + µ0 hγτ q, Vτ (γτ h∗ )iτ ,Γ − µ0 hλ, Bqi1/2,Γ .

74

Cap´ıtulo 4

Luego, usando la primera identidad de (4.2.10), b∗, q b b) ∀b ıω (µh, R−1 q)0,R3 = a(h q := (q, η) ∈ X(Ω) . En consecuencia, la primera ecuaci´on del Problema 4.2.2 garantiza que b∗ , q b b) + c(h∗ , q) + d(r, q) = l(q) ∀b a(h q := (q, η) ∈ X(Ω) , que es la primera ecuaci´on del Problema 4.2.3. b∗ ∈ X(Ω), b Por lo tanto, h r ∈ M(Ωd ) es una soluci´on del Problema 4.2.3. Para demostrar la unicidad de soluci´on del Problema 4.2.3, consideramos b = (h, λ) ∈ X(Ω), b h r ∈ M(Ωd ) que sea soluci´on del Problema 4.2.3 homog´eneo. Testeando la primera ecuaci´on del Problema 4.2.3 homog´eneo con funciones de b b := (0, η) ∈ X(Ω), la forma q deducimos que λ = −V −1 (Bh) .

(4.2.11)

Testeando nuevamente la primera ecuaci´on del Problema 4.2.3 homog´eneo con b b := (q, 0) ∈ X(Ω) funciones de la forma q y substituyendo λ por la expresi´on (4.2.11), ıω (µh, q)0,Ω + ıωµ0 hγτ q, Vτ (γτ h)iτ ,Γ + ıωµ0 hV −1 (Bh), Bqi1/2,Γ + d(r, q) + c(h, q) = 0 ∀q ∈ X(Ω) . Razonando como en la demostraci´on de la Proposici´on 4.2.4, podemos asociar al campo h ∈ X(Ω) un u ´nico potencial arm´onico ψ ∈ W 1 (Ωe ) caracterizado por la condici´on rotΓ (γψ) = γτ h. An´alogamente, podemos asociar a cada funci´on q ∈ X(Ω) un u ´nico potencial arm´onico ϕ ∈ W 1 (Ωe ) tal que rotΓ (γϕ) = γτ q. Como ψ es arm´onico, cumple las ecuaciones integrales (2.7.8) y por tanto la u ´ltima ecuaci´on es equivalente a 1 ∂ψ ıω (µh, q)0,Ω + ıωµ0 hN γψ, γϕi1/2,Γ − ıωµ0 h( I − K∗ ) , γϕi1/2,Γ 2 ∂n + c(h, q) + d(r, q) = 0 ∀q ∈ X(Ω) .

Problema arm´onico en conductores no simplemente conexos

75

Aplicando nuevamente que ψ satisface (2.7.8) e integrando por partes, ıω (µh, q)0,Ω + ıωµ0 (∇ψ, ∇ϕ)0,Ωe + c(h, q) + d(r, q) = 0 ∀q ∈ X(Ω) . Es decir, ıω (µR−1 h, q)0,R3 + c(R−1 h, q) + d(r, q) = 0

∀q ∈ X(R3 ) .

Teniendo en cuenta esta identidad y la segunda ecuaci´on del Problema 4.2.3 homog´eneo, resulta que R−1 h ∈ X(R3 ), r ∈ M(Ωd ) es una soluci´on del Problema 4.2.2 homog´eneo. Entonces el Teorema 4.2.1 garantiza que R−1 h = 0 en R3 y r = 0 en Ωd , y concluimos que el Problema 4.2.3 homog´eneo s´olo admite la soluci´on trivial. ¥ Introducimos los operadores lineales continuos A : X(Ω) → X0 (Ω) ,

C : X(Ω) → X0 (Ω)

y D : X(Ω) → M0 (Ωd ) ,

definidos por hAq1 , q2 iX0 (Ω)×X(Ω) := a0 (q1 , q2 ) ,

hCq1 , q2 iX0 (Ω)×X(Ω) := c(q1 , q2 ) ,

hDq, miM0 (Ωd )×M(Ωd ) := d(m, q) , para cualesquiera q1 , q2 , q ∈ X(Ω) y m ∈ M(Ωd ). Adem´as, denotamos por −1/2

B ∗ : H0

(Γ) → X0 (Ω) y D ∗ : M(Ωd ) → X0 (Ω)

los operadores duales de B y D, respectivamente. Entonces reformulamos el Problema 4.2.3 como      A + C −ıωµ0 B∗ D ∗ h `           (4.2.12) µ0 V 0  λ =  0  ,  ıωµ0 B      D 0 0 r `d siendo h`, qiX0 (Ω)×X(Ω) := l(q)

y

h`d , miM0 (Ωd )×M(Ωd ) := (jd , m)0,Ωd

76

Cap´ıtulo 4

para cualesquiera q ∈ X(Ω) y m ∈ M(Ωd ). Utilizando la expresi´on (4.2.11), podemos eliminar λ del sistema (4.2.12) y reescribirlo de forma reducida como      A + ıωµ0 B∗ V −1 B + C D ∗ h `    =  . (4.2.13) D 0 r `d

4.3.

Esquema discreto

En lo que sigue, suponemos que Ω y Ωc son poliedros de Lipschitz. Consideramos {Th (Ω); h > 0} una familia regular de mallas del dominio Ω tales que cada o en Ωd . Denotamos por h elemento T ∈ Th (Ω) es un tetraedro contenido en Ωc ´ el m´aximo di´ametro de los tetraedros de Th (Ω). Adem´as, definimos Th (Ωd ) como la restricci´on de Th (Ω) al dominio Ωd , ª © Th (Ωd ) := T ∈ Th (Ω) ; T ⊂ Ωd , y Th (Γ) como la triangulaci´on que induce Th (Ω) sobre la superficie Γ, Th (Γ) := {F cara de alg´ un T ∈ Th (Ω) ; F ⊂ Γ} .

4.3.1.

Descripci´ on del problema discreto

En este apartado, proponemos una versi´ on discreta de (4.2.3) combinando elementos finitos y elementos de contorno para aproximar conformemente el b espacio X(Ω). M´as concretamente, consideramos el espacio de elementos finitos de arista © ª N D h (Ω) := q ∈ H(rot, Ω) ; q|T ∈ N D(T ) ∀T ∈ Th (Ω) ,

(4.3.1)

y aproximamos las funciones de X(Ω) usando las del espacio de dimensi´on finita Xh (Ω) := N D h (Ω) ∩ X(Ω) ;

Problema arm´onico en conductores no simplemente conexos

77

v´eanse las definiciones (3.2.1) y (3.2.4). −1/2

Adem´as, aproximamos las funciones de H0 (Γ) con funciones que son constantes a trozos con respecto a la triangulaci´on Th (Γ). De hecho, utilizamos el espacio de dimensi´on finita © ª Λh (Γ) := η ∈ L2 (Γ) ; hη, 1iΓ = 0 , η|F ∈ C ∀F ∈ Th (Γ) ; v´ease la definici´on (3.2.5). Finalmente, hacemos una discretizaci´on conforme del espacio H(div, Ωd ) utilizando elementos finitos de Raviart–Thomas de primer orden: RT h (Ωd ) := {q ∈ H(div, Ω) ; q|T ∈ RT (T ) ∀T ∈ Th (Ωd )} , © ª siendo RT (T ) := a x + b ; a ∈ C, b ∈ C3 ; v´ease [45, Secci´on 5.4]. Entonces aproximamos las funciones de M(Ωd ) usando las del espacio de dimensi´on finita Mh (Ωd ) := RT h (Ωd ) ∩ M(Ωd ) . Proponemos la siguiente versi´ on discreta del Problema 4.2.3, donde denotab mos Xh (Ω) := Xh (Ω) × Λh (Γ). b h := (hh , λh ) ∈ X b h (Ω) y rh ∈ Mh (Ωd ) con Problema 4.3.1 Buscar h bh, q b) + c(hh , q) + d(rh (t), q) = l(q) a(h

b h (Ω) , ∀b q = (q, η) ∈ X

d(m, hh ) = (jd , m)0,Ωd ∀m ∈ Mh (Ωd ) . Nota 4.3.1 Tal y como se˜ nalamos en la Nota 3.2.1, la inversa del operador tangencial rotΓ est´a involucrada en la definici´on de la forma sesquilineal a(·, ·) y, en consecuencia, nuestro m´etodo num´erico parece dif´ıcil de implementar y computacionalmente costoso. De nuevo, esto no sucede en la pr´actica porque dicho operador puede eliminarse en el c´alculo efectivo de la matriz del sistema lineal asociado al Problema 4.3.1. Al respecto, v´ease el P´arrafo 4.4.

78

Cap´ıtulo 4

4.3.2.

An´ alisis del esquema discreto

Introducimos los espacios X0 (Ω) := {q ∈ X(Ω) ; rot q = 0 en Ωd } , X0,h (Ω) := {q ∈ Xh (Ω) ; d(m, q) = 0 ∀m ∈ Mh (Ωd )} , que son subespacios cerrados de X(Ω) y de Xh (Ω), respectivamente. Adem´as, en [45, Lema 5.40] se muestra que {rot q|Ωd ; q ∈ Xh (Ω)} ⊂ Mh (Ωd ); en consecuencia, tenemos que X0,h (Ω) ⊂ X0 (Ω) . Denotamos por I h el operador de interpolaci´on asociado al espacio de elementos finitos de N´ed´elec N D h (Ω). Recordemos que en el P´arrafo 3.2.1 definimos expl´ıcitamente este operador y se˜ nalamos sus propiedades fundamentales. Por otra parte, sea I RT el operador de interpolaci´on asociado al espacio h de elementos finitos de Raviart–Thomas RT h (Ωd ). En [45, Teorema 5.25] se demuestra que, para cada s ∈ (1/2, 1], existe una constante C > 0 tal que s kq − I RT h qk0,Ωd ≤ C h kqks,Ωd

∀q ∈ Hs (Ωd ) .

(4.3.2)

Adem´as, en [45, Lema 5.40] se prueba que, para todo s > 1/2, s rot(I h q)|Ωd = I RT h (rot q|Ωd ) ∀q ∈ H (rot, Ω) .

(4.3.3)

El siguiente resultado es una herramienta fundamental para el an´alisis del esquema discreto (4.3.1). Lema 4.3.1. Existe una constante β > 0 tal que sup q∈Xh (Ω), q6=0

para todo m ∈ Mh (Ωd ).

|d(m, q)| ≥ β kmk0,Ωd , kqkrot,Ω

Problema arm´onico en conductores no simplemente conexos

79

Demostraci´ on. Seg´ un el Lema 4.1.3, la funci´on Mm|Ωd ∈ H1 (Ωd ) satisface que rot(Mm) = m en Ωd . En particular, para cada tetraedro T ∈ Th (Ωd ) tenemos que Mh m|T ∈ H1 (T ) y, en consecuencia, la traza de Mm sobre ∂T pertenece al espacio Lp (∂T ) para todo p ≥ 2. En estas condiciones, [7, Lema 4.7] permite definir Mh m := I dh (Mm|Ωd ), donde I dh denota el operador de interpolaci´on de N´ed´elec asociado al espacio © ª N D h (Ωd ) := q ∈ H(rot, Ω) ; q|T ∈ N D(T ) ∀T ∈ Th (Ωd ) . Por un lado, razonando como en la demostraci´on de [7, Proposici´on 4.6], deducimos que ³ ´ kMh mk0,Ωd ≤ C1 h krot(Mm)k0,Ωd + kMmk1,Ωd . (4.3.4) Por otro lado, la propiedad (4.3.3) garantiza que rot(Mh m) = I RT h (rot(Mm)) = m en Ωd .

(4.3.5)

Luego, utilizando el Lema 4.1.3 junto con (4.3.4–4.3.5), tenemos que kMh mkrot,Ωd ≤ C2 kmk0,Ωd . Adem´as, en [2] se muestra que existe un operador de extensi´on discreto N D h (Ωd ) → N D h (Ω) que es uniformemente continuo con respecto al par´ametro h > 0. As´ı que existe f h m ∈ Xh (Ω) que extiende Mh m a todo Ω y que cumple una funci´on M f h mkrot,Ω ≤ C3 kMh mkrot,Ω ≤ C3 C2 kmk0,Ω . kM d d

(4.3.6)

Por lo tanto, combinando (4.3.5) y (4.3.6), concluimos que sup q∈Xh (Ω), q6=0

f h m)| |d(m, M |d(m, q)| 1 ≥ kmk0,Ωd . ≥ f h mkrot,Ω kqkrot,Ω C2 C3 kM ¥

80

Cap´ıtulo 4

Teorema 4.3.2. El Problema 4.3.1 tiene una u ´nica soluci´ on. Demostraci´ on. En primer lugar, el Lema 4.3.1 junto con [29, Teorema 1.1] garantizan la existencia de una funci´on hdh ∈ Xh (Ω) tal que d(m, hdh ) = (jd , m)0,Ωd

∀m ∈ Mh (Ωd ).

Entonces definimos h∗h := hh − hdh ∈ X0,h (Ω) y reescribimos el Problema 4.3.1 en b ∗ := (h∗ , λh ) ∈ X0,h (Ω) × Λh (Γ) como t´erminos de h h h ∗ ∗ b∗ , q a(h q ∈ X0,h (Ω) × Λh (Γ) , h b ) + c(hh , q) = l (q) ∀b

(4.3.7)

donde b) − c(hdh , q) . l∗ (q) := l(q) − a((hdh , 0), q Introducimos los operadores lineales y acotados −1/2

Ah : X(Ω) → X0h (Ω) ,

Vh : H0

D ∗h : M(Ωd ) → X0h (Ω) ,

D h : X(Ω) → M0h (Ωd ) , Bh : X(Ω) → Λ0h (Γ) ,

(Γ) → Λ0h (Γ) ,

−1/2

B∗h : H0

(Γ) → X0h (Ω) ,

C h : X(Ω) → X0h (Ω) , caracterizados por hAh q, qh iX0h (Ω)×Xh (Ω) = a0 (q, qh ) ,

hVh η, ηh iΛ0h (Γ)×Λh (Γ) = hVη, η h i1/2,Γ ,

hD h q, mh iM0h (Ωd )×Mh (Ωd ) = d(mh , q) ,

hD ∗h m, qh iX0h (Ω)×Xh (Ω) = d(m, qh ) ,

hBh q, ηh iΛ0h (Γ)×Λh (Γ) = hBq, η h i1/2,Γ ,

hB∗h η, qh iX0h (Ω)×Xh (Ω) = hBqh , ηi1/2,Γ ,

hC h q, qh iX0h (Ω)×Xh (Ω) = c(q, qh ) , −1/2

respectivamente, para cualesquiera q ∈ X(Ω), m ∈ M(Ωd ), η ∈ H0 (Γ) y para todo qh ∈ Xh (Ω), mh ∈ Mh (Ωd ), ηh ∈ Λh (Γ). Adem´as, seg´ un el Teorema 2.7.5, hVh ηh , η h iΛ0h (Γ)×Λh (Γ) = hVηh , ηh i1/2,Γ ≥ α1 kηh k2−1/2,Γ

∀ηh ∈ Λh (Γ) ,

Problema arm´onico en conductores no simplemente conexos

81

de modo que Vh : Λh (Γ) → Λ0h (Γ) es un isomorfismo. Esta propiedad nos permite desacoplar el sistema (4.3.7) como      Ah + ıωµ0 B∗h Vh−1 Bh + C h 0 h∗h `∗    =  , (4.3.8) Vh−1 Bh I λh −Vh−1 Bh hdh donde `∗ ∈ X00,h (Ω) est´a definido por h`∗ , qiX0h (Ω)×Xh (Ω) = l∗ (q) − ıωµ0 hB∗h Vh−1 Bh hdh , qiX0h (Ω)×Xh (Ω)

∀q ∈ X0,h (Ω) .

Tambi´en tenemos que · ¸ 1−ı ∗ −1 0 Re h(Ah + ıωµ0 Bh Vh Bh + C h ) q, qiXh (Ω)×Xh (Ω) ≥ ω µ0 µ 1 1 ≥ ( q, q)0,Ω + ( rot q, rot q)0,Ωc , µ0 ωµ0 σ para cada q ∈ X0,h (Ω), as´ı que el operador Ah + ıωµ0 B∗h Vh−1 Bh + C h : X0,h (Ω) → X00,h (Ω) es invertible. En consecuencia, queda garantizada la existencia y unicidad de la inc´ognita h∗h ∈ X0,h (Ω) y, por lo tanto, tambi´en de la funci´on λh = −Vh−1 Bh (h∗h + hdh ) .

(4.3.9)

Para estudiar la existencia y unicidad del multiplicador a nivel discreto, rh ∈ Mh (Ωd ), introducimos el operador lineal Fh : Xh (Ω) → C b h , (q, 0)) . q 7→ Fh (q) := l(q) − c(hh , q) − a(h

(4.3.10)

b h es equivalente a Con esta notaci´on, la ecuaci´on que satisface h Fh (q) = 0 ∀q ∈ X0,h (Ω) .

(4.3.11)

Luego, recordando la definici´on de X0,h (Ω) y aplicando el Lema 4.3.1, concluimos que existe un u ´nico rh ∈ Mh (Ωd ) con d(rh , q) = Fh (q) ∀q ∈ Xh (Ω) .

(4.3.12) ¥

82

Cap´ıtulo 4

4.3.3.

An´ alisis de la convergencia

Teniendo en cuenta el Teorema 2.7.5, podemos definir el operador lineal −1/2

ρ : X(Ω) → H0

(Γ)

q 7→ ρq := V −1 Bq , que es continuo. Recordemos que el operador V : H −1/2 (Γ) → H 1/2 (Γ) es el´ıptico 1/2 (v´ease el Teorema 2.7.5) y que el operador B : X(Ω) → H0 (Γ) es continuo; −1/2 as´ı que el Lema de Lax–Milgram garantiza que la funci´on ρq ∈ H0 (Γ) se caracteriza de forma u ´nica por cumplir hη, Vρqi1/2,Γ = hη, Bqi1/2,Γ

−1/2

∀η ∈ H0

(Γ) .

(4.3.13)

Adem´as, en virtud de (4.2.6), tenemos que kρqk−1/2,Γ ≤ C1 kγτ qkH−1/2 (Γ) .

(4.3.14)

×

Ahora definimos el operador ρh : X(Ω) → Λh (Γ) utilizando una versi´ on discreta de la ecuaci´on variacional anterior. Es decir, ρh q ∈ Λh (Γ) se caracteriza de forma u ´nica por hη, Vρh qiΓ = hη, BqiΓ

∀η ∈ Λh (Γ) .

(4.3.15)

En efecto, utilizando nuevamente la elipticidad de V y la propiedad (4.2.6), deducimos que el operador ρh est´a bien definido y que kρh qk−1/2,Γ ≤ C2 kγτ qkH−1/2 (Γ) .

(4.3.16)

×

Adem´as, se cumple la siguiente desigualdad tipo Lema de C´ea: kρq − ρh qk−1/2,Γ ≤ C3 kρq − ηk−1/2,Γ

∀η ∈ Λh (Γ) ,

para cualquier q ∈ X(Ω). Lema 4.3.3. Dado s ∈ (1/2, 1], existe una constante C > 0 tal que kρq − ρh qk−1/2,Γ ≤ C h1/2 kqkHs (rot,Ω) , para cada q ∈ X(Ω) ∩ Hs (rot, Ω).

(4.3.17)

Problema arm´onico en conductores no simplemente conexos

83

Demostraci´ on. Fijado s > 1/2, sabemos que γτ : Hs (rot, Ω) → L2τ (Γ) est´a acotado; cf. [21, Teorema 3]. Utilizando el isomorfismo definido en (4.2.5), deducimos que existe un potencial η ∈ H 1/2 (Γ) con rotΓ η = γτ q ∈ L2τ (Γ) . N´otese que este potencial se caracteriza por ser la u ´nica soluci´on del problema variacional   buscar η ∈ H 1 (Γ)/C tal que (4.3.18) 1  hrot η, rot λi = hγ q, rot λi Γ Γ Γ Γ Γ ∀λ ∈ H (Γ)/C . τ Como η ∈ H 1 (Γ), tenemos que 1 Vρq = Bq = ( I − K)η ∈ H 1 (Γ)/C . 2 Entonces el Teorema 2.7.5 garantiza que ρq ∈ L20 (Γ) con ³ ´ 1 kρqk0,Γ ≤ C kρqk−1/2,Γ + k( I − K)ηk1,Γ ≤ C kqkHs (rot,Ω) . 2 Finalmente, teniendo en cuenta esta estimaci´on, la desigualdad (4.3.17) y la propiedad de aproximaci´on ´ınf

kλ − λh k−1/2,Γ ≤ C h1/2 kλk0,Γ

λh ∈Λh (Γ)

∀λ ∈ L20 (Γ) ,

concluimos que kρq − ρh qk−1/2,Γ ≤ C h1/2 kρqk0,Γ ≤ C h1/2 kqkHs (rot,Ω) . ¥ A continuaci´on estudiaremos el orden de convergencia del esquema discreto dado por el Problema 4.3.1 cuando hacemos hip´otesis de regularidad sobre la soluci´on exacta del Problema 4.2.3.

84

Cap´ıtulo 4

Teorema 4.3.4. Sean h y hh las primeras componentes de la soluciones de los Problemas 4.2.3 y 4.3.1, respectivamente. Si h ∈ X(Ω) ∩ Hs (rot, Ω) para alg´ un s > 1/2, entonces tenemos la siguiente estimaci´ on del error: kh − hh krot,Ω ≤ C h1/2 khkHs (rot,Ω) . Demostraci´ on. En esta demostraci´on, seguimos el razonamiento expuesto en [29, Teorema 1.1]. Comenzamos notando que el Lema 4.3.1 permite definir ph ∈ Xh (Ω) con D h ph = D h (h − I h h) .

(4.3.19)

Entonces la funci´on hdh := ph + I h h ∈ Xh (Ω) cumple d(m, hdh ) = (jd , m)0,Ωd

∀m ∈ Mh (Ωd ) .

Descomponemos h − hh como h − hh = e1h − e2h , donde e1h := h − hdh

y

e2h := hh − hdh .

A continuaci´on acotamos la norma de los campos e1h y eh2 en H(rot, Ω). Por un lado, aplicando la desigualdad triangular para e1h = h − I h h − ph , ke1h krot,Ω ≤ (1 +

kdk ) kh − I h hkrot,Ω , β

(4.3.20)

donde kdk > 0 es la norma de la forma sesquilineal d(·, ·) y β > 0 es la constante involucrada en el Lema 4.3.1. Por otro lado, la primera ecuaci´on del sistema (4.2.13) implica que a0 (h, q) + ıωµ0 hρh, Bqi1/2,Γ + c(h, q) = l(q) ∀q ∈ X0 (Ω) . An´alogamente, la primera ecuaci´on de (4.3.8) garantiza que a0 (hh , q) + ıωµ0 hρh hh , Bh qi1/2,Γ + c(hh , q) = l(q) ∀q ∈ X0,h (Ω) .

Problema arm´onico en conductores no simplemente conexos

85

Puesto que X0,h (Ω) ⊂ X0 (Ω), deducimos de las dos u ´ltimas identidades que a0 (e2h , q) + c(e2h , q) + ıωµ0 hVh−1 Bh e2h , Bh qiΛ0h (Γ)×Λh (Γ) = = a0 (e1h , q) + c(e1h , q) + ıωµ0 hρh − ρh h, Bqi1/2,Γ + ıωµ0 hρh e1h , Bqi1/2,Γ para cualquier q ∈ X0,h . En particular, se cumple para q = e2h ∈ X0,h (Ω) y entonces, aplicando (4.2.7), (4.2.6) y (4.3.16), |a0 (e2h , e2h ) + µ0 hB∗h Vh−1 Bh e2h , e2h iX0h (Ω)×Xh (Ω) + c(e2h , e2h )| ≤ ´ ³ ≤ C1 ke1h krot,Ω + kρh − ρh hk−1/2,Γ ke2h krot,Ω . Utilizando la hip´otesis (1.1.6), la propiedad (4.2.8) y que el operador B∗h Vh−1 Bh es no negativo, deducimos que ³ ´ ke2h k2rot,Ω ≤ C2 ke1h krot,Ω + kρh − ρh hk−1/2,Γ ke2h krot,Ω . Simplificando y aplicando el Lema 4.3.3, ³ ´ ke2h krot,Ω ≤ C3 ke1h krot,Ω + h1/2 khkHs (rot,Ω) .

(4.3.21)

Finalmente, aplicando la desigualdad triangular para h − hh = e1h − e2h junto con las cotas (4.3.20) y (4.3.21), resulta que ³ ´ kh − hh krot,Ω ≤ C4 kh − I h hkrot,Ω + h1/2 khkHs (rot,Ω) , y concluimos utilizando la estimaci´on del error de interpolaci´on del Teorema 3.2.2. ¥ Corolario 4.3.5. Bajo las hip´ otesis del Teorema 4.3.4, tenemos que kλ − λh k−1/2,Γ ≤ C h1/2 khkHs (rot,Ω) . Demostraci´ on. Partiendo de la primera ecuaci´on del Problema 4.2.3, dedujimos la ecuaci´on (4.2.11), que es equivalente a λ = −ρh. An´alogamente, partiendo de

86

Cap´ıtulo 4

la primera ecuaci´on del Problema 4.3.1 obtenemos que λh = −ρh hh . Luego, utilizando la desigualdad triangular, kλ − λh k−1/2,Γ ≤ kρh − ρh hk−1/2,Γ + kρh (h − hh )k−1/2,Γ . Adem´as, la propiedad (4.3.16) garantiza que kρh (h − hh )k−1/2,Γ ≤ C1 kh − hh krot,Ω . As´ı que concluimos el resultado aplicando el Lema 4.3.3 y el Teorema 4.3.4.

¥

Corolario 4.3.6. Supongamos que se cumplen las hip´ otesis del Teorema 4.3.4 y r que r ∈ M(Ωd ) ∩ H (Ωd ) para alg´ un r > 1/2. Entonces ³ ´ kr − rh k0,Ωd ≤ C h1/2 khkHs (rot,Ω) + krkr,Ωd . Demostraci´ on. Introducimos el operador lineal continuo F : X(Ω) → C b (q, 0)) , q 7→ F(q) := l(q) − c(h, q) − a(h, y recordemos que definimos una versi´ on discreta Fh : Xh (Ω) → C en (4.3.10). N´otese que los multiplicadores de Lagrange r ∈ M(Ωd ) y rh ∈ Mh (Ωd ) satisfacen d(r, q) = F(q) ∀q ∈ X(Ω)

y

d(rh , q) = Fh (q) ∀q ∈ Xh (Ω) .

Estas propiedades sugieren tomar una funci´on auxiliar e rh ∈ Mh (Ωd ) tal que d(e rh , q) = F(q) ∀q ∈ Xh (Ω) . En efecto, el Lema 4.3.1 garantiza que esta funci´on existe y es u ´nica, porque el 0 operador F ∈ X (Ω) se anula sobre X0 (Ω) y, en particular, sobre X0,h (Ω). Descomponemos r − rh como r − rh = e1h + e2h , donde e1h := r − e rh

y e2h := e rh − rh .

Problema arm´onico en conductores no simplemente conexos

87

A continuaci´on acotamos la norma de los campos e1h y e2h en L2 (Ωd ). Por un lado, utilizando la desigualdad triangular para e1h = r − s + s − e rh y aplicando el Lema 4.3.1, ke1h k0,Ωd ≤ kr − sk0,Ωd +

1 β

d(e rh − s, q) q∈Xh (Ω), q6=0 kqkrot,Ω sup

∀s ∈ Mh (Ωd ) .

Como d(e rh , q) = d(r, q) para todo q ∈ Xh (Ω), deducimos que ke1h k0,Ωd ≤ (1 +

kdk ) β

´ınf

kr − sk0,Ωd .

s∈Mh (Ωd )

Adem´as, bajo la hip´otesis de regularidad r ∈ M(Ωd ) ∩ Hr (Ωd ), la propiedad (4.3.3) y el Teorema 4.1.3 implican que RT I RT h r = I h (rot(Mr)) = rot(Mh r) ∈ Mh (Ωd ) ,

siendo Mh r := I dh (Mr|Ωd ). Luego ke1h k0,Ωd ≤ (1 +

kdk ) kr − I RT h rk0,Ωd . β

(4.3.22)

Por otro lado, seg´ un el Lema 4.3.1, ke2h k0,Ωd ≤

1 β

|d(e rh − rh , q)| 1 ≤ kqkrot,Ω β q∈Xh (Ω), q6=0 sup

|F(q) − Fh (q)| . kqkrot,Ω q∈Xh (Ω), q6=0 sup

Observemos que

³ ´ |F(q) − Fh (q)| ≤ C1 k(h − hh )krot,Ω + kλ − λh k−1/2,Γ kqkrot,Ω

∀q ∈ Xh (Ω) .

As´ı que, aplicando el Teorema 4.3.4 y el Corolario 4.3.5, ke2h k0,Ωd ≤

1 β

|F(q) − Fh (q)| ≤ C2 h1/2 khkHs (rot,Ω) . kqk rot,Ω q∈Xh (Ω), q6=0 sup

(4.3.23)

Finalmente, utilizando la desigualdad triangular para r − rh = e1h + e2h junto con las cotas (4.3.22) y (4.3.23), tenemos que kr − rh k0,Ωd ≤ (1 +

kdk 1/2 ) kr − I RT khkHs (rot,Ω) . h rk0,Ωd + C2 h β

Por lo tanto, concluimos el resultado aplicando la estimaci´on del error de interpolaci´on (4.3.2). ¥

88

Cap´ıtulo 4

4.4.

Forma matricial del esquema discreto

En este apartado reformulamos el Problema 4.3.1 como un sistema lineal. b h (Ω) Para ello comenzamos determinando expl´ıcitamente una base del espacio X y un sistema generador de Mh (Ωd ). Sea Ah (Ω) el conjunto de las aristas de la malla Th (Ω). Definimos los conjuntos A0h (Ω) y Ah (Ωd ) de aristas interiores a Ω y Ωd , respectivamente; es decir, © ª A0h (Ω) := {E ∈ Ah (Ω) ; E ∩ Ω 6= ∅} y Ah (Ωd ) := E ∈ Ah (Ω) ; E ⊆ Ωd . Tambi´en introducimos el conjunto Vh (Ω) de v´ertices de la malla Th (Ω), as´ı como Vh (Γ) := {v ∈ Vh (Ω) ; v ∈ Γ} . En la Proposici´on 3.3.1 vimos que, si fijamos un v´ertice v0 ∈ Vh (Γ), entonces el siguiente conjunto determina una base de Xh (Ω): © ª Bh (Ω) := qE ; E ∈ A0h (Ω) ∪ {∇ϕv ; v ∈ Vh (Γ) \ {v0 }} . Adem´as, en el P´arrafo 3.3.1 se˜ nalamos que, si fijamos una cara F0 ∈ Th (Γ), entonces el siguiente conjunto constituye una base de Λh (Γ): ¾ ½ 1 1 χF − χF ; F ∈ Th (Γ) \ {F0 } . ρF := |F | |F0 | 0 En lo que respecta al espacio M(Ωd ), n´otese que {(rot q)|Ωd ; q ∈ Xh (Ω)} ⊆ Mh (Ωd ) . Por otro lado, la inclusi´on rec´ıproca es consecuencia de la condici´on inf–sup del Lema 4.3.1. En efecto, para cada rh ∈ Mh (Ωd ), dicha condici´on inf–sup garantiza la existencia de alg´ un qh ∈ Xh (Ω) tal que d(qh , m) = (rh , m)0,Ωd

∀m ∈ Mh (Ωd ) ;

es decir, rh = (rot qh )|Ωd . En consecuencia, {(rot q)|Ωd ; q ∈ Xh (Ω)} = Mh (Ωd )

Problema arm´onico en conductores no simplemente conexos

89

y, por lo tanto, el conjunto {(rot qE )|Ωd ; E ∈ Ah (Ωd )} es un sistema generador de Mh (Ωd ). Expresamos las inc´ognitas del Problema 3.2.1 en t´erminos de las bases de los espacios Xh (Ω) y Λh (Γ): X X X hh = hE qE + hv ∇ϕv y λh = λF ρF ; E∈A0h (Ω)

v∈Vh (Γ)\{v0 }

F ∈Th (Γ)\{F0 }

y del sistema generador de Mh (Ωd ): X rh = rE (rot qE )|Ωd . E∈Ah (Ωd )

Recordemos que en (3.3.2) introdujimos las matrices complejas µ ı qE , qE 0 )0,Ω , CE,E 0 := − c(qE , qE 0 ) , µ0 ωµ0 µ µ qE , ∇ϕv )0,Ω , AΓv,v0 := ( ∇ϕv , ∇ϕv0 )0,Ω , AΩΓ E,v := ( µ0 µ0 1 := (rotΓ ϕv0 , Vτ (rotΓ ϕv ))0,Γ , BF,v := (ρF , ( I − K)γϕv )0,Γ , 2 AΩ E,E 0 := (

Rv,v0

VF,F 0 := (ρF 0 , VρF )0,Γ , para cualesquiera aristas interiores E, E 0 ∈ A0h (Ω) ; v´ertices sobre la frontera v, v 0 ∈ Vh (Γ) \ {v0 } y caras de la frontera F, F 0 ∈ Th (Γ) \ {F0 } . Ahora, definimos adem´as la matriz DE,E 0 := −

ı (rot qE , rot qE 0 )0,Ωd , ωµ0

para cualesquiera E ∈ A0h (Ω) y E 0 ∈ Ah (Ωd ). Con esta notaci´on, el Problema 4.3.1 es equivalente a determinar los vectores hΩ := (hE )E∈A0 (Ω) , h

λ := (λF )F ∈Th (Γ)\{F0 } ,

hΓ := (hv )v∈Vh (Γ)\{v0 } , r := (rE )E∈A0 (Ωd ) , h

90

Cap´ıtulo 4

soluci´on del sistema lineal      AΩ + C AΩΓ 0 D hΩ fΩ       ΩΓ >  Σ    Γ + R B> (A ) A 0 h 0        = ,       0     B V 0 −λ 0       ⊥ Ω d D 0 0 0 r f

(4.4.1)

donde el t´ermino independiente viene dado por f Ω := −

4.4.1.

ı (l(qE ))E∈A0 (Ω) h ωµ0

y

f Ωd := −

´ ı ³ d (j , rot qE )0,Ωd . ωµ0 E∈Ah (Ωd )

Implementaci´ on del esquema discreto

Para implementar nuestro m´etodo num´erico, necesitamos aproximar las integrales asociadas a la matriz y al segundo miembro del sistema lineal (4.4.1). Con este objetivo, seguimos la estrategia descrita en el P´arrafo 3.3.3 y el Ap´endice A. Es importante recordar que dedujimos el sistema (4.4.1) utilizando, no una base del espacio Mh (Ωd ), sino un sistema que es generador pero no libre. En consecuencia, dicho sistema lineal es indeterminado. Esta dificultad puede solucionarse introduciendo un multiplicador de Lagrange adicional para imponer la restricci´on de divergencia nula asociada al espacio M(Ωd ); v´ease [5, Secci´on 5.2]. Ahora bien, esta estrategia aumenta el tama˜ no del sistema lineal y, en consecuencia, el coste computacional de su ensamblaje y su resoluci´on. Por esta raz´on, pensamos que en la pr´actica es preferible resolver el sistema indeterminado (4.4.1) con un m´etodo iterativo del tipo gradiente conjugado; v´ease [5, Nota 5.1].

4.5.

Resultados num´ ericos

Hemos implementado el m´etodo descrito en este cap´ıtulo mediante un c´odigo MATLAB. A continuaci´on presentamos algunos resultados num´ericos que hemos

Problema arm´onico en conductores no simplemente conexos

91

obtenido con dicho c´odigo para el problema test del P´arrafo 3.4.1 cuando Ωc := [(−1, 1) × (−1, 1) × (−0,5, 0,5)] \ (−0,5, 0,5)3 ,  (x2 − |x|2 )α si |x| < x0 , 0 f (x) := 0 si no , siendo x0 > 0 y α ∈ [3, +∞) dos par´ametros reales. En lo que respecta a los datos f´ısicos, de nuevo consideramos  1 si x ∈ Ωc , σ(x) := 0 si x ∈ Ωe , c

µ := 1

en R3 ,

ω := 1

en R3 .

Adem´as, tomamos como dominio computacional Ω := (−1,5, 1,5) × (−1,5, 1,5) × (−1, 1) . Para estudiar el orden de convergencia del esquema discreto, hemos resuelto el Problema 4.3.1 para diferentes mallas, refinadas sucesivamente. En la Tabla 4.1 mostramos los resultados que hemos obtenido cuando x0 = 1 y α = 4. M´as concretamente, identificamos cada malla con su di´ametro, que denotamos por h, y desglosamos el error total en las partes asociadas al campo magn´etico h, al potencial λ y al multiplicador de Lagrange r ∈ M(Ωd ). Tambi´en mostramos el orden de convergencia experimental para el campo magn´etico, definido por ²=

logkh − hh krot,Ω − logkh − heh krot,Ω , log h − log e h

donde h y e h denotan los di´ametros de dos mallas consecutivas. Por otra parte, la Figura 4.2 muestra en escala logar´ıtmica el error correspondiente al campo magn´etico h frente al di´ametro h de la malla. Observamos que este error presenta una dependencia lineal con respecto a dicho di´ametro

92

Cap´ıtulo 4

h

kh − hh krot,Ω

kλ − λh k0,Γ

kr − rh k0,Ωd

²

0.9299

56.1051

0.7107

52.0686

—–

0.6601

39.3998

0.1413

24.4914

1.0315

0.4572

28.6757

0.0620

21.4884

0.8651

0.3653

22.0725

0.0173

11.7229

1.1663

Tabla 4.1: Resultados num´ericos h, mejorando el resultado previsto por el Teorema 4.3.4. De hecho, para este problema test, nuestro m´etodo se comporta como un m´etodo FEM mixto (sin elementos de contorno); v´ease la nota a continuaci´ on. Nota 4.5.1 Hemos comparado los resultados num´ericos obtenidos utilizando nuestro m´etodo BEM–FEM mixto con los resultados num´ericos que se obtienen aplicando el m´etodo num´erico propuesto en [5]. M´as concretamente, dicho m´etodo es de tipo FEM mixto y consiste en resolver el sistema lineal (4.4.1) bloqueando a cero las inc´ognitas hΓ y λ. En [5, P´arrafo 6], se se˜ nala que este esquema num´erico presenta el orden de convergencia O(h) siempre que la soluci´on exacta sea suficientemente regular. Para el problema test anterior y empleando las mismas mallas, hemos observado que los resultados num´ericos obtenidos con ambos m´etodos son similares para las inc´ognitas hΩ y r. Esto sucede porque la inc´ognita de frontera λ se anula, y adem´as el campo magn´etico h y el multiplicador r tienen sus soportes contenidos en el dominio computacional Ω.

Problema arm´onico en conductores no simplemente conexos

93

1.8

10

1.7

Error absoluto

10

1.6

10

1.5

10

error absoluto y=Ch 1.4

10

−0.4

10

−0.3

10

−0.2

10

−0.1

10

Diámetro de la malla

Figura 4.2: Error absoluto para hh en H(rot, Ωc ) vs. h (escala logar´ıtmica)

Cap´ıtulo 5

Eddy currents de evoluci´ on En este cap´ıtulo estudiamos el problema de eddy currents de evoluci´ on (1.1.1–1.1.3) bajo las condiciones de compatibilidad (1.1.7). Para ello, utilizamos una estrategia an´aloga a la del Cap´ıtulo 4. M´as concretamente, suponemos que J y div H0 tienen soporte compacto. Entonces podemos introducir un dominio computacional Ω conexo y simplemente conexo tal que Ωc ⊂ Ω ,

soporte(J) ⊂ [0, T ] × Ω

y

soporte(div H0 ) ⊂ Ω .

Adem´as, tomamos Ω de forma que su frontera Γ := ∂Ω sea conexa, simplemente conexa y Lipschitz continua. Subdividimos el diel´ectrico Ωec en el dominio acotado Ωd := Ω\Ωc y la regi´on no acotada Ωe := R3 \ Ω; v´ease la Figura 4.1. Es importante se˜ nalar que los campos involucrados en el problema de evoluci´on (1.1.1–1.1.3) son reales, a diferencia de lo que sucede para el problema en r´egimen arm´onico (1.2.1–1.2.2). Por lo tanto, en este cap´ıtulo trabajamos con funciones reales en vez de funciones complejas. Adem´as, para cada funci´on f : (0, T ) × R3 → R, denotamos f (t) : R3 → R x 7→ f (t)(x) := f (t, x) . 95

96

5.1.

Cap´ıtulo 5

Formulaci´ on variacional

En este p´arrafo deducimos una formulaci´ on variacional del problema (1.1.1) con un razonamiento an´alogo al del Cap´ıtulo 4.

5.1.1.

Una formulaci´ on variacional global

En primer lugar, estudiamos el espacio funcional en el que buscaremos el campo magn´etico H. Con este fin, notamos que, seg´ un la primera ecuaci´on de (1.1.1), rot H(t, x) = J(t, x) en (0, T ) × Ωe . Adem´as, tomando el operador divergencia en la segunda ecuaci´on de (1.1.1) y utilizando la hip´otesis (1.1.5), ∂t (div H(t, x)) = 0 en (0, T ) × Ωe . Recordemos que hemos introducido el dominio computacional Ω ⊂ R3 de modo que J(t, x) = 0 en (0, T ) × Ωe y div H(0, x) = div H0 = 0 en Ωe , luego rot H(t, x) = 0 y

div H(t, x) = 0

en (0, T ) × Ωe .

Partiendo de esta propiedad junto con la condici´on asint´ otica (1.1.3) y aplicando la tercera identidad de (2.0.1), deducimos que   ∆H(t, x) = 0 en (0, T ) × (R3 \ BR ) ,     div H(t, x) = 0 en (0, T ) × (R3 \ BR ) ,   1    H(t, x) = O( ) cuando |x| → ∞ , |x| siendo BR una bola centrada en el origen con radio R suficientemente grande y donde el comportamiento asint´otico es uniforme en casi todas las direcciones x (t, |x| ) (con t ∈ (0, T ), x ∈ R3 \ {0}). En esta situaci´on, en [6, Proposici´on 3.1] se hace una expansi´on del campo H(t) en polinomios esf´ericos arm´onicos para mostrar que H(t, x) = O(

1 ) uniformemente cuando |x| → ∞ . |x|2

Problema de evoluci´on

97

En particular, tenemos que H(t) ∈ L2 (R3 ) en c.t.p. t ∈ (0, T ). Teniendo en cuenta esta propiedad, as´ı como la primera ecuaci´on de (1.1.1) y el Lema 4.1.3, buscaremos el campo magn´etico H(t) en el espacio af´ın M(J(t)|Ωd ) + X0 (R3 ) en c.t.p. t ∈ (0, T ). En lo que sigue, denotamos por simplicidad Jd := J|(0,T )×Ωd y MJd (t) := M(Jd (t)) en c.t.p. t ∈ (0, T ). Para deducir una formulaci´ on variacional del problema (1.1.1–1.1.3), testeamos la segunda ecuaci´on de (1.1.1) con una funci´on q ∈ X0 (R3 ) e integramos por partes con la f´ormula (2.3.3): d (µH(t), q)0,R3 + (E(t), rot q)0,Ωc = 0 dt

en D0 (0, T ).

(5.1.1)

Seg´ un la primera ecuaci´on de (1.1.1), E=

1 (rot H − Jc ) σ

en (0, T ) × Ωc ,

donde denotamos Jc := J|(0,T )×Ωc . As´ı que podemos eliminar el campo el´ectrico de (5.1.1) como sigue: d (µH(t), q)0,R3 + c(t, H(t), q) = l(t, q) en D0 (0, T ) , dt siendo c(t, q1 , q2 ) := (

1 rot q1 , rot q2 )0,Ωc σ(t)

y

l(t, q) := (

1 Jc (t), rot q)0,Ωc , σ(t)

para cualesquiera funciones q1 , q2 , q ∈ X(R3 ). Introducimos el espacio W 1 ((0, T ); X0 (R3 ), X00 (R3 )) , que consta de las funciones q ∈ L2 ((0, T ); X0 (R3 )) con ∂t q ∈ L2 ((0, T ); X00 (R3 )). Entonces proponemos la siguiente formulaci´ on d´ebil del problema de eddy currents de evoluci´on (1.1.1–1.1.3).

98

Cap´ıtulo 5

Problema 5.1.1 Buscar H ∈ MJd + W 1 ((0, T ); X0 (R3 ), X00 (R3 )) tal que H(0) = H0 , y que satisfaga la siguiente ecuaci´on en D0 (0, T ): d (µ H(t), q)0,R3 + c(t, H(t), q) = l(t, q) ∀q ∈ X0 (R3 ) . dt N´otese que la condici´on inicial del problema anterior tiene sentido siempre que Jd sea suficientemente regular. De hecho, el Lema 4.1.1 junto con [58, Teorema 41.15] garantizan que el espacio W 1 ((0, T ); X0 (R3 ), X00 (R3 )) se incluye con continuidad en C 0 ([0, T ]; U0 (R3 )). Teorema 5.1.1. Supongamos que la densidad de corriente J ∈ L2 ((0, T ) × R3 ) presenta la regularidad Jd ∈ H 1 ((0, T ); M(Ωd )) y que el campo magn´etico inicial H0 ∈ U(R3 ) satisface la condici´ on de compatibilidad rot H0 |Ωd = Jd (0). Entonces el Problema 5.1.1 tiene una soluci´ on H ∈ MJd + W 1 ((0, T ); X0 (R3 ), X00 (R3 )), que es u ´nica y se caracteriza por Z − 0

Z

T

(µ H(t), ∂t q(t))0,R3 dt +

T

c(t, H(t), q(t)) dt = Z T = (µ H0 , q(0))0,R3 + l(t, q(t)) dt , 0

(5.1.2)

0

para cada q ∈ C ∞ ([0, T ]; X0 (R3 )) con q(T ) = 0. Demostraci´ on. Introducimos la variable auxiliar H∗ := H − MJd y reescribimos el Problema 5.1.1 en t´erminos de esta nueva inc´ognita. Con este fin, notamos que la hip´otesis de regularidad Jd ∈ H 1 ((0, T ); M(Ωd )) junto con el Lema 4.1.2 implican que ∂t E(Jd (t)) = E(∂t Jd (t)) en c.t.p. t ∈ (0, T ). Luego ∂t (MJd (t)) = ∂t (Φ ∗ E(Jd (t))) = Φ ∗ E(∂t Jd (t)) = M(∂t Jd (t)) ∈ X(R3 ) ,

Problema de evoluci´on

99

en c.t.p. t ∈ (0, T ) y, aplicando el Lema 4.1.3, resulta que MJd ∈ H 1 (0, T ; X(R3 )). Como consecuencia de esta propiedad, deducimos que el operador l∗ (t, q) := l(t, q) − (µ ∂t (MJd (t)), q)0,R3 − c(t, MJd (t), q) pertenece al espacio L2 ((0, T ); X00 (R3 )). Tambi´en deducimos que tiene sentido evaluar MJd (t) en t = 0, porque el espacio H 1 ((0, T ); X(R3 )∩H1 (R3 )) se incluye con continuidad en C 0 ([0, T ]; X(R3 ) ∩ H1 (R3 )); cf. [58, Teorema 41.15]. Estas propiedades nos permiten reescribir el Problema 5.1.1 como  buscar H∗ ∈ W 1 ((0, T ); X0 (R3 ), X00 (R3 )) tal que        H∗ (0) = H0 − MJd (0) ,  y que satisface la siguiente ecuaci´on en D0 (0, T ) :     d   (µH∗ (t), q)0,R3 + c(t, H∗ (t), q) = l∗ (t, q) ∀q ∈ X0 (R3 ) . dt

(5.1.3)

N´otese que la forma bilineal X0 (R3 ) × X0 (R3 ) → R (q1 , q2 ) 7→ c(t, q1 , q2 ) est´a acotada y es coerciva uniformemente con respecto a t ∈ (0, T ): c(t, q1 , q2 ) = ( c(t, q, q) =

1 rot q1 , rot q2 )0,Ωc σ(t) 1 rot q, rot q)0,Ωc ( σ(t)

≤ ≥

1 kq1 krot,R3 kq2 krot,R3 , σ0 ´ 1 ³ kqk2rot,R3 − kqk20,R3 , σ1

para cualesquiera q1 , q2 , q ∈ X0 (R3 ) y en c.t.p. t ∈ (0, T ). En estas condiciones, [57, Teorema 23A] garantiza que el problema (5.1.3) est´a bien planteado. Finalmente, concluimos notando que (5.1.3) es equivalente al Problema 5.1.2. ¥

100

5.1.2.

Cap´ıtulo 5

Una formulaci´ on variacional mixta

Siguiendo la estrategia del P´arrafo 4.2.2, en este apartado imponemos la restricci´on rot H|(0,T )×Ωd = Jd introduciendo un multiplicador de Lagrange. Con este fin, definimos Υ : L2 ((0, T ); X(R3 )) → L2 ((0, T ); M0 (Ωd )) tal que Z T Z T hΥq(t), m(t)iM0 (Ωd )×M(Ωd ) dt = d(m(t), q(t)) dt , 0

0

para cualesquiera q ∈ L2 ((0, T ); X(R3 )) y m ∈ L2 ((0, T ); M(Ωd )). N´otese que este operador est´a acotado; en efecto, para cada q ∈ L2 ((0, T ); X(R3 )), kΥqk2L2 ((0,T );M0 (Ωd ))

Z TÃ = 0

d(m, q(t)) sup m∈M(Ωd ), m6=0 kmk0,Ωd

!2

Z dt ≤ 0

T

krot q(t)k20,Ωd dt .

Denotamos su dual por Υ∗ : L2 ((0, T ); M(Ωd )) → L2 ((0, T ); X0 (R3 )); es decir, Z T Z T ∗ hΥ m(t), q(t)iX0 (R3 )×X(R3 ) dt = d(m(t), q(t)) dt, 0

0

para cualesquiera q ∈ L2 ((0, T ); X(R3 )) y m ∈ L2 ((0, T ); M(Ωd )). Razonando como en el Lema 4.2.2, deducimos el siguiente resultado. Lema 5.1.2. El operador lineal Υ∗ : L2 ((0, T ); M(Ωd )) → L2 ((0, T ); X0 (R3 ))0 define un isomorfismo, donde L2 ((0, T ); X0 (R3 ))0 representa el espacio de las funciones L ∈ L2 ((0, T ); X0 (R3 )) tales que Z T hL(t), q(t)iX0 (R3 )×X(R3 ) dt = 0 0

para cada q ∈ L2 ((0, T ); X0 (R3 )). Problema 5.1.2 Buscar H ∈ W 1 ((0, T ); X(R3 ), X0 (R3 ))

y

r ∈ L2 ((0, T ); M(Ωd ))

Problema de evoluci´on

101

tales que H(0) = H0 , y que cumplan las siguientes ecuaciones en D0 (0, T ): ´ d³ (µ H(t), q)0,R3 + d(r(t), q) + c(t, H(t), q) = l(t, q) ∀q ∈ X(R3 ) , dt d(m, H(t)) = (Jd (t), m)0,Ωd

∀m ∈ M(Ωd ) .

Teorema 5.1.3. Bajo las hip´ otesis del Teorema 5.1.1, el Problema 5.1.2 tiene una u ´nica soluci´ on. Demostraci´ on. La segunda ecuaci´on variacional del Problema 5.1.2 significa que rot H(t)|Ωd = Jd (t) en c.t.p. t ∈ (0, T ) , luego H ∈ MJd + W 1 ((0, T ); X0 (R3 ), X00 (R3 )). Entonces, teniendo en cuenta la condici´on inicial junto con la primera ecuaci´on variacional del Problema 5.1.2 para funciones q ∈ X0 (R3 ), deducimos que H coincide con la u ´nica soluci´on del Problema 5.1.1. Para estudiar la existencia y unicidad del multiplicador de Lagrange r asociado al Problema 5.1.2, introducimos G ∈ L2 ((0, T ); X0 (R3 )) definida por Z 0

T

hG(t), q(t)iX0 (R3 )×X(R3 ) dt := =−

Z T³ ´ (µ H(t)−µ H0 , q(t))0,R3 + cˆ(t, H(t), q(t)) − hˆ`(t), q(t)iX0 (R3 )×X(R3 ) dt , 0

para cada q ∈ L2 ((0, T ); X(R3 )). Aqu´ı hemos denotado Z t Z t ˆ cˆ(t, H(t), q(t)) := c(s, H(s), q(t)) ds y `(t) := `(s) ds , 0

0

donde ` ∈ L2 ((0, T ); X0 (R3 )) se caracteriza por Z T Z T h`(t), qiX0 (R3 )×X(R3 ) dt := l(t, q) dt ∀q ∈ X(R3 ) . 0

0

102

Cap´ıtulo 5 Dados q ∈ D((0, T ); X0 (R3 )) y m ∈ L2 ((0, T ); M(Ωd )), consideramos Z T Z t ˇ (t) := ˆ q(s) ds y m(t) := m(s) ds ∀t ∈ [0, T ] . q t

0

Observemos que, aplicando el Lema 2.8.1, Z

Z T³Z

T

ˇ (t)) dt = d(m(t), q 0

0

Z T³Z s ´ ´ d(m(t), q(s)) ds dt = d(m(t), q(s)) dt ds =

T

t

0

Z

T

ˆ d(m(t), q(t)) dt = 0

= 0

0

∀m ∈ L2 ((0, T ); M(Ωd )) .

ˇ ∈ C ∞ ([0, T ]; X0 (R3 )) y, puesto que q ˇ (T ) = 0, la ecuaci´on En consecuencia q (5.1.2) garantiza que Z − 0

Z

T

ˇ (t))0,R3 dt + (µ H(t), ∂t q

T

ˇ (t)) dt = c(t, H(t), q 0

Z ˇ (0))0,R3 + = (µ H0 , q

T

ˇ (t)) dt . l(t, q 0

ˇ y el Lema 2.8.1, reescribimos esta ecuaci´on como Utilizando la definici´on de q Z T³ ´ (µH(t) − µH0 , q(t))0,R3 + cˆ(t, H(t), q(t)) − hˆ`(t), q(t)iX0 (R3 )×X(R3 ) dt = 0 , 0

o equivalentemente, Z 0

T

hG(t), q(t)iX0 (R3 )×X(R3 ) dt = 0 .

Esta propiedad, junto con la densidad de D((0, T ); X0 (R3 )) en L2 ((0, T ); X0 (R3 )), garantiza que G ∈ L2 ((0, T ); X0 (R3 ))0 . Luego, en virtud del Lema 5.1.2, existe una u ´nica funci´on r ∈ L2 ((0, T ); M(Ωd )) con Υ∗ r(t) = G(t) .

(5.1.4)

Concluimos notando que r ∈ L2 ((0, T ); M(Ωd )) satisface la primera ecuaci´on variacional del Problema 5.1.2. En efecto, consideramos q ∈ X(R3 ) y ϕ ∈ D(0, T );

Problema de evoluci´on

103

dϕ tomando el producto de dualidad de (5.1.4) con q(x) (t) ∈ L2 ((0, T ); X(R3 )), dt d d b y utilizando que cˆ(t, H(t), q) = c(t, H(t), q) y ˆ`(t) = `(t), resulta que dt dt h

d d(r(t), q) , ϕ iD0 (0,T )×D(0,T ) = dt d = h (µH(t), q)0,R3 + c(t, H(t), q) + l(t, q) , ϕ iD0 (0,T )×D(0,T ) , dt

que es la primera ecuaci´on variacional del Problema 5.1.2 en el sentido de D0 (0, T ). ¥

5.1.3.

Una formulaci´ on variacional BEM–FEM

Definimos las formas bilineales a0 (q1 , q2 ) := (µq1 , q2 )0,Ω + µ0 hγτ q2 , Vτ (γτ q1 )iτ ,Γ , ³ ´ b2 ) := a0 (q1 , q2 ) + µ0 hη2 , Vη1 i1/2,Γ + hη2 , Bq1 i1/2,Γ − hη1 , Bq2 i1/2,Γ , a(b q1 , q b b1 := (q1 , η1 ) y q b2 := (q2 , η2 ) en X(Ω). para cualesquiera q N´otese que, bajo la hip´otesis (1.1.5) sobre el coeficiente µ, los Teoremas 2.7.5 y 2.7.6 y la propiedad (4.2.6) garantizan que ³ b2 )| ≤ C1 kq1 k0,Ω kq2 k0,Ω + kγτ q1 kH−1/2 (Γ) kγτ q2 kH−1/2 (Γ) |a(b q1 , q ×

×

+ kη1 k−1/2,Γ kη2 k−1/2,Γ + kη2 k−1/2,Γ kγτ q1 kH−1/2 (Γ) ×

´ + kη1 k−1/2,Γ kγτ q2 kH−1/2 (Γ) , (5.1.5) ×

b) ∈ R con y que adem´as a(b q, q ³ b) ≥ C2 kqk20,Ω + kγτ qk2 a(b q, q

−1/2

H×

(Γ)

´ + kηk2−1/2,Γ ,

(5.1.6)

104

Cap´ıtulo 5

b b1 := (q1 , η1 ), q b2 := (q2 , η2 ) y q b = (q, η) en X(Ω). para cualesquiera q b 0 = (H0 |Ω , γ ext H0 ), donde H0 ∈ U(R3 ) representa el campo magn´etico Sea H n inicial. Entonces proponemos la siguiente formulaci´ on variacional BEM–FEM del problema (1.1.1–1.1.3). Problema 5.1.3 Buscar −1/2 b := (H, λ) ∈ L2 ((0, T ); X(Ω)) b H ∩ C 0 ([0, T ]; L2 (Ω) × H0 (Γ)) ,

r ∈ L2 ((0, T ); M(Ωd )) , tales que b b0, H(0) =H y que satisfagan las siguientes ecuaciones en D0 (0, T ): ´ d³ b b) + d(r(t), q) + c(t, H(t), q) = l(t, q) a(H(t), q dt d(m, H(t)) = (Jd (t), m)0,Ωd

b ∀b q ∈ X(Ω) , ∀m ∈ M(Ωd ) .

Teorema 5.1.4. Bajo las hip´ otesis del Teorema 5.1.1, el Problema 5.1.3 tiene una u ´nica soluci´ on. Adem´ as, se relaciona con la soluci´ on H del Problema 5.1.2 ext b a trav´es de la relaci´ on H := (H|Ω , γn H) y el multiplicador de Lagrange r es el mismo para ambos problemas. Demostraci´ on. Sea H ∈ W 1 ((0, T ); X(R3 ), X0 (R3 )), r ∈ L2 ((0, T ); M(Ωd )) la b ∗ := (H∗ , λ), donde H∗ (t) := H(t)|Ω y soluci´on del Problema 5.1.2. Definimos H λ(t) := γnext H(t). Entonces −1/2 b ∗ ∈ L2 ((0, T ); X(Ω)) b H ∩ C 0 ([0, T ]; L2 (Ω) × H0 (Γ)) ,

r ∈ L2 ((0, T ); M(Ωd )) es una soluci´on del Problema 5.1.3. En efecto, utilizando la condici´on inicial del Problema 5.1.2, tenemos que H∗ (0) = H(0)|Ω = H0 |Ω

y λ(0) = γnext H(0) = γnext H0 .

Problema de evoluci´on

105

Adem´as, aplicando la segunda ecuaci´on del Problema 5.1.2, d(m, H∗ (t)) = d(m, H(t)) = (Jd (t), m)0,Ωd

∀m ∈ M(Ωd ) .

b ∗ y r satisface la primera ecuaci´on Por tanto lo que tenemos que estudiar es si H del Problema 5.1.3. Para ello, consideramos el potencial ψ ∈ L2 ((0, T ); W 1 (Ωe )) asociado a la funci´on H(t) ∈ X(R3 ) en c.t.p. t ∈ (0, T ) como en el Lema 4.2.4. An´alogamente, introducimos el potencial ϕ ∈ W 1 (Ωe ) correspondiente a una funci´on gen´erica q ∈ X(R3 ). Recordemos que ψ(t) ∈ W 1 (Ωe ) es arm´onico y que γτ H∗ (t) = rotΓ ψ(t) ,

(5.1.7)

∂ψ en c.t.p. t ∈ (0, T ). En consecuencia, definiendo λ(t) := γnext H(t) = (t) y ∂n utilizando las propiedades (2.7.8), resulta que V(λ(t)) + B(H∗ (t)) = 0

y

1 λ(t) = −N (γψ(t)) + ( I − K∗ )(λ(t)) 2

(5.1.8)

en c.t.p. t ∈ (0, T ). Ahora, integramos por partes en (µ H(t), q)0,R3 = (µ H∗ (t), q)0,Ω + µ0 (∇ψ(t), ∇ϕ)0,Ωe , de modo que (µ H(t), q)0,R3 = (µ H∗ (t), q)0,Ω − µ0 hλ(t), γϕi1/2,Γ . Combinando esta ecuaci´on con la segunda de (5.1.8) y con (2.7.7), (µ H(t), q)0,R3 = (µ H∗ (t), q)0,Ω 1 + µ0 hrotΓ ϕ, Vτ (rotΓ ψ(t))iτ ,Γ − µ0 hλ(t), ( I − K)γϕi1/2,Γ . 2 Aplicando la propiedad (5.1.7), (µ H(t), q)0,R3 = (µ H∗ (t), q)0,Ω + µ0 hγτ q, Vτ (γτ H∗ (t))iτ ,Γ − µ0 hλ(t), Bqi1/2,Γ .

106

Cap´ıtulo 5

A˜ nadimos a esta ecuaci´on la siguiente versi´ on de la primera identidad de (5.1.8): ³ ´ −1/2 µ0 hη, V(λ(t))i1/2,Γ + hη, B(H∗ (t))i1/2,Γ = 0 ∀η ∈ H0 (Γ), y as´ı obtenemos que b ∗ (t), q b b) ∀b (µ H(t), R−1 q)0,R3 = a(H q := (q, η) ∈ X(Ω) , en c.t.p. t ∈ (0, T ). Por lo tanto, la primera ecuaci´on del Problema 5.1.2 garantiza que se cumple la siguiente identidad en D0 ((0, T )): ´ d ³ b∗ b b) + d(q, r(t)) + c(t, H∗ (t), q) = l(t, q) ∀b a(H (t), q q := (q, η) ∈ X(Ω) , dt que es la primera ecuaci´on del Problema 5.1.3. Nos resta estudiar la unicidad de soluci´on del Problema 5.1.3. Con este fin, consideramos −1/2 b ∗ = (H∗ , λ) ∈ L2 ((0, T ); X(Ω)) b H ∩ C 0 ([0, T ]; L2 (Ω) × H0 (Γ)) ,

r ∈ L2 ((0, T ); M(Ωd )) que sea una soluci´on del Problema 5.1.3 homog´eneo. A partir de la primera b b = (0, η) ∈ X(Ω), ecuaci´on del Problema 5.1.3 para funciones de la forma q deducimos que λ(t) = −V −1 B(H∗ (t)) en c.t.p. t ∈ (0, T ). Tambi´en, a partir de b b := (q, 0) ∈ X(Ω) dicha ecuaci´on para funciones de la forma q y eliminando la inc´ognita λ con la expresi´on anterior, tenemos que ´ d d³ (µ H∗ (t), q)0,Ω +µ0 hγτ q, Vτ (γτ H∗ (t))iτ ,Γ + hV −1 B(H∗ (t)), Bqi1/2,Γ + dt dt d + d(r(t), q) + c(t, H∗ (t), q) = 0 , dt en D0 ((0, T )) y para todo q ∈ X(Ω). Denotamos por ψ ∈ L2 ((0, T ); W 1 (Ωe )) el potencial asociado a la funci´on H∗ (t) ∈ X(R3 ) en c.t.p. t ∈ (0, T ) como en el Lema 4.2.4. An´alogamente, ϕ ∈ W 1 (Ωe ) representa el potencial correspondiente

Problema de evoluci´on

107

a una funci´on gen´erica q ∈ X(R3 ). Utilizando que ψ(t) satisface (5.1.7) y (5.1.8), reescribimos la ecuaci´on anterior como ´ d d³ 1 ∂ψ (µH∗ (t), q)0,Ω + µ0 hN (γψ(t)), γϕiΓ − h( I − K∗ )( (t)), γϕi1/2,Γ + dt dt 2 ∂n d + d(r(t), q) + c(t, H∗ (t), q) = 0 , dt en D0 ((0, T )) y para todo q ∈ X(Ω). Aplicando las identidades (5.1.8) e integrando por partes, tenemos que ´ d³ (µH∗ (t), q)0,Ω + µ0 (∇ψ(t), ∇ϕ)0,Ωe + d(r(t), q) + c(t, H∗ (t), q) = 0 , dt en D0 ((0, T )) y para todo q ∈ X(Ω). Es decir, ´ d³ (µ R−1 H∗ (t), q)0,R3 + d(r(t), q) + c(t, R−1 H∗ (t), q) = 0 , dt en D0 ((0, T )) y para todo q ∈ X(Ω). Esta ecuaci´on, junto con la condici´on inicial y la segunda identidad del Problema 5.1.3, garantizan que (R−1 H∗ , r) es soluci´on del Problema 5.1.2 homog´eneo y, en consecuencia, es id´enticamente nula. Por lo tanto, el Problema 5.1.3 homog´eneo s´olo admite la soluci´on trivial; as´ı que, por superposici´on, concluimos que el Problema 5.1.3 admite a lo sumo una soluci´on. ¥ Reescribimos el Problema 5.1.3 como         A −µ0 B∗ D ∗ H(t) C(t) 0 0 H(t) `(t)                d µ0 B  λ(t)  +  0  λ(t)  =  0 , (5.1.9) µ0 V 0  0 0           dt         0 0 0 r(t) D 0 0 r(t) `d (t) donde utilizamos la notaci´on del P´arrafo 4.2.3 y adem´as hC(t)q1 , q2 iX0 (Ω)×X(Ω) := c(t, q1 , q2 ) h`(t), qiX0 (Ω)×X(Ω) := l(t, q) h`d (t), miM0 (Ωd )×M(Ωd ) := (Jd (t), m)0,Ωd

∀q1 , q2 ∈ X(Ω) , ∀q ∈ X(Ω) , ∀m ∈ M(Ωd ) .

108

Cap´ıtulo 5

Dado que λ(t) = −V −1 B(H(t)) en c.t.p. t ∈ (0, T ), podemos eliminar la inc´ognita λ del sistema (5.1.9). As´ı obtenemos el sistema reducido         ∗ −1 ∗ A + µ B V B D H(t) C(t) 0 H(t) `(t) 0 d  +  = . (5.1.10) dt 0 0 r(t) D 0 r(t) `d (t)

5.2.

Esquema semidiscreto

En lo que sigue, suponemos que Ω y Ωc son poliedros de Lipschitz y utilizamos los espacios discretos introducidos en la Secci´on 4.3.

5.2.1.

Descripci´ on y an´ alisis del problema semidiscreto

Proponemos la siguiente versi´on semidiscreta del Problema 5.1.3. Problema 5.2.1 Buscar b h (t) := (Hh (t), λh (t)) ∈ X b h (Ω) H

y

rh (t) ∈ Mh (Ωd )

tales que b h (0) ≈ H b0, H y que satisfagan las siguientes ecuaciones en todo t ∈ (0, T ]: ´ d³ b b) + d(rh (t), q) + c(t, Hh (t), q) = l(t, q) a(Hh (t), q dt d(m, Hh (t)) = (Jd (t), m)0,Ωd

b h (Ω) , ∀b q∈X ∀m ∈ Mh (Ωd ) .

Teorema 5.2.1. El Problema 5.2.1 tiene una u ´nica soluci´ on. Demostraci´ on. Deducimos del Lema 4.3.1 y de [29, Teorema 1.1] que para todo t ∈ [0, T ] existe alg´ un Hdh (t) ∈ Xh (Ω) con d(m, Hdh (t)) = (Jd (t), m)0,Ωd

∀m ∈ Mh (Ωd ) .

Problema de evoluci´on

109

Introducimos H∗h (t) := Hh (t)−Hdh (t) ∈ X0,h (Ω) y reescribimos el Problema 5.2.1 b ∗ (t) := (H∗ (t), λh (t)) como en t´erminos de la nueva inc´ognita H h h d b∗ b ∗ (t), q) = l∗ (t, q) ∀b b) + c(t, H a(Hh (t), q q ∈ X0,h (Ω) × Λh (Γ) , h dt

(5.2.1)

para todo t ∈ (0, T ] y denotando l∗ (t, q) := l(t, q) −

d b) − c(t, Hdh (t), q) . a((Hdh (t), 0), q dt

La ecuaci´on (5.2.1) es equivalente a  Ah + µ0 B∗h Vh−1 Bh  µ0 Vh−1 Bh

     ∂t H∗h (t) C h (t) 0 H∗h (t)  +  = µ0 I ∂t λh (t) 0 0 λh (t)   `∗h (t) , (5.2.2) = −1 d −µ0 Vh Bh ∂t Hh (t) 0

donde utilizamos la notaci´on de la demostraci´on del Teorema 4.3.2 y adem´as hC h (t)q, qh iX0h (Ω)×Xh (Ω) = c(t, q, qh ) , h`∗h (t), qh iX0h (Ω)×Xh (Ω) = l∗ (t, qh ) − µ0 hB∗h Vh−1 Bh ∂t Hdh (t), qh iX0h (Ω)×Xh (Ω) , para cada t ∈ (0, T ] y para cualesquiera q ∈ X(Ω) y qh ∈ Xh (Ω). Recordemos que el operador Ah + B∗h Vh−1 Bh es inversible; v´ease la demostraci´on del Teorema 4.3.2. Entonces la teor´ıa cl´asica de sistemas lineales de ecuaciones diferenciales garantiza que existe una u ´nica funci´on H∗h (t) ∈ X0,h (Ω) solub h (t) soluci´on del Proci´on del sistema (5.2.2). Luego existe una u ´nica funci´on H blema 5.2.1 y ´esta tiene la forma Hh (t) = Hdh (t) + H∗h (t)

y

λh (t) = λh (0) − Vh−1 Bh (Hh (t) − Hh (0)) . (5.2.3)

Para estudiar la existencia y unicidad del multiplicador de Lagrange, introducimos el operador Fh : [0, T ] × Xh (Ω) → R caracterizado por Z t³ ´ b h (t) − H b h (0), (q, 0)) . (5.2.4) Fh (t, q) := l(s, q) − c(t, Hh (s), q) ds − a(H 0

110

Cap´ıtulo 5

En particular, para cada q ∈ X0,h (Ω) tenemos que Fh (0, q) = 0 y que ∂t Fh (t, q) = l(t, q) −

d b a(Hh (t), (q, 0)) − c(t, Hh (t), q) = 0 ∀t ∈ (0, T ) , dt

luego Fh (t, q) = 0 ∀t ∈ [0, T ] . Recordando la definici´on del espacio X0,h (Ω) y aplicando el Lema 4.3.1, concluimos que para cada t ∈ [0, T ] existe un u ´nico rh (t) ∈ Mh (Ωd ) con d(rh (t), q) = Fh (t, q) ∀q ∈ Xh (Ω) .

(5.2.5) ¥

Nota 5.2.1 Como se˜ nalamos en las Notas 3.2.1 y 4.3.1, la inversa del operador tangencial rotΓ est´a involucrada en la definici´on de la forma bilineal a(·, ·) y, en consecuencia, nuestro m´etodo num´erico parece dif´ıcil de implementar y computacionalmente costoso. De nuevo, esto no sucede en la pr´actica porque dicho operador puede eliminarse en el c´alculo efectivo de la matriz del sistema de ecuaciones diferenciales lineales ordinarias asociado al Problema 5.2.1. Al respecto, v´ease el P´arrafo 5.3.

5.2.2.

An´ alisis de la convergencia

En este apartado estudiamos el orden de convergencia del esquema semidiscreto con respecto al Problema 5.1.3. Teorema 5.2.2. Sean H y Hh las primeras componentes de las soluciones de los Problemas 5.1.3 y 5.2.1, respectivamente. Supongamos que H pertenece al espacio C 1 ([0, T ]; Hs (rot, Ω) ∩ X(Ω)) para alg´ un s > 1/2 y que Hh (0) = I h H(0). Entonces tenemos la siguiente estimaci´ on de error: Z T 2 2 sup kH(t)−Hh (t)k0,Ω + sup kγτ (H(t)−Hh (t))k −1/2 + kH−Hh k2rot,Ωc dt ≤ t∈[0,T ]

H×

t∈[0,T ]

h ≤Ch

(Γ)

0

i sup kH(t)k2Hs (rot,Ω) + sup k∂t H(t)k2Hs (rot,Ω) .

t∈[0,T ]

t∈[0,T ]

Problema de evoluci´on

111

Demostraci´ on. A continuaci´ on, seguimos el razonamiento expuesto en [29, Teorema 1.1]; de hecho, esta demostraci´on es an´aloga a la que hemos utilizado para probar el Teorema 4.3.4. En primer lugar, n´otese que el Lema 4.3.1 nos permite introducir una funci´on ph (t) ∈ Xh (Ω) tal que D h ph (t) = D h (H(t) − I h H(t)) ,

(5.2.6)

para todo t ∈ [0, T ]. Entonces Hdh (t) := ph (t) + I h H(t) cumple d(m, Hdh (t)) = (Jd (t), m)0,Ωd

∀m ∈ Mh (Ωd ) ,

para cada t ∈ [0, T ]. Descomponemos H − Hh como H − Hh = e1h − e2h , siendo e1h (t) := H(t) − Hdh (t) ∈ X(Ω) y e2h (t) := Hh (t) − Hdh (t) ∈ X0,h (Ω) , para todo t ∈ [0, T ]. A continuaci´ on, acotamos las normas de los campos e1h (t) y e2h (t) en el espacio H(rot, Ω). Por un lado, escribiendo e1h (t) = H(t)−I h H(t)+I h H(t)−Hdh (t) y aplicando la desigualdad triangular, tenemos que ke1h (t)krot,Ω ≤ (1 +

kdk ) kH(t) − I h H(t)krot,Ω , β

(5.2.7)

para todo t ∈ [0, T ]. Hemos denotado por kdk > 0 la norma de la forma bilineal d(·, ·) y por β > 0 la constante que aparece en el Lema 4.3.1. Adem´as, bajo las hip´otesis de regularidad del teorema, podemos derivar en la ecuaci´on (5.2.6) y as´ı obtenemos que D h (∂t ph (t)) = D h (∂t H(t) − I h ∂t H(t)) ∀t ∈ [0, T ] . Luego con un razonamiento an´alogo al anterior tenemos que ∂t e1h (t) cumple estimaci´on kdk k∂t e1h (t)krot,Ω ≤ (1 + ) k∂t H(t) − I h ∂t H(t)krot,Ω , (5.2.8) β en todo t ∈ [0, T ].

112

Cap´ıtulo 5 Por otro lado, la primera ecuaci´on del sistema (5.1.10) garantiza que a0 (∂t H(t), q) + µ0 hρ(∂t H(t)), BqiΓ + c(t, H(t), q) = l(t, q) .

para cada funci´on q ∈ X0 (Ω) y todo t ∈ (0, T ). An´alogamente, la primera ecuaci´on de (5.2.2) garantiza que a0 (∂t Hh (t), q) + µ0 hρh (∂t Hh (t)), BqiΓ + c(t, Hh (t), q) = l(t, q) , para cada funci´on q ∈ X0,h (Ω) y todo t ∈ (0, T ). Dado que X0,h (Ω) ⊂ X0 (Ω), deducimos de las dos u ´ltimas identidades que a0 (∂t e2h (t), q) + c(t, e2h (t), q) + µ0 hVh−1 Bh ∂t e2h (t), Bh qiΛ0h (Γ)×Λh (Γ) = = a0 (∂t e1h (t), q) + c(t, e1h (t), q)+ + µ0 hρ(∂t H(t)) − ρh (∂t H(t)), BqiΓ + µ0 hρh (∂t e1h (t)), BqiΓ , para cada funci´on q ∈ X0,h (Ω) y todo t ∈ (0, T ). En particular, se cumple para q := e2h (t) ∈ X0,h (Ω): ´ 1d ³ a0 (e2h (t), e2h (t)) + µ0 hB∗h Vh−1 Bh e2h (t), e2h (t)iΛ0h (Γ)×Λh (Γ) + 2 dt + c(t, e2h (t), e2h (t)) = a0 (∂t e1h (t), e2h (t)) + c(t, e1h (t), e2h (t))+ + µ0 hρ(∂t H(t)) − ρh (∂t H(t)), Be2h (t)iΓ + µ0 hρh (∂t e1h (t)), BH∗h (t)iΓ , para cada t ∈ (0, T ). Aplicando las propiedades (5.1.5), (4.2.6) y (4.3.16), ´ 1d ³ a0 (e2h (t), e2h (t)) + µ0 hB∗h Vh−1 Bh e2h (t), e2h (t)iM0h (Ωd )×Mh (Ωd ) + 2 dt ³ + c(t, e2h (t), e2h (t)) ≤ C1 k∂t e1h (t)k0,Ω ke2h (t)k0,Ω + + krot e1h (t)k0,Ωc krot e2h (t)k0,Ωc + kγτ ∂t e1h (t)kH−1/2 (Γ) kγτ e2h (t)kH−1/2 (Γ) + ×

×

´ + kρ(∂t H(t)) − ρh (∂t H(t))k−1/2,Γ kγτ e2h (t)kH−1/2 (Γ) . ×

Problema de evoluci´on

113

Si integramos esta ecuaci´on en (0, t) ⊂ (0, T ) y utilizamos (5.1.6) junto con la no negatividad del operador B∗h Vh−1 Bh , resulta que

ke2h (t)k20,Ω

+

kγτ e2h (t)k2 −1/2 H (Γ) ×

Z T³ ≤ 2 C1

2 + σ1

Z 0

t

krot e2h k20,Ωc ds ≤

k∂t e1h k0,Ω ke2h k0,Ω + kγτ (∂t e1h )kH−1/2 (Γ) kγτ e2h kH−1/2 (Γ) + ×

0

×

´ + kρ(∂t H) − ρh (∂t H)k−1/2,Γ kγτ e2h kH−1/2 (Γ)+krot e1h k0,Ωc krot e2h k0,Ωc ds + ×

+ ke2h (0)k20,Ω + kγτ e2h (0)k2

−1/2

H×

(Γ)

.

Razonando de forma est´andar, obtenemos la siguiente cota: Z sup ke2h (t)k20,Ω t∈[0,T ] Z ≤ C2

T

0

0

−1/2 H× (Γ)

0

T

T

+

³ k∂t e1h k20,Ω + kγτ ∂t e1h k2

³Z + C2

+

sup kγτ e2h (t)k2 −1/2 H× (Γ) t∈[0,T ]

krot e2h k20,Ωc ds ≤

´ + krot e1h k20,Ωc ds+

kρ(∂t H)−ρh (∂t H)k2−1/2,Γ ds+ke2h (0)k20,Ω +kγτ e2h (0)k2

−1/2 H× (Γ)

´ .

Por lo tanto, utilizando la desigualdad triangular y el Teorema 2.3.5, sup kH(t) − Hh (t)k20,Ω + sup kγτ (H(t) − Hh (t))k2

t∈[0,T ]

Z

T

+ 0

Z + 0

krot(H(t) − Hh (t))k20,Ωc dt ≤ C3 T

−1/2

H×

t∈[0,T ]

³

(Γ)

+

sup ke1h k2rot,Ω + sup k∂t e1h k2rot,Ω

t∈[0,T ]

t∈[0,T ]

´ kρ(∂t H) − ρh (∂t H)k2−1/2,Γ dt + ke1h (0)k2rot,Ω + kH(0) − Hh (0)k2rot,Ω .

114

Cap´ıtulo 5

Luego, aplicando (5.2.7) y (5.2.8), deducimos que sup kH(t) − Hh (t)k20,Ω + sup kγτ (H(t) − Hh (t))k2

t∈[0,T ]

Z + 0

T

−1/2

H×

t∈[0,T ]

krot(H − Hh )k20,Ωc dt ≤ C4

³

t∈[0,T ]

+

sup kH(t) − I h H(t)k2rot,Ω +

t∈[0,T ]

Z + sup k∂t H(t) − I h ∂t H(t)k2rot,Ω +

(Γ)

0

T

´ kρ(∂t H) − ρh (∂t H)k2−1/2,Γ dt .

Ahora concluimos utilizando el Lema 4.3.3 y la estimaci´on del error de interpolaci´on del Teorema 3.2.2. ¥

Corolario 5.2.3. Supongamos que se cumplen las hip´ otesis del Teorema 4.3.4 y que adem´ as λ0 = ρh Hh (0). Entonces tenemos la siguiente estimaci´ on de error: sup kλ(t) − λh (t)k2−1/2,Γ ≤ C h

t∈[0,T ]

³

´ sup kHk2Hs (rot,Ω) + sup k∂t Hk2Hs (rot,Ω) .

t∈[0,T ]

t∈[0,T ]

Demostraci´ on. La primera ecuaci´on (5.1.8) garantiza que λ(t) = −ρH(t) en cada t ∈ [0, T ]. Adem´as, la hip´otesis sobre λ0 junto con la propiedad (5.2.3) implican que λh (t) = −ρh Hh (t) en todo t ∈ [0, T ]. Entonces, aplicando la desigualdad triangular, kλ(t) − λh (t)k−1/2,Γ ≤ kρH(t) − ρh H(t)k−1/2,Γ + kρh (H(t) − Hh (t))k−1/2,Γ , para cada t ∈ [0, T ]. N´otese que, utilizando (4.3.16), kρh (H(t) − Hh (t))k−1/2,Γ ≤ C1 kH(t) − Hh (t)krot,Ω , para todo t ∈ [0, T ]. As´ı que el resultado es consecuencia del Lema 4.3.3 y del Teorema 5.2.2. ¥

Problema de evoluci´on

115

Corolario 5.2.4. Supongamos que se cumplen las hip´ otesis del Teorema 5.2.2 2 y del Corolario 5.2.3, y que adem´ as r ∈ L ((0, T ); M(Ωd ) ∩ Hr (Ωd )) para alg´ un r > 1/2. Entonces tenemos la siguiente estimaci´ on del error: Z T kr(t) − rh (t)k20,Ωd dt ≤ 0

³ ≤Ch

Z sup kH(t)k2Hs (rot,Ω) t∈[0,T ]

+ sup k∂t H(t)k2Hs (rot,Ω) t∈[0,T ]

T

+ 0

´ kr(t)k2r,Ωd dt .

Demostraci´ on. Introducimos el operador continuo F : [0, T ] × X(Ω) → R con Z t³ ´ b b F(t, q) := l(s, q) − c(t, h(s), q) ds − a(H(t) − H(0), (q, 0)) , 0

y recordemos que definimos una versi´ on discreta Fh : [0, T ] × Xh (Ω) → R en (5.2.4). N´otese que los multiplicadores r(t) ∈ M(Ωd ) y rh (t) ∈ Mh (Ωd ) satisfacen las ecuaciones d(r(t), q) = F(t, q)

∀q ∈ X(Ω) ,

d(rh (t), q) = Fh (t, q)

∀q ∈ Xh (Ω) ,

en c.t.p. t ∈ (0, T ). Estas propiedades sugieren introducir una funci´on auxiliar e rh (t) ∈ Mh (Ωd ) de modo que d(e rh (t), q) = F(t, q) ∀q ∈ Xh (Ω) , para todo t ∈ [0, T ]. En efecto, el Lema 4.3.1 garantiza que esta funci´on existe y es u ´nica para cada t ∈ [0, T ], porque el operador F(t, ·) ∈ X0 (Ω) se anula sobre X0 (Ω) y, en particular, sobre X0,h (Ω). Descomponemos r − rh como r − rh = e1h + e2h , siendo e1h (t) := r(t) − e rh (t)

y

e2h (t) := e rh (t) − rh (t) ,

en c.t.p. t ∈ (0, T ). A continuaci´ on estudiamos las normas de e1h y de e2h en L2 ((0, T ); L2 (Ωd )). Por un lado, utilizando el Lema 4.3.1 y la desigualdad triangular, ke1h (t)k0,Ωd ≤ (1 +

kdk ) ´ınf kr(t) − sk0,Ωd , β s∈Mh (Ωd )

116

Cap´ıtulo 5

en c.t.p. t ∈ (0, T ). Adem´as, bajo la hip´otesis de regularidad sobre r, tenemos que r(t) ∈ M(Ωd ) ∩ Hr (Ωd ) en c.t.p. t ∈ (0, T ) y entonces la propiedad (4.3.3) y el Lema 4.1.3 implican que RT rot Mh r(t) = I RT h (rot Mr(t)) = I h (r(t)) ∈ Mh (Ωd ) ,

en c.t.p. t ∈ (0, T ). Luego ke1h (t)k0,Ωd ≤ (1 +

1 ) kr(t) − I RT h (r(t))k0,Ωd , β

(5.2.9)

en c.t.p. t ∈ (0, T ). Por otro lado, seg´ un el Lema 4.3.1, ke2h (t)k0,Ωd ≤

1 β

|d(e rh (t) − rh (t), q)| ≤ kqkrot,Ω q∈Xh (Ω), q6=0 sup

≤

1 β

|F(t, q) − Fh (t, q)| , kqkrot,Ω q∈Xh (Ω), q6=0 sup

para todo t ∈ [0, T ]. Observemos que λ |F(t, q) − Fh (t, q)| ≤ C (εH h + εh ) kqkrot,Ω ,

para cualesquiera t ∈ [0, T ], q ∈ Xh (Ω), y donde εH := h

ελh :=

sup kH(t) − Hh (t)k0,Ω + sup kγτ (H(t) − Hh (t))kH−1/2 (Γ) + × t∈[0,T ] ´1/2 ³Z T krot(H(t) − Hh (t)k20,Ωc dt , +

t∈[0,T ]

0

sup kλ(t) − λh (t)k−1/2,Γ .

t∈[0,T ]

As´ı que, aplicando el Teorema 5.2.2 y el Corolario 5.2.3, resulta que kdk |F(t, q) − Fh (t, q)| sup ≤ β q∈Xh (Ω), q6=0 kqkrot,Ω ³ ´ ≤ C2 h1/2 sup (kH(t0 )kHs (rot,Ω) + sup (k∂t H(t0 )kHs (rot,Ω) , (5.2.10)

ke2h (t)k0,Ωd ≤

t0 ∈[0,T ]

t0 ∈[0,T ]

Problema de evoluci´on

117

para cada t ∈ [0, T ]. Finalmente, utilizando la desigualdad triangular para r − rh = e1h + e2h junto con las cotas (5.2.9) y (5.2.10), tenemos que 1 ) kr(t) − I RT h (r(t))k0,Ωd + β ³ ´ + C2 h1/2 sup kH(t0 )kHs (rot,Ω) + sup k∂t H(t0 )kHs (rot,Ω) ,

kr(t) − rh (t)k0,Ωd ≤ (1 +

t0 ∈[0,T ]

t0 ∈[0,T ]

en c.t.p. t ∈ (0, T ). Por tanto concluimos el resultado aplicando la estimaci´on del error de interpolaci´on (4.3.2). ¥

5.3.

Forma matricial del esquema semidiscreto

El objetivo de este p´arrafo es reescribir el Problema semidiscreto 5.2.1 en b h (Ω) y el sistema generador forma matricial. Para ello utilizamos las bases de X de Mh (Ωd ) introducidas en el P´arrafo 4.4. De hecho, seguimos un razonamiento an´alogo al all´ı expuesto, si bien ahora todos los t´erminos son reales y adem´as la matriz C, las inc´ognitas y el segundo miembro dependen del tiempo. M´as concretamente, expresamos las inc´ognitas del Problema 5.2.1 en t´erminos de las bases de los espacios Xh (Ω) y Λh (Γ): Hh (t) =

X

hE (t) qE +

E∈A0h (Ω)

X

hv (t) ∇ϕv

y λh (t) =

y del sistema generador de Mh (Ωd ): X

λF (t) ρF ;

F ∈Fh (Γ)\{F0 }

v∈Vh (Γ)\{v0 }

rh (t) =

X

rE (t) (rot qE )|Ωd .

E∈Ah (Ωd )

118

Cap´ıtulo 5

Consideramos las matrices reales µ qE , qE 0 )0,Ω , µ0 µ := ( qE , ∇ϕv )0,Ω , µ0

AΩ E,E 0 := (

CE,E 0 (t) :=

AΩΓ E,v

AΓv,v0 := (

Rv,v0 := (rotΓ ϕv0 , Vτ (rotΓ ϕv ))0,Γ ,

BF,v

1 c(t, qE , qE 0 ) , µ0

µ ∇ϕv , ∇ϕv0 )0,Ω , µ0 1 := (ρF , ( I − K)γϕv )0,Γ , 2

VF,F 0 := (ρF 0 , VρF )0,Γ , para cualesquiera aristas interiores E, E 0 ∈ A0h (Ω) ; v´ertices sobre la frontera v, v 0 ∈ Vh (Γ) \ {v0 } y caras de la frontera F, F 0 ∈ Fh (Γ) \ {F0 } . Adem´as, introducimos la matriz DE,E 0 := (rot qE , rot qE 0 )0,Ωd , para cualesquiera E ∈ A0h (Ω) y E 0 ∈ Ah (Ωd ). Con esta notaci´on, el Problema 5.2.1 consiste en determinar los vectores hΩ (t) := (hE (t))E∈A0 (Ω) , h

λ(t) := (λF (t))F ∈Fh (Γ)\{F0 } ,

hΓ (t) := (hv (t))v∈Vh (Γ)\{v0 } , r(t) := (rE (t))E∈A0 (Ωd ) , h

soluci´on del sistema de ecuaciones diferenciales     AΩ AΩΓ 0 D hΩ C 0      ΩΓ > Γ    (A ) A +R B> 0 hΣ  0 0 d    +     dt   0 B V 0−λ  0 0     0 0 0 0 r D> 0

ordinarias     0 0 hΩ fΩ         0 0 hΣ  0    =  ,     0 0−λ  0      0 0 r f Ωd

donde el t´ermino independiente viene dado por f Ω (t) :=

1 (l(t, qE ))E∈A0 (Ω) h µ0

y f Ωd(t) := (Jd (t), rot qE )0,Ωd .

(5.3.1)

Ap´ endice A

Integraci´ on num´ erica Tal y como se˜ nalamos en el P´arrafo 3.3.3, a la hora de implementar efectivamente el sistema lineal (3.3.3), necesitamos estimar las integrales involucradas en el ensamblaje de la matriz y del segundo miembro del sistema. Lo mismo sucede con los sistemas (4.4.1) y (5.3.1). En este ap´endice estudiamos las integrales de contorno involucradas en el ensamblaje del sistema (3.3.3). M´as concretamente, las evaluamos de forma exacta o las aproximamos num´ericamente con f´ ormulas de cuadratura que dise˜ namos al efecto. Por un lado, tenemos BF,v

1 = 2

µ

(1, ϕv )0,F (1, ϕv )0,F0 − |F | |F0 |

¶

µ −

(K∗ 1, ϕv )0,F (K∗ 1, ϕv )0,F0 − |F | |F0 |

¶ ,

para cualesquiera cara F ∈ Th (Σ) \ {F0 } y v´ertice v ∈ Vh (Σ) \ {v0 } sobre la frontera. N´otese que podemos calcular (1, ϕv )0,F de forma exacta. Adem´as, podemos descomponer (K∗ 1, ϕv )0,F como ∗

(K 1, ϕv )0,F = (1, Kϕv )0,F =

X

Z µZ

F 0 ∈Th (Σ) F

119

F0

¶ (x − y) · ny dSx ϕv (y) dSy . 4π |x − y|3

120

Ap´endice A Por otro lado, tenemos X

Rv,v0 =

Z Z (rotΣ ϕv0 )|F · (rotΣ ϕv )|F 0

F,F 0 ∈Th (Σ)

F

F0

dSx dSy , 4π |x − y|

para cualesquiera v´ertices v, v 0 ∈ Vh (Σ) \ {v0 } sobre la frontera, y adem´as VF,F 0

1 = |F | |F 0 |

Z Z

Z Z dSx dSy dSx dSy 1 + 2 |F0 | F0 F0 4π |x − y| F F 0 4π |x − y| Z Z Z Z dSx dSy dSx dSy 1 1 − − , 0 |F0 | |F | F0 F 0 4π |x − y| |F0 | |F | F0 F 4π |x − y|

para cualesquiera caras F, F 0 ∈ Th (Σ) \ {F0 } sobre la frontera. En consecuencia, para ensamblar las matrices B, R y V necesitamos aproximar integrales de la forma ¶ Z µZ Z Z (y − x) · ny dSx dSy dSx p(y) dSy y , (A.0.1) 3 |x − y| F F0 F F 0 |x − y| para F, F 0 ∈ Th (Σ) y donde p ∈ P1 (F ). Es importante observar que estas integrales presentan singularidades cuando F y F 0 tienen puntos en com´ un. Por esta raz´on, aproximamos las integrales (A.0.1) con f´ormulas de cuadratura est´andar u ´nicamente cuando las caras F y F 0 est´en suficientemente alejadas. En caso contrario, hacemos una integraci´ on semianal´ıtica; es decir, realizamos las dos etapas siguientes: 1. En primer lugar, estudiamos las integrales interiores de (A.0.1). Si es posible, las calculamos de forma exacta. En caso contrario, eliminamos sus singularidades y entonces las aproximamos con f´ormulas de cuadratura est´andar. 2. En segundo lugar, estudiamos las integrales exteriores de (A.0.1). Como en el paso anterior eliminamos sus posibles singularidades, ya podemos aproximarlas utilizando la f´ormula de Gauss de orden 3 para tri´angulos: Z ¢ |F | ¡ f (z F1 ) + f (z F2 ) + f (z F3 ) ∀f ∈ C(F ) , f (x) dSx ' 3 F

Integraci´on num´erica con los nodos

121   6 z F1    6 z F2     6 zF 3

= 4 v F1 + v F2 + v F3 , = v F1 + 4 v F2 + v F3 , = v F1 + v F2 + 4 v F3 ,

donde v Fi (i = 1, 2, 3) denotan los v´ertices del tri´angulo F .

Aproximaci´ on de las integrales singulares El objetivo de este p´arrafo es resolver la primera etapa del m´etodo de integraci´on semianal´ıtica. Es decir, estudiamos Z Z (y − x) · ny 1 dSx , I1 (y) := dSx y I2 (y) := 3 |y − x| |y − x| F F para cada vector y ∈ R3 y cualquier tri´angulo F ⊂ R3 . Por simplicidad de notaci´on, identificamos los vectores de R3 con los puntos del espacio af´ın asociado. Sea hy := (py − y) · nF ∈ R y definimos py := y + hy nF , que representa la proyecci´on de y sobre el plano que contiene al tri´angulo F ; v´ease la Figura A.1. y nF

F1

F3 py F2

Figura A.1: Proyecci´on py

122

Ap´endice A

N´otese que y = −hy nF + py , luego Z Z hy nF · ny (py − x) · ny I1 (y) = − ¡ ¢3/2 dSx + ¡ ¢ dSx , 2 2 3/2 F h2 F h2 y + |py − x| y + |py − x| Z 1 I2 (y) = ¡ ¢1/2 dSx . 2 F hy + |py − x|2

(A.0.2) (A.0.3)

En [30, 27] se calcula (A.0.3) de forma exacta; damos su resultado al final del presente p´arrafo. Por otro lado, en la literatura no hemos encontrado aproximaciones de (A.0.2). Por esta raz´on, a continuaci´ on proponemos una aproximaci´ on de (A.0.2) siguiendo el razonamiento utilizado en [30, 27] para calcular (A.0.3). Cuando y est´a en el plano que contiene a F , es claro que hy = 0 y que ny · (y − x) = 0, de forma que I1 (y) = 0. Por lo tanto, en lo que sigue, suponemos sin p´erdida de generalidad que y 6= py . Comenzamos estudiando el segundo sumando de (A.0.2). Para ello, descomponemos la integral utilizando los tres tri´angulos F1 , F2 y F3 determinados por py y los v´ertices de F tomados de dos en dos; v´ease la Figura A.1. De esta forma, Z F

Z 3 X (py − x) · ny (py − x) · ny ²m ¡ ¢3/2 dSx = ¡ ¢ dSx , 2 3/2 Fm h 2 h2y + |py − x|2 m=1 y + |py − x|

donde ²m = ±1. En consecuencia, es suficiente aproximar Z (py − x) · ny ¡ ¢ dSx 2 3/2 F h2 y + |py − x| para un tri´angulo F que tiene a py como v´ertice. Descomponemos nuevamente esta integral proyectando py sobre el lado opuesto de F . Es decir, si F tiene por v´ertices py , v 1 y v 2 , entonces definimos q y como la proyecci´ on de py sobre el lado {s v 1 + (1 − s) v 2 ; s ∈ [0, 1]}. Para m = 1, 2, representamos por Fem el tri´angulo de v´ertices py , q y y v m ; v´ease la Figura A.2 a continuaci´ on. Con esta notaci´on, resulta que Z Z 2 X (py − x) · ny (py − x) · ny e ²m ¡ ¡ ¢3/2 dSx = ¢ dSx , 2 2 3/2 F h2 Fem h2 m=1 y + |py − x| y + |py − x|

Integraci´on num´erica

123 v1 qy v2 Fe1

Fe2

py Figura A.2: Descomposici´on de F siendo e ²m = ±1. Por lo tanto, basta aproximar Z (py − x) · ny ¡ ¢ dSx 2 3/2 F h2 + |p − x| y y cuando F es un tri´angulo rect´angulo con v´ertices py , qy y v tales que (py − qy ) · (v − qy ) = 0

y

py − y = hy nF .

(A.0.4)

Tomando coordenadas polares con centro en py , obtenemos que ! Z Z φ ÃZ ρ cos θ (py − x) · ny (py − x) · ny r dr dθ , ¡ ¢ dSx = 2 3/2 (h2 + r2 )3/2 0 F h2 0 y + |py − x| donde denotamos ρ := |qy − py | y φ := ∠(qy − py , v − py ); v´ease la Figura A.3. Descomponemos el vector ny en sus componentes normal y tangencial como b y ; v´ease la Figura A.4. ny = (ny · nF ) nF + n Si tomamos θb := ∠(b ny , q − py ), entonces ny · (py − x) = −|b ny | |py − x| cos(b ny , x − py ) = −|b ny | r cos(θb + θ) , as´ı que Z

(py − x) · ny ny | ¡ ¢3/2 dSx = −|b 2 F hy + |py − x|2

Z φÃZ 0

0

ρ cos θ

! r2 dr cos(θb + θ) dθ . (h2y + r2 )3/2

124

Ap´endice A x2 v

x φ py

θ

r qy

x2

ρ Figura A.3: Coordenadas polares con centro en py

Adem´as, integrando por partes, tenemos que Z 0

ρ cos θ

#¯ " r2 r r ¯¯ dr = − ¡ ¢1/2 + arg sinh |h | ¯¯ (h2y + r2 )3/2 y h2y + r2 r=

. ρ cos θ

Luego Z

(py − x) · ny ¡ ¢ dSx = 2 3/2 F h2 y + |py − x| Z Ã φ

= |b ny |

0

ρ

¡ ¢1/2 h2y cos2 θ + ρ2

ρ − arg sinh cos θ |hy |

! cos(θb + θ) dθ .

N´otese que con este razonamiento hemos eliminado la singularidad del integrando. En consecuencia, podemos aproximar la integral en la variable θ con una f´ormula de cuadratura unidimensional est´andar. Por ejemplo, podemos aplicar la siguiente f´ormula de Gauss: Z 1 8 5 5 f (θ) dθ ' f (−0,7745) + f (0) + f (0,7745) ∀f ∈ C([−1, 1]) . 9 9 9 −1

Integraci´on num´erica

125

ny

(ny · nF ) nF

py

by n

Figura A.4: Descomposici´on del vector ny Haciendo una descomposici´on del tri´angulo F an´aloga a la anterior, reducimos el c´alculo del primer sumando de (A.0.2) al caso de un tri´angulo rect´angulo F de v´ertices py , qy y v que satisfagan (A.0.4). Entonces, tomando coordenadas polares, resulta que Z F

hy

Z

¡ ¢3/2 dSx = h2y + |py − x|2

φ

ÃZ

0

ρ cos θ

0

! hy r dr dθ . (h2y + r2 )3/2

Ahora integramos de forma expl´ıcita en la variable r: Z 0

ρ cos θ

(h2y

donde

hy r hy dr = − ³ ´1/2 + signo(hy ) , 3/2 2 2 +r ) h2y + cosρ 2 θ  h   y signo(hy ) := |hy |   0

si hy 6= 0 , si hy = 0 .

Entonces tenemos que Z F

³ hy dS = signo(h ) φ− x y (h2y + |p − x|2 )3/2

Z 0

φ

´ |hy | dθ . (h2y + ρ2 + ρ2 tan2 θ)1/2

126

Ap´endice A

Para calcular la integral en la variable θ, seguimos la estrategia propuesta en [30] y hacemos el cambio de variable ρ tan u = ¡ ¢1/2 tan θ , 2 hy + ρ2 cuyo jacobiano es

¡ 2 ¢1/2 hy + ρ2 1 + tan2 u dθ = du . ρ 1 + tan2 θ De esta forma, resulta que ¡ ¢1/2 Z φ Z |hy | |hy | υ 1 + tan2 u du , ¡ ¢1/2 dθ = ρ h2 +ρ2 0 0 1+ y h2 + ρ2 + ρ2 tan2 θ tan2 u y

ρ2

ρ tan φ siendo tan υ := tan ¡ ¢1/2 . Reescribimos este resultado como h2y + ρ2 Z Z φ |hy | |hy | υ cos u du . dθ = ³ ´2 ¡ ¢ 1/2 ρ 0 hy 2 0 h2y + ρ2 + ρ2 tan2 θ 1+ sin u ρ

Haciendo el nuevo cambio de variable v = sin u y definiendo ν := sin υ, Z Z φ |hy | |hy | ν |hy | ν 1 . dθ = ³ ´2 dv = arg tan ¡ ¢ 1/2 ρ 0 ρ |h | 0 h2y + ρ2 + ρ2 tan2 θ 1+ y v ρ

Ahora aplicamos la relaci´on ρ tan φ ν = sin υ = ¡ ¢1/2 , h2y + ρ2 + ρ2 tan2 φ y as´ı deducimos que Z φ 0

|hy | |hy | tan φ dθ = arg tan ¡ ¡ ¢ ¢1/2 ; 1/2 h2y + ρ2 + ρ2 tan2 θ h2y + ρ2 + ρ2 tan2 φ

equivalentemente, Z φ 0

|hy | |hy | sin φ ¡ ¢1/2 dθ = arg sin ¡ ¢1/2 . 2 h2y + ρ2 + ρ2 tan θ h2y + ρ2

Integraci´on num´erica

127

Por lo tanto, tenemos que Z

hy

F

Ã

¡ ¢3/2 dSx h2y + |p − x|2

! |hy | sin φ = signo(hy ) φ − arg sin ¡ ¢1/2 . h2y + ρ2

Concluimos este ap´endice se˜ nalando las l´ıneas generales seguidas en [30] para calcular de forma exacta (A.0.3). De nuevo, descomponemos el tri´angulo F y reducimos el c´alculo de (A.0.3) al caso en el que F es un tri´angulo rect´angulo cuyos v´ertices py , qy y v cumplen (A.0.4). En esta situaci´on, en [30] se muestra que Z F

1

Z

¡ ¢1/2 dSx = h2y + |py − x|2

0

φ

ÃZ

ρ cos θ

r

!

¡ ¢1/2 dr dθ = h2y + r2 ¶ Z φµ ρ2 1/2 2 = (hy + ) − hy dθ . cos2 θ 0 0

Entonces, haciendo cambios de variable an´alogos a los anteriores, resulta que Z F

1

¡ ¢1/2 h2y + |py − x|2

¡ 2 ¢1/2 hy + ρ2 + ρ2 tan2 φ + ρ tan φ dSx = ρ lg ¡ ¢1/2 h2y + ρ2 hy sin φ + hy arg sin ¡ ¢1/2 − hy φ . h2y + ρ2 + ρ2 tan2 φ

Bibliograf´ıa [1] A. Alonso, A mathematical justification of the low-frequency heterogeneous time-harmonic Maxwell equations, Math. Models Methods Appl. Sci. 9 (1999), no. 3, 475–489. [2] A. Alonso, A. Valli, An optimal domain decomposition preconditioner for low-frequency time-harmonic Maxwell equations, Math. Comp. 68 (1999), 607–631. [3] A. Alonso, P. Fernandes, A. Valli, The time-harmonic eddy current problem in general domains: solvability via scalar potentials. UTM–611, Dipartimento di Matematica, Universita degli Studi di Trento, 2001. [4] A. Alonso, A. Valli, A domain decomposition approach for heterogeneus time-harmonic Maxwell equations, Comput. Methods Appl. Mech. Engrg. 143 (1997), 97–112. [5] A. Alonso, R. Hiptmair, A. Valli, Mixed finite element approximation of eddy current problems, IMA J. Numer. Anal. 24 (2004), 255–271. ´de ´lec, A justification of eddy currents [6] H. Ammari, A. Buffa, J.-C. Ne model for the Maxwell equations, SIAM J. Appl. Math. 60 (2000), 1805– 1823. [7] C. Amrouche, C. Bernardi, M. Dauge, V. Girault, Vector potentials in three-dimensional nonsmooth domains, Math. Methods Appl. Sci. 21 (1998), 823–864. 129

130

Bibliograf´ıa

´ dez, R. Rodr´ıguez, P. Salgado, A finite element method with [8] A. Bermu Lagrange multipliers for low-frequency harmonic Maxwell equations, SIAM J. Numer. Anal. 40 (2002), 1823–1849. ´ dez, R. Rodr´ıguez, P. Salgado, Modeling and numerical [9] A. Bermu treatment of boundary data in an eddy currents problem, C. R. Acad. Sci. Paris, Serie I, 335 (2002), 633–638. ´rite ´, The TRIFOU code: Solving the 3-D eddy-currents [10] A. Bossavit, J. Ve problem by using H as state variable, IEEE Trans. MAG-19 6 (1983), 2465– 70. [11] A. Bossavit, Two dual formulations of the 3D eddy-currents problem, COMPEL 4 (1985), 103–116. [12] A. Bossavit, A rationale for edge elements in 3D field computations, IEEE Trans. Mag., 24 (1988), 74–79. [13] A. Bossavit, The computation of eddy-currents in dimension 3 by using mixed finite elements and boundary elements in association, Math. Comput. Modelling, 15 (1991), 33–42. ´ [14] A. Bossavit, Electromagn´ etisme, en vue de la mod´elisation, SpringerVerlag, Paris Berlin Heidelberg, 1993. [15] A. Bossavit, Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements, Academic Press, SanDiego, 1997. [16] A. Buffa, P. Ciarlet, Jr., On traces for functional spaces related to Maxwell’s equation. Part I: An integration by parts formula in Lipschitz polyhedra, Math. Meth. Appl. Sci. 24(1) (2001), 9–30. [17] A. Buffa, Hodge decompositions on the boundary of a polyhedron: the multi-connected case, Math. Meth. Model. Appl. Sci. (2001) Vol. 11, No. 9, pp. 1491–1504.

Bibliograf´ıa

131

[18] A. Buffa, Traces for functional spaces related to Maxwell equations: an overview, Proceedings of GAMM-Workshop, Kiel, 2001. [19] A. Buffa, M. Costabel, Ch. Schwab, Boundary element methods for Maxwell’s equations on non-smooth domains, Numer. Math. (2002) Vol. 92, No. 4, pp. 679–710. [20] A. Buffa, M. Costabel, D. Sheen, On traces for H(rot, Ω) in Lipschitz domains, IAN-CNR 1185 (2000), Universidad de Pavia. [21] A. Buffa, R. Hiptmair, Galerkin boundary element methods for electromagnetic scattering. Topics in computational wave propagation, Lect. Notes Comput. Sci. Eng. 31 Springer, Berlin (2003) 83–124. ´zis, An´ [22] H. Bre alisis Funcional, Alianza Universidad Textos, Madrid, 1984. [23] E. Casas, Introducci´ on a las ecuaciones en derivadas parciales, Servicio de Publicaciones, Universidad de Cantabria, 1992. [24] M. Cessenat, Mathematical methods in electromagnetism: linear theory and applications, World Scientific, Singapur, 1996. [25] 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. [26] R. Courant, D. Hilbert, Methods of mathematical physics, Vol. I, Interscience Publishers, Nueva York, 1953. [27] R. Dautray, J.L. Lions, Analyse math´ematique et calcul num´erique pour les sciences et les techniques, Volume 5, Masson, Paris Milan Barcelone, 1988. [28] L. H. Dixon, Jr., Eddy current losses in transformer windings and circuit wiring, Unitrode Coorporation.

132

Bibliograf´ıa

[29] V. Girault, P. A. Raviart, Finite element approximation of the NavierStokes equations: theory and algorithms, Springer, Berlin Heidelberg New York, 1986. [30] J. Giroire, Mise en œvre de methodes d’elements finis de frontiere, D.E.A. d’Analyse Num´erique, Universidad Pierre et Marie Curie, 1987–1988. [31] J.L. Guermond, P.D. Minev, Mixed finite element approximation of an MHD problem involving conducting and isolating regions: the 3D case, Numer. Methods Partial Differential Equations, 19(6) (2003) 709–731. [32] R. Hiptmair, Symmetric coupling for eddy current problems, Technical Report 148, Soderforschungsbereich 382, Universidad de T¨ ubigen, 2001. [33] M. T. Holmberg, Three-dimensional finite element computation of eddy currents in synchronous machines, Technical Report 350, Universidad de G¨otenborg, 1998. ´de ´lec, On the coupling of boundary integral and [34] C. Johnson, J. C. Ne finite element methods, Math. of Comp. 35 (1980), 1063–1079. [35] W. McLean, Strongly elliptic systems and boundary integral equations, Cambridge University Press, 2000. [36] S. Meddahi, An optimal iterative process for the Johnson-Nedelec method of coupling boundary and finite elements, SIAM Journal on Numerical Analysis 35 (1998), 1393–1415. [37] S. Meddahi, F.J. Sayas, A fully discrete BEM-FEM for the exterior Stokes problem in the plane, SIAM Journal on Numerical Analysis 37 (2000), 2082– 2102. ´s, O. Mene ´ndez and P. Pe ´rez, On the coupling [38] S. Meddahi, J. Valde of boundary integral and mixed finite element methods, J. Comput. Appl. Math. 69 (1996), 127–141.

Bibliograf´ıa

133

[39] S. Meddahi, A mixed-FEM and BEM coupling for a two-dimensional eddy current problem, SIAM Numer. Funct. Anal. and Optimiz. 22 (5&6) (2001), 675–696. ´ rquez, A combination of spectral and finite elements [40] S. Meddahi, A. Ma methods for an exterior problem in the plane, Applied Numerical Mathematics 43 (3) (2002), 275–295. [41] S. Meddahi, V. Selgas, A mixed–FEM and BEM coupling for a three– dimensional eddy current problem, Matematical Modelling and Numerical Analysis, 37(2) (2003), 291–318. ´ rquez, V. Selgas, Computing acoustic waves in an [42] S. Meddahi, A. Ma inhomogeneous medium of the plane by a coupling of spectral and finite elements, SIAM J. Numer. Anal., 41 (5) (2003), 1729–1750. ´ rquez, V. Selgas, A new BEM-FEM coupling strat[43] S. Meddahi, A. Ma egy for two dimensional fluid-solid interaction problems, J. Comput. Phys. 199 (1) (2004), 205–220. [44] S. Meddahi, V. Selgas, An H-based FEM-BEM formulation for a time dependent eddy current problem, Technical Report 05–07 (2005), Departamento de Matem´aticas, Universidad de Oviedo. [45] P. Monk, Finite element methods for Maxwell’s equations, Clarendon Press, Oxford, 2003. ´de ´lec, Acoustic and electromagnetic equations. Integral represen[46] J. C. Ne tations for harmonic problems, Springer-Verlag, New York, 2001. ´de ´lec, Mixed finite elements in R3 , Numer. Math. 35 (1980), [47] J. C. Ne 315–341. [48] W. R. Smythe, Static and dynamic electricity, McGraw-Hill, New York, 1950.

134

Bibliograf´ıa

[49] E. P. Stephan, M. Maischak, A posteriori error estimates for fem-bem couplings of three-dimensional electromagnetic problems, Comput. Methods Appl. Mech. Engrg. 194 (2005), no. 2-5, 441–452. [50] M. Teltscher, M. Maischak, E. P. Stephan, A residual error estimator for an electromagnetic fem-bem coupling problem in R3 , Technical Report 52 (2003) Inst. Angew. Mathematik, Universit¨ at Hannover. [51] M. Teltscher, M. Maischak, E. P. Stephan, A hierarchical error estimator for an electromagnetic scattering problem coupling fem and bem in R3 , Technical Report 53 (2003) Inst. Angew. Mathematik, Universit¨ at Hannover. [52] M. Teltscher, M. Maischak, E. P. Stephan, A hierarchical error estimator for an eddy current fem-bem coupling problem, Technical Report 56 (2003) Inst. Angew. Mathematik, Universit¨ at Hannover. [53] F. Treves, Basic linear partial differential equations, Academic Press, New York, 1975. [54] H.A. van der Vorst, J.B.M. Melissen, A Petrov–Galerkin type method for solving Ax = b, where A is symmetric complex, IEEE Transactions on Magnetics, 26(2) (1990), 706-708. [55] H. Whitney, Geometric Integration Theory, Princeton University Press, Princeton, 1957. [56] K. Yosida, Functional Analysis, Springer–Verlag, Berlin Heidelberg, 1980. [57] E. Zeidler, Nonlinear functional analysis and its applications. II/A. Linear monotone operators, Springer–Verlag, New York, 1990. ˇ ´ıˇ [58] A. Zen sek, Nonlinear elliptic and evolution problems and their finite element approximations, Academic Press, London, 1990.