Modelo de Elementos Finitos No-Hidrostático Integrado en la Vertical para la Propagación de Olas

DEPTHINTEGRATED NON HYDROSTATIC FINITE ELEMENT MODEL FOR WAVE PROPAGATION

Lucas E. Calvo Gobbetti^1*

Centro de Investigaciones Hidráulicas e Hidrotécnicas
Universidad Tecnológica de Panamá

[email protected]

Paulo C. Colonna Rosman²

COPPE-PEnO, Área de Engenharia Costeira e Oceanográfica
Universidade Federal do Rio de Janeiro

[email protected]

Resumen–Un modelo de elementos finitos no hidrostático integrado en la profundidad para la propagación y transformación del oleaje en zonas costeras fue desarrollado con suceso a partir del modelo hidrostático SisBahia®. El modelo utiliza elementos finitos cuadrilaterales cuadráticos para la aproximación de las velocidades horizontales y elementos finitos cuadrilaterales lineales para la aproximación de las elevaciones de la superficie del agua y las presiones no hidrostáticas. El modelo es verificado con una solución analítica y validado usando datos experimentales. Al no requerir del uso de mallas intercaladas el presente modelo puede ser usado en mallas no estructuradas de elementos finitos.

Palabras Claves–Propagación de olas, no hidrostático, método de elementos finitos, SisBahia®.

Abstract–A depth integrated non hydrostatic finite element model for the propagation and transformation of waves in coastal areas was developed with success from the hydrostatic model SisBahia®. The model use quadratic quadrilateral finite element for horizontal velocities approximation and linear quadrilateral finite elements for water surface elevations and non hydrostatic pressures approximations. Because the model does not require the use of staggered grids it can be used on non structured finite element meshes.

Keywords–Waves propagation, non hydrostatic, finite element method, SisBahia®.

1. Introducción

Las ecuaciones de aguas poco profundas, con aproximación hidrostática de presiones, tienen muchas aplicaciones en la modelación de las ondas largas, tales como mareas y olas de tormenta. Sin embargo, la falta de dispersión (conservación del momento vertical) en los modelos de aguas poco profundas impide su uso en la modelación de la propagación de las olas desde mar adentro hasta la costa.

Casulli y Stelling [1] y Stansby y Zhou [2] propusieron independientemente simular el movimiento de las olas usando las ecuaciones de Navier Stokes con presión no hidrostática. Stelling y Zijlema [3] formularon un modelo de diferencias finitas que toma en cuenta la dispersión a través de un término de presión no hidrostática. Tanto en versiones integradas en la profundidad como en formulaciones de varias capas, ellos descomponen la presión en componentes hidrostáticas y no hidrostáticas, como en [4]. Yamazaki et al. [5] propusieron un modelo no hidrostático integrado en la profundidad de diferencias finitas con un esquema advectivo de conservación del momento que permite la modelación del quiebre de las olas y del run up.

Las ecuaciones integradas en la profundidad son análogas a las ecuaciones de aguas poco profundas con la adición de una ecuación de momento vertical y un término de presión no hidrostático en las ecuaciones del momento horizontal. Walters [6] adaptó este enfoque no hidrostático a un modelo mixto de elementos finitos y volúmenes finitos. Wei y Jia [7] desarrollaron sobre un modelo existente, el CCHE2D, una versión no hidrostática donde mallas intercaladas de elementos finitos son usadas para la determinación de las derivadas espaciales.

En este trabajo un modelo de elementos finitos no hidrostático integrado en la profundidad para la propagación de olas es desarrollado partiendo de un modelo hidrostático existente, el SisBahia® [8]. SisBahia® (Sistema Base de Hidrodinámica Ambiental) es un sistema profesional de modelos computacionales registrado por la Fundación Coppetec, ligada al Instituto Alberto Luiz Coimbra de Pós-Graduação e Pesquisa de Engenharia (COPPE) de la Universidade Federal do Rio de Janeiro (UFRJ). El modelo no hidrostático es verificado con una solución analítica y validado usando datos experimentales.

2. Modelo Matemático

Los modelos no hidrostáticos para la dinámica de olas son generalmente desarrollados basados en sus predecesores hidrostáticos. El modelo hidrostático SisBahia® viene siendo utilizado ampliamente desde 1987 y actualmente incluye módulos hidrodinámicos 2DH y 3D, de calidad de agua, transporte de sedimentos, transporte Euleriano y Lagrangeano y generación de oleaje por viento, entre otros. Las ecuaciones gobernantes en el modelo hidrodinámico no hidrostático integrado en la profundidad, sin incluir los términos turbulentos y de Coriolis, son:

Donde U, V y W son componentes del promedio de las velocidades en las direcciones x, y, y z, p es la densidad del agua,Ch es el coeficiente de Chezy, tradicionalmente utilizado para parametrizar la fricción de fondo, yg es la aceleración de la gravedad. La profundidad de flujo se define como H = + h, donde es el nivel de la superficie del agua medido desde un nivel de referencia fijo y h es la profundidad medida desde ese mismo nivel. En estas ecuaciones se asume una distribución lineal en la vertical tanto para las presiones no hidrostáticas como para las velocidades verticales. La presión no hidrostática en la superficie libre es tomada como cero y en el fondo como siendo qb. La velocidad vertical promedio, W, es (w( + wb)/2, donde w( es la velocidad vertical en la superficie y wb es la velocidad vertical en el fondo.

Las condiciones de frontera cinemáticas de la superficie libre y del fondo son:

Con Uξ, Vξ, Ub y Vb siendo las componentes en x y y de las velocidades junto a la superficie libre y el fondo.

3. Formulación de Elementos Finitos y Solución Numérica

3.1 Elementos Finitos

La discretización espacial en el modelo hidrostático SisBahia® se efectúa utilizando elementos cuadrilaterales cuadráticos (9 nodos) tanto para las velocidades como para los niveles de la superficie del agua. Para garantizar que los sistemas lineales resultantes de la formulación numérica del modelo no hidrostático estén bien planteados, es necesario usar elementos cuadrilaterales mixtos: cuadráticos (9 nodos) para las velocidades y lineales (4 nodos) para las elevaciones de la superficie del agua y las presiones no hidrostáticas (figura 1).

Figura 1. Malla y elementos finitos cuadrilaterales usados en el modelo: Cuadráticos (9 nodos, círculos rellenos) para velocidades; Lineales (4 nodos, círculos) para elevación de la superficie y presión no hidrostática.

Las velocidades U, V, el nivel de la superficie del agua ξ y la presión no hidrostática en el fondo son aproximados en un elemento por:

Donde son, respectivamente, las funciones de aproximación cuadráticas (9 nodos) y lineales (4 nodos) de las variables de solución.

3.2 Solución Numérica

3.2.1 Primer Paso

El procedimiento de solución comienza resolviendo las siguientes ecuaciones del momento horizontal (1) y (2), sin los términos de presión no hidrostática, de forma acoplada hasta encontrar las soluciones provisionales para el paso de tiempo n+1, Ũ^n+1,Ṽⁿ⁺¹:

Las ecuaciones acopladas (9) y (10) son tratadas con el esquema numérico del SisBahia®, basado en el método de los elementos finitos de Galerkin [8]. El esquema de discretización temporal adoptado en el modelo SisBahia® es de 2ª orden, centrado en el instante t +Δ/2, y usa la formulación de Crank-Nicholson para los términos lineales y la factorización implícita para los términos no lineales. Suponiendo la siguiente forma general:

Con L representando una función lineal y L1 L2 una función no lineal que puede ser escrita como el producto de dos funciones lineales, el esquema de discretización temporal adoptado sería:

En la documentación del modelo SisBahia® se utiliza la siguiente notación:

El valor extrapolado en el instante t +Δt es calculado por:

Un modelo de elementos finitos de Galerkin centrado en el instante t +Δ/2 para resolver las ecuaciones del momento horizontal (9) y (10) acopladas sobre un elemento Ωe es:

Las integrales en el domino del elemento Ωe son definidas como:

Para controlar las oscilaciones espurias provenientes de los términos advectivos, el modelo SisBahia® usa la interfase anti-disipativa explícita desarrollada por Rosman [8]. Esta interfase alisa las oscilaciones modificando los valores calculados de las variables usando valores calculados de puntos adyacentes, sin provocar una difusión excesiva de la solución. En el caso de mallas no estructuradas, con muchos elementos compartiendo el mismo nodo, el resultado de la aplicación de la interfase depende de cuales de estos elementos sean considerados.

En el presente modelo no-hidrostático las oscilaciones espurias de los términos advectivos son tratadas con el ya tradicional método Streamline Upwind Petrov-Galerkin (SUPG), cuya aplicación no depende del tipo de malla utilizada. En el método SUPG [9] son incluidas en las ecuaciones del momento términos estabilizadores, que en el caso del modelo en (15) son dados por:

con:

Luego de ensamblados todos los elementos del dominio se resuelve el sistema lineal resultante para conseguir las velocidades horizontales provisionales Ũ^n+1,Ṽⁿ⁺¹.

3.2.2 Segundo Paso

En el segundo paso se construye una ecuación de Poisson que se resuelve implícitamente para obtener las presiones no hidrostáticas. Una aproximación de la ecuación del momento vertical en (3) es:

La velocidad vertical en el fondo se estima a partir de la condición de borde cinemática (6) como:

Las velocidades horizontales finales, influenciadas por los términos de presión no hidrostática, pueden ser expresadas de la siguiente forma:

Para obtener una solución correcta entre el campo de velocidades y las presiones no hidrostáticas, la ecuación de la continuidad se aplica directamente sobre la columna de agua:

Al substituir las ecuaciones (28) a (32) en la ecuación (33) se establecen una ecuación de Poisson dela cual se obtiene la presión hidrostática en el fondo q_bⁿ⁺¹.

Un modelo de elementos finitos de Galerkin para resolver la ecuación de Poisson sobre un elemento Ωe es:

Con las integrales en el dominio del elemento Ωe definidas por:

Las integrales sobre el contorno del elemento гe pueden ser escritas como:

Con nx y ny siendo las componentes en x y y del vector normal unitario en el nodo i y qb un valor de qb extrapolado para el paso de tiempo n+1. Una vez ensamblados todos los elementos del dominio se resuelve el sistema lineal resultante para obtener las presiones no hidrostáticas q_bⁿ⁺¹.

3.2.3 Tercer Paso

Finalmente, una vez conocidas las presiones no hidrostáticas q_bⁿ⁺¹ y las soluciones provisionales Ũⁿ⁺¹ y Ṽⁿ⁺¹,se resuelven las siguientes ecuaciones a través del método de los elementos finitos de Galerkin, para obtener las soluciones finales Ũⁿ⁺¹,Ṽⁿ⁺¹,ξⁿ⁺¹.

Un modelo de elementos finitos de Galerkin para resolver de forma acoplada las ecuaciones (43) a (45) sobre un elemento Ωe es:

Las integrales en el dominio del elemento Ωe son definidas como:

Luego de ensamblados los elementos se resuelve el sistema para alcanzar las soluciones finales. El tiempo de cálculo del modelo es aproximadamente 1.5 veces el tiempo requerido para el modelo hidrostático.

4. Verificación y Validación del Modelo

En esta sección se verifica y valida la exactitud y precisión del modelo de ola no hidrostático integrado en la profundidad usando tres casos de referencia. El primer caso es la propagación de una onda solitaria a lo largo de un canal de profundidad de agua constante. Debido a que la onda solitaria no es una solución de las ecuaciones de aguas poco profundas, este caso verifica la correcta propiedad no hidrostática del modelo. El segundo caso trata de la propagación de olas regulares sobre una barra sumergida, validando la capacidad del modelo para el manejo de la no linealidad de las olas. El tercer caso es el estudio del run up de una onda solitaria sobre una pared vertical.

4.1 Propagación de una onda solitaria a lo largo de un canal de profundidad de agua constante

La onda solitaria es una onda no lineal con amplitud finita; si el fluido es no viscoso, y el fondo horizontal es sin fricción, la onda debe mantener su forma y su velocidad durante todo el proceso de propagación. Este caso se ha utilizado para verificar muchos modelos no hidrostáticos [3] [10] [6] [5].

En este estudio, consideramos un canal sin fricción de 500 m de largo y profundidad del agua h de 10 m. Una condición de radiación de Sommerfeld basada en la elevación de la superficie del agua se especifica en la entrada y en la salida del canal. La onda solitaria se encuentra inicialmente en X₀= 150 m con una altura de ola A = 2m. Los tamaños de malla son Δx=Δy=1 y el paso de tiempo es Δt= 0.025 seg. El número de Courant en cuanto a la celeridad de la ola es C_r= ≈ 0.25, pero números de Courant > 1 pueden ser soportados. Las soluciones analíticas de la elevación de la superficie libre del agua y la velocidad son dadas por:

Donde la celeridad de la onda es C =y el número de onnda k =.

En la figura 2, se presenta la onda solitaria inicial y las formas de onda analíticas y simuladas a lo largo del canal a los 5, 10 y 15 seg. Existe una ligera reducción de la altura de las ondas en el comienzo de la simulación debido a la condición inicial aproximada por la solución analítica; fenómenos similares fueron reportados por Walters [6] y Yamazaki et al. [5]. Aunque pequeñas ondas residuales fueron observadas, la forma de la ola y su amplitud están bien conservadas durante la simulación.

4.2 Propagación de olas sinusoidales sobre una barra sumergida

Beji y Battjes [11] y Luth et al. [12] llevaron a cabo experimentos físicos acerca de la propagación de olas regulares sobre una barra trapezoidal sumergida en un canal de 37.7 m de largo, 0.8 m de ancho y 0.75 m de profundidad.

Figura 2. Propagación de ondas simuladas y solución analítica para diferentes tiempos de simulación.

La figura 3 muestra la configuración numérica del experimento. La profundidad del agua quieta es 0.4 m. Una barra trapezoidal de 0.3 m de alto, con una pendiente frontal de 1:20 y una pendiente posterior de 1:10, está situada entre los 6 y los 17 m en el canal. Olas sinusoidales incidentes con una amplitud de 1.0 cm y un período de 2.02 seg, correspondientes al parámetrokH � 0.67, se generan en el lado izquierdo en base a la teoría de las olas lineal. La frontera del experimento en el lado derecho se modela como un área de flujo abierto donde se impone una condición de radiación de Sommerfeld. La presión no hidrostáticaqb se asume cero en los lados izquierdo y derecho El dominio computacional de 35 m de largo de la figura es discretizado con una malla de elementos quadrilaterales de lados ∆x = ∆y = 1.25 cm. Al inicio de la simulación la velocidad y la elevación de la superficie del agua es 0 en todos los nodos y un intervalo pequeño de tiempo ∆t= 0.002 seg se utiliza para desarrollar la ola suavemente.

Figura 3. Configuración numérica para la propagación de olas sinusoidales sobre una barra sumergida (Fuente: Yamazaki et al. [5]).

Comparaciones de las elevaciones simuladas y medidas de la superficie del agua se presentan en las figuras 4-a (puntos 4 a 6) y 4-b (puntos 7 a 9). Se observa que el proceso de asomeramiento en el lado frontal de la barra (puntos 4 y 5) y la transformación de las olas de una zona de dispersión de baja frecuencia (punto 6) a una zona de dispersión de alta frecuencia (punto 8) son correctamente predichas por el modelo. Sin embargo, discrepancias notables aparecen entre las elevaciones de la superficie libre simuladas y medidas detrás de la barra (punto 9), donde se liberan los armónicos más altos. Un análisis espectral [13] ha demostrado que el parámetrokH oscila desde 6 hasta 10 en esta zona; estas olas altamente dispersivas están fuera del rango de aplicación de los modelos no hidrostáticos integrados en la profundidad [6] [5]. Un modelo no hidrostático con múltiples capas podría mejorar las simulaciones en el punto 9.

Figura 4-a. Comparación de elevaciones de la superficie simuladas y medidas. Simuladas (líneas sólidas), datos de Beji y Battjes [11] (círculos).

Figura 4-b. Comparación de elevaciones de la superficie simuladas y medidas. Simuladas (líneas sólidas), datos de Beji y Battjes [11] (círculos).

4.3 Run Up de una onda solitaria sobre una pared vertical

Para estudiar los parámetros físicos envueltos en el run up de una onda de tsunami, fueron realizados una serie de experimentos a escala de ondas de tsunami financiados por la National Science Foundation [14] [15]. Estos experimentos se han convertido en testes padrón para la validación de modelos numéricos de ondas [16] [17]. Uno de estos experimentos fue simulado para evaluar el desempeño del modelo no hidrostático desarrollado.

En la figura 5 se muestra el esquema del experimento del run up de una onda solitaria sobre una pared vertical. El canal tiene 23.23 m de largo y 0.45 m de ancho, mientras que las pendientes compuestas sobre la batimetría fija del fondo consisten de tres diferentes pendientes (1:53, 1:150 y 1:13) y una sección plana al final del canal. La pared vertical está localizada al final de la pendiente de 1:13. La profundidad del agua en reposo en la sección plana del canal es h = 21.8 cm. Para la medición de la elevación de la superficie del agua fueron usados 10 medidores de ondas de capacitancia ubicados a lo largo de la línea central del canal. Se consideraron dos casos, denominados como I y II, con alturas de onda relativa de A/h = 0.039 y 0.246 respectivamente. El Caso I corresponde a una onda pequeña que no crece lo suficiente como para quebrarse antes de llegar a la pared. El Caso II es una onda mayor que quiebra justo antes de golpear la pared vertical. Este caso ha sido usado como test por muchos modelos numéricos hidrostáticos [15]. Debido a la falta de dispersión, los modelos hidrostáticos no lineales predijeron frentes de onda más empinados que los medidos, llevando eventualmente a un quiebre de onda prematuro. Adicionalmente, diferencias de fase entre los resultados numéricos y los medidos fueron también observados [18]. En este estudio, la resolución de la malla es determinada por la longitud de la onda para asegurar que la onda sea aproximada por lo menos por cien elementos; debido a ello se utilizaron mallas de elementos cuadrilaterales de 0.05 y 0.02 m de lado para los Casos I y II respectivamente. Correspondientemente, el intervalo de tiempo utilizado en los cálculos fue de ∆t = 0.01 y 0.005 segundos para los Casos I y II respectivamente.

Figura 5. Esquema del experimento del run up de una onda solitaria sobre una pared vertical.

Las figuras 6 y 7 muestran las comparaciones entre los perfiles de agua simulados y medidos. Los datos medidos en el punto 5 son usados como referencia para el ajuste temporal de las formas de onda simuladas. En el Caso I (figura 6), con una altura de onda relativa pequeña, se observa una buena coincidencia entre los valores simulados y medidos tanto para la onda entrando como saliendo luego de reflejarse en la pared vertical. Pequeñas diferencias con respecto al nivel 0 (≈2 mm) de los valores medidos en los puntos 7, 8 y 9 pueden ser debidos al montaje del experimento [19]. En el Caso II (figura 7) la onda se quiebra antes de llegar a la pared y el modelo numérico fue capaz de continuar sin el uso de parámetros empíricos, obteniendo predicciones satisfactorias.

Figura 6. Comparación de elevaciones de las superficies simuladas y medidas (Caso I). Simuladas (líneas sólidas), datos experimentales (círculos).

Figura 7. Comparación de elevaciones de las superficies simuladas y medidas (Caso II). Simuladas (líneas sólidas), datos experimentales (círculos).

5.Conclusiones

El desarrollo de modelos no-hidrostáticos tiene aplicación potencial en áreas de interés como: energía oceánica, ambiente, infraestructura costera, educación e investigación.

Un modelo de elementos finitos no-hidrostático integrado en la profundidad para la propagación del oleaje en zonas costeras fue desarrollado a partir del modelo hidrostático SisBahia®. Varias pruebas estandarizadas fueron llevadas a cabo para verificar y validar el modelo no hidrostático. La inclusión de un módulo de oleaje dinámico en el sistema SisBahia® será de bastante utilidad, ayudando a difundir aún más su uso y a propiciar su desarrollo.

El modelo no hidrostático fue capaz de mantener la forma y la velocidad de propagación de una onda solitaria a lo largo de un canal de profundidad de agua constante. Esto es debido a que las distribuciones verticales de presión y velocidad del modelo no hidrostático son equivalentes a las de las ecuaciones de Boussinesq clásicas, una solución de las cuales son precisamente las ecuaciones de la onda solitaria (59) y (60).

También, el modelo consiguió describir la transformación de olas sinosoidales sobre una barra sumergida. Específicamente, se reprodujeron con éxito los procesos de asomeramiento y transformación de las olas desde zonas de dispersión de baja frecuencia a zonas de dispersión de alta frecuencia.

Finalmente, el modelo fue capaz de describir el run up y la reflexión de una onda solitaria sobre una pared vertical.

Al no requerir el uso de mallas intercaladas, a diferencia de los modelos en Wei y Jia [7] y en Wei y Jia [19], el presente modelo puede ser usado en mallas no estructuradas de elementos finitos.

En el futuro, este modelo de ola será reforzado mediante un algoritmo de mojado y secado para la simulación del run up de la onda rompiente sobre superficies inclinadas, las más encontradas en la práctica, y una técnica de múltiples capas para el modelado de olas altamente dispersivas.

6. Referencias

Casulli V., Stelling G. S. “Numerical simulation of 3D quasi- hydrostatic, free-surface flows”. Journal of Hydraulic Engineering, Vol. 124, No. 7, pp. 678–686, 1998.

Stansby P. K., Zhou J. G. “Shallow-water flow solver with non-hydrostatic pressure: 2D vertical plane problems”. International Journal for Numerical Methods in Fluids, Vol. 28, No. 3, pp. 541–563, 1998.

Stelling G., Zijlema M. “An accurate and efficient finite difference algorithm for non-hydrostatic free surface flow with application to wave propagation”. International Journal for Numerical Methods in Fluids, Vol. 43, No. 1, pp. 1–23, 2003.

Casulli V. A. “Semi-implicit finite difference method for non- hydrostatic free surface flows”. International Journal for Numerical Methods in Fluids, Vol. 30, No. 4, pp. 425–440, 1999.

Yamazaki Y., Kowalik Z., Cheung K. F. “Depth-integrated, non-hydrostatic model for wave breaking and run-up”. International Journal for Numerical Methods in Fluids, Vol. 61, No. 5, pp. 473–497, 2008.

Walters R. A. “A semi implicit finite element model for non- hydrostatic (dispersive) surface waves”. International Journal for Numerical Methods in Fluids, Vol. 49, No. 7, pp.721–737, 2005.

Wei Z., Jia Y. “A depth-integrated non-hydrostatic finite element model for wave propagation”. International Journal for Numerical Methods in Fluids, Vol. 73, July, 2013.

Rosman P. C. Referência Técnica do SisBahia©.(http://www.sisbahia.coppe.ufrj.br/SisBAHIA_RefTec_V95.p df), 2015.

Brooks A., Hughes T. “Stream-line upwind/Petrov Galerkin formulation for Convection dominated flows with particular emphasis on the incompressible Navier-Stokes equation”. Comp. Meth. Appl. Mech. Eng., 32:199–259, 1982.

Zijlema M., Stelling G. S. “Further experiences with computing non-hydrostatic free-surface flows involving water waves.” International Journal for Numerical Methods in Fluids, Vol. 48, No. 2, pp.169–197, 2005.

Beji S., Battjes J. A. “Experimental investigation of wave propagation over a bar”. Coastal Engineering, Vol. 19, No. 1, pp. 151–162, 1993.

Luth H. R., Klopman G., Kitou N. “Project 13G: kinematics of waves breaking partially on an offshore bar; LDV measurements for waves with and without a net onshore current.” Technical Report H1573, Delft Hydraulics, Delft, The Netherlands, 1994.

Roeber V, Cheung K. F., Kobayashi M. H. “Shock-capturing Boussinesq-type model for nearshore wave processes”. Coastal Engineering, Vol. 57, No. 4, pp. 407–423, 2010.

Briggs, M.J., Synolakis, C.E., Harkins, G.S., Green, D.R., 1995. Laboratory experiments of tsunami runup on a circular island. Tsunamis: 1992–1994, Springer, pp. 569–593, 1995.

Yeh, H.H.J., Liu, P.L.F., Synolakis, C.E. Long-wave Runup Models, World Scientific, Singapore, 1996.

Fuhrman, D.R., Madsen, P.A. “Simulation of nonlinear wave run-up with a high order Boussinesq model”. Coastal Engineering 55, 139–154, 2008.

Zijlema, M., Stelling, G.S., Smit, P. “SWASH: an operational public domain code for simulating wavefields and rapidly variedflows in coastal waters”. Coastal Engineering 58, 992– 1012, 2011.

Watson, G., Barnes, T.C.D. Peregrine, D.H “Numerical modelling of solitary wave propagation and breaking on a beach and runup on a vertical wall”. In: Yeh, H.H.J., Liu, P.L.F., Synolakis, C.E. (Eds.), Long-wave Runup Models. World Scientific, Singapore, pp. 291–297, 1996.

Wei Z., Jia Y., “Simulation of nearshore wave processes by a depth-integrated non-hydrostatic finite element model”. Coastal Engineering 83, 93-107, 2014.