Un modelo numérico de refracción-difracción de ondas en zonas costeras

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 6, 3, 409436( 1990) UN MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS EN ZONAS COSTERAS P. ORTIZ Y M. PASTOR División de Cálculo y Análisis Matemático, Cedex, MOPU, Alfonso XII, 3, 28014 Madrid. RESUMEN Se estudia en el presente trabajo el tratamiento numérico, mediante el método de elementos finitos (MEF), de la ecuación aproximada de la onda para variaciones graduales de profundidad. Es planteada una nueva aproximación a la condición de borde "mar afueran basada en el método de las características que, además, permite definir, bajo los mismos criterios, las restantes condiciones de contorno. El sistema discreto es resuelto mediante el método del gradiente conjugado precondicionado (PCG), lo que se traduce en menores requerimientos de memoria. Para demostrar la efectividad del modelo, se han incluido sus predicciones en diversos ejemplos teóricos y en el análisis de la respuesta a ondas largas de recintos portuarios. SUMMARY In this paper, the numerical solution by the finite element method (FEM) for the inild slope wave equation is considered. A new approach to the "open sea" boundary condition, based on the method of characteristics, is introduced. In the same way, the others boundary conditions can be derived. Preconditioned Conjugate Gradients (PCG) is applied to solve the discrete system, resulting in important storage saving. The numerical examples included show the effectiveness of the predictions made by the model in various theoretical probleiiis and in tlie study of harbour response to long waves. INTRODUCCION En el proyecto de obras portuarias, se debe asegurar la correcta ejecución de movimientos y tareas para las cuales un puerto ha sido concebido, por lo que es preciso conocer el efecto que las "excitaciones exteriores" provocan en él. Así, las "ondas largas" (con períodos que oscilan en general entre 30 segundos y 5 minutos) pueden producir amplificaciones de movimientos horizontales31 que deriven en atrasos en las actividades y, ocasionalmente, en accidentes por roturas de amarras y colisiones entre buques. Recibido: Enero 1989 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 409 P. ORTIZ Y M. PASTOR El fenómeno "resonante"32 se producirá cuando la respuesta del sistema, sometido a excitaciones como las citadas, presente una amplificación máxima. En consecuencia, el problema depende, además, de la geometría del recinto, mientras que su magnitud está influida por los mecanismos de disipación energética presentes en forma simultánea a él. Tanto la predicción de la respuesta de un recinto portuario a ondas largas, como el estudio de su agitación debida al oleaje, son problemas que, salvo algunas excepcione~"~~~~~~' l , no poseen soluciones analíticas, por lo que en casos prácticos es necesario recurrir a la modelización físicalg y numérica. Las dos vías más generales de establecer los principios que rigen la propagación de una onda superficial de gravedad, en una región en que existen obstáculos y variaciones de profundidad que modifican sus características son (teniendo presente su posterior tratamiento numérico): Las ecuaciones de Boussinesql, válidas para números de Ursell: U = 0(1), con U = H X2/h3, siendo H y A, respectivamente, medidas de la amplitud y de una longitud horizontal característica de la onda y h una medida de la profundidad. Bajo la hipótesis de profundidades reducidas, se tratará de la ecuación no lineal de ondas largasZ4J6 (U >> l), mientras que en el caso lineal, se referirá a la teoría de ondas de pequeña amplitudz4. La formulación potencial, que es empleada en este trabajo en el caso de aproximación lineal, al recurrir a la ecuación aproximada de ondas para variaciones graduales de profundidad6 (caracterizada por U << 1). Sin estar contemplados en forma directa en estos grupos, merecen citarse los modelos de propagación basados en una formulación parabólicaa3 que se adaptan mejor a problemas costeros que a portuarios, al ignorar el campo de ondas reflejado. En la ingeniería práctica suelen utilizarse los modelos de refracción3', que, fundados en principios de óptica geométrica, permiten efectuar una aproximación al estudio del oleaje en grandes extensiones en las cuales la refracción no sea relevante. Ahora bien, la simulación numérica de las ecuaciones de Boussinesq y de las ecuaciones no lineales de ondas largas lleva a modelos que frecuentemente son de elevado costo, lo que impide su aplicación generalizada. En esta línea se han irnplementado modelos tanto por el método de diferencias finitas1, como por el de elementos finitos4s,lz,36 , existiendo esquemas explícitos de elevada economía y estabilidad13. La obtención de soluciones numéricas a la variante "elíptica", planteada en la segunda opción, son numerosas. Debido a las ventajas que otorga el MEF al tratar zonas acotadas de geometrías complejas, los modelos existentes se pueden clasificar unificadamente'' en función de las aproximaciones propuestas al problema "exterior", es decir, a la condición de radiación, que debe ser acoplada a la solución interior. De esta forma, pueden enumerarse: - Amortiguadores planos, cilíndricos, esféricos y de orden superior3*", que presentan la ventaja de su simplicidad de aplicación a un modelo MEF. Sin embargo, su limitación se basa en que su cumplimiento, en rigor, sólo existe en inñnito. La hipótesis de profundidad exterior constante, asumida en esta opción, es también MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS adoptada en todos los modelos de similar potencia. - Soluciones exteriores analíti~aP*'~*'~, que acoplan funciones que satisfacen la condición de radiación de Somrnerfeld y la ecuación de Helmholtz (profundidad constante), destruyendo parte del carácter bandeado de la matriz de coeficientes. - Integrales de donde se mantiene la simetría del problema interior, conservando el inconveniente del método anterior. - Elementos infinito^^*^^^, que no destruyen la simetríani el carácter de matriz banda del sistema, afectando en cambio su condicionamiento, característica que influye sustancialmente si se utilizan métodos indirectos de resolución. Estas soluciones tienen, pues, como objetivo principal luego de la discretización interior, salvar la dificultad de modelar la disipación por radiación, que es la más relevante en problemas de ondas largas. En el presente estudio se formula una condición de borde basada en conceptos del método de las características, extendido al problema armónico, que permite englobar tanto la condición de radiación como la condición de bordes absorbentes de energía. Cuando el modelo es aplicado a problemas de ondas largas, la condición propuesta implica dividir el dominio total en dos zonas: la primera supuesta como "mar abierto7>, que no es incluida en el análisis, y la segunda que es la que da lugar a la discretización. Bajo la hipótesis de generación de ondas largas por transferencia de energía desde ondas de alta frecuencia con extenso desarrollo, se supone, pues, que el proceso se produce en la primera región, considerada como de "generación", mientras que en la segunda se admite que dicha transferencia ha cesado o es despreciable. En ciertos as os^^,^', el origen parece más orientado hacia fenómenos de fluctuación barométrica. Por otra parte, la discretización de la región interior depende de la longitud de onda "menor", pudiéndose veriñcar mediante experimentación numéricaag, que es necesario emplear un mínimo de 8 a 12 puntos por longitud de onda para aproximar la función solución con errores relativos menores a 0.1%. Esta condición es crítica en numerosos casos prácticos, pues lleva a sistemas con elevado número de grados de libertad, que limita fuertemente el uso de algoritmos de resolución frecuentes en el MEF, como el método frontal2'. Los métodos indirectos de resolución se presentan como una eficaz alternativa, debido a sus drásticos ahorros de memoria central requerida para almacenar los coeficientes de la matriz que surge de la discretización. El presente problema tiene la particularidad de estar deñnido por matrices de coeficientes no definidas positivas, por lo que muchos métodos iterativos clásicos no son convergentes. Así, se presenta en este trabajo una extensión del método del gradiente conjugado precondicionado a sistemas con matrices de coeficientes complejos no deñnidas positivas. De esta forma se consigue aplicar el modelo a problemas de elevado número de grados de libertad con equipos medios e implementarlo en ordenadores personales con facilidad, permitiendo acceder al ingeniero proyectista a una potente herramienta de diseño. P. ORTIZ Y M. PASTOR La forma simple y compacta de concebir las condiciones de borde del problema otorga un mejor condicionamiento de la matriz de coeficientes en comparación con los métodos citados previamente, lo que influye en forma decisiva en la eficiencia del método indirecto de solución. Con el fin de demostrar el correcto funcionamiento del modelo propuesto, se han ejecutado pruebas de casos con solución teórica conocida o con resultados experimentales fiables, mientras que los casos prácticos incluidos tienen especial énfasis en el estudio de problemas de ondas largas. TEORIA DE ONDAS. FORMULACION POTENCIAL A continuación se considerarán ondas de gravedad superficiales que se propagan en un medio homogéneo e isótropo. Tanto el fondo como los bordes sólidos se admitirán como impermeables, por lo que no existirán corrientes netas de entrada o salida a través de ellos, excepto cuando se especifiquen contornos absorbentes o con velocidades normales prescritas. El fluido es incompresible y el flujo no viscoso. Las singularidades en el campo de velocidades, tales como las producidas en entornos de bordes agudos, no son aceptadas en la teoría elegida25. La adopción de la teoría lineal implica que las sobreelevaciones a partir del nivel medio del agua, q, son mucho menores que la profundidad h y, además, que la relación entre altura H y longitud de onda X debe ser pequeña, equivalente a afirmar que quedará invalidada para peraltes de onda considerables. La ecuación aproximada de ondas para variaciones graduales de profundidad exige que el parámetro p cumpla: con k el número de onda. Acepta además cualquier rango de relación profundidadlongitud de onda kh, siendo válida, por lo tanto, desde ondas cortas a ondas en profundidades reducidas. La ecuación aproximada, obtenida por BerkhofF, y luego tratada por diversos auto re^^^-^*^^, está basada en la teoría del potencial de velocidades. Por lo tanto, es asumido que el flujo es irrotacional, pudiendo definirse un potencial de velocidades 4 tal que: siendo ü el vector de velocidad y V el operador gradiente. Así, se cumplirá la ecuación de Laplace: El problema queda definido si a (1) se agrega la condición de borde en el fondo impermeable: MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS a4 ah a4 -.- +- = O, en 23 = -h(xi) (i = 1,2) azi azi ax3 y la condición de contorno a superficie libre (linealizada): a4 a24 +2= 1 apa 9- --- en 23 = O az3 at p at considerando las coordenadas cartesianas (zi, i = 1,2,3) o (2, y, z) indistintamente, siendo el plano x - y o xi - 22 coincidente con el nivel medio del agua, g la aceleración de gravedad terrestre, p la densidad del agua y pa la presión atmosférica. Si se emplean estas expresiones, sumándose la hipótesis de periodicidad de las variables que intervienen en el problema, es decir, que la dependencia de ellas del tiempo puede separarse de la dependencia de las variables espaciales xi, como por ejemplo: 4(x, t) = 4(x) . exp(.iwt) con w la frecuencia angular, y x el vector posición, se obtiene la ecuación aproximada de ondas. Para ello se supone que con variaciones graduales de profundidad, la relación de dispersión: w2 = gvkstanhkh y la solución al problema a profundidad constante: -ig. ~(x, y) cosh[k(z + h)] ~o(z,Y, z) = cosh kh (2) tienen aún validez, adoptándose los valores locales de k y h. En (2), q(x, y) satisface la ecuación de Helmhotz (a profundidad constante): v2~ + k2 -7 = O mientras que representa el potencial de velocidades tridimensional. Luego de algunas operacione~~~, la ecuación de la onda que se obtiene es, en función del potencial: V(c. cgV4) + u"-) o c (3) donde c es la celeridad de fase local, definida como: W c=- k mientras que cg es la velocidad de grupo, dada por: 1 cg = c-n; n = ?(l + 2kh sinh 2kh 1 P. ORTIZ Y M. PASTOR En el caso de profundidades reducidas, el factor kh es pequeño, por lo que: sinh2kh x 2kh, tanhkh FZ kh -t c x c, x Jgh reduciéndose la expresión (3) a la ecuación de Helmholtz para profundidad variable: El problema descrito es elíptico", por lo que deben especificarse condiciones de borde en el contorno I' que limita la región R donde rige el mismo, las cuales son: - Condición de Dirichlet: cuando es considerado el valor de la sobreelevación por, por ejemplo, mediciones in situ. - Condición de Cauchy: dentro de este tipo pueden incluirse la condición de reflexión perfecta: con n la normal al contorno reflectante rn, la condición de velocidades normales al contorno rq prescritas (Neurnann): donde q es el valor prescrito, y la condición de absorción (total o parcial): si se introduce un coeficiente de absorción a definido como energía absorbida a = energía total La determinación práctica de a es complicada, debido a que los mecanismos físicos que pretende simular son, en realidad, altamente no lineales (como, por ejemplo, fenómenos de rotura de oleaje). Por lo tanto su uso se restringe a la obtención de valores medios en una zona dada, o como filtro de ruidos. Debe observarse, además, que la expresión (5) sólo se cumple cuando la incidencia es normal al contorno. Por Último, hay que tener en cuenta que parte del dominio Cl es no definido o, a fines de su estudio, infinito. En esta zona debe imponerse la condición de que las ondas radiadas al exterior no retornen, tendiendo a anularse cuando la distancia tiende a infinito. La condición de radiación, propuesta por Sommerfeld (1949)" para casos periódicos, se puede expresar como: MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS donde r es la distancia medida a partir de un punto fijo y 4 es solución a la ecuación de la onda, mientras n es la dimensión del problema. Esta condición supone una cantidad positiva de flujo saliente de energía en un período TI1, representando la disipación energética más importante en ondas de baja frecuencia. En la siguiente sección serán tratadas las diferentes condiciones vistas en forma unificada, introduciendo conceptos del método de las características. CONDICIONES DE CONTORNO Si en la ecuación de ondas (3) se elimina la condición de periodicidad, considerando, por simplicidad, una dimensión espacial, ésta quedará expresada como siendo una ecuación de carácter hiperbólico de segundo orden. De esta forma, es posible obtener la pendiente de la curva característica correspondiente, dada por Si se considera la hipótesis de que la profundidad se mantiene constante en un entorno a la línea característica, se puede llegar a que a lo largo de ella se cumplirá 84 - * ~~"93 = cte at ax En el caso de profundidades reducidas (6) se transforma en Al introducir 71 y 72, que corresponden a cantidades que se propagan en direcciones opuestas, la expresión (6) puede escribirse también como El principio aplicado para determinar las condiciones de contorno es el de tratar las cantidades descritas prescribiendo o no su valor, obteniendo, al combinarlas, las diferentes condiciones físicas que se desean simular3'. De esta forma, sólo pueden ser prescritas aquellas magnitudes "entrantes", es decir, cuya característica penetra en el dominio. Las ecuaciones (7) pueden expresarse en función de los invariantes 71 y 72 como: P. ORTIZ Y M. PASTOR Si se prescribe la cantidad TI, que ingresa al recinto, como una onda progresiva plana 4~ dada por: q51(z, y, t) = A exp(i(k T - wt)) siendo r la dirección de incidencia, al escoger cualquiera de las expresiones (7) u (8) y teniendo presente la condición de periodicidad, se obtiene que que es la expresión, en el caso unidimensional, que define el contorno "cargado", es decir, el que conecta el dominio en estudio con el indefinido. Luego de algunas operaciones, (9) queda finalmente como válida para cualquier relación longitud de onda-profundidad. Si se trata de ondas largas, la ecuación (10) se convierte en Para obtener una extensión a dos dimensiones, se asumirá el cumplimiento de la condición (9) en la dirección de incidencia. En este caso, la dirección de incidencia no coincide, en general, con la dirección R, normal al contorno, por lo que se introduce p como el ángulo formado entre estas direcciones, quedando, en consecuencia, a partir de (9): 1 84 6--- + iwqi = Jccs.i.kq51 + iw41 cos p dn La expresión final es: y en aguas poco profundas: La condición de absorción es simple de obtener mediante los mismos argumentos, con la simplificación de admitir la incidencia normal al contorno. Para absorción total, se cumple que: 72: libre (saliente), TI = 0, por lo que, reemplazando en la segunda expresión de (8): MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS luego, en la primera de ellas, queda Esta condición, con incidencia normal, es válida para el problema bidimensional, o sea: luego, si se introduce la condición de periodicidad, se llega finalmente a para absorción total normal que, en el caso equivalente a absorción parcial es: La condición de Neumann (4) se puede concebir mediante la prescripción de las componentes de velocidad de (7), que corresponden a la onda que se aleja, lo que no se contradice con lo expuesto, ya que este caso es equivalente a conocer el comportamiento de una de las variables en la región exterior. Para reflexión total, su obtención es inmediata, ya sea anulando el coeficiente de absorción a en (12), o haciendo TI = 72 y 72 =libre en (8), considerando velocidades normales al contorno. FORMULACION NUMERICA La ecuación (3), junto con las condiciones de borde cargado (ll), de absorción (12), de Neumann (4) y de Dirichlet, han sido discretizadas mediante la técnica de BoubnovGalerkin. La formulación débil del problema puede expresarse como: donde 7, E, (, p son coeficientes de proporcionalidad, R es el dominio definido, r,, r,, I'+, r, son, respectivamente, los contornos absorbente, cargado, con valores prescritos y con velocidades prescritas y 3 representa los valores prescritos de 4. P. ORTIZ Y M. PASTOR evitando el almacenamiento de A. Las estrategias posibles son: 1) Almacenar A ensamblada, efectuando el producto A. pr, en cada iteración, 2) Almacenar las ~natrices elementales sin ensamblar Ae y efectuar xfZ1 Ae pk en cada iteración y 3) Calcular Ae en cada iteración y efectuar el producto de la misma forma que en la opción 2). La primera Enea es la de menor tiempo de ejecución y máxima memoria o entradasalida, poseyendo escasas ventajas respecto a un método directo. La segunda alternativa es intermedia y la tercera es la de menor almacenamiento y mayor tiempo de proceso. Como resultado de experimentación numérica, pueden compararse los órdenes del número de coeficientes necesarios de almacenamiento en cada caso (sin contar con los datos topológicos), siendo N el número de grados de libertad: Caso 1): Ñ N*; Caso 2): Ñ ION; Caso 3): 0. Estos órdenes se obtienen suponiendo que en el caso 1) ha sido optimizada la numeración nodal, mientras que en el segundo se adopta una relación número de nodos/número de elementos de 3.5 (media de todos los casos estudiados), siendo guardada sólo la matriz triangular superior. Se ha comprobado además que para casos con grados de libertad entre 81 y 509 se obtienen tiempos de proceso de alrededor de 7 veces mayores en el caso 3) respecto al 2), lo que muestra que, salvo en problemas críticos de capacidad, la opción práctica es la segunda. Es conveniente añadir que si se compara cualquier método indirecto con uno directo, en el tiempo de ejecución de éste Último, estrictamente, debe tenerse en cuenta el tiempo invertido en la optirnización de frente o de banda. APLICACIONES El programa desarrollado en base a los conceptos previos emplea elementos de ocho nodos (cuadriláteros), de la familia serendíptica y funciones de forma de continuidad CO. Se ha optado por 16 puntos de integración por elemento y 4 en el caso de integrales curvilíneas. La solución del sistema proporciona sobreelevaciones y función potencial en nodos, así como agitaciones medias en el área de estudio y velocidades en puntos de Gauss. La excitación exterior es impuesta especificando el contorno en "mar abierto", r,, la frecuencia angular de la onda incidente, su dirección y amplitud máxima. Dentro de las distintas pruebas efectuadas al modelo, se incluyen en este trabajo 3 problemas típicos, presentándose posteriormente un caso real referente a puertos del Mar Cantábrico. Los problemas referidos son los puertos rectangulares, rompeolas circular y refracción por variación parabólica de profundidad. Puertos rectangulares El problema del puerto rectangular, estudiado en forma teórica, numérica y e~perimental"~~~~ es analizado en esta ocasión mediante dos configuraciones diferentes, representadas en las Figuras 1.a y l.b, siendo en ambos casos las paredes totalmente reflectantes. En la primera configuración, el ancho de la dársena coincide con el de la bocana, mientras que en la segunda la bocana está parcialmente cerrada por dos rompeolas rectos. MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS C/2d = 0.1 b Figura 1. Puertos Rectangulares. Para incluir correctamente la influencia de la costa recta lindante a la dársena, pueden considerarse dos opciones: i) considerarla totalmente reflejante o ii) totalmente absorbente. Si se elige la primera, se debe componer a la onda incidente la reflejada en la costa exterior. Para el primer ejemplo se empleó una malla de 88 elementos y 317 nodos, mientras que en el segundo se recurrió a 160 elementos y 537 nodos. A continuación se hizo incidir una onda plana de amplitud unitaria perpendicular a la costa recta, obteniendo resultados equivalentes en las opciones i) e ii) en ambas dársenas, siendo resumidos en la Figura 2. En esta Figura se representan los factores de amplificación en los puntos P representados en la Figura 1, definidos como la relación entre la amplificación máxima en el punto considerado y la amplificación de la onda incidente si no estuviese implantado el puerto. En el primer análisis se intentó reproducir el primer pico resonante, que es el de mayor interés, obteniendo resultados satisfactorios, superiores a los que se obtienen en modelo físico, al no contabilizar pérdidas locales y por fricción16. Para la segunda configuración, se estudiaron los cuatro primeros picos de la curva de amplificación en función de la longitud de onda relativa, comparados con valores obtenidos en modelo físico16, siendo esta vez los valores muy cercanos entre ambas P. ORTIZ Y M PASTOR 20 18 16 u 14 L t'2 U 2 10 A -8I -62Ili7i"'1 - - - - * - F :: - i :, ,' \ ,,-\ - .-N .-?--- ,*' '- -____ .' -- --_ _ _-- I'(lI 00~234567a910 K L Caso A: - - - : [16]. o: Modelo nuiiiérico. 1 o 8 L u 26 18 - u - - U Q - LL UI 012345678013 Caso B: -: [16]. m: Modelo numérico. Figura 2. Puertos Rectangulares. Factores de amplificación vs. longitudes de ondas relativas kl. técnicas. Rompeolas circular Esta prueba consiste en una planta ideal formada por una costa recta infinita donde, adyacente a ella, se emplaza un rompeolas semicircular con dos aberturas, haciendo incidir sobre la región por él limitada una onda unidad paralela a la costa. La malla utilizada (Figura 3) consta de 24 elementos y 79 nodos siendo la profundidad constante. Figura 3. Rompeolas Circular. r,: Contorno "cargadon. r,: Contorno absorbente. MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS Este test, al ser usado en el modelo de elementos infinitoss, exigía contemplar una región exterior como transición al ensamble de elementos infinitos, lo que implica mayor número de grados de libertad. Los resultados, comparados con el estudio citado y con la solución analítica se exponen en la Figura 4, donde se reproduce la respuesta en función de la longitud de onda relativa kr, que esta vez relaciona la longitud de onda y el radio de la planta semicircular. Por otra parte, la respuesta de la dársena es evaluada mediante: siendo A el área total encerrada por el rompeolas. Luego, en la Figura 4 se representa el cociente entre esta sobreelevación "media" y la de la bocana orientada hacia la dirección de incidencia de la onda. O 0.4 0.8 1.2 1.6 2 2.4 2.8 L onda rel. KR Figura 4. Rompeolas Circular. Respuesta. -: Solución Analítica. o: Elementos Infinitos. O: Modelo Actual. La coincidencia con los valores analíticos es buena, incluso en zonas alejadas del segundo pico, a pesar de lo gruesa de la malla. La distribución de sobreelevaciones también posee un buen ajuste, como puede verse en la Figura 5 (kr = 3.0). Por las características geométricas que presenta este ejemplo, la influencia que puede tener el optar por un contorno absorbente recto en la bocana "no cargada" es poca, lo que es confirmado por los buenos resultados obtenidos. Sin embargo, el uso de la condición de absorción debe ser hecho con cautela, dada sus limitaciones. Cuando la batimetría de la zona costera permite un efecto de refracción importante, la absorción normal puede ser admisible, al tender los frentes de onda a ser paralelos a las líneas batimétricas. P. ORTIZ Y M. PASTOR Figura 5. Rompeolas Circular. Factores de Amplificación. Ba = 3, a: radio. -: Analítico [8]. - - - : Modelo Numérico. Caso A. Ondas Largas. Caso B. Ondas Cortas. Figura 6. Variación parabólica de profundidad. Refracción por variación parabólica de profundidad. Para contemplar un caso donde pueda comprobarse la efectividad de la aproximación para variaciones importantes de profundidad, se ha estudiado la refracción de una onda plana que se propaga por una zona circular donde la profundidad varía en forma parabólica (Figura 6). Se han efectuado dos casos, representados como A y B en la ~igura 6, correspondiendo el primero a un problema en que no es respetada la condición de MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS variación gradual de la profundidad, mientras que en el segundo esta limitación no es transgredida. En el caso A (Figura 7), los resultados, que se presentan como sobreelevaciones en la dirección de incidencia, no son muy ajustados, siendo similares a los obtenidos mediante el uso de elementos finitos8. DlRECClON DE ONDA INCIDENTE Figura 7. Variación parabólica de profundidad. Amplitud relativa. Caso A. -: Analíticase. o: Elementos Infinitos'. o: Modelo actual. El número de elementos empleados en esta ocasión fue de 168 y el de nodos 505. El problema B (Figura 8) (ondas cortas) ha necesitado de 576 elementos y 1729 nodos, dada la reducida longitud de onda. Los resultados, esta vez, presentan una aproximación aceptable siendo nuevamente similares al modelo en elementos infinitos, con un ahorro de un 10% de grados de libertad respecto al mismo. dist. del centro (cm.) Figura 8. Variación parabólica de profundidad. Amplitud relativa. Caso B. - o - Analíticase. +: Elementos Infinitos8. O: Modelo actual. P. ORTIZ Y M. PASTOR Mar Cantábrico Los mecanismos de generación de ondas de períodos en la banda de 1 a 30 minutos no están aún completamente determinados, al contrario de otras ondas largas como por ejemplo tsunamis. La hipótesis antes citada de transferencia de energía desde frecuencias altas podría interpretarse como la modificación de un tren inicial de ondas con poca energía en la zona de baja frecuencia que, al propagarse hacia la costa, con variaciones de profundidad importantes y extensos desarrollos, muestra un pico en la distribución de frecuencias en una región de valores menores. El pico de la nueva distribución corresponde al período dominante de la onda larga. En el Mar Cantábrico, donde se estima que las ondas son producidas por mecanismos de este tipo, presentan un período preponderante cercano a los 4 minutos. Este período, relacionado con el de los grupos de ondas del espectro inicial, puede ser estimado inicialmente como el del asociado al grupo de ondas. Una metodología de estudio de ondas largas en un recinto portuario puede ser, en primer lugar, la obtención de frecuencias propias principales del sistema no a~nortiguado~~, detectando así las zonas de máximas amplitudes de movimientos horizontales y verticales3'. El primer caso permite decidir emplazamientos de medición de ondas largas, mientras el segundo indica sectores de atraque desfavorables. Los modos propios pueden calcularse mediante la ecuación de autovalores generalizada, oytándose por las frecuencias más bajas. Una forma alternativa de obtener la frecuencia propia fundamental de un recinto semicerrado de agua es mediante la ecuación lineal de ondas largaP, excitando el sistema con una función armónica e imponiendo condiciones iniciales nulas. En la respuesta aparecerán oscilaciones espúreas, mostrando un máximo de energia en la frecuencia propia fundamental del sistema3'. Al tratar la respuesta mediante la transformada discreta de Fourier, en el contenido en frecuencias surge el pico citado y, en consecuencia, la frecuencia buscada. La influencia del "mar abierto" en el modelo simple no amortiguado puede asumirse mediante dos opciones: reflexión perfecta en la bocana (lo que equivale a calcular los períodos propios de un cuerpo de agua cerrado), o suponer amplitud vertical nula en la bocana y, por lo tanto, una línea nodal en ella. Aunque la realidad es una situación intermedia, la hipótesis más razonable es la segunda, confirmada al compararse con modelos físicos32. La parte final de un estudio es la aplicación del modelo disipativo descrito en este trabajo, que permite evaluar respuestas cuantitativamente más acordes con la realidad, limitando su uso a la zona del espectro acotada por el estudio previo. Dentro de esta Última etapa, se incluyen en este trabajo dos puertos pesqueros dC la costa Vasca: Motrico y Ondárroa, cuyo comportamiento fue analizado junto con el de puertos de dimensiones mayoresas. El interés de estos pequeños puertos reside en su comportamiento totalmente opuesto frente a una excitación similar. Los pequeños puertos pesqueros de bajura del Mar Cantábrico, tales como los nombrados, poseen en algunos casos características geométricas que ocasionan problemas resonantes ante las excitaciones de períodos cercanos a 4 minutos. Los primeros estudios sistemáticos de estos fenómenos se deben a Iribarren17*18, MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS cuya explicación intuitiva del fenómeno es validada por el modelo numérico. El puerto de Motrico es el que presentaba este efecto en forma más pronunciada, pues era suficiente cualquier marejada exterior para que la agitación interior fuese elevada. Así, en la Figura 9 se representan las líneas de igual elevación para una excitación unidad en la situación resonante (245 seg.). Puede observarse el elevado gradiente en la zona de fondeo, indicada con A. En el modelo reducido" se define el factor de amplificación como el cociente entre velocidades en A y en la zona exterior, obteniéndose una curva de amplificación (Figura 10) poco pronunciada. Esto puede deberse a la "ambigua" definición del coeficiente de amplificación, al ser las velocidades exteriores elegidas poco representativas de la excitación. Con el presente modelo se ha optado por definir el factor comparando sobreelevaciones en el extremo de la dársena interior con amplitud de la excitación, otorgando la respuesta un pico más definido. Puede observarse que la coincidencia de los picos entre modelo físico y numérico es buena. SEP. ISOL.= 0.5 Figura 9. Puerto de Motrico. Líneas de igual elevación. T = 245 s. A: Zona de fondeo. B: Zona exterior. C: Zona eliminada. Para el modelo se empleó una malla de 345 nodos y 92 elementos, habiéndose incluido una zona exterior B (Figura 9) relativamente pequeña. Los problemas fueron parcialmente resueltos eliminando la zona C (Figura 9), y dragando, es decir modificando la geometría del recinto. A pocos kilómetros al oeste de Motrico se encuentra Ondárroa, cuyo puerto consta de un antepuerto y una segunda dársena conectados por una bocana interior. Los problemas del puerto de Motrico exigía en muchas ocasiones que los pesqueros allí fondeados tuviesen que trasladarse a Ondárroa, con las consiguientes pérdidas económicas y de capacidad de ambos puertos17. Pero ¿qué razón producía la calma en el puerto vecino al de Motrico cuando éste era impracticable para las actividades de los pesqueros?. La contestación a esta situación surge de la curva de respuestas de P. ORTIZ Y M. PASTOR Figura 10. Puerto de Motrico. Factores de amplificación. o : Amplificación de velocidades17. +: Modelo numérico (elevaciones). frecuencia angular Figura 11. Puerto de Ondárroa. Curva de respuesta. Ondárroa (Figura 11): para la banda crítica de 230 a 250 seg. (0.027 a 0.025 rad/seg. de frecuencia angular) la amplificación es mínima, hecho detectado visualmente por los pescadores. En la Figura 11 se representa el barrido en frecuencias realizado para el intervalo 170 s.-600 s. El factor de amplificación ha sido calculado por el cociente entre sobreelevaciones en el nodo de control P (Figura 12), donde se producen máximas elevaciones para los dos primeros picos resonantes y la onda incidente unidad exterior. En la misma Figura 12 se representan las líneas de igual elevación para T = 245 s., donde puede observarse las condiciones favorables del puerto para esta frecuencia. MODELO NUMERICO DE REFRACCION-DIFRACCION DE ONDAS SEP. ISOL.=0.5 Figura 12. Puerto de Ondárroa. Líneas de igual elevación. T = 245 s. A: Zona exterior de playas. Para realizar el estudio se construyó una malla de 509 nodos y 142 elementos, contemplando en la discretización la zona exterior de playas (zona A en la Figura 12) con contornos reflejantes a onda larga, dado que la absorción a éstas es muy baja. CONCLUSIONES El tratamiento numérico de la ecuación aproximada de difracción y refracción exige resolver el problema '(exterior" adecuadamente, lo que es propuesto mediante diversas vías, de las cuales puede destacarse los modelos en elementos infinitos. En este artículo se plantea una forma diferente de encarar el tema, en base al método de las características. A efectos de la eficiencia resolutiva del modelo, se obtiene menor número de grados de libertad que en caso de usar elementos infinitos, así como un mejor condicionamiento de la matriz de coeficientes que conserva, a su vez, el carácter bandeado simétrico. El modelo se presenta como un serio competidor a los de similar potencia, poseyendo similares limitaciones, derivadas de la forma de imponer el comportamiento de la onda radiada. A este respecto, una variante de las condiciones aquí derivadas pueden obtenerse en base a la dirección normal al contorno "cargado", habiéndose conseguido resultados más pobres y más dependientes de la forma del contorno. Las limitaciones de la condición de absorción normal implica su uso a zonas reducidas del dominio. Si el problema se formula en su forma hiperbólica y se plantea un esquema explícito de solución, una vía de mejorar el problema de absorción es la de corregir el módulo y dirección de la velocidad que incide en el contorno absorbente en cada salto de tiempo. La aplicación de estos conceptos en el caso armónico lleva a un problema no lineal y, por lo tanto, mucho más costoso. Para salvar el problema del elevado número de grados de libertad presente en ciertos