Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 Modeling and Parameter Identification of Shipping Water Events: A Numerical Evaluation with Non-Integer Derivatives ⋆ Marcos A. González-Olvera ∗ Lizeth Torres ∗∗ Jassiel V. Hernández-Fontes ∗∗∗ ∗ Universidad Autónoma de la Ciudad de México, Ciudad de México, México. (e-mail: marcos.angel.gonzalez@uacm.edu.mx). ∗∗ Instituto de Ingeniería, Universidad Nacional Autónoma de la México, Ciudad de México, México. (e-mail: ftorreso@iingen.unam.mx) ∗∗∗ Departamento de Engenharia Naval, Escola Superior de Tecnologia, Universidade do Estado do Amazonas, Manaus 69050-020, Brasil (e-mail: jvfontes@uea.edu.br) Abstract: In this article, we explore the use of a generalized diffusion equation (GDE) for modeling shipping water (SW) events, which are bodies of water that flood the deck of marine structures during extreme weather conditions. The GDE is obtained by substituting the first-order time derivative of the standard diffusion equation by a generalized derivative of non-integer order. The GDE, which in this contribution is formulated using the Caputo fractional derivative, may permit to approximate the behavior between diffusion and wave propagation observed in SW events. We present numerical results to demonstrate the agreement between the response of the GDE and the standard advection-difussion equation, which has been previously used to describe SW events. Keywords: Shipping Water Events, Fractional-Time Difussion Equation, System Identification 1. INTRODUCCIÓN Los embarques de agua (EE-A), también conocidos como eventos de embarque de agua en cubierta, son un fenómenos hidrodinámicos en los que una masa de agua invade la cubierta de una estructura marítima (Goda and Miyamoto, 1976; Greco et al., 2005; Buchner, 2002). Se ha demostrado que estos eventos representan un riesgo importante en la operación de estas estructuras, ya que el peso agregado puede superar las tolerancias de diseño de la estructura (Pakozdi et al., 2016; Silva et al., 2017; Hernández-Fontes et al., 2020a), la tripulación puede resultar gravemente herida y los equipos pueden resultar dañados por las olas de agua que se propagan por la cubierta (Buchner, 2002). En la literatura se han propuesto diversas formas de modelar este tipo de fenómenos y de esta forma predecir su dinámica (Chuang et al., 2018; Zhang et al., 2019). La ecuación de difusión advección (EDA) se ha propuesto recientemente para modelarlos EE-A, la cual, además de proporcionar una aceptable descripción de la dinámica de estos fenómenos, es simple y puede resolverse analí⋆ Los autores desean agradecer a CONACyT por su apoyo a este trabajo. M.A.G.-O. desea agradecer al proyecto UACM PI-CCyT2019-09 por el apoyo dado a este trabajo. Número Especial 2020 ticamente para ciertas condiciones, lo que significa una ventaja ya que las soluciones analíticas de esta ecuación pueden utilizarse en aplicaciones numéricas de bajo costo computacional (Hernández-Fontes et al., 2019a, 2020b, 2019b). Recientemente, Hernández-Fontes et al. (2019a) encontraron una solución analítica para la EDA con coeficientes de advección y difusión constantes, la cual utilizaron para expresar la dinámica de un sólo evento de embarque en una estructura fija. Después, utilizaron esta misma solución para estimar las cargas en cubierta (HernándezFontes et al., 2019b) y para modelar múltiples eventos (Hernández-Fontes et al., 2020b). Posteriormente, utilizaron las solución analítica de la EDA, pero con coeficientes variables, en el modelado de embarques de agua consecutivos generados con olas más largas (Hernández-Fontes et al., 2020c). Los resultados mostraron el gran potencial de la EDA. Sin embargo, aún se requiere más investigación para representar mejor la evolución de eventos generados con olas de diferente inclinación, particularmente olas más empinadas. Por lo anterior, en este artículo se explora el uso de la ecuación de difusión generalizada para modelar embarques de agua en cubierta, que es una extensión de la ecuación clásica de difusión de orden entero, ya que el 1 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 orden de la derivada temporal puede tomar cualquier valor del conjunto de los números reales. En particular, cuando el orden de la derivada temporal varía entre uno y dos, es decir, α ∈ [1, 2], la EDG permite describir un proceso de evolución que es intermedio entre la difusión clásica y la propagación de onda clásica (Huang and Liu, 2005). Este proceso intermedio se ha denominado anteriormente fenómeno de ondas de difusión fraccionarias y sus efectos son ondas de difusión fraccionarias, que son una especie de ondas amortiguadas que involucran oscilaciones (Mainardi, 2010). Asimismo, se evalúa numéricamente el uso de la ecuación de difusión generalizada para aproximar la dinámica de los EE-A, cuyo orden de la derivada temporal puede ser cualquier valor real en lugar de un entero. En particular, se propone que la derivada temporal se ubique en el intervalo (1,2), por lo que describe la dinámica de los EEA con un comportamiento intermedio entre la difusión clásica y la propagación de onda, con menos términos. Por lo tanto, la dinámica estaría ubicada en el llamado fenómeno de onda de difusión fraccionaria que da como resultado ondas de difusión fraccionaria (Mainardi, 2010). Para evaluar la capacidad de EDG de predecir eventos de agua, se llevaron a cabo diferentes simulaciones numéricas con la EDA, considerando a ésta como el modelo patrón o modelo de referencia. De estas simulaciones se obtuvieron series de tiempo de elevaciones de agua en varias posiciones a lo largo de la cubierta de una estructura fija. Parte de los datos obtenidos durante las simulaciones se utilizaron para calibrar la EDG, i.e., para estimar tanto el orden fraccionario como el coeficiente de pseudo-difusión que hacen que la EDG se ajuste mejor a los datos de calibración. La calibración se abordó como un problema de optimización, el cual se llevó a cabo utilizando un algoritmo genético y la transformada pulso. Esta última facilitó el proceso de calibración. El presente artículo se encuentra organizado de la siguiente forma: en la Sección 2 se presenta el modelo que se propone para modelar los embarques de agua, es decir, la EDG. En la Sección 4 se muestran y discuten los resultados numéricos de aplicar la EDG como los eventos EEA, así como el algoritmo de identificación basado en un algoritmo genético. Finalmente, se discuten las conclusiones y el trabajo futuro en la Sección 5. 2. ANTECEDENTES 2.1 Derivada de orden fraccionario La integral o derivada de orden fraccionario puede verse como la generalización de la integral o derivada de orden entero de una determinada función f (t). En este sentido, entre varias definiciones, dos de los enfoques principales en los sistemas de control son el operador fraccionario de Riemann-Liouville y el operador de Caputo (Gorenflo and Mainardi, 1997). Número Especial 2020 El operador diferencial fraccionario de Caputo (1967), que se define como Dα f (t) :=  Z t 1 f (n)(τ )   dτ ; n − 1 < α < n ∈ N, Γ(n − α) a (t − τ )α+1−n n   d f (t); α = n ∈ N, dtn es el más utilizado en control automático, ya que puede incorporar condiciones iniciales en ecuaciones diferenciales de orden fraccionario de manera similar a su contraparte de orden entero. Sin embargo, los valores de inicialización no se pueden dar con tanta fluidez como en el caso del orden de números enteros, y se debe tener especial cuidado en la definición de las mismas (Lorenzo and Hartley, 2001). Por lo anterior, en este trabajo se emplea esta definición de la derivada de orden fraccionario, considerando que al inicio de los fenómenos estudiados la condición inicial es nula y constante para t ≤ 0. 2.2 Ecuación de advección-difusión La ecuación de difusión advección, que previamente se ha utilizado como modelo de eventos de agua en embarques, está dada por ∂η(x, t) ∂ 2 η(x, t) ∂η(x, t) + A(x) = B(x) , (1) ∂t ∂x ∂x2 donde η(x, t) es la elevación del agua, x denota la coordenada del dominio espacial por el que se propaga el flujo de agua (esto es, la cubierta de una estructura marina) y t es la coordenada temporal. A(x) y B(x) son los coeficientes de advección y difusión, respectivamente, los cuales dependen de la coordenada espacial. Estos coeficientes describen la disipación de energía y dependen de diversos factores, como la velocidad horizontal del líquido, la rugosidad de la superficie y la viscosidad del fluido. Más detalles pueden ser encontrados en el trabajo de Hernández-Fontes et al. (2019a). En este artículo la EDA es utilizada como modelo patrón o modelo de referencia, es decir, que se usa para validar el modelo que se propone, es decir, la ecuación de difusión generalizada que a continuación se presenta. 3. ECUACIÓN DE DIFUSIÓN GENERALIZADA En este artículo, consideramos la ecuación de difusión generalizada como un nuevo modelo candidato para eventos de agua embarcada. Esta ecuación se obtiene reemplazando la derivada temporal de primer orden en la ecuación de difusión estándar por una derivada generalizada de orden α ∈ R: ∂ 2 η(x, t) ∂ α η(x, t) = A(x) , (2) α ∂t ∂x2 donde, para EE-A, la variable η(x, t) es la elevación del agua sobre la superficie horizontal, x y t son las coordenadas espacial y temporal, respectivamente. Cuando α = 1, el parámetro A(x) actua como coeficiente de difusión dependiente de la variable espacial, mientras que cuando α = 2 actua como la velocidad de la propagación de onda. En este trabajo este coeficiente se considera constante en 2 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 el tiempo y variante en el espacio. Sin embargo, en trabajos futuros se llevará a cabo un análisis de la EDG con un coeficiente dependiente del tiempo o incluso dependiente de la propia variable η(x, t). Fα =       La versión espacial-discreta de la EDG se obtuvo empleando el método de diferencias finitas, que implicó aproximar ∂ 2 η(x, t)/∂x2 de la siguiente forma: ∂ 2 η(x, t) ηi+1 (t) − 2ηi (t) + ηi−1 (t) ≈ , 2 ∂x (∆x)2 (3) de tal manera que la versión discreta en el espacio de (2) resultó (αi ) ηi (4) = A′i (ηi+1 − 2ηi + ηi−1 ). donde A′i = A(xi )/∆x2i , xi = i∆x and i = 1, 2, ..., N . g2 g1 0 .. . 0 g3 g2 g1 .. . 0 (10) y g1 = 1, gj = j α+1 −2(j−1)α+1 +(j−2)α+1 , j = 2, . . . , N . Esto resulta en ηi (t) − ai J α ηi (t) = ( bi,1 , bi,2 ) J α ui (t). (11) Considerando que
T T ψ(N ) η̄i+1 ψ(N ) , ηi (t) ≈ η̄iT ψ(N ) (t), ui (t) ≈ η̄i−1 entonces (11) se aproxima mediante η̄iT (IN ×N − ai Fαi )ψN (t) = (bi,1 ūi,1 + bi,2 ūi,1 )Fαi ψN (t). (12) 3.1 Transformación pulso En un marco experimental, todas las señales se obtendrían de sensores con un tiempo de muestreo dado y se ubicarían en ubicaciones específicas de la variable espacial. Por lo tanto, como se muestra en Tang et al. (2015); GonzálezOlvera and Tang (2016), en (4), ηi−1 y ηi+1 son tomadas como señales externas con ηi como el i-ésimo pseudoestado: Dα ηi (t) = ai ηi (t) + ( bi,1 , bi,2 ) ui (t), ∆ ∆ g1 0 0 .. . 0 α T 1 × N Γ(α + 2)  . . . gN . . . gN −1  . . . gN −2   , ..  .. . .  . . . g1 N ×N  (5) ∆ ∆ donde ai = −2A′i , bi,1 = A′i , bi,2 = A′i , ui (t) = T ( ηi−1 (t) ηi+1 (t) ) . Considerando un tiempo de muestreo presente en experimentos Ts en el intervalo t ∈ [0, T ], se emplean funciones pulso, que están definidas como el conjunto de N funciones continuas a tramos constantes en los subintervalos j t ∈ [ j−1 N T, N T ]:  1, t ∈ Ωi ψi (t) = (6) 0, otherwise. Finalmente, η̄iT = L(αi , ai )(bi,1 ūi,1 + bi,2 ūi,1 )Fαi , donde L(αi , ai ) = (IN ×N − ai Fαi ), resultando así en una ecuación algebraica en lugar de un problema diferencial. 4. RESULTADOS NUMÉRICOS 4.1 Identificación basada en un Algoritmo Genético Con el objetivo de determinar los parámetros que minimizan el error cuadrático del proceso de identificación, y dado que se manifiestan de forma no lineal en (3.1), se emplea un algoritmo genético (GA) para determinar un valor sub-óptimo de α y Ai que minimice a la función de costo V (η̄i ) = η̄iT η̄i . (13) (7) Para este trabajo, los individuos de cada generación son los parámetros codificados (αi , Ai ) para cada pesudoestado ηi mediante números reales en genes, y se clasifican por el valor que resulte de aplicarlos en la función (13). La siguiente generación se crea a partir de la selección de los mejores miembros de una iteración dada (elitismo), y el resto es generado mediante una recombinación de dos genes aleatorios. Finalmente, los individuos resultantes son modificados de forma aleatoria con un determinado porcentaje de mutación para definir la siguiente generación que seá evaluada. los valores muestreados y retenidos en cada subintervalo ηi (t) se aproximan mediante Se considera el siguiente algoritmo para cada i-ésima discretización ηi (t): Estas funciones son ortogonales bajo el producto punto ∆ RT hw(t), z(t)i = 0 x(t)z(t)dt, i.e., hψj (t), ψk (t)i = 0 ∀j 6= k. Entonces, al agrupar las muestras mediante η̄i = (ηi (0), ηi (T /N ), ηi (2T /N ), . . . , T ηi (jT /N ), . . . , ηi (T (N − 1)/N )) , j = 0, . . . , N − 1, ηi (t) ≈ η̃i (t) = η̄iT ψ(N ) (t), (8) 1. Selección del número de genes para cada iteración (o generación) Ng . 2. Generar valores iniciales aleatorios para αi y Ai dentro de un intervalo de búsqueda definido. 3. Obtener los mejores %E individuos por elitismo. 4. Generar el resto de los individuos mediante recombinación. 5. Mutar el %P M de la generación resultante por recombinación. No mutar a los individuos élite. donde ψ(N ) (t) = ( ψ1 (t) ψ2 (t) . . . ψN (t) ). R t f (τ ) 1−α 1 dτ a ηi (t) de (8) y Aplicando J α f (t) = Γ(α) 0 t−τ asumiendo ηi (t) = 0∀t < 0, se obtiene: J α ηi (t) ≈ η̄iT Fα ψ(N ) (t), (9) donde Número Especial 2020 3 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 0.03 Elevation [m] 0.025 0.02 0.015 0.01 0.005 0 0 1 2 3 4 5 6 7 t[s] Figura 2. Condiciones de frontera empleados en la simulación 0.105 1.55 0.1 1.5 0.095 ∆x αi Ai RMS [mm] 0.04 0.08 0.12 0.16 0.20 0.24 1.6900 1.6001 1.6036 1.6000 1.6000 1.6002 0.2151 0.1478 0.1539 0.1590 0.1464 0.0858 0.8535 0.4542 0.4114 0.3891 0.3989 0.6522 Tabla 1. Parámetros obtenidos y el error cuadrático medio (RMS) para cada discretización de la variable espacial. 6. Repetir y regresar a 3) hasta obtener un valor mínimo determinado en (13) o alcanzar un número determinado de iteraciones. Por lo tanto, los parámetros de diseño del GA son Ng , %E and %P M . Con el fin de obtener los valores discretizados de ηi (t) para la identificación, primero se obtuvo la solución numérica de (1) empleando las condiciones de frontera mostradas en la Fig. 2, las cuales se obtuvieron a partir de los datos publicados en Hernández-Fontes et al. (2020b). Posteriormente, se tomaron seis valores cada ∆x = 0.04 m de un dominio espacial con longitud L = 0.24 m para alimentar el algoritmo genético y estimar (αi , Ai ). Note que x = i∆x, donde i = 0, 1, 2, ..., 7. Note también que i = 0 y i = 7 corresponden a las condiciones de frontera η0 y η7 . 4.2 Resultados La estimación paramétrica se obtuvo mediante la aplicación del algoritmo genético descrito en la Sección 4.1 con Ng = 100, Ne = 5 % y %P M = 1 %, y el espacio de búsqueda se consdieró como Ai = [0, 0.5] y αi = [1.2, 1.7]. Los resultados tras la identificación se muestran en la Fig. 4.2 y la Tabla 1. Cada término ηi (t) se obtuvo empleando una implementación del método PECE de tipo AdamsBashforth-Moulton Diethelm and Freed (1998). Número Especial 2020 1.45 0.09 (x) D(x) Figura 1. Algoritmo de recombinación 0.085 1.4 0.08 1.35 0.075 0.07 1.3 0.04 0.08 0.12 0.16 0.2 0.24 x [m] Figura 3. Resultados de la identificación paramétrica. Se puede observar que el uso del EDG puede aproximarse al de la EDA cuando se usa una discretización de la variable espacial en una restricción 1-D, y los parámetros identificados aparecen consistentemente entre[1.3, 1.5] for α y Ai = [0.07, 0.102] para Ai . Asimismo, el error RMS se encuentra por debajo de 1 mm, logrando una reproducción del fenómeno de embarques de agua mediante la EDG. 5. CONCLUSIONES En este trabajo se ha presentado una demostración numérica de las capacidades de la EDG para modelar eventos de embarque. A través de resultados numéricos y simulaciones, se encontró que si la derivada fraccionaria de la EDG se situaba entre la dinámica clásica de difusión y propagación de ondas, puede aproximarse efectivamente a su dinámica reduciendo un término en la EDA y agregando un parámetro adicional asociado con el tiempo término fraccionario. El trabajo futuro se refiere a la adaptación de este método a los datos experimentales obtenidos de embarques de agua en un entorno confinado, para que los resultados puedan validarse en condiciones del mundo real. AGRADECIMIENTOS Los autores desean agradecer a CONACyT por su apoyo a este trabajo. M.A.G.-O. desea agradecer al proyecto UACM PI-CCyT-2019-09 por el apoyo dado a este trabajo. 4 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 , x=0.04 [m]; =1.4721 0.015 0.01 0.005 0 0 1 2 3 4 5 6 Elev. agua [m] Elev. agua [m] 0.01 0.005 0 0 1 2 4 5 6 Elev. agua [m] 0.005 0 0 1 2 3 4 1 2 5 6 4 5 6 7 5 6 7 5 6 7 Aprox. EDG Sol. EDA 0.005 0 1 2 3 4 Tiempo [s] , x=0.24 [m]; =1.5191 6 Aprox. EDG Sol. EDA 0.01 0.005 0 7 Tiempo [s] 3 Tiempo [s] , x=0.16 [m]; =1.3907 4 0.015 Aprox. EDG Sol. EDA 0.01 0 0.01 0 7 Tiempo [s] , x=0.2 [m]; =1.3474 5 0.015 Elev. agua [m] 3 0.01 0.005 0.015 Aprox. EDG Sol. EDA , x=0.08 [m]; =1.4611 Aprox. EDG Sol. EDA 0.015 0 7 Tiempo [s] , x=0.12 [m]; =1.4394 3 0.015 2 0.02 Aprox. EDG Sol. EDA Elev. agua [m] Elev. agua [m] 1 0.02 0 1 2 3 4 Tiempo [s] Figura 4. Resultados de la identificación y de la TFDE. Se muestra en cada figura el resultado de la identificación para diferentes valores de la variable espacial. En línea azul continua se muestra el resultado de la aproximación empleando la EDG. En línea roja discontinua, el resultado de la solución numérica de la EDA. REFERENCIAS Buchner, B. (2002). Green water on ship-type offshore structures. PhD Thesis. Delft University of Technology, Delft, The Netherlands. Ph.D. thesis, Delft University of Technology Delft, The Netherlands. Caputo, M. (1967). Linear models of dissipation whose q is almost frequency independent—ii. Geophysical Journal International, 13(5), 529–539. Chuang, W.L., Chang, K.A., and Mercier, R. (2018). Kinematics and dynamics of green water on a fixed platform in a large wave basin in focusing wave and random wave conditions. Experiments in Fluids, 59(6), 100. Diethelm, K. and Freed, A.D. (1998). The fracpece subroutine for the numerical solution of differential equations of fractional order. Forschung und wissenschaftliches Rechnen, 1999, 57–71. Goda, K. and Miyamoto, T. (1976). A study of shipping water pressure on deck by two-dimensional ship model tests. Journal of The Society of Naval Architects of Japan, 1976(140), 16–22. González-Olvera, M.A. and Tang, Y. (2016). Identification of incommensurate state-space fractional systems. In 2016 IEEE 13th International Conference on Networking, Sensing, and Control (ICNSC), 1–4. IEEE. Gorenflo, R. and Mainardi, F. (1997). Fractional calculus. In Fractals and fractional calculus in continuum mechanics, 223–276. Springer. Greco, M., Faltinsen, O., and Landrini, M. (2005). Shipping of water on a two-dimensional structure. Journal of Fluid Mechanics, 525, 309. Hernández-Fontes, J.V., Esperança, P.T., Silva, R., Mendoza, E., and Sphaier, S.H. (2020a). Violent waterstructure interaction: Overtopping features and vertical Número Especial 2020 loads on a fixed structure due to broken incident flows. Marine Structures, 74, 102816. Hernández-Fontes, J.V., Hernández, I.D., Mendoza, E., and Silva, R. (2020b). Green water evolution on a fixed structure induced by incoming wave trains. Mechanics Based Design of Structures and Machines, 1–29. Hernández-Fontes, J.V., Vitola, M.A., Paulo de Tarso, T.E., and Sphaier, S.H. (2019a). Analytical convolution model for shipping water evolution on a fixed structure. Applied Ocean Research, 82, 415–429. Hernández-Fontes, J.V., Vitola, M.A., Paulo de Tarso, T.E., and Sphaier, S.H. (2019b). Assessing shipping water vertical loads on a fixed structure by convolution model and wet dam-break tests. Applied Ocean Research, 82, 63–73. Hernández-Fontes, J., Torres, L., Mendoza, E., and Escudero, M. (2020c). Identification of the advectiondiffusion equation for predicting green water propagation. Ocean Engineering, 214, 107658. Huang, F. and Liu, F. (2005). The time fractional diffusion equation and the advection-dispersion equation. The ANZIAM Journal, 46(3), 317–330. Lorenzo, C.F. and Hartley, T.T. (2001). Initialization in fractional order systems. In 2001 European Control Conference (ECC), 1471–1476. IEEE. Mainardi, F. (2010). Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models. World Scientific. Pakozdi, C., Östman, A., Ji, G., Stansberg, C.T., Reum, O., Ovrebo, S., Vestbøstad, T., Sorvaag, C., and Ersland, J. (2016). Estimation of wave loads due to green water events in 10000-year conditions on a tlp deck structure. In International Conference on Offshore Mechanics and Arctic Engineering, volume 49927, V001T01A050. American Society of Mechanical Engi- 5 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx Memorias del Congreso Nacional de Control Automático ISSN: 2594-2492 neers. Silva, D.F., Esperança, P.T., and Coutinho, A.L. (2017). Green water loads on fpsos exposed to beam and quartering seas, part ii: Cfd simulations. Ocean Engineering, 140, 434–452. Tang, Y., Liu, H., Wang, W., Lian, Q., and Guan, X. (2015). Parameter identification of fractional order systems using block pulse functions. Signal Processing, 107, 272–281. Zhang, X., Draper, S., Wolgamot, H., Zhao, W., and Cheng, L. (2019). Eliciting features of 2d greenwater overtopping of a fixed box using modified dam break models. Applied Ocean Research, 84, 74–91. Número Especial 2020 6 Copyright©AMCA. Todos los Derechos Reservados www.amca.mx