Algoritmo de Mallado Autoadaptativo

Vista previa del material en texto

UNIVERSIDAD POLITÉCNICA DE MADRID
ESCUELA TÉCNICA SUPERIOR DE INGENIEROS DE
TELECOMUNICACIÓN
PROYECTO FIN DE CARRERA
ALGORITMO DE MALLADO AUTOADAPTATIVO PARA LA
RESOLUCIÓN DE PROBLEMAS ELECTROMAGNÉTICOS EN
TRES DIMENSIONES MEDIANTE EL MÉTODO DE LOS
ELEMENTOS FINITOS
Juan Córcoles Ortega
Septiembre de 2004
TÍTULO: ALGORITMO DE MALLADO AUTOADAPTATIVO PARA LA RESOLU-
CIÓN DE PROBLEMAS ELECTROMAGNÉTICOS EN TRES DIMENSIO-
NES MEDIANTE EL MÉTODO DE LOS ELEMENTOS FINITOS
AUTOR: Juan Córcoles Ortega
TUTOR: Sergio Llorente Romano
PONENTE: Magdalena Salazar Palma
DEPARTAMENTO: Señales, Sistemas y Radiocomunicaciones
Miembros del tribunal calificador
PRESIDENTE:
VOCAL:
SECRETARIO:
FECHA DE LECTURA:
CALIFICACIÓN:
UNIVERSIDAD POLITÉCNICA DE MADRID
ESCUELA TÉCNICA SUPERIOR DE INGENIEROS DE
TELECOMUNICACIÓN
DEPARTAMENTO DE SEÑALES, SISTEMAS Y RADIOCOMUNICACIONES
PROYECTO FIN DE CARRERA
ALGORITMO DE MALLADO AUTOADAPTATIVO PARA LA
RESOLUCIÓN DE PROBLEMAS ELECTROMAGNÉTICOS EN
TRES DIMENSIONES MEDIANTE EL MÉTODO DE LOS
ELEMENTOS FINITOS
Autor:
Juan Córcoles Ortega
Tutor:
Sergio Llorente Romano
Ponente:
Magdalena Salazar Palma
Madrid, Septiembre de 2004
RESUMEN
El presente proyecto realiza el análisis electromagnético mediante elementos finitos de es-
tructuras resonantes tridimensionales (eventualmente multidieléctricas) con la ayuda de
un algoritmo de mallado autoadaptativo. Un algoritmo de mallado autoadaptativo es aquél
capaz de primeramente detectar las zonas del dominio bajo estudio donde se comete un
mayor error en la solución aproximada por el método de elementos finitos, y posteriormen-
te refinar dichas regiones, o añadir más grados de libertad a las mismas, de tal manera que
en la siguiente iteración la aplicación del método proporcione una solución más precisa.
Dicho algoritmo se aplica al problema de autovalores/autovectores (frecuencia de resonan-
cia/campo ~E o ~H) generado por las ecuaciones de Maxwell en problemas electromagnéticos
cerrados y se analiza la convergencia del autovalor.
Este proyecto ha sido financiado con una beca de colaboración del Departamento de Seña-
les, Sistemas y Radiocomunicaciones, dentro del Grupo de Microondas y Radar, otorgada
por el Ministerio de Educación.
PALABRAS CLAVE
Elementos Finitos, Mallado Autoadaptativo, Tetraedros curl-conformes, Indicador de Error
Residual, Indicador de Error Zienkiewicz-Zhu, Algoritmos de Refinamiento Local, Conver-
gencia de Autovalores, Cavidades Resonantes, Resonadores Dieléctricos.
A mis padres.
Agradecimientos
El primer agradecimiento, que siento absolutamente necesario, es para mis padres,
Juan e Isabel, y mi hermana Isabel, por haber estado siempre ah́ı, para lo bueno y para
lo malo, por ayudarme, por aguantarme, por enseñarme cosas que jamás se aprenderán de
ningún libro ni en ninguna escuela, y por cien millones de motivos más. Muchas gracias.
Mil gracias, como no, a mis yayos, Simón, Juana y Ascensión, y a mi difunto abuelo
Fidel, que sin duda estaŕıa orgulloso de ver donde he llegado. Sé que todos se alegran por
mı́ más que si fuesen ellos mı́smos.
Y por supuesto, gracias a mis t́ıas Mari y Reme, que nunca han dejado de apoyarme
y darme consejos.
Me gustaŕıa agradecer a Sergio, mi tutor, su ayuda, que sin duda va much́ısimo más
allá de las meras cuestiones técnicas que hayan surgido.
Quiero aprovechar la oportunidad que tengo aqúı de agradecer a Magdalena su
confianza y el hecho de haberme ofrecido la oportunidad de colaborar con el departamento.
Dar las gracias a Luis Emilio, pues en ocasiones, cuando daba por agotadas todas
las opciones, él siempre surǵıa con una respuesta.
Gracias también a todos los integrantes del grupo por el excepcional ambiente de
trabajo, y en especial becarios y doctorandos, por hacer que muchas tardes frente al
ordenador fuesen más amenas de lo que jamás me habŕıa imaginado.
No puedo olvidar aqúı a todos mis compañeros y amigos del C.M.U. San Juan
Evangelista, no sólo por hacer que estos cinco años en Madrid hayan estado repletos de
ilusiones y experiencias, sino por convertir el Johnny en mi segunda casa, y a ellos en mi
segunda familia.
Y por supuesto, a mis amigos y compañeros de siempre, los de Albacete, con los
que compart́ı quizá los mejores años de mi vida, los del colegio y el instituto, y que en
estos años fuera de mi ciudad natal me han hecho comprender el verdadero sentido de la
amistad. Gracias a todos.
Índice general
Índice general I
Índice de figuras V
Índice de tablas IX
1. Introducción 1
1.1. El Método de los Elementos Finitos y su aplicación al Electromagnetismo . 1
1.2. Algoritmos de mallado autoadaptativo y trabajo desarrollado . . . . . . . . 2
2. Análisis de cavidades resonantes mediante el MEF 7
2.1. Formulación mediante elementos finitos . . . . . . . . . . . . . . . . . . . . 7
2.1.1. Problema original, formulación débil y discretización mediante ele-
mentos curl-conformes . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.2. Formulación Normalizada . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2. Convergencia de autovalores en el MEF . . . . . . . . . . . . . . . . . . . . 10
2.3. Postproceso de la solución . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1. Cálculo de ~V (x, y, z) . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.2. Cálculo de ∇× ~V (x, y, z) . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.3. Cálculo de ∇× ν−1∇× ~V (x, y, z) . . . . . . . . . . . . . . . . . . . . 15
2.3.4. Visualización gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3. Indicadores de Error para problemas electromagnéticos de autovalores
en 3D 19
3.1. Errores a posteriori en el MEF . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1. Error Real VS Error Estimado en problemas electromagnéticos . . . 19
3.1.2. Medida del error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.3. Indicador VS Estimador del Error . . . . . . . . . . . . . . . . . . . 20
3.2. Indicador de Error Residual . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.1. Error Residual . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.2. Error Residual para la formulación normalizada . . . . . . . . . . . . 22
3.2.3. Residuo Interior o Volumétrico . . . . . . . . . . . . . . . . . . . . . 22
3.2.4. Residuo Singular en las Caras . . . . . . . . . . . . . . . . . . . . . . 24
3.3. Indicador de Error de tipo Zienkiewicz-Zhu . . . . . . . . . . . . . . . . . . 28
3.3.1. Indicadores de error de tipo Zienkiewicz-Zhu mediante recuperación
por parches con ajuste por mı́nimos cuadrados . . . . . . . . . . . . 28
3.3.2. Indicador basado en ∇× ~V . . . . . . . . . . . . . . . . . . . . . . . 31
3.3.3. Indicador basado en ∇× ν−1∇× ~V . . . . . . . . . . . . . . . . . . 32
i
ii ÍNDICE GENERAL
3.3.4. Indicador basado en ϑ~V . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3.5. Consideración de las condiciones de salto . . . . . . . . . . . . . . . 35
4. Algoritmos de refinamiento local para mallas tetraédricas 39
4.1. Algoritmo no 1 basado en la bisección . . . . . . . . . . . . . . . . . . . . . 40
4.1.1. Descripción del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . 40
4.1.2. Implementación computacional . . . . . . . . . . . . . . . . . . . . . 42
4.2. Algoritmo no 2 basado en la octasección . . . . . . . . . . . . . . . . . . . . 47
4.2.1. Descripción del algoritmo . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.2. Implementación computacional . . . . . . . . . . . . . . . . . . . . . 57
4.2.3. Una solución intuitiva a la posible aparición de tetraedros degenerados 65
4.3. Medidas de calidad en tetraedros y mallados . . . . . . . . . . . . . . . . . . 68
4.3.1. Medidas de calidad de un tetraedro . . . . . . . . . . . . . . . . . . . 69
4.3.2. Medida de calidad de la malla . . . . . . . . .. . . . . . . . . . . . . 70
4.4. Ejemplos y resultados numéricos . . . . . . . . . . . . . . . . . . . . . . . . 71
4.4.1. Validación de la solución a la degeneración de los tetraedros pro-
puesta para el algoritmo no 2 . . . . . . . . . . . . . . . . . . . . . . 71
4.4.2. Dependencia de los algoritmos con la calidad de las mallas iniciales
y coste computacional . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5. Resultados numéricos: aplicación a cavidades y resonadores dieléctricos 85
5.1. Algoritmos de mallado autoadaptativo . . . . . . . . . . . . . . . . . . . . . 86
5.1.1. Especificación de la condición de finalización . . . . . . . . . . . . . 86
5.1.2. Elección de los elementos a refinar . . . . . . . . . . . . . . . . . . . 87
5.1.3. Factores de ponderación del indicador de error residual y parámetro
de refinamiento γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2. Cavidad rectangular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.2.1. Algoritmo autoadaptativo con indicadores de error de tipo residual . 88
5.2.2. Algoritmo autoadaptativo con indicadores de error de tipo
Zienkiewicz-Zhu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.2.2.1. Algoritmo de refinamiento no 1 . . . . . . . . . . . . . . . . 92
5.2.2.2. Algoritmo de refinamiento no 2 . . . . . . . . . . . . . . . . 96
5.3. Cavidad tipo ridge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.3.1. Análisis del primer modo . . . . . . . . . . . . . . . . . . . . . . . . 100
5.3.2. Análisis del tercer modo . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.4. Cavidad con dieléctrico en el centro . . . . . . . . . . . . . . . . . . . . . . 116
5.4.1. Análisis del primer modo . . . . . . . . . . . . . . . . . . . . . . . . 116
5.4.2. Análisis del primer modo superior o segundo modo . . . . . . . . . . 123
5.5. Cavidad rectangular cargada con dieléctrico . . . . . . . . . . . . . . . . . . 130
6. Conclusiones y trabajo futuro 153
A. Elementos tetraédricos curl-conformes de 2o orden subparamétricos 157
A.1. Elementos curl-conformes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
A.2. Definición del elemento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
A.3. El elemento de referencia o canónico . . . . . . . . . . . . . . . . . . . . . . 160
A.4. Transformación geométrica . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
B. Integración de funciones tridimensionales polinómicas en tetraedros 165
ÍNDICE GENERAL iii
C. Glosario de śımbolos y terminoloǵıa 169
Bibliograf́ıa 171
iv ÍNDICE GENERAL
Índice de figuras
1.1. Estructura general de un algoritmo de mallado autoadaptativo . . . . . . . 4
4.1. Analoǵıa 2D mostrando la diferencia entre una malla no conforme y una
conforme. Como se ve, en la primera existe lo que se denomina vértice colgante 39
4.2. Diagrama del algoritmo de refinamiento no 1 . . . . . . . . . . . . . . . . . 40
4.3. Tipos de tetraedros según sus marcas. La arista de refinamiento está se
muestra en rojo con la letra r mientras que las aristas marcadas se muestran
en azul con la letra m . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4. Pasos para la bisección de un tetraedro. La arista de refinamiento está se
muestra en rojo con la letra r para el tetraedro padre y r1 y r2 para ambos
hijos mientras que las aristas marcadas se muestran en azul con la letra m,
m1 y m2 igualmente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.5. Transiciones entre los diferentes tipos de tetraedros producidos por el algo-
ritmo básico de bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.6. Tetraedro con su correspondiente numeración en vértices y aristas . . . . . 44
4.7. Diagrama de flujo del algoritmo de refinamiento no 1 . . . . . . . . . . . . . 54
4.8. Ejemplo que muestra el criterio a cumplir por las caras de un mallado a
refinar por el algoritmo no 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.9. Clases de tetraedros según el número de aristas marcadas que poseen y su
correspondiente división en subtetraedros . . . . . . . . . . . . . . . . . . . 56
4.10. Numeración en las divisiones de cada clase de tetraedro . . . . . . . . . . . 58
4.11. Diagrama de flujo del algoritmo no 2 básico . . . . . . . . . . . . . . . . . . 61
4.12. Analoǵıa en 2D (frente de caras) de la degeneración de los tetraedros. En
gris, los elementos inicialmente seleccionados. . . . . . . . . . . . . . . . . . 65
4.13. Esquema de mallas padre-hijo del algoritmo no 2 . . . . . . . . . . . . . . . 65
4.14. Diagrama de flujo del algoritmo no 2 completo. Las dos flechas a izquierda
◭◭ indican refinamiento siguiendo el algoritmo básico de la figura 4.11 . . . 67
4.15. Analoǵıa en 2D (frente de caras) de la aplicación de la solución del algoritmo
no 2 para evitar la degeneración de los tetraedros. En gris, los elementos
inicialmente seleccionados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.16. Zona de refinamiento elegida en el cubo para la validación de la solución a
la degeneración en el algoritmo no 2 . . . . . . . . . . . . . . . . . . . . . . 71
4.17. Malla inicial para el cubo 10 × 10 × 10. 12 elementos. η(T ) = 0,7472191. . . 72
4.18. Malla final para el cubo 10 × 10 × 10 aplicando el algoritmo no 2 con la
solución para la degeneración. 461620 elementos. η(T ) = 0,1676262. . . . . 72
4.19. Malla final para el cubo 10×10×10 aplicando la versión del algoritmo no 2
sin la solución para la degeneración. 448598 elementos. η(T ) = 0,0080865053. 73
v
vi ÍNDICE DE FIGURAS
4.20. Distribución acumulada de calidades de la malla para las mallas original y
finales con el algoritmo no 2 con y sin solución. . . . . . . . . . . . . . . . . 74
4.21. Evolución de la calidad de la malla η(T ) en función del número de elementos
para el algoritmo no 2 con y sin solución. . . . . . . . . . . . . . . . . . . . 75
4.22. Zonas de refinamiento elegidas para la pirámide degenerada y la pirámide
regular. El valor de r se reduce en un factor de 0.9 para cada iteración del
algoritmo no 2 y cada tres del no 1. . . . . . . . . . . . . . . . . . . . . . . . 77
4.23. Malla final para la pirámide degenerada con el algoritmo no 1. 26975 ele-
mentos. η(T ) = 0,1865965. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.24. Malla final para la pirámide degenerada con el algoritmo no 2. 26766 ele-
mentos. η(T ) = 0,042878415. . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.25. Malla final para la pirámide regular con el algoritmo no 1. 30292 elementos.
η(T ) = 0,4280632. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.26. Malla final para la pirámide regular con el algoritmo no 2. 25838 elementos.
η(T ) = 0,2721259. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.27. Distribución acumulada de calidades de la malla para las mallas original y
finales con el algoritmo no 1 y no 2 de la pirámide degenerada. . . . . . . . 80
4.28. Distribución acumulada de calidades de la malla para las mallas original y
finales con el algoritmo no 1 y no 2 de la pirámide regular. . . . . . . . . . . 81
4.29. Coste computacional de los algoritmos no 1 y no 2 en el caso de la pirámide
degenerada. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.30. Coste computacional de los algoritmos no 1 y no 2 en el caso de la pirámide
regular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.1. Mallado Inicial de la cavidad rectangular. 18 elementos. . . . . . . . . . . . 89
5.2. Convergencia del autovalor k̃2
0 del primer modo de la cavidad rectangular
con mallado uniforme y adaptado mediante el indicador de error residual. . 91
5.3. Convergencia del autovalor k̃2
0 para los dos primeros modos de la cavidad
rectangular con una malla adaptada al primeromediante el indicador de
error residual y el algoritmo de refinamiento no 1. . . . . . . . . . . . . . . 92
5.4. Convergencia del autovalor k̃2
0 para los dos primeros modos de la cavidad
rectangular con una malla adaptada al primero mediante el indicador de
error residual y el algoritmo de refinamiento no 2. . . . . . . . . . . . . . . 93
5.5. Convergencia del autovalor k̃2
0 del primer modo de la cavidad rectangular
con mallado uniforme y adaptado con los diferentes indicadores de error con
el algoritmo de refinamiento no 1. . . . . . . . . . . . . . . . . . . . . . . . . 95
5.6. Convergencia del autovalor k̃2
0 del primer modo de la cavidad rectangular
con mallado uniforme y adaptado con los diferentes indicadores de error con
el algoritmo de refinamiento no 2. . . . . . . . . . . . . . . . . . . . . . . . 99
5.7. Mallado denso para la obtención de valores precisos de números de onda en
la cavidad ridge. 3833 elementos. η(T ) = 0,3822556 . . . . . . . . . . . . . . 100
5.8. Mallado inicial para la cavidad ridge. 90 elementos. η(T ) = 0,3813936 . . . 101
5.9. Convergencia del autovalor k̃2
0 del primer modo de la cavidad ridge según
los distintos métodos autoadaptativos. Formulación ~H. . . . . . . . . . . . . 103
5.10. Mallado inicial para el análisis con formulación ~E ce la cavidad ridge. 171
elementos. η(T ) = 0,3270671 . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.11. Convergencia del autovalor k̃2
0 del primer modo de la cavidad ridge según
el algoritmo de refinamiento usado con indicador residual. Formulación ~E. . 105
ÍNDICE DE FIGURAS vii
5.12. Campos ~E y ~H del primer modo de la cavidad tipo ridge. . . . . . . . . . . 106
5.13. Cortes en X = 0,4, Y = 0,1 y Z = 0,375 para el campo ~E del primer modo
de la cavidad tipo ridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.14. Cortes en X = 0,4, Y = 0,1 y Z = 0,375 para el campo ~H del primer modo
de la cavidad tipo ridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.15. Convergencia del autovalor k̃2
0 para el tercer modo de la cavidad ridge según
los distintos métodos autoadaptativos. . . . . . . . . . . . . . . . . . . . . . 109
5.16. Malla final para el proceso con indicador de error Zienkiewicz-Zhu clase 3
y algoritmo de refinamiento no 1 en cavidad ridge. 404 elementos. η(T ) =
0,3396284 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.17. Malla final para el proceso con indicador de error de tipo residual y algorit-
mo de refinamiento no 2 en cavidad ridge. 465 elementos. η(T ) = 0,2169389 112
5.18. Campos ~E y ~H del tercer modo de la cavidad tipo ridge. . . . . . . . . . . . 113
5.19. Cortes en X = 0,4, Y = 0,1 y Z = 0,375 para el campo ~E del tercer modo
de la cavidad tipo ridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.20. Cortes en X = 0,4, Y = 0,1 y Z = 0,375 para el campo ~H del tercer modo
de la cavidad tipo ridge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
5.21. Vistas de la cavidad con dieléctrico en el centro . . . . . . . . . . . . . . . 116
5.22. Mallado inicial de la cavidad con dieléctrico en el centro . 96 elementos.
η(T ) = 0,4651330 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.23. Campos ~E y ~H del primer modo de la cavidad con dieléctrico en el centro
calculados con sus respectivas formulaciones . . . . . . . . . . . . . . . . . . 120
5.24. Cortes en X = 0,5, Y = 0,5 y Z = 0,5 para el campo ~E del primer modo de
la cavidad con dieléctrico en el centro con formulación ~E. . . . . . . . . . . 121
5.25. Cortes en X = 0,5, Y = 0,5 y Z = 0,5 para el campo ~H del primer modo
de la cavidad con dieléctrico en el centro con formulación ~H. . . . . . . . . 122
5.26. Malla final obtenida tras el proceso autoadaptativo con formulación ~E. Cor-
te en X = 0,5. 1256 elementos. η(T ) = 0,4045892 . . . . . . . . . . . . . . . 124
5.27. Malla final obtenida tras el proceso autoadaptativo con formulación ~H.
Corte en X = 0,5. 1102 elementos. η(T ) = 0,4018410 . . . . . . . . . . . . . 125
5.28. Campos ~E y ~H del segundo modo de la cavidad con dieléctrico en el centro
con Formulación ~H . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.29. Cortes en X = 0,5, Y = 0,5 y Z = 0,5 para el campo ~E del segundo modo
de la cavidad con dieléctrico en el centro con Formulación ~H. . . . . . . . . 127
5.30. Cortes en X = 0,5, Y = 0,5 y Z = 0,5 para el campo ~H del segundo modo
de la cavidad con dieléctrico en el centro con Formulación ~H. . . . . . . . . 128
5.31. Vistas de la cavidad cargada con dieléctrico . . . . . . . . . . . . . . . . . . 130
5.32. Mallado inicial de la cuarta parte de la cavidad cargada con dieléctrico. 46
elementos. η(T ) = 0,4149189 . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.33. Detalle interior del mallado inicial de la cuarta parte de la cavidad cargada
con dieléctrico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.34. k̃2
0 del primer modo de la cavidad cargada con dieléctrico según la formula-
ción para los método autoadaptativos RES Alg. 1, RES Alg. 2, ZZ Alg. 1 y
ZZ Alg. 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
5.35. k̃2
0 del cuarto modo de la cavidad cargada con dieléctrico según la formula-
ción para los método autoadaptativos RES Alg. 1, RES Alg. 2, ZZ Alg. 1 y
ZZ Alg. 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
viii ÍNDICE DE FIGURAS
5.36. Mallado final para el primer modo de la cavidad cargada con dieléctrico con
formulación ~E, indicador de error residual y algoritmo de refinamiento no
2. 485 elementos. η(T ) = 0,2181919 . . . . . . . . . . . . . . . . . . . . . . . 139
5.37. Mallado final para el segundo modo de la cavidad cargada con dieléctri-
co con formulación ~E, indicador de error Zienkiewicz-Zhu y algoritmo de
refinamiento no 2. 313 elementos. η(T ) = 0,2588586 . . . . . . . . . . . . . . 139
5.38. Mallado final para el tercer modo de la cavidad cargada con dieléctrico con
formulación ~H, indicador de error residual y algoritmo de refinamiento no
1. 494 elementos. η(T ) = 0,4004977 . . . . . . . . . . . . . . . . . . . . . . . 140
5.39. Mallado final para el cuarto modo de la cavidad cargada con dieléctrico con
formulación ~H, indicador de error Zienkiewicz-Zhu y algoritmo de refina-
miento no 1. 398 elementos. η(T ) = 0,3840674 . . . . . . . . . . . . . . . . . 140
5.40. Campos ~E y ~H del primer modo en la cuarta parte de la cavidad cargada
con dieléctrico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.41. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~E del primer modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . . . 142
5.42. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~H del primer modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . . . 143
5.43. Campos ~E y ~H del segundo modo en la cuarta parte de la cavidad cargada
con dieléctrico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
5.44. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~E del segundo modo
en la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . 145
5.45. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~H del segundo modo
en la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . 146
5.46. Campos ~E y ~H del tercer modo en la cuarta parte de la cavidad cargada
con dieléctrico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
5.47. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~E del tercer modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . . . 148
5.48. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~H del tercer modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . .. . . . 149
5.49. Campos ~E y ~H del cuarto modo en la cuarta parte de la cavidad cargada
con dieléctrico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
5.50. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~E del cuarto modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . . . 151
5.51. Cortes en X = 0, Y = 0 y Z = 0,875 para el campo ~H del cuarto modo en
la cuarta parte de la cavidad cargada con dieléctrico. . . . . . . . . . . . . . 152
6.1. Analoǵıa en dos dimensiones de la necesidad de adaptación a la geometŕıa
original de un algoritmo de refinamiento . . . . . . . . . . . . . . . . . . . 154
A.1. Tetraedro de Grado 2, elemento de referencia . . . . . . . . . . . . . . . . . 161
A.2. Transformación geométrica entre un tetraedro real y el tetraedro canóncio . 162
B.1. Productos cruzados y no cruzados al elevar al cuadrado un polinomio . . . 167
Índice de tablas
2.1. Correspondencias para la ecuación del doble rotacional . . . . . . . . . . . . 8
3.1. Vértices que definen las caras de un tetraedro (numeración local) . . . . . . 25
3.2. Aristas que definen las caras de un tetraedro (numeración local) . . . . . . 25
3.3. Coordenadas del punto medio de cada arista para la cuadratura por Gauss-
Legendre en una cara de un tetraedro . . . . . . . . . . . . . . . . . . . . . 27
3.4. Validez de los indicadores de tipo Zienkiewicz-Zhu según los casos dados
por ϑ y ν. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.1. Definición de las aristas del tetraedro en función de los vértices . . . . . . . 44
4.2. Definición de las caras del tetraedro en función de las aristas (aristas copla-
nares) y los vértices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.3. Definición de vértices para los tetraedros hijos según la arista de refinamien-
to del padre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.4. Caras heredada, cortadas y nueva de los tetraedros hijo según la arista de
refinamiento del padre definiendo los vértices en los hijos tal como muestra
la tabla 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 1 . 48
4.6. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 2 . 49
4.7. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 3 . 50
4.8. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 4 . 51
4.9. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 5 . 52
4.10. Vértices, condiciones de contorno y marcas para los tetraedros hijos proce-
dentes de la bisección cuando la arista de refinamiento del padre r(τ) = 6 . 53
4.11. Vértices y condiciones de contorno para los cuatro primeros hijos (hijos
exteriores) de un tetraedro primario de aristas base α3 y α5 . . . . . . . . . 59
4.12. Vértices y condiciones de contorno para los cuatro segundos hijos (hijos
interiores) de un tetraedro primario de aristas base α3 y α5 . . . . . . . . . 60
4.13. Vértices y condiciones de contorno para los cuatro hijos de un tetraedro
secundario tipo A cuando sus aristas marcadas son α1, α2 y α3 . . . . . . . 62
4.14. Vértices y condiciones de contorno para los cuatro hijos de un tetraedro
secundario tipo B cuando sus aristas marcadas son α2 y α4 . . . . . . . . . 63
ix
x ÍNDICE DE TABLAS
4.15. Rotaciones necesarias en los tetraedros secundarios de tipo A para conseguir
que sus aristas marcadas sean α1,α2 y α3 . . . . . . . . . . . . . . . . . . . 64
4.16. Rotaciones necesarias en los tetraedros secundarios de tipo B para conseguir
que sus aristas marcadas sean α2 y α4 . . . . . . . . . . . . . . . . . . . . . 64
4.17. Calidades y proporciones en las mallas originales y obtenidas con ambos
algoritmos para las pirámides degenerada y regular. T0: malla inicial. T1:
malla final algoritmo no 1. T2: malla final algoritmo no 2. . . . . . . . . . . 80
5.1. Resultados para el primer modo de la cavidad rectangular con indicador de
error residual y algoritmo de refinamiento no 1. Valor exacto del número de
onda k0 = 3,088572 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2. Resultados para el primer modo de la cavidad rectangular con indicador de
error residual y algoritmo de refinamiento no 2. Valor exacto del número de
onda k0 = 3,088572 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.3. Indicadores de error residual para el primer modo de la cavidad rectangular
usando el algoritmo de refinamiento no 1. . . . . . . . . . . . . . . . . . . . 90
5.4. Indicadores de error residual para el primer modo de la cavidad rectangular
usando el algoritmo de refinamiento no 2. . . . . . . . . . . . . . . . . . . . 90
5.5. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 1 y algoritmo de refinamiento no 1. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 93
5.6. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 2 y algoritmo de refinamiento no 1. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 94
5.7. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 3 y algoritmo de refinamiento no 1. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 94
5.8. Indicadores de error Zienkiewicz-Zhu clase 1 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 1. γ = 0,5. . . . 96
5.9. Indicadores de error Zienkiewicz-Zhu clase 2 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 1. γ = 0,5. . . . 97
5.10. Indicadores de error Zienkiewicz-Zhu clase 3 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 1. γ = 0,5. . . . 97
5.11. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 1 y algoritmo de refinamiento no 2. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 97
5.12. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 2 y algoritmo de refinamiento no 2. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 98
5.13. Resultados para el primer modo de la cavidad rectangular con indicador de
error Zienkiewicz-Zhu clase 3 y algoritmo de refinamiento no 2. γ = 0,5.
Valor exacto del número de onda k0 = 3,088572 . . . . . . . . . . . . . . . . 98
5.14. Indicadores de error Zienkiewicz-Zhu clase 1 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 2. γ = 0,5. . . . 98
5.15. Indicadores de error Zienkiewicz-Zhu clase 2 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 2. γ = 0,5. . . . 99
5.16. Indicadores de error Zienkiewicz-Zhu clase 3 para el primer modo de la
cavidad rectangular usando el algoritmo de refinamiento no 2. γ = 0,5. . . . 99
ÍNDICE DE TABLAS xi
5.17. Número de onda k̃0 del primer modo de la cavidad ridge para los distintos
procesos autoadaptativos. Formulación ~H. . . . . . . . . . . . . . . . . . . . 102
5.18. Números de onda k̃0 del primer modo de la cavidad ridge para ambos al-
goritmos de refinamiento con el indicador de error residual. Formulación
~E . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . 105
5.19. Número de onda k̃0 del tercer modo de la cavidad ridge para los distintos
procesos autoadaptativos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
5.20. Números de onda k̃0 para el proceso autoadaptativo del primer modo de
la cavidad con dieléctrico en el centro según ambas formulaciones. Valor
obtenido con Ansoft HFSS c©: k0 = 2,747427 . . . . . . . . . . . . . . . . . . 118
5.21. Indicadores de error residual para el primer modo de la cavidad con dieléc-
trico en el centro usando el algoritmo de refinamiento no 1. . . . . . . . . . 119
5.22. Valor del número de onda del primer modo de la cavidad con dieléctrico en
el centro obtenido en cada iteración según la formulación y error relativo
de la media de los mismos k̃0 =
(
k̃
( ~H)
0 + k̃
(~E)
0
)
/2 con respecto al valor
proporcionado por Ansoft HFSS c©: k0 = 2,747427 . . . . . . . . . . . . . . 119
5.23. Números de onda k̃0 para el proceso autoadaptativo del segundo modo de
la cavidad con dieléctrico en el centro no 1. Formulación ~H. Valor obtenido
con Ansoft HFSS c©: k0 = 2,859478 . . . . . . . . . . . . . . . . . . . . . . . 129
5.24. Combinaciones de pared magnética y eléctrica para los plano de simetŕıa
de la cavidad cargada con dieléctrico . . . . . . . . . . . . . . . . . . . . . . 130
5.25. Número de onda calculado k̃0 para los diferentes procesos autoadaptativos
del primer modo de la cavidad cargada con dieléctrico . . . . . . . . . . . . 133
5.26. Número de onda calculado k̃0 para los diferentes procesos autoadaptativos
del segundo modo de la cavidad cargada con dieléctrico . . . . . . . . . . . 134
5.27. Número de onda calculado k̃0 para los diferentes procesos autoadaptativos
del tercer modo de la cavidad cargada con dieléctrico . . . . . . . . . . . . . 135
5.28. Número de onda calculado k̃0 para los diferentes procesos autoadaptativos
del cuarto modo de la cavidad cargada con dieléctrico . . . . . . . . . . . . 136
5.29. Número de onda k̃0 del primer y segundo modo para ambas formulaciones
según el método adaptativo seguido . . . . . . . . . . . . . . . . . . . . . . 138
xii ÍNDICE DE TABLAS
Caṕıtulo 1
Introducción
1.1. El Método de los Elementos Finitos y su aplicación al
Electromagnetismo
Los sistemas de radiocomunicaciones que soportan las nuevas y futuras aplicaciones
multimedia (sistemas de acceso radio, sistemas y redes por satélite , redes móviles de banda
ancha, etc.) se están desarrollando en la parte alta de la banda de microondas y en la banda
de onda milimétricas. La razón fundamental es que dichas aplicaciones requieren anchos de
banda elevados, por lo que se opta por frecuencias portadoras de valor alto donde además
el espectro electromagnético está menos explotado. El aumento de la frecuencia de trabajo
implica sin embargo un diseño del subsistema de RF más complejo y costoso.
En la parte baja de la banda de microondas las herramientas de diseño asistido por
ordenador (CAD) habituales son simuladores circuitales que utilizan modelos equivalentes
de los diferentes componentes. Sin embargo los modelos simples usados por dichos simu-
ladores impiden que esas herramientas puedan utilizarse a frecuencias altas. De hecho,
la complejidad de los dispositivos de alta frecuencia requiere herramientas de diseño que
efectúen un análisis electromagnético de onda completa que permita caracterizar todos los
efectos presentes en la estructura para aśı evitar técnicas de diseño basadas en “prueba y
error”.
Actualmente, un solo dispositivo puede llegar a presentar una estructura intrinca-
da que conste de varios conductores, dieléctricos e incluso semiconductores de tamaño
arbitrario y caracteŕısticas f́ısicas complejas. La cada vez mayor variedad de geometŕıas
y materiales de dichas estructuras, descartan, en general, el uso exclusivo de técnicas de
cálculo anaĺıticas o semianaĺıticas, tales como Transformación Conforme, Separación de
Variables, Ajuste Modal, etc. Por tanto, en general, se hace preciso el empleo de un método
numérico, bien en frecuencia o en tiempo, y sobre problemas con ecuaciones diferenciales,
integrales o integro-diferenciales, como el Método de las Diferencias Finitas, el Método de
los Momentos, el Método de Elementos Finitos, etc . . .
El Método de los Elementos Finitos (en adelante, MEF ) data de los primeros años
40, y aunque inicialmente fue concebido para el campo del cálculo de estructuras, pronto
se llegó a la evidencia de que era una técnica de uso general para la resolución numérica de
problemas formulados mediante ecuaciones en derivadas parciales. En el campo del electro-
magnetismo, y más en concreto en lo concerniente a microondas y ondas milimétricas, el
método de los elementos finitos no fue utilizado hasta finales de los años 60; desde entonces
su popularidad en la resolución de problemas en esta área ha crecido considerablemente
1
2 Introducción
creando una especial actividad investigadora al respecto en los últimos años.
Es preciso mencionar que en esta memoria no se ha incluido ningún caṕıtulo ni
apéndice espećıfico sobre el propio MEF. Se supone de antemano que el lector posee una
cierta familiaridad en cuanto a la mecánica de aplicación del método (generación del
mallado, discretización mediante funciones de base locales a cada elemento, ensambla-
do, etc. . . ). En cualquier caso, para consultar dichas generalidades existen multitud de
libros entre los cuales [Dhat and Touzot, 1981, Reddy, 1993] se suelen recomendar para
una primera lectura sobre el MEF. Entre los libros que tratan el MEF en el contexto
de su aplicación a problemas electromagnéticos de una forma genércia se pueden citar
[Jin, 1993, Volakis et al., 1998].
De dichas generalidades, las más relevantes con respecto a este trabajo son las concer-
nientes a la transformación geométrica y la numeración local/global de elementos, nodos,
etc. . . En la primera, se pasa del sistema real (x, y, z) a un sistema de coordenadas deno-
minado canónico (ξ, η, ζ), de tal forma que los vértices de un tetraedro real (perteneciente
al primer sistema) sean (0, 0, 0), (1, 0, 0), (0, 1, 0) y (0, 0, 1) en el sistema canónico. Hay que
destacar que todas las caras en un elemento finito, vistas desde el exterior, han de estar
descritas por sus vértices en sentido anti-horario; de esta manera se asegura que al calcular
el volumen del elemento, éste siempre sea positivo. Con respecto a la numeración se ha de
discernir entre la de tipo local, a nivel de elemento, y la de tipo global, a nivel de malla.
Este criterio se puede aplicar a cualquier parte del elemento o la malla, como aristas, caras,
vértices, etc. . . Por ejemplo, el vértice no 1 (numeración local) de un tetraedro puede ser el
vértice no 1250 (numeración global) de la malla. En el caso de nodos y grados de libertad
esta numeración es de suma importancia a la hora de realizar el ensamblado para obte-
ner matrices dispersas; diversos métodos de renumeración se emplean para opitmizar esta
cuestión, habiéndose escogido aqúı el método de Gibbs (proporcionado por MODULEF1).
1.2. Algoritmos de mallado autoadaptativo y trabajo desa-
rrollado
Una de las caracteŕısticas más potentes del MEF reside en el hecho de que la con-
tribución de cada elemento es calculada de forma independiente y después ensamblada
para formar la solución total. Aśı pues, una alta densidad de grados de libertad puede ser
usada en las regiones donde la variación del campo sea abrupta, mientras que un menor
número de los mismos puede ser asignado a regiones con una variación del campo relati-
vamente pequeña. Diversas técnicas basadas en diferentes indicadores o estimadores del
error de aproximación se usan para decidir que partes del dominio han de ser refinadas, y
aplicando recursivamente un algoritmo de refinamiento seleccionado se obtiene un mallado
adaptado al comportamiento del campo en la estructura. Éstaes la esencia del mallado
autoadaptativo.
La estructura básica de un algoritmo o método de mallado autoadaptativo se mues-
tra en la figura 1.1. Primeramente se genera una malla que en principio debeŕıa constar
de pocos grados de libertad, a partir de la cual se resuelve el problema mediante el MEF.
Posteriormente se realiza una indicación o estimación del error cometido, y se comprueba
si se ha alcanzado cierta condición de finalización, que puede ser del tipo estimación por
debajo de un umbral, no de iteraciones máximo alcanzado, etc. . . Si es aśı, se procede al
postproceso de la solución, que generalmente suele venir acompañado de alguna visua-
1al final del caṕıtulo se explicará qué es MODULEF
1.2 Algoritmos de mallado autoadaptativo y trabajo desarrollado 3
lización gráfica. Si la condición no se cumple, se ha de seleccionar, atendiendo a algún
criterio, los elementos del mallado cuyo indicador de error sea mayor para su posterior
refinamiento. Tras este refinamiento se obtiene una nueva malla (que dependiendo del re-
finamiento empleado tendrá bien un número distinto de elementos, ó nuevos grados de
libertad añadidos, etc. . . con respecto a la de la iteración anterior), que vuelve a ser usada
para resolver el problema mediante el MEF.
El trabajo desarrollado en este proyecto incluye las partes que no se muestran en
ĺınea discontinua en la figura 1.1.
Las mallas usadas contienen tan sólo elementos de tipo tetraédrico. Este tipo de
elementos se prefiere, entre otros muchos motivos, por su capacidad de adaptación a geo-
metŕıas complejas.
Los indicadores de error aqúı desarrollados responden a dos tipos (posiblemente
los más ampliamente usados por los investigadores): uno de corte más clásico que tra-
dicionalmente ha sido usado en el MEF, denominado residual, y otro desarrollado más
recientemente conocido por los nombres de sus precursores, Zienkiewicz-Zhu. Aunque en
la figura 1.1 no se especifica, según se han implementado, para el cálculo de los mismos
se requiere una parte de lo que se podŕıa denominar postproceso, aunque de una manera
muy básica. Lógicamente, en la literatura se encuentran una gran variedad de tipos de
indicadores de error.
Cuando se han seleccionado los tetraedros con mayor indicador de error, se puede
proceder a varios tipos de refinamiento; entre los más conocidos se encuentran:
Refinamiento tipo h Se dividen los elementos seleccionados en dos o más nuevos ele-
mentos.
Refinamiento tipo p Se incrementa el orden de las funciones base de interpolación en
los elementos seleccionados.
Refinamiento tipo r Se recolocan los nodos de la malla en una mejor distribución de
los mismos.
Técnicas h́ıbridas r − h, h − p Se trata de combinaciones de los algoritmos anteriores.
En este proyecto se llevará a cabo un refinamiento tipo h, por lo que los tetraedros se-
leccionados que formen parte de la malla en una determinada iteración se dividirán en
tetraedros más pequeños que formarán parte de la malla de la siguiente iteración; de esta
manera se añaden más grados de libertad al problema allá donde lo necesite.
Para comprobar la eficacia de las técnicas de refinamiento y de indicación del error,
éstas se aplicarán al análisis mediante el MEF de cavidades resonantes, donde se espera
obtener un grado de refinamiento mayor en las zonas de campo con variación más abrupta
(y eventualmente, singularidades del campo).
Dentro de la parte concerniente al postproceso se ha desarrollado un módulo capaz
de dibujar los campos de mayor interés que presentan las estructuras analizadas, tanto en
forma de vector (con sentido y magnitud), como en módulo, y tanto para la geometŕıa
completa en 3D como para cortes de los planos más interesantes.
Sin duda el lector encontrará interesante el glosario que se encuentra en el apéndice C,
pues posiblemente en algún momento de la lectura tenga que recurrir a él para descifrar
el significado de algún śımbolo que se haya quedado en el tintero.
La mayoŕıa del trabajo se ha implementado haciendo uso de las funciones y rutinas
proporcionadas por MODULEF ([INRIA, 1999]), que es una libreŕıa de elementos finitos
4 Introducción
Figura 1.1: Estructura general de un algoritmo de mallado autoadaptativo
1.2 Algoritmos de mallado autoadaptativo y trabajo desarrollado 5
de propósito general creada por el instituto de investigación francés INRIA. Dos ejem-
plos básicos de su uso que se emplean a modo de gúıa de iniciación se encuentran en
[Arnold and Narasimhan, 1994a, Arnold and Narasimhan, 1994b]. El lenguaje de progra-
mación empleado por la misma y, por ende, en la implementación de los distintos módulos
es Fortran ([Intel, 2004]).
Además, para la creación inicial de mallados se habrá de emplear el software GID c©,
desarrollado en el instituto internacional de métodos numéricos CIMNE, aśı como una
modificación del mismo (Gid2Modulef, desarrollada por el departamento de Teoŕıa de la
Señal de la Universidad de Alcalá de Henares) para adaptar los mallados resultantes a las
estructuras de datos que maneja la biblioteca MODULEF. Para el módulo de dibujo se ha
empleado la potencia gráfica proporcionada por MATLAB c©([Mathworks, 2001]).
6 Introducción
Caṕıtulo 2
Análisis de cavidades resonantes
mediante el MEF
2.1. Formulación mediante elementos finitos
En esta sección se resumen los pasos seguidos para resolver el problema de las cavi-
dades resonantes mediante el método de los elementos finitos. No se pretende aqúı hacer
un desarrollo exhaustivo de todos, pues muchos no son fáciles y llevan un desarrollo más
intrincado que se ha obviado; para un mayor detalle consultar [Salazar-Palma et al., 1998,
Garćıa-Castillo, 1998].
Es importante remarcar que en el tipo de problemas abordados los tensores de
permitividad y permeabilidad se consideran reales y simétricos. Con este tipo de tensores,
las componentes del campo están en fase, por lo que la solución del MEF puede elegirse
como real sin pérdida de generalidad. Aún aśı, se hace hincapié en que el método de los
elementos finitos podŕıa resolver problemas con dichos tensores complejos y no simétricos.
2.1.1. Problema original, formulación débil y discretización mediante
elementos curl-conformes
Dado un dominio Ω con una frontera Γ las ecuaciones de Maxwell expresadas me-
diante la ecuación del doble rotacional tienen la forma del problema de valores propios
siguiente:
∇× (f−1∇× ~W ) = ω2p ~W en Ω
n̂ × ~W = ~0 en ΓDirichlet
n̂ × f−1∇× ~W en ΓNeumman
(2.1)
Las expresiones involucradas en la ecuación 2.1 toman su valor de la tabla 2.1 según
la formulación usada en el problema.
Para una aproximación de la solución ~̃W , tenemos que que el residuo no será exac-
tamente nulo:
R( ~̃W ) = ∇× (f−1∇× ~̃W ) − ω2p ~̃W ≈ 0 (2.2)
donde ~̃W consiste en una combinación lineal de un número dado (número de grados de
libertad Nd) de funciones vectoriales (funciones base):
7
8 Análisis de cavidades resonantes mediante el MEF
Tabla 2.1: Correspondencias para la ecuación del doble rotacional
Formulación ~E Formulación ~H
~W ~E ~H
f µ ε
p ε µ
ΓDirichlet ΓElectrica ΓMagnetica
ΓNeumman ΓMagnetica ΓElectrica
~̃W =
Nd
∑
i=1
~Nidi (2.3)
Teniendo en cuenta que la mayoŕıa de las veces se opera a nivel de elemento, obten-
dŕıamos, donde el ı́ndice e se refiere a las mismas cantidades anteriores en un sólo elemento,
y Ne es el número total de los mismos en la malla:
~̃W =
Ne
∑
e=1
Ne
d
∑
i=1
~N e
i de
i (2.4)
Multiplicando 2.2 por una función test ~̃W ′ (que ha de cumplir ciertas especificacio-
nes contempladas en [Garćıa-Castillo, 1998]), integrando e igualando a 0 por el Teorema
de Mejor Aproximación (consultar, por ejemplo, [Debnath and Mikusinski, 1990]), obte-
nemos:
< ~̃W ′;R( ~̃W ) >=
∫∫∫
Ω
(
~̃W ′ · ∇ × (f−1∇× ~̃W ) − ω2 ~̃W ′ · p ~̃W
)
dΩ = 0 (2.5)
Teniendo en cuenta ahora las condiciones de contorno, aśı como identidades en
el cálculovectorial como el teorema de Green (consultar [Marsden and Tromba, 1991,
Salas and Hille, 1986]), la ecuación 2.5 puede reescribirse como:
< ~̃W ′;R( ~̃W ) >=
∫∫∫
Ω
(
(∇× ~̃W ′) · (f−1∇× ~̃W ) − ω2 ~̃W ′ · p ~̃W
)
dΩ = 0 (2.6)
La ecuación 2.6 es la formulación débil del problema original (ecuación 2.1).
El siguiente paso es la discretización para obtener un sistema de ecuaciones alge-
braico. El vector ~̃W debe cumplir 2.6 para un número finito M de funciones test ~̃W ′
independientes. Normalmente, este número es igual al de grados de libertad, M = Nd. Si
usamos el método de Galerkin donde las funciones test se eligen iguales a las funciones
vectoriales base ~Ni, obtenemos el siguiente sistema de ecuaciones:
< ~Ni;R( ~̃W ) >=
∫∫∫
Ω
(
(∇× ~Ni) · (f−1∇× ~̃W ) − ω2 ~Ni · p ~̃W
)
dΩ = 0 i = 1, . . . , Nd
(2.7)
2.1 Formulación mediante elementos finitos 9
Si la ecuación 2.7 se aplica sobre el elemento e de la malla y se tiene en cuenta 2.4
obtenemos:
{W e} = ([ke] − ω2[me]){de} (2.8)
ke
ij =
∫∫∫
Ωe
(∇× ~N e
i ) · (fe
−1∇× ~N e
j )dΩ, i, j = 1, . . . , N e
d (2.9)
me
ij =
∫∫∫
Ωe
~N e
i pe
~N e
j dΩ, i, j = 1, . . . , N e
d (2.10)
donde N e
d son los grados de libertad del elemento e, ~N e
i las funciones base de dicho ele-
mento, {de} el vector que agrupa los mencionados grados de libertad en el elemento, y fe
y pe son los tensores permeabilidad y permitividad (según formulación) en el elemento e
Realizando el proceso de ensamblado de las matrices [ke] y [me] en [K] y [M ] repec-
tivamente, el sistema final puede escribirse como:
([K] − ω2[M ]){D} = 0 (2.11)
En terminoloǵıa del MEF, [K] es la matŕız de rigidez global, [M ] la matriz de masa
global, y {D} el vector que agrupa los grados de libertad globales del problema.
2.1.2. Formulación Normalizada
Si realizamos la conversión
~Wn =
√
(p0) ~W (2.12)
teniendo en cuenta que
f = f0fr (2.13)
p = p0pr (2.14)
y, finalmente llevando a cabo un cambio de notación para mayor sencillez:
~V = ~̃Wn (2.15)
ν = fr (2.16)
ϑ = pr (2.17)
la ecuación 2.6 puede expresarse como
< ~V ′;R(~V ) >=
∫∫∫
Ω
(
(∇× ~V ′) · (ν−1∇× ~V ) − k2
0
~V ′ · ϑ~V
)
dΩ = 0 (2.18)
donde ahora el autovalor del problema ha pasado a ser k2
0 en vez de ω2. Siguiendo
los mismos pasos anteriores, se llega a un sistema final:
([K] − k2
0[M ]){D} = 0 (2.19)
10 Análisis de cavidades resonantes mediante el MEF
que es un sistema de autovalores / autovectores que tendrá que ser resuelto. En este
caso, el autovalor k2
0 si se corresponde con el autovalor del problema original, mientras que
el autovector {D} corresponde al vector de los grados de libertad (de los que algunos serán
forzados a tener valor nulo para cumplir condiciones de contorno) que formarán parte de
la interpolación para hallar el autovector del problema ~V . En [Garćıa-Castillo, 1998] se
aborda este tema con mayor detalle.
2.2. Convergencia de autovalores en el MEF
Dado que el problema que se va a analizar constituye inherentemente un problema
de autovalores, es interesante comprobar la convergencia del error a priori que proporciona
el método de los elementos finitos para este tipo de problemas.
Supongamos un problema de valores propios de orden 2s en un dominio Ω con solu-
ción exacta perteneciente a un espacio de Hilbert1 de orden r. Dicho dominio se discretiza
con una malla cuasiuniforme: siendo he el diámetro de la esfera circunscrita en el elemento
e, tenemos que h = máx(he); como cada elemento ha de ser de similar tamaño se habrá
de cumplir h/he < τ , donde τ es un número finito real. En concreto, definiendo el factor
de uniformidad de una malla como:
FU =
h
mı́n(he)
(2.20)
En una malla totalmente uniforme, este factor seŕıa igual a 1. Lógicamente, en
geometŕıas complejas esto es dif́ıcil de conseguir. Suponemos que el orden de los elementos
usados para la interpolación es p. La solución exacta puede ser regular y suficientemente
suave (r > p), o no suficientemente suave (r < p); esto se puede deber a varias razones:
contornos no suaves, cambios bruscos en las constantes f́ısicas que caracterizan los medios
del dominio, excitaciones singulares y discontinuidades en las condiciones de contorno
[Babuska, 1977]. En problemas electromagnéticos, acoplos y efectos de proximidad entre
conductores también producen soluciones no suaves. Según [Salazar-Palma et al., 1998],
asumiendo que h < 1 el error en el autovalor i-ésimo vendrá dado por la expresión:
|λi − λ̃i| ≤ C1h
2(k+1−s)λ
k+1
s
i , k = mı́n(p, r) (2.21)
donde λi es el valor exacto de dicho autovalor y λ̃i el computado.
Para problemas de segundo orden se obtiene:
|λi − λ̃i| ≤ C2h
2kλk+1
i (2.22)
Como se observa, la convergencia para autovalores mayores en magnitud es mayor
que para los más pequeños. Es importante también ver que si la solución que se está inten-
tando alcanzar es no suave, a pesar de usar elementos de orden p, el error en el autovalor
viene determinado como si usasen elementos de orden r. En dichos casos, se entiende que
una malla cuasiuniforme el error de interpolación es mayor en aquellos elementos localiza-
dos en las regiones del dominio donde el comportamiento de la solución es más abrupto.
Aśı pues, si se generase una malla no uniforme con elementos más pequeños en dichas
regiones, los errores locales estaŕıan todos en el mismo orden de magnitud. A tales mallas
1consultar [Salazar-Palma et al., 1998] apéndice A para las cuestiones matemáticas sobre espacios vec-
toriales en relación al método de los elementos finitos
2.3 Postproceso de la solución 11
se les denomina consecuentemente mallas equilibradas. Para poder comparar en términos
de convergencia tales mallas con las uniformes, se habrá de expresar el error en términos
del número de grados de libertad Nd, de tal manera que la expresión 2.22 queda como:
|λi − λ̃i| ≤ C3N
− 2
n
mı́n(p,r)
d (2.23)
donde n es la dimensión del problema; en nuestro caso n = 3. Aśı pues, si se usan mallas
equilibradas, obtendremos que:
|λi − λ̃i| ≤ C4N
− 2
n
p
d (2.24)
De esta última expresión se deduce que para elementos de segundo orden, como
los que se usarán en adelante, el error en el autovalor ha de converger al menos con una
pendiente de −4/3. Esto será fácil de lograr para problemas que presenten una solución
suave; el algoritmo autoadaptativo desarrollado intentará que también se alcance dicha
cota para cualquier problema, aún presentando éste una solución no suave.
2.3. Postproceso de la solución
A continuación se llevará a cabo una presentación del desarrollo realizado para obte-
ner el vector solución y otras cantidades relacionadas con el mismo. Se ha de mantener pre-
sente en todo momento que el tipo de elementos usados son los tetraedros curl-conformes
subparamétricos de 2o orden descritos en el apéndice A; además, se recuerda que en to-
do caso que los tensores de permitividad eléctrica y permeabilidad magnética son reales,
simétricos y constantes en el elemento.
2.3.1. Cálculo de ~V (x, y, z)
Del Apéndice A sabemos que la solución interpolada en cada elemento toma la forma:
~V (x, y, z) = ~Υ(ξ, η, ζ) = ([Je]T )−1~V(ξ, η, ζ)





Vx(x, y, z)
Vy(x, y, z)
Vz(x, y, z)





=





Υx(ξ, η, ζ)
Υy(ξ, η, ζ)
Υz(ξ, η, ζ)





= ([Je]T )−1





Vξ(ξ, η, ζ)
Vη(ξ, η, ζ)
Vζ(ξ, η, ζ)





(2.25)
Expresando de forma matricial la función vectorial respecto de las coordenadas ca-
12 Análisis de cavidades resonantes mediante el MEF
nónicas, la ecuación 2.25 se puede escribir como:
~Υ(ξ, η, ζ) = [i]{Π2}

















Υx(ξ, η, ζ)
Υy(ξ, η, ζ)
Υz(ξ, η, ζ)

















=









i1,1 i1,2 i1,3 i1,4 i1,5 i1,6 i1,7 i1,8 i1,9 i1,10
i2,1 i2,2 i2,3 i2,4 i2,5 i2,6 i2,7 i2,8 i2,9 i2,10
i3,1 i3,2 i3,3 i3,4 i3,5 i3,6 i3,7 i3,8 i3,9 i3,10

























































































1
ξ
η
ζ
ξ2
ξη
ξζ
η2
ηζ
ζ2

















































































(2.26)
donde [i] se define como la matriz de interpolación y {Π2} como el vector de coordenadas
de segundo orden. A su vez, teniendo en cuenta que al estar utilizando una transforma-
ción geométrica subparamétrica (de primer orden) ([Je]T )−1 es una matriz de elementos
constantes, la matriz de interpolación puede escribirse como:
[i] = ([Je]T )−1[ι] (2.27)
Si expresamos matricialmente las funciones de interpolación en el elemento canónico
(ver Apéndice A) obtenemos:







~N
(k)
ξ
~N
(k)
η
~N
(k)
ζ







=
=



a1
(k) a2
(k) a3
(k) a4
(k) 0 −F (k) −G(k) D(k) J(k) H(k)
b1
(k) b2
(k) b3
(k) b4
(k) F (k) −D(k) −J+K(k) 0 −E(k) I(k)
c1(k) c2(k) c3(k) c4(k) G(k) −K(k) −H(k) E(k) −I(k) 0
































1
ξ
η
ζ
ξ2
ξη
ξζ
η2
ηζ
ζ2





























(2.28)
De donde podemos definir la matriz de interpolación en el elemento canónico [ι]
como la suma de las matrices de coeficientes de cada nodo por el correspondiente valor de
su(s) grado(s) de libertad:
[ι] =
20
∑
k=1



a1
(k) a2
(k) a3
(k) a4
(k) 0 −F (k) −G(k) D(k) J(k) H(k)
b1
(k) b2
(k) b3
(k) b4
(k) F (k) −D(k) −J+K(k) 0 −E(k) I(k)
c1(k) c2(k) c3(k) c4(k) G(k) −K(k) −H(k) E(k) −I(k) 0



dk (2.29)
De esta forma obtenemos:
~V (x, y, z) = ~Υ(ξ, η, ζ) = ([Je]T )−1
20
∑
k=1
{ ~N (k)}dk = [i]{Π2} (2.30)
2.3 Postproceso de la solución 13
2.3.2. Cálculo de ∇× ~V (x, y, z)
Para el cálculo de ∇× ~V (x, y, z) partimos de la definición de rotacional:
∇× ~V (x, y, z) =
∣
∣
∣
∣
∣
∣
∣
∣
x̂ ŷ ẑ
∂
∂x
∂
∂y
∂
∂z
Vx Vy Vz
∣
∣
∣
∣
∣
∣
∣
∣
=







∂Vz
∂y − ∂Vy
∂z
∂Vx
∂z − ∂Vz
∂x
∂Vy
∂x − ∂Vx
∂y







(2.31)
Para calcular las derivadas parciales de ~V aplicamos la regla de la cadena en la ecua-
ción 2.25 (consultar las referencias [Marsden and Tromba, 1991], [Salas and Hille, 1986]
para las cuestiones referentes al cálculo vectorial en tres dimensiones), obteniendo la igual-
dad:




∂Υx
∂ξ
∂Υx
∂η
∂Υx
∂ζ
∂Υy
∂ξ
∂Υy
∂η
∂Υy
∂ζ
∂Υz
∂ξ
∂Υz
∂η
∂Υz
∂ζ




=




∂Vx
∂x
∂Vx
∂y
∂Vx
∂z
∂Vy
∂x
∂Vy
∂y
∂Vy
∂z
∂Vz
∂x
∂Vz
∂y
∂Vz
∂z








∂x
∂ξ
∂x
∂η
∂x
∂ζ
∂y
∂ξ
∂y
∂η
∂y
∂ζ
∂z
∂ξ
∂z
∂η
∂z
∂ζ




(2.32)
donde el segundo multiplicando del término de la derecha resulta ser la matriz jacobiana
de la transformación geométrica.
De la ecuación 2.32 se obtiene inmediatamente:
∂ξ,η,ζ
~Υ = ∂x,y,z
~V [Je]
∂x,y,z
~V = ∂ξ,η,ζ
~Υ[Je]−1
∂x,y,z
~V = ([Je]T )−1∂ξ,η,ζ
~V[Je]−1
(2.33)
Expresando ahora ~Υ como en la ecuación 2.30 y recordando que la única dependencia
con (ξ, η, ζ) la presenta el vector de coordenadas de segundo orden {Π2} podemos escribir:
∂ξ,η,ζ
~Υ = [i]∂ξ,η,ζ{Π2}




∂Υx
∂ξ
∂Υx
∂η
∂Υx
∂ζ
∂Υy
∂ξ
∂Υy
∂η
∂Υy
∂ζ
∂Υz
∂ξ
∂Υz
∂η
∂Υz
∂ζ




= [i]





















0 0 0
1 0 0
0 1 0
0 0 1
2ξ 0 0
η ξ 0
ζ 0 ξ
0 2η 0
0 ζ η
0 0 2ζ





















(2.34)
Si introducimos la ecuación anterior en la ecuación 2.33, y teniendo nuevamente en
cuenta que al estar utilizando una transformación geométrica subparamétrica (de primer
orden) [Je]−1 es una matriz de elementos constantes, que definiremos como:
[j] = [Je]−1 (2.35)
podemos expresar las derivadas parciales de ~V como:
14 Análisis de cavidades resonantes mediante el MEF
∂x,y,z
~V = [i]∂ξ,η,ζ{Π2}[Je]−1




∂Vx
∂x
∂Vx
∂y
∂Vx
∂z
∂Vy
∂x
∂Vy
∂y
∂Vy
∂z
∂Vz
∂x
∂Vz
∂y
∂Vz
∂z




= [i]





















0 0 0
1 0 0
0 1 0
0 0 1
2ξ 0 0
η ξ 0
ζ 0 ξ
0 2η 0
0 ζ η
0 0 2ζ
























j1,1 j1,2 j1,3
j2,1 j2,2 j2,3
j3,1 j3,2 j3,3



(2.36)
Desarrollando el segundo término de la ecuación anterior llegamos al producto ma-
tricial:




i1,1 i1,2 i1,3 i1,4 i1,5 i1,6 i1,7 i1,8 i1,9 i1,10
i2,1 i2,2 i2,3 i2,4 i2,5 i2,6 i2,7 i2,8 i2,9 i2,10
i3,1 i3,2 i3,3 i3,4 i3,5 i3,6 i3,7 i3,8 i3,9 i3,10

























0 0 0
j1,1 j1,2 j1,3
j2,1 j2,2 j2,3
j3,1 j3,2 j3,3
j1,12ξ j1,22ξ j1,32ξ
j1,1η + j2,1ξ j1,2η + j2,2ξ j1,3η + j2,3ξ
j1,1ζ + j3,1ξ j1,2ζ + j3,2ξ j1,3ζ + j3,3ξ
j2,12η j2,22η j2,32η
j2,1ζ + j3,1η j2,2ζ + j3,2η j2,3ζ + j3,3η
j3,12ζ j3,22ζ j3,32ζ





















(2.37)
Realizando dicho producto, podemos expresar finalmente las derivadas parciales de
~V como:
∂Vtn
∂τm
= (d(n,m)){Π1}
∂Vtn
∂τm
=
(
d
(n,m)
0 d
(n,m)
1 d
(n,m)
2 d
(n,m)
3
)











1
ξ
η
ζ











n = 1 . . . 3
m = 1 . . . 3
t1 = x t2 = y t3 = z
τ1 = ξ τ2 = η τ3 = ζ
(2.38)
donde:
d
(n,m)
0 = in,2j1,m + in,3j2,m + in,4j3,m
d
(n,m)
1 = 2in,5j1,m + in,6j2,m + in,7j3,m
d
(n,m)
2 = in,6j1,m + 2in,8j2,m + in,9j3,m
d
(n,m)
3 = in,7j1,m + in,9j2,m + 2in,10j3,m
(2.39)
2.3 Postproceso de la solución 15
Introduciendo las derivadas parciales obtenidas mediante la ecuación 2.38 en la de-
finición de rotacional (ecuación 2.31) obtenemos finalmente:
~Ψ(ξ, η, ζ)
.
= ∇× ~V (x, y, z) = [r]{Π1}





Ψx(ξ, η, ζ)
Ψy(ξ, η, ζ)
Ψz(ξ, η, ζ)





.
=





(∇× ~V )x(x, y, z)
(∇× ~V )y(x, y, z)
(∇× ~V )z(x, y, z)





=



r1,1 r1,2 r1,3 r14
r2,1 r2,2 r2,3 r24
r3,1 r3,2 r3,3 r34














1
ξ
η
ζ











(2.40)
donde {Π1} se define como vector de coordenadas de primer orden, mientras que [r] es la
matriz de interpolación del rotacional definida por:
[r] =



(d(3,2)) − (d(2,3))
(d(1,3)) − (d(3,1))
(d(2,1)) − (d(1,2))



(2.41)
2.3.3. Cálculo de ∇× ν−1∇× ~V (x, y, z)
Análogamente al procedimiento seguido en el apartado anterior, se calcularán pri-
mero las derivadas parciales de la función ν−1∇× ~V (x, y, z) necesarias para el cálculo del
rotacional de dicha función
Teniendo en cuenta que∇×~V (x, y, z) puede ser expresado como una función vectorial
dependiente de las coordenadas canónicas ξ, η, ζ, tal y como muestra la ecuación 2.40,
podemos escribir:
ν−1∇× ~V (x, y, z) = ν−1~Ψ(ξ, η, ζ) (2.42)
Aplicando como en la seccion 2.3.2 la regla de la cadena sobre esta ecuación obte-
nemos:





∂(ν−1~Ψ)x
∂ξ
∂(ν−1~Ψ)x
∂η
∂(ν−1~Ψ)x
∂ζ
∂(ν−1~Ψ)y
∂ξ
∂(ν−1~Ψ)y
∂η
∂(ν−1~Ψ)y
∂ζ
∂(ν−1~Ψ)z
∂ξ
∂(ν−1~Ψ)z
∂η
∂(ν−1~Ψ)z
∂ζ





=
=





∂(ν−1∇×~V )x
∂x
∂(ν−1∇×~V )x
∂y
∂(ν−1∇×~V )x
∂z
∂(ν−1∇×~V )y
∂x
∂(ν−1∇×~V )y
∂y
∂(ν−1∇×~V )y
∂z
∂(ν−1∇×~V )z
∂x
∂(ν−1∇×~V )z
∂y
∂(ν−1∇×~V )z
∂z









∂x
∂ξ
∂x
∂η
∂x
∂ζ
∂y
∂ξ
∂y
∂η
∂y
∂ζ
∂z
∂ξ
∂z
∂η
∂z
∂ζ




(2.43)
lo que inmediatamente nos conduce a:
∂x,y,z(ν
−1∇× ~V ) = ∂ξ,η,ζ(ν
−1~Ψ)[Je]−1 (2.44)
Expresando ahora de forma matricial la función ~Ψ según la ecuación 2.40 y recordan-
do que la única dependencia con (ξ, η, ζ) la presenta el vector de coordenadas de primer
orden {Π1} tenemos:
∂ξ,η,ζ(ν
−1~Ψ) = [ν]−1[r]∂ξ,η,ζ{Π1} (2.45)
16 Análisis de cavidades resonantes mediante el MEF
que insertado en la ecuación 2.44 nos da:
∂x,y,z(ν
−1∇× ~V ) = [ν]−1[r]∂ξ,η,ζ{Π1}[Je]−1




∂(ν−1∇×~V )x
∂x
∂(ν−1∇×~V )x
∂y
∂(ν−1∇×~V )x
∂z
∂(ν−1∇×~V )y
∂x
∂(ν−1∇×~V )y
∂y
∂(ν−1∇×~V )y
∂z
∂(ν−1∇×~V )z
∂x
∂(ν−1∇×~V )z
∂y
∂(ν−1∇×~V )z
∂z





= [ν]−1[r]






0 0 0
1 0 0
0 1 0
0 0 1






[Je]−1
(2.46)
Finalmente, insertando las derivadas parciales obtenidas mediante la ecuación ante-
rior en la definición de rotacional (ecuación 2.31) llegamos a:
∇× ν−1∇× ~V =









∂(ν−1∇×~V )z
∂y − ∂(ν−1∇×~V )y
∂z
∂(ν−1∇×~V )x
∂z − ∂(ν−1∇×~V )z
∂x
∂(ν−1∇×~V )y
∂x − ∂(ν−1∇×~V )x
∂y









=





p1
p2
p3





= {p} (2.47)
donde {p} se define, por similitud con los casos anteriores, como la columna de interpolación
del rotacional doble. Por motivos que se verán más adelante (consultar sección 3.2.3)
definimos por conveniencia la matriz de interpolación del rotacional doble como:
[p] = [{p} [0]3×9] =



p1 0 0 0 0 0 0 0 0 0
p2 0 0 0 0 0 0 0 0 0
p3 0 0 0 0 0 0 0 0 0



(2.48)
De tal manera, ∇× ν−1∇× ~V (x, y, z) puede expresarse como:
∇× ν−1∇× ~V = [p]{Π2} (2.49)
2.3.4. Visualización gráfica
Para poder visualizar cualquier cantidad de las calculadas anteriormente, habremos
de especificar una serie de puntos xi, yi, zi y calcular el valor de dicha cantidad en esos
puntos. Para ello habremos de encuadrar cada punto en el elemento que pertenezca, realizar
la transformación a las coordenadas canónicas y, finalmente, llevar a cabo la interpolación
a nivel local. La transformación geométrica en los elementos usados viene dada por:





x
y
z





=



x1 x2 x3 x4
y1 y2 y3 y4
z1 z2 z3 z4














1 − ξ − η − ζ
ξ
η
ζ











(2.50)
La ecuación 2.50 puede ser fácilmente expresada como:





x
y
z





=



x2 − x1 x3 − x1 x4 − x1
y2 − y1 y3 − y1 y4 − y1
z2 − z1 z3 − z1 z4 − z1








ξ
η
ζ





+





x1
y1
z1





{v} = [Je]





ξ
η
ζ





+ {v1}
(2.51)
2.3 Postproceso de la solución 17
Aśı pues, las coordenadas canónicas en el elemento de referencia del punto x, y, z
vienen expresadas como:





ξ
η
ζ





= [Je]−1 ({v} − {v1}) (2.52)
Lógicamente, a priori no se sabe a qué elemento e pertenece el punto en cuestión.
Para averiguarlo, una vez tenemos las coordenadas en el elemento canónico, están han de
cumplir obligatoriamente:
0 ≤ ξ ≤ 1
0 ≤ η ≤ 1 − ξ
0 ≤ ζ ≤ 1 − ξ − η
(2.53)
Si no es aśı, el punto x, y, z pertenecerá a otro elemento.
Una vez tenemos las coordenadas ξ, η, ζ basta aplicar las mismas a los vectores de
coordenadas {Π1} ó {Π2} y multiplicar por la matriz correspondiente para calcular el
valor del vector deseado en el punto x, y, z. Con una cantidad suficiente de puntos se
puede realizar una representación gráfica de cualquiera de los vectores calculados en el
postproceso. Bajo esta idea se ha desarrollado un módulo de dibujo de que se hará uso en
el capitulo 5 para visualizar los campos ~E y ~H.
18 Análisis de cavidades resonantes mediante el MEF
Caṕıtulo 3
Indicadores de Error para
problemas electromagnéticos de
autovalores en 3D
Para saber en qué regiones la solución calculada del campo modela de una manera
menos exacta la solución verdadera, se habrá de estimar o dar una indicación del error
cometido tanto de manera global como a nivel de elemento. Como ya se ha comentado
en la introducción, en este caṕıtulo se presentan dos tipos de indicadores de error que
actualmente son ampliamente usados.
3.1. Errores a posteriori en el MEF
3.1.1. Error Real VS Error Estimado en problemas electromagnéticos
El error cometido al resolver mediante el método de los elementos finitos un problema
de tipo vectorial (sea campo ~E o ~H) toma la forma:
~ǫ = ~V − ~V (3.1)
donde ~V representa la solución real del problema y ~V la obtenida mediante el MEF.
Existen diferentes tipos de error, y el error total cometido es la suma de las contri-
buciones de los mismos. No es el propósito de este trabajo dar aqúı una lista completa
y detallada de los mismos, que puede fácilmente encontrarse en la literatura (consultar
[Salazar-Palma et al., 1998], [D́ıaz-Morcillo, 2000], [Nédélec, 1992]), aunque si cabe des-
tacar que algunos, como los de discretización e interpolación, son inherentes al propio
método de elementos finitos, mientras que otros, como los de truncamiento, dependen de
las capacidades computacionales de la máquina donde se lleve a cabo la implementación
de dicho método.
El propio hecho de hacer uso de un método numérico para resolver un problema
implica que su solución exacta no se conoce, por lo que su error real ~ǫ también es imposible
conocerlo. Sin embargo, se puede conseguir un error estimado a posteriori ~e. Una vez más,
no es propósito mencionar aqúı todos los tipos, pues en la bibliograf́ıa se pueden encontrar
multitud de clases de estimadores (residuales, complementarios, etc . . . ), destacándose
aqúı la lista proporcionada en [Salazar-Palma et al., 1998].
19
20 Indicadores de Error para problemas electromagnéticos de autovalores en 3D
Bien es cierto que, como ya sucedió con el propio MEF en sus inicios, la mayor parte
de la teoŕıa de estimación de errores se ha desarrollado en el campo de la ingenieŕıa civil,
mecánica o de fluidos. Sin embargo, en los últimos años se han logrado un gran número de
avances en lo referente al campo del electromagnetismo, estando esta parcela en continua
evolución.
3.1.2. Medida del error
Para realizar una medida fiable, una estimación o indicación del error cometido,
recurrimos a calcular la norma sobre el vector ~e. De esta manera, particularizando ya
sobre los elementos que forman nuestro mallado, podemos obtener que la estimación o
indicación del error total cometido viene dada por:
‖~eT ‖ =
Ne
∑
e=1
‖~ee‖ (3.2)
donde Ne es el número de elementos que componen la malla.
La norma elegida puede ser simplemente la norma L2 o la norma enerǵıa:
‖~ee‖L2 =
√
∫∫∫
Ωe
~ee∗ · ~eedΩ (3.3)
‖~ee‖E =
√
∫∫∫
Ωe
~ee∗ · L · ~eedΩ (3.4)
donde L es el operador diferencial propio de la formulación inicial del problema. En los
problemas que nos atañen, L = ∇× ν−1∇×−k2
0ϑ.
3.1.3. Indicador VS Estimador del Error
Hasta ahora se ha hablado de manera aparentemente indiferente sobre estimación e
indicación del error. Muchos autores, de hecho, no realizan distinción entre ambas definicio-
nes. Sin embargo aqúı se adoptará la postura mencionada en [Salazar-Palma et al., 1998],
donde por estimador de error se entiende una cantidad ‖~eT ‖ que conforme aumentan las
iteraciones en el proceso de mallado autoadaptativo se aproxima al error real cometido
‖~ǫ‖, de tal manera que:
‖~eT ‖/‖~ǫ‖ → 1 (3.5)
La cantidad presentada en 3.5 se denomina ı́ndice de efectividad y se aplica en pro-
blemas de tipo determinista, como el análisis cuasiestático de ĺıneas de transmisión, en
[Salazar-Palma et al., 1998]. En la misma referencia, para problemas de autovalores (en
2D), se introduce el concepto de indicador de error, que no constituye una verdadera esti-
mación del mismo, pero que como su propio nombre refleja, indica las zonas con una mayor
error respecto al resto (elementos con ‖~ee‖ mayores), y si el error total ‖~eT ‖ disminuye
con el paso de las iteraciones. Aśı pues, en los problemas que se discuten en lo siguiente, al
tratarse de problemas de autovalores en 3D, se hablará siempre en términos del indicador
de error, aunque quizá a veces se haga uso de la palabra estimador, bien por relajación
dada la ambigüedad que envuelve ambos términos, bien por evitar repetición continuada
de una misma palabra.
3.2 Indicador de Error Residual 21
3.2. Indicador de Error Residual
Este indicador de error ha sido tradicionalmente usado desde los inicios de la estima-
ciónde error a posteriori en el MEF. En esta sección primeramente se aplica la teoŕıa sobre
los problemas que se están tratando con el tipo de elementos usados, para posteriormente,
en las secciones 3.2.3 y 3.2.4, afrontar el cálculo de los residuos implicados.
3.2.1. Error Residual
Una medida de la norma enerǵıa del error cometido en cada elemento se pue-
de calcular a través de los residuos en cada tetraedro mediante la expresión dada en
[D́ıaz-Morcillo, 2000] que constituye una extensión de la que se puede encontrar en la
referencia [Salazar-Palma et al., 1998]:
‖~ee‖2 = fv
h2
e
λf−1
e,min
p
∫∫∫
Ωe
~re∗
v · ~re
vdΩ + fs
he
λf−1
e,min
p
ne
c
∑
k=1
∫∫
Γk /∈ΓD
~re∗
sk
· ~re
sk
dΓk (3.6)
donde
λf−1
e,min
es el menor autovalor del tensor fe
−1, esto es:
λf−1
e,min
= mı́n(λ) (3.7)
sujeto a
[fe]
−1 − λ[I] = [0] (3.8)
p es el orden de las funciones de interpolación utilizadas, en nuestro caso:
p = 2 (3.9)
he es una medida del tamaño del elemento, definida por el diámetro de la esfera
circunscrita en el tetraedro, o lo que es lo mismo, la longitud de su arista más larga
o mayor distancia entre sus vértices.
he = máx(lij), 1 ≤ i ≤ j ≤ 4 (3.10)
ne
c es el número de caras del elemento, que en nuestro caso al tratarse de un tetraedro
resulta:
ne
c = 4 (3.11)
El primer término de la expresión 3.6 constituye una medida del residuo interno o
residuo volumétrico, que viene dado por:
~re
v = ∇× fe
−1∇× ~̃W e − ω2ge
~̃W e (3.12)
El segundo término de la expresión 3.6 representa una medida del residuo singular
que se evalúa en las caras del elemento que no presentan condiciones de contorno de
Dirichlet (donde es igual a 0). En este caso hay que considerar dos situaciones:
22 Indicadores de Error para problemas electromagnéticos de autovalores en 3D
Caras interiores de la malla compartidas por dos elementos. El residuo es:
~re
sk
|int = n̂ × (fe1
−1∇× ~̃W e1 − fe2
−1∇× ~̃W e2) (3.13)
Caras pertenecientes a la frontera bajo condiciones de contorno de Neumann. En
este caso el residuo será:
~re
sk
|Neumann = n̂ × fe
−1∇× ~̃W e (3.14)
Finalmente, fv y fs son factores de ponderación de la medida de estos residuos que
generalmente suelen cumplir:
fv + fs = 1, fv ≥ 0, fs ≥ 0 (3.15)
3.2.2. Error Residual para la formulación normalizada
En este caso se puede demostrar fácilmente que (consultar sección 2.1.2 para más
detalles), elidiendo el ı́ndice e, el indicador de error en el elemento queda como:
‖~e‖2 = fv
h2
λν−1
min
p
∫∫∫
Ω
~r∗v · ~rvdΩ + fs
h
λν−1
min
p
nc
∑
k=1
∫∫
Γk /∈ΓD
~r∗sk
· ~rsk
dΓk (3.16)
El residuo interior o volumétrico viene ahora dado por:
~rv = ∇× ν−1∇× ~V − k2
0ϑ
~V (3.17)
Mientras que la expresión para el residuo singular en las caras es:
Caras interiores de la malla compartidas por dos elementos.
~rsk
|int = n̂ × (ν1
−1∇× ~V 1 − ν2
−1∇× ~V 2) (3.18)
Caras pertenecientes a la frontera bajo condiciones de contorno de Neumann.
~rsk
|Neumann = n̂ × ν−1∇× ~V (3.19)
Caras pertenecientes a la frontera bajo condiciones de contorno de Dirichlet.
~rsk
|Dirichlet = 0 (3.20)
El resto de los términos (p, h, nc, fv, fs) permanecen inalterados.
3.2.3. Residuo Interior o Volumétrico
De las ecuaciones 2.30 y 2.49 sabemos que el residuo interior puede ser expresado
como:
~rv = ∇× ν−1∇× ~V − k2
0ϑ~V = ([p] − k2
0 [ϑ][i]){Π2} (3.21)
Si definimos la matriz de residuo volumétrico como:
[rv] = ([p] − k2
0[ϑ][i]) (3.22)
3.2 Indicador de Error Residual 23
tenemos que la ecuación 3.21 se puede expresar de la siguiente manera:
~rv = [rv]{Π2} =
=




rv
x,c rv
x,ξ rv
x,η rv
x,ζ rv
x,ξ2 rv
x,ξη rv
x,ξζ rv
x,η2 rv
x,ηζ rv
x,ζ2
rv
y,c rv
y,ξ rv
y,η rv
y,ζ rv
y,ξ2 rv
y,ξη rv
y,ξζ rv
y,η2 rv
y,ηζ rv
y,ζ2
rv
z,c rv
z,ξ rv
z,η rv
z,ζ rv
z,ξ2 rv
z,ξη rv
z,ξζ rv
z,η2 rv
z,ηζ rv
z,ζ2













































1
ξ
η
ζ
ξ2
ξη
ξζ
η2
ηζ
ζ2









































(3.23)
Para realizar la medida de este residuo hemos de integrar en el elemento canónico,
pues dicho residuo está expresado en coordenadas canónicas. Aśı pues, si queremos tener
una medida en nuestro dominio habremos de aplicar la fórmula del cambio de variables
para la integración:
∫∫∫
Ωe
|~rv(x, y, z)|2dxdydz =
∫∫∫
Ω̂
|~rv(ξ, η, ζ)|2|Je|dξdηdζ (3.24)
Si utilizamos, como es nuestro caso, elementos de 2o orden subparamétricos, el de-
terminante de la matriz jacobiana será constante, por lo que finalmente podremos expresar
la medida del residuo como:
∫∫∫
Ωe
|~rv|2dΩ = |Je|
(
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0
|rv,x|2dξdηdζ +
+
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0
|rv,y|2dξdηdζ +
+
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0
|rv,z|2dξdηdζ
)
=
= |Je|
3
∑
n=1
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0
|rv,tn |2dξdηdζ (3.25)
t1 = x t2 = y t3 = z
De la ecuación 3.23 se obtiene que:
|rv,tn
|2 = (rv
tk,c+rv
tk,ξξ+rv
tk,ηη+rv
tk,ζζ+rv
tk,ξ2ξ2+rv
tk,ξηξη+rv
tk,ξζξζ+rv
tk,η2η2+rv
tk,ηζηζ+rv
tk,ζ2ζ2)2
(3.26)
Refiriendo al lector al apéndice B para realizar la cuadratura numérica de la integral
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0 |rv,tn |2dξdηdζ, queda finalmente:
24 Indicadores de Error para problemas electromagnéticos de autovalores en 3D
∫ 1
ξ=0
∫ 1−ξ
η=0
∫ 1−ξ−η
ζ=0
|rv,tn |2dξdηdζ =
=
(rv
tn,c)
2
6
+
+
2(rv
tn,cr
v
tn,ξ + rv
tn,cr
v
tn,η + rv
tn,cr
v
tn,ζ)
24
+
+
2(rv
tn,cr
v
tn,ξη + rv
tn,ξr
v
tn,η + rv
tn,cr
v
tn,ξζ + rv
tn,ξr
v
tn,ζ + rv
tn,cr
v
tn,ηζ + rv
tn,ηr
v
tn,ζ)
120
+
+
2(rv
tn,ξr
v
tn,ηζ + rv
tn,ηr
v
tn,ξζ + rv
tn,ζr
v
tn,ξη)
720
+
+
2rv
tn,c(r
v
tn,ξ2 + rv
tn,η2 + rv
tn,ζ2) + (rv
tn,ξ)
2 + (rv
tn,η)
2 + (rv
tn,ζ)
2
60
+
+
2(rv
tn,ξr
v
tn,ξη + rv
tn,ξ2r
v
tn,η + rv
tn,ξr
v
tn,ξζ + rv
tn,ξ2r
v
tn,ζ)
360
+
+
2(rv
tn,ηr
v
tn,ξη + rv
tn,η2r
v
tn,ξ + rv
tn,ηr
v
tn,ηζ + rv
tn,η2r
v
tn,ζ)
360
+
+
2(rv
tn,ζr
v
tn,ξζ + rv
tn,ζ2r
v
tn,ξ + rv
tn,ζr
v
tn,ηζ + rv
tn,ζ2r
v
tn,η)
360
+
+
2(rv
tn,ξηr
v
tn,ξζ + rv
tn,ξ2r
v
tn,ηζ + rv
tn,ξηr
v
tn,ηζ + rv
tn,η2r
v
tn,ξζ + rv
tn,ξζr
v
tn,ηζ + rv
tn,ζ2r
v
tn,ξη)
2520
+
+
2rv
tn,ξ2r
v
tn,η2 + (rv
tn,ξη)
2 + 2rv
tn,ξ2r
v
tn,ζ2 + (rv
tn,ξζ)
2 + 2rv
tn,η2r
v
tn,ζ2 + (rv
tn,ηζ)
2
1260
+
+
2(rv
tn,ξ2r
v
tn,ξ + rv
tn,η2r
v
tn,η + rv
tn,ζ2r
v
tn,ζ)
120
+
+
2(rv
tn,ξ2r
v
tn,ξη + rv
tn,ξ2r
v
tn,ξζ + rv
tn,η2r
v
tn,ξη + rv
tn,η2r
v
tn,ηζ + rv
tn,ζ2r
v
tn,ξζ + rv
tn,ζ2r
v
tn,ηζ)
840
+
+
(rv
tn,ξ2)
2 + (rv
tn,η2)
2 + (rv
tn,ζ2)
2
210
(3.27)
3.2.4. Residuo Singular en las Caras
La evaluación de la medida del residuo singular en cada cara
∫∫
Γk /∈ΓD
|~rsk
|2dΓk (3.28)
presenta un nivel mayor de complejidad que la del residuo interior ya que si se trata de
una cara interior del mallado se ha de realizar la interpolación a nivel local sobre dos ele-
mentos distintos, tal y como indica la ecuación 3.18 cuyas transformaciones geométricas
son distintas y, por tanto, cuyas coordenadas canónicas no pueden agruparse dentro del
mismo polinomio. Cabe la posibilidad de realizar una doble transformación geométrica
de coordenadas del elemento de referencia al real (de 3 dimensiones a 3 dimensiones), y,
posteriormente, una transformación de la cara implicada a un triángulo canónico (de 3
dimensiones a 2 dimensiones), pudiendo utilizar una cuadratura similar a la del apartado
anterior, esta vez sobre un triángulo. Sin embargo, tal operación se antoja dif́ıcil si sólo
3.2 Indicador de Error Residual 25
Tabla 3.1: Vértices que definen las caras de un tetraedro (numeración local)
Ci V
(Ci)
1 V
(Ci)
2 V
(Ci)
3
1 1 3 2
2 1 4 3
3 1 2 4
4 2 3 4
Tabla 3.2: Aristas que definen las caras de un tetraedro (numeración local)
Ci A
(Ci)
1 A
(Ci)
2 A
(Ci)
3
1 3 1 2
2 4 3 6
3 1 4 5
4 2 5 6
se quiere hacer uso de la interpolación local. Por ello se ha optado por una cuadratura de
Gauss-Legendre