TFG Pedro Palomo

Trabajo Fin de Grado
Grado en Ingeniería Civil
Programación en Matlab del Método de los

Elementos de Contorno aplicado a problemas
acústicos
Autor: Pedro Palomo Cantador
Tutor: Luis Rodríguez de Tembleque Solano
Equation Chapter 1 Section 1
Departamento de Mecánica de los Medios

Continuos y Teoría de Estructuras
Universidad de Sevilla
Escuela Técnica Superior de Ingeniería
Sevilla, Junio de 2015

Trabajo Fin de Grado
Grado en Ingeniería Civil
Programación en Matlab del Método de los

Elementos de Contorno aplicado a problemas
acústicos
Departamento de Mecánica de los Medios Continuos y

Teoría de Estructuras
Universidad de Sevilla
Escuela Técnica Superior de Ingeniería
Sevilla, Junio de 2015
ii
Trabajo Fin de Grado: Programación en Matlab del Método de los Elementos de Contorno aplicado a
problemas acústicos.
El tribunal nombrado para juzgar el Trabajo arriba indicado, compuesto por los siguientes miembros:
Presidente:
Vocales:
Secretario:
Acuerdan otorgarle la calificación de:
Sevilla, 2015
El Secretario del Tribunal

Agradecimientos
Quisiera dedicar esta sección a agradecer a mi tutor Luis Rodrı́guez de Tembleque

Solano, por darme la oportunidad de realizar este trabajo y por guiarme, dándome además
un gran apoyo durante toda la realización del mismo. Te considero un referente tanto a
nivel académico como humano.
Me gustarı́a dedicarle este trabajo también a mis padres, por su apoyo incondicional y
su entrega durante toda la carrera, para facilitármelo todo y que solo tuviera que centrarme
en los estudios. Nunca me faltó vuestro consuelo en los momentos de debilidad y tristeza,
gracias.
Al resto de mis familiares: mi hermano, mis tı́os, mis abuelos . . . por colaborar en la
medida de lo posible a que nada me faltara y me fuera todo más fácil, os lo agradezco
también profundamente.
Agradecer a mis compañeros de carrera su ayuda e interés por mi persona, no solo en

lo relativo a la universidad, sino también en la vida que hay fuera de ella, me ayudasteis
no solo a mejorar como alumno y compañero, sino también como persona. También a los
que comenzaron y no pudieron terminar, todavı́a hoy los recuerdo y les deseo lo mejor en
sus nuevos desafı́os y oportunidades.
Tampoco quisiera olvidarme de darle las gracias a los profesores que he tenido durante
estos años, por ayudarme a crecer en el plano intelectual, pero sobre todo en el de la cultura
del esfuerzo, el sacrificio y la superación diaria. Estos valores son imprescindibles en la
vida y sin vosotros no los habrá alcanzado del mismo modo.
Por último, quisiera condensar toda mi etapa en la universidad en una frase célebre que
me marcó y me ayudó a querer dar siempre el máximo: El éxito tiene una simple fórmula:
da lo mejor de ti y puede que a la gente le guste-Sam Ewing.
2
3
.
.
.
.
.
.
.
.
.
.
.
.
A mis padres,
Antonio y Cándida.
Índice general
Agradecimientos 2
Índice general 4
Índice de figuras 7
Índice de cuadros 11
1. Introducción 12
1.1. Estudio y motivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.2. Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3. Estructuración del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2. Conceptos fundamentales de acústica 14

2.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2. Ecuaciones básicas del problema acústico . . . . . . . . . . . . . . . . . 14
2.2.1. La ecuación de onda linealizada . . . . . . . . . . . . . . . . . . 14
2.2.2. Ondas armónicas planas: la ecuación de Helmholtz . . . . . . . . 15
2.2.3. Condiciones de contorno . . . . . . . . . . . . . . . . . . . . . . 15
2.3. Magnitudes para medir el sonido . . . . . . . . . . . . . . . . . . . . . . 16
2.3.1. Intensidad acústica . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.2. Impedancia acústica especı́fica . . . . . . . . . . . . . . . . . . . 16
2.3.3. Nivel de presión sonora . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.4. Coeficiente de pérdida por inserción . . . . . . . . . . . . . . . . 17
3. Representación integral para problemas escalares 19

3.1. Relación de reciprocidad . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4
Índice general 5
3.2. Solución fundamental . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.3. Igualdad integral en el contorno . . . . . . . . . . . . . . . . . . . . . . 22
3.4. Condición de radiación en el infinito . . . . . . . . . . . . . . . . . . . . 25
3.5. Igualdad para el problema exterior . . . . . . . . . . . . . . . . . . . . . 26
3.6. Problemas de scatter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.7. No unicidad de la solución . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.8. Uso del Método de los Chief Points . . . . . . . . . . . . . . . . . . . . . 31
3.9. Solución fundamental para el semiespacio . . . . . . . . . . . . . . . . . 31
4. Formulación e implementación del MEC en problemas acústicos 2D 35

4.1. Introducción y descripción general de la metodologı́a . . . . . . . . . . . 35
4.2. Ecuación integral de Helmholtz para problemas bidimensionales . . . . . 35
4.3. Discretización del contorno . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.3.1. Elementos constantes . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3.2. Elementos lineales . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.3. Elementos cuadráticos . . . . . . . . . . . . . . . . . . . . . . . 39
4.3.4. Desarrollo del plantemiento numérico de la identidad integral . . 39
4.4. Utilización del Método de los Chief Points y solución por mı́nimos cuadrados 41
4.5. Aplicación al estudio de pantallas de atenuación acústica . . . . . . . . . 41
4.6. Funciones implementadas . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5. Formulación e implementación del MEC en problemas acústicos 3D 53

5.1. Introducción y descripción general de la metdologı́a . . . . . . . . . . . . 53
5.2. Ecuación integral de Helmholtz para problemas tridimensionales . . . . . 54
5.3. Discretización del contorno . . . . . . . . . . . . . . . . . . . . . . . . . 55
5.3.1. Elementos constantes . . . . . . . . . . . . . . . . . . . . . . . . 56
5.3.2. Elementos lineales . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.3.3. Desarrollo del plantemiento numérico de la identidad integral . . 57
5.4. Utilización del Método de los Chief Points y solución por mı́nimos cuadrados 60
5.5. Funciones implementadas . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6. Ejemplos numéricos 67
6.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.2. Ejemplos de validación . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.2.1. Problema 1 interior 2D . . . . . . . . . . . . . . . . . . . . . . . 67
Índice general 6
6.2.2. Problema 2 interior 2D . . . . . . . . . . . . . . . . . . . . . . . 72

6.2.3. Problema exterior 2D . . . . . . . . . . . . . . . . . . . . . . . . 75
6.2.4. Problema de scatter 2D . . . . . . . . . . . . . . . . . . . . . . . 79
6.2.5. Problema interior 3D . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2.6. Problema exterior 3D . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3. Aplicación al cálculo de pantallas de atenuación acústica . . . . . . . . . 94
6.3.1. Pantalla con contorno reflejante . . . . . . . . . . . . . . . . . . 94
6.3.2. Pantalla con contorno absorbente . . . . . . . . . . . . . . . . . 96
6.4. Análisis de resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
7. Resumen y conclusiones 100

7.1. Resumen de los contenidos desarrollados . . . . . . . . . . . . . . . . . 100
7.2. Conclusiones del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.3. Trabajo futuro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Anexos 104
A. Interacción LaTeX-Matlab para la incorporación de códigos a documentos 105
B. Código acústico interior 2D 108
C. Código acústico exterior 2D 128
D. Código para problemas de scatter 2D 134
E. Código para pantallas acústicas 2D 144
F. Código acústico interior 3D 166
G. Código acústico exterior 3D 209
H. Aplicación del MEC a problemas de potencial 211
I. Código para problemas de potencial 2D 215
Bibliografı́a 240
Índice de figuras
3.1. Ecuación Integral y Solución Fundamental. . . . . . . . . . . . . . . . . . . . . 21

3.2. Igualdad integral en el contorno. . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.3. Término libre. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.4. (a)Problema interior.(b)Problema Exterior. . . . . . . . . . . . . . . . . . . . . 26
3.5. Formulación integral del problema exterior. . . . . . . . . . . . . . . . . . . . . 27
3.6. Problema de scatter con una onda incidente Pi . . . . . . . . . . . . . . . . . . . 28
3.7. Semiespacio 2D con la pantalla. . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.8. Representación de fuente e imagen para la solución fundamental. . . . . . . . . . . 33
4.1. Tipologı́as de elementos empleadas para discretizar. . . . . . . . . . . . . . . . . 37

4.2. Función de forma constante. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.3. Funciones de forma lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4. Función de forma cuadráticas. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.5. Diagrama de flujo del problema interior para elementos constantes. . . . . . . . . . . 50
4.6. Diagrama de flujo del problema exterior para elementos constantes. . . . . . . . . . 50
4.7. Diagrama de flujo del problema de pantallas para elementos cuadráticos. . . . . . . . 51
4.8. Desglose de la función GHMATPQ. . . . . . . . . . . . . . . . . . . . . . . . 51
4.9. Desglose de la función INTERPQ. . . . . . . . . . . . . . . . . . . . . . . . . 52
5.1. Angulo interno en problemas 3D. . . . . . . . . . . . . . . . . . . . . . . . . . 55

5.2. Función de forma constante del problema 3D. . . . . . . . . . . . . . . . . . . . 56
5.3. Ejemplo de función de forma lineal 3D(φ3 ). . . . . . . . . . . . . . . . . . . . . 57
5.4. Elemento j que actúa sobre nodo i. . . . . . . . . . . . . . . . . . . . . . . . . 62
5.5. Elemento j que actúa sobre nodo i. . . . . . . . . . . . . . . . . . . . . . . . . 63
5.6. Diagrama de flujo para problema interior 3D con elementos lineales. . . . . . . . . . 65
5.7. Diagrama de flujo para problema exterior 3D con elementos lineales. . . . . . . . . . 66
7
6.1. Geometrı́a de la cavidad 2D(condiciones de contorno del problema 1). . . . . . . . . 68

6.2. Solución en los nodos del contorno con elementos constantes. . . . . . . . . . . . . 69
6.3. Solución a lo largo de una recta del dominio con elementos constantes. . . . . . . . . 69
6.4. Solución en los nodos del contorno con elementos lineales. . . . . . . . . . . . . . 70
6.5. Solución a lo largo de una recta del dominio con elementos lineales. . . . . . . . . . 70
6.6. Solución en los nodos del contorno con elementos cuadráticos. . . . . . . . . . . . . 71
6.7. Solución a lo largo de una recta del dominio con elementos cuadráticos. . . . . . . . 71
6.8. Solución en los nodos del contorno con elementos constantes. . . . . . . . . . . . . 72
6.9. Solución a lo largo de una recta del dominio con elementos constantes. . . . . . . . . 73
6.10. Solución en los nodos del contorno con elementos lineales. . . . . . . . . . . . . . 73
6.11. Solución a lo largo de una recta del dominio con elementos lineales. . . . . . . . . . 74
6.12. Solución en los nodos del contorno con elementos cuadráticos. . . . . . . . . . . . . 74
6.13. Solución a lo largo de una recta del dominio con elementos cuadráticos. . . . . . . . 75
6.14. Solución en los nodos del contorno, para la primera frecuencia con elementos lineales. . 76
6.15. Solución a lo largo de una recta del dominio infinito, para la primera frecuencia con
elementos lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.16. Solución sin usar chief points (desestabilizada) en los nodos del contorno, para la segunda
frecuencia con elementos lineales. . . . . . . . . . . . . . . . . . . . . . . . . 77
6.17. Solución sin usar chief points (desestabilizada) a lo largo de una recta del dominio infi-
nito, para la segunda frecuencia con elementos lineales. . . . . . . . . . . . . . . . 77
6.18. Solución estabilizada en el contorno, (empleo del Método de los Chief Points ) para la
segunda frecuencia con elementos lineales. . . . . . . . . . . . . . . . . . . . . 78
6.19. Solución estabilizada en el dominio, (empleo del Método de los Chief Points) para la
segunda frecuencia con elementos lineales. . . . . . . . . . . . . . . . . . . . . 78
6.20. Solución en el contorno cuando hay una onda plana incidente. . . . . . . . . . . . . 79
6.21. Solución en el dominio considerando solo la difracción, cuando hay una onda plana in-
cidente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.22. Solución en el dominio considerando reflexión y difracción, cuando hay una onda plana
incidente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.23. Solución en el contorno cuando hay una onda esférica incidente. . . . . . . . . . . . 81
6.24. Solución en el dominio considerando solo la difracción, cuando hay una onda esférica
incidente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.25. Solución en el dominio considerando reflexión y difracción, cuando hay una onda esféri-
ca incidente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.26. Geometrı́a de la cavidad 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6.27. Solución en los nodos del contorno del prisma rectangular con elementos lineales. . . . 84
6.28. Solución a lo largo de una recta del dominio, paralelela al eje y con elementos lineales. . 84
6.29. Vista 3D de la solución a lo largo de rectas del dominio, paralelelas al eje y con elementos
lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.30. Solución en los nodos del contorno del prisma rectangular con elementos lineales. . . . 85
6.31. Solución a lo largo de una recta del dominio, paralelela al eje y con elementos lineales. . 86
6.32. Vista 3D de la solución a lo largo de rectas del dominio, paralelelas al eje y con elementos
lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.33. Solución en los nodos del contorno de la circunferencia, para la primera frecuencia. . . 87
6.34. Solución a lo largo de una recta del dominio infinito, para la primera frecuencia. . . . . 88
6.35. Vista 3D de la solución a lo largo de una recta del dominio infinito, para la primera
frecuencia. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.36. Solución en el contorno, para la segunda frecuencia sin emplear chief points. . . . . . 89
6.37. Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia sin
emplear chief points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.38. Vista 3D de la solución a lo largo de una recta del dominio infinito, para la segunda
frecuencia sin emplear chief points. . . . . . . . . . . . . . . . . . . . . . . . . 90
6.39. Solución en el contorno, para la segunda frecuencia empleando 1 chief point. . . . . . 90
6.40. Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia emplean-
do 1 chief point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.41. Vista 3D de la solución a lo largo de una recta del dominio infinito, para la segunda
frecuencia empleando 1 chief point. . . . . . . . . . . . . . . . . . . . . . . . . 91
6.42. Solución en el contorno, para la segunda frecuencia empleando 10 chief points. . . . . 92
6.43. Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia emplean-
do 10 chief points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.46. Geometrı́a de la primera pantalla(reflejante). . . . . . . . . . . . . . . . . . . . . 94
6.47. Espectro de pérdida por inserción de la primera pantalla, obtenido mediante los códigos
adjuntados en los anexos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.48. Espectro de pérdida por inserción de la primera pantalla, sacado del libro de Maeso y
Aznárez [4] (2004). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.49. Geometrı́a de la segunda pantalla(absorbente). . . . . . . . . . . . . . . . . . . . 96
6.50. Espectro de pérdida por inserción de la segunda pantalla, obtenido mediante los códigos
adjuntados en los anexos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.51. Espectro de pérdida por inserción de la segunda pantalla, sacado del libro de Maeso y
Aznárez [4] (2004). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
H.1. Geometrı́a del problema de potencial estudiado. . . . . . . . . . . . . . . . . . . 212

Índice de cuadros
H.1. Resultados en el dominio con elementos constantes. . . . . . . . . . . . . 213

H.2. Resultados en el dominio con elementos lineales. . . . . . . . . . . . . . 213
H.3. Resultados en el dominio con elementos cuadráticos. . . . . . . . . . . . 213
H.4. Solución analı́tica de los resultados en el dominio. . . . . . . . . . . . . . 214
11
Capı́tulo 1
Introducción
1.1. Estudio y motivación

Este trabajo se inició con el objetivo de profundizar en los conocimientos de los méto-
dos matemáticos avanzados y la programación para la compresión de una de las principales
herramientas empleadas para la modelización en ingenierı́a como es el Método de los Ele-
mentos de Contorno y para la posterior implementación en Matlab de los planteamientos
teóricos estudiados. Se enmarca por tanto este trabajo, dentro del contexto de un proyecto
más amplio en el que el MEC no es la única herramienta empleada para la modelización,
sino que también se incluyen otras como puede ser el MEF. En cualquier caso, este docu-
mento está exlusivamente orientado a la formulación e implementación del MEC, con el
objetivo de abarcar la mayor cantidad posible de aplicaciones que ofrece dicho método.
1.2. Antecedentes
Este trabajo tiene como referencia fundamental cinco libros, incluidos en la biblio-
grafı́a, en los que se recogen los principios fundamentales para la implementación de los
problemas acústico y de potencial mediante el MEC. En algunos de dichos libros se reco-
gen además códigos en lenguage Fortran, que se adaptaron para poderse traducir a lenguaje
de Matlab. A partir de esos códigos ya en el nuevo lenguaje, se elaboraron otros códigos
más complejos que permitı́an un funcionamiento más general de los programas, ası́ co-
mo integraban de una forma más compacta y vinculada las distintas subrutinas, evitando
la necesidad de que el usuario del programa debiera realizar varias interacciones con el
12
Capı́tulo 1. Introducción 13
mismo, siendo tan solo necesaria la introducción de un fichero inicial con los datos de los
parámetros a emplear en la correspondiente ejecución del programa.
Dentro de los libros que este trabajo tiene como referencia, me gustarı́a remarcar la
importancia que ha tenido para el desarrollo de este documento, el libro de Dominguez
[2] (1993) en el que se aplica el Método de los Elementos de Contorno a problemas de
potencial, pues fueron sus planteamientos teóricos junto con el código adjuntado en len-
guaje Fortran nuestro punto de partida, para poder luego desarrollar la implementación
del problema acústico. Se han incluido dos anexos al final de este documento en los que
se ha hecho un resumen de las ideas para la formulación del MEC para los problemas de
potencial y se ha adaptado y traducido el código en lenguaje Fortran a lenguaje de Matlab.
1.3. Estructuración del trabajo

El trabajo se ha estructurado de la siguiente forma: tras este primer capı́tulo introduc-
torio, los dos siguientes capı́tulos se han dedicado a explicar los conceptos básicos y la
representación integral del problema acústico de una forma general, común a los espacios
bidimensional y tridimensional; los dos siguientes, es decir los capı́tulos 4 y 5, se han
dedicado a desarrollar la formulación e implementación del MEC para el problema 2D y
para el 3D respectivamente; el capı́tulo 6 recoge todos los ejemplos numéricos de aplica-
ción de los códigos descritos en los dos anteriores capı́tulos, incluyendo tanto ejemplos
de validación como ejemplos de aplicación práctica en la realidad; el capı́tulo 7 incluye
un resumen de los contenidos desarrollados a lo largo del trabajo, ası́ como las principales
conclusiones obtenidas en el mismo y una breve sección destinada a ubicar este documen-
to dentro de un marco de trabajos futuros. Por último se adjuntan una serie de anexos, en
los que se recogen los códigos explicados anteriormente ası́ como los de otra aplicación
del MEC como son los problemas de potencial, que como se ha señalado fue el punto de
partida de este trabajo.
Capı́tulo 2
Conceptos fundamentales de acústica
2.1. Introducción
En este capı́tulo se introducirán los fundamentos de la acústica, aplicables a problemas
con pequeñas perturbaciones de amplitud en fluidos homogéneos, no viscosos y en reposo.
Sobre estas hipótesis de acústica lineal, serán sobre las que se apoyará el desarrollo de todo
este trabajo. Por contra, los problemas con perturbaciones de amplitud finita, se incluyen
dentro de la acústica no lineal y quedan fuera del alcance de este documento.
2.2. Ecuaciones básicas del problema acústico
2.2.1. La ecuación de onda linealizada

Para el tratamiento de pequeñas perturbaciones de amplitud en problemas aplicados a
medios homogéneos, no viscosos y en reposo, las ecuaciones de la acústica lineal son las
siguientes:
∂P
+ ρO c2 ∇ ·V = 0
∂t
∂V
ρ0 = −∇P (2.1)
∂t
∂P 2 ∂P

+V · ∇Po = c +V · ∇Po
∂t ∂t
14
Capı́tulo 2. Conceptos fundamentales de acústica 15
La primera es la ecuación de conservación de la masa, la segunda es la ecuación de

Euler(ecuación de momento linealizado) y la tercer se obtiene de asumir que el proceso
es adiabático. Estas tres ecuaciones se combinan para dar lugar a la ecuación de onda
linealizada:
1 ∂ 2P

1
ρ0 ∇ · ∇P − 2 2 = 0 (2.2)
ρ0 c ∂t
Para medios homogéneos, la ecuación (2.2) se simplifica en:
1 ∂ 2P
∇2 P − =0 (2.3)
c2 ∂ t 2
2.2.2. Ondas armónicas planas: la ecuación de Helmholtz
Si la variación espacial de todas las variables acústicas del problema depende de una
sola coordenada, nos encontramos ante un caso de ondas planas. Si además las condicio-
nes de contorno y las fuentes internas de presión(en el caso de haberlas) presentan una
variación temporal de tipo de armónico con pulso ω , todas las variables acústicas serán
también armónicas con el mismo pulso. Bajo estas condiciones, se obtiene la ecuación de
onda en el dominio de la frecuencia:
d2 p 1
2
+ k2 p + 2 b = 0 (2.4)
dx c
donde p = p(x, ω ) es la función de respuesta compleja, x es una coordenada genérica
ω
y k es el número de onda k = . En ausencia de fuentes internas, la ecuación homogénea
c
que se obtiene es la conocida como ecuación de Helmholtz:
d2 p
+ k2 p = 0 (2.5)
dx2
2.2.3. Condiciones de contorno
Los problemas acústicos en el dominio de la frecuencia deben cumplir con las condi-
ciones impuestas en los contornos, que pueden ser de tres tipos:
p = p̄ : presión conocida(condición tipo Dirichlet)

∂p
= q̄ : flujo conocido(condición tipo Neumann)
∂n
∂p
= −ikβ p : condición mixta(condición tipo Robin)
∂n
En la última condición de contorno, β representa la admitancia del contorno, de forma
que si β = 0, el contorno es acústicamente rı́gido y si β = 1, toda la energı́a que incide
en el contorno, sale del dominio, por tanto implica que ambos medios tienen idénticas
propiedades.
2.3. Magnitudes para medir el sonido
2.3.1. Intensidad acústica

La intensidad acústica I de una onda, se define como el promedio de la potencia acústi-
ca transmitida a lo largo del tiempo a través de un área unitaria normal a la dirección
de propagación, en una posición fija del espacio. En el sistema internacional se mide en
(W /m2 ). Para nuestro caso particular de ondas planas en régimen armónico, tiene la si-
guiente expresión:
Z T
1 1
I= p2 dt (2.6)
ρ0 c T 0
2.3.2. Impedancia acústica especı́fica

La impedancia acústica especı́fica se define como el cociente entre la presión acústica
en un medio y la velocidad de partı́cula asociada. Para ondas planas, este cociente es un
valor real, cuyo signo depende del sentido de la propagación:
z = ±ρ0 c (2.7)
2.3.3. Nivel de presión sonora

El oı́do humano responde de manera lineal a la energı́a que le llega, la cual es propor-
cional al cuadrado de la presión acústica.
El rango de presión detectable por el oı́do humano es bastante amplio,(entre 20µ Pa y

60 Pa)por lo que se debe usar una escala comprimida. Se utiliza la escala logarı́tmica con
un factor de 10, definiéndose de esta forma el decibelio(dB).
El nivel de presión sonora p se dice que es L p decibelios(dB) mayor o menor que una
determinada presión de referencia pre f según la siguiente expresión:
< p2 >
L p = 10 · log10 (dB) (2.8)
< p2re f >
donde hp2 i se refiere a la media cuadrática temporal.
2.3.4. Coeficiente de pérdida por inserción

En este documento se tratará el uso de pantallas acústicas para la atenuación de sonido.
El medir en decibelios la reducción obtenida por interponer una de dichas pantallas, es una
buena forma de evaluar la eficacia de dicha estrategia. Se denomina pérdida por inserción
de la pantalla IL a la diferencia de presión sonora entre antes y después de introducir la
medida correctora:
IL = SPLdespues − SPLantes
< p2despues > < p2antes >

IL = 10 · log10 − 10 · log10
< p2re f > < p2re f > (2.9)
< p2despues >
IL = 10 · log10 (dB)
< p2antes >
Si se pretende calcular el coeficiente de inserción global, conocidos los coeficientes

de inserción para cada frecuencia y el espectro en frecuencia de una fuente genérica, el
coeficiente de inserción para un tono puro, puede expresarse de la siguiente forma:
|Pi|2
IL = 10 · log10 (2.10)
|Pi∗ |
2
Introducimos ahora el coeficinete αi como el porcentaje de presión sonora aportado en

la fuente por la frecuencia i a la presión sonora total:
ILi
|Pi |2 10 10
αi = 2
= ILi (2.11)
|P| ∑i 10 10
Sustituı́mos lo anterior en la definición de IL total y operando llegamos a:
!
ILi
IL = 10 · log10 ∑ αi · 10 10 (dB) (2.12)
i
Esta expresión permite obtener el coeficiente IL medio espectral a partir de los ILi de
cada componente del espectro.
Capı́tulo 3
Representación integral para problemas

escalares
3.1. Relación de reciprocidad

Partimos de la ecuación de Helmholtz, que gobierna la propagación de ondas de pre-
sión en un medio homogéneo, no viscoso y lineal en el deominio de la frecuencia:
1
∇2 p + k2 p +
b=0 (3.1)
c2
Ahora buscamos una solución del problema que anule la integral de volumen del resi-
duo ponderado por un campo virtual de presiones p∗ , a la que por el momento solo se le
exige que pueda ser derivada.Obtenemos:
Z
2 1 2∗
∇ p + k p + 2 b p dΩ = 0 (3.2)
Ω c
Integramos por partes dos veces y obtenemos lo siguiente:
Z Z Z
∗ 2 ∗ ∗ ∗
p ∇ p dΩ = (p p,i − p,i p),i dΩ + p∇2 p dΩ (3.3)
Ω Ω Ω
Ahora aplicamos el teorema de la divergencia sobre la primera integral de dominio del
segundo miembro:
Z Z Z
∗ 2 ∗ ∗ ∗
p ∇ p dΩ = (p p,i − p,i p)ni dΓ + p∇2 p dΩ (3.4)
Ω Γ Ω
donde Γ representa el contorno del dominio bajo estudio y ni las componentes del
19
Capı́tulo 3. Representación integral para problemas escalares 20
vector normal unitario en puntos del mismo. Teniendo en cuenta además que:
∗
∂p ∂p
= p,i ni y = p,i ni (3.5)
∂n ∂n
podemos escribir (3.4) del siguiente modo:
∗
∂p ∗ ∂p
Z Z Z Z
∗ 2 ∗
p ∇ p dΩ = p dΓ − p dΓ + p∇2 p dΩ (3.6)
Ω Γ ∂n Γ ∂n Ω
Sustituyendo esta expresión en el planteamiento residual de 3.2:
∗
∂p ∗ ∂p 1 ∗
Z Z Z Z
2 ∗ 2 ∗
p dΓ − p dΓ + ∇ p +k p pdΩ + 2
bp dΩ = 0 (3.7)
Γ ∂n Γ ∂n Ω Ωc
∗
Ahora introducimos una segunda exigencia a la función p : deberá verificar la ecuación
de gobierno del problema real en su formulación genérica. De esta forma:
∗ ∗ 1 ∗
∇2 p + k2 p = − b (3.8)
c2
Introducimos la ecuación anterior en 3.7 y tras reordenar, obtenemos finalmente:
∗
∂p 1 ∗ ∂p ∗ 1 ∗
Z Z Z Z
p dΓ + 2
b pdΩ = p dΓ + 2
bp dΩ (3.9)
Γ ∂n Ωc Γ ∂n Ωc
Esta ecuación es la relación de reciprocidad buscada y constituye el punto de partida
para la aplicación del MEC.
3.2. Solución fundamental

En esta sección el objetivo es convertir la relación de reciprocidad 3.9 en una igual-
dad integral que afecte sólo a los contornos del dominio en estudio. Para ello, la función
de ponderación elegida deberá verificar (3.8). A dicha función de ponderación la deno-
minaremos solución fundamental. Partiremos por tanto de la ecuacin 3.9, para el caso
simplificado en que no existen fuentes internas:
∗
∂p 1 ∗ ∂p ∗
Z Z Z
p dΓ + 2
b pdΩ = p dΓ (3.10)
Γ ∂n Ωc Γ ∂n
∗
Consideramos la función de ponderación p como el campo de presiones generado por
una función delta de Dirac aplicada en una región virtual de iguales propiedades que la
analizada en el problema real, en la que no se tienen en cuenta las condiciones de contorno
∗
de esta. En consecuencia, según (3.8) p solo tendrá que verificar la ecuación de gobierno
del problema:
∗ ∗
∇2 p + k2 p + δ (x − xi ) = 0 (3.11)
La función delta de Dirac es singular en x = xi y vale 0 en el resto de puntos del

dominio. Su integral de volumen extendida a todo el dominio es igual a la unidad. Ma-
temáticamente, puede interpretarse como el lı́mite de una función impulso rectangular
cuando la dimensión de su base tiende a cero. Fı́sicamente, representa una fuente puntual
armónica que pulsa a la misma frecuencia que el problema real, aplicada en un punto de
coordenadas xi , que en un principio consideramos que es interior al dominio de estudio.
Teniendo en cuenta estas propiedades, la integral de volumen de 3.10 puede escribirse as:
1 ∗
Z Z
2
b pdΩ = δ (x − xi )pdΩ = pi (3.12)
Ωc Ω
Con lo que la relación de reciprocidad queda de la siguiente forma:
∗
∂p ∂p ∗
Z Z
pi + p dΓ = p dΓ (3.13)
Γ ∂n Γ ∂n
De esta forma, (3.13) permite obtener, mediante la solución fundamental, el valor de
la presión en cualquier punto interno del dominio, conocido su valor y el de su derivada
en el contorno del mismo.
Γ
r r’
∞ Ω i ∞
i’
∞
Figura 3.1: Ecuación Integral y Solución Fundamental.
La solución fundamental (correspondiente a la solución de 3.11) para problemas armóni-

cos escalares bi y tridimensionales, se muestra a continuación:
∗ 1
p (k, r) = K0 (ikr) problemas 2D
2π
∗ 1 −ikr (3.14)
p (k, r) = e problemas 3D
4π r
donde r es la distancia entre el punto i de aplicación de la fuente puntual(punto de

∗
colocación) y el punto donde se quiere calcular el valor de la variable p (punto de obser-
vación). En el problema bidimensional, K0 representa la función modificada de Bessel de
segunda especie y orden cero, cuyo argumento es ikr.
También podemos colocar la fuente de presión en un punto exterior al dominio en es-

tudio. En este caso se anula la integral de volumen correspondiente de (3.10) y la relación
de reciprocidad queda ası́:
∗
∂p ∂p ∗
Z Z
p dΓ = p dΓ (3.15)
Γ ∂n Γ ∂n
Teniendo en cuenta todo lo anterior e incluyendo además la posibilidad de una fuente
de presión interna (delta de Dirac) en un punto x0 , la relación de reciprocidad queda de la
siguiente forma:
∗
∂p ∂p ∗
Z Z
∗
pi + p dΓ = p0 + p dΓ (3.16)
Γ ∂n Γ ∂n
Esta será la expresión que emplearemos como punto de partida en la formulación del
MEC.
3.3. Igualdad integral en el contorno

En esta sección se pretende obtener ahora una igualdad integral que afecte solo a pun-
tos del contornos, lo cual requiere de la aplicación de la fuente en puntos del mismo. La
dificultad se encuentra en las singularidades que aparecen en los integrandos de (3.16)
cuando situamos el punto xi sobre el contorno Γ del problema.
Γε n
ε
i Γ−Γ
ε
Ω
Figura 3.2: Igualdad integral en el contorno.
La estrategia que permite esquivar la singularidad a la hora de evaluar las integrales en

el contorno, consiste en recrecer infinitisimalmente el dominio en el contorno del punto
de colocación mediante un contorno circular de radio ε . Con ello, el punto de colocación
es interno al dominio y por tanto se puede aplicar (3.16). Se presentan a continuación las
integrales de contorno como la suma de otras dos extendidas a los contornos Γ − Γε y Γε
respectivamente.
∗ ∗
∂p ∂p ∂p ∗ ∂p ∗
Z Z Z Z
∗
pi + p dΓ + p dΓ = p0 + p dΓ + p dΓ (3.17)
Γ−Γε ∂n Γε ∂n Γ−Γε ∂n Γε ∂n
Ahora hay que conseguir que esta igualdad integral implique sólo variables del con-
torno. Para ello será necesario estudiar el comportamiento de estas integrales cuando ε
tiende a cero. Las primeras integrales no presentan problemas ya que el contorno sobre el
que se extienden no incluye singularidad. Por otra parte, las extendidas al contorno infini-
tesimal Γε están perfectamente acotadas en el lı́mite.
Teniendo en cuenta lo anterior, para las integrales extendidas al contorno diferencial,

distinguiremos entre aquella que incorpora en su integrando la solución fundamental en
∗
presiones p y aquella que tiene que ver con su derivada respecto de la normal.
Z α
∂p ∗ ∂p 1
Z
lı́m p dΓ = lı́m K0 (ikε )ε d α (3.18)
ε →0 Γε ∂n ε →0 0 ∂ n 2π
Ω ε dα
ε
dα
i
θ
Figura 3.3: Término libre.
Teniendo en cuenta que ε es constante, y que en el lı́mite la variable de campo es

constante y con el valor que adopta en el punto de colocación, podemos reescribir 3.18
como sigue:
∂p ∗ ∂p α
Z
lı́m p dΓ = ( )i lı́m [K0 (ikε )ε ] (3.19)
ε →0 Γε ∂n ∂ n 2π ε →0
donde α es el ángulo exterior de la esquina en radianes. Teniendo además en cuenta la
expresión de K0 en el lı́mite, podemos concluir.

∂p ∗ ∂p α 1
Z
lı́m p dΓ = lı́m Ln( )ε = 0 (3.20)
ε →0 Γε ∂n ∂n i 2π ε →0 ε
Por otro lado, la integral que incorpora la derivada de la solución fundamental tiene un
comportamiento diferente. Realizando un proceso análogo podemos escribir:
∗ Z α
∂p (−ik) ∂r
Z
lı́m p dΓ = lı́m p K1 (ikε ) ε dα (3.21)
ε →0 Γε ∂n ε →0 0 2π ∂n
Para el contorno cirular Γε se cumple además:
∂r
= r,i ni = 1 (3.22)
∂n
Teniendo en cuenta la expresión de K1 (ikε ) junto con lo anterior, podemos desarrollar

(3.21) del siguiente modo:
∗
∂p α α

1
Z
lı́m p dΓ = pi (−ik) lı́m ε = − pi (3.23)
ε →0 Γε ∂n 2π ε →0 ikε 2π
Para terminar, sustituı́mos (3.20) y (3.23) en (3.17) y obtenemos:
∗
∂p ∂p ∗
Z Z
∗
ci pi + p dΓ = p0 + p dΓ (3.24)
Γ ∂n Γ ∂n
donde:
α θ
=
ci = 1 − (3.25)
2π 2π
En la expresión anterior θ representa el ángulo interno de la esquina en radianes. El
término adicional (ci ) a considerar en la fomulación integral en el contorno, recibe el
nombre de término libre. Depende exclusivamente de la geometrı́a del contorno y para
contornos suaves su valor es ci = 0,5.
3.4. Condición de radiación en el infinito

Para el problema exterior, la solución tiene que cumplir las condiciones de radiación
en el infinito, condición de radiación de Sommerfeld, que para problemas n-dimensionales
se expresa matemáticamente como sigue:
∂p

n−1
lı́m r 2 + ikp = 0 (3.26)
r→∞ ∂r
Particularizando la expresión anterior para el problema 3-D en un dominio infinito y
con una fuente puntual, obtenemos:
∗
!
∂p ∗ ∗
lı́m r + ikp = lı́m (p ) = 0 (3.27)
r→∞ ∂r r→∞
Estas condiciones suelen expresarse como dos ecuaciones por separado, de forma que
se obtiene:
∗
!
∂p ∗
lı́m r + ikp =0
r→∞ ∂r
(3.28)
∗
lı́m (p ) = 0
r→∞
donde la primera subecuación representa la condición de radiación y la segunda es

la denominada ecuación de regularidad. Esta última está vinculada al amortiguamiento
por radiación propio de ondas bi y tridimensionales. Se puede obeservar por tanto que la
solución fundamental propuesta del problema cumple las condiciones de radiación y re-
gularidad en el infinito.
Concluiremos este capı́tulo reiterando la importancia de la consideración de las con-

diciones de radiación para la obtención de la única solución posible en problemas con
dominio infinito.
3.5. Igualdad para el problema exterior

El problema exterior es muy habitual en acústica. En él, la región estudiada es in-
finita y los contornos son interiores a la misma: son huecos en el interior del dominio
con condiciones de contorno conocidas. En la formulación integral del problema exterior,
las integrales de superficies se extienden a los contornos interiores del modelo. Para su
tratamiento, consideraremos el problema exterior como un caso particular del problema
interior y aplicaremos la igualdad integral (3.24). Cerraremos artificialmente el dominio
con una circunferencia para los problemas 2-D y con una esfera para los problemas 3-D,
de radio r → ∞ y centrada en el punto de colocación donde se quiere aplicar la igualdad
integral.
n
a b
Γ Γ
Ω Ω
n
Figura 3.4: (a)Problema interior.(b)Problema Exterior.

Γ∞
Γ i
n
*x0
Ω
Figura 3.5: Formulación integral del problema exterior.
La igualdad integral en el dominio encerrado por el contorno Γ∞ podemos escribirla

ası́:
∗ ∗
∂p ∂p ∂p ∗ ∂p ∗
Z Z Z Z
∗
ci pi + p dΓ + p dΓ = p0 + p dΓ + p dΓ (3.29)
Γ ∂n Γ∞ ∂r Γ ∂n Γ∞ ∂ r
donde se reemplaza ∂ /∂ n por ∂ /∂ r en el contorno Γ∞ . La solución buscada para el

problema tendrá que cumplir las condiciones de radiación en puntos de este contorno,
por lo que dicha solución se aproxima a una onda cilı́ndrica(2-D) o esférica (3-D) que se
propaga en la dirección de los r crecientes. Por tanto en puntos de Γ, se puede escribir:
∂p
= −ikp (3.30)
∂r
Sustituyendo esta última expresión en (3.29), se puede escribir:
∗ ∗
!
∂p ∂p ∂p ∗
Z Z Z
∗
ci pi + p dΓ + p + ikp dΓ = p0 + p dΓ (3.31)
Γ ∂n Γ∞ ∂r Γ ∂n
Concluiremos remarcando el interés en que la solución fundamental del problema ex-

terior verifique las condiciones de radiación en el infinito, ya que la integral de contorno
extendida a Γ∞ se anula y la formulación integral se extiende tan solo a un contorno Γ
interior. En cambion en el problema interior, no se necesita este requerimiento. Puede ob-
servarse por tanto la importante ventaja que supone la aplicación del MEC frente a los
métodos del dominio como el MEF, ya que al tratar con dominios de naturaleza infinita,
los métodos del dominio presentan el gran inconveniente de que no se puede extender la
discretización del medio hasta el infinito.
3.6. Problemas de scatter

Los problemas de scatter, también conocidos como scattering problems o problemas
con obstáculos constituyen un problema tı́pico en el que una onda incidente pI impacta
con un obstáculo V ′ que constituye el hueco que queda fuera del dominio de un problema
exterior. La onda dispersada la denotaremos por pS , de forma que la onda de presión total
será la suma de ambas ondas:
p = pI + pS (3.32)
V
Pi
V’ Γ
Ps
Figura 3.6: Problema de scatter con una onda incidente Pi.
En cuanto a pI , puede interpretarse también como la presión en ausencia del obstáculo

V ′ . Si la onda incidente es una onda plana:
pI = Ae−ik(α x+β y+γ z) (3.33)
donde A es la amplitud de la onda plana incidente y α , β y γ son los cosenos directores

de la dirección de incidencia. Si la onda incidente es una onda esférica:
e−ikR
pI = B (3.34)
R
donde la constante B está relacionada con la fuerza de la fuente puntual y R es la dis-

tancia desde la fuente hasta el punto considerado. Cabe destacar que tanto pI como pS
deben satisfacer la ecuación de Helmholtz, pero solo pS tiene que satisfacer la condición
de radiación de Sommerfeld.
Partimos de la ecuación integral de Helmholtz para el problema exterior, sin incluir

el término asociado a las fuentes internas de presión, por simplicidad de las expresiones,
aunque de considerarlas, el planteamiento serı́a el mismo incluyendo dicho término. Apli-
camos la ecuación a la onda pS :
∗
S∂p ∂ pS ∗
Z Z
ci pSi + p dΓ = + p dΓ (3.35)
Γ ∂n Γ ∂n
A continuación, aplicamos la ecuación integral de Helmholtz para el problema interior,

a la onda incidente pI , en el dominio interior V ′ , donde no es necesario cumplir la ecuación
de radiación de Sommerfeld. Si mantenemos el mismo sentido de la normal que en (3.35),
obtenemos:
∗
I∂p ∂ pI ∗
Z Z
c0i pIi + p dΓ = − p dΓ (3.36)
Γ ∂n Γ ∂n
Si restamos (3.39) a (3.38), hacemos uso de (3.31) y tenemos en cuenta que ci +c0i = 1,
llegamos a:
∗
∂p ∂p ∗
Z Z
ci pi + dΓ = + p dΓ + pIi (3.37)
Γ ∂n Γ ∂n
Esta es la ecuación integral para problemas con obstáculos. Puede observarse que la
única diferencia con respecto a la ecuación integral original es la incorporación del término
libre asociado a la onda incidente.
3.7. No unicidad de la solución
A pesar de las ventajas señaladas en los apartados anteriores, que presenta el MEC
para problemas exteriores, la aplicación directa del método presenta también algún incon-
veniente que requiere ser tratado para poder subsanarse.
El principal incoveniente es que la ecuación integral de Helmholtz no tiene solución

única para determinadas frecuencias asociadas con el correspondiente problema interior.
Explicaremos este fenómeno partiendo de dicha ecuación integral de Helmholtz y teniendo
como referencia el libro de Wu [5] (2000). Por simplicidad en la formulación asumiremos
que el contorno es suave, que no hay fuentes de presión internas y que no hay una onda
incidente, aunque lo que se va a desarrollar a continuación es perfectamente extrapolable
a problemas con cualquier tipo de contorno, con fuentes internas y con obstáculos.
La ecuación integral de Helmholtz para problemas exteriores (3.31), se simplica de

la siguiente forma de acuerdo a las consideraciones hechas y después de reordenar los
términos:
∗
∂p ∗ 1 ∂p
Z Z
− p dΓ = − pi − p dΓ (3.38)
Γ ∂n 2 Γ ∂n
Si consideramos además un problema auxiliar interior definido en el hueco que ocupa

la parte que está fuera del dominio del problema exterior y mantenemos el mismo sentido
de la normal que en dicho problema exterior, la ecuación integral para el problema auxiliar
es la siguiente:
∗
∂p ∗ 1 ∂p
Z Z
− p dΓ = pi − p dΓ (3.39)
Γ ∂n 2 Γ ∂n
Observando (3.38) y (3.39), nos damos cuenta de que ambas ecuaciones tienen el mis-
mo miembro de la izquierda. Debido a que (3.39) es para un problema exterior, se pro-
ducirá un fen’meno de resonancia para un conjunto de frecuencias (frecuencias naturales
del problema en estudio), es decir, las matrices de (3.38) se volverán singulares, lo cual
afectará a las matrices de (3.39), que también experimentarán singularidades.
3.8. Uso del Método de los Chief Points

Para solventar el problema planteado en el apartado anterior, sobre la no unicidad de
la solución para determinadas frecuencias, se propone un método conocido con el nombre
de Método de los Chief Points, cuyo objetivo es precisamente el de estabilizar la solución
para todas las frecuencias.
Para ello, este método utiliza la ecuación integral de Helmholtz en algunos puntos del
hueco interior que queda fuera del dominio del problema exterior en cuestión, además de
para su uso habitual en el contorno. Estos puntos son los que se conocen como chief points.
Lo que se pretende es generar ecuaciones adicionales de restricción. Estas ecuaciones,
conocidas como chief equations tienen la misión de hacer cumplir el valor nulo de la
presión en los chief points, lo que supone dichas restricciones adicionales,(a las ecuaciones
asociadas a las integrales en el contorno) que ayudan a estabilizar la solución. Como ci vale
cero fuera del dominio en estudio, las chief equations, adoptan la siguiente forma:
∗
∂p ∂p ∗
Z Z
p dΓ − p dΓ = 0 (3.40)
Γ ∂n Γ ∂n
Al incluir estas ecuaciones adicionales, se forma un sistema de ecuaciones sobredeter-
minado, que se resuelve por mı́nimos cuadrados.
3.9. Solución fundamental para el semiespacio

En este apartado presentaremos una nueva modificaicón de la solución fundamental
de forma que también cumpla alguna de las condiciones de contorno del problema real,
con el objetivo de no considerar las integrales extendidas al contorno cuya condición es
satisfecha por la solución fundamental.
Esta necesidad está vinculada a una nueva forma de discretizar el contorno que se
utilizará en el apartado dedicado al cálculo de panatallas acústicas. Esta nueva forma de
discretizar tiene como objetivo el realizar dicha discretización tan solo en el contorno que
delimita la pantalla.
En un principio, la aplicación de la igualdad integral en los términos conocidos, im-

plica integrales a los tres contornos del problema (Γ∞ , Γs y Γ p ). Aunque sabemos que las
integrales extendidas a Γ∞ desaparecen si la solución fundamental cumple las condiciones
de radiación en el infinito. Teniendo esto en cuenta, podemos escribir (3.16) en términos de
cada uno de los contornos que definen el problema, para un punto de colocación genérico
i interior al dominio, de la siguiente forma:
∗ ∗
∂p ∂p ∂p ∗ ∂p ∗
Z Z Z Z
∗
pi + p dΓ + p dΓ = p + p dΓ + p dΓ (3.41)
Γs ∂ n Γp ∂ n Γs ∂ n Γp ∂ n
Ω
Γp(βp)
fuente(x0)
Γs(βs)
*
Figura 3.7: Semiespacio 2D con la pantalla.

∗
Elegiremos una función p que verifique la condición de contorno Γs
∗
∂p ∂p ∗
= −ikβs p y = −ikβs p en Γs (3.42)
∂n ∂n
El incluir la condición anterior, permite la cancelación mutua de las integrales exten-
didas al contorno del semiespacio Γs . Se obtiene:
∗
∂p ∂p ∗
Z Z
∗
pi + p dΓ = p + p dΓ (3.43)
Γp ∂ n Γp ∂ n
De esta forma la identidad integral se extiende únicamente a los contornos del obstácu-
lo(pantalla) en cuestión, con lo que se ha conseguido el objetivo que se pretendı́a de reali-
zar la discretización tan solo en el contorno que delimita la pantalla.
Para un punto de colocación sobre el contorno del obstáculo, la identidad integral es la

siguiente:
∗
!
∂p
Z
∗ ∗
ci pi = p0 − + ikβ p p pdΓ (3.44)
Γp ∂n
Mostraremos a continuación la expresión matemática de la solución fundamental a la

que nos hemos estado refiriendo. El caso más sencillo es aquel en el que el contorno del
semiespacio es completamente reflejante (βs = 0). Esta condición se consigue consideran-
do además de la fuente aplicada en el punto de colocación donde se escribe la identidad
integral, otra fuente simétrica respecto del contorno del semiespacio, pulsando a la misma
frecuencia. Considerando esto, llegamos a:
∗ 1 1
p (k, r) = K0 (ikr) + K0 (ikr̄) problemas 2D
2π 2π
∗

1 1 −ikr 1 −ikr̄
(3.45)
p (k, r) = e + e problemas 3D
4π r r̄
donde r y r̄ representan la distancias de la fuente y la imagen al punto de observación

respectivamente.
(x,y)
punto de r punto de
colocación observación
(xi,yi)
y r
x Γs(βs)
(xi,−yi)
imagen
Figura 3.8: Representación de fuente e imagen para la solución fundamental.
Para el caso en que βs 6= 0 la expresión de la solución fundamental se complica más.

Mostraremos la forma de abordarla para el caso bidimensional, que es el que se utili-
zará después en las aplicaciones. Se basa para dicho caso en la ecuación de propagación
en campo lejano y se expresa de la siguiente manera:
∗ 1
p (k, r) = [K0 (ikr) + K0 (ikr̄) + Φ(ikr̄)] (3.46)
2π
El término adicional añadido respecto a (3.37), es el siguiente:
√
r
a
Φ(ikr̄) = π ω (i ikr̄a)[V (γ ) − 1]e−ikr̄ (3.47)
2
donde
y + yi
q
a = 1 + γβs − (1 − βs2 )(1 − γ 2 ) ; γ = cosθ =
r̄
√ √ (3.48)
γ − βs
V (γ ) = ; ω (i ikr̄a) = eikr̄a er f c( ikr̄a)
γ + βs
donde a su vez er f c(z) es la función complementaria de error con argumento comple-

joy V (γ ) es un término que representa el coeficiente de refleción para ondas planas que
inciden con ángulo θ sobre el contorno en estudio.
El primero término de (3.46) representa la componente directa de la propagación desde

el emisor y el segundo y el tercer término engloban la componente debida a las reflexiones
en el contorno del semiespacio.
Capı́tulo 4
Formulación e implementación del

MEC en problemas acústicos 2D
4.1. Introducción y descripción general de la metodologı́a

En este apartado se desarrollarán los procedimientos necesarios para realizar el plan-
teamiento numérico de la identidad integral (3.44) en problemas bidimensionales. Se
usarán estrategias análogas a las que se utilizan en la fomulación del Método de los Ele-
mentos Finitos y se aproximará el comportamiento de las variables de forma localizada
en una serie de elementos en que se divide el contorno del problema. Dentro de cada ele-
mento la variable se interpolará en función del valor que adopta en los nodos del mismo.
Se emplearán elementos constantes, lineales y cuadráticos, cuyas caracterı́sticas se ex-
pondrán más adelante. Una de las diferencias fundamentales en la discretización del MEC
frente a la del MEF es la de ser necesario aproximar tanto la variable fundamental como
su derivada, ya que ambas se encuentran en la formulación integral inicial (3.44).
4.2. Ecuación integral de Helmholtz para problemas bi-

dimensionales
Tal como se ha citado en la sección anterior, el punto de partida para desarrollar el
planteamiento numérico es la identidad integral (3.44). Aunque es cierto que dicha ecua-
ción hace referencia al caso particular que se aplicará al cálculo de pantallas, donde el
contorno a discretizar se reduce a Γ p . Por ello, a continuación expondremos un caso más
35
Capı́tulo 4. Formulación e implementación del MEC en problemas acústicos 2D 36
general(partiendo del problema interior y pasando luego al exterior) en el que partiremos

de todo el contorno Γ.
Para el problema interior definido en un dominio finito 2D, Ω la ecuación integral en

el contorno es la siguiente:
∗
!
∂p
Z
∗ ∗
c0i pi = p0 − + ikβ p p pdΓ (4.1)
Γ ∂n
En esta expresión, la normal está apuntando hacia fuera del dominio interior y el coe-
ficiente c0i vale 1 en el interior del dominio, 1/2 en los puntos i que se encuentren en
contornos suaves y 0, fuera del dominio. Para cualquier otro caso, habrá que emplear la
expresión (3.25).
Para el problema exterior definido en un dominio Ω infinito 2D la ecuación integral en

el contorno es la siguiente:
∗
!
∂p
Z
∗ ∗
ci pi = p0 − + ikβ p p p + pI dΓ (4.2)
Γ ∂n
En esta expresión pI representa una onda incidente en el caso de haberla, y la normal

apunta en sentido contrario al caso anterior, es decir, hacia fuera del dominio exterior. El
valor del coeficiente ci es el complementario sobre 1 del de c0i , de forma que:
ci = 1 − c0i (4.3)
4.3. Discretización del contorno
Para resolver la ecuación integral de Helmholtz, en primer lugar habrá que discretizar
el contorno en un número determinado de elementos curvilı́neos. Dependiendo del tipo de
elemento, este tendrá un número dado de nodos. En cualquier caso la geometrı́a de cada
elemento se representa interpolando entre los nodos, es decir:
n
x = ∑ xi Ni (ξ )
i=1
n (4.4)
y = ∑ yi Ni (ξ )
i=1
donde xi e yi son las coordenadas de los nodos, Ni (ξ ) son las funciones de forma
definidas en un elemento con coordenadas locales −1 ≤ ξ ≤ 1 y n es el número de nodos
del elemento.
elementos elementos elementos

constantes lineales cuadráticos
elemento elemento elemento
nodo nodos nodos
Figura 4.1: Tipologı́as de elementos empleadas para discretizar.
4.3.1. Elementos constantes

Para el caso de elementos constantes, se emplea un único nodo, situado en el centro
del elemento. Hay por tanto una única función de forma que será constante a lo largo de
todo el elemento:
φ1 = 1 (4.5)
nodo
φ1=1
ξ=−1 ξ=0 ξ=+1
Figura 4.2: Función de forma constante.
4.3.2. Elementos lineales

Para el caso de elementos lineales, se emplean dos nodos, uno situado en cada ex-
tremo del elemento. Hay por tanto dos funciones de forma, que adoptan las siguientes
expresiones:
φ1 = 21 (1 − ξ )
(4.6)
φ2 = 21 (1 + ξ )
φ1=0.5(1−ξ)
ξ=−1 ξ=0 ξ=+1
nodos
φ2=0.5(1+ξ)
ξ=−1 ξ=0 ξ=+1
Figura 4.3: Funciones de forma lineales.

4.3.3. Elementos cuadráticos

Para el caso de elementos cuadráticos, se emplean tres nodos, uno en cada extremo
del elemento y otro en el centro del mismo. Hay en consecuencia tres funciones de forma,
según las siguientes expresiones:
N1 = 21 ξ (ξ − 1)
N2 = 21 (1 + ξ )(1 − ξ ) (4.7)
N3 = 21 ξ (ξ + 1)
φ1=0.5ξ(ξ−1)
ξ=−1 ξ=0 ξ=+1
φ2=(1+ξ)(1−ξ)
ξ=−1 ξ=0 ξ=+1
φ3=0.5ξ(ξ+1)
ξ=0 ξ=+1
Figura 4.4: Funciones de forma cuadráticas.
4.3.4. Desarrollo del plantemiento numérico de la identidad integral

Independientemente del tipo de elementos seleccionado para discretizar, (y en conse-
cuencia del número de nodos de los mismos) la aplicación de la expresión de la identidad
integral en un nodo i genérico de la discretización se realiza de la misma forma, de mo-
do que se introduce la aproximación de la variable de acuerdo con el tipo de elementos
elegido. Se puede escribir por tanto:
NE N
∗ ij ij j
ci pi = p0 − ∑ ∑ (hk + ikβ j gk )pk (4.8)
j=1 k=1
donde NE representa el número de elementos en que se ha dividido el contorno y β j

la admitancia del contorno en cada elemento, siendo a su vez:
∗
∂p
Z Z
ij ∗
hk = φk dΓ j ; gikj = p φk dΓ j (4.9)
Γj ∂ n Γj
los denominados núcleos de integración de la solución fundamental y cuya evaluación

se realiza mediante procesos de cuadratura numérica estándar. Con este preceso hemos
transformado la igualdad integral en una ecuación algebraica donde las incógnitas son los
valores de presión en los nodos. Si aplicamos esta discretización utilizando como puntos
de colocación todos los nodos del contorno, se plantea un sistema de ecuaciones del tipo:
∗
(H + ikβ G)P = P0 (4.10)
donde H y G son matrices cudradas(NxN) siendo N el número de nodos de la dis-

cretización, al igual que en la expresión (5.8). Además, P es un vector que contiene los
∗
valores de la presión en los nodos(incógnitas) y P0 es otro vector donde se almacena el va-
lor obtenido de la solución fundamental en la fuente interna real de presión (x0 ) utilizando
como punto de colocación cada nodo del contorno. La resolución de este sistema permite
obtener por tanto los valores de presión y flujo en todos los nodos del contorno.
Una vez conocidos los valores en el contorno, se pueden obtener los correspondientes
valores en puntos interiores del dominio, aplicando la siguiente ecuación algebraica:
NE N
∗
pi = p0 − ∑ ∑ (hikj + ikβ j gikj )pkj (4.11)
j=1 k=1
donde ahora i pasa a ser un punto genérico del dominio donde se quiere calcular la
j
presión y pk la presión en los nodos del contorno, obtenida tras aplicar (4.10).
Cabe destacar por último que los términos de la matriz de coeficientes del sistema de
ecuaciones (4.10), son todos distintos de cero, es decir, la matriz es llena, y además es una
matriz no simétrica. Esto es un incoveniente que presenta el MEC frente al MEF, donde la
matriz de coeficientes del sistema de ecuaciones es una matriz en banda y simétrica, con
lo que se ahorran tanto recursos de memoria como tiempo de computación.
4.4. Utilización del Método de los Chief Points y solución

por mı́nimos cuadrados
Tal como se dijo en el capı́tulo dedicado a la representación integral del problema
acústico, se utilizará el Método de los Chief Points para generar una serie de ecuaciones
que impongan restricciones adicionales sobre el valor de la presión.
Aunque para que este método resulte eficaz, hay que seleccionar adecuadamente los
chief points, ya que si se colocan estos en puntos que caen en el interior del problema
interior asociado, carecen de utilidad, puesto que la presión en esa zona ya es nula de por
sı́ y no proporcionan por tanto un efecto restrictivo adicional.
Para nuestro problema exterior 2D, será necesario con la colocación de un solo chief
point en el interior del hueco que queda fuera del dominio, para poder estabilizar la solu-
ción. Este resultado se mostrará en el capı́tulo 6, donde se comparará la solución obtenida
al utilizar chief points con la obtenida al no utilizarlos.
4.5. Aplicación al estudio de pantallas de atenuación acústi-

ca
Un ejemplo de aplicación del MEC en Ingenierı́a Civil es el diseño de pantallas de
atenuación acústica para su implantación en carreteras. En este documento se incluirá un
código en Matlab traducido de Fortran del libro de (Maeso y Aznárez) [4] (2004) y adap-
tado a la estructura de código que se venı́a desarrollando anteriormente. Para ello, se par-
tirá del código exterior 2D programado anteriormente y se añadirán las modificaciones
necesarias para conseguir el resultado final obtenido en el libro citado.
Las modificaciones fundamentales son: la adición ya explicada de nuevos términos a

la solución fundamental, para conseguir una discretización más reducida del contorno; la
implementación de una nueva función que realice de forma independiente la integración
de los nodos de colocación coincidentes con los nodos del elemento integrado, de forma
que se aborde de una forma especı́fica la singularidad y por último la incorporación de una
serie funciones relaciones con un modelo para estudiar la capacidad de absorción de las
pantallas y del semiespacio conocido como Modelo de Delany y Bazley. Todo esto se es-
tudiará con más detalle al explicar las distintas funciones implementadas. A continuación
explicaremos las ideas básicas en que se fundamenta el modelo de Delany y Bazley.
El Modelo de Delany y Bazley propone una relación empı́rica para determinar la ad-
mitancia β de un contorno en el caso de materiales porosos, que son los que se utilizan
habitualmente como materiales absorbentes. Dicha relación es:
−0,75 −0,73
1 3 f 3 f
= 1 + 9,08 · 10 − i · 11,9 10 (4.12)
β σ σ
donde f representa la frecuencia de excitación (Hz) y σ la resistividad al flujo de aire
del material.
Aunque en muchas aplicaciones, se utilizan como recubrimiento de superficies refle-

jantes, materiales absorbentes de espesor relativamente pequeño. En este caso:
βe = β tanh(γ e) (4.13)
donde β es la admitancia obtenida mediante (4.12), e el espesor del recubrimiento (m)

y γ el número de onda asociado a la propagación en el medio absorbente. Existen también
en la bibliografı́a distintas fuentes que proporcionan relaciones empı́ricas para obtener es-
te último valor.
Para tener en cuenta que la energı́a puede incidir en el contorno en todas las direccio-
nes y no solo en la normal, se utiliza un coeficiente de absorción estadı́stico, propuesto
por Morse y Bolt(1944), que caracteriza el comportamiento absorbente de los materiales
utilizados como recubrimiento.La expresión de dicho coeficiente de absorción estadı́stico
es la siguiente:
8cosθ cosθ cos2θ ξ senθ

αst = 1− ln(1 + 2ξ cosθ + ξ 2 ) + tan−1 (4.14)
ξ ξ ξ senθ 1 + ξ cosθ
donde Z/(ρ + c) = R + iX , ξ 2 = R2 + X 2 yθ = tan−1(X /R)

Una vez obtenido el coeficiente estadı́stico de absorción de energı́a, se puede calcular
el valor de la admitancia a través de la siguiente expresión:
√
1 − 1 − α st
β= √ (4.15)
1 + 1 − α st
El valor obtenido mediante la expresión(4.15) se caracteriza por ser un valor real,

frente al obtenido mediante las expresiones (4.12) y (4.13) que en general será un valor
complejo.
4.6. Funciones implementadas

Para el problema acústico 2D, se han desarrollado una serie de programas que permi-
ten calcular los valores de presión y el flujo tanto en el contorno como en el dominio, de
acuerdo con las ecuaciones recogidas en las secciones anteriores de este capı́tulo. Se han
implementado 4 tipos de programas para el caso 2D: problema interior, problema exte-
rior, problema de scatter y problema de aplicación al cálculo de pantallas de atenuación
acústica. Todos ellos se han programado a su vez mediante elementos constantes, lineales
y cuadráticos.
Cada uno de los programas se compone a su vez de una serie de funciones y aunque
cada una de ellas tiene una serie de particularidades según el programa al que pertenece,
todos los programas se estructuran en una serie de funciones básicas que permiten reali-
zar las operaciones necesarias para obtener los resultados tanto en el contorno como en el
dominio. Dependiendo del programa,(de cada uno de los 4 tipos señalados en el párrafo
anterior) este podrá contener además otra serie de funciones complementarias que ayuden
a realizar cálculos auxiliares. La estructura básica de todos los programas es la siguiente:
fichero de entrada de datos, programa principal que gestiona al resto de funciones y las
llama en el orden adecuado para que realicen las operaciones correspondientes; función de
lectura de los datos del fichero de entrada(INPUT); función que realiza las integraciones
para los nodos del contorno y calcula los valores de presión y flujo en el mismo(GHMAT)
y por último una función que realiza las integraciones para puntos interiores al dominio,
calculando también los valores de presión y flujo(INTER).
A continuación empezaremos describiendo más en detalle esas funciones para el pro-

grama correspondiente al problema interior(son las únicas que forman dicho programa) y
luego pasaremos al resto de programas, describiendo las nuevas funciones que aparezcan
en cada uno de ellos, y comentando las diferencias que presentan las funciones básicas
con respecto a los programas anteriores.
Problema interior 2D
Fichero de entrada de datos
En él se introducen los valores de los parámetros que se utilizarán para la ejecución del
problema. Los datos son: número de nodos que se emplearán para la discretización, núme-
ro de puntos del dominio para el que se calcularán sus valores de presión y flujo, coorde-
nadas (x,y) de los nodos, una pareja de valores para cada nodo, de forma que el primer
valor será un 0 si la condición de contorno que se conoce del nodo es una presión y un
1, si la condición de contorno conocida es un flujo y el segundo valor de la pareja será el
valor para dicha condición de contorno y por último se incluyen también las coordenadas
(x,y) de los puntos del dominio.Cabe destacar que para el problema interior, el orden de
introducción de los nodos debe ser en sentido antihorario, de forma que la normal apunte
siempre hacia el exterior del dominio.
Función principal
Su misión es la de coordinar las distintas funciones; las va llamando en orden de forma que
cada una realice sus cálculos en el momento adecuado. En primer lugar llama a la función
de lectura de los datos del fichero de entrada, luego a la que realiza las integraciones en el
contorno y tras realizar una reordenación de las matrices G y H obtenidas, resuelve el sis-
tema de ecuaciones y dibuja los resultados obtenidos en el contorno. Posteriormente llama
a la función que realiza la integración para los puntos del dominio y dibuja también los
resultados obtenidos en el dominio. Esta función será muy parecida para los distintos tipos
de elementos(constantes, en los que la función se ha denominado ACCONBE; lineales, en
los que se ha denominado ACLINBE y cuadráticos, a la que se ha llamado ACQUABE.
La mayor precaución habrá que tenerla a la hora de introducir los nodos de los elementos,
pues como ya se ha explicado, a cada tipo de elemento le corresponde un número diferente
de nodos.
Función de lectura de datos del fichero de entrada
Esta función es la primera en ser llamada por el programa principal. La primera operación
que realiza es volcar todos los datos del fichero de entrada a un vector. Posteriormente
irá asignando el valor adecuado a cada una de las variables para después proporcionárse-
las al programa principal. Es por ello imprescindible introducir los valores en el orden
establecido en el fichero de entrada, para que la posterior asignación de esos valores a las
variables, sea correcta. Esta función también será muy parecida para los tres tipos de ele-
mentos. Los nombres de las funciones, siguiendo la misma nomenclatura que antes, serán:
INPTPC para elementos constantes, INPUTPL para elementes lineales INPUTQ para ele-
mentos cuadráticos.
Función de integración en el contorno
Esta función realizará la integración para todos los nodos del contorno, empleando una
cuadratura de Gauss de 6 puntos para los elementos constantes y lineales y de 10 puntos
para los elementos cuadráticos. En este programa se realizará de manera conjunta la in-
tegración para los nodos en que el punto de colocación coincide con el elemento. Esto es
debido a que la singularidad que se produce para el problema acústico en estos casos, no
es tan pronunciada como ocurre en otros problemas, como los de potencial. De esta forma
es posible utilizar el mismo tipo de integración para todos los nodos, sin que esto apenas
afecte al tiempo de computación(siempre que las mallas generadas no sean excesivamente
densas). Las funciones asociadas a los elementos constantes y lineales(GHMATPC y GH-
MATPL respectivamente), son muy parecidas entre sı́ y sin embargo difieren notablemente
de la asociada a elementos cuadráticos(GHMATPQ). Esta última tiene una mayor comple-
jidad tal como puede verse en los códigos adjuntados en los anexos, aunque es cierto que
los resultados que ofrece son también más precisos.
Función de integración en el dominio
Esta función es análoga a la anterior, de forma que pueden aplicarse las mismas ideas:
también se realiza la integración mediante 6 puntos de Gauss para elementos constantes
y lineales (funciones INTERPC e INTERPL respectivamente), y mediante 10 puntos de
Gauss (función INTERPQ). Se hará además de una forma conjunta para todos los nodos y
la complejidad para los distintos tipos de elementos es la misma que en el caso anterior.
La diferencia está en que para este caso no hay que calcular una matriz C, puesto que los
puntos del dominio tienen todos valor unidad para su término libre. Además, una vez ob-
tenidas las matrices G y H, en lugar de resolver un sistema de ecuaciones, lo que hay que
hacer es multiplicar las anteriores matrices por los valores de presión y flujo calculados en
el contorno.
Problema exterior 2D
En el problema exterior se emplean todas las funciones anteriores y además se añade

una nueva función para la integración en los chief points, que permite incluir ecuaciones
adicionales que estabilicen la solución a determinadas frecuencias. La otra gran diferencia
con respecto al problema interior, es que a la hora de introducir las coordenadas de los no-
dos en el fichero de entrada, estos se incluyen en sentido horario, de forma que la normal
también apunte hacia el exterior del dominio.
Función de integración de los chiefpoints
El objetivo de esta función(llamada CHIEFPOINTS) es realizar la integración en los ya

explicados chief points de forma que se generan ecuaciones adicionales que estabilicen la
solución. Estas ecuaciones se añaden apropiadamente a las matrices G y H .Dichas ecua-
ciones adicionales, no aportan términos distintos de 0 a la matriz C. Una forma de resolver
esto es añadir los términos distintos de cero a la matriz H después de haberle sumado C.
También se pueden añadir previamente los términos distintos de cero a H, las filas de ceros
correspondientes a C y posteriormente sumar ambas matrices. Esta última es la solución
por la que se ha optado.
Es importante hacer notar que el uso de esta función en el problema exterior es opcio-
nal. Esto está controlado por una pareja de valores k y nchi f que se dan al final del fichero
de entrada. Cuando k = 0, no se hace la llamada a la función CHIEFPOINTS, mientras
que en el caso k = 1 sı́ se hace. El número de puntos que se utilizarán entonces, será el que
indique la variable nchi f .
Problema de scatter
En este caso las funciones son análogas a las del problema exterior. Las variaciones se
encuentran en que tal como se explicaba en la sección dedicada a problemas con obstácu-
los, en este caso se incluirá la presencia de ondas incidentes, lo que implica un término
adicional en la identidad integral. Este término adopta una forma distinta según la onda
sea plana o esférica.
Se han considerado 4 argumentos de entrada en el programa principal(también puede

hacerse a través del fichero de entrada, pero se ha hecho en el programa principal para ex-
plicar dentro del mismo mediante comentarios la utilización de los argumentos de entrada)
que adquieren durante la ejecución del programa significados diferentes según la onda sea
plana o esférica. Para que la onda sea considerada como plana, se indica con un 1 en el
primer argumento de entrada. Esto implica que el segundo argumento será el ángulo que
forma el frente de avance de la onda con respecto al eje X, el tercero será el ángulo que
forma con respecto al eje Y(el incluir ambos permite controlar con más comodidad los
sentidos de propagación de la onda) y el cuarto será la amplitud de la onda. En cambio,
si se quiere que la onda sea esférica, se debe indicar con un 2 en el primer argumento de
entrada, lo que implica que el segundo sea considerado como la coordenada X de la fuente
de emisión, el tercero como su coordenada Y y el cuarto está asociado al parámetro B que
depende de la potencia de la fuente.
Por último, tras realizar la integración para los puntos del dominio (con la respectiva
función INTERP), se han considerado dos tipos de soluciones. Por un lado se ha consi-
derado una solución en la que solo se tiene en cuenta el efecto de la difracción debida al
obstáculo y por otro se ha considerado otra solución en la que se tienen en cuenta tanto
la difracción como la reflexión. Para ello, en el primero de los casos no se ha incluido el
término asociado a la onda incidente, (solo para la integración de los puntos en el dominio,
sı́ se ha incluido para la integración en el contorno) mientras que en el segundo sı́ se ha
incluido.
Aplicación al cálculo de pantallas acústicas
Este es el último en desarrollarse de los programas para la implementación del proble-

ma acústico 2D, es por ello también el más complejo. Aunque como se indicó al inicio
de esta sección, la estructura básica de funcionamiento que sigue, es la misma que la de
los programas anteriores, solo que incluye otra serie de funciones para realizar algunos
cálculos adicionales que son necesarios para esta aplicación. Se indican a continuación
cuáles son esas nuevas funciones y qué modificaciones han sufrido las que ya estaban.
Funciones MORSE y DELANYBAZLEY
Estas funciones están directamente relacionadas entre sı́, pues dentro de la primera se lla-
ma a la segunda. El objetivo conjunto es el de calcular el valor de la admitancia(asociada a
la variable BETA del programa) en el contorno. Esta admitancia es un valor real obtenido a
partir del coeficiente de absorción estadı́stico, para el cálculo del cual se utiliza la función
MORSE. Este cofieciente se obtiene a su vez a partir de la admitancia compleja del modelo
empı́rico de Delany y Bazley (asociado a la correspondiente función DELANYBAZLEY).
Básicamente lo que se ha hecho es programar las fórmulas (4.12) a (4.15) desarrolladas
en el apartado anterior de este mismo capı́tulo.
Funciones LOCIN y EXTIN
El uso de estas funciones parte de la idea de separar la integración de los nodos de colo-
cación coincidentes con los nodos de los elementos, de la integración del resto de nodos
de colocación, para esquivar la singularidad que se produce en el primero de los casos(que
para problemas acústicos es menos pronunciada que para otros tipos de problema como
los de potencial). Por ello, la función GHMAT en la que aparecen ambos tipos de integra-
ción, (en la INTERP no aparece el tipo señalado como particular, ya que todos los puntos
que se integran son interiores al dominio) llama a su vez a dos funciones: LOCIN, que es
la que realiza una integración para el caso particular de nodos de colocación coincidentes
con nodos de los elementos y EXTIN, que es la que realiza la integración general para
el resto de casos. En los programas anteriores desarrollados para el problema interior, ex-
terior y de scatter, todas las integraciones se realizan mediante el procedimiento general
que se utiliza en EXTIN, por lo que la función GHMAT incorpora directamente este tipo
de integración y no llama a ninguna función. Tal como se ha indicado antes, la función
INTERP no realiza el caso particular de integración, por lo que en este programa llama
solo a la función EXTIN y en el resto, incopora directamente el procedimiento general al
igual que lo hace GHMAT.
Funciones PREKAW y FLUKAW
La última diferencia que presenta este programa con respecto a los anteriores es la modi-
ficación de la solución fundamental. Esta presenta dos variaciones: por un lado incorpora
un término asociado a la integración en los puntos de la imagen(además de en los de la
fuente) y por otro presenta otro término adicional más (φ ) para los casos en que βs 6= 0.
Para el cálculo de este último término es para lo que se emplean las funciones PREKAW
y FLUKAW. La primera función calcula el valor de φ , que se incluye dentro de la matriz
de coeficientes G, mientras que la segunda, calcula el valor de su derivada respecto de la
normal, que se incluye dentro de la matriz H.
El resto de funciones de este programa, son análogas a las desarrolladas en los an-
teriores(salvo adaptaciones al propio contexto). Tan solo en el programa principal se ha
incluido un bucle tras la llamada a la función de lectura de datos del fichero de entrada,
para poder hacer un barrido de frecuencias, recogiendo los resultados para cada una de
ellas.
Por último, se recogen una serie de diagramas de flujo con los que se pretende resumir
y clarificar las ideas expuestas en esta sección. Se ha decidido incorporarlos todos juntos,
empezando por el más sencillo, para que pueden apreciarse de un simple vistazo las mo-
dificaciones que sufren los programas al ir aumentando la complejidad de los problemas
estudiados. Se ha decidido además incluir los diagramas de los problemas interior y exte-
rior, para elementos constantes, mientras que para el programa de pantallas acústicas, se
ha incluido el de elementos cuadráticos. Con ello puede observarse que todos ellos inclu-
yen las funciones básicas, diferenciándose estas tan solo en la nomenclatura utilizada para
distinguir los distintos tipos de elementos.
Figura 4.5: Diagrama de flujo del problema interior para elementos constantes.
Figura 4.6: Diagrama de flujo del problema exterior para elementos constantes.
Figura 4.7: Diagrama de flujo del problema de pantallas para elementos cuadráticos.
Figura 4.8: Desglose de la función GHMATPQ.

Figura 4.9: Desglose de la función INTERPQ.

Capı́tulo 5
Formulación e implementación del

MEC en problemas acústicos 3D
5.1. Introducción y descripción general de la metdologı́a
En este apartado se desarrollarán los procedimientos necesarios para realizar el plan-

teamiento numérico de la identidad integral (3.44) en problemas tridimensionales. Tendre-
mos como referencia para este capı́tulo las tesis de Davidsson [1] (2004) y Holmstrom [3]
(2001).
Las estrategias utilizadas para el problema 3D son completamente análogas a las utili-
zadas en el 2D, por lo que se omitirán algunos de los desarrollos. Al igual que ya se hizo
para dicho caso 2D, ahora también se aproximará el comportamiento de las variables, de
forma localizada en una serie de elementos en que se divide el contorno del problema.
Aunque para el problema 3D, nos centraremos en el uso de elementos constantes y linea-
les, no empleando por tanto elementos cuadráticos. Tal y como se ha indicado, la identidad
integral es la misma que para el caso 2D, tan solo cambia la solución fundamental, que
para el caso 3D será la segunda expresión de (3.14).
53
5.2. Ecuación integral de Helmholtz para problemas tri-

dimensionales
Para el planteamiento numérico partiremos de la identidad integral (3.44). Aunque co-
mo ya se dijo, dicha ecuación hace referencia al caso particular que se aplicará al cálculo
de pantallas, donde el contorno a discretizar se reduce a Γ p . Mostraremos ahora el caso
general en el que partamos de todo el contorno Γ, ya que como se ha dicho además, no
aplicaremos los códigos 3D al cálculo de pantallas.
Para los problemas interior y exterior, las consideraciones son las mismas que para el
caso 2D. Por lo tanto, podemos decir que para un problema definido en un dominio finito
3D, Ω la ecuación integral en el contorno es la siguiente:
∗
!
∂p
Z
∗ ∗
c0i pi = p0 − + ikβ p p pdΓ (5.1)
Γ ∂n
En esta expresión, la normal está apuntando hacia fuera del dominio interior y el coe-
ficiente c0i vale 1 en el interior del dominio, 1/2 en los puntos i que se encuentren en
contornos suaves y 0, fuera del dominio. La principal diferencia con el problema 2D se
encuentra en que para el resto de casos(para un caso general), en lugar de utilizar la ex-
presión (3.25), se utiliza la siguiente:
θ
ci = (5.2)
4π
siendo θ el ángulo interno de la esquina en radianes.
Ω
y
Figura 5.1: Angulo interno en problemas 3D.
Para el problema exterior, es todo análogo al caso 2D, es decir, para un problema
exterior definido en un dominio (Ω) infinito 3D la ecuación integral en el contorno es la
siguiente:
∗
!
∂p
Z
∗ ∗
ci pi = p0 − + ikβ p p p + pI dΓ (5.3)
Γ ∂n
donde pI representa una onda incidente en el caso de haberla, y la normal apunta en

sentido contrario al caso anterior, es decir, hacia fuera del dominio exterior. El valor del
coeficiente ci es el complementario sobre 1 del de c0i , de forma que:
ci = 1 − c0i (5.4)
5.3. Discretización del contorno
Para resolver la ecuación integral de Helmholtz, en primer lugar habrá que discretizar
el contorno en un número determinado de elementos curvilı́neos. Para el problema 3D
emplearemos elementos constantes y lineales. Para ambos casos la geometrı́a de cada
elemento se representa interpolando entre los nodos, es decir:
n
x = ∑ xi Ni (ξ , η )
i=1
n
y = ∑ yi Ni (ξ , η ) (5.5)
i=1
n
z = ∑ zi Ni (ξ , η )
i=1
donde xi , yi y z son las coordenadas de los nodos, Ni (ξ , η ) son las funciones de forma
definidas en un elemento con coordenadas locales −1 ≤ ξ ≤ 1; −1 ≤ η ≤ 1 y n es el
número de nodos del elemento.
5.3.1. Elementos constantes

Para el caso de elementos constantes, se necesita un único nodo para describir las
variables fı́sicas en el elemento. Dicho elemento puede tener una forma cuadrangular ar-
bitraria y el nodo se sitúa en el centro del elemento. Se usa una única función de forma
que será constante en todo el elemento:
φ1 = 1 (5.6)
nodo
φ =1
1
Figura 5.2: Función de forma constante del problema 3D.

5.3.2. Elementos lineales
Para el caso de elementos lineales, se emplearán elementos cuadrangulares de 4 nodos

y la presión y el flujo en dichos elementos se obtendrán por interpolación, utilizando las
siguientes funciones de forma:
φ1 = 41 (ξ − 1)(η − 1)
φ2 = − 14 (ξ + 1)(η − 1)
(5.7)
φ3 = 41 (ξ + 1)(η + 1)
φ4 = − 14 (ξ − 1)(η + 1)
en el sistema de coordenadas (ξ , η ).
4:(−1,1)
1
η 3:(1,1)
ξ
1:(−1,−1)
2:(1,−1)
Figura 5.3: Ejemplo de función de forma lineal 3D(φ3 ).
5.3.3. Desarrollo del plantemiento numérico de la identidad integral
Al igual que en el caso 2D, independientemente del tipo de elementos seleccionado

para discretizar, (y en consecuencia del número de nodos de los mismos) la aplicación de
la expresión de la identidad integral en un nodo i genérico de la discretización se realiza
de la misma forma, de modo que se introduce la aproximación de la variable de acuerdo
con el tipo de elementos elegido. Se puede escribir por tanto:
NE N
∗
ci pi = p0 − ∑ ∑ (hikj + ikβ j gikj )pkj (5.8)
j=1 k=1
donde NE representa el número de elementos en que se ha dividido el contorno y β j

la admitancia del contorno en cada elemento.
Para el caso de emplear elementos constantes, la ecuación anterior puede reescribirse

de la siguiente forma, teniendo en cuenta que los contornos serán siempre suaves (ci =
1/2), usando n como el producto de NxNE y dejando de utilizar el término asociado a la
fuente de presión interna para simplicar las expresiones siguientes:
n n
1
pi − ∑ H̄i j p j = iρ0 ω ∑ Gi j q j (5.9)
2 j=1 j=1
Introduciendo:
1
Hi j = δi j − H̄i j (5.10)
2
donde δi j es la delta de Kronecker. Introduciendo ahora (5.10) en (5.9), obtenemos:
n n
∑ Hi j p j = iρ0ω ∑ Gi j q j (5.11)
j=1 j=1
donde G y H son 2 matrices nxn. Desarrollaremos a continuación Gi j y Hi j . En un

sistema de coordenadas cartesiano, se obtiene:
√ 2 2 2
e−ik (x j −xi ) +(y j −yi ) +(z j −zi )
Z
Gi j = p dΓ (5.12)
Γ j 4π (x j − xi )2 + (y j − yi )2 + (z j − zi )2
Z
H̄i j = (∇g)T n̂dΓ (5.13)
Γj
siendo:
 
nx
 
n̂ = 
 ny 
 (5.14)
nz
y
 √ !
2 2 2
xe−ik x +y +z 1
− ik + p
 4π (x2 + y2 + z2 ) x2 + y2 + z2 

 √ !
 ye−ik x2 +y2 +z2
 
1 
∇g = 
 − 4π (x2 + y2 + z2 ) ik + p  (5.15)
 2 2
x +y +z  2 
 √ !
2 2 2
 ye−ik x +y +z
 
1
− ik + p

4π (x2 + y2 + z2 ) x2 + y2 + z2
Hi j es igual a H̄i j cuando i 6= j. Cuando i = j los términos de H̄ii son cero, ya que la
normal n̂ y el elemento son perpendiculares entre sı́.
Por tanto H̄ii = 0 → Hii = 21 .
Una vez conocidas la presión y el flujo en los nodos del contorno, se pueden obtener los
correspondientes valores en puntos interiores del dominio, aplicando la siguiente ecuación
algebraica:
n n
pi = ∑ H̄i j p j + iρ0ω ∑ Gi j q j (5.16)
j=1 j=1
donde i vuelve a ser un punto genérico del dominio, en el que se quiere calcular la
presión.
Para el caso de emplear elementos lineales, el planteamiento es completamente análo-

go, pero haciendo uso en este caso de las funciones de forma(ahora no se empleará una
única función de forma con valor unidad, que se simplifica en los cálculos) para efectuar
la interpolación. No se incluye el desarrollo matemático en este documento, ya que no
aporta nada conceptual aparte de lo ya citado y en cambio las expresiones se complican
enormemente. En cualquier caso, como se ha dicho, el objetivo sigue siendo el mismo, es
decir, crear un sistema de ecuaciones de la forma:
H p = iρω Gq (5.17)
Una vez obtenidos los valores de la presión en el contorno, puede calcularse en puntos
interiores del dominio, usando la expresión (5.16).
5.4. Utilización del Método de los Chief Points y solución

por mı́nimos cuadrados
De igual forma a como se empleó en el problema 2D, utilizaremos el Método de los
Chief points para estabilizar la solución del problema exterior 3D. En este caso también
colocaremos chief points en puntos del interior del hueco que queda fuera del dominio.
Compararemos las soluciones obtenidas al emplear dicho método y al no emplearlo.
5.5. Funciones implementadas

Al igual que ocurrı́a para los programas del problema 2D, se han implementado distin-
tos tipos. En este caso, nos hemos centrado en 2 distintos: problema interior y problema
exterior. Ambos se han programado mediante elementos constantes y lineales.
También se utiliza la misma idea de que cada programa se compone a su vez de una
serie de funciones que comparten una estructura común basada en una serie de funciones
básicas que permiten realizar las operaciones necesarias para obtener los resultados tanto
en el contorno como en el dominio. Aunque el problema exterior incorpora algunas modi-
ficaciones respecto al interior. La estructura básica de todos los programas es la siguiente:
fichero de entrada de datos, programa principal que gestiona al resto de funciones; fun-
ción de lectura de los datos del fichero de entrada; función que gestiona la generación de
las matrices H, G y C para la formación del sistema de ecuaciones con los valores del
contorno; función que realiza las integraciones necesarias, a partir de los puntos que se
le suministren, calculando los términos de las matrices H y G(usada tanto para los nodos
del contorno como para los puntos del dominio), función que calcula la matriz C a partir
de la geometrı́a de la discretización, función que reordena las matrices H y G en función
de las condiciones de contorno y resuelve el sistema de ecuaciones que proporciona los
resultados en el contorno y por último una función que gestiona la generación las matrices
H y G para los puntos del dominio y calcula posteriormente los valores de presión y flujo
en el mismo.
Puede observarse que la estructura es muy parecida a la que habı́a en el problema 2D,
salvo que la complejidad de añadir una dimensión más obliga a organizar los datos de una
forma más elaborada. Por ello, para el problema 3D lo que se hará será explicar algunos
detalles de cada una de las funciones que aparecen nuevas con respecto al 2D. Comenza-
remos con el problema interior:
Problema interior 3D
Fichero de entrada de datos

La idea es la misma que para el caso 2D, salvo que la complejidad y dimensión del fi-
chero para el caso 3D es mucho mayor. Por ello se ha incluido una función adicional que
permite la generación automática del fichero de entrada según una serie de parámetros
que están adaptados a los problemas que vamos a resolver en los que la geometrı́a es un
prisma rectangular de dimensiones 6x1x1. La función auxiliar Generación fichero, tiene 6
argumentos de entrada: el primero está asociado a la densidad de la malla, concretamente
al número de nodos por unidad de longitud que se desea que tenga la misma; el segundo
está asociado a la generación de las condiciones de contorno, de manera que un valor 1
generará las correspondientes al problema interior 1 que se explicará en el capı́tulo 6 y un
valor 2 hará lo mismo con el problema interior 2; el tercero es el número de puntos inter-
nos que se quieren calcular; el cuarto es el que decide si utilizar chief points(valor 1) o no
hacerlo(valor 0), el quinto es el número de ellos en caso de emplearlos y el último permite
introducir la geometrı́a aprovechando la simetrı́a con respecto a alguno de los ejes(X,Y o
Z).
En cuanto a la introducción de los nodos, estos también se incluirán para cada elemen-
to en sentido antihorario para el problema interior, de forma que la normal apunte siempre
hacia el exterior del dominio.
Función principal
Esta función también llama al resto de funciones en el orden adecuado y dibuja los resul-
tados en el contorno y en el dominio. Salvo por el hecho de llamar a una función adicional
que ensambla los valores en el contorno, no añade nada nuevo respecto al caso 2D. Los
nombres de las funciones son ACCONBE3D para elementos constantes y ACLINBE3D
para elementos lineales.
Función que genera H, G y C para el contorno
Esta función lo que hace es gestionar la formación de las matrices G, H y C. Ahora es nece-
sario calcular esta última por separado debido a la mayor complejidad que puede presentar
la geometrı́a en problemas tridimensionales. El mayor interés de esta función (nos referi-
mos tanto a GHMATPC3D como a GHMATPL3D) es el de ensamblar adecuadamente los
valores obtenidos de la integración para cada nodo de colocación en las matrices globales
H y G y sobre todo el de ensamblar el valor de c asociado a cada nodo, en la matriz global
C, de forma que no haya que interactuar con el programa durante la ejecución del mismo
para decirle cómo es la discretización de la geometrı́a de nuestro problema.
Función de integración
Esta función es claramente distinta para elementos constantes y para elementos lineales,
por lo que trataremos ambas por separado:
EXTINPC3D es la función para elementos constantes, utiliza elementos cuadrangu-

lares de un único nodo i que se ve influenciado por un elemento j(a través de su también
único nodo), como puede verse en la figura (5.4)
(xj3,yi3,zj3
i3 j3 i3
)
(xj4,yi4,zj4
i4 j4 i4
)
i (xj2,yi2,zj2
i2 j2 i2
)
y (xi1,yi1,zi1)
(xj3,yj3,zj3)
(xj4,yj4,zj4)
x
z
j
(xj1,yj1,zj1)
(xj2,yj2,zj2)
Figura 5.4: Elemento j que actúa sobre nodo i.
Por tanto las variables de entrada para esta función son las coordenadas de los nodos
del elemento j y las del elemento i que contienen al nodo i(además de las propiedades del
medio y el valor de n asociado a la normal, que valdrá 1 en caso de introducir los nodos
del elemento en el sentido recomendado, que es el antihorario para el problema interior, y
-1 en caso contrario).
EXTINPL3D es la función para elementos lineales, utiliza elementos cuadrangulares

de cuatro nodos. Cada uno de estos nodos de un elemento genérico j influencia a su vez
a todos los nodos,(incluido él mismo) representados por un nodo genérico i, como puede
verse en la figura (5.5).
(xi,yi,zi)
y
(xj3,yj3,zj3)
(xj4,yj4,zj4)
x
z
(xj1,yj1,zj1)
(xj2,yj2,zj2)
Figura 5.5: Elemento j que actúa sobre nodo i.
Las variables de entrada para esta función son las coordenadas de los nodos del ele-
mento j y las del nodo i sobre el que actúa el elemento anterior(además de las propiedades
del medio y el valor de n asociado a la normal).
Función para el cálculo de la matriz C

Para el problema 3D esta matriz no es tan fácil de generar pues pueden aparecer tanto
puntos cóncavos y convexos como puntos de silla y es necesario crear una función solo
para su obtención. Dicha función recibe como argumentos de entrada, las coordenadas del
nodo para el que se va a calcular el término libre y las coordenadas de los nodos de los
elementos a los que pertenece dicho nodo. La misión de la función GHMAT3D, dentro de
la que se ubica CMAT3D, será la de proporcionarle a esta última las coordenadas de los
otros nodos que sean necesarias para calcular el término libre en cada uno de los nodos de
la discretización.
Función de ensamblaje y resolución del sistema de ecuaciones

Esta función recibe las matrices H y G generadas mediante GHMAT3D, ası́ como 3 ma-
trices que contienen los 3 tipos de condiciones de contorno posibles (presión conocida,
flujo conocido e impedancia conocida). Cada una de estas 3 matrices(de dimensiones nx2)
contiene el número de los nodos(según el orden de introducción de los mismos) a los que
se va a aplicar la respectiva condición de contorno, en una columna, y el valor de dicha
condición de contorno asociado a cada nodo, en otra columna. Teniendo en cuenta las
condiciones de contorno de los nodos se reordenan las matrices H y G, obteniéndose el
sistema de ecuaciones definitivo que tras resolverse proporciona los resultados de presión
y flujo en cada nodo del contorno.
Función que genera H y G para el dominio

La misión de esta función es análoga a la que realiza GHMAT3D, pero para el dominio
en lugar de para el contorno. Concretamente lo que se hace es gestionar la formación de
las matrices G y H (no es necesario generar C para el dominio ya que todos sus términos
valen 1). Para ello se ensamblan adecuadamente los valores obtenidos de la integración
para cada punto de colocación del interior del dominio en las matrices globales H y G.
Los nombres de las funciones son INTERPC para elementos constantes e INTERPL para
elementos lineales.
Problema exterior 3D
Las dos principales diferencias que presenta el programa del problema exterior respec-
to al del interior, son las mismas que en el caso 2D:
La primera es que al introducir los nodos de cada elemento, se hará recorriendo el

mismo en sentido horario, para que la normal apunte hacia el exterior del dominio.
La segunda será la utilización de una función con la que se implementa el método de

los chief points. Dicha función lleva el mismo nombre CHIEFPOINTS y se empleará para
la estabilización de las soluciones. Se podrá comprobar en el capı́tulo dedicado a los ejem-
plos numéricos, que la estabilización para el caso 3D es más compleja que para el 2D, si
bien la filosofı́a de la función es exactamente la misma en ambos casos, por lo que no se
realizará ningún comentario adicional.
Se incluirá también por último en esta sección un diagrama de flujo para el problema
interior y otro para el exterior, empleando elementos lineales en ambos programas:
Figura 5.6: Diagrama de flujo para problema interior 3D con elementos lineales.
Figura 5.7: Diagrama de flujo para problema exterior 3D con elementos lineales.
Capı́tulo 6
Ejemplos numéricos
6.1. Introducción
En este capı́tulo se expondrán los resultados obtenidos de aplicar los distintos códigos
que se han explicado anteriormente. Lo hemos dividido en tres partes. La primera, está de-
dicada a los ejemplos de validación de los distintos códigos programados, que será la más
extensa, puesto que este trabajo está orientado a abarcar una gran cantidad de formula-
ciones del problema acústico, todas ellas con sus respectivos códigos los cuales se debe
comprobar que funcionan adecuadamente. La segunda está dedicada a una aplicación con-
creta del problema acústico 2D y de uso frecuente en ingenierı́a civil, como es el cálculo
de pantallas de atenuación acústica. Finalmente se dedicará la última parte a comparar y
analizar los resultados obtenidos.
En todos los casos se ha seleccionado el aire como el fluido a través del cual se produce
la propagación de las ondas. Se han incluido por tanto en los códigos las prociedades
fı́sicas de este fluido: ρ = 1,21g/l y c = 343m/s
6.2. Ejemplos de validación
6.2.1. Problema 1 interior 2D

En esta subsección se incluirán los resultados del problema correspondiente a una
cavidad rectangular de dimensiones 6x1, en la que la pared de la izquierda tiene como
condición de contorno conocida, la presión(valor 1) y en el resto de paredes es conocido
67
Capı́tulo 6. Ejemplos numéricos 68
el flujo(valor 0).
Para este problema se han incluido los resultados de aplicar elementos constantes,
lineales y cuadráticos para la integración numérica, ya que todos ellos proporcionan re-
sultados algo diferentes, aunque como se verá, no tanto como en el caso del siguiente
problema. Para todos ellos se presentará el valor de la presión en los nodos del contorno,
ası́ como el valor de la presión a lo largo de una recta del interior del dominio, paralela a
las paredes exteriores. El resultado es el mismo para cualquiera de dichas lı́neas paralelas,
independientemente de a qué distancia se encuentre de los paredes, ya que el valor de la
presión solo depende de la coordenada x(paralela a las anteriores) y es independiente de
la coordenada y(perpendicular a la anterior).
q=0
1
y[m]
p=1 q=0
0
q=0
−1
−2
−2 0 2 4 6 8
x[m]
Figura 6.1: Geometrı́a de la cavidad 2D(condiciones de contorno del problema 1).
En primer lugar se presentan los resultados obtenidos mediante elementos constantes.

Como puede observarse, los elementos constantes proporcionan una muy buena aproxima-
ción. La amplitud de la onda que proporcionan es ligeramente superior a la de la solución
exacta, pero una cantidad casi imperceptible.
Término libre C(p) y presión en el contorno ps
0.5
0.4
C(p)
0.3
0.2
MEC
0.1 Solución analítica analítica
0
10 20 30 40 50 60 70
Nodos
1
Presión (Pa)
−1 MEC
Solución analítica
−2
0 10 20 30 40 50 60 70
Nodos
Figura 6.2: Solución en los nodos del contorno con elementos constantes.
Presión en el dominio de la cavidad rectangular

1.5
MEC
0.5
Presión [Pa]
−0.5
−1
−1.5
0 1 2 3 4 5 6
x [m]
Figura 6.3: Solución a lo largo de una recta del dominio con elementos constantes.
A continuación se presentan los resultados obtenidos mediante elementos lineales. Al

contrario de lo que cabrı́a esperar, los elementos lineales proporcionan una solución más
inexacta que los constantes. Aunque como se puede observar, la aproximación no deja de
ser buena.
0.5
0.4
C(p)
0.3
0.2
MEC
0
10 20 30 40 50 60 70
Nodos
1
Presión (Pa)
−1 MEC
−2
0 10 20 30 40 50 60 70 80
Nodos
Figura 6.4: Solución en los nodos del contorno con elementos lineales.

1.5
MEC
0.5
Presión [Pa]
−0.5
−1
−1.5
0 1 2 3 4 5 6
x [m]
Figura 6.5: Solución a lo largo de una recta del dominio con elementos lineales.
Por último se presentan los resultados obtenidos mediante elementos cuadráticos, que
tal como se puede apreciar, proporcionan la solución más parecida a la exacta, siendo casi
imperceptible la diferencia con respecto a la misma(más aún que con los elementos cons-
tantes que ya vimos que era muy buena).
0.5
0.4
C(p)
0.3
0.2
MEC
0
10 20 30 40 50 60 70
Nodos
1
Presión (Pa)
−1 MEC
−2
0 10 20 30 40 50 60 70 80
Nodos
Figura 6.6: Solución en los nodos del contorno con elementos cuadráticos.

1.5
MEC
0.5
Presión [Pa]
−0.5
−1
−1.5
0 1 2 3 4 5 6
x [m]
Figura 6.7: Solución a lo largo de una recta del dominio con elementos cuadráticos.
6.2.2. Problema 2 interior 2D

En esta subsección se incluirán los resultados del problema correspondiente a la cavi-
dad rectangular anterior, en la que ahora la pared de la izquierda y la de la derecha tienen
como condición de contorno conocida, la presión(valor 1) y en la pared de arriba y en la
de abajo es conocido el flujo(valor 0).
Para este problema también se han incluido los resultados de aplicar elementos cons-
tantes, lineales y cuadráticos para la integración numérica, ya que todos ellos proporcionan
resultados considerablemente diferentes.
En primer lugar se presentan los resultados obtenidos mediante elementos constantes.

Como puede observarse, los elementos constantes proporcionan una buena aproximación,
aunque la amplitud de la onda es algo menor a la de la solución exacta.
0.5
0.4
C(p)
0.3
0.2
MEC
0
10 20 30 40 50 60 70
Nodos
2
Presión (Pa)
−2 MEC
−4
0 10 20 30 40 50 60 70
Nodos
Figura 6.8: Solución en los nodos del contorno con elementos constantes.

4
MEC
3
1
Presión [Pa]
−1
−2
−3
−4
0 1 2 3 4 5 6
x [m]
Figura 6.9: Solución a lo largo de una recta del dominio con elementos constantes.
A continuación se presentan los resultados obtenidos mediante elementos lineales. Es-

tos vuelven a proporcionar una solución menos exacta que los elementos constantes. Aun-
que al contrario que ocurre con estos últimos, la amplitud de onda proporcionada por los
elementos lineales es mayor que la de la solución exacta.
0.5
0.4
C(p)
0.3
0.2
MEC
0
10 20 30 40 50 60 70
Nodos
2
Presión (Pa)
−2 MEC
−4
0 10 20 30 40 50 60 70 80
Nodos
Figura 6.10: Solución en los nodos del contorno con elementos lineales.

4
MEC
3
1
Presión [Pa]
−1
−2
−3
−4
0 1 2 3 4 5 6
x [m]
Figura 6.11: Solución a lo largo de una recta del dominio con elementos lineales.
Por último se presentan los resultados obtenidos mediante elementos cuadráticos, que
como era de esperar, siguen proporcionando la solución más parecida a la exacta.
0.5
0.4
C(p)
0.3
0.2
MEC
0.1
Solución analítica analítica
0
10 20 30 40 50 60 70
Nodos
2
Presión (Pa)
−2 MEC
−4
0 10 20 30 40 50 60 70 80
Nodos
Figura 6.12: Solución en los nodos del contorno con elementos cuadráticos.

4
MEC
3
1
Presión [Pa]
−1
−2
−3
−4
0 1 2 3 4 5 6
x [m]
Figura 6.13: Solución a lo largo de una recta del dominio con elementos cuadráticos.
6.2.3. Problema exterior 2D

En esta subsección se incluirán los resultados del problema correspondiente a una cir-
cunferencia pulsante de radio unidad, en la que el flujo es conocido en todo el contorno de
la misma, siendo igual a uno.
Para este problema se ha incluı́do una única gráfica de los valores de presión calcu-
lados en el contorno y otra única gráfica para los valores calculados en el dominio, para
cada una de las dos frecuencias estudiadas. Esto es debido a que para este problema, no
solo los elementos cuadráticos, sino también los constantes y los lineales proporcionan un
resultado prácticamente exacto en todos los casos. Mostraremos los resultados empleando
por ejemplo elementos lineales. Como se ha dicho, se han estudiado dos frecuencias dis-
tintas, una que no produce resonancia y otra que sı́ la produce, para poder comprobar cómo
en el primero de los casos la solución no necesita ser estabilizada, mientras que en el se-
gundo caso, es necesario recurrir al método de los chief points para lograr la estabilización.
En primer lugar se presentan los resultados obtenidos cuando la frecuencia de pulsa-

ción es f rec = c/(2π ) siendo c la velocidad de propagaciñ en el aire(Esta frecuencia se ha
obtenido de igualar a 1 el número de onda k).
Para este caso, no se produce resonancia, por lo tanto la solución es estable y el resul-
tado obtenido es el mismo tanto si se aplica como si no, el método de los chief points. Se
presenta a continuación dicho resultado:
Término libre C(p) y presión en el contorno

0.5
C(p) MEC
0
5 10 15 20 25 30 35 40
Nodos
Presión(p. imaginaria) Presión (p. real)
200 MEC
0
5 10 15 20 25 30 35 40
Nodos
100
MEC
0
5 10 15 20 25 30 35 40
Nodos
Figura 6.14: Solución en los nodos del contorno, para la primera frecuencia con elementos lineales.
Presión en el dominio
Sol. analítica(parte real)
MEC(parte real)
300
Sol. analítica(parte imaginaria)
Parte real e imaginaria de la presión (Pa)
MEC(parte imaginaria)
200
100
−100
−200
−300
2 4 6 8 10 12 14 16 18 20
Distancia a la superficie (m)
Figura 6.15: Solución a lo largo de una recta del dominio infinito, para la primera frecuencia con elementos
lineales.
A continuación se muestran los resultados para una frecuencia de pulsación f rec =

131Hz que sı́ produce resonancia. Se mostrarán por tanto los resultados para la solución
desestabilizada y para la estabilizada.
En primer lugar se muestran los resultados de la solución desestabilizada ante el

fenómeno de resonancia:
0.5
C(p)
MEC
0
5 10 15 20 25 30 35 40
Nodos
Presión(p. imaginaria) Presión (p. real)
500
MEC
0
5 10 15 20 25 30 35 40
Nodos
100
50 MEC
0
5 10 15 20 25 30 35 40
Nodos
Figura 6.16: Solución sin usar chief points (desestabilizada) en los nodos del contorno, para la segunda
frecuencia con elementos lineales.
400
MEC(parte real)
300 MEC(parte imaginaria)
200
100
−100
−200
−300
−400
2 4 6 8 10 12 14 16 18 20
Figura 6.17: Solución sin usar chief points (desestabilizada) a lo largo de una recta del dominio infinito,
para la segunda frecuencia con elementos lineales.
Ahora se mostrará la solución estabilizada mediante el método de los chief points. Para
conseguirlo, se ha colocado un único punto en el interior del hueco de la circunferencia,
que queda fuera del dominio.

0.5
C(p)
MEC
0
5 10 15 20 25 30 35 40
Nodos
Presión (p. real)
400
200 MEC
0
5 10 15 20 25 30 35 40
Presión(p. imaginaria)
Nodos
50
MEC
0
5 10 15 20 25 30 35 40
Nodos
Figura 6.18: Solución estabilizada en el contorno, (empleo del Método de los Chief Points ) para la
segunda frecuencia con elementos lineales.
400
MEC(parte real)
300 MEC(parte imaginaria)
200
100
−100
−200
−300
−400
2 4 6 8 10 12 14 16 18 20
Figura 6.19: Solución estabilizada en el dominio, (empleo del Método de los Chief Points) para la segunda
frecuencia con elementos lineales.
6.2.4. Problema de scatter 2D

En esta subsección se incluirán los resultados de un problema exterior, en el que la cir-
cunferencia de radio unidad, en lugar de ser pulsante, va a ser tratada como un obstáculo
ante la llegada de los distintos frentes de onda. Se recogerán los resultados de dos frentes
distintos: uno correspondiente a una onda plana incidente y el otro correspondiente a una
onda esférica.
Los distintos tipos de elementos, al igual que ocurre con los diferentes problemas ex-
teriores, proporcionan resultados casi idénticos, por lo que se presentarán tan solo los
resultados utilizando uno de los tipos(por ejemplo elementos cuadráticos, que darán resul-
tados ligeramente más exactos).
En primer lugar se muestran las gráficas correspondientes a una onda plana que sigue
una trayectoria horizontal, con sentido de izquierda a derecha. Se le ha asignado amplitud
unidad a la onda.
Término libre y presión en el contorno

0.5
C(p)
MEC
Sol. analítica
0
50 100 150 200 250 300 350 400
Nodos
Presión(p.real)
0
MEC
−2
0 50 100 150 200 250 300 350 400
Nodos
Presión(p.imaginaria)
0
MEC
−2
0 50 100 150 200 250 300 350 400
Nodos
Figura 6.20: Solución en el contorno cuando hay una onda plana incidente.
Se muestran además los resultados en el dominio teniendo en cuenta solo el fenómeno

de difracción de la onda y teniendo en cuenta tanto la difracción como la reflexión.
Solución considerando solo difracción

5
4
1
3
2 0.5
1
y[m]
0 0
−1
−0.5
−2
−3
−1
−4
−5
−5 −4 −3 −2 −1 0 1 2 3 4 5
x[m]
Figura 6.21: Solución en el dominio considerando solo la difracción, cuando hay una onda plana incidente.
Solución considerando reflexión y difracción

5
1.5
4
3 1
2
0.5
1
y[m]
0
0
−1
−2 −0.5
−3
−1
−4
−5
−5 −4 −3 −2 −1 0 1 2 3 4 5
x[m]
Figura 6.22: Solución en el dominio considerando reflexión y difracción, cuando hay una onda plana
incidente.
A continuación se muestran los resultados correspondientes a una onda esférica cuya

fuente de emisión se ha situado en el punto de coordenadas (-4,0) y a la que se le ha
asignado un valor unidad a su parámetro de potencia de la fuente.
Término libre y presión en el contorno

C(p) 0.5
MEC
Sol. analítica
0
50 100 150 200 250 300 350 400
Nodos
Presión(p.real)
0
MEC
−1
0 50 100 150 200 250 300 350 400
Nodos
Presión(p.imaginaria)
0
MEC
−1
0 50 100 150 200 250 300 350 400
Nodos
Figura 6.23: Solución en el contorno cuando hay una onda esférica incidente.
Se muestran también los resultados en el dominio teniendo en cuenta solo el fenómeno

de difracción de la onda y teniendo en cuenta tanto la difracción como la reflexión.
Solución considerando solo difracción

5
4 0.3
3
0.2
2
0.1
1
y[m]
0 0
−1
−0.1
−2
−0.2
−3
−4 −0.3
−5
−5 −4 −3 −2 −1 0 1 2 3 4 5
x[m]
Figura 6.24: Solución en el dominio considerando solo la difracción, cuando hay una onda esférica inci-
dente.
Solución considerando reflexión y difracción

5 3
4 2.5
3
2
2
1.5
1
1
y[m] 0
0.5
−1
0
−2
−0.5
−3
−1
−4
−1.5
−5
−5 −4 −3 −2 −1 0 1 2 3 4 5
x[m]
Figura 6.25: Solución en el dominio considerando reflexión y difracción, cuando hay una onda esférica
incidente.
6.2.5. Problema interior 3D

En esta subsección se incluyen los resultados de los problemas de dos cavidades con
forma de prisma rectangular de dimensiones 6x1x1. Se han recogido ambos en una única
sección, ya que tanto las ideas como los resultados son los mismos que para el caso 2D.
Esto quiere decir que las condiciones de contorno son las mismas para cada uno de los
problemas, salvo que ahora en lugar de imponerse en lados de un rectángulo, se imponen
en caras de un prisma rectangular. Además la presión sigue variando a lo largo de una
única coordenada y,(esta es la misma que definimos como x para el problema 2D) que
para el caso 3D es la que une la cara con condición de contorno impuesta en presión con
la cara enfrentada, para el problema 1 y ambas caras con condición de contorno impuesta
en presión para el problema 2.
Mostramos a continuación un esquema de la geometrı́a de la cavidad para el problema

3D:
1 6
5
4
z[m]
0.5
3
2
0 1
0 0.5 y[m]
1 0
x[m]
Figura 6.26: Geometrı́a de la cavidad 3D.
En este caso se emplean elementos constantes y lineales para el cálculo de los campos
de presiones. Dichos campos a lo largo de cualquier recta que recorra la definida como
coordenada y para cualquier pareja de valores x, z constantes(recta paralela al eje y) son
los mismos que los mostrados para los dos problemas interiores 2D, en las figuras respec-
tivas de las subsecciones 1 y 2 de este apartado.
Se ha representado para ambos problemas el perfil de presiones en el dominio para 9

rectas que recorren la coordenada y, asociadas a 9 parejas de valores de las coordenadas
x, z, constantes.
Se muestra en primer lugar una gráfica con el valor del término libre calculado en cada
nodo. Puede observarse que pasa de haber dos grupos de valores distintos como habı́a en
el caso 2D(asociados a nodos del interior de tramos rectos y a nodos de esquinas) a haber
3 grupos de valores, asociados a nodos del interior de las caras, nodos de las aristas y
nodos de las esquinas. En dicha gráfica aparece también el valor calculado de la presión
en los nodos del contorno para el primer problema interior. Solo se representa la solución
obtenida mediante el MEC y no la solución analı́tica, ya que la ordenación de los nodos no
siguen en todos los casos la dirección de la coordenada y. En cualquier caso, dichos valores
de presión en el contorno son correctos, ya que los valores en el dominio lo son, tal como
se verá después. Se han recogido los resultados empleando solo elementos lineales para
no incluir gráficas repetidas, ya que tal como se ha comentado, los resultados obtenidos
para cada recta paralela al eje y son los mismos que para el caso 2D.

0.5
0.4
C(p)
0.3
0.2
MEC
0.1
0 50 100 150 200 250 300 350 400 450
Nodos
1
Presión (Pa)
−1
MEC
−2
0 50 100 150 200 250 300 350 400 450
Nodos
Figura 6.27: Solución en los nodos del contorno del prisma rectangular con elementos lineales.
Pueden observarse en la gráfica anterior los 3 valores distintos que toma el término
libre en función de la situación de cada nodo en la cavidad: caras, aristas y vértices. Se
pueden distinguir fácilmente los 8 valores menores(0.125), asociados a los 8 vértices de
la cavidad. En segundo lugar se muestran para el primer problema interior, los perfiles
de presión calculados a lo largo de 9 rectas del dominio paralelas al eje y. Presentaremos
tanto una vista 2D, correspondiente a una sola de las rectas, como una vista 3D incluyendo
todas las rectas.
Presión en el dominio de la cavidad
MEC
2
1.5
0.5
Presión [Pa]
−0.5
−1
−1.5
−2
0 1 2 3 4 5 6
y [m]
Figura 6.28: Solución a lo largo de una recta del dominio, paralelela al eje y con elementos lineales.
Vista 3D de la presión en el dominio de la cavidad

1
0.8
0.5
Presión [Pa]
0.6
0
0.4
−0.5
0.2
−1
0
5 −0.2
4 −0.4
3 −0.6
2
−0.8
MEC
y [m] 1
Solución analítica 0.5 −1
Cavidad rectangular
x [m]
Figura 6.29: Vista 3D de la solución a lo largo de rectas del dominio, paralelelas al eje y con elementos
lineales.
A continuación se muestran los resultados análogos a los anteriores, pero para el segun-
do problema interior. Puede observarse que los valores del término libre son los mismos
para todos los nodos, puesto que la geometrı́a del prisma es la misma. Se presenta también
el valor de la presión en los nodos del contorno de este segundo problema.

0.5
0.4
C(p)
0.3
0.2
MEC
0.1
0 50 100 150 200 250 300 350 400 450
Nodos
2
Presión (Pa)
−2
MEC
−4
0 50 100 150 200 250 300 350 400 450
Nodos
Figura 6.30: Solución en los nodos del contorno del prisma rectangular con elementos lineales.
Como últimas gráficas de esta subsección, se incluyen las correspondientes a los per-
files de presión a lo largo de las 9 rectas del dominio, paralelas al eje y, para el segundo
problema interior, (vista 3D) ası́ como una vista 2D asociada a una sola de las rectas.
4
MEC
3
1
Presión [Pa]
−1
−2
−3
−4
0 1 2 3 4 5 6
x [m]
Figura 6.31: Solución a lo largo de una recta del dominio, paralelela al eje y con elementos lineales.
Vista 3D de la presión en el dominio de la cavidad

4 3
3
2
2
1 1
Presión [Pa]
0
0
−1
−1
−2
−3 −2
−4
0 −3
2 1
MEC 0.5
Solución analítica 4
0
y [m] 6
Cavidad rectangular
x [m]
Figura 6.32: Vista 3D de la solución a lo largo de rectas del dominio, paralelelas al eje y con elementos
lineales.
6.2.6. Problema exterior 3D

En esta subsección se incluirán los resultados del problema correspondiente a una es-
fera pulsante de radio uno, en la que el flujo es conocido en todo el contorno con valor
unidad para todo el mismo.
Para este problema, debido a que los elementos constantes y los lineales producen re-
sultados casi idénticos, se ha decidido incluir las gráficas para uno solo de los tipos(se ha
escogido los elementos lineales, aunque prácticamente no se distinguen los resultados con
respecto a los que proporcionan los elementos constantes).
Al igual que en el problema exterior 2D, se han estudiado dos frecuencias distintas,
una que no produce resonancia y otra que sı́ la produce, para poder comprobar también la
eficacia del Método de los Chief Points para estabilizar la solución.
En primer lugar se presentan los resultados obtenidos cuando la frecuencia de pulsa-

ción es f rec = c/(2π )Hz siendo c la velocidad de propagación en el aire.
Para este caso, no se produce resonancia, por lo tanto la solución es estable y el resul-
tado obtenido es el mismo tanto si se aplica como si no, el Método de los Chief Points. Se
presenta a continuación dicho resultado:

MEC
C(p)
0.5
0.4
0 5 10 15 20
Nodos
MEC
Presión(p.real)
1.5
1 Solución analítica
0.5
0
−0.5
0 5 10 15 20
Nodos
Presión(p.imag.)
1.5 MEC
0.5
0
−0.5
0 5 10 15 20
Nodos
Figura 6.33: Solución en los nodos del contorno de la circunferencia, para la primera frecuencia.
0.5
MEC(parte real)
0.4
0.3
P. real e imag. de la presión

0.2
0.1
−0.1
−0.2
−0.3
1 2 3 4 5 6 7 8 9 10
x[m]
Figura 6.34: Solución a lo largo de una recta del dominio infinito, para la primera frecuencia.
Se muestra también una representación 3D de la misma solución anterior, en la que se

puede ver la esfera pulsante junto al campo de presiones que genera en el dominio, hasta
una distancia 10 veces mayor que su radio(la esfera tiene radio unidad):
Vista 3D de la presión presión en el dominio

0
10
−1 8
6
0.5 4
0 2
−0.5
0 x[m] Sol. analítica(parte real)
y[m] MEC(parte real)
Figura 6.35: Vista 3D de la solución a lo largo de una recta del dominio infinito, para la primera frecuencia.
A continuación se muestran los resultados para una frecuencia de pulsación f rec =

6c/(2π )Hz que sı́ produce resonancia. Se mostrarán por tanto los resultados para la solu-
ción desestabilizada y para la estabilizada. Para esta frecuencia, colocaremos chief points,
pero la solución no es independiente del número de los mismos que utilicemos. Veremos
cómo por el hecho de utilizar alguno(que sea válido y produzca restricción) la solución se
estabiliza mucho. Pero al aumentar el número, la solución en el dominio(luego veremos

por qué nos referimos especı́ficamente al mismo) se estabiliza cada vez más, hasta llegar
al número óptimo de chief points. A partir de dicho número, la solución comenzará de
nuevo a desestabilizarse hasta llegar a ser incluso menos exacta que la que obtenida sin
emplear el método.
En primer lugar se muestran los resultados de la solución desestabilizada ante el

fenómeno de resonancia:
C(p)
0.5
MEC
0.4
0 5 10 15 20
Nodos
Presión(p.real)
1.5
1
0.5
0 MEC
−0.5 Solución analítica
0 5 10 15 20
Nodos
Presión(p.imag.)
1.5 MEC
0.5
0
−0.5
0 5 10 15 20
Nodos
Figura 6.36: Solución en el contorno, para la segunda frecuencia sin emplear chief points.
1
MEC(parte real)
0.5 MEC(parte imaginaria)
−0.5
−1
−1.5
1 2 3 4 5 6 7 8 9 10
x[m]
Figura 6.37: Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia sin emplear
chief points.
Mostramos también la representación 3D de la solución anterior.
1
10
8
0
6
−1 4
2
0.5
0 Sol. analítica(parte real)
−0.5 0 x[m]
MEC(parte real)
y[m]
Figura 6.38: Vista 3D de la solución a lo largo de una recta del dominio infinito, para la segunda frecuencia
sin emplear chief points.
Ahora mostramos la solución obtenida al colocar un chief point, que como se anti-
cipó estabiliza mucho la solución tanto en el dominio como en el contorno:

C(p)
0.5
MEC
0.4 Solución analítica
0 5 10 15 20
Nodos
Presión(p.real)
1.5
1
0.5
0 MEC
−0.5
0 5 10 15 20
Nodos
Presión(p.imag.)
1.5 MEC
0.5
0
−0.5
0 5 10 15 20
Nodos
Figura 6.39: Solución en el contorno, para la segunda frecuencia empleando 1 chief point.
1
0.8 MEC(parte real)
0.6

0.4
0.2
−0.2
−0.4
−0.6
−0.8
1 2 3 4 5 6 7 8 9 10
x[m]
Figura 6.40: Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia empleando
1 chief point.
Se presenta también la representación 3D:

1
10
0 8
6
−1
4
0.5 2
0 Sol. analítica(parte real)
−0.5 0 x[m]
MEC(parte real)
y[m]
Figura 6.41: Vista 3D de la solución a lo largo de una recta del dominio infinito, para la segunda frecuencia
empleando 1 chief point.
Tras esto, mostramos la solución obtenida al colocar 10 chief points, que es el número
óptimo de los mismos para estabilizar la solución en el dominio. Aunque cabe destacar que
la solución en el contorno es menos exacta que en el caso de emplear un solo chief point.
Esto es debido a que al aumentar el número de restricciones, la solución por mı́nimos
cuadrado del sistema de ecuaciones se ve alterada. En cualquier caso, la estabilización de
la solución es mayor que la alteración que sufre debido a este efecto. Puede apreciarse en
las siguientes gráficas:
C(p)
0.5
MEC
0.4
0 5 10 15 20
Nodos
Presión(p.real)
1.5
1
0.5
0 MEC
−0.5
0 5 10 15 20
Nodos
Presión(p.imag.)
1.5 MEC
0.5
0
−0.5
0 5 10 15 20
Nodos
Figura 6.42: Solución en el contorno, para la segunda frecuencia empleando 10 chief points.
1
0.8 MEC(parte real)
0.6
0.4
0.2
−0.2
−0.4
−0.6
−0.8
1 2 3 4 5 6 7 8 9 10
x[m]
Figura 6.43: Solución a lo largo de una recta del dominio infinito, para la segunda frecuencia empleando
10 chief points.
No incluimos la representación 3D para este caso, ya que es muy parecida a la que se

mostró anteriormente para el caso de colocar un solo chief point (6.41).
Por último, mostramos la solución obtenida al colocar 20 chief points; un número

considerablemente mayor que el óptimo, que era 10. Puede apreciarse que tanto la solución
en el contorno como en el dominio, vuelven a estar desestabilizadas, siendo el resultado
proporcionado más inexacto que el que tenı́amos antes de aplicar el método de los chief
points. Esto se explica teniendo en cuenta que el número de ecuaciones de restricción
es tan alto, que altera enormemente la solución por mı́nimos cuadrados del sistema de
ecuaciones. En este caso, dicho efecto es mucho mayor que el efecto estabilizador, por
ello es por lo que podemos apreciar en las siguientes gráficas cómo se desestabiliza la
solución tanto en el contorno como en el dominio:
C(p)
0.5
MEC
0.4
0 5 10 15 20
Nodos
Presión(p.real)
1.5
1
0.5
0 MEC
−0.5
0 5 10 15 20
Nodos
Presión(p.imag.)
1.5 MEC
0.5
0
−0.5
0 5 10 15 20
Nodos
1
0.8 MEC(parte real)
0.6 MEC(parte imaginaria)
0.4
0.2
−0.2
−0.4
−0.6
−0.8
−1
1 2 3 4 5 6 7 8 9 10
x[m]
Tampoco incluimos la representación 3D para este caso, ya que se distingue muy poco
de la proporcionada para el caso en que no se colocaban chief points (6.38).
6.3. Aplicación al cálculo de pantallas de atenuación acústi-

ca
A continuación aplicaremos el código 2D desarrollado al cálculo de pantallas de ate-
nuación acústica y mostraremos los resultados obtenidos para los distintos casos estudia-
dos.
6.3.1. Pantalla con contorno reflejante

En el primer caso se ha colocado una pantalla de dimensiones 5m x 5m en un semi-
espacio reflejante. Las superficies de la pantalla son también reflejantes. Además hay una
fuente puntual situada 2.5m por encima del contorno del semiespacio y a 5m en horizontal
del contorno de la pantalla. El objetivo es calcular los resultados para 3 receptores situados
5m a la derecha(en horizontal) de la pantalla y a unas alturas sobre el semiespacio de 0,2.5
y 5m respectivamente. De esta forma mostraremos el efecto de atenuación que realiza la
pantalla a través del coeficiente de pérdida por inserción (IL).
6 Receptores
β=0
3
4
fuente
2
2
βs=0
0 1
0 5 10 15
Figura 6.46: Geometrı́a de la primera pantalla(reflejante).

Se muestran a continuación los resultados obtenidos para las distintas frecuencias y

para cada uno de los receptores:
Espectro de pérdida por inserción

10
Receptor 1
Receptor 2
Receptor 3
0
pérdida por inserción (IL)(dB)
−10
−20
−30
−40
1 2 3 4
10 10 10 10
frecuencia(Hz)
Figura 6.47: Espectro de pérdida por inserción de la primera pantalla, obtenido mediante los códigos
adjuntados en los anexos.
Este resultado, puede también validarse al compararlo con el obtenido en el libro de

Maeso y Aznárez [4] (2004) para el mismo problema. Los valores son prácticamente
idénticos.
Figura 6.48: Espectro de pérdida por inserción de la primera pantalla, sacado del libro de Maeso y Aznárez
[4] (2004).
En ambas gráficas, puede observarse a pesar de las oscilaciones asociadas a las posicio-
nes de la fuente y los receptores que la pérdida por inserción es creciente lı́neas generales
al disminuir la cota del receptor.
6.3.2. Pantalla con contorno absorbente

El segundo caso consiste en colocar una pantalla triangular con propiedades absorben-
tes. El semiespacio sigue siendo reflejante, pero la pantalla en cambio se va a estudiar
para 3 recubrimientos absorbentes diferentes: σ1 =300000 Nsm−4 , σ2 =100000 Nsm−4 y
σ1 =10000 Nsm−4 . En los tres casos el espesor del tratamiento será de 0.2m. La forma de
la pantalla será la de un triángulo is’sceles de 5m de base y 3m de altura. La fuente se co-
locará a la altura del semiespacio y 5m a la izquierda de la pantalla y el receptor, también
está a la altura del semiespacio, pero 5m a la derecha de la pantalla. Al igual que en el
caso anterior, el objetivo es comprobar la eficacia de la pantalla a través del coeficiente de
pérdida por inserción.
2
β 0
fuente βs=0 receptor
0
−2 0 2 4 6 8 10 12 14 16 18
Figura 6.49: Geometrı́a de la segunda pantalla(absorbente).

Presentamos a continuación los resultados obtenidos para los tres tratamientos y para
las distintas frecuencias, incluyendo además el caso en que la pantalla sea completamente
reflejante.
Espectro de pérdida por inserción

10
σ1=300000 Nsm−4
σ2=100000 Nsm−4
0 σ1=10000 Nsm−4
pérdida por inserción (IL)(dB)
reflejante
−10
−20
−30
−40
1 2 3 4
10 10 10 10
frecuencia(Hz)
Figura 6.50: Espectro de pérdida por inserción de la segunda pantalla, obtenido mediante los códigos
adjuntados en los anexos.
Al comparar con los resultados sacados del libro de Maeso y Aznárez [4] (2004), se
aprecia también que los resultados son prácticamente idénticos. Presentamos por último
la gráfica sacada de dicho libro.
Figura 6.51: Espectro de pérdida por inserción de la segunda pantalla, sacado del libro de Maeso y Aznárez
[4] (2004).
Puede observarse que para cada frecuencia, al reducir el valor de σ , aumenta la pérdida
por inserción (IL). Ello se justifica teniendo en cuenta que disminuir σ implica aumentar
el coeficiente de absorción del contorno.
Otra conclusión un poco más profunda es que al aumentar la frecuencia, se reduce el

efecto de la difracción y debido a ello, el efecto de la reflexión supera al de la absorción.
Ası́ se explica que tan solo para el tratamiento σ3 se consiga un aumento significativo de la
eficacia de la pantalla(con respecto a la situación totalmente reflejante) para las frecuencias
medias.
6.4. Análisis de resultados

En este apartado, haremos un breve resumen de los resultados obtenidos a lo largo
del capı́tulo, no basándonos en los valores numéricos, ya que la mayorı́a de los ejemplos
son de validación, sino en las ideas utilizadas, que son extrapolables a cualquier problema
estudiado mediante el MEC.
En primer lugar cabe destacar la eficacia de los distintos tipos de elementos utilizados,
pues aunque solo los elementos cuadráticos proporcionan un resultado casi exacto, los
constantes y los lineales proporcionan también un muy buen resultado en general. Tan so-
lo para el segundo problema interior, (tanto para el problema 2D como para el 3D) son los
elementos cuadráticos los únicos que dan un muy buen resultado. En cambio para todos
los problemas exteriores, los resultados proporcionados por los tres métodos son práctica-
mente idénticos.
Se ha comprobado también la eficacia del método de los chief points para estabilizar
las soluciones de los problemas exteriores para aquellas frecuencias de pulsación que pro-
ducen resonancia. También ha quedado de manifiesto la diferente complejidad que hay
entre utilizar este método para problemas 2D y utilizarlo para problemas 3D. Para geo-
metrı́as análogas entre ambos espacios dimensionales , como son una circunferencia en el
caso 2D y una esfera en el caso 3D, en el primero de los casos, basta con colocar un único
chief point en el interior de la circunferencia, para estabilizar la solución. En cambio para
el caso 3D, se producen distintos grados de estabilización, alcanzándose el resultado más
parecido a la solución exacta, para un número óptimo de chief points a partir del cual la
solución vuelve a desestabilizarse llegando a obtenerse peores resultados que en el caso

de no emplear dicho método de estabilización.
Capı́tulo 7
Resumen y conclusiones
7.1. Resumen de los contenidos desarrollados

A lo largo de este trabajo se han recogido los principales planteamientos teóricos del
problema acústico aplicado al MEC, tanto para el caso 2D como para el 3D. Dichos plan-
teamientos se han programado a su vez mediante códigos en lenguage de Matlab y se han
acompañado en todos los casos de ejemplos numéricos de validación para comprobar que
los códigos funcionan correctamente. Además para el caso 2D se ha anádido un ejemplo
de aplicación práctica en el mundo de la Ingenierı́a Civil, como es el cálculo de pantallas
acústicas.
7.2. Conclusiones del trabajo

El Método de los Elementos de Contorno, es una herramienta muy útil cuando se trata
de modelizar problemas donde hay una razón superficie/volumen baja, ya que se consu-
men menos recursos computacionales que en otros método como por ejemplo los métodos
de discretización de volumen(Método de los Elementos Finitos, Método de Diferencias Fi-
nitas, Método de Volúmenes Finitos). También lo es para problemas en los que el dominio
es infinito, ya que si la solución fundamental cumple la condición de radiación de Som-
merfeld, es posible discretizar todo el dominio, lo cual no es posible en dichos métodos
de discretización de volumen. En tales condiciones se ha utilizado el MEC en este trabajo
para modelizar problemas de potencial y acústicos, pero también es ampliamente utili-
zado en disciplinas como el electromagnetismo,la mecánica de fluidos y la mecánica de
100
Capı́tulo 7. Resumen y conclusiones 101
la fractura. Concretamente, este trabajo está fundamentalmente enfocado hacia problemas

acústicos, en los cuales el MEC es muy superior al MEF en todos los aspectos comentados.
Aunque en cualquier caso, generalmente el principal incoveniente de este tipo de for-

mulaciones es que dan como resultado matrices llenas, lo que implica que los recursos de
almacenamiento y tiempo computacional tienden a crecer proporcionalmente al cuadrado
del tamaño del problema. Esto se pudo comprobar al ejecutar el código correspondiente
a los problemas interiores 3D, que involucran una mayor cantidad de nodos y el tiempo
de computación asociado no es despreciable. Para una densidad de la malla de 3 nodos
por unidad de longitud, el tiempo de ejecución en el ordenador utilizado era de unos 30
segundos. Al aumentar a 4 nodos por unidad de longitud, el tiempo se elevaba a 2 minutos
y para una densidad de la malla de 6 nodos por unidad de longitud, el tiempo de compu-
tación alcanzaba los 5 minutos.
En cambio, las matrices de elementos finitos son tı́picamente bandeadas (los elemen-
tos son solo localmente conexos) y las necesidades de almacenamiento para este tipo de
matrices crece tı́picamente linealmente con el tamaño del problema. Para aliviar este in-
coveniente se pueden usar técnicas de compresión (como expansiones multipolares o al-
macenamiento jerárquico de matrices), aunque el costo es la complejidad que se añade y
la utilidad está muy ligada al tipo de problemas que se resuelve y a su geometrı́a. En los
códigos desarrollados se propusieron funciones adicionales para realizar una integración
independiente en aquellos puntos de colocación que coincidieran con nodos dentro del
elemento integrado, de forma que se esquivara la singularidad. En cualquier caso tal como
recoge la bibliografı́a y se ha podido comprobar, dicha singularidad es menos pronunciada
para el problema acústico que para el de potencial y es posible para el primero de ellos
realizar una integración conjunta para todos los puntos de colocación sin que esto práctica-
mente afecte al tiempo de computación para mallas no demasiado densas. De hecho, entre
los códigos de este documento para el problema acústico 2D, se incluye un programa en
el que se realiza la integración conjunta para todos los puntos de colocación(Anexos A
y B, para la validación de los ejemplos numéricos del problema interior y exterior res-
pectivamente) y otro en el que se realiza por separado la integración de los puntos con la
singularidad(Anexo E, para el cálculo de pantallas acústicas).
También es importante destacar que para la mayorı́a de los problemas, no solo los ele-
mentos cuadráticos, sino también los constantes y los lineales, producen resultados muy
precisos, distinguiéndose muy poco entre ellos y con respecto a la solución exacta(analı́ti-
ca, para los problemas de validación). Si bien es cierto que para algunos problemas, como
el segundo problema interior estudiado en este documento, tan solo los elementos cuadráti-
cos proporcionan resultados muy similares a la solución analı́tica conocida. Aunque la
mayor eficacia que ofrecen este tipo de elementos, implica también un salto cualitativo en
la complejidad que implica su programación e implementación.
7.3. Trabajo futuro

Si bien este trabajo recoge una gran cantidad de códigos de aplicación del MEC, to-
davı́a quedan muchos otros que se pueden implementar para lograr una biblioteca más
completa. Ası́, podrı́a por ejemplo considerarse la posibilidad de usar elementos cuadráti-
cos de 8 nodos para el problema 3D. Tal como se ha indicado en la sección anterior, la
implementación de estos elementos implica una mayor complejidad, pero también propor-
ciona unos resultados notablemente mejores en algunos casos. Es por ello por lo que se
considera de interés incorporarlos. Dando un paso más, podrı́a plantearse la posibilidad de
desarrollar los siguientes planteamientos del MEC:
✐ Extensión a Elementos de Contorno isogeométricos.
✐ Aplicación de técnicas de resolución jerárquicas.
✐ Elementos de Contorno espectrales.
Además, el campo de la modelización mediante el uso de métodos matemáticos avan-

zados, tal como se indicó en el capı́tulo introductorio, ha desarrollado otras herramientas
además del MEC para el estudio de problemas, como es por ejemplo el MEF. Podrı́a con-
siderarse su estudio e implementación e incluso hacerlo para la interacción entre ambos
métodos, conocida mayoritariamente como BEM/FEM, debido a que estas son las siglas
de ambos métodos en inglés.
Por último, como trabajo futuro, se pretende recoger los principales códigos y resul-
tados de este trabajo, ası́ como los que se pretenden realizar, citados en los párrafos an-
teriores de este apartado y junto con otros que vayan surgiendo durante la realización de
dicho trabajo, en un documento escrito en inglés, con el objetivo de conseguir una mayor
difusión de los conocimientos incorporados.

TFG Pedro Palomo

Cargado por

Copyright:

Formatos disponibles

TFG Pedro Palomo

Cargado por

Información del documento

Descripción original:

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

TFG Pedro Palomo

Cargado por

Copyright:

Formatos disponibles

Trabajo Fin de Grado

Grado en Ingeniería Civil

Programación en Matlab del Método de los

Autor: Pedro Palomo Cantador

Tutor: Luis Rodríguez de Tembleque Solano

Equation Chapter 1 Section 1

Departamento de Mecánica de los Medios

Sevilla, Junio de 2015

Programación en Matlab del Método de los

Autor: Pedro Palomo Cantador

Tutor: Luis Rodríguez de Tembleque Solano

Departamento de Mecánica de los Medios Continuos y

Sevilla, Junio de 2015

Autor: Pedro Palomo Cantador

Tutor: Luis Rodríguez de Tembleque Solano

Acuerdan otorgarle la calificación de:

El Secretario del Tribunal

Quisiera dedicar esta sección a agradecer a mi tutor Luis Rodrı́guez de Tembleque

Agradecer a mis compañeros de carrera su ayuda e interés por mi persona, no solo en

2. Conceptos fundamentales de acústica 14

3. Representación integral para problemas escalares 19

3.2. Solución fundamental . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4. Formulación e implementación del MEC en problemas acústicos 2D 35

5. Formulación e implementación del MEC en problemas acústicos 3D 53

6.2.2. Problema 2 interior 2D . . . . . . . . . . . . . . . . . . . . . . . 72

7. Resumen y conclusiones 100

A. Interacción LaTeX-Matlab para la incorporación de códigos a documentos 105

B. Código acústico interior 2D 108

C. Código acústico exterior 2D 128

D. Código para problemas de scatter 2D 134

E. Código para pantallas acústicas 2D 144

F. Código acústico interior 3D 166

G. Código acústico exterior 3D 209

H. Aplicación del MEC a problemas de potencial 211

I. Código para problemas de potencial 2D 215

3.1. Ecuación Integral y Solución Fundamental. . . . . . . . . . . . . . . . . . . . . 21

4.1. Tipologı́as de elementos empleadas para discretizar. . . . . . . . . . . . . . . . . 37

5.1. Angulo interno en problemas 3D. . . . . . . . . . . . . . . . . . . . . . . . . . 55

6.1. Geometrı́a de la cavidad 2D(condiciones de contorno del problema 1). . . . . . . . . 68

6.26. Geometrı́a de la cavidad 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

H.1. Geometrı́a del problema de potencial estudiado. . . . . . . . . . . . . . . . . . . 212

H.1. Resultados en el dominio con elementos constantes. . . . . . . . . . . . . 213

1.1. Estudio y motivación

1.3. Estructuración del trabajo

Conceptos fundamentales de acústica

2.2. Ecuaciones básicas del problema acústico

2.2.1. La ecuación de onda linealizada

La primera es la ecuación de conservación de la masa, la segunda es la ecuación de

2.2.2. Ondas armónicas planas: la ecuación de Helmholtz

2.2.3. Condiciones de contorno

p = p̄ : presión conocida(condición tipo Dirichlet)

2.3. Magnitudes para medir el sonido

2.3.1. Intensidad acústica

2.3.2. Impedancia acústica especı́fica

2.3.3. Nivel de presión sonora

El rango de presión detectable por el oı́do humano es bastante amplio,(entre 20µ Pa y

donde hp2 i se refiere a la media cuadrática temporal.

2.3.4. Coeficiente de pérdida por inserción

< p2despues > < p2antes >