Nuevas formulaciones numéricas para el estudio de propagación de ondas en fluidos

Tesis para la obtención del grado Doctor por la Universidad de Sevilla Nuevas formulaciones numéricas para el estudio de propagación de ondas en fluidos Autor: Antonio Cerrato Casado Director: José Ángel González Pérez Profesor Titular Codirector: Luis Rodríguez-Tembleque Solano Profesor Contratado Doctor Departamento de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2016 Tesis Doctoral: Nuevas formulaciones numéricas para el estudio de propagación de ondas en fluidos Autor: Antonio Cerrato Casado Director: José Ángel González Pérez Codirector: Luis Rodríguez-Tembleque Solano El tribunal nombrado para juzgar la Tesis arriba indicada, compuesto por los siguientes doctores: Presidente: Vocales: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: A mis padres, ellos se merecen esta línea más que nadie. Agradecimientos En primer lugar, deseo expresar mi agradecimiento al director de esta tesis doctoral, José Ángel González, por su tiempo, dedicación, apoyo y cercanía. Igualmente agradecer al codirector de la misma, Luis Rodríguez de Tembleque, su disponibilidad en todo momento, consejos y, cómo no, también el apoyo prestado. Esta andadura la comencé gracias a Ramón Abascal, al cual agradezco a título póstumo la ayuda en los primeros pasos del camino. He de expresar mi más sincero agradecimiento al Profesor M.H. Ferri Aliabadi por haberme acogido en el Department of Aeronautics of Imperial College (London). Sin duda, una experiencia que siempre recordaré con mucho cariño. Quisiera expresar mi gratitud a Andrés Sáez, Director del Departamento, por su disponibilidad, accesibilidad y por haberme permitido hacer uso del departamento en todo momento. Este agradecimiento es extensible a todos aquellos con los que he compartido lugar de trabajo durante estos años. A Antonio y a Gabriel, con los que he compartido muchas horas en el despacho. A los profesores Héctor, Pedro, Mario, Antonio y Fernando. También a Reyes, a M ª José, a Manolo y al resto de personas con las que he coincidido en el departamento. Gracias a mis amigos, por todo su apoyo y cariño. Incluido Hamdy, el cual, por suerte, apareció un día para arreglar un entuerto que hice. Pero sobre todo, quisiera agradecer a mis padres, Antonio y M ª Manuela, su inmenso apoyo incondicional, al que he recurrido en innumerables ocasiones y sin el cual esta tesis nunca se habría escrito. Por supuesto también quiero dar las gracias a mis hermanos, Manuel y Ana, pues siempre he podido contar con ellos. Antonio Cerrato Casado Sevilla, 2016 III Resumen En esta tesis se estudian diferentes problemas de propagación de ondas lineales en fluidos desde el punto de vista numérico. La tesis está compuesta por tres artículos donde se abordan dos temas fundamentalmente: la propagación de oleaje sobre batimetrías variables aplicando el Método de los Elementos de Contorno (MEC), el Método de los Elementos Finitos (MEF) y el Método de los Elementos Espectrales (MEE). También se tratan problemas de interacción fluido-estructura en acústica mediante formulaciones particionadas no-simétricas MEC-MEF. Para el tratamiento del problema de propagación de oleaje mediante el MEC, se presenta en esta tesis una solución fundamental que permite resolver el problema para profundidad variable en una dirección. Esta morfología de batimetría es especialmente relevante, pues está presente en un gran número de problemas reales en ingeniería. Además, se propone un acoplamiento MEC-MEF, incorporando una discretización en ambos métodos basada en elementos espectrales. Mediante dicha discretización, problemas con batimetrías complejas y con condiciones de contorno de dominio abierto pueden ser resueltos de forma precisa y eficiente, evitando reflexiones indeseadas en el contorno exterior y reduciendo extremadamente el error de polución en la solución. Debido a diversos motivos, en el problema de interacción fluido-estructura, las mallas no tienen por qué coincidir en la interfase. Además, en ocasiones resulta de interés resolver ambos dominios mediante métodos diferentes y por computación en paralelo. El problema acústico de interacción fluido-estructura es tratado en esta tesis mediante un acoplamiento MEC-MEF; donde el MEC es utilizado para modelar el fluido y el MEF para la estructura. Para el acoplamiento se estudian y comparan dos estrategias diferentes que permiten una resolución particionada del problema: una basada en el método Mortar y otra en los Multiplicadores de Lagrange Localizados. Esta última es por primera vez utilizada este tipo de problemas. V XII Índice 2.3.5. Propagación de ondas armónicas y dispersión 19 2.3.6. Derivación de la ecuación de pendiente suave 21 2.3.7. Cambio de variable: ecuación de Helmholtz 26 3. Aproximaciones numéricas 27 3.1. Método de los Elementos Finitos 27 3.1.1. Discretización del dominio 28 3.2. Método de los Elementos de Contorno 29 3.2.1. Introducción 29 3.2.2. Aplicación del Método de los Elementos de Contorno (MEC) a la ecuación de Helmholtz bidimensional 30 Ecuación integral de contorno 30 Discretización de la ecuación integral de contorno 31 3.2.3. Solución fundamental de la ecuación de Helmholtz 32 3.2.4. Solución fundamental para número de onda variable 32 El problema transformado en el espacio de Fourier 34 Integración en el plano complejo 35 Análisis del comportamiento de la función de Green en el espacio transformado 38 Solución en el espacio real 39 3.2.5. Resumen 40 3.3. Aproximación mediante elementos espectrales 44 3.3.1. Introducción 44 3.3.2. Polinomios de Legendre 46 3.3.3. Cuadratura de Legendre-Gauss-Lobatto (LGL) 47 3.3.4. Aproximación polinomial óptima 49 3.3.5. Aproximación pseudoespectral: Elementos de referencia 50 Elemento espectral 1-D 51 Elemento espectral 2-D 52 3.3.6. Integración numérica en los elementos espectrales 52 3.3.7. Resumen 53 3.4. Métodos de acoplamiento 54 3.4.1. Introducción 54 3.4.2. Interacción Fluido-Estructura. Problema variacional 54 3.4.3. Interacción fluido-estructura mediante el Método Mortar 56 Formulación variacional continua del problema de acoplamiento 56 Formulación variacional discreta 56 Ecuaciones del problema de acoplamiento fluido-estructura 58 3.4.4. Método de los Mulplicadores de Lagrange Localizados (MLL) 59 Formulación variacional contínua 59 Formulación variacional discreta 60 Resolución del problema de interacción fluido-estructura 61 Localización de los nodos del marco de acoplamiento: Regla del momento nulo 62 Índice XIII 3.4.5. Resumen 63 4. Conclusiones y trabajos futuros 65 4.1. Resumen y conclusiones 65 4.2. Trabajos futuros 68 Bibliografía 69 Apéndice A.Artículo I 83 Apéndice B.Artículo II 123 Apéndice C.Artículo III 159 Índice de Figuras 193 Índice alfabético 195 Glosario 195 1 Introducción 1.1 Introducción Esta tesis aborda el problema de la propagación de ondas lineales en fluidos desde el punto de vista numérico. El estudio de la propagación de ondas armónicas es de primordial interés en un gran número de problemas prácticos en campos tan diversos como la elastomecánica, la sismología, el electromagnetismo, técnicas de radar y sonar, espectrometría atómica, la acústica o la propagación de oleaje. Esta tesis se centra especialmente en estas dos últimas disciplinas, donde el fenómeno de propagación se puede modelar de forma linealizada en el dominio de la frecuencia mediante la ecuación de Helmholtz. 1.2 Motivación y revisión histórica 1.2.1 El problema acústico El interés por la creación y la propagación del sonido ha acompañado a la humanidad desde el inicio de su existencia, pues el sonido forma parte del entorno en el que vivimos. Ya en la antigua Grecia, matemáticos como Pitágoras o Aristóteles se preocuparon por la física del problema desde el punto de vista de la concepción musical. En el siglo I a.C. aparece el primer trabajo conocido en cuyo contenido se trata la propagación del sonido desde un punto de vista físico e ingenieril. En él, aparecen recomendaciones a tener en cuenta para el diseño de anfiteatros con el fin de conseguir una buena acústica en el recinto. La obra se titulaba De Architectura libri decem y su autor fue el arquitecto romano Marco Vitruvio Polión. Fue en el siglo XVIII cuando la acústica empezó a desarrollarse como ciencia en sí misma. Joseph Sauveur (1653-1716), físico francés, fue el que acuñó el término Acústica para designar a todo lo referente al sentido del oído. Durante dicho siglo, grandes físicos y matemáticos como Newton (1642-1727), Leibniz (1646-1716) Euler (1707-1783), 1 2Capítulo 1. Introducción d ´ Alembert (1717-1783) o Lagrange (1736-1813), por citar algunos, estudiaron el movimiento oscilatorio en problemas de la mecánica clásica, aportando los fundamentos necesarios para el posterior desarrollo teórico de la acústica. El siglo XIX estuvo marcado por los avances en la comprensión de la electricidad y el electromagnetísmo, destacando Faraday, Maxwell o Hertz; pero también por la teoría de la elasticidad con Navier, Cauchy, Stokes etc.. Ambas áreas de estudio contribuyeron de manera esencial a la comprensión de la física de la propagación del sonido. Este fenómeno fue ampliamente estudiado por Hermann von Helmholtz (1824-1894), de quien, entre otros trabajos, cabe destacar sus teoría sobre resonancia y su obra: On the sensations of tone as a physiological basis for the theory of music. También son de enorme importancia en esta disciplina las aportaciones de Lord Rayleigh (1842-1919), en su obra: Theory of Sound. Para una revisión histórica más detallada sobre la acústica el lector puede consultar [147, 110]. En la actualidad, la Ingeniería Acústica tiene una notable relevancia por su carácter multidisciplinar y sus importantes aplicaciones en muchos ámbitos de la ciencia y la ingeniería. Los diferentes campos de investigación y de interés práctico son muchos: ruido ambiental, ruido industrial y laboral, acústica arquitectónica, tecnología sonar, ultrasonidos, ecografía, acústica submarina, bioacústica, instrumentación acústica, etc. En general, la acústica está presente en infinidad de ámbitos de la actividad humana, lo que explica el interés de su estudio. 1.2.2 Propagación del oleaje El estudio del comportamiento del oleaje es uno de los principales pilares de la ingeniería marítima y la oceanografía. El conocimiento sobre su mecanismo de propagación juega un importante papel como factor de diseño de actuaciones costeras, mantenimiento de puertos y playas y el estudio de transporte de sedimentos. Durante muchos siglos, el fenómeno físico de propagación del oleaje ha interesado a muchos e importantes autores. Entre todos ellos se podrían destacar los trabajos de Airy [ 2 ] y Stokes [ 151 ], que describieron modelos para la propagación de oleaje lineal y no lineal respectivamente, o los de Boussinesq [ 26 ] y Lamb [ 97 ]. En su trabajo, Boussinesq desarrolló una teoría de propagación no lineal en aguas someras [ 26 ], y Lamb (discípulo de Stokes y Maxwell), propuso el modelo de long-wave para propagación de olas en su obra Hydrodynamics . La modelización numérica del problema de propagación puede ser complicada y requerir un gran coste computacional dependiendo de las hipótesis consideradas en la formulación. Sin embargo, son muchos los casos en los que los modelos lineales pueden ser considerados como una buena aproximación del fenómeno real o, en cualquier caso, ser considerados el punto de partida para estudios más complejos. Los modelos lineales pueden proporcionar la obtención de una estimación a priori que nos permita acotar las zonas de estudio que requieran el empleo de modelos no-lineales. En general, los modelos lineales pueden ser válidos por sí solos cuando la altura de ola ( H ) es pequeña en comparación con la longitud de la onda (L), es decir: H/L<< 1. El modelo lineal de propagación de oleaje más extendido por su simplicidad y fidelidad con la física real del problema fue presentado por W. Berkhoff en los años 70 [ 18 , 19 ]. Este modelo, conocido como mild-slope equation o Ecuación de Pendiente Suave (EPS), 1.3 Solución numérica de la ecuación de onda 3 es capaz de reproducir los efectos de asomeramiento, refracción, difracción y reflexión y, como su propio nombre indica, es sólo conveniente para pendientes pequeñas, es decir, cuando ∇h/kh 1 , donde h es la profundidad y k es el número de onda. Históricamente, fue Eckart a principios de los 50 quien propuso un modelo muy parecido a la EPS [ 50 ], sin embargo, el área de aplicación de este modelo estaba limitada a aguas someras. En la década de los 70 aparecen también otras aproximaciones similares a la de Berkhoff, como la de Ito y Tanimoto [ 87 ] que, a diferencia de la EPS, fueron formuladas en el dominio del tiempo. Schönfeld [ 140 ], de manera contemporánea a Berkhoff, también deriva la EPS, pero de forma diferente y menos compacta. Smith y Sprink derivan también la EPS en 1975 en [148], donde hacen una comparación entre esta teoría y la de Eckart. Muchos trabajos han seguido ampliando y analizando la formulación presentada por Berkhoff, ya sea en su forma elíptica o parabólica, como el trabajo publicado por Radder [ 131 ]. Booij, presentó en su tesis, a principios de los ochenta, una formulación basada también en la versión parabólica de la EPS en la que se incluía el efecto de las corrientes [ 23 ]. Más tarde, publicó un trabajo complementario a la EPS donde concluía que el modelo era válido para pendientes en la batimetría de hasta 1:3. Algo después Kirby también presentó una extensión de la EPS en el dominio del tiempo [ 91 ], donde resuelve problemas de propagación sobre fondos ondulados. Otros de sus trabajos relacionados con la forma parabólica de la EPS son [93, 92]. Por otro lado, Massel propuso en [ 112 ] una extensión del modelo de Berkhoff que incorpora términos de orden superior de la forma del fondo, incluyendo los proporcionales a la curvatura, y los modos evanescentes. Chamberlain y Porter, presentaron a mediados de los noventa la modified mild-slope equation [ 35 ]. Esta formulación es probablemente la que más ha prevalecido sobre otras extensiones, en ella se conservan los términos de segundo orden que en un primer momento fueron despreciados por Berkhoff. Otras extensiones posteriores de la EPS, parecidas a la anterior, pueden ser encontradas en [ 152 , 36 ]. Lee et al. et al. [ 100 ], sigiendo los trabajos de Copeland [ 43 ] y Massel [ 112 ], presentaron una versión hiperbólica de la EPS extendida. Esta formulación es esencialmente como la de Copeland, salvo que retiene los términos de orden superior de la variación del fondo. Otros trabajos interesantes que incluyen variaciones sobre la EPS orginal son [103, 79, 80, 78]. 1.3 Solución numérica de la ecuación de onda Determinados fenómenos físicos de propagación de ondas en fluidos, como es el caso de la acústica lineal o la propagación de oleaje de pequeña amplitud, pueden ser reproducidos matemáticamente mediante la ecuación de Helmholtz. La ecuación de Helmholtz puede tratarse de forma analítica o numérica. En el campo analítico encontramos métodos como el de las imágenes [ 119 ], el de separación de variables, expansión en series de autovalores, transformadas de Fourier, etc. [ 53 ]. Sin embargo las soluciones analíticas son por lo general difíciles de hallar y están reservadas a casos específicos con geometrías simples. Los métodos numéricos son, en cambio, mucho más versátiles y permiten resolver gran cantidad geometrías, inhomogeneidades y condiciones de contorno diferentes. 4Capítulo 1. Introducción El desarrollo de los métodos numéricos ha venido de la mano del desarrollo de los ordenadores modernos desde la segunda mitad del siglo pasado. En general, los métodos numéricos convierten el problema diferencial inicial en un problema discreto, obteniendo así un valor estimado de la solución en un número finito de puntos del dominio de estudio. A pesar de no ser exactos llegan a conseguir una aproximación de la solución lo suficientemente buena como para satisfacer las necesidades requeridas para la mayor parte de los problemas físicos o de ingeniería actuales. Dos han sido los principales métodos numéricos usados para la solución de las ecuaciones elípticas o problemas de contorno: el Método de las Diferencias Finitas (MDF) y el Método de los Elementos Finitos (MEF). En ambos métodos los sistemas de ecuaciones lineales resultantes de su aplicación están formados por matrices dispersas. El MEC, aunque ha recibido menos atención en la literatura [ 41 ], también ha sido ampliamente utilizado en este tipo de problemas, pues ofrece ventajas como la reducción de la malla necesaria o el cumplimiento de la condición de radiación en el infinito. Otro importante grupo de métodos de resolución de ecuaciones en derivadas parciales que ha tenido un relevante papel son los Métodos Espectrales (MMEE), de los que hablaremos más adelante. 1.3.1 El Método de las Diferencias Finitas En el MDF el problema diferencial es tratado en su forma fuerte. En este método, las derivadas de las variables son aproximadas mediante operadores diferenciales basados en desarrollos en serie y los valores de la solución son obtenidos en una serie de puntos colocados en forma de rejilla sobre el dominio. El MDF es, en general, fácil de implementar y ofrece una buena convergencia en la solución tanto al refinar la malla como al aumentar el orden de los operadores diferenciales. Además, las matrices que contienen los coeficientes de los sistemas de ecuaciones a los que da lugar son dispersas. Algunos trabajos destacables, en cuanto a la ecuación de Helmholtz se refiere, son los de Harari y Turkel [ 71 ] o Turkel et al. [ 156 ]. Un aspecto negativo del MDF es que, debido a la topología de la malla, el modelo no se puede adaptar a cualquier contorno, dificultando una correcta aproximación geométrica del dominio. 1.3.2 El Método de los Elementos Finitos Originalmente el MEF fue una herramienta desarrollada por ingenieros que se planteó para la resolución de problemas estructurales en la década de los 50 [ 157 ], sin embargo, rápidamente se convirtió en el método general para resolver ecuaciones diferenciales, aplicándose a todos los ámbitos de la ciencia y la ingeniería. El MEF se deriva a partir de la forma variacional del problema, dividiéndose el dominio en un número finito de subregiones donde la solución se interpola mediante un conjunto de funciones de aproximación a partir de los valores nodales. Al contrario que con el MDF, el MEF, se puede adaptar fácilmente a la geometría del dominio en el que se requiere econtrar la solución. Tanto el MEF como el MDF han sido ampliamente usados para resolver el problema acústico o de propagación de oleaje. Sin embargo, ambos comparten algunas deficiencias, como son el alto número de grados de libertad requeridos en comparación con el MEC, el error de polución o la correcta aproximación de las condiciones de contorno, especialmente en dominios abiertos. El error de polución será tratado con más detalle en el Capítulo 3. 1.3 Solución numérica de la ecuación de onda 5 Muchos han sido los trabajos que han abordado la condición de contorno necesaria para modelizar dominios abiertos, donde se ha de cumplir la condición de radiación en el infinito (o condición de radiación de Sommerfeld), desde la década de los 70. Bettess y Zienkiewicz [ 20 ], y posteriormente Lau y Ji [ 99 ], implementaron los elementos infinitos para lidiar con este problema. Berkhoff [ 19 ] y Tsay y Liu [ 154 ], se sirvieron del Método de los Elementos Híbridos que originalmente propuso Mei [ 115 ]. En esta línea continuaron Houston [ 77 ] Chen [ 37 ] y Tsay et al. [ 155 ]. Chen y Mei [ 38 ], propusieron también una solución de la condición de contorno abierta a través de una expansión analítica de autovalores. Posteriormente este trabajo fue continuado por Houston [ 77 ] entre otros. Algo posterior son la condiciones de contorno de tipo Dirichlet to Neumann (DtN) presentadas por Keller y Givoli [ 89 ], esta técnica consiste en un procedimiento analítico que reproduce de forma exacta la condición de radiación de Sommerfeld para algunos caso particulares de geometría sencilla [ 60 , 59 ]. El método fue aplicado en el problema acústico por Harari y Huges [ 70 ] y Harari [ 70 ]. Continuando en esta línea, Bonet y Givoli presentaron la condición de contorno tipo Discrete Non Local (DNL), basadas en la DtN. Otra técnica para simular dominios abiertos, y una de las más relevantes, es el uso de los Perfectly Matched Layers (PMLs), presentados en [ 14 ] para el modelado bidimensional de propagación de ondas electromagnéticas mediante el MEF. Desde entonces esta técnica ha sido ampliamente utilizada para problemas de propagación de ondas en el dominio de la frecuencia en todos los ámbitos; en acústica y oleaje podemos destacar [ 113 , 10 , 11 , 118 ]. Además de todos lo métodos mencionados anteriormente, existen otros menos sofisticados basados en los métodos iterativos. Algunos trabajos en los que estos métodos son utilizados son [13, 12, 150, 39, 108]. 1.3.3 El Método de los Elementos de Contorno Se puede considerar el MEC como una técnica numérica basada en el método de los residuos ponderados [ 41 ], como el MEF, exceptuando que como función de ponderación se elige la función de Green del problema con el fin de eliminar la necesidad de discretizar el interior del dominio [ 28 ]. O también puede ser visto como una implementación numérica de la Ecuación Integral de Contorno, donde el concepto de elemento proveniente del MEF es introducido [96]. Es difícil definir el origen del MEC, y ni mucho menos es un objetivo de esta tesis, pero uno de los trabajos que se pueden destacar es de principios de los 60, cuya autoría pertenece a Jawson [ 88 ]. En este trabajo se utiliza la aproximación más simple posible en el MEC, que son los elementos constantes. Cuatro años después Rizzo extendió la formulación de Jawson a la teoría de la elasticidad [ 134 ], iniciando un periodo fecundo para el MEC. Una aportación que puede considerarse clave es la de Lachat y Watson [ 96 ], que incorporaron los métodos de discretización y las funciones de aproximación que se venían utilizando en el MEF, dando lugar a una versión del MEC más actual. Otros trabajos destacables de los inicios del método son los llevados a cabo por Hess [ 75 ], extendiendo la aplicación a problemas tridimensionales de fluidos, o Banerjee y Daviles [ 7 ], que para problemas de cimentaciones por pilotaje aplicó lo que más adelante se llamaría el método indirecto. Para una extensa revisión histórica del método se refiere al lector a [41]. 6Capítulo 1. Introducción La ecuación de Helmholtz ha sido ampliamente tratada mediante el MEC. Algunos trabajos destacables en el ámbito de la acústica son los de Seznec [ 141 ], De Lacerda et al. [ 45 , 46 ], donde se trata la difracción producida por barreras acústicas mediante el MEC. Li et al. [ 104 ] utilizan el MEC para resolver el problema acústico en el semi-espacio teniendo en cuenta un plano reflejante de impedancia arbitraria. Harari y Hughes presentan en [ 68 ] un estudio comparativo del coste computacional requerido por el MEC y el MEF para el problema acústico. El MEC también ha sido empleado para resolver problemas de Control Activo de Sonido (CAS), como en los trabajos de Bai y Chang [6] o Bracanti et al. [27]. El MEC también ha sido utilizado para resolver problemas de propagación de oleaje, sin embargo, su aplicación es más difícil que en otros campos como la acústica o el electromagnetismo. Esto es debido a que la propagación del oleaje se produce sobre un medio no homogéneo (debido a los cambios de profundiad). Es por esto que, lo modelos basados en el MEC se ven limitados a zonas off-shore o interior de puertos, donde en muchos casos la batimetría puede ser considerada constante. Algunos trabajos destacables en este campo son los de Hauguel [ 72 ] y Shaw y Falby [ 142 ], a finales de los 70, donde el MEC es acoplado al MEF para satisfacer la condición de contorno de dominio abierto. Issacson y Qu [ 86 ] y Hamanaka [ 82 ] propusieron sendos modelos basados en el MEC con el fin de reproducir la propagación del oleaje en puertos, donde la profundidad se puede asumir constante, teniendo en cuenta diferentes coeficientes de reflexión en la línea de costa. Más tarde, Lee y Williams [ 101 ] tratan el mismo problema teniendo en cuenta la variabilidad del espectro de oleaje. A pesar de esta limitación, el problema de propagación del oleaje puede ser tratado mediante el MEC recurriendo a técnicas especiales, como el Método de Elementos de Contorno de Reciprocidad Dual (DRBEM), utilizado por Zhu [ 169 ], Zhu et al. [ 170 ] o Liu [ 107 ]. En este método es necesario recurrir a integrales de dominio con el fin de incluir la inhomogeneidad del medio. Otra posible alternativa es utilizar una solución fundamental que permita modelar mediante el MEC problemas que incluyan batimetría variable sin necesidad de recurrir a ninguna otra técnica, como es el caso de Cerrato et al. [ 32 ], donde se resuelven problemas para batimetría variable en una dirección. La solución fundamental propuesta en dicho trabajo está basada en la función de Green propuesta por Belibassakis [9] para este tipo de problemas. 1.3.4 Métodos Espectrales Los MMEE están basados en una aproximación de la solución de la ecuación diferencial mediante una familia de funciones ortogonales, lo que confiere al método una rápida convergencia comparado con otras técnicas numéricas. Aunque se pueden encontrar trabajos relacionados con el tema ya en los años treinta [ 98 ] los MMEE recibieron una especial atención a finales de los 60 y durante los 70, en paralelo con el desarrollo de las transformadas rápidas. Algunas referencias destacables son las de Fox y Parker [ 95 ] y Gottlieb y Orszag [ 44 ]. Dentro de la familia de los MMEE encontramos también los pseudo-espectrales, que reciben este nombre cuando la interpolación de la solución se realiza a través de valores nodales y no de los coeficientes que acompañan al espectro de funciones ortogonales [145]. 1.3 Solución numérica de la ecuación de onda 7 Los MMEE y pseudo-espectrales constituyen una importante alternativa a los métodos anteriormente descritos a la hora de resolver ecuaciones en derivadas parciales elípticas como la de Helmholtz o la EPS. Sus principales ventajas son su rápida convergencia y la derivabilidad de sus soluciones. En los MMEE, las técnicas más utilizadas son el Método de Colocación, el Método de Galerkin y el Método Tau. El Método de Colocación se deriva directamente a partir de la formulación fuerte del problema, como ocurre en el MDF. Los operadores diferenciales se obtienen derivando las funciones de aproximación nodales. El Método Galerkin Espectral, consiste al igual que en el MEF, en forzar la ortogonalidad del residuo con respecto a todas las funciones de aproximación de la solución. A diferencia del MEF, el espacio no es divido en subregiones y la solución y la función de ponderación se aproximan mediante la familia de polinomios ortogonales elegida y sus respectivos coeficientes. El Método Tau se deriva de la formulación débil del problema y puede ser visto como un caso particular del Método de Petrov-Galerkin. Su uso es adecuado para problemas con condiciones de contorno no periódicas sin embargo, en problemas con coeficientes no constantes no es tan eficiente como el método de colocación. En general es un método poco utilizado. Una explicación detallada sobre estos métodos, su evolución, su funcionamiento e implementación puede ser encontrada en los trabajos de Gottlieb [ 44 ] o Canuto et al. [30, 34, 31] entre otros. La aproximación mediante elementos espectrales fue propuesta por Patera en 1984 [ 127 ] para la solución de las ecuaciones de Navier-Stokes para fluidos incompresibles. El Método de los Elementos Espectrales (MEE), como Patera lo denominó, consiste en una combinación del MEF con los MMEE. Las funciones de aproximación que propuso originalmente Patera estaban asociadas a valores nodales, por lo que se trata de un método pseudo-espectral, sin embargo, existen otras funciones de aproximación de tipo espectral aptas para el MEE como son las modal boundary-adapted basis [ 34 , 143 ]. Este método híbrido cuenta una convergencia exponencial a medida que aumenta el orden de la aproximación, al mismo tiempo que se beneficia de la flexibilidad del MEF para adaptarse a geometrías complicadas. El MEE se ha utilizado para resolver la ecuación de Helmholtz en combinación con las PMLs por Mehdizadeh y Paraschivoiu en [ 113 ], estudiando con detenimiento la convergencia del método. Le ha seguido trabajos como el de Rong y Xu [ 135 ] o el de Shen y Wang [144] para números de onda elevados. Recientemente, el MEE ha despertado un fuerte interés en diversas disciplinas para la resolución de problemas de propagación de ondas. Uno de los trabajos más recientes es el de He et al. y [73], donde se utiliza el método para resolver la ecuación de Helmholtz en medios estratificados y se combina con técnicas similares a las condiciones de contorno DtN para conseguir condiciones de contorno transparentes. 1.3.5 Estrategias particionadas de acoplamiento fluido-estructura para el problema acústico La interacción dinámica entre fluido y estructura es un problema presente en muchas y variadas disciplinas de la ingeniería, por ejemplo la acústica. Encontrar métodos eficientes que permitan calcular este tipo de problemas es de suma importancia y actualidad, pues 14 Capítulo 2. Propagación de ondas Derivando (2.16) con respecto al tiempo y haciendo uso de la ecuación (2.15) , se obtiene la ecuación clásica de onda: ∂2P0 ∂t2−c2 o∇2P0=0.(2.16) La solución de esta expresión puede ser obtenida mediante el método de separación de variables, comprobándose que es de tipo armónica. Considerando que la perturbación se producen con una frecuencia ω , las variaciones de presión y velocidad pueden ser representadas de la siguiente manera: P0=peiωt(2.17) v0=veiωt(2.18) donde p y v son las distribuciones (complejas) espaciales de la presión y velocidad. Por simplicidad en la notación se ha reutilizado el término v. Introduciendo (2.17) en la ecuación (2.16) se obtiene la conocida ecuación de Helmholtz: ∇2p+k2p=0,(2.19) donde k es el número de onda k=ω/c0 y el operador ∇2 el operador Laplaciano. El campo de velocidades puede ser extraído a partir de las presiones, utilizando la ecuación de conservación de la cantidad de movimiento (2.15): v=−1 iωρ0 ∇p,(2.20) pudiendo aproximarse la velocidad normal al contorno como: vn=−1 iωρ0 ∇p·n,(2.21) donde n es el vector normal a la superficie. Esta relación permite aplicar condiciones de contorno de tipo Neumann o Robin a la ecuación (2.19). 2.3 Propagación del oleaje En esta sección se deriva la EPS, presentada por Berkhoff en 1972 en [ 18 ] y que cuatro años más tarde detalló en su tesis doctoral [ 19 ]. Este modelo lineal asume las hipótesis de fluido homogéneo, no viscoso, incompresible y flujo irrotacional. Esta ecuación es válida para olas de bajo peralte propagándose en zonas dispersivas con pendientes en la batimetría de hasta 1:3, como expuso Booij en 1983 [ 24 ]. La EPS también es conocida con el nombre con el que se presentó inicialmente: "Modelo de Refracción-Diffracción", pero no sólo es capaz de captar estos efectos, pues también capta el efecto de asomeramiento y la reflexión producida por obstáculos. Para la obtención de la EPS en esta sección se ha seguido la tesis de Berkhoff [ 19 ]. Existe sin embargo otra forma de derivar la ecuación de pendiente suave, en la que se 2.3 Propagación del oleaje 15 h(x,y) x z H ζ(x,y,t) L Figura 2.1 Sección esquemática en el plano x−z de olas propagándose sobre un fondo variable. procede de una forma variacional o Hamiltoniana de la que se pueden encontrar más detalles en [48]. El esquema del problema de propagación sobre un fondo variable se muestra en la Figura 2.1, donde además se definen algunas de las variables que se utilizarán a continuación para la derivación de la EPS. La elevación de la superficie del agua viene dada por ζ(x,y,t) , siendo H y L la amplitud y la longitud de onda, respectivamente. Por último, la cota cero del sistema de referencia se sitúa sobre la superficie media del agua y la profundidad se denomina como h(x,y). 2.3.1 Ecuaciones de gobierno para el fenómeno del oleaje Bajo la hipótesis de irrotacionalidad, la velocidad de las partículas en un fluido se puede expresar como el gradiente de un potencial de velocidad, v=∇Φ,(2.22) convirtiendo la ecuación de conservación de la masa (2.3) en la ecuación de Laplace si consideramos un fluido incompresible: ∇2Φ=0.(2.23) 16 Capítulo 2. Propagación de ondas La ecuación de la conservación de la cantidad de movimiento es, bajo estas hipótesis, la ecuación de Bernouilli, Φ,t+1 2(∇Φ ·∇Φ) = 1 ρP+C,(2.24) donde P(x,y,z,t) es la presión en el fluido, C es una función escalar de cuyo gradiente son obtenidas las fuerzas de volumen y Φ,t es la derivada temporal del potencial. Por simplicidad en la notación, de aquí en adelante de escribirá así la derivada. 2.3.2 Condiciones de contorno La solución del problema tridimensional descrito por la ecuación de Laplace (2.23) y la de Bernouilli (2.24) en términos del potencial de velocidad y la presión, requieren de la aplicación de las condiciones de contorno tanto en el fondo como en la superficie libre. Estas condiciones de contorno son de dos tipos: las cinemáticas, que afectan al movimiento de las partículas, y las dinámicas, donde la restricción se ejerce mediante la distribución de presiones. En la superficie del fondo y en la superficie libre, se ha de exigir que la trayectoria de las partículas del fluido no tengan una componente normal a dicha superficie, es decir, su trayectoria sea tangente. El incumplimiento de esta condición implicaría la transferencia de masa a través de los contornos. Estas restricciones son conocidas como condiciones de contorno cinemáticas. Por otro lado, la presión que actúa sobre la superficie libre es la presión atmosférica ( Pa ), esta condición de contorno dinámica permitirá establecer una relación directa entre la energía cinética y potencial en la superficie. Condiciones de contorno cinemáticas Si se considera un fondo impermeable tal que podamos considerar que la velocidad de las partículas en la dirección normal a su superficie sea nula, el producto escalar de la normal a la superficie y la velocidad de las partículas ha de ser nulo: Φ,xh,x+Φ,yh,y+Φ,z=0; en z=−h(x,y),(2.25) proporcionándose así la expresión que representa la condición cinemática en el fondo. La condición cinemática en la superficie se obtiene imponiendo la no transferencia de masa a través de la interfase agua-aire. Esto implica que las partículas localizadas en la superficie habrán de seguir su movimiento, permaneciendo en la superficie en todo momento. Siguiendo el esquema presentado en la Figura 2.1, la superficie libre viene dada por z=ζ(x,y,t) . Por lo tanto, para cualquier partícula localizada en la superficie, su velocidad en la dirección vertical vendrá dada por: z,t=ζ,t+ζ,xx,t+ζ,yy,ten z=ζ,(2.26) expresión en la que introduciendo (2.22) se obtiene: ζ,t+Φ,xζ,x+Φ,yζ,y=Φ,zen z=ζ,(2.27) 2.3 Propagación del oleaje 17 que es la condición de contorno cinemática en la superficie. Condición de contorno dinámica La condición de contorno dinámica se obtiene introduciendo en la ecuación de Bernouilli (2.24) la presión atmosférica en la superficie, por lo que: Φ,t+1 2(∇Φ)2+gζ=Pa ρ,(2.28) donde el potencial C toma el valor gζ . Si además consideramos Pa la presión de referencia de forma que en la superficie podamos considerar presión nula, la ecuación anterior puede reescribirse como: Φ,t+1 2(∇Φ)2+gζ=0.(2.29) Resumen del problema no lineal de propagación de oleaje Recapitulando las ecuaciones obtenidas anteriormente, el problema de propagación se puede resumir en que la ecuación de gobierno del problema es la ecuación de Laplace (2.23) , a la que es necesario imponer las condiciones de contorno cinemáticas en el fondo (2.25) y en la superficie (2.27) y la condición de contorno dinámica en la superficie. El problema completo es, por tanto: ∇2Φ=0; −h(x,y)≤z≤ζ(x,y,t);(2.30) ζ,t+Φ,xζ,x+Φ,yζ,y−Φ,z=0; z=ζ(x,y,t),(2.31) Φ,t+1 2(∇Φ)2+gζ=0; z=ζ(x,y,t),(2.32) Φ,z+Φ,xh,x+Φ,yh,y=0; z=−h(x,y),(2.33) que es un problema no lineal en el que el contorno z=ζ(x,y,t) está en continuo movimiento. 2.3.3 Linealización de las condiciones de contorno A pesar de que el conjunto de ecuaciones (2.30-2.33) constituye un problema no lineal, se pueden aplicar ciertas restricciones adicionales que simplifiquen las condiciones de contorno, permitiendo así la linealización del problema. Para simplificar este proceso, las variables son adimensionalizadas utilizando las siguentes transformaciones [114]:                Φ=HωL 2πΦ0, (x,y,z,h) = L 2π(x0,y0,z0,h0), t=1 ωt0, ζ=Hζ0; (2.34) donde H=|ζ| , L y ω la longitud de onda y frecuencia angular respectivamente. Introduciendo estas nuevas variables el conjunto de ecuaciones (2.30) - (2.33) se reescribe de 18 Capítulo 2. Propagación de ondas manera adimensional como: ∇02Φ0=0; −h(x,y)≤z≤ζ(x,y,t);(2.35) ζ0 ,t0+εΦ0 ,x0ζ0 ,x0+Φ0 ,y0ζ0 ,y0=Φ0 ,z0;z0=εζ0(x0,y0,t0),(2.36) Φ0 ,t0+1 2ε(∇0Φ0)2+κζ0=0; z0=εζ0(x0,y0,t0),(2.37) Φ0 ,z0+Φ0 ,x0h0 ,x0+Φ0 ,y0h0 ,y0=0; z0=−h0(x,y),(2.38) siendo ε=2πH/L y κ=2πg/(ω2L) . Además, el potencial de velocidad Φ0 en la superficie z0=εζ0 puede ser aproximado mediante un desarrollo en serie de Taylor en torno a z0=0 : Φ0(x0,y0,εζ0,t0) = Φ0z0=0+εζ0Φ0 ,z0z0=0+O(ε2),(2.39) lo que permite despreciar aquellos términos de la serie multiplicados por ε si las olas poseen un peralte bajo ( H/L1 ). Si introducimos esta aproximación en (2.36) y en (2.37) las ecuaciones pueden ser reescritas como: ζ0 ,t0=Φ0 ,z0;z0=0,(2.40) Φ0 ,t0+κζ0=0; z0=0,(2.41) donde ahora las condiciones de contorno quedan impuestas sobre el plano medio de la superficie, convirtiendo ambas ecuaciones en ecuaciones lineales. Derivando (2.40) con respecto al tiempo podemos combinarla con (2.41), resultando la expresión: Φ0 ,t0t0+κΦ0 ,z0=0,en z0=0.(2.42) La elevación de la superficie queda ligada al potencial de velocidad de manera proporcional en z=0como: ζ(x,y,t) = −1 gΦ,t,(2.43) relación que proviene de la versión linealizada de la ecuación de contorno dinámica en la superficie del agua (2.41). 2.3.4 El problema lineal Recapitulando las ecuaciones (2.35-2.43) y volviendo a las variables originales del problema de potencial, el sistema no lineal (2.30-2.33) queda simplificado de la forma: ∇2 HΦ+Φ,zz =0; −h(x,y)≤z≤0,(2.44) Φ,z+1 gΦ,tt =0; z=0,(2.45) Φ,z+Φ,xh,x+Φ,yh,y=0; z=−h(x,y),(2.46) 2.3 Propagación del oleaje 19 donde ∇2 H=∂,xx +∂,yy . La introducción de este operador servirá para facilitar la notación en los apartados siguientes. 2.3.5 Propagación de ondas armónicas y dispersión A pesar de ser un problema tridimensional, la propagación del oleaje se produce en el plano horizontal preferencialmente. A partir del problema de propagación linealizado (2.44-2.46) es posible eliminar la dependencia de la variable espacial z , quedando un problema bidimensional en el plano x−y . Para llevar a cabo tal propósito se haya la solución del sistema de ecuaciones mediante el método de separación de variables, de forma que: Φ(x,y,z,t) = Z(h,z)φ(x,y)T(t),(2.47) donde Z(h,z) representa el perfil que adopta el potencial a lo largo de dicha coordenada, φ(x,y)es el potencial de velocidad en el plano x−yyT(t)es la función temporal. Si introducimos esta aproximación de Φ en (2.44) y asumimos que T(t)6=0 se obtiene: 1 φ∇2 Hφ=−1 ZZ,zz;(2.48) siendo ambos lados de la ecuación igual a una constante, −k2 . La solución para la ecuación unidimensional 1 ZZ,zz =k2(2.49) será de la forma Z(z) = Aekz +Be−kz,(2.50) siempre y cuando k2>0 . Si suponemos puntualmente un fondo plano, la condición de contorno (2.46) pasará a ser Z,z=0en z=−h,(2.51) que aplicada a la solución (2.50) resulta que: Ae−kh −Bekh =0.(2.52) Asignando Γ/2=Ae−kh =Bekh e introduciéndolo en (2.49) , obtenemos, tras tener en cuenta algunas relaciones trigonométricas: Z(h,z) = Γcosh(k(z+h)).(2.53) Introduciendo la aproximación del potencial descrita en (2.47) , incluyendo la solución obtenida para Z(h,z) , dentro de la condición de contorno linealizada de la superficie libre (2.42) se obtiene: T,tt +gktanh(kh)T=0,(2.54) ecuación que presenta una solución del tipo: T(t) = Qeiωt,(2.55) 20 Capítulo 2. Propagación de ondas donde Q es una constante que por conveniencia tomará valor unidad. Introduciendo (2.55) en (2.54) resulta: ω2=gktanh(kh),(2.56) que es la conocida como Ecuación de Dispersión. Esta expresión implica que olas con diferente periodo viajarán a diferente velocidad para una misma profundidad confiriéndole un carácter dispersivo a la propagación. La Ecuación de Dispersión tiene como solución una sola raíz real y positiva e infinitas complejas. Sin embargo, en esta tesis se considera únicamente la raíz real al igual que Berkhoff [ 19 ]. Las soluciones complejas de la ecuación componen el conjunto de los llamados modos evanescentes (para más detalles véase [ 48 ]). Estos modos pueden ser usados para obtener una mejor aproximación del problema, como por ejemplo la Ecuación de Pendiente Suave Extendida (EPSE) [ 112 , 130 ] o la versión de Belibassakis en [ 10 ], donde se presenta un modelo de ecuaciones acopladas que incluyen además un modo asociado a la pendiente del fondo. La velocidad de propagación del oleaje ( c ) puede ser obtenida de la Ecuación de Dispersión: c=ω k=rg ktanh(kh)(2.57) que se aproxima en los límites de kh como: c= (gh)1/2,kh 1,(2.58) c= (g/k)1/2,kh 1.(2.59) Analizando estas expresiones se observa que la velocidad de propagación del oleaje depende de la gravedad y de la profundidad cuando la longitud de onda es grande y la profundidad es pequeña ( kh 1 ). Lo cual quiere decir que en estas circunstancias la propagación pierde su carácter dispersivo. Por otro lado, cuando se trata de mar de fondo u olaje en aguas profundas ( kh 1 ), la velocidad de propagación depende del número de onda y por lo tanto del periodo, adquiriendo el carácter dispersivo mencionado. Un concepto importante que conviene introducir es la velocidad de grupo ( cg ), que se define como cg=dω/dk . Para fondo plano puede ser obtenida derivando la Ecuación de Dispersión: cg=c 21+2kh sinh2kh.(2.60) La importancia de la velocidad de grupo radica en que indica la velocidad a la que se propaga la energía. El comportamiento dispersivo de la velocidad y de la velocidad de grupo se muestra en la Figura 2.2, donde se aprecia como la razón entre ambas es la unidad en aguas profundas ( kh 1 ), mientras que en aguas someras ( kh 1 ) la velocidad de grupo tiende a la mitad de la velocidad de fase. Para una explicación más extensa sobre el tema se refiere al lector a [116]. Puesto que la solución de (2.49) se cumple para cualquier valor de su constante Γesta se selecciona de forma que Z(0) = 1, recayendo así el valor de la amplitud de la onda en 2.3 Propagación del oleaje 21 Figura 2.2 Curvas de dispersión para una onda propagándose sobre un fondo de profundidad variable. la superficie z=0sobre el valor de la amplitud de la función φ(x,y), por lo tanto: Γ=1 cosh(kh),(2.61) resultando que Φ(x,y,z,t) = cosh(k(z+h)) cosh(kh)φ(x,y)eiωt.(2.62) 2.3.6 Derivación de la ecuación de pendiente suave La EPS es una aproximación del fenómeno de propagación de ondas siempre y cuando los efectos no lineales no sean significativos y que la propagación se produzca en zonas donde la reflexión producida por los cambios de profundidad sea igualmente despreciable. Esto último se traduce cuantitativamente en que la relación entre la profundidad y el número de onda ha de ser tal que: ∇Hh kh 1.(2.63) Las ecuaciones que reproducen el fenómeno de propagación, incluyendo todos los efectos son: ∇2ϕ=0; −h(x,y)≤z≤0; (2.64) ϕ,z−ω2 gϕ=0; z=0,(2.65) ϕ,z+ϕ,xh,x+ϕ,yh,y=0; z=−h(x,y),(2.66) 22 Capítulo 2. Propagación de ondas D α σ=tan(α) D Figura 2.3 Definición de los parámetros utilizados para la adimensionalización del problema. que son las ecuaciones lineales (2.44) - (2.46) considerando la solución armónica obtenida en (2.62) y donde ϕ=Z(h,z)φ(x,y) . En esta sección se siguen los pasos descritos por Berkhoff en su tesis, sin embargo, existen otras derivaciones como las de Smith y Spring [148] o la de Dingemans [48]. Las magnitudes espaciales se pueden adimensionalizar mediante la longitud característica lo=g/ω2 , que posee un valor del orden de la longitud de onda en la superficie para el periodo estudiado. Este cambio de variable permite eliminar los coeficientes presentes en la condición de contorno en la superficie. Las nuevas variables consideradas son: (x0,y0,z0,h0) = 1 lo (x,y,z,h),(2.67) y por lo tanto, ∇ϕ=1 lo ∇0ϕ.(2.68) La EPS se plantea para pequeñas variaciones en la pendiente del fondo, como su propio nombre indica, por lo que el gradiente de la profundidad con respecto a las coordenadas globales ( x,y ) será ,bajo estas condiciones, despreciable. Para adimensionalizar y al mismo tiempo expresar este gradiente como un valor en torno a la unidad, se introducen las nuevas variables adimensionales ˆxeˆy, con valores: (ˆx,ˆy) = σ D(x,y);(2.69) donde D es la profundidad media y σ es la pendiente media en una zona de longitud D (véase Figura 2.3). Teniendo en cuenta estas nuevas variables el gradiente de la profundidad adimensional se podrá expresar como: ∇h0=σ Dloˆ ∇h0.(2.70) 2.3 Propagación del oleaje 23 Introduciendo los nuevos operadores adimensionales en las ecuaciones (2.64) - (2.66) se obtiene: ∇02 Hϕ+ϕ,z0z0=0; −h0(x,y)≤z0≤0,(2.71) ϕ,z0−ϕ=0; z0=0,(2.72) ϕ,z0+εϕ,x0h0 ,ˆx+ϕ,y0h0 ,ˆy=0; z0=−h0(x,y),(2.73) donde ε=σ/µ y µ=D/lo y de nuevo se utiliza el subíndice H para denotar el plano horizontal. Por simplicidad y para una mejor comprensión, se reescribe el conjunto de ecuaciones anterior eliminando las primas: ∇2 Hϕ+ϕ,zz =0; −h(x,y)≤z≤0,(2.74) ϕ,z−ϕ=0; z=0,(2.75) ϕ,z+εϕ,xh,ˆx+ϕ,yh,ˆy=0; z=−h(x,y).(2.76) Siguiendo el mismo procedimiento que en la Sección 2.3.5, la solución del potencial de velocidades ϕ se expresa nuevamente mediante el producto de dos funciones, en este caso dependientes de las variables adimensionales introducidas: ϕ=Z(h,z)φ(x,y,z)(2.77) con Z(z) = cosh(κ(h+z)) cosh(κh),(2.78) donde la variable κ , o número de onda adimensional, puede ser obtenida de la ecuación de dispersión 1=κtanh(κh).(2.79) Es en este punto que la derivación de la EPS diverge del modelo de refracción-asomeramiento y del de difracción puro. En el modelo de refracción el potencial se aproxima como ϕ=HeiS , donde H(z,y,z) es la función asociada a la amplitud del potencial y S(x,y) la fase. Este procedimiento da lugar a prescindir del gradiente horizontal del perfil vertical de velocidades Z(h,z) , lo que conduce a dos ecuaciones desacopladas. Por otro lado, el modelo de difracción puro está basado en la asunción de fondo plano, dando lugar a la ecuación de Helmholtz. Sustituyendo (2.77) en el sistema (2.74-2.76) resulta: Z∇2φ+2∇Z·∇φ+φ∇2Z+φZ,zz +2Z,zφ,z+Zφ,zz =0; −h(x,y)≤z≤0,(2.80) φ,z=0; z=0,(2.81) Zφ,z+εh(∇Z·ˆ ∇h)φ+Z(∇φ·ˆ ∇h)i=0; z=−h,(2.82) donde se ha suprimido el subíndice H del operador ∇ por simplicidad. En el conjunto de ecuaciones anterior, el gradiente en el plano x−y de la función Z(h,z) puede reescribirse 30 Capítulo 3. Aproximaciones numéricas Sin embargo también son patentes algunas importantes desventajas frente a otros métodos. La modelización mediante el MEC de un problema regido por una ecuación diferencial determinada requiere el uso de la solución fundamental para dicho problema, aunque desafortunadamente las soluciones fundamentales conocidas son limitadas y en su mayoría, se restringen a ecuaciones diferenciales lineales con coeficientes constantes. Consecuentemente problemas no lineales o con inhomogeneidades suelen quedar fuera del ámbito de aplicación del MEC a no ser que se use algún tipo de tratamiento especial. Otra desventaja del MEC es que el sistema de ecuaciones resultante de su aplicación es no-simétrico y lleno, lo cual aumenta considerablemente el coste computacional del proceso de resolución. Sin embargo, en muchas ocasiones esto se ve sobradamente compensado con el menor esfuerzo requerido en la generación de la malla y el reducido número de incógnitas resultantes. 3.2.2 Aplicación del MEC a la ecuación de Helmholtz bidimensional El problema de Helmholtz es uno de los principales marcos de aplicación del MEC, puesto que se trata de un problema de condición de frontera o contorno. Su descripción puede ser encontrada en cualquier libro de referencia, como lo son [ 29 , 49 ], existiendo libros especializados como [166, 94], dedicados a la aplicación del MEC a la acústica lineal. Ecuación integral de contorno Si la ecuación de Helmholtz se cumple en un dominio Ω , entonces, integrando en dicho dominio la ecuación completa multiplicada por una función ψ(x) y definida en el mismo dominio resultará que ZΩ(∇2φ+k2φ)ψdΩ=0,(3.11) relajando las exigencias del problema. Suponiendo ahora una función de ponderación ψ(x) que también cumple la ecuación diferencial y aplicando dos veces el teorema de Gauss convenientemente, se obtiene la Segunda identidad de Green oTeorema de Reciprocidad, punto de partida del MEC: ZΩ(∇2ψ+k2ψ)φdΩ+ZΓ ψ∇φ·ndΓ=ZΩ(∇2φ+k2φ)ψdΩ+ZΓ ∇ψ·ndΓ,(3.12) expresión para la que de momento sólo se ha exigido que la función de ponderación cumpla la ecuación diferencial. Sin embargo, si se considera que dicha función es la función de Green del problema, la ecuación anterior podrá reescribirse en el contorno como C(x0)φ(x0)+ZΓ ψ(x,x0)∇φ(x)·ndΓ−ZΓ(∇ψ(x,x0)·n)φ(x)dΓ=0(3.13) en ausencia de fuentes internas o sumideros. La ecuación anterior contiene el término libre C(x0) , que depende de la localización del punto de colocación x0 , siendo x el punto de observación. Para más detalles sobre la obtención del término libre se refiere al lector a la basta bibliografía existente sobre el MEC, por ejemplo [ 29 ]. Los valores que toma el término libre en el problema bidimensional son: C(x0) = 1 para x0 situado en el interior del dominio, C(x0) = 1/2 cuando x0 está situado sobre el contorno Γ y éste tiene una forma 3.2 Método de los Elementos de Contorno 31 suave y C(x0) = θ/(2π) cuando x0 además de estar situado sobre el contorno está sobre un vértice del mismo, siendo θel ángulo interior. En el caso de que existan fuentes internas la ecuación integral (3.12) quedará de la forma: C(x0)φ(x0)+ZΓ ψ(x,x0)∇φ(x)·ndΓ−ZΓ(∇ψ(x,x0)·n)φ(x)dΓ=ψ(x0,y)(3.14) donde ψ(x0,y) es la función de Green simulando una fuente posicionada en y y evaluada en el punto de colocación x0 . En el caso de tratarse de un problema de scattering en un dominio abierto, el término ψ(x0,y)tomará el valor de la onda incidente [166]. Discretización de la ecuación integral de contorno Para la resolución numérica de la ecuación integral de contorno (3.14) el dominio es dividido en un número finito de elementos de forma que Γ= ne [ e=1 Γj, ne \ e=1 Γj=/0,(3.15) siendo ne el número de elementos. Las variables φ y q son aproximadas mediante unas funciones de interpolación a partir de los valores nodales: φ= nn ∑ i=1 Niφi,q= nn ∑ i=1 Niqi,(3.16) donde Ni es la función interpolante correspondiente a cada nodo i y nn es el número de nodos existente el contorno. Introduciendo (3.15) y (3.16) en la ecuación (3.14) , ésta puede ser reescrita en su forma discreta: Ciδi jφj+ ne ∑ e=1ZΓe ∂ψ(x,x0) ∂nNjφjdΓe− ne ∑ e=1ZΓe ψ(x,x0)NjqjdΓe=ψi,(3.17) siendo δi j la delta de Kronecker. De forma más compacta, la expresión anterior puede ser representada en su forma matricial: Hφ φ φ−Gq =φ φ φin,(3.18) donde los coeficientes de las matrices HyGvienen dados por: Hi j =Ciδi j + ne ∑ e=1ZΓe ∂ψ(x,x0) ∂nNjφjdΓe,(3.19) Gi j = ne ∑ e=1ZΓe ψ(x,x0)NjqjdΓe,(3.20) y el vector φ φ φin contiene los valores del potencial incidente sobre el contorno. 32 Capítulo 3. Aproximaciones numéricas 3.2.3 Solución fundamental de la ecuación de Helmholtz La solución fundamental para la ecuación de Helmholtz en dos dimensiones se obtiene de la solución de la ecuación: ∇2ψ+k2ψ+δ(x,x0) = 0,x∈R2,(3.21) donde δ(x,x0) es la función Delta de Dirac aplicada en el punto x0 . La condición de contorno del problema es la condición de radiación en el infinito o condición de radiación de Sommerfeld: l´ ım r→∞r∂φ ∂r−ikφ=0.(3.22) La solución de (3.21) puede ser encontrada en una extensa bibliografía, como por ejemplo [ 29 , 166 ]o[ 102 , 163 ], donde se obtiene con más detalle. La función de Green viene dada por: ψ(x,x0) = i 4H(1) 0(kr),(3.23) donde r=|x−x0| y H(1) 0 es la función de Hankel de primera especie y de orden zero. La derivada de esta función de Green en dirección radial es ψ,r(x,x0) = −i 4kH(1) 1(kr),(3.24) siendo H(1) 1 la función de Hankel de primera especie de orden uno. El conjunto de las dos expresiones anteriores, (3.23) y (3.24) , conforman la solución fundamental necesaria para la formulación MEC del problema de Helmholtz en un medio homogéneo. Para más detalles sobre su obtención nos referimos aquí a [166]. 3.2.4 Solución fundamental para número de onda variable En el Artículo I anejo a esta tesis se ha propuesto una solución fundamental para el problema de propagación de oleaje sobre zonas con batimetría variable, ampliando así la aplicación del MEC a problemas más complejos. Esta solución fundamental está basada en la función de Green propuesta por Belibassakis en [ 9 ] y es válida no sólo para la ecuación mild-slope, sino también para la mild-slope modificada [ 35 ] o para cualquier ecuación diferencial que pueda ser reescrita como la ecuación de Helmholtz con número de onda variable. Además, esta solución fundamental se ha introducido en una formulación MEC y se han resuelto varios ejemplos para validarla. Las principales contribuciones de Belibassakis en [ 9 ] fueron la obtención de una solución numérica de la función de Green ψ(x,x0) para batimetría variable en una dirección y el estudio de su comportamiento asintótico en el campo lejano. Su desarrollo estaba enfocado a ampliar la gama de problemas que podían ser resueltos mediante formulaciones como la de Berkhoff [ 16 , 19 ], en las que se trata el problema de difracción del oleaje provocado por una región con una batimetría irregular sobre un fondo plano, o como la de Xu [ 167 ] donde se trata el problema de difracción en acústica cuando existe una zona inhomogénea no-estratificada en un medio estratificado. El interés por una función de Green de estas 3.2 Método de los Elementos de Contorno 33 características es patente y se pueden encontrar algunas soluciones para casos particulares para dos y tres dimensiones en [102, 76, 105, 111]. Al ser la función de Green propuesta por Belibassakis válida para cualquier perfil batimétrico, el trabajo se sitúa como un excelente punto de partida para una formulación MEC del problema de propagación de oleaje. Sin embargo, Belibassakis en su trabajo no obtuvo las derivadas de la función, dejando incompleta la solución fundamental necesaria para la formulación MEC del problema. El principal aporte del Artículo I es la obtención de dichas derivadas, siendo casi inmediata su obtención en la dirección en la que cambia la batimetría pero complicándose en muchos aspectos en la dirección paralela a las isóbatas. Si bien es cierto que la solución fundamental propuesta está limitada a batimetrías compuestas por isóbatas rectas y paralelas, una formulación MEC que permita resolver problemas con este tipo de inhomogeneidad sigue siendo de enorme importancia, puesto que se está resolviendo una forma batimétrica de lo más común en problemas de propagación de oleaje en zonas costeras [ 121 ]. Gracias a esta aportación, el MEC puede ser usado por sí solo como una herramienta precisa y eficiente para la realización de modelos de propagación de oleaje. Para problemas más complejos, en los que el interés se centrara en zonas con una batimetría más compleja, el MEC puede ser acoplado con un modelo basado en el MEF que reproduzca esta zona de manera más precisa, aportando una mejora significativa en las condiciones de contorno. Considerando que la dirección en la que cambia el número de onda es x , la función de Green ψ(x,x0)ha de satisfacer: ∇2ψ+ˆ k(x)2ψ+δ(x−x0) = 0(3.25) donde el número de onda es obtenido mediante la ecuación de dispersión (2.56) a partir de la profundidad y del periodo del oleaje. El perfil a lo lago de xvaría de la forma: ˆ k(x) =    ˆ ka,x<a, ˆ k(x),a≤x≤b, ˆ kb,x>b, (3.26) siendo ˆ ka y ˆ kb constantes, lo cual permitirá que la condición de Sommerfeld sea satisfecha exactamente, como veremos más adelante. El perfil tipo utilizado, descrito en (3.26) se muestra en la Figura 3.1. La ecuación (3.25) puede ser transformada en una familia de problemas unidimensionales mediante la transformada de Fourier en el eje y . La relación entre el potencial original y su transformado viene dada por: Ψ(x,x0;ξ) = Z∞ −∞ ψ(x,x0;ˆ k)e−iyξdy (3.27) y ψ(x,x0;ˆ k) = 1 2πZ∞ −∞ Ψ(x,x0;ξ)eiyξdξ,(3.28) quedando definida la transformada de Fourier (F(•)) en (3.27) y la antitransformada (F−1(•)) en (3.28) y donde el potencial transformado está representado por Ψ(x,x0;ξ) . A partir de la expresión integral del potencial dada en la ecuación anterior, las expresiones 34 Capítulo 3. Aproximaciones numéricas x ˆ k(x) xaxb ˆ kb ˆ ka ˆ k Figura 3.1 Variación del número de onda con la distancia x. de las derivadas en ambas direcciones físicas, xey, son fácilmente obtenidas, quedando: ψ,x(x,x0;ˆ k) = 1 2πZ∞ −∞ Ψ,x(x,x0;ξ)eiyξdξ,(3.29) ψ,y(x,x0;ˆ k) = 1 2πZ∞ −∞ iξΨ(x,x0;ξ)eiyξdξ.(3.30) La función de Green (3.28) y sus derivadas completan la solución fundamental de forma que pueda ser introducida en una formulación MEC. El problema transformado en el espacio de Fourier Introduciendo (3.28) en (3.25) y teniendo en cuenta la condición de radiación en el infinito (3.22), se obtiene la familia de problemas unidimensionales: Ψ,xx +κ2(x)Ψ+δ(x−x0) = 0parax∈[a,b],(3.31) Ψ,x+iα(ξ)Ψ=0en x=a, Ψ,x−iβ(ξ)Ψ=0en x=b,(3.32) donde κ2(x) = ˆ k2 c(x)−ξ2(3.33) y los coeficientes α y β son α(ξ) = (ˆ k2 a−ξ2)1 2 y β(ξ) = (ˆ k2 b−ξ2)1 2 respectivamente. La solución en los intervalos semi-infinitos (−∞,a] y [b,∞) puede ser obtenida mediante la solución del problema para número de onda constante una vez ha sido calculada la 3.2 Método de los Elementos de Contorno 35 solución en el intervalo (a,b): Ψ(x,x0,ξ) = Ψ(xa,x0,ξ)e−iα|x−xa|,x≤a,(3.34) Ψ(x,x0,ξ) = Ψ(xb,x0,ξ)eiβ|x−xb|,x≥b.(3.35) Con el fin de obtener de nuevo el potencial original ψ , la familia de problemas (3.31, 3.32) ha de ser resuelta para todo valor de ξ∈R . Sin embargo, es conocido que la solución presenta cuatro polos situados sobre el eje real, simétricos respecto al origen, localizados en ±ˆ ka , y ±ˆ kb [ 9 ], que obligan a recurrir a la integración el plano complejo para poder evaluar la integral de la antitransformada. En vista de lo anterior, existen dos cuestiones principales a las que hacer frente: •Cómo resolver el problema unidimensional descrito en (3.31) y (3.32). • Cuál es el camino de integración en el plano complejo más adecuado para llevar a cabo la antitransformada de Fourier. La primera cuestión no será abordada en este capítulo. La solución del problema unidimensional puede ser hallada como lo hizo Belibassakis [ 9 ], mediante un esquema MDF, o como se ha resuelto en el Artículo I , donde la solución se obtiene mediante el MEF. La elección del camino de integración se basa fundamentalmente en tres criterios: la simetría y el comportamiento asintótico de la solución del problema unidimensional, la posibilidad de usar el algoritmo de la Transformada Rápida de Fourier (FFT) y el criterio de signos elegido para el sentido positivo de la propagación. Integración en el plano complejo El camino de integración de las expresiones (3.28-3.30) en el plano complejo ( ξ∈C|ξ= ξ1+iξ2 ) puede ser expresado como la diferencia entre un contorno cerrado y parte del mismo, de forma que ZCS Ψ(x,x0;ξ)eiyξdξ=Iγ Ψ(x,x0;ξ)eiyξdξ−ZCR Ψ(x,x0;ξ)eiyξdξ,(3.36) siendo el contorno γ=CS∪CR , donde CS≡{ξ= [−R,R]} y CR≡ξ=Reiθ,∀θ∈[0,π] , tal y como se muestra en la Figura 3.2. Según el lema de Jordan, haciendo el límite cuando Rtiende a infinito, la integral sobre el contorno semicircular se hace nula: l´ ım R→∞ZCR Ψ(x,x0;ξ)eiyξdξ=0,(3.37) por lo que Iγ Ψ(x,x0;ξ)eiyξdξ=l´ ım R→∞ZCS Ψ(x,x0;ξ)eiyξdξ.(3.38) La importancia que tiene escribirlo de esta manera es que se puede aplicar el teorema de Cauchy a la integral en el contorno cerrado γ permitiendo desviar el camino de integración, situado en el eje real, al plano complejo. 36 Capítulo 3. Aproximaciones numéricas ˆ kb ˆ ka −ˆ ka −ˆ kbξ1 ξ2 R CR CS γ=CS∪CR Figura 3.2 Camino de integración cerrado en el plano complejo. En este punto es importante introducir algunas propiedades de la solución en el espacio transformado de Fourier. De (3.31) y (3.33) se deduce que la solución es simétrica con respecto a ξ, resultando Ψ(x,x0;ξ1+iξ2) = Ψ(x,x0;−ξ1−iξ2),(3.39) por lo que, en caso de considerar un camino de integración simétrico, será suficiente con calcular la solución del problema unidimensional en el semiplano correspondiente a la parte positiva real de ξ . Por otro lado se observa que cuando ξˆ k(x) la solución del problema unidimensional con número de onda variable tiende a la solución del problema con coeficiente constante, siendo por tanto, Ψ(x,x0;ξ)≃i 2κ(x0)eiκ(x0)|x−x0|,ξˆ k(x0),(3.40) definiendo el comportamiento asintótico de la solución. Teniendo en cuenta (3.39) , es razonable escoger un camino de integración simétrico con respecto al origen como el mostrado en la Figura 3.3. La razón de bordear los polos por el segundo y cuarto cuadrante se debe a la física del problema y a la consecuente elección del criterio de signos para el sentido positivo de la propagación. Tomando C como camino de integración, el intervalo [ˆ ka,ˆ kb] queda dentro del contorno γ y el resultado al llevar a cabo la antitransformada es una onda irradiada desde la fuente. Sin embargo, si se integrara dejando el intervalo del eje real [−ˆ kb,−ˆ ka] dentro del contorno γ , el resultado es una onda que se dirige desde el infinito hasta el punto de la fuente, o sumidero en este caso. La diferencia entre incluir una pareja de polos o la otra proviene de la elección en el criterio de signos para la solución del problema unidimensional. En este trabajo un argumento positivo de la función exponencial en (3.40) supone una onda irradiada desde la fuente. Para una descripción más detallada el lector puede consultar [102, 163]. 3.2 Método de los Elementos de Contorno 37 ˆ kb ˆ ka −ˆ ka −ˆ kbξ1 ξ2 C Figura 3.3 Camino de integración simétrico respecto al origen. Los intervalos del eje real en los que la solución del problema unidimensional es singular son evitados pasando por el segundo y el cuarto cuadrante. ˆ kb ˆ ka −ˆ ka −ˆ kbξ1 ξ2 τ ξ1=−Ξ ξ1=Ξ C− 3C− 2 C− 1 C+ 1 C+ 2C+ 3 Figura 3.4 Camino de integración por tramos para resolver de forma numérica. Para la evaluación numérica de la inversa de la transformada necesaria para recuperar el potencial original, la curva de integración C se puede descomponer en los seis tramos rectos que aparecen en la Figura 3.4: C− 3={ξ|ξ1∈(−∞,−Ξ),ξ2=τ}, C− 2={ξ|ξ1∈[−Ξ,0),ξ2=τ}, C− 1={ξ|ξ1=0,ξ2∈[τ,0]}, C+ 1={ξ|ξ1=0,ξ2∈(0,−τ]}, C+ 2={ξ|ξ1∈(0,Ξ],ξ2=−τ}, C+ 3={ξ|ξ1∈(Ξ,∞),ξ2=−τ}, facilitando así las operaciones necesarias. Para valores de τˆ ka y Ξˆ kb , el comportamiento asintótico de la solución descrito en (3.40) decae exponencialmente, pudiendo ser aproximada la solución como: Ψ(x,x0;ξ)≃1 2ξe−ξ|x−x0|,ξ−→∞,(3.41) 38 Capítulo 3. Aproximaciones numéricas x -5 0 5 10 15 ˆ k(x) 0 5 10 Figura 3.5 Función ˆ k(x)utilizada para el caso práctico. y su derivada como Ψ,x(x,x0;ξ)≃−1 2e−ξ|x−x0|,ξ−→∞,(3.42) lo que permite un tratamiento analítico de las integrales sobre los tramos C− 3yC+ 3. Análisis del comportamiento de la función de Green en el espacio transformado Con el fin de ilustrar el comportamiento de la función de Green en el espacio transformado de Fourier se resuelve aquí la familia de problemas unidimensionales (3.31, 3.32) para la función ˆ k(x): ˆ k(x) =    1,x<0, 5.5cos(π(1−x)/10),0≤x≤10, 10,x>10, (3.43) cuyo perfil está representado en la Figura 3.5. La fuente se sitúa en el punto medio del perfil, x=5 . El esquema utilizado para hallar la solución del problema es el MEF descrito en el Artículo I, utilizando 20 elementos lineales por longitud de onda. En primera instancia se resuelve para ξ≡{ξ1∈[0,12],ξ2=0} . El resultado se muestra en la Figura 3.6, donde se puede ver el valor absoluto que toma el potencial transformado Ψ(x,x0,ξ) (a) y su fase (b). Los polos, situados en ξ=1 y ξ=10 , están marcados con sendas líneas rojas. Una de las dificultades que afectan a la integración en la dirección ξ1 es que la función Ψ(x,x0;ξ1−iτ) se vuelve altamente oscilatoria en la zona cercana a los polos para valores de x alejados de la fuente, tal y como se aprecia en la Figura 3.6(b). Para valores de ξ1>ˆ k(x) se observa como el potencial transformado decae exponencialmente para todo x menos en x0 , de ahí la importancia del estudio asintótico de la solución, puesto que esta zona la integral sobre los caminos C− 3 y C+ 3 no será en ningún caso despreciable. Con el fin de evitar polos y discontinuidades en la derivada, se toma un camino por el plano complejo como el mostrado en la Figura 3.4. El camino de integración se adentra por el plano complejo una distancia ξ2=−τ , siendo τ=6ˆ k(x0)/1024 . La elección de este valor no es arbitraria, si no que está elegida acorde a los criterios presentes en [ 9 ] y [ 32 ]. El valor del potencial transformado Ψ(x,x0;ξ1−iτ) para valores de ξ1∈[0,12] se muestra en la Figura 3.7, donde en (a) se representa su valor absoluto y en (b) la fase. Como se puede apreciar, al recorrer la solución por el camino C+ 2 y C+ 3 con un valor de τ pequeño, la solución tiene una forma parecida, pero sin embargo se suaviza y se evitan las singularidades, facilitando así la integración numérica. 3.2 Método de los Elementos de Contorno 39 (a) |Ψ(x,x0;ξ1)|.(b) arg(Ψ(x,x0;ξ1)). Figura 3.6 Solución del potencial transformado Ψ(x,x0,ξ) para el conjunto de puntos ξ≡ {ξ1∈[0,12],ξ2=0} . En (a) se muestra el valor absoluto del potencial, siendo su fase representada en (b). Los polos situados en ξ=1,10 están marcados con lineas rojas. Solución en el espacio real Una vez que la solución del potencial en el espacio transformado ha sido obtenida, la función de Green y sus derivadas en las direcciones x e y pueden ser recuperadas mediante el procedimiento de integración descrito en el Artículo I , donde las antitransformadas (3.28-3.30) son llevadas a cabo con la ayuda de la Transformada Rápida de Fourier Inversa (IFFT), la FFT y las propiedades asintóticas del problema unidimensional. La descripción detallada del proceso necesario para realizar las transformadas y antitransformadas está explicado con detalle en el Artículo I , centrándose esta sección en aportar simplemente una descripción gráfica de la solución del campo potencial y una breve discusión sobre ella. La función de Green del problema de Helmholtz para una distribución del número de onda según (3.43) es representada en la Figura 3.8. La amplitud del campo potencial generado por una fuente puntual situada en (x=5,y=0) se muestra en (a), apreciándose la acumulación de energía en unas ramas con forma parabólica. La Figura 3.8(b) muestra la fase del campo, dejando patente el efecto de refracción. Si bien en un medio homogéneo la isolíneas de fase hubieran sido concéntricas, en este caso la dirección de propagación cambia tanto que prácticamente llegan a ponerse perpendiculares a la dirección en la que cambia el número de onda. Finalmente en la Figura 3.8(c) se representa una perspectiva tridimensional de la parte real del campo potencial. Esta perspectiva ofrece una visión más comprensible del comportamiento de la solución. La derivada en la dirección x del campo potencial ψ(x,x0;ˆ k) se muestra en la Figura 3.9. 46 Capítulo 3. Aproximaciones numéricas para distribuciones de tipo CGL o LGL, conlleva a un requerimiento menor de número de grados de libertad para la solución del problema, las matrices que contienten los coeficientes de los sistemas de ecuaciones generados se van volviendo cada vez más llenas, por lo que la eficiencia en lo que a coste computacional se refiere se ve comprometida. En este sentido Vos et al. estudian en [ 159 ] el grado de aproximación óptimo, concluyendo que en general, aproximaciones del entorno de p=6 alcanzan un buen punto de acuerdo entre precisión y eficiencia. En esta tesis, véase el Artículo II , se han introducido elementos espectrales de alto orden con una distribución nodal LGL tanto en la formulación MEC como en la de tipo MEF para la resolución de la ecuación mild-slope. Los elementos espectrales usados utilizan polinomios de Lagrange como funciones de interpolación de las variables consideradas a partir de los valores nodales. Cuatro son las principales aportaciones de este trabajo: • Los elementos espectrales han sido utilizados para resolver por primera vez la ecuación elíptica mild-slope, considerando batimetrías bariables. • Una formulación MEC para la ecuación mild-slope con una discretización basada en elementos espectrales ha sido presentada y probada en varios ejemplos de referencia. • Una formulación acoplada MEC-MEF basada en una discretización con elementos espectrales ha sido usada para resolver problemas de contorno abierto incluyendo batimetría variable en la formulación MEC. • Por primera vez se ha estudiado la convergencia de un método en los problemas propuestos por Ito y Tanimoto (tren de ondas propagándose sobre un montículo parabólico) en [ 87 ] y por Berkhoff en [ 17 ] (tren de ondas propagándose sobre un montículo elíptico superpuesto en una rampa). El error de las soluciones presentadas es de varios ordenes de magnitud menor que los obtenidos hasta ahora. A continuación, en las subsecciones 3.3.2 y 3.3.3 se presentan algunas de las características básicas de los polinomios de Legendre y se introducen los detalles básicos de la cuadratura de interés en este trabajo, la de LGL. Para más detalles sobre esta cuadratura, los polinomios de Legendre y otras familias de polinomios clásicas en los MMEE se refiere al lector a [ 44 , 30 , 143 ]. En la subsección 3.3.4 se justifica la convergencia de las aproximaciones pseudoespectrales. Para terminar, en la subsección 3.3.5 se incluye una descripción más detallada de los elementos de referencia que en el Artículo II. 3.3.2 Polinomios de Legendre Los polinomios de Legendre ( Pk(x),k=0,1,2,... ) son las autofunciones del problema particular de Sturm-Liouville (1−x2)P0 k(x)0+k(k+1)Pk(x) = 0(3.46) y pueden ser expresados mediante la fórmula de Rodrigues: Pn(x) = 1 2nn! dn dxn(x2−1)n.(3.47) Los polinomios de Legendre poseen las siguientes propiedades: 3.3 Aproximación mediante elementos espectrales 47 •Ortogonalidad: hPm,Pniw=Z1 −1 Pm(x)Pn(x)w(x)dx =γnδmn,γn=2 2n+1,con w(x) = 1, (3.48) hP0 m,P0 niw=Z1 −1 P0 m(x)P0 n(x)dx =n(n+1)γnδmn,con w(x) = (1−x2), (3.49) donde δmn es la delta de Kroneker. •Simetría: Pn(−x) = (−1)nPn(x),Pn(±1) = (±1)n,(3.50) •Recursividad: Pk+1=2k+1 k+1xPk(x)−k k+1Pk−1(x),(3.51) •Recursividad en la derivada: P0 n(x) = n−1 ∑ k=0 k+n impar (2k+1)Pk(x).(3.52) Los cinco primeros polinomios de Legendre se pueden ver en la Figura 3.11, siendo sus derivadas representadas en la Figura 3.12. Los elementos espectrales más comunes están basados en los polinomios de Chebyshev y en los polinomios de Legendre, ambos familias particulares de los polinomios de Jacobi. La principal característica de estos elementos de alto orden es que sus nodos están localizados en los puntos de integración de la cuadratura en cuestión: Legendre-Gauss (LG), Legendre- Gauss-Lobatto o Chebyshev-Gauss-Lobatto. En esta tesis se han incorporado los elementos espectrales a la formulación MEF y MEC descritas anteriormente. 3.3.3 Cuadratura de LGL La integral de una función f(x) en el intervalo x∈[−1,1] puede ser aproximada mediante la cuadratura de LGL, de la forma: Z1 −1 f(x) = n ∑ j=0 f(xj)wj,(3.53) donde los puntos de la cuadratura xjson las raíces del polinomio: (1−x2)P0 n(x),(3.54) 48 Capítulo 3. Aproximaciones numéricas Figura 3.11 Polinomios de Legendre de orden p=0,1,...,4. Figura 3.12 Derivada de los polinomios de Legendre de orden p=0,1,...,4. siendo Pn el polinomio de Legendre de orden n . Los pesos ( wj ) vienen dados por la expresión: wj=2 p(p+1) 1 [Pn(xj)]2,j=0,1,...,n.(3.55) Esta cuadratura es exacta para polinomios de orden hasta 2n−1. 3.3 Aproximación mediante elementos espectrales 49 3.3.4 Aproximación polinomial óptima La aproximación polinomial es crucial en los métodos computacionales, puesto que en realidad son las únicas funciones que un ordenador puede evaluar de forma exacta. La aproximación de una functión f(x) en el intervalo x∈[−1,1] puede llevarse a cabo mediante polinomios de Lagrange, de forma que: IX nf(x) = n ∑ k=0 f(xk)lk(x),(3.56) siendo X={x0,x1,...,xn} el conjunto de n+1 puntos donde f(xk) = IX nf(xk) . El polinomio aproximante de Lagrange lk(x) correspondiente a cada punto del conjunto X viene dado por: lk(x) = n ∏ j=0,j6=k x−xj xk−xj .(3.57) El error que se produce al usar este tipo de interpolación para una función continua hasta su derivada (n+1) -ésima puede ser acotado superiormente mediante el teorema del error de interpolación de Cauchy: e(x) = kf(x)−IX nf(x)k∞≤kf(n+1)(x)k∞ (n+1)!kWX n+1(x)k∞,(3.58) donde kg(x)k∞=m´ ax x∈[−1,1]|g(x)| , g(n)(x) indica la derivada enésima de g(x) y WX n+1(x) es un polinomio de Lagrange de orden n+1 con raíces situadas en el conjunto de puntos X : WX n+1(x) = n ∏ k=0 (x−xk),xk∈X.(3.59) De la ecuación (3.58) , que es una expresión exacta, se deduce que la única manera posible de disminuir el límite superior por el que viene acotado el error e(x) es minimizando kWX n+1(x)k∞ , lo que se consigue con una correcta elección del conjunto de puntos contenidos en X . La opción más inmediata a la hora de escoger el conjunto de puntos X es mediante una distribución equiespaciada, como ocurren en la versión p del MEF. sin embargo, es conocida la existencia del fenómeno de Runge [ 161 ], que pone de manifiesto la falta de convergencia de este tipo de distribución nodal para funciones del tipo f(x) = 1/(1+cx2) , siendo c una constante. La distribución nodal óptima que minimiza kWX n+1(x)k∞es la de Chebyshev-Gauss (CG) [52] que proporciona el límite superior kWX n+1(x)k∞=1 2n,(3.60) En la Figura 3.13 se muestra la evolución del valor de kWX n+1(x)k∞ en función del grado de los polinomios interpolantes n para las distribuciones nodales de LGL, CGL, CG y equiespaciada. Como se puede ver en la figura la convergencia de la norma-infinito para 50 Capítulo 3. Aproximaciones numéricas Figura 3.13 Convergencia de kWX n+1(x)k∞para diferentes distribuciones nodales. las distribuciones LGL, CGL y CG es exponencial, mientras que para la distribución equiespaciada es una expresión algebraica dependiente de la distancia entre los puntos. 3.3.5 Aproximación pseudoespectral: Elementos de referencia El MEE, tal y como lo presentó Patera en [ 127 ], consiste en un método híbrido entre entre el MEF y los MMEE que combina las virtudes de ambos: la flexibilidad y la versatilidad del primero con la precisión espectral del segundo. Al igual que en MEF, el MEE proviene de la aplicación del Método de los Residuos Ponderados (MRP) en el dominio de interés Ω , sin embargo el MEE utiliza una aproximación pseudoespectral dentro de los elementos. El dominio es pues, dividido en un número finito Kde subregiones, de forma que: Ω= K [ j=1 Ωj, K \ j=1 Ωj=/0,(3.61) donde la diferencia fundamental con el MEF clásico se encuentra en la distribución nodal y el grado de las funciones interpolantes en cada subdominio Ωj. En esta tesis se realiza el mismo procedimiento para el MEC en el ámbito de la ecuación mild-slope. Para ello el contorno de interés Γ se divide en un número finito N de subdominios, de forma que: Γ= N [ j=1 Γj, N \ j=1 Γj=/0.(3.62) 3.3 Aproximación mediante elementos espectrales 51 Figura 3.14 Nodos y funciones de forma del elemento de referencia unidimensional ( n=4 ). siendo Γjcada uno de los subcontornos. Para poder llevar a cabo los procesos de interpolación, integración y diferenciación, los subdominios son llevados a unos espacios normalizados. En estas coordenadas los elementos tanto del MEC como los del MEE toman la forma de los elementos de referencia. A continuación se presentan los dos elementos espectrales necesarios para la solución de la ecuación mild-slope, el unidimensional y el bidimensional. Elemento espectral 1-D El elemento unidimensional de referencia está compuesto por un conjunto de nodos {(ξ0,ξ1,...,ξn)ξ∈[−1,1]} y sus respectivas funciones de aproximación asociadas (L0(ξ), L1(ξ),..., Ln(ξ)) . La posición de los nodos se corresponde con los puntos de la cuadratura de LGL, obtenidos mediante la solución del polinomio (3.54) , y sus funciones de aproximación vienen dadas por el conjunto de polinomios de Lagrange: Lk(ξ) = n ∏ j=0,j6=k ξ−ξj ξk−ξj ,(3.63) tal y como se representa en la Figura 3.14 para un elemento de orden n=4. La distribución nodal según la cuadratura LGL proporciona una aproximación cuya precisión aumenta de forma exponencial con el grado de las funciones nodales Lk(ξ) , al igual que ocurre en los métodos espectrales, de ahí su denominación de pseudoespectral, pues demuestran una convergencia similar y su distribución nodal está basada en una familia de polinomios típica de los MMEE. Sin embargo, la aproximación se realiza de manera espacial mediante valores nodales y sus funciones interpolantes. 52 Capítulo 3. Aproximaciones numéricas El elemento unidimensional es requerido en este trabajo para dos propósitos, la aproximación del potencial y el flujo en la expresión discreta de la Ecuación Integral de Contorno (EIC) y la aproximación del flujo en el contorno del dominio modelado con el MEE. Elemento espectral 2-D Para la implementación del MEE es necesaria la definición de un elemento de referencia bidimensional sobre el que realizar de manera simplificada las operaciones necesarias. Podemos definir la región que ocupa el elemento en el espacio normalizado como {(ξ,ζ)|ξ∈[−1,1],ζ∈[−1,1]} , donde la posición de sus nodos en el interior vendrá dada por una distribución LGL en las direcciones ξyζ: X≡(ξk,ζk) = (ξi,ζj),   i=0,1,...,n, j=0,1,...,n, k= (i+1)+ j(n+1), (3.64) Al igual que en (3.56) , una función f(ξ,ζ) puede ser aproximada en el elemento a partir de una serie de polinomios y valores nodales: IInf(ξ,η) = n ∑ k=1 f(ξk,ζk)Lk(ξ,ζ),(3.65) donde IIn denota el operador de interpolación de grado n sobre el conjunto de (n+1)2 puntos internos en el elemento ¯ Ω2 . La función de aproximación correspondiente a cada nodo viene dada como un producto de sus funciones unidimensionales escritas en (3.63) en cada una de las direcciones. El conjunto de todas las funciones de aproximación viene dado entonces por el producto tensorial: Lk(ξ,ζ) = Li(ξ)Lj(ζ),   i=0,1,...,n, j=0,1,...,n, k= (i+1)+ j(n+1), (3.66) En la Figura 3.15(a) se muestra una descripción de la distribución nodal en un elemento de orden n=4 y como los nodos son numerados en su interior. Las subfiguras 3.15(b), 3.15(c) y 3.15(d) muestran las funciones de aproximación correspondientes a los nodos 21, 14 y 22 respectivamente, que son representativos de las diferentes posiciones en las que se puede encontrar un punto: esquina, contorno o interior. 3.3.6 Integración numérica en los elementos espectrales Como se mencionaba en el punto 3.3.3 la cuadratura de LGL es exacta para polinomios de orden hasta 2n−1 , siendo el número de puntos de la cuadratura n+1 . A pesar de ser exacta para polinomios de hasta un orden de 2n−1 y que los polinomios que se están integrando son de orden 2pla precisión se mantiente [42]. 3.3 Aproximación mediante elementos espectrales 53 Figura 3.15 Elemento de referencia bidimensional. (a) elemento y numeración de los nodos, (b) función de forma del nodo 21, (c) función de forma del nodo 14 y (d) función de forma del nodo 22. 3.3.7 Resumen • Si se utilizan los puntos de cuadratura Gaussiana como nodos en la discretización espacial y se les asocia una base de polinomios de Lagrange, se obtiene una base discreta ortogonal, deseable para métodos computacionales. • La interpolación mediante una base de polinomios de Lagrange asociada a los puntos de cuadratura Gaussiana está estrechamente relacionada con la interpolación mediante el espectro de polinomios asociados a la cuadratura. Este tipo de aproximación se denomina pseudo-espectral. • Utlizando una aproximación espectral o pseudo-espectral se consigue convergencia exponencial o espectral. • Los polinomios de Legendre son ortogonales respecto al producto escalar hPm,Pniw con peso w(x) = 1 . Esto hace que distribuciones nodales asociadas a estos polinomios sean especialmente atractivas frente al uso de otras distribuciones nodales, como aquellas basadas en los polinomios de Chebyshev. Esta característica en su ortogonalidad permite que las integrales de las funciones de forma nodales pueden ser llevadas a cabo mediante la propia cuadratura LGL manteniendo su exactitud en la integración hasta polinomios de orden 2n−1. 54 Capítulo 3. Aproximaciones numéricas 3.4 Métodos de acoplamiento 3.4.1 Introducción El acoplamiento superficial entre diferentes dominios tiene un gran interés para la resolución numérica de problemas no convencionales, como por ejemplo aquellos en los que hay un gran número de grados de libertad y se hace necesaria su resolución de manera particionada y computación en paralelo; o también para problemas de contacto [ 149 ], o de interacción fluido estructura en acústica como es el caso de esta tesis. La posibilidad de utilizar además diferentes métodos numéricos para los diferentes dominios, según las características y necesidades de éste, abre la puerta a la obtención de resultados más precisos y con tiempos de computación menores. En el caso de los problemas de acústica lineal, en los que el dominio fluido puede ser infinito, el MEC es más indicado que el MEF, mientras que para la estructura no tiene por qué ser así. Mediante estas técnicas de acoplamiento ambos dominios pueden ser acoplados de manera eficiente y precisa incluso con mallas no conformes. El principal objetivo de esta sección es servir de introducción al Artículo III anejo a esta tesis y aportar algunos detalles más sobre la formulación allí presentada. Este artículo, titulado "Partitioned solution strategies for coupled BEM-FEM acoustic fluid-structure interaction problems", trata específicamente sobre el acoplamiento fluido-estructura mediante el Método Mortar y Multiplicadores de Lagrange Localizados. En este trabajo el dominio fluido acústico es modelado con el MEC, por la gran cantidad de ventajas que presenta para este tipo de problemas (véase la sección 3.2), mientras que para el dominio sólido se utiliza el MEF. En este trabajo, el Método Mortar, basado en multiplicadores clásicos, es comparado con los Multiplicadores de Lagrange Localizados , utilizados aquí por primera vez en este tipo de problemas. Mediante ambos, se resuelven problemas con mallas altamente no conformes y sus sistemas de ecuaciones resultantes son resueltos mediante el método del bi-gradiente conjugado, el cual presenta una buena escalabilidad en este tipo de problemas. En el artículo se han resuelto varios ejemplos de especial relevancia, como son la cavidad acústica propuesta por Sandberg et al. en [ 138 ], donde originalmente se resolvía el conjunto mediante acoplamiento con mallas conformes MEF-MEF. Este problema se estudia con detenimiento en el artículo para varias frecuencias, además, se realiza un barrido para encontrar los modos de vibración naturales del sistema. Otro ejemplo que se resuelve en el artículo es un conducto por el que se transmite una onda de presión. En este ejemplo el dominio fluido es subdivido en varias partes y se estudia como el hecho de particionar el dominio afecta a la convergencia. Para terminar, se estudia una cavidad inmersa en un medio fluido infinito que cuenta con una pared flexible, mostrando, como era de esperar, una convergencia similar al problema de la cavidad interna. 3.4.2 Interacción Fluido-Estructura. Problema variacional La variación del potencial total en el problema de interacción fluido-estructura ( δΠtotal ) puede ser expresada mediante el Principio de los Trabajos Virtuales (PTV): δΠtotal(us,uf,tc) = δΠs(us) +δΠf(uf)+ δΠc(uc s,uc f,tc),(3.67) 3.4 Métodos de acoplamiento 55 donde δΠs , δΠf y δΠc representan la variación del trabajo en el dominio sólido, en el fluido y en la interfase respectivamente. Estos potenciales dependen de los desplazamientos producidos en la estructura (us ), en el fluido ( uf ) y de las fuerzas que surgen en la interfase entre ambos ( tc ). El superíndice c se utiliza para denotar los valores de los desplazamientos en la interfase. La variación del potencial en la estructura viene dada por: δΠs(us) = δuT s{Ks(ω)us−f},(3.68) siendo Ks(ω) la matriz de rigidez dinámica de la estructura en el dominio de la frecuencia. La variación del potencial en el fluido modelado con el MEC puede expresarse como una integral de contorno mediante la fórmula de Clapeyron [22, 62, 66], resultando: δΠf=ZΓf (p−tf)·δun f dΓ,(3.69) donde un f son los desplazamientos normales al contorno del dominio fluido Γf , p es la presión del fluido y tf son tracciones externas normales a la superficie. La ecuación anterior se puede representar de manera discreta como: δΠf=δdT fZΓf NTNdΓ(p−¯ t),(3.70) donde los desplazamientos normales a la superficie, la presión y las tracciones externas son aproximados como un f ≃Ndf,(3.71) p≃Np,(3.72) tf≃Ntf(3.73) respectivamente, siendo N el vector que contiene las funciones de interpolación y d , p y tf los vectores que contienen los valores nodales de los desplazamientos normales a la superficie, de la presión y de las tracciones externas. En la expresión anterior, las presiones pueden ser expresadas en función de los desplazamientos recordando la forma matricial del sistema de ecuaciones del MEC descrito en (3.18). Según la cual: p=iωH−1Gdf,(3.74) donde la relación entre la velocidad y el desplazamiento normal a la superficie vn f =iωun f ha sido utilizada. Aplicando esta expresión en (3.70), resulta: δΠf=δdT fKf(ω)df−M¯ t=0,(3.75) donde Kf(ω) = iωMH−1G , siendo M=RΓfNTNdΓ . Finalmente, de aquí se deduce el sistema de ecuaciones proveniente del MEC que será utilizado más adelante para acoplar 62 Capítulo 3. Aproximaciones numéricas Recordando las expresiones del trabajo virtual producido en la estructura (3.68) y el fluido (3.75) podemos reescribir la expresión anterior como: δΠtotal =δds[K(ω)ds−fs] +δdfKf(ω)df−ff +δλ λ λT sdc s−δλ λ λT sLsdc+δdcT sλ λ λs−δdT cLT sλ λ λs +δλ λ λT fdc f−δλ λ λT fLfdc+δdcT fλ λ λf−δdT cLT fλ λ λf=0, (3.108) que ha de cumplirse para cualquier variación arbitraria (pero compatible con las condiciones de contorno) de los desplazamientos de los dominios, del marco y de las tracciones de la interfase. Por lo tanto, se puede reescribir la ecuación anterior reagrupando sus términos en función del operador al que multiplican como: δds[Ks(ω)ds−fs+Bsλ λ λs] = 0(3.109) δdfKf(ω)df−ffBfλ λ λf=0(3.110) δλ λ λT shBT sds−δλ λ λT sLsdci=0(3.111) δλ λ λT fhBT fdf−δλ λ λT fLfdci=0(3.112) δdT cLT sλ λ λs+LT fλ λ λf=0(3.113) Los términos Bα(α=s,f) representan sendas matrices booleanas que sirven para seleccionar los grados de libertad del dominio α situados en la interfase, es decir, dc α=Bαdα . Obviando la solución trivial, este sistema de ecuaciones nos permite obtener directamente la solución del campo de desplazamientos de ambos dominios, de los desplazamientos del marco y de las tracciones de la interfase. De manera matricial, el sistema se escribe:       Ks(ω)0 Bs0 0 0 Kf(ω)0 Bf0 BT s0 0 0 −Ls 0 BT f0 0 −Lf 0 0 −LT s−LT f0                  ds df λ λ λs λ λ λf dc            =           fs ff 0 0 0 0 0 0 0 0 0            ,(3.114) donde el grupo de ecuaciones (3.109-3.113) está montado de manera monolítica. Para más detalles sobre la resolución de este sistema de forma particionada e iterativa véase el Artículo III. Localización de los nodos del marco de acoplamiento: Regla del momento nulo Hasta ahora no se ha hecho ninguna especificación en cuanto a la malla del marco, sin embargo en el método de los MLL los nodos del marco no pueden estar situados de forma arbitraria, sino que su disposición ha de cumplir los requerimientos de Babuska-Brezi [ 51 , 90 , 171 ]. Para establecer la localización de los nodos del marco de acoplamiento se sigue la Regla del momento nulo, estudiada por Park et al. [ 125 ]. Esta regla consiste en extraer los puntos del marco en los que la resultante de los momentos producidos por los 3.4 Métodos de acoplamiento 63 Ω1: Sólido Marco Ω2: Fluido 1/3 2/3 2/3 1/3 11/2 1/2 Ley de momentos 2L Figura 3.19 Aplicación de la regla del momento nulo. Representación del diagrama de momentos que producen los MLL del fluido y la estructura sobre el marco. Los nodos de la malla del marco han de situarse sobre los puntos de valor nulo del diagrama. MLL de la interfase son nulos. En la Figura 3.19 se representa un ejemplo de una ley de momentos sobre el marco para un conjunto de MLL dados. 3.4.5 Resumen En el Artículo III se han estudiado dos estrategias de acoplamiento MEC-MEF para problemas de interacción fluido-estructura en acústica. Los métodos empleados han sido el Método Mortar y el Método de los Multiplicadores de Lagrange Localizados, de los cuales este último ha sido por primera vez utilizado para este tipo de problemas. En este artículo se propone además el uso del Método del Bi-Gradiente Conjugado para la resolución de los sistemas de ecuaciones resultantes, aportando escalabilidad al problema. Para probar las formulaciones y el método de resolución se han resuelto varios ejemplos, en los que se comprueba que el Método de los Multiplicadores de Lagrange Localizados ofrece mejores resultados que el Método Mortar en General. Por otro lado también se ha comprobado la escalabilidad de los problemas resueltos. 4 Conclusiones y trabajos futuros 4.1 Resumen y conclusiones En esta tesis se han desarrollado diferentes formulaciones numéricas para el estudio de problemas de propagación de ondas lineales en fluidos. Además, se ha abordado el problema acústico de interacción fluido-estructura, modelando el medio fluido y el sólido mediante diferentes métodos numéricos con mallas no conformes en la interfase. El trabajo llevado a cabo durante la realización de esta tesis ha dado lugar a los tres artículos anejos y de los que se hace una revisión a continuación; destacando sus resultados y conclusiones más relevantes. •Artículo I En el artículo contenido en el Apéndice A, se desarrolla una solución fundamental para el problema de Helmholtz en dos dimensiones con el número de onda variable en una dirección. El trabajo se ha desarrollado en el contexto de la propagación de oleaje modelado mediante la EPS, que reproduce el fenómeno de propagación de olas de pequeña amplitud, incluyendo los efectos de difracción, refracción y asomeramiento. Dicha solución fundamental, que está basada en la inversión mediante el algoritmo IFFT de un problema unidimensional en el espacio de Fourier, es además introducida en una formulación MEC que posteriormente es utilizada para resolver algunos ejemplos de referencia. Entre esos ejemplos se encuentran dos especialmente destacables. El primero de ellos consiste en un problema de propagación de olaje sobre una rampa en la que se superpone un montículo con forma elíptica. Esta variación en la geometría del fondo provoca una perturbación en el campo de la altura de ola. El problema ha sido resuelto mediante un acoplamiento MEC-MEF, utilizando el MEF para modelar la región interna con batimetría irregular y el MEC para la región externa con batimetría variable en una dirección. Este problema es especialmente relevante pues ha servido como prueba para numerosas formulaciones sobre la propagación de oleaje desde los años 70. 65 66 Capítulo 4. Conclusiones y trabajos futuros El último problema resuelto el artículo es el modelo de propagación de olaje en el puerto de Chipiona (Cádiz). Este caso ha sido tratado exclusivamente mediante el MEC, aproximando la forma real de la batimetría mediante líneas rectas y paralelas a la línea de costa. El modelo permite hallar la solución del problema cubriendo un gran área alrededor del puerto con tan sólo discretizar el contorno de este. Los resultados contenidos en el Artículo I muestran que el MEC provisto de la solución fundamental para profundidad variable es una herramienta muy útil, pues en muchos problemas reales la batimetría puede ser aproximada mediante líneas rectas y paralelas a la linea de costa. Además, la solución fundamental puede ser incluida dentro de un modelo acoplado MEC-MEF, donde el MEF es usado para modelar la región interna y el MEC la externa. Este acoplamiento permite resolver problemas con batimetrías complejas, cumpliéndose la condición de radiación en el infinito para condiciones de contorno de dominio abierto incluso con batimetría variable. •Artículo II El artículo contenido en el Apéndice B presenta una formulación de acoplamiento MEC-MEF para la EPS mediante elementos espectrales. El tipo de elemento utilizado aproxima las variables de forma espacial mediante polinomios de Lagrange con posiciones nodales localizadas en los puntos de cuadratura LGL. Esta clase de aproximación es usualmente denominada como pseudo-espectral y confiere al MEC y al MEF una convergencia exponencial de la solución con el incremento del orden de los elementos. Las principal aportación de este artículo es la formulación del Método de los Elementos de Contorno Espectrales (MECE) para el problema de Helmholtz con número de onda constante y variable. Dentro de esta formulación, se detalla también un método semi-analítico para la evaluación de las integrales impropias inherentes al MEC. Además, se lleva a cabo el acoplamiento MECE-MEE para la resolución de problemas de propagación de oleaje en dominios abiertos con fondo variable, combinando las ventajas de ambos métodos. La convergencia del MECE ha sido comprobada y comparada con la del MEE en un problema de propagación sobre fondo plano de solución conocida. El ejemplo se resuelve incrementando el orden de aproximación elemental, demostrando que el MECE presenta una convergencia similar a la del MEE. El acoplamiento MECE-MEE ha sido examinado mediante la solución de dos problemas de referencia para modelos basados en la EPS. El primero de ellos consiste en un problema de propagación de oleaje sobre una región con fondo plano en la que existe un montículo con forma parabólica. En el segundo ejemplo, el oleaje se propaga sobre una rampa en la que hay un montículo con forma elíptica; para hallar la solución de este último ejemplo es necesario dotar al MECE de la solución fundamental presentada en el Artículo I . En ambos casos se hace uso de la versatilidad del acoplamiento MECE-MEE, que permite modelar regiones de batimetría variable en dominios abiertos, satisfaciendo la condición de radiación en el infinito. El método demuestra tener una convergencia exponencial, reduciendo en 4.1 Resumen y conclusiones 67 gran medida el número de grados de libertad requeridos comparado con un método de acoplamiento MEC-MEF clásico. •Artículo III El tercer artículo que forma parte de esta tesis, incluido en el Apéndice C, estudia diferentes estrategias de acoplamiento MEC-MEF para problemas de interacción fluido-estructura en acústica. En este trabajo, el MEC es utilizado para modelar el dominio fluido aprovechando sus ventajas en cuanto a precisión, número de grados de libertad requeridos y cumplimiento de la condición de radiación de Sommerfeld en el infinito; mientras que el MEF es utilizado para modelar la estructura. Los métodos de acoplamiento utilizados son el método Mortar y el método de los MLL, siendo este último utilizado por primera vez para el problema acústico. Ambas formulaciones permiten llevar a cabo acomplamientos con mallas no conformes en la interfase. Además, los sistemas de ecuaciones lineales resultantes de la aplicación de estas formulaciones particionadas son resueltos mediante el método iterativo BiCGSTAB, demostrando la escalabilidad del proceso de solución del sistema particionado. En el artículo se resuelven tres ejemplos en los que se analiza el comportamiento del método Mortar y del método de los MLL. En el primero de ellos se analizan los modos de vibración de una cavidad acústica mediante un barrido en frecuencia, presentándose además, los desplazamientos y presiones en la estructura, el fluido y la interfase de forma detallada. Este ejemplo ha sido resuelto para mallas altamente no conformes entre el fluido y la estructura. El segundo ejemplo resuelve un problema en el que una onda de presión viaja a través de un conducto. Este problema es resuelto aumentando el número particiones, comprobándose que ante este hecho el número de iteraciones requerido para convergencia no se ve afectado de manera significativa. Finalmente, se resuelve un problema externo en el que un tren de ondas incide sobre una cavidad hueca con una pared flexible. En todos los ejemplos estudiados los sistemas resultan ser escalables, dependiendo el número de iteraciones necesarias para el BiCGSTAB principalmente del número de grados de libertad en la interfase y el grado de no conformidad existente entre ambas mallas. En conclusión, en los tres artículos que componen esta tesis se han abordado problemas de transmisión de ondas lineales en fluidos, aportando cada uno de ellos soluciones a tres problemas de actual interés en ingeniería. En el primero de ellos se propone una solución fundamental que permite resolver problemas de propagación de ondas en medios no homogéneos mediante el MEC, en el segundo se aborda la dispersión numérica y el error de polución a través del uso de elementos espectrales de alto orden y, por último, en el tercer artículo se estudian estrategias de acoplamiento para problemas de interacción fluido-estructura en acústica. 68 Capítulo 4. Conclusiones y trabajos futuros 4.2 Trabajos futuros Como resultado de esta tesis quedan abiertas muchas e interesantes líneas de investigación para el futuro. A partir del primero de los artículos anejos podría continuarse con: • La extensión de la formulación MEC a modelos de propagación de oleaje más complejos, como la Ecuación de Pendiente Suave Modificada (EPSM) [ 35 ] o la EPSE [ 112 ], lo que permitiría ampliar el espectro de problemas reales que pueden ser resueltos mediante el MEC. • Por otro lado, se podría desarrollar una solución fundamental con el número de onda variable en una dirección para el problema tridimensional siguiendo un procedimiento parecido. Esta solución fundamental permitiría resolver el problema acústico en medios no homogéneos mediante el MEC. A partir el segundo artículo, en el que se conjugan el MEF, el MEC y los elementos espectrales, surgen las siguientes ideas que pueden ser llevadas a cabo en el futuro: • La formulación mediante el MECE del problema acústico tridimensional. Esta técnica aportaría una reducción significativa del coste computacional, además de incrementar la precisión, especialmente para altas frecuencias. • La aplicación del MEE y el MECE en teorías de propagación de oleaje no lineales en el dominio del tiempo, como el modelo de Boussinesq. Discretizaciones basadas en elementos espectrales están empezando a incorporarse a estas formulaciones, mostrando muy buenos resultados. Por último, como consecuencia del tercer artículo y continuando las propuestas anteriores, las siguientes líneas de trabajo podrían ser abordadas en el futuro: • La extensión al espacio tridimensional de las formulaciones MEC-MEF acopladas mediante el método Mortar y el método de los MLL en problemas de interacción fluido estructura en acústica. Sin duda esta línea es de gran interés, pues está ligada numerosas aplicaciones en la industria y la construcción. Tras los buenos resultados demostrados en dominios bidimensionales, es de esperar que el comportamiento de estas formulaciones sea parecido para el caso tridimensional. • Otro aspecto interesante sería llevar a cabo a cabo una formulación MECE-MEE para problemas de interacción fluido-estructura mediante el Método Mortar o el método de los MLL. A partir de mallas conformes en la interfase el refinamiento en cualquiera de ellas podría llevarse a cabo aumentando el orden del elemento. Bibliografía [1] M. Ainsworth. Discrete dispersion relation for hp-version finite element approximation at high wave number. SIAM Journal on Numerical Analysis, 42(2):553–575, 2004. [2] G. B. Airy. Encyclopedia Metropolitana, chapter Tides and waves, page 289. 1845. [3] I. Babu ˘ ska and M. Suri. The p- and h-p versions of the finite element method, an overview. Computer Methods in Applied Mechanics and Engineering, 80(1):5 – 26, 1990. [4] I. Babu ˘ ska, B. A. Szabo, and I. N. Katz. The p-version of the finite element method. SIAM Journal on Numerical Analysis, 18(3):515–545, 1981. [5] I. M. Babu ˘ ska and S. A. Sauter. Is the pollution effect of the fem avoidable for the helmholtz equation considering high wave numbers? SIAM Journal on Numerical Analysis, 34(6):2392–2423, 1997. [6] M. R. Bai and S. Chang. Active noise control of enclosed harmonic fields by using bem-based optimization techniques. Applied Acoustics, 48(1):15 – 32, 1996. [7] P. K. Banerjee and T. G. Davies. The behaviour of axially and laterally loaded single piles embedded in nonhomogeneous soils. Géotechnique, 28(3):309–326, 1978. [8] A. Bayliss, C. Goldstein, and E. Turkel. On accuracy conditions for the numerical computation of waves. Journal of Computational Physics, 59(3):396 – 404, 1985. [9] K. Belibassakis. The Green’s function of the mild-slope equation: The case of a monotonic bed profile. Wave Motion, 32:339–361, 2000. [10] K. Belibassakis, G. Athanassoulis, and T. Gerostathis. A coupled-mode model for the refraction-diffraction of linear waves over steep three-dimensional bathymetry. Applied Ocean Research, 23(6):319 – 336, 2001. 69 70 Bibliografía [11] K. Belibassakis, T. Gerostathis, and G. Athanassoulis. A coupled-mode model for water wave scattering by horizontal, non-homogeneous current in general bottom topography. Applied Ocean Research, 33(4):384 – 397, 2011. [12] G. Bellotti, G. M. Beltrami, and P. D. Girolamo. Internal generation of waves in 2d fully elliptic mild-slope equation FEM models. Coastal Engineering, 49:71 – 81, 2003. [13] G. M. Beltrami, G. Bellotti, P. D. Girolamo, and P. Sammarco. Treatment of wave breaking and total absorption in a mild-slope equation fem model. Journal of Waterway, Port, Coastal, and Ocean Engineering, 127(5):263–271, 2001. [14] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2):185 – 200, 1994. [15] P. G. Bergmann. The equation in a medium with a variable index of refraction. The Journal of the Acoustical Society of America, 17(4):329 – 333, 1946. [16] J. Berkhoff. Finite Element in Fluids, volume 1, chapter Linear wave propagation problem and the finite element method, pages 251–264. Wiley New York, 1975. [17] J. Berkhoff, N. Booy, and A. Radder. Verification of numerical wave propagation models for simple harmonic linear water waves. Coastal Engineering, 6(3):255 – 279, 1982. [18] W. Berkhoff. Computation of combined refraction-diffraction. In Proceedings of 13th International Conference on Coastal Engineering. ASCE, 1972. [19] W. Berkhoff. Mathematical Models for Simple Harmonic Linear Water Waves. Wave Diffraction and Refraction. PhD thesis, Delft Hydraulics Laboratory, 1976. [20] P. Bettess and O. C. Zienkiewicz. Diffraction and refraction of surface waves using finite and infinite elements. International Journal for Numerical Methods in Engineering, 11(8):1271–1290, 1977. [21] D. Blackstock. Fundamentals of Physical Acoustics. A Wiley-Interscience publication. Wiley, 2000. [22] M. Bonnet. Encyclopedia of Computational Mechanics: Boundary integral equation methods for elastic and platic methods, volume 2. John Wile and Sons, 2004. [23] N. Booij. Gravity Waves on Water with Non-uniform Depth and Current. PhD thesis, 1981. [24] N. Booij. A note on the accuracy of the mild-slope equation. Coastal Engineering, 7(3):191 – 203, 1983. [25] P. Bouillard. Influence of the pollution on the admissible field error estimation for fe solutions of the helmholtz equation. International Journal for Numerical Methods in Engineering, 45(7):783–800, 1999. Bibliografía 71 [26] J. V. Boussinesq. Thèorie des ondes et des remous qui se propagent le long dún canal rectangulaire horizontal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. Journal de Mathèmatiques Pures et Appliqueès, 2(17):55 – 108, 1872. [27] A. Brancati, M. H. Aliabadi, and V. Mallardo. A bem sensitivity formulation for three-dimensional active noise control. International Journal for Numerical Methods in Engineering, 90(9):1183–1206, 2012. [28] C. Brebbia and S. Walker. Boundary Element Techniques in Engineering. Newnes- Butterworths, 1980. [29] C. A. Brebbia and J. Domínguez. Boundary elements: an introductory course. Computational Mechanics, 2 edition, 1992. [30] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Z. Jr. Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer-Verlag Berlin Heidelberg, 1 edition, 1988. [31] C. G. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral methods: fundamentals in single domains. Scientific Computation. Springer, 2nd printing. 3rd printing. edition, 2010. [32] A. Cerrato, J. A. González, and L. Rodríguez-Tembleque. Boundary element formulation of the mild-slope equation for harmonic water waves propagating over unidirectional variable bathymetries. Engineering Analysis with Boundary Elements, 62:22 – 34, 2016. [33] O. Cessenat and B. Despres. Application of an ultra weak variational formulation of elliptic pdes to the two-dimensional helmholtz problem. SIAM Journal on Numerical Analysis, 35(1):255–299, 1998. [34] C. C.G., H. M.Y., and Z. T.A. Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. 2007. [35] P. Chamberlain and D. Porter. The modified mild-slope equation. Journal Fluid Mechanics, 291:393 – 407, 1995. [36] C. N. Chandrasekera and K. F. Cheung. Extended linear refraction-diffraction model. Journal of Waterway, Port, Coastal, and Ocean Engineering, 123(5):280 – 286, 1997. [37] H. Chen. Effects of bottom friction and boundary absorption on water wave scattering. Applied Ocean Research, 8(2):99 – 104, 1986. [38] H. Chen and C. Mei. Oscillations and wave forces in a man-made harbor in the open sea. In Symposium on Naval Hydrodynamics 10th. Cambridge, June 1974. 78 Bibliografía [123] K. C. Park, C. A. Felippa, and U. Gumaste. A localized version of the method of lagrange multipliers and its applications. Computational Mechanics, 24:476–490, 2000. [124] K. C. Park, C. A. Felippa, and R. Ohayon. Partitioned formulation of internal fluid-structure interaction problems by localized lagrange multipliers. Computer Methods in Applied Mechanics and Engineering, 190(24-25):2989 – 3007, 2001. Advances in Computational Methods for Fluid-Structure Interaction. [125] K. C. Park, C. A. Felippa, and G. Rebel. A simple algorithm for localized construction of non-matching structural interfaces. International Journal for Numerical Methods in Engineering, 53:2117–2142, 2002. [126] K. C. Park, R. Ohayon, C. A. Felippa, and J. A. González. Partitioned formulation of internal and gravity waves interacting with flexible structures. Computer Methods in Applied Mechanics and Engineering, 199(9-12):723 – 733, 2010. [127] A. T. Patera. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. Journal of Computational Physics, 54:468 – 488, 1984. [128] S. Petersen, D. Dreyer, and O. von Estorff. Assessment of finite and spectral element shape functions for efficient iterative simulations of interior acoustics. Computer Methods in Applied Mechanics and Engineering, 195(44â € “47):6463 – 6478, 2006. [129] A. Pierce. Acoustics: An Introduction to Its Physical Principles and Applications. Acoustical Society of America, 1989. [130] D. Porter and D. Staziker. Extensions of the mild-slope equation. Journal of Fluid Mechanics, 300:367–382, 1995. [131] A. C. Radder. On the parabolic equation method for water-wave propagation. Journal of Fluid Mechanics, 95:179–176, 1979. [132] C. Rajakumar and A. Ali. Boundary element-finite element coupled eigenalalysis of fluid-structure systems. International Journal for Numerical Methods in Engineering, 39(10):1625–1634, 1996. [133] G. Rebel, K. C. Park, and C. A. Felippa. A contact formulation based on localized lagrange multipliers: formulation and application to two-dimensional problems. International Journal for Numerical Methods in Engineering, 54:263–297, 2002. [134] F. J. Rizzo. An integral equation approach to boundary value problems of classical esatostatics. Quarterly of Applied Mathematics, 25(1):83–95, 1967. [135] Z. Rong and C. Xu. Numerical approximation of acoustic waves by spectral element methods. Appl. Numer. Math., 58(7):999–1016, July 2008. [136] M. R. Ross, C. A. Felippa, K. C. Park, and M. A. Sprague. Treatment of acoustic fluid-structure interaction by localized lagrange multipliers: Formulation. Computer Methods in Applied Mechanics and Engineering, 197(33-40):3057 – 3079, 2008. Bibliografía 79 [137] M. R. Ross, M. A. Sprague, C. A. Felippa, and K. C. Park. Treatment of acoustic fluid-structure interaction by localized lagrange multipliers and comparison to alternative interface-coupling methods. Computer Methods in Applied Mechanics and Engineering, 198(9-12):986 – 1005, 2009. [138] G. Sandberg, P. Wernberg, and P. Davidson. Fundamentals of fluid-structure interaction. In G. Sandberg and R. Ohayon, editors, G. Sandberg and R. Ohayon (Eds.) Computational Aspects of structural acoustics and vibratinon, pages 23–102, New York, 2008. Springer Wien New York. [139] S. Schneider. Fe/fmbe coupling to model fluid-structure interaction. International Journal for Numerical Methods in Engineering, 76(13):2137–2156, 2008. [140] J. C. Schönfeld. Propagation of two dimensional short waves. Delft University of Technology Manuscript (in Dutch), 1972. [141] R. Seznec. Diffraction of sound around barriers: Use of the boundary elements technique. Journal of Sound and Vibration, 73(2):195 – 209, 1980. [142] R. Shaw and W. Falby. FEBIE - A combined finite element-boundary integral equation method. Computers & Fluids, 6(3):153 – 160, 1978. [143] J. Shen, T. Tang, and L.-L. Wang. Spectral Methods: Algorithms, Analysis and Applications. Springer Series in Computational Mathematics 41. Springer-Verlag Berlin Heidelberg, 1 edition, 2011. [144] J. Shen and L.-L. Wang. Spectral approximation of the helmholtz equation with high wave numbers. SIAM Journal on Numerical Analysis, 43(2):623–644, 2005. [145] B. Shizgal. Spectral Methods in Chemistry and Physics: Applications to Kinetic Theory and Quantum Mechanics. Scientific Computation. Springer Netherlands, 1 edition, 2015. [146] J. Simo, P. Wriggers, and R. Taylor. A perturbed lagrangian formulation for the finite element solution of contact problems. Computer Methods in Applied Mechanics and Engineering, 50:163–180, 1985. [147] E. Skudrzyk. The Foundations of Acoustics: Basic Mathematics and Basic Acoustics. Springer-Verlag Wien, 1 edition, 1971. [148] R. Smith and T. Sprinks. Scattering of surface waves by a conical island. Journal of Fluid Mechanics, 72:373–384, 1975. [149] L. R.-T. Solano. Formulación numérica de la interacción mecánica entre superficies de sólidos 3D. PhD thesis, Universidad de Sevilla, Sevilla, mayo 2009. [150] D. Steward and V. Panchang. Improved coastal boundary condition for surface water waves. Ocean Engineering, 28(1):139 – 157, 2001. 80 Bibliografía [151] G. G. Stokes. On the theory of oscillatory waves. In Transactions of the Cambridge Philosophical Society, volume 8, pages 441–455. 1847. [152] K. D. Suh, C. Lee, and W. S. Park. Time-dependent equations for wave propagation on rapidly varying topography. Coastal Engineering, 32:91 – 117, 1997. [153] L. L. Thompson and P. M. Pinsky. A galerkin least-squares finite element method for the two-dimensional helmholtz equation. International Journal for Numerical Methods in Engineering, 38(3):371–397, 1995. [154] T.-K. Tsay andP. L.-F. Liu. A finite element model for wave refraction and diffraction. Applied Ocean Research, 5(1):30 – 37, 1983. [155] T.-K. Tsay, W. Zhu, and P.-F. Liu. A finite element model for wave refraction, diffraction, reflection and dissipation. Applied Ocean Research, 11(1):33 – 38, 1989. [156] E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov. Compact 2d and 3d sixth order schemes for the helmholtz equation with variable wave number. Journal of Computational Physics, 232(1):272 – 287, 2013. [157] M. J. Turner, R. W. Clough, H. C. Martin, and L. J. Topp. Stiffness and deflection analysis of complex structures. Journal of the Aeronautical Sciences, 23:805–823, 1956. [158] G. van de Sande. Literature study for the msc thesis. the mild-slope equation and its numerical implementation. Master’s thesis, Delft University of Technology, 12 2011. [159] P. E. Vos, S. J. Sherwin, and R. M. Kirby. From h to p efficiently: Implementing finite and spectral/hp element methods to achieve optimal performance for low- and high-order discretisations. Journal of Computational Physics, 229(13):5161 – 5181, 2010. [160] D. Wang, R. Tezaur, J. Toivanen, and C. Farhat. Overview of the discontinuous enrichment method, the ultra-weak variational formulation, and the partition of unity method for acoustic scattering in the medium frequency regime and performance comparisons. International Journal for Numerical Methods in Engineering, 89(4):403–417, 2012. [161] H. Wang and S. Xiang. On the convergence rates of legendre approximation. Mathematics of Computation, 81(278):861–877, 2012. [162] K. Wang and Y. S. Wong. Pollution-free finite difference schemes for nonhomogeneous Helmholtz equation. International Journal of Numerical Analysis and Modeling, 4(11):787–815, December 2014. [163] K. Watanabe. Integral Transform Techniques for Green’s Function. Lecture Notes in Applied and Computational Mechanics 76. Springer International Publishing, 2 edition, 2015. Bibliografía 81 [164] D. T. Wilton. Acoustic radiation and scattering from elastic structures. International Journal for Numerical Methods in Engineering, 13(1):123–138, 1978. [165] B. I. Wohlmuth. Discretization Methods and Iterative Solvers Based on Domain Decomposition. Lecture Notes in Computational Science and Engineering 17. Springer-Verlag Berlin Heidelberg, 1 edition, 2001. [166] T. Wu. Boundary Element Acoustics, Fundamentals and Computer Codes, volume 3. WIT Press, Boston, 5 edition, 2000. [167] Y. Xu. Radiation condition and scattering problem for time-harmonic acoustic waves in a stratified medium with a nonstratified inhomogeneity. IMA Journal of Applied Mathematics, 54(1):9–29, 1995. [168] L. Yao, Y. Li, and L. Li. Dispersion error reduction for acoustic problems using the smoothed finite element method (sfem). International Journal for Numerical Methods in Fluids, 80(6):343–357, 2016. FLD-14-0421.R1. [169] S. Zhu. A new DRBEM model for wave refraction and diffraction. Engineering Analysis with Boundary Elements, 12(4):261–274, 1993. [170] S.-P. Zhu, H.-W. Liu, and K. Chen. A general DRBEM model for wave refraction and diffraction. Engineering Analysis with Boundary Elements, 24(5):377 – 390, 2000. [171] O. Zienkiewicz, S. Qu, R. Taylor, and S. Nakazawa. Finite element modeling of rolling contact. International Journal for Numerical Methods in Engineering, 23:1873–1883, 1986. Apéndice A Artículo I: Boundary element formulation of the Mild-Slope Equation for harmonic water waves propagating over unidirectional variable bathymetries The original version of this paper can be found in www.sciencedirect.com •DOI: 10.1016/j.enganabound.2015.09.006 •Journal Name: Engineering Analysis with Boundary Elements •ISSN: 0955-7997 •Journal Citation Reports (2015): Impact factor: 1.862 –Classification (Engineering, Multidisciplinary): Q1(20/85) –Classification (Mathematics, Interdisciplinary Applications): Q1(23/101) •SCIMAGO (2015): Impact factor: 1.251 –Classification (Applied Mathematics): Q1(86/461) –Classification (Computational Mathematics): Q1(22/134) –Classification (Engineering (miscellaneous)): Q1(26/544) –Classification (Analysis): Q2(34/118) 83 85 ANTONIO CERRATO CASADO BOUNDARY ELEMENT FORMULATION OF THE MILD-SLOPE EQUATION FOR HARMONIC WATER WAVES PROPAGATING OVER UNIDIRECTIONAL VARIABLE BATHYMETRIES ENGINEERING ANALYSIS WITH BOUNDARY ELEMENTS - 62 FEATURED IN THE SCIENCEDIRECT TOP 25 LIST OF MOST DOWNLOADED ARTICLES RANKED 4 TH ON THE TOP 25 FOR ENGINEERING ANALYSIS WITH BOUNDARY ELEMENTS - OCTOBER TO DECEMBER 2015 A.1 Introduction 87 Boundary element formulation of the Mild-Slope Equation for harmonic water waves propagating over unidirectional variable bathymetries Authors: Antonio Cerrato, José A. González and Luis Rodríguez-Tembleque Escuela Ténica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain Abstract This paper presents a boundary element formulation for the solution of the Mild-Slope equation in wave propagation problems with variable water depth in one direction. Based on the Green’s function approximation proposed by Belibassakis [ 3 ], a complete fundamental-solution kernel is developed and combined with a boundary element scheme for the solution of water wave propagation problems in closed and open domains where the bathymetry changes arbitrarily and smoothly in a preferential direction. The ability of the proposed formulation to accurately represent wave phenomena like refraction, reflection, diffraction and shoaling, is demonstrated with the solution of some example problems, in which arbitrary geometries and variable seabed profiles with slopes up to 1:3 are considered. The obtained results are also compared with theoretical solutions, showing an excellent agreement that demonstrates its potential. Keywords: Wave propagation; Mild-slope equation; Helmholtz equation; Boundary element method A.1 Introduction Wave propagation in variable-depth waters is a problem of significant importance in coastal engineering with applications in the design and maintenance of harbors, coastal defense works and hydrodynamic and sediment transportation studies. It is well known that the transmission of linear waves in intermediate and deep waters can be reproduced by the elliptical Mild-Slope Equation (MSE) which was derived by Berkhoff [ 9 ] in the early 70s. The MSE considers simultaneously the effects of diffraction, refraction, reflection and shoaling of linear water surface waves and it is formally valid for slowly varying sea bed slopes, i.e. ∇h<< kh , being h the water depth and k the wave number. The validity of the MSE has been evaluated by Tsay and Liu [ 51 ] demonstrating that it produces accurate results for bottom slopes up to 1:1 when waves are propagating perpendicularly to the bathymetry contour lines. Nevertheless, Booij [ 13 ] verified that, for general directions of wave propagation, the MSE is able to provide acceptable accuracy for bottom profiles with slopes up to 1:3, enough for practical applications. 94 Artículo I ξ2 ξ1 ˆ k1ˆ k3 −ˆ k1 −ˆ k3 ξ2=−τC C+ 1 C+ 2C+ 3 C− 1 C− 2 C− 3 ξ2=τ −Ξ Ξ Figure A.4 Approximation and decomposition into linear paths of the integration path in the complex ξ -plane used to compute the inverse Fourier transform of the function Ψ. velocity potential in the following way: ψ(r,ro;ˆ k) = 1 πZ−τ 0 Ψ(x,xo;iξ2)cosh(yξ2)dξ2+ 1 2πeτyZΞ −Ξ Ψ(x,xo;|ξ1|−iτ)eiyξ1dξ1+ 1 πsinh(τy)ZΞ 0 Ψ(x,xo;ξ1−iτ)e−iyξ1dξ1+ 1 πZ∞ Ξ Ψ(x,xo;ξ1−iτ)cos((ξ1−iτ)y)dξ1,(A.23) where the first integral, along C+ 1 , can be disregarded if τ is considered sufficiently small. The second and the third integrals correspond to the integration along the second path C+ 2 , with the third one containing an hyperbolic sine of small argument that can also be neglected without an important loss of accuracy. The last integral corresponds to the third path C+ 3 , where an asymptotic behavior of the integrand, defined by (A.12) , can be assumed. Under these assumptions, the approximation of the velocity potential can be finally reduced to: ψ(r,ro;ˆ k)≃1 2πeτyZΞ −Ξ Ψ(x,xo;|ξ1|−iτ)eiyξ1dξ1+ 1 2πcosh(τy)ℜ{E1((|x−xo|+iy)Ξ)},(A.24) expression previously proposed by Belibassakis in [ 3 ] for the evaluation of the Green’s function. However, to complete the fundamental solution kernel, we also need to compute the spatial derivatives. A.3 Fundamental solution for variable wave number 95 The x -derivative of the velocity potential is obtained using the same procedure. Substituting the integration path C+in (A.21), we have: ψ,x(r,ro;ˆ k) = 1 πZ−τ 0 Ψ,x(x,xo;iξ2)cosh(yξ2)dξ2 +1 2πeτyZΞ −Ξ Ψ,x(x,xo;|ξ1|−iτ)eiyξ1dξ1 +1 πsinh(τy)ZΞ 0 Ψ,x(x,xo;ξ1−iτ)e−iyξ1dξ1 +1 πZ∞ Ξ Ψ,x(x,xo;ξ1−iτ)cos((ξ1−iτ)y)dξ1,(A.25) where, again, the first and the third integrals can be neglected for a small value of τ . Using the asymptotic behavior of the velocity potential derivative, the integral along the semi-infinite interval C+ 3can be evaluated analytically: Z∞ Ξ Ψ,x(x,xo;ξ1−iτ)cos((ξ1−iτ)y)dξ1= 1 2 e−|x−xo|(Ξ−iτ) (x−xo)2+y2[ysin((Ξ−iτ)y)−|x−xo|cos((Ξ−iτ)y)],(A.26) an substituting back in (A.25) we arrive to the final approximation for the x -derivative of the velocity potential: ψ,x(r,ro;ˆ k)≃1 2πeτyZΞ −Ξ Ψ,x(x,xo;|ξ1|−iτ)eiyξ1dξ1+ +1 2π e−|x−xo|(Ξ−iτ) (x−xo)2+y2[ysin((Ξ−iτ)y)−|x−xo|cos((Ξ−iτ)y)].(A.27) However, obtaining the y -derivative is more involved. We see from its definition (A.22) that the integrand is antisymmetric, so integration in C+ can be carried out considering only the antisymmetric part of the exponential complex function as follows: ψ,y(r,ro;ˆ k) = 1 πZ−τ 0 ξ2Ψ(x,xo;ξ2)sinh(yξ2)dξ2+ 1 2πeτyZΞ −Ξ i(|ξ1|−iτ)Ψ(x,xo;|ξ1|−iτ)eiyξ1dξ1+ 1 πcosh(τy)ZΞ 0 i(ξ1−iτ)Ψ(x,xo;ξ1−iτ)e−iyξ1dξ1+ 1 πZ∞ Ξ i(ξ1−iτ)Ψ(x,xo;ξ1−iτ)sin((ξ1−iτ)y)dξ1,(A.28) where, once more, the first integral can be neglected, but the hyperbolic cosine that appears now in the third integral should be retained. To facilitate its numerical computation, this 96 Artículo I third integral will be extended to a symmetrical integration interval using the symmetry properties of the integrand in the following way: ZΞ 0 i(ξ1−iτ)Ψ(x,xo;ξ1−iτ)e−iyξ1dξ1= 1 2ZΞ −Ξ i(|ξ1|−iτ)Ψ(x,xo;|ξ1|−iτ)e−iyξ1dξ1+ 1 2ZΞ −Ξ iξΨ(x,xo;|ξ1|−iτ)e−iyξ1dξ1,(A.29) and the fourth integral, along the path C+ 3 , is evaluated analytically using the asymptotic value of the transformed velocity potential: Z∞ Ξ i(ξ1−iτ)Ψ(x,xo;ξ1−iτ)sin((ξ1−iτ)y)dξ1= 1 2 e−|x−xo|(Ξ−iτ) (x−xo)2+y2[|x−xo|sin((Ξ−iτ)y)+ycos((Ξ−iτ)y)],(A.30) results that are substituted back in (A.28) to find the final approximation for the y-derivative of the velocity potential: ψ,y(r,ro;ˆ k)≃1 2πeτyZΞ −Ξ i(|ξ1|−iτ)Ψ(x,xo;|ξ1|−iτ)eiyξ1dξ1+ 1 2πcosh(τy)ZΞ −Ξ i(|ξ1|−iτ)Ψ(x,xo;|ξ1|−iτ)e−iyξ1dξ1+ 1 2πcosh(τy)ZΞ −Ξ iξΨ(x,xo;|ξ1|−iτ)e−iyξ1dξ1+ 1 2π e−|x−xo|(Ξ−iτ) (x−xo)2+y2[|x−xo|sin((Ξ−iτ)y)+ycos((Ξ−iτ)y)],(A.31) an expression where all the integration paths are now symmetric. As we will see, this symmetry is needed for an efficient numerical evaluation using FFT. In conclusion, equations (A.24) , (A.27) and (A.31) constitute the complete kernel of the fundamental solution back-transformed to the space domain. Next step is to devise an efficient and accurate numerical procedure to evaluate these integrals. Numerical evaluation of the Fourier integrals The final expression of the velocity potential (A.24) and its derivatives, (A.27) and (A.31) , can be calculated very efficiently by means of the FFT algorithm, as proposed in [3]. Starting with the velocity potential (A.24), it can be expressed in a compact form as: ψ(r,ro;ˆ k)≃ψ(1)+ψ(2) an (A.32) where ψ(1) is the first definite integral and ψ(2) an represents the second analytic term on the A.3 Fundamental solution for variable wave number 97 right hand side. To compute ψ(1) , the domain of integration is discretized using a uniform mesh of N elements, where ξl= (l−1)∆ξ are the nodal locations of the l=1,···,N+1 sampling points and ∆ξ=Ξ/N is the element length; in addition, the physical space y∈[0,Y] is discretized with a similar uniform distribution of nodes yj= (j−1)∆y , for j=1,···,N+1 , separated a distance ∆y=Y/N using a fixed value of ∆y=π/Ξ . The symmetry property expressed in (A.10) assures that the discrete values of the transformed velocity potential Ψl=Ψ(x,xo;ξl−iτ) , for l=1,...,N+1 , are symmetric with respect to ξ1 and hence Ψl=Ψ2N−l+2 , for l=2,...,N . Based on these discretizations and applying the IFFT algorithm, the value of ψ(1)can then be approximated by the finite series: ψ(1)=1 2πeτyj"M ∑ l=1 Ψlei2π M(j−1)(l−1)#(A.33) where the total number of points M=2N is selected as a power of two to be efficiently evaluated. In a similar way, the x -derivative given by (A.27) can be written as the addition of two terms: ψ,x(r,ro;ˆ k)≃ψ(1) ,x+ψ(2) ,x an (A.34) with a definite integral ψ(1) ,x that needs to be evaluated at the same interval. Noting that the discrete values of the velocity potential derivative Ψ,x l =Ψ,x(x,xo;ξl−iτ) , for l=1,...,N+1 , are symmetric with respect to ξ1 and using the same discretization, we can sample Ψ,x l =Ψ,x(x,xo;ξl−iτ) , for l=1,...,N+1 , and apply the symmetry property Ψ,x l =Ψ,x2N−l+2, for l=2,...,N, to approximate Ψ(1) ,xas the IFFT sequence: ψ(1) ,x=1 2πeτyj"M ∑ l=1 Ψ,x l ei2π M(j−1)(l−1)#.(A.35) Finally, according to equation (A.31) , the y -derivative of the velocity potential can be decomposed into four different terms: ψ,y(r,ro;ˆ k)≃ψ(1) ,y+ψ(2) ,y+ψ(3) ,y+ψ(4) ,y an (A.36) with three definite integrals and one analytical term. The complex parameter present in the integrands of ψ(1) ,y and ψ(2) ,y is ξs=|ξ1|−iτ , a symmetric function with respect to variable ξ1 . This variable can be defined in a discrete form as ξs l= (l−1)∆ξ−iτ , for l=1,...,N+1 , and ξs 2N−l+2=ξs l , for l=2,...,N . By doing that, the first integral ψ(1) ,y , can be evaluated applying the IFFT algorithm as we did before: ψ(1) ,y=1 2πeτyj"M ∑ l=1 iξs lΨlei2π M(j−1)(l−1)#(A.37) 98 Artículo I and, on the contrary, the second integral ψ(2) ,yis approximated using the FFT algorithm: ψ(2) ,y=1 2πcosh(τyj)"M ∑ l=1 iξs lΨle−i2π M(j−1)(l−1)#(A.38) due to its negative exponential term. In the last integral ψ(3) ,y , the complex parameter ξ is defined in discrete form as ξa l= (l−1)∆ξ−iτ , for l=1,...,N+1 , and ξa 2N−l+2=−ξa l , for l=2,...,N . Substituting and performing the FFT, the last definite integral is approximated as follows: ψ(3) ,y=1 2πcosh(τyj)"M ∑ l=1 iξa lΨle−i2π M(j−1)(l−1)#(A.39) closing the derivation of a complete fundamental solution kernel. The accuracy of this approach highly depends on a proper selection of parameters τ , Ξ and N . As explained by Belibassakis [ 3 ], for the calculation of ψ(r,ro;ˆ k) using this technique, the value of τ must be small enough to make it possible to neglect the contribution of the first and third integrals of (A.23) and (A.25) , and the first term of (A.28) , but at the same time, it can not be too small because the aliasing effect is attenuated at least by factor of exp(−2τY) . On the other hand, experience demonstrates that a value of Ξ≈4−6ˆ k∗ , being ˆ k∗ the maximum value of ˆ k(x) , is large enough to make the asymptotic expression (A.12) valid, and consequently the approximation of integrals along the path C+ 3 . In our calculations, we have used a sampling of M=212 =4096 points inside the interval [−Ξ,Ξ], fixing the other two parameters to Ξ=6ˆ k∗and τ=∆ξ. FEM solution of the transformed velocity potential As we have seen, in order to evaluate the fundamental solution, it is necessary to solve the transformed one-dimensional wave transmission problem defined by equations (A.13-A.14) for different values of the Fourier parameter ξ . For this task, Belibassakis [ 3 ] proposes a second-order central finite difference scheme. In our experience, the use of the finite element method improves the solution near the source point xo , increasing this way the final accuracy of the fundamental solution. Applying the method of weighted residuals, with a test function w(x) defined in the domain [a,b], the weak form of equation (A.13) can be expressed: Zb a w(Ψ,xx +κ2(x)Ψ+δ(x−xo))dx =0,(A.40) and integrating by parts: Zb a [w,xΨ,x−wκ2(x)Ψ]dx =w(xo)+ iβ(ξ)w(b)Ψ(b)+iα(ξ)w(a)Ψ(a),(A.41) where we have substituted the matching conditions (A.14) that close the transformed domain. A.3 Fundamental solution for variable wave number 99 Using the classical finite element Galerkin formulation, the transformed velocity potential and the weighted residual function are approximated as: Ψ(x) = n ∑ j=1 Nj(x)Ψj,w(x) = n ∑ j=1 Nj(x)wj(A.42) where Ψj and wj are the corresponding nodal values, Nj(x) are linear shape functions and n is the number of nodes distributed in the domain. Substituting the discretization (A.42) into (A.41), the following FEM system is obtained:         a11 −iα(ξ)a12 ··· ··· 0 a21 a22 ··· ··· 0 . . .. . .ai j .... . . . . .. . ........ . . 0 0 ··· ··· ann −iβ(ξ)                        Ψ1 . . . Ψj . . . Ψn                =                0 . . . 1 . . . 0                ,(A.43) where the components of the tridiagonal matrix: ai j =Zb aNi,xNj,x−κ2(x)NiNjdx.(A.44) are computed with a two-point Gauss quadrature. The right hand side of system (A.43) contains the contribution of the load, where a node j is purposely located at xo . As a practical rule, we use at least 20 elements per wave-length to discretize the transformed one-dimensional transmission problems. A.3.3 Validation of the fundamental solution for constant water depth To check the accuracy of the numerical integration process described in Section A.3.2 to approximate the MSE fundamental solution, we solve the problem for constant water-depth and compare the numerical results with the analytical solution of the equivalent Helmholtz problem for constant wave number given by: ψH=i 4H(1) 0(kr),ψH,r=−i 4kH(1) 1(kr),(A.45) where H(1) 0 and H(1) 1 are Hankel functions of the first kind of order zero and one, r=|r−ro| is the distance in the radial direction and subindex H refers to Helmholtz solution. Using the series expansion of the Bessel’s function for small argument [ 1 ] the fundamental solution ψHnear the source can be expressed as: ψp H=−1 2π(ln kr 2+γ)+ i 4+O(kr)2ln(kr)(A.46) where γ is the Euler-Mascheroni constant and the superindex p denotes the polynomial approximation. The expansion contains a weak singularity of the the real part and a 100 Artículo I 0 50 100 150 200 −0.1 0 0.1 0.2 0.3 x [ m ] ψ[m2/s] Real(ψ) Imag(ψ) Real(ψH) Imag(ψH) 0 50 100 150 200 −0.1 0 0.1 0.2 0.3 y [ m ] ψ[m2/s] Real(ψ) Imag(ψ) Real(ψH) Imag(ψH) Figure A.5 Comparison of analytical and numerical fundamental solution ψfor constant water depth. Cross section of the 2D fundamental solution along the lines y=0(top) and x=0(bottom). constant imaginary value. Figures A.5 and A.6 show a comparison between the numerical fundamental solution ψ and the analytical solution ψH for constant water depth. The wave-period and water-depth used are T=5s and h=14m respectively. The solution profiles at y=0 and x=0 shown in Figure A.5 are in excellent agreement. The same degree of approximation is obtained for xand yderivatives, represented in Figure A.6. A.4 Boundary element formulation In this section we present the basis of the BEM for the Helmholtz problem, originated from an integral equation expressing a reciprocity relation between the unknown velocity A.4 Boundary element formulation 101 0 50 100 150 200 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 x [ m ] ψ,x[m/s] Real(ψ,x) Imag(ψ,x) Real(ψH,x) Imag(ψH,x) 0 50 100 150 200 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 y [ m ] ψ,y[m/s] Real(ψ,y) Imag(ψ,y) Real(ψH,y) Imag(ψH,y) Figure A.6 Comparison of analytical and numerical fundamental solution derivatives ψ,x and ψ,y for constant water depth. Cross section of the x -derivative along the line y=0(top) and the y-derivative along x=0(bottom). potential field (ˆ φ,∇ˆ φ)and the fundamental solution kernel (ψ,∇ψ). A.4.1 Direct boundary integral equation The BEM formulation for an acoustic medium is well known and can be found in different texts [ 53 , 2 ]. Let us consider a domain Ω closed by a boundary Γ of outward normal n . Multiplying both sides of Helmholtz equation (A.5) by the fundamental solution and applying Green’s second identity and Sommerfeld’s radiation condition, the following Boundary Integral Equation (BIE) is obtained in absence of internal loads: C(ro)ˆ φ(ro)+ ZΓ ∇ψ(r,ro;ˆ k)·nˆ φ(r)dΓ−ZΓ ψ(r,ro;ˆ k)∇ˆ φ(r)·ndΓ=0,(A.47) 102 Artículo I where ∇ˆ φ·n=∂ˆ φ/∂n=ˆq is the normal flux, ro is the collocation point and C(ro) is a geometrical coefficient, function of the regularity of the boundary, that takes the value C(ro) = θ(ro)/2πwith θ(ro)being the internal angle of the boundary at point ro. This BIE can be solved numerically discretizing the boundary into ne elements, Γ= ne S i=1 Γe and approximating the fields of modified velocity potential ˆ φ(r) and normal flux ˆq(r)using isoparametric linear elements: ˆ φ(ζ) = n ∑ i=1 Ni(ζ)ˆ φi,ˆq(ζ) = n ∑ i=1 Ni(ζ)ˆqi,(A.48) where ˆ φi and ˆqi represent the nodal values of velocity potential and flux, Nj(ζ) are the element shape functions and n is the number of nodes per element. Introducing this approximation in (A.47), the discretized form of the BIE can be written as: Ciˆ φi+ ne ∑ e=1ZΓe ∇ψ(r,ri;ˆ k)·nNjˆ φjdΓe= ne ∑ e=1ZΓe ψ(r,ri;ˆ k)NjˆqjdΓe,(A.49) an equation that can be expressed in matrix form as: ˆ Hˆ φ φ φ=ˆ Gˆ q,(A.50) being ˆ φ φ φ the nodal vector of modified velocity potentials ˆ φj , ˆ q the nodal vector of fluxes ˆqj, together with the matrix coefficients: ˆ Hi j =Ciδi j + ne ∑ e=1ZΓe ∇ψ(r,ri;ˆ k)·nNjdΓe(A.51) ˆ Gi j = ne ∑ e=1ZΓe ψ(r,ri;ˆ k)NjdΓe(A.52) where δi j is the Kronecker δ-function and rithe position of node i. In general, the boundary element integrals present in (A.51) and (A.52) can be computed using a standard Gauss quadrature formula. But when the collocation point is located in one of the element nodes, integral (A.52) becomes weakly singular and needs a special treatment. In that particular case, it is possible to express the integrand as the sum of two terms: ˆ Gi j =ZΓe [ψ(r,ri;ˆ k)−℘(ˆ kir)]NjdΓe+ZΓe ℘(ˆ kir)NjdΓe(A.53) one completely regular treated by a standard Gauss quadrature and another term of order O(lnr) that can be evaluated numerically using a special quadrature. We have used the real part of the first term of the series expansion of ψH around ri , obtained in (A.46) , to define the kernel: ℘(ˆ kir) = −1 2π(ln ˆ kir 2+γ)(A.54) A.4 Boundary element formulation 103 x[m] 50 60 70 ψ[m2/s] -0.4 -0.2 0 0.2 0.4 0.6 y[m] 0 10 20 -0.4 -0.2 0 0.2 0.4 0.6 Figure A.7 Behavior of the Green’s function ψ(r,ro;ˆ k) (thin solid line), the kernel ℘(ˆ kor) (dashed line) and the regularized integrand (thick solid line). The source point is located at xo=60m and the bathymetry is defined in Section A.5.1 with a wave period T=5s. needed to regularize the first integral, see Figure A.7. Details of the numerical treatment of the second integral can be found in [20]. Finally, back-substituting the change of variable introduced in (A.4) , the BEM system (A.50) can be rewritten in terms of the original velocity potential as: Hφ φ φ=Gq,(A.55) a complex non-symmetrical linear system that provides the solution of the problem. A.4.2 Integral formulation for scattering problems In wave transmission problems, it is usually interesting to consider the effect on the object under study of an incident wave radiated from a distant source. In these scattering problems [53], we divide the total velocity potential into a scattered wave and an incident wave: ˆ φ(r) = ˆ φsc(r) + ˆ φin(r)(A.56) where ˆ φsc is the scattered field radiated by the object that should satisfy the Helmholtz equation and ˆ φin is the incident field that would exist in the absence of obstacles. If we substitute this decomposition in (A.47), the BIE adopts the new form: C(ro)ˆ φ(ro)+ZΓ ∇ψ(r,ro;ˆ k)·nˆ φ(r)dΓ−ZΓ ψ(r,ro;ˆ k)∇ˆ φ(r)·ndΓ=ˆ φin(ro).(A.57) 110 Artículo I −50 0 50 100 0 0.5 1 1.5 2 x[m] |φ|/|φo| BEM Analytical solution Figure A.12 Wave scattering by a cylinder in constant water depth. Comparison with the analytical solution of the normalized magnitude of the velocity potential |φ|/|φo|at y=0obtained using BEM. waves with a shallowness ratio varying between 0.30 and 0.074 from the deepest area to the shallow region. In Figure A.16 (top) we represent the contours of equal-amplitude lines of normalized wave height obtained over the elliptical shoal using the BEM-FEM coupling technique. Two sections are defined in the domain to compare the solution with experimental results obtained by Berkhoff et. al. [ 8 ] and numerical results from Belibassakis et. al [ 5 ] solving the MMSE for this case, see Figure A.17. Equal phase-range zones of the solution are shown at the bottom of Figure A.16, where it can be appreciated its continuity and absence of spurious reflections at the BEM-FEM interface. Sections A-A’ and B-B’ of the solution are represented in Figure A.17, demonstrating the good behavior of the fundamental solution and presenting a very reasonable agreement with the experimental data. The discrepancies between numerical and experimental results, as the overestimation of the focal peak in the transversal section B-B’ and the lack of prediction along the longitudinal section A-A’, are mainly attributed to the non-linear effects [ 5 ], not considered in the MSE linear theory. The phase results represented in Figure A.16 (bottom) are in very good agreement with those measured by Berkhoff et al. [8], capturing the presence of two low amplitude points behind the shoal observed in the experiments. A.5.4 Harbor resonance study Finally, in order to explore the possibilities of the proposed BEM formulation in applications with more complex geometries, it is studied the wave transmission problem in a small A.5 Numerical examples 111 −50 0 50 100 −100 −80 −60 −40 −20 0 20 40 60 80 100 x[m] y[m] −2−1 0 1 2 Figure A.13 BEM solution for a decreasing water depth in the x -direction and incidence angle θ=π/2 . Contour field of the real part of the normalized velocity potential. Shadow region behind the cylinder is now affected by a strong shoaling effect due to the variable bathymetry. harbor and its near-by coastal region with variable bathymetry. The harbor is located in the coastal village of Chipiona, at coordinates 36◦4405700 latitude and 6◦2504200 longitude on the Atlantic shoreline of the Cádiz province, southwest of Spain. The objective of the simulation is to estimate the wave amplification in the interior of the harbor and reproduce the diffraction effects in front of the breakwater. The geometry and boundary conditions of the numerical model are shown in Figure A.18, with open boundaries to represent the shoreline and completely reflecting boundaries in the dike and quay walls of the harbor. The open boundaries are modeled assuming a complete absorption of the incoming wave. Bathymetric lines are defined parallel to the shoreline, varying with the polynomial depth function (A.60) with coefficients a0=100m , a1=0m−1 , a2=−4.0816327×10−4m−2 and a3=3.239391×10−7m−3 . In the region, the tidal range is 3.4m with a lowest tide of 4m that is considered as mean water-depth inside the harbor. The boundary of the harbor is discretized using 1234 linear boundary elements with at least 20 elements per wave length. As external source, it is considered an incident wave of heigh 1m , angle of incidence θ=π/6 and a common period value for this area, T=10s . As part of the post-process and after solving the boundary problem, we computed the solution in the exterior domain at 72539 collocation points that are evenly distributed with a distance between two adjacent points of 1/20 times the local wave length. The WAF 112 Artículo I −50 0 50 100 −100 −80 −60 −40 −20 0 20 40 60 80 100 x[m] y[m] −2−1 0 1 2 Figure A.14 BEM solution for a decreasing water depth in the x -direction and incidence angle θ=π/6 . Contour field of the real part of the normalized velocity potential. Combined shoaling and refraction effects behind the cylinder due to the variable bathymetry. obtained for this example is shown in Figure A.19. It can be observed that some points inside the harbor exhibit a significant amplification factor because dissipation effects, like friction with the sea floor and partial-reflecting boundaries, are not considered in the analysis. However, the use total reflection boundaries is a common practice in harbor resonance studies. A.6 Summary and Conclusions A complete fundamental solution and BEM formulation for the elliptic Mild-Slope equation in waters of variable depth in one direction has been presented. The Green’s function proposed by Belibassakis [ 3 ] has been used as the starting point of our BEM formulation and different examples have been solved to validate the approximation. The main conclusions and findings of this work are the following. • The Green’s function of Belibassakis [ 3 ] for the MSE with one-directional variable bathymetries has been extended and combined with a boundary element formulation to simulate surface water-wave transmission problems in medium to shallow transition waters. A.6 Summary and Conclusions 113 −50 0 50 100 0 0.5 1 1.5 2 2.5 x[m] |φ|/|φo| BEM solution for θ=π/2 BEM solution for θ=π/6 Analytical solution (h=14m,θ=π/2) Figure A.15 Wave scattering by a cylinder in variable water depth. Normalized magnitude of the velocity potential computed with the BEM along the profile located at y=0for two different incidence angles. • The one-dimensional wave equations in the transformed domain associated with the evaluation of the fundamental solution are solved using a classical Galerkin finite element approximation. It has been observed in our numerical experiments that similar accuracy is obtained using this approach, compared to the second order FD scheme proposed by Belibasakis [3] for this task. • The formulation is able to correctly reproduce the phenomena appearing in waterwave transmission problems: shoaling, diffraction, refraction and the result of their combined effects. • Bathymetries with slopes up to 1:3 and contour lines parallel to the shoreline are very common in real problems. The proposed BEM technique allows to simulate these conditions, providing accurate solutions for practical coastal engineering problems. • This BEM formulation can be coupled with classical design techniques, like FEM or FDM formulations of the MSE or MMSE [ 43 ], to model open sea conditions of variable bathymetry. Finally, it is important to mention that the proposed BEM formulation can also be combined with more advanced partial-reflection boundary conditions and FEM-BEM coupling techniques to extend the range of practical applications. Work in this direction is under way. 114 Artículo I −5 0 5 10 15 −8 −6 −4 −2 0 2 4 6 8 x[m] y[m] A B B’ A’ 0 0.5 1 1.5 2 2.5 Figure A.16 Scattering produced by an elliptic bank on a sloping bottom. Absolute value of the normalized wave height (top) and zones of equal phase within an interval of π radians (bottom). Two thick lines indicate the location of the sections represented in Figure A.17. A.6 Summary and Conclusions 115 0 2 4 6 8 10 12 0 0.5 1 1.5 2 2.5 3 [m] Normalized wave height Lab. data MMSE B BEM−FEM MSE −5 0 5 0 0.5 1 1.5 2 2.5 3 [m] Normalized wave height Lab. data MMSE B BEM−FEM MSE Figure A.17 Normalized wave height along a longitudinal section A-A’ (top) and a transversal section B-B’ (bottom), see Figure A.16. Comparison of experimental data from Berkhoff et al. [ 8 ], numerical results based on the MMSE obtained by Belibassakis et al. [ 5 ] and numerical solution of the MSE obtained by BEM-FEM coupling. 116 Artículo I y x φin θ Reflecting boundaries Open boundaries Water depth Figure A.18 Nearshore bathymetric mapping and boundary conditions used for the harbor simulation. Water depth changes along the x -coordinate from 100m to 4m between x=0m and x=840m . Water depth remains constant for x≥840m . A.7 Acknowledgements This work was supported by the Ministerio de Economía y Competitividad of Spain, under the research projects DPI2010-19331and DPI2013-43267-P, which were co-funded by European Regional Development Funds (ERDF). A.7 Acknowledgements 117 Figure A.19 WAF diagram calculated using BEM for the harbor of Chipiona example under inclined incident waves and 10-second period wave condition. Bibliography [1] M. Abramowitz and I.A.Stegun. Handbook of mathematical functions. NBS, 10 edition, 1972. [2] N. Atalla and F. Sgard. Finite Element and Boundary Methods in Structural Acoustics and Vibration. CRC Press, 2015. [3] K. Belibassakis. The Green’s function of the mild-slope equation: The case of a monotonic bed profile. Wave Motion, 32:339–361, 2000. [4] K. Belibassakis and G. Athanassoulis. Three-dimensional Green’s function for harmonic water waves over a bottom topography with different depths at infinity. Journal of Fluid Mechanics, 510:267 – 302, 2004. [5] K. Belibassakis, G. Athanassoulis, and T. Gerostathis. A coupled-mode model for the refraction-diffraction of linear waves over steep three-dimensional bathymetry. Applied Ocean Research, 23(6):319 – 336, 2001. [6] G. Beltrami, G. Bellotti, P. Girolamo, and P. Sammarco. Treatment of wave breaking and total absorption in a mild-slope equation fem model. Journal of Waterway, Port, Coastal, and Ocean Engineering, 127(5):263–271, 2001. [7] P. G. Bergmann. The equation in a medium with a variable index of refraction. The Journal of the Acoustical Society of America, 17(4):329 – 333, 1946. [8] J. Berkhoff, N. Booy, and A. Radder. Verification of numerical wave propagation models for simple harmonic linear water waves. Coastal Engineering, 6(3):255 – 279, 1982. [9] W. Berkhoff. Computation of combined refraction-diffraction. In Proceedings of 13th International Conference on Coastal Engineering. ASCE, 1972. [10] W. Berkhoff. Mathematical Models for Simple Harmonic Linear Water Waves. Wave Diffraction and Refraction. PhD thesis, Delft Hydraulics Laboratory, 1976. 119 126 Artículo II made for non-lineal equations like the Boussinesk-type equations [ 19 , 20 ]. Nevertheless, the application of pseudo-spectral methodologies to the mild-slope equation has not been treated in the literature, despite of the fact that these type of models represent the basic framework for the simulation of surface wave propagation problems in variable water depths. The elliptical mild-slope equation (MSE) [ 9 , 10 ] allows to consider simultaneously the effects of diffraction, refraction, reflection and shoaling of linear water surface waves, being formally only valid for slowly varying sea-bed slopes, i.e. ∇h<< kh , being h the water depth and k the wave number [ 60 , 13 ]. That condition converts the MSE in a powerful tool to study large coastal areas with less computational effort compared to non-lineal models. In order to extend the applicability of the MSE, some modifications has been proposed for more complicated bathymetries, as in [ 54 ] [ 48 ], or to include energy dissipation effects, such as wave breaking and bottom friction [ 46 ]. The modified mildslope equation (MMSE) presented in [ 15 ] and later improved in [ 54 ], retains the second order terms discarded by Berkhoff in the original formulation of the MSE. More extensions of the MSE can be found in [58, 16, 38, 41, 28, 29]. The MSE has been traditionally solved using the finite element method (FEM) [ 10 ], and the finite difference method (FDM) [ 52 , 42 ]. However, classical finite difference and finite element methods present common deficiencies, as the pollution effect and difficulties to reproduce properly natural boundary conditions. The pollution effect is directly related to the numerical dispersion, which is an error due to the loss of the ellipticity of the Helmhotz’s equation as the wave number k increases. This error deteriorates the solution even if the coefficient kh , where h is the element size, is kept small. Hence the number of nodes per wave length is not sufficient to determine the accuracy of the solution. This matter has been object of many works [ 3 , 31 , 32 , 18 ], where the pollution error is in general estimated as a rational function of kh with the coefficients depending on the order of the interpolation function p . It is found that with higher p , the number of points per wave length required for a given accuracy presents a slower increasing rate. In order to overcome this deficiency several formulations have been proposed, as the Galerkin Least-squares (GLS) methods [ 59 ], spectral element methods [49] or the smoothed finite element methods [43]. One of the most important problems of the FEM and FDM formulations is in the difficulty to reproduce unbounded domains properly. This deficiency has been studied by many authors. Bettess and Zienkiewicz [ 11 ] and Lau and Ji [ 37 ] used infinite elements in the outer regions. Dirichlet to Neumann (DtN) boundary conditions were proposed by Givoli et al. [ 23 , 35 , 22 ] as an analytical procedure to reproduce exact non-reflecting boundary conditions in some particular cases. This idea was followed by Bonet [ 12 ] to derive the discrete non-local (DNL) boundary condition. Other methods, such perfectly matched layers (PML) has been used in [ 21 ] and [ 50 ]. More rudimentary iterative methods have also been proposed to define absorbing boundary conditions in [6, 57, 17, 45]. Boundary element techniques prove to be very accurate in wave-diffraction problems in unbounded domains, presenting the additional benefit that the radiation condition to infinity is automatically satisfied. In order to improve the solution of the FEM schemes, Hauguel [ 25 ] and Shaw and Falby [ 56 ] proposed coupling FEM and BEM at the end of the seventies. Hamanaka [ 30 ] and later Isaacson and Qu [ 33 ] introduced a boundary integral formulation B.1 Introduction 127 to reproduce the wave field in harbors with partial reflecting boundaries and Lee et al. [ 40 , 39 ] included the effect of incoming random waves. Boundary integral techniques have also been used for non-constant bathymetries by means of the dual reciprocity boundary element method (DRBEM). Some of these works also include wave run-ups, see Zhu [ 63 ]. Later, this technique was extended to model internal regions with variable depth surrounded by exterior regions with constant bathymetry [ 44 , 64 , 65 , 27 ]. More recently, Naserizabeh et al. [ 51 ] proposed a coupled BEM-FDM formulation to solve the MSE in unbounded domains. Later, Cerrato et al. [ 14 ] derived a complete kernel of a fundamental solution for variable water depths, based on the Green’s function presented in [ 4 ], and introduced it in a standard BEM formulation. In the context of SEM, we have the contributions of Mehdizadeh and Paralchivoiu (2003) [ 49 ] or, more recently, He et al. (2016) [ 26 ]. These works presented a SEM formulation based on the Legendre polynomials for the two dimensional Helmholtz’s equation using special techniques to reproduce unbounded domains. To model the open boundaries, Mehdizadeh and Peralchivoiu use an absorbing layer, or PML, surrounded the inner region while He et al. used the Dirichlet to Neumann (DtN) boundary conditions. This last work is also concerned with the modelling of waves propagating over layered media. An specific work about water wave modelling with BEM combined with spectral approximation based on Chebyshev polynomials was presented recently by Kumar et al. [ 36 ], but their formulation is limited to constant bathymetries. The combination of spectral elements and BEM formulations has been also recently applied to other fields, such as elastodynamic promblems [ 66 ]. A final interesting work is [ 61 ], where the authors study the hp -convergence of the spectral element method for the two dimensional Helmholtz’s equation. In order to deal with the issues mentioned above, this work presents a spectral boundary element-finite element formulation for linear water-wave propagation problems. The SEM is used to model inner regions, allowing to model bathymetries with arbitrary shape, while the BSEM is used as a natural boundary condition which fulfils automatically the Sommerfeld radiation condition. An important feature of the BSEM formulation is that it uses the fundamental solution derived by Cerrato et al. [ 14 ] to consider a variable bathymetry in the outer region. The approach of the solution is made by considering a pseudospectral approximation inside the elements, using nodal basis functions with nodes located at the LGL points. The paper is organized as follows. First, Section B.2 presents the MSE for the linear water-wave propagation problem. In Section B.3, the boundary spectral element formulation is explained in detail. In Section B.4, special attention is paid to the finite spectral element formulation. The coupled boundary spectral element–finite spectral element formulation is presented in Section B.5. Section B.6 is dedicated to the validation of the proposed formulation through the solution of different wave propagation problems for variable water depth, where mathematical and numerical convergence studies are presented. Finally, the paper concludes with the summary and conclusions. 128 Artículo II B.2 The mild-slope equation To formulate the MSE, according to [ 9 , 10 ], we consider a Cartesian coordinate system with the (x,y) -plane located on the quiescent water surface and the z direction pointing upwards. The still water depth is given by hw(x,y) and ∇= (∂x,∂y) is used to represent the gradient operator. Under the assumption of potential flow and integrating the velocity potential in the vertical direction with appropriated boundary conditions, the velocity potential of the water surface can be assumed to be of the form: Φ(x,y,t) = φ(x,y)e−iωt,(B.1) where i is the imaginary unit and t is the time variable. This potential has to satisfy the homogeneous MSE, that may be written as: ∇·(ccg∇φ)+ ωcg cφ=0.(B.2) where c is the wave velocity and cg is the group velocity. The water depth function hw(x,y) , the wave number k and the angular frequency ω of the waves are related by the dispersion equation: ω2=gktanh(khw) , where g is the gravitational acceleration ( g=9.81m/s2 ). So, for a fixed frequency and variable bathymetry, the wave number k(x,y) is a function of the local water depth. In the linear wave theory, the wave height ( H ), is linearly related to the velocity potential on the water surface by the following relation: H=2ω|φ|/g,(B.3) an important design variable for practical problems. The MSE can be simplified introducing the following change of variable due to Bergmann [7]: φ=ˆ φ/√ccg. This relation transforms (B.2) into a Helmholtz equation: ∇2ˆ φ+ˆ k2ˆ φ=0,(B.4) with a modified wave number ˆ k(x,y)given by: ˆ k2(x,y) = k2−∇2√ccg √ccg ,(B.5) that is a known function of the wave characteristics and the local water depth. In this work, we discard the second term and the modified wave number is approached as ˆ k=k(x,y). Note that a similar change of variable can be done for treating the same problem under the framework of the MMSE. Simply by modifying accordingly the expression of the wave number (B.5) , to account for additional effects associated with higher-order contributions of bottom slope and curvature, we obtain a MMSE model that extends the applicability of the mild-slope equation. B.3 Boundary spectral element method 129 B.3 Boundary spectral element method The Boundary Integral Equation (BIE) is a powerful tool to construct accurate and efficient formulations for linear water wave problems. The BIE presents some important advantages as a the fullfilment of the Sommerfeld radiation condition at infinity, or the strong reduction on the number of degrees of freedom required to construct the discrete model of the domain. The main disadvantage of the application of the BIE to water wave propagation problems is the difficulty to include variable bathymetries. The prefered technique to deal with this difficulty is the DRBEM. However it requires to compute domain integrals and the modelisation of infinite domains continues being only applicable to constant-depth areas. Nevertheless, in order to model infinity domains with non-homogeneus water-depth, we can use the standard BIE but equipped with an appropriated fundamental solution [14]. In the BSEM, the BIE is discretized by spectral elements. As in other spectral element techniques, the nodal basis functions based on a family of orthonormal polynomials provide to the BIE an exponencial convergence and highly accurate description of the boundary shape. B.3.1 Boundary integral equation We start writing the classical BIE of the BEM in order to fix the notation. Considering a boundary Γ of the region of interest, which is governed by equation (B.4) , the BIE can be expressed as: C(x0)ˆ φ(x0)+ ZΓ ∇ψ(x,x0;ˆ k)·nˆ φ(x)dΓ=ZΓ ψ(x,x0;ˆ k)ˆq(x)dΓ+ˆ φin,(B.6) where x0 is the collocation point, ˆq is the normal flux through the boundary considering an outward normal and ψ(x,x0;ˆ k) is the Green’s function of the Helmholtz problem. The free term is, in general, given by C(x0) = θ(x0)/2π , being θ(x0) the internal angle of the boundary at the collocation point. In order to consider also scattering problems, the incident potential ˆ φin is included in the above expression [62]. B.3.2 Spectral boundary elements In order to evaluate the integrals of equation (B.6) the boundary is divided into ne elements, being the discrete boundary Γ=Sne i=1Γe . In the BSEM formulation presented here, the geometry and the variables on the boundary are approximated using nodal basis functions based on the Lagrange polynomials with nodes located on the Legendre-Gauss-Lobatto (LGL) quadrature points. Inside each spectral element, the modified velocity potential and the normal flux through the boundary is approximated by: ˆ φ(ξ) = p ∑ m=0 Lm(ξ)ˆ φ(e) m,ˆq(ξ) = p ∑ m=0 Lm(ξ)ˆq(e) m,(B.7) where ˆ φ(e) m and ˆq(e) m are the nodal values of the element e and Lm(ξ) the nodal basis functions of degree p on the interval [−1,1] . This approximation functions can be constructed by 130 Artículo II n Γe x y ξ=−1ξ=1ξ Lm(ξ) Figure B.1 Example of a spectral boundary element. An element of order p=7is represented in the global x−y -axis (left) and in the local coordinate ξ (rigth). Also the nodal function corresponding to the fifth node is shown. means of Lagrange polynomials with roots at the Gauss-Lobatto points: Lm(ξ) = p ∏ l=0(l6=m) ξ−ξl ξm−ξl ,m=0,1,...,p,(B.8) being ξm the local coordinate of the node under consideration and ξl the coordinate of the rest of the nodes inside the element. In the above expressions the index 0 and p are used to denote the nodes located at the extremes of the interval ( ξ=−1 and ξ=1 ). A representation of this spectral boundary element is given in Figure B.1. It is possible to use the same high order polynomials to approximate the geometry, specially when the real geometry is curved and complex. B.3.3 Discrete boundary integral equation Introducing in the BIE given by (B.6) the approximations of the variables (B.7) , the discrete form of the BIE can be written as Ciˆ φi+ ne ∑ e=1Z1 −1 ∇ψ(x(ξ),xi;ˆ k)·nLm(ξ)J(e)(ξ)dξˆ φ(e) m= ne ∑ e=1Z1 −1 ψ(x(ξ),xi;ˆ k)Lm(ξ)J(e)(ξ)dξˆq(e) m+ˆ φin i (B.9) where the sub-index i refers to the collocation point and the super-index e to the elements. The jacobian of the transformation from a global coordinate to the local coordinate system is denoted by J(e)(ξ) . Assembling all the element contributions included in expression (B.9), the discrete form of the BIE can be finally written in a matrix form as: ˆ Hˆ φ φ φ=ˆ Gˆ q+ˆ φ φ φin,(B.10) B.3 Boundary spectral element method 131 where the vectors ˆ φ φ φ and ˆ q contain the unknown nodal values of the modified potential and flux respectively, and ˆ φ φ φin is the load vector with the values of the incident wave. The matrices ˆ Hand ˆ Gare composed of the following terms: ˆ Hi j =Ciδi j + E(j) ∑ e=1Z1 −1 ∇ψ(x(ξ),xi;ˆ k)·nLj(ξ)J(e)(ξ)dξ,(B.11) ˆ Gi j = E(j) ∑ e=1Z1 −1 ψ(x(ξ),xi;ˆ k)Lj(ξ)J(e)(ξ)dξ,(B.12) where E(j)represents the group of elements sharing node j. B.3.4 Numerical integration The integrals of the ˆ Hi j and ˆ Gi j terms are evaluated numerically. One of the principal advantages of using a Lagrangian approach combined with a LGL distribution of the nodes within the element, is that the numerical integration can be carried out by means of its associated LGL quadrature rule, which is exact for polynomials up to 2p−1 . This is specially accurate for those elements that are not located close to the collocation point, over a region with smooth variation of the wave number and without a complex geometry. In this case, the values of the Lagrange interpolants at the quadrature points are given by the relation Lm(ξk) = δmk,(B.13) being δ the Kronecker delta. Using this property, the integral of the product of a function f(ξ) with the nodal basis functions Lm(ξ) reduced to a single evaluation of the integrand at the quadrature point as follows: Z1 −1 f(ξ)Lm(ξ)dξ=f(ξm)wm,(B.14) where the weight is given by the expression: wm=2 p(p+1) 1 P2 p(ξm),(B.15) and where Pp(x) is the Lagrange polynomial of order p . This way, the values of the integrals inside (B.11) and (B.12) can be computed fast and efficiently. In other cases, like for example integration on elements close to the collocation point, over regions where the geometry contains a sharp corner, or the wave number changes with a complex function, a standard Gaussian quadrature is preferred. B.3.5 Regularization of the singular integrals During the collocation process, when the element under integration contains the source load, the integrand of (B.12) becomes weakly singular. This problem can be treated by a regularisation technique; subtracting the singularity and integrating it separately [ 2 , 1 , 24 ]. 132 Artículo II ξ=−1ξ=1 ξ=ξ0 ε SεS V(x,x0) ρ=ˆ ρ+ ρ=ˆ ρ−ρ Figure B.2 Singular integral around Sε when the collocation point is located inside the element. The distance to the singularity ρ=|ξ−ξ0| is defined to extract the singular kernel. Let us define the S and the Sε regions inside the element as shown in Figure B.2, where Sε surrounds the singularity and S covers the entire element. Taking the limit when ε tends to zero, the singular integrals in (B.12) can be calculated by means of their Cauchy Principal Value: I=lim ε−→0ZS−Sε V(x,x0)−V1(x,x0)dξ+ZS−Sε V1(x,x0)dξ(B.16) where the integrands are V(x,x0) = ψ(x(ξ),x(ξ0))Li(x(ξ))J(x(ξ)) and V1(x,x0) is the first term of the Taylor series expansion of V(x,x0)at ξ=ξ0. To subtract the singularity, we use the fundamental solution of the 2D-Helmholtz equation for constant wave number in terms of its series expansion: ψ(x,x0) =℘(kρ)+ i 4+O((kρ)2log(kρ)),(B.17) where ℘(kρ) = −1 2π(log(kρ/2)+ γ)(B.18) and ρ=|x−x0| is the distance to the load. Then, in order to express the weakly singular kernel of (B.16) when ρ goes to zero, the function containing the singularity (B.18) is written in terms of the local variable of the element ρ=|ξ−ξ0| , as shown in Figure B.2. Then, the distance from x0 in any direction xi can be approximated using the first order B.3 Boundary spectral element method 133 term of its Taylor series expansion: xi−x0 i=dxi dξξ=ξ0 ρ+O(ρ2).(B.19) Therefore, the distance between x0and any point xcan be finally expressed as ρ=ρA+O(ρ2),(B.20) where A2=∑2 i=1xi,ξξ=ξ02 contains the Jacobian of the transformation. Substituting (B.20) into (B.18) we obtain the kernel of the singularity in terms of the local distance ρ: ℘(kρ) =℘ o+℘ 1logρ+O(ρ),(B.21) where℘ o=−1 2π(γ+log kA 2)and℘ 1=−1 2π. The same procedure can be applied to all the terms appearing in the integrand, like: L(ρ) = Lo+O(ρ)(B.22) and J(ρ) = Jo+O(ρ),(B.23) being Lo=Li(ξ0)and Jo=J(ξ0). Collecting all the terms, we can write: V(x,x0) = To+T1logρ+O(ρ)(B.24) where To=℘ oLoJo and T1=℘ 1LoJo . Therefore, the first term of the series expansion of the integrand V(x,x0)can be written as: V1(x,x0) = To+T1logρ.(B.25) Introducing this in equation (B.16) we regularize integral in the domain Sε. To compute the first regular integral of (B.16) we can use a standard Gaussian quadrature. The second integral, which contains the singularity, can be computed analytically dividing the domain Sinto two regions: ZS V1(x,x0)dξ=Z0 ˆ ρ−To+T1logρdρ+Zˆ ρ+ 0 To+T1logρdρ,(B.26) and integrating, ZS V1(x,x0)dξ=2To+T1ˆ ρ−log ˆ ρ−+ˆ ρ+log ˆ ρ+−2,(B.27) being ˆ ρ−=|−1−ξ0|and ˆ ρ+=|1−ξ0|. 134 Artículo II B.4 Spectral element method The SEM is a high-order finite element technique that combines the geometrical flexibility of finite elements with the high accuracy of spectral methods. This method was pioneered in the mid 80’s by Patera [ 53 ] for the Navier-Stokes equation and exhibits several favorable computational properties, such as the use of naturally diagonal mass matrices that facilitate the iterarive solution process. Moreover, in the context of Helmholtz’s equation, the SEM reduces the pollution effect and improves the efficiency of the classical Galerkin finite element approximation, i.e., it requires a lower number of nodes per wavelength compared to FEM in order attain the same accuracy [49]. B.4.1 Weak form of the MSE The weighted residual form of the Helmholtz’s equation (B.4) for a domain Ω bounded by Γafter integrating by parts is expressed as: ZΩ ∇η∇ˆ φdΩ−ZΓ ηˆq dΓ−ZΩ ˆ k2ηˆ φdΩ=0,(B.28) where η is the weight function. This variational formulation is the starting point of the SEM. B.4.2 Two dimensional spectral elements A discretization of the domain Ωis carried out by using nespectral elements, fulfilling Ω= ne [ e=1 Ωeand ne \ e=1 Ωe=/0,(B.29) where the modified velocity potential is approached inside the elements by means of Lagrangian interpolants through the LGL grid points in the form: ˆ φ(ξ,ζ) = p ∑ m=0 p ∑ n=0 Lm(ξ)Ln(ζ)ˆ φ(ξm,ζn).(B.30) where Lm is the one dimensional shape function associated to node m as in (B.8) , p is the polynomial degree and (ξ,ζ)∈[−1,1]2 are the intrinsic coordinates of the reference element. Figure B.3 shows an example of a two-dimensional spectral element and its transformation to the normalized coordinates (ξ,ζ) . Also two the nodal basis functions corresponding to L21(ξ,η) = L0(ξ)L4(ζ)and L6(ξ,η) = L0(ξ)L1(ζ)are represented. Equation (B.30) can be written in a more compact form as ˆ φ(ξ ξ ξ) = Nc ∑ j=1 Lj(ξ ξ ξ)ˆ φj,(B.31) B.4 Spectral element method 135 x y ζ ξ ζξζξ L0(ξ)L4(ζ)L0(ξ)L1(ζ) Figure B.3 A finite spectral element of order p=4 is represented in the physical ( x−y )- coordinate system and in the normalized (ξ,ζ)∈[−1,1]2 reference coordinate system (top). Two representative nodal basis functions of the element are also showed (bottom). where ξ ξ ξ= (ξ,ζ) , Nc is the total number of element nodes ( (p+1)2 for quadrilateral elements), Lj(ξ ξ ξ) the two dimensional shape functions and ˆ φj the nodal values of the modified velocity potential. B.4.3 Discrete finite element equations Applying the Spectral-Galerkin method to the weighted residual form (B.28) and considering the domain discretization (B.29), the following system of equation is obtained: ˆ A−ˆ Cˆ φ φ φ ˆ q=0,(B.32) where vectors ˆ φ φ φ and ˆ q group all the nodal values of the modified velocity potentials and boundary fluxes respectively. The matrix ˆ A and ˆ C are obtained by assembling the elemental matrices ˆ Ae and ˆ Ce as in the standard finite element formulation. The domain integrals of (B.28) over each element eprovide the coefficients of the ˆ Aematrix: ˆ Ae i j =Ke i j −Me i j,(B.33)