Numerical simulation of sloshing in liquid storing tanks by an arbitrary Lagrangian-Eulerian strategy

==Resumen==

La agitación de fluidos contenidos en tanques de almacenamiento con superficie libre es numéricamente simulada mediante una formulación lagrangiana-euleriana arbitraria. El fluido es considerado viscoso y de comportamiento newtoniano, mientras que el flujo se asume laminar e incompresible. Se emplea un código computacional particionado y distribuido, que resuelve tres instancias en cada paso de tiempo: (i) la determinación del estado en el fluido, representado por las ecuaciones de Navier–Stokes; (ii) el desplazamiento de la superficie libre; y (iii) la actualización de la posición de los nodos interiores de la malla de elementos finitos, que es deformada como consecuencia del desplazamiento de la superficie libre. El propósito del trabajo es verificar la aplicabilidad del método en problemas de agitación de solución conocida, así como también resolver algunos ejemplos prácticos. Los ejemplos numéricos incluyen la validación con soluciones analíticas, en las cuales el período de la onda y la tasa de amortiguamiento viscoso son bien capturadas, comparaciones con soluciones de referencia tomadas de otros autores y un caso de agitación inducida por acción sísmica.

==Abstract==

Sloshing of fluids with a free surface contained in liquid storage tanks is numerically simulated by an arbitrary Lagrangian-Eulerian formulation. The fluid is considered viscous and Newtonian, while the flow is assumed laminar and incompressible. A partitioned and distributed computational code is employed, which solves three instances each time step: (i) the determination of the fluid state, given by the Navier–Stokes equations; (ii) the displacement of the free surface; and (iii) the update of the position of the internal nodes of the mesh, that is deformed as a consequence of the free surface displacement. The purpose of the work is verifying the applicability of the method to sloshing problems with known solutions, as well as the resolution of some practical examples. Numerical examples include validations against analytical solutions, where the wave period and damping rate are well captured, comparisons with reference results from other authors and a sample of sloshing induced by seismic actions.

==Palabras clave==

Mecánica de fluidos ; Superficie libre ; Agitación ; Método de los elementos finitos ; Estrategia lagrangiana-euleriana arbitraria ; Movimiento de malla

==Keywords==

Fluid mechanics ; Free surface ; Sloshing ; Finite element method ; Arbitrary Lagrangian Eulerian strategy ; Mesh movement

==1. Introducción==

Es conocida la gran variedad de casos en los cuales se presenta el flujo con Superficie Libre (SL) en problemas de ingeniería, entre los cuales pueden citarse el flujo en canales abiertos [[#bib0070|[14]]]  u otras estructuras civiles [[#bib0215|[43]]] , la hidrodinámica naval [[#bib0045|[9]]]  and [[#bib0210|[42]]]  y la agitación de líquidos contenidos en tanques [[#bib0140|[28]]]  que se emplean para su transporte o simple almacenamiento.      

En particular, la agitación de fluido en tanques reviste interés en distintas circunstancias. Entre ellas, se pueden mencionar casos de recipientes de almacenamiento de fluidos que sufren una excitación producto de, por ejemplo, acciones sísmicas, induciendo la variación de las aceleraciones en distintas direcciones [[#bib0080|[16]]] . Otros casos incluyen la detención súbita de un vehículo de transporte de líquidos, o la agitación que se produce en el interior de las bodegas de los buques cisterna como consecuencia del movimiento de la nave inducido por el oleaje, en los cuales se corre el riesgo de que la frecuencia de agitación sea cercana a alguno de los modos propios más bajos del sistema, por lo que este pueda entrar en resonancia; en esta eventualidad, podría producirse un movimiento del fluido en el interior del tanque suficientemente violento como para afectar la estabilidad del vehículo o nave [[#bib0055|[11]]] .      

Los métodos analíticos o semianalíticos para la resolución de problemas de agitación resultan limitados por la complejidad que pueden presentar tanto la geometría involucrada como las excitaciones antes presentadas, razón por la cual una amplia variedad de métodos numéricos son empleados para la resolución de situaciones reales, entre ellos los métodos de Elementos Finitos (MEF), y los de diferencias finitas y de volúmenes finitos.

Existen distintas propuestas de resolución para problemas de flujos con interfases móviles, entre los cuales el flujo con SL es un caso particular: las técnicas de captura de interfase y las de seguimiento de interfase [[#bib0185|[37]]] . En el primer grupo, los métodos más difundidos hoy día son el de ''Volume of Fluid''  (VOF)  [[#bib0085|[17]]]  and [[#bib0175|[35]]]  y el de pseudoconcentración, o de ''Level Set''  (LS)  [[#bib0135|[27]]]  and [[#bib0180|[36]]] . A grandes rasgos, estos métodos consisten en representar el dominio del problema con una malla fija, capturando la posición de la interfase mediante una fracción de fluido, en el caso de VOF, o mediante con una función marcadora, en LS, lo cual exige modelar tanto la fase líquida como la fase gaseosa. Además, estas formulaciones cuentan con la ventaja de admitir la rotura y/o el pliegue de la SL. Por el otro lado, las técnicas de seguimiento de interfase consisten en definir a la SL mediante entidades tales como partículas, nudos, aristas o, en el caso del MEF, caras de elementos, de manera tal que al evolucionar el problema en el tiempo se realiza un seguimiento preciso de lo que ocurre con la interfase. Entre estos métodos pueden citarse el ''Marker-and-Cell''  (MAC)  [[#bib0075|[15]]]  and [[#bib0130|[26]]] , los métodos de partículas, combinaciones de elementos finitos con partículas [[#bib0100|[20]]] , y la estrategia denominada lagrangiana euleriana arbitraria (ALE, por ''Arbitrary Lagrangian-Eulerian'' )  [[#bib0095|[19]]]  and [[#bib0090|[18]]] , que es la adoptada en este trabajo como metodología de resolución. En esta, el dominio del problema tiene como parte de su frontera a la SL, que es desplazada como consecuencia del movimiento del fluido, de manera que la malla cambia de forma al evolucionar temporalmente el problema. Para ello, se determina el valor de los desplazamientos de las entidades de la malla sobre la cual está definida la SL en función de las velocidades del fluido, y se adopta algún criterio para que la grilla sobre la cual se plantea el método numérico se adapte a la nueva forma del dominio, o sea regenerada. El propósito de este trabajo es resolver numéricamente problemas de agitación transientes en tanques que contienen líquidos viscosos, incompresibles y de comportamiento newtoniano, cuyo flujo es regido por las ecuaciones de Navier–Stokes (NS), aplicando una técnica ALE en la cual la deformación instantánea del dominio es resuelta con una técnica de movimiento de malla, en dos y tres dimensiones espaciales. Para tal fin, se resuelven mediante el MEF tres instancias distintas en cada paso de tiempo: el estado en el fluido, el desplazamiento de las entidades que definen la posición de la SL y el movimiento de la malla. Los ejemplos incluyen validaciones con soluciones analíticas, tanto viscosas como no viscosas, y con resultados numéricos de referencia.      

==2. Ecuaciones de gobierno==

La presente metodología consiste en resolver en cada paso de tiempo las siguientes instancias, una a continuación de la otra: (i) el estado en el dominio del flujo, obteniendo el campo de velocidades y el de presión; (ii) el desplazamiento de las entidades que definen la posición instantánea de la SL, con las velocidades calculadas en la instancia anterior; y (iii) la actualización de la posición de los nodos interiores del dominio mediante una estrategia de movimiento de malla. Cada una de estas instancias es resuelta por separado, conformando así un esquema de acoplamiento débil entre ellas e implementado en el programa PETSc-FEM [[#bib0150|[30]]] , que usa herramientas de ''Portable Extensible Toolkit for Scientific Computations''  (PETSc)  [[#bib0010|[2]]]  y de ''Message Passing Interface''  (MPI)  [[#bib0120|[24]]]  para efectuar el cálculo en paralelo de cada una de las etapas antes indicadas.      

===2.1. Descripción del movimiento mediante un paradigma ALE===

Hay 3 formas de representar el movimiento del fluido: lagrangiana, euleriana y arbitraria o ALE. La configuración material o lagrangiana ''Ω''<sub>'''''X'''''</sub>  está constituida por el conjunto o sistema de partículas materiales asociadas a las coordenadas '''''X''''' . Por otro lado, la descripción espacial o euleriana ''Ω''<sub>'''''x'''''</sub>  se compone de puntos correspondientes a las coordenadas espaciales '''''x''''' . Por último, la configuración de referencia ''Ω''<sub>'''''χ'''''</sub> , que es arbitraria, es la asociada a los datos de la discretización espacial del método numérico, y cuyas coordenadas referenciales '''''χ'''''  son las empleadas en la resolución numérica.      

La velocidad absoluta de una partícula con respecto al marco euleriano, o velocidad material, es definida como:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}\left(\boldsymbol{X},t\right)={\frac{\partial \boldsymbol{x}}{\partial t}}_\boldsymbol{X};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

y coincide con la velocidad de una partícula material '''''X'''''  que en el instante ''t''  se encuentra en la posición espacial '''''x'''''  en la cual se evalúa la ecuación [[#eq0005|(1)]] . Esta derivada parcial, calculada «a '''''X'''''  fija», es la que aparece como efecto convectivo al plantear las ecuaciones sobre un marco euleriano e indica la velocidad de desplazamiento entre el fluido y ese referencial euleriano. Un procedimiento mediante el cual puede arribarse a las expresiones consideradas puede seguirse en Donea y Huerta [[#bib0050|[10]]]  o en Folch Durán [[#bib0060|[12]]] .      

Por otro lado, al emplear para el movimiento una formulación ALE, en la cual el dominio de referencia no coincide con el dominio espacial, la velocidad relativa de la malla con respecto al marco euleriano es dada por la expresión [[#bib0050|[10]]]  and [[#bib0060|[12]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\mbox{ˆ}}{\boldsymbol{v}}\left(\boldsymbol{\chi },t\right)=</math><math>{\frac{\partial \boldsymbol{x}}{\partial t}}_{\boldsymbol{\chi }}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

La velocidad relativa de la partícula con respecto al marco de referencia es definida como [[#bib0050|[10]]]  and [[#bib0060|[12]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{w}\left(\boldsymbol{X},t\right)={\frac{\partial \boldsymbol{\chi }}{\partial t}}_\boldsymbol{X}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

Por último, la velocidad convectiva '''''c'''''  se define de la siguiente manera [[#bib0050|[10]]]  and [[#bib0060|[12]]] :

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{c}:=\boldsymbol{v}-\overset{\mbox{ˆ}}{\boldsymbol{v}}=</math><math>\frac{\partial \boldsymbol{x}}{\partial \boldsymbol{\chi }}\boldsymbol{w};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

y corresponde a la velocidad relativa entre las configuraciones material y de referencia. Recordando que <math display="inline">\boldsymbol{w}</math>  es la velocidad relativa de la partícula con respecto a la configuración de referencia, la ecuación [[#eq0020|(4)]]  implica que la velocidad convectiva coincide con <math display="inline">\boldsymbol{w}</math> , es decir <math display="inline">\boldsymbol{c}=</math><math>\boldsymbol{w}</math> , cuando el movimiento del fluido es exclusivamente de traslación, esto es, si ∂'''''x''''' /∂'''''χ'''''  = '''''I''''' .      

El uso de un paradigma ALE requiere la expresión de las ecuaciones de gobierno en la configuración de referencia, para lo cual se necesita establecer las relaciones entre esta y las cantidades intervinientes. La forma a emplear con este fin es la denominada «Ecuación Fundamental ALE»[[#bib0050|[10]]] , que para una cantidad física escalar ''f''  es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\frac{\partial f}{\partial t}}_\boldsymbol{X}={\frac{\partial f}{\partial t}}_{\boldsymbol{\chi }}+</math><math>\boldsymbol{c}\cdot \nabla f.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

===2.2. Ecuaciones en el dominio del fluido===

Las ecuaciones en derivadas parciales que representan el flujo de un fluido viscoso, homogéneo e incompresible son las de NS incompresibles, esto es, la expresión diferencial de conservación de la cantidad de movimiento, junto con la condición de incompresibilidad del fluido. En un referencial euleriano, estas pueden escribirse como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho \left({\partial }_t\boldsymbol{v}+\boldsymbol{v}\cdot \nabla \boldsymbol{v}-\right. </math><math>\left. \boldsymbol{f}\right)-\nabla \cdot \boldsymbol{\sigma }=</math><math>0;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

sobre el dominio del flujo ''Ω''  para el tiempo ''t''  ∈ [0, ''T  '' ], en la cual <math display="inline">\boldsymbol{v}</math>  es la velocidad del fluido, '''''f'''''  es la fuerza de cuerpo, o fuerza por unidad de masa, ''ρ''  la densidad del fluido, ''T''  cierto tiempo final, ∂<sub>''t''</sub>  indica derivación parcial con respecto al tiempo y ∇ = ∇ <sub>'''''x'''''</sub>  es el operador de derivación con respecto a las coordenadas espaciales '''''x'''''  del marco euleriano. El tensor <math display="inline">\boldsymbol{\sigma }=</math><math>\boldsymbol{\sigma }(\boldsymbol{v},p)</math>  es el de tensiones en el fluido, que puede descomponerse aditivamente en una parte isotrópica −''p'''''''I'''''  y otra desviadora '''''T''''' , en la forma:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\sigma }=-p\boldsymbol{I}+\boldsymbol{T};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

donde ''p''  es la presión e '''''I'''''  es el tensor identidad. Dado que únicamente se considerarán fluidos de comportamiento newtoniano, la parte desviadora '''''T'''''  se relaciona linealmente con la tasa de deformación <math display="inline">\boldsymbol{\epsilon }=</math><math>\boldsymbol{\epsilon }(\boldsymbol{v})</math>  como sigue,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{T}=2\mu \boldsymbol{\epsilon };</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

en la cual

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\epsilon }=\frac{1}{2}\left[\nabla \boldsymbol{v}+\right. </math><math>\left. {\left(\nabla \boldsymbol{v}\right)}^T\right];</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

con ''μ''  la viscosidad dinámica del fluido y donde (. . .)<sup>''T''</sup>  indica transposición.      

Las expresiones dadas son referidas a coordenadas espaciales '''''x''''' , y no contemplan la posibilidad de la deformación del dominio ''Ω'' . Debido a que el problema del flujo debe plantearse en un dominio cuya forma geométrica varía como consecuencia del desplazamiento y la deformación de la SL, se establece una relación entre las ecuaciones antes consideradas, para el dominio espacial fijo ''Ω''<sub>'''''x'''''</sub> , y las que relacionan estas con un dominio de referencia arbitrario ''Ω''<sub>'''''χ'''''</sub> . Para ello, se utilizan las formas ALE de las ecuaciones de conservación de cantidad de movimiento y de continuidad [[#bib0050|[10]]] , dadas por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho \left({\partial }_t\boldsymbol{v}+\boldsymbol{c}\cdot \nabla \boldsymbol{v}-\right. </math><math>\left. \boldsymbol{f}\right)-\nabla \cdot \boldsymbol{\sigma }=</math><math>0;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

Nótese que solo la ecuación de cantidad de movimiento se ve afectada por la velocidad del dominio de referencia <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{v}}</math> , que interviene a través de la velocidad convectiva '''''c'''''  dada por la ecuación [[#eq0020|(4)]] .      

Las condiciones de contorno en las fronteras ''Γ''  del dominio son las siguientes:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}=\boldsymbol{v}_D\quad en\quad {\Gamma }_D;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

donde las condiciones Dirichlet están dadas sobre ''Γ''<sub>D</sub> , que corresponde a los contornos sólidos cuya velocidad <math display="inline">\boldsymbol{v}_D</math>  es conocida, y ''Γ''<sub>'''''t'''''</sub>  es una superficie de tracción o interfase entre 2 fluidos, donde actúan las fuerzas de tracción '''''t''''' . Además, se verifica que ''Γ''  = ''Γ''<sub>D</sub>  ∪ ''Γ''<sub>'''''t'''''</sub>  y que ''Γ''<sub>D</sub> ∩ ''Γ''<sub>'''''t'''''</sub>  = ∅.      

En el caso de SL, y para una aproximación lagrangiana con representación homogénea del fluido, la condición de contorno sobre la SL, ''Γ''<sub>''SL''</sub>  ⊆ ''Γ''<sub>'''''t'''''</sub> , consiste en el equilibrio entre las proyecciones en dirección normal a la interfase de los tensores de tensiones en el líquido y en el gas, '''''σ'''''<sub>''l''</sub>  y '''''σ'''''<sub>''g''</sub> , respectivamente,

<span id='eq0090'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\sigma }}_l\cdot \boldsymbol{n}={\boldsymbol{\sigma }}_g\cdot \boldsymbol{n}\quad \quad en\quad {\Gamma }_{SL}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

Debido a que la viscosidad y la densidad de la fase gaseosa son consideradas despreciables, tenemos que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{T}\cdot \boldsymbol{n}=\boldsymbol{0}\quad en\quad {\Gamma }_{SL};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

en las cuales ''P''<sub>atm</sub>  es la presión que ejerce el gas sobre el líquido, normalmente la atmosférica, y '''''T'''''  representa las tensiones producidas por el gas sobre la interfase. Entonces, las fuerzas de tracción sobre ''Γ''<sub>''SL''</sub>  valen '''''t'''''  = − ''P''<sub>atm</sub>  '''''n''''' , o bien:

<span id='eq0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\sigma }}_g\cdot \boldsymbol{n}=-P_{atm}\quad \boldsymbol{n}\quad sobre\quad {\Gamma }_{SL}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

En el caso de adoptarse ''P''<sub>atm</sub>  = 0, el término de contorno sobre ''Γ''<sub>''SL''</sub>  desaparece. Notar que los efectos de la tensión superficial no son tenidos en cuenta en el modelo, y que la presión del lado del gas sobre la SL no es igual a la que se produce del lado del fluido, dada la ecuación [[#eq0090|(15)]] , que impone la igualdad de tracciones.      

Las condiciones iniciales en los problemas a ser resueltos por las ecuaciones de NS, en sus variantes para fluidos homogéneos dadas por las ecuaciones   and  , o bien las   and  , son indicadas sobre la velocidad:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}(\boldsymbol{x},0)=\boldsymbol{v}_0(\boldsymbol{x});</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

esto es, <math display="inline">\boldsymbol{v}_0</math>  es el campo de velocidades inicial definido sobre todo el dominio en tiempo ''t''  = 0.      

La discretización se realiza sobre las ecuaciones   and [[#bib0190|[38]]] . Para ello, el dominio ''Ω''  se divide en ''n''<sub>el</sub>  elementos finitos ''Ω''<sup>''e''</sup> , siendo <math display="inline">E</math>  el conjunto de todos los elementos y ''H''<sup>1''h''</sup>  el espacio de dimensión finita definido como sigue,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>H^{1h}=\lbrace {\psi }^h\vert {\psi }^h\in C^0(\overline{\Omega }),{\psi }^h\vert _{{\Omega }^e}\in P^1,\forall {\Omega }^e\in E\rbrace ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

donde ''P''<sup>1</sup>  es el espacio de polinomios de primer orden. Los espacios de funciones de peso e interpolación son

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>S_v^h=\lbrace \boldsymbol{v}^h\vert \boldsymbol{v}^h\in (H^{1h})^{n_d},\mbox{con}\mbox{ }\boldsymbol{v}^h\doteq \boldsymbol{v}_D\quad sobre\quad {\Gamma }_D\rbrace ;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

siendo ''n''<sub>''d''</sub>  el número de dimensiones espaciales. El método estabilizado con ''streamline upwind/Petrov-Galerkin''  (SUPG)  [[#bib0035|[7]]]  y ''pressure stabilizing/Petrov-Galerkin''  (PSPG)  [[#bib0220|[44]]]  and [[#bib0225|[45]]]  se escribe de la siguiente forma: ''hallar  ''<math display="inline">\boldsymbol{v}^h\in S_v^h</math>''y  ''<math display="inline">p^h\in S_p^h</math>''tal que:''

<span id='eq0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{\Omega }\boldsymbol{w}^h\cdot \rho \left({\partial }_t\boldsymbol{v}^h+\boldsymbol{v}^h\cdot \nabla \boldsymbol{v}^h-\boldsymbol{f}\right)\quad d\Omega </math>
|-
|<math>+{\int }_{\Omega }\boldsymbol{\epsilon }(\boldsymbol{w}^h):{\boldsymbol{\sigma }}^hd\Omega </math>
|-
|<math>+\sum_{e=1}^{n_{el}}{\int }_{{\Omega }^e}{\boldsymbol{\delta }}^h\cdot \left[\rho \left({\partial }_t\boldsymbol{v}^h+\boldsymbol{v}^h\cdot \nabla \boldsymbol{v}^h-\boldsymbol{f}\right)-\nabla \cdot {\boldsymbol{\sigma }}^h\right]d{\Omega }^e</math>
|-
|<math>+\sum_{e=1}^{n_{el}}{\int }_{{\Omega }^e}{\boldsymbol{\iota }}^h\cdot \left[\rho \left({\partial }_t\boldsymbol{v}^h+\boldsymbol{v}^h\cdot \nabla \boldsymbol{v}^h-\boldsymbol{f}\right)-\nabla \cdot {\boldsymbol{\sigma }}^h\right]d{\Omega }^e</math>
|-
|<math>+{\int }_{\Omega }q^h\quad \nabla \cdot \boldsymbol{v}\quad d\Omega </math>
|-
|<math>+\sum_{e=1}^{n_{el}}{\int }_{{\Omega }^e}{\boldsymbol{\kappa }}^h\quad \nabla \cdot \boldsymbol{v}^hd{\Omega }^e</math>
|-
|<math>={\int }_{{\Gamma }_\boldsymbol{t}}\boldsymbol{w}^h\boldsymbol{t}^h\quad d\Gamma \quad \forall \boldsymbol{w}^h\in V_v^h,\forall q^h\in S_p^h;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 24)
|}

en la cual el término que contiene '''''δ'''''<sup>''h''</sup>  es el de SUPG, el de PSPG es el afectado por '''''ι'''''<sup>''h''</sup> , y el de estabilización de mínimos cuadrados sobre la restricción de incompresibilidad (''Least-Squares on Incompressibility Constraint'' , LSIC)  [[#bib0230|[46]]]  es el de '''''κ'''''<sup>''h''</sup> , con <math display="inline">{\boldsymbol{\sigma }}^h=</math><math>\boldsymbol{\sigma }(\boldsymbol{v}^h,p^h)</math> . Los parámetros de estabilización son dados por las siguientes expresiones:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\delta }}^h={\tau }_{SUPG}\quad (\boldsymbol{v}^h\cdot \nabla )\boldsymbol{w}^h;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 25)
|}

con el tiempo característico de SUPG dado como [[#bib0230|[46]]] :

<span id='eq0170'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau _{SUPG}={\left(\frac{1}{{\tau }_{S1}^2}+\frac{1}{{\tau }_{S2}^2}+\right. </math><math>\left. \frac{1}{\tau _{S3}^2}\right)}^{-\frac{1}{2}};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 28)
|}

siendo

<span id='eq0175'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau _{S1}=\frac{h}{2\parallel \boldsymbol{v}\parallel };\quad \tau _{S2}=</math><math>\frac{\Delta t}{2};\quad \tau _{S3}=\frac{h^2}{4\nu }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 29)
|}

El denominado tiempo característico para PSPG es [[#bib0230|[46]]] :

<span id='eq0180'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tau }_{PSPG}={\left(\frac{1}{{\tau }_{P1}^2}+\frac{1}{{\tau }_{P2}^2}+\right. </math><math>\left. \frac{1}{{\tau }_{P3}^2}\right)}^{-\frac{1}{2}};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 30)
|}

con

<span id='eq0185'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tau }_{P1}=\frac{h^\#}{2\parallel \boldsymbol{v}\parallel };\quad {\tau }_{P2}=</math><math>\frac{\Delta t}{2};\quad {\tau }_{P3}=\frac{(h^\#)^2}{4\nu }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 31)
|}

El parámetro restante, la viscosidad cinemática característica, es:

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\nu }_{LSIC}=\frac{\parallel \boldsymbol{v}\parallel h}{2}\quad z(Pe);</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 32)
|}

en la cual el número de Péclet se calcula como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Pe=\frac{\parallel \boldsymbol{v}\parallel h}{2\nu }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 33)
|}

Nótese que si <math display="inline">\parallel \boldsymbol{v}\parallel =</math><math>0</math> , algunas de las expresiones en las ecuaciones [[#eq0175|(29)]]  o [[#eq0185|(31)]]  son singulares; sin embargo, las restantes no lo son, y los parámetros de estabilización de las ecuaciones [[#eq0170|(28)]]  y [[#eq0180|(30)]]  resultan siempre regulares.      

En un paradigma ALE, la velocidad con la cual se determinan los coeficientes de estabilización de las ecuaciones [[#eq0170|(28)]] , [[#eq0180|(30)]]  y [[#eq0190|(32)]]  es la velocidad convectiva '''''c''''' . Además, el tamaño de elemento ''h''  se determina mediante:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h=2{\left(\sum_{a=1}^{n_{en}}\vert \boldsymbol{s}\cdot \nabla w_a\vert \right)}^{-1};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 34)
|}

donde las <math display="inline">w_a</math>  son las funciones asociadas al nodo ''a'' , ''n''<sub>en</sub>  es el número de nodos en el elemento considerado y '''''s'''''  es un versor orientado según las líneas de corriente. Así, la longitud ''h''<sup>#</sup>  para las expresiones en la ecuación [[#eq0185|(31)]]  viene dada por el diámetro del círculo de igual área que el elemento en 2D, o de la esfera de igual volumen en 3D. Por último, la función ''z'' (Pe) se calcula como sigue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>z(Pe)=\begin{array}{ll}
Pe/3, & 0\leq Pe<3;\\
1, & 3\leq Pe.
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 35)
|}

Una vez realizada la discretización espacial de la ecuación [[#eq0145|(24)]]  se obtiene el sistema de ecuaciones algebraicas que se resuelve numéricamente. La integración temporal consiste en un esquema en diferencias, mediante un método tipo ''θ''[[#bib0050|[10]]] , tal que 0 ≤ ''θ''  ≤ 1. En este esquema, ''θ''  > 1/2 es incondicionalmente estable, correspondiendo ''θ''  = 1 a Retro-Euler y ''θ''  = 0 a Euler hacia adelante.      

Las condiciones de contorno sobre ''Γ''<sub>D</sub> , esto es, sobre la velocidad <math display="inline">\boldsymbol{v}_D</math> , se adoptan de acuerdo al tipo de problema a resolver numéricamente. Las condiciones más usuales son las siguientes.
* ''Pared inmóvil:''  toma la forma

<span id='eq0210'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}_D=\boldsymbol{0},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 36)
|}

dada la hipótesis de fluido viscoso adoptada para el modelo.                          
* ''Puntos de contacto fluido-fluido-sólido:''  en las zonas de frontera sólida cercanas a la SL, la imposición de la ecuación  [[#eq0210|(36)]]  en el problema discreto produce inconvenientes en el movimiento de la malla, ya que las entidades que definen la línea de contacto se verían sometidos a grandes gradientes en la elevación de la SL. Es por este motivo que la condición de no deslizamiento es relajada en la línea de contacto, siendo remplazada, por ejemplo, por la denominada «condición de deslizamiento de Navier»[[#bib0070|[14]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\cdot (\boldsymbol{n}\cdot \boldsymbol{\sigma })=</math><math>-\frac{1}{\beta }(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\cdot (\boldsymbol{v}_S-</math><math>\boldsymbol{v}_D);</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 37)
|}

en la cual <math display="inline">\boldsymbol{v}_S</math>  es la velocidad del fluido en la línea de contacto, mientras que '''''I'''''  − '''''nn'''''  proyecta la componente de velocidades sobre el plano tangente y ''β''  es un coeficiente de deslizamiento de valor arbitrario. Los valores límite de ''β''  son ''β''  = 0 para una condición de no deslizamiento, en tanto que ''β'' → ∞ lleva a considerar un deslizamiento perfecto.                          
* ''Fluidos de baja viscosidad:''  en problemas en los que intervienen líquidos de muy baja viscosidad cinemática, por ejemplo 1 × 10<sup>−6</sup>  m<sup>2</sup> /s. En estos casos, puede emplearse la condición de deslizamiento perfecto <math display="inline">\boldsymbol{v}_D\cdot \boldsymbol{n}=</math><math>0</math>[[#bib0090|[18]]] , [[#bib0195|[39]]]  and [[#bib0160|[32]]] , que requiere un esfuerzo computacional menor que las condiciones antes mencionadas.                          
* ''Planos de simetría u otras fronteras artificiales:  ''  cuando el fluido no circula a través de los contornos, siendo estos ficticios; para estos casos, la condición de borde a aplicar es la de deslizamiento perfecto, esto es, <math display="inline">\boldsymbol{v}_D\cdot \boldsymbol{n}=</math><math>0</math> .

===2.3. Desplazamiento de la superficie libre===

Los desplazamientos de los nodos de la SL se proponen restringidos a una dirección dada por un versor <math display="inline">{\overset{\mbox{ˆ}}{\boldsymbol{s}}}_j</math> , fija, de manera que la nueva posición del nodo ''j''  sobre la SL en el tiempo ''t''  es  [[#bib0020|[4]]] :

<span id='eq0220'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{x}_j(t)=\boldsymbol{x}_{0,j}+{\eta }_j(t)\quad {\overset{\mbox{ˆ}}{\boldsymbol{s}}}_j;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 38)
|}

en la cual <math display="inline">{\overset{\mbox{ˆ}}{\boldsymbol{s}}}_j</math>  da la dirección fija o «espina»''(spine)'' , '''''x'''''<sub>0,''j''</sub>  es la posición inicial del nodo ''j''  y ''η''<sub>''j''</sub>  la coordenada escalar, todos ellos indicados en la [[#fig0005|figura 1]] . Usualmente, <math display="inline">{\overset{\mbox{ˆ}}{\boldsymbol{s}}}_j</math>  es normal a la SL en reposo para el problema a resolver.

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr1.jpg|center|150px|Desplazamiento de los nodos de la SL.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Desplazamiento de los nodos de la SL.

</span>
|}

Es importante remarcar que las espinas se usan exclusivamente para los deplazamientos de la SL, y que en el interior del dominio las posiciones nodales son determinadas mediante alguno de los métodos presentados en la sección siguiente. De esta forma, la distorsión de la malla es reducida en comparación con casos en los cuales la totalidad de los nodos se mueven sobre direcciones fijas.

El movimiento de la SL está regido por la denominada condición cinemática. Esta condición contiene el significado físico de una interfase material, a través de la cual no hay intercambio de fluido [[#bib0065|[13]]] , [[#bib0125|[25]]] , [[#bib0145|[29]]]  and [[#bib0205|[41]]] , y puede expresarse como:

<span id='eq0225'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}\cdot \boldsymbol{n}={\partial }_t\boldsymbol{\eta }\cdot \boldsymbol{n}\quad sobre\quad {\Gamma }_{SL};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 39)
|}

donde '''''n'''''  es la dirección normal a la SL y '''''η'''''  es el desplazamiento o elevación de la interfase sobre la dirección de la espina <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{s}}</math> , es decir <math display="inline">\boldsymbol{\eta }=</math><math>\eta \quad \overset{\mbox{ˆ}}{\boldsymbol{s}}</math> , con ''η''  representando el desplazamiento escalar, esquematizado en la  [[#fig0005|figura 1]] . Remplazando en la ecuación [[#eq0225|(39)]] ,

<span id='eq0230'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\partial }_t\eta =\frac{\boldsymbol{v}\cdot \boldsymbol{n}}{\overset{\mbox{ˆ}}{\boldsymbol{s}}\cdot \boldsymbol{n}}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 40)
|}

Nótese que los desplazamientos considerados son aquellos que se producen en dirección normal a la SL, mientras que los tangenciales son irrelevantes.      

Teniendo en cuenta la ecuación [[#eq0220|(38)]]  y que el campo escalar en 3 dimensiones es ''η''  = ''η'' (''x''<sub>1</sub> , ''x''<sub>2</sub> , ''t'' ), la SL puede expresarse de manera implícita como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F(\boldsymbol{x},t)=x_3-\eta =0;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 41)
|}

asumiendo que la dirección de las espinas se adopta vertical, i.e. <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{s}}\parallel x_3</math> , la normal a la SL puede calcularse como

<span id='eq0240'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{n}=\nabla F(\boldsymbol{x})={\left[-\frac{\partial \eta }{\partial x_1};\quad -\right. </math><math>\left. \frac{\partial \eta }{\partial x_2};\quad +\right. </math><math>\left. 1\right]}^T.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 42)
|}

Remplazando en la ecuación [[#eq0230|(40)]]  con [[#eq0240|(42)]] , e introduciendo el parámetro de proyección

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>H=\overset{\mbox{ˆ}}{\boldsymbol{s}}\cdot \boldsymbol{n};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 43)
|}

esto es, siendo ''H''  la proyección de la normal en la dirección de la espina, resulta:

<span id='eq0250'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\partial }_t\eta =\frac{1}{H}\boldsymbol{v}\cdot \boldsymbol{n}=</math><math>\frac{1}{H}\left[-v_1\frac{\partial \eta }{\partial x_1}-\right. </math><math>\left. v_2\frac{\partial \eta }{\partial x_2}+v_3\right];</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 44)
|}

con lo cual se arriba a una expresión para la elevación ''η''  en la forma de una ecuación de advección. La reducción al caso de SL en dominio bidimensional es inmediata. A partir de la ecuación  [[#eq0250|(44)]] , y con el fin de mantener una notación uniforme en las expresiones subsiguientes, el problema a resolver para obtener el desplazamiento de la SL ''η''  se escribe como el siguiente sistema advectivo,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\partial }_t\eta +\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\nabla }\eta =</math><math>s\quad en\quad {\Omega }_{SL}^{{'}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 45)
|}

en la que la velocidad normal a las espinas es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{v}_{\parallel }=\frac{1}{H}\quad {\left[v_1;\quad v_2\right]}^T;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 47)
|}

mientras que el gradiente bidimensional de la elevación ''η''  es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\mbox{ˆ}}{\nabla }\eta ={\left[\frac{\partial \eta }{\partial x_1};\quad \frac{\partial \eta }{\partial x_2}\right]}^T;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 48)
|}

y <math display="inline">s=v_3/H</math>  es el término fuente dado por la velocidad del fluido en la dirección vertical para cada punto de la SL. Además, <math display="inline">{\Omega }_{SL}^{{'}}=</math><math>{\Gamma }_{SL}</math>  en la ecuación [[#eq0110|(18)]]  es el dominio de interfase de dimensión ''n''<sub>dim</sub>  − 1, correspondiente a un dominio ''n''<sub>dim</sub>  para el problema del fluido con SL, véase la [[#fig0010|figura 2]] . Por último, siendo <math display="inline">{\Gamma }_{SL}^{{'}}=</math><math>\partial {\Omega }_{SL}^{{'}}</math> , sobre <math display="inline">{\Gamma }_D^{{'}}</math>  se impone la condición de contorno Dirichlet para la ecuación de advección dada. Debe tenerse presente que la ecuación [[#eq0260|(45)]]  representa un problema hiperbólico, y por lo tanto <math display="inline">{\Gamma }_D^{{'}}</math>  es el contorno de ingreso en el que se verifica <math display="inline">\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\boldsymbol{n}}<0</math> , con <math display="inline">\overset{\mbox{ˆ}}{\boldsymbol{n}}=</math><math>[n_1;n_2]^T</math> , es decir, <math display="inline">{\Gamma }_D^{{'}}=</math><math>\left\{{\Gamma }_{SL}^{{'}}\vert \boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\boldsymbol{n}}<0\right\}</math> .

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr2.jpg|center|150px|Dominios y contornos para el problema de advección de la SL.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

Dominios y contornos para el problema de advección de la SL.

</span>
|}

La formulación explícita del desplazamiento de la SL dada por la ecuación [[#eq0260|(45)]] , al igual que la dada en la ecuación [[#eq0230|(40)]] , resulta numéricamente inestable para ondas de superficie, dado su carácter advectivo [[#bib0015|[3]]] . En este sentido, la estabilización es necesaria en aquellos problemas en los cuales el término de transporte es dominante frente al término fuente ''s'' . Este hecho ya ha sido indicado por diversos autores, quienes apelaron a diferentes métodos de estabilización para obtener soluciones satisfactorias. Entre estos métodos se encuentran el de ''streamline upwind/Petrov-Galerkin''  (SUPG)  [[#bib0035|[7]]]  aplicado a la SL por Soulaïmani et al. [[#bib0195|[39]]]  y Güler et al. [[#bib0070|[14]]] , y el denominado ''Galerkin/Least-Squares''  (GLS), elegido por Behr y Abraham  [[#bib0025|[5]]] .      

El método de estabilización SUPG, introducido por Brooks y Hughes [[#bib0035|[7]]] , es empleado en este trabajo para la resolución del problema de advección de la ecuación [[#eq0260|(45)]]  mediante Galerkin, como paso previo a la actualización de la malla, realizado sobre los resultados del NS.      

La formulación variacional adoptada como guía es tomada de Donea y Huerta [[#bib0050|[10]]] , y se basa en los espacios funcionales <math display="inline">S_t:=</math><math>\lbrace \eta \vert \eta (\cdot ,t)\in H^1({\Omega }_{SL}^{{'}}),t\in [0,T]\vert \eta =</math><math>{\eta }_D\quad sobre\quad {\Gamma }_D^{{'}}\rbrace </math>  y <math display="inline">W:=</math><math>\lbrace w\in H^1({\Omega }_{SL}^{{'}})\vert w=0\quad sobre\quad {\Gamma }_D^{{'}}\rbrace </math> . El problema se propone de la siguiente forma: ''para cualquier tiempo t''  ∈ [0, ''T'' ]'', hallar  ''<math display="inline">\eta \in S</math>''tal que''

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_{{\Omega }_{SL}^{{'}}}w\quad \frac{\partial \eta }{\partial t}d{\Omega }_{SL}^{{'}}+</math><math>{\int }_{{\Omega }_{SL}^{{'}}}w\left(\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\nabla }\eta \right)d{\Omega }_{SL}^{{'}}=</math><math>{\int }_{{\Omega }_{SL}^{{'}}}w\quad s\quad d{\Omega }_{SL}^{{'}};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 49)
|}

o, expresado en forma compacta,

<span id='eq0285'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(w,\eta \right)+c\left(\boldsymbol{v}_{\parallel };w,\eta \right)=</math><math>\left(w,s\right)\quad en\quad {\Omega }_{SL}^{{'}};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 50)
|}

con condición inicial

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left(w,\eta (\cdot ,0)\right)=\left(w,{\eta }_0\right).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 51)
|}

En la ecuación [[#eq0285|(50)]] , los términos se identifican de la siguiente manera,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(w,\eta \right)={\int }_{{\Omega }_{SL}^{{'}}}w\quad \frac{\partial \eta }{\partial t}\quad d{\Omega }_{SL}^{{'}};</math>
|-
|<math>c\left(\boldsymbol{v}_{\parallel };w,\eta \right)={\int }_{{\Omega }_{SL}^{{'}}}w\left(\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\nabla }\eta \right)d{\Omega }_{SL}^{{'}};</math>
|-
|<math>\left(w,s\right)={\int }_{{\Omega }_{SL}^{{'}}}w\quad s\quad d{\Omega }_{SL}^{{'}}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 52)
|}

Una estabilización consistente para la formulación de la ecuación [[#eq0285|(50)]]  puede escribirse como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(w,\eta \right)+c\left(\boldsymbol{v}_{\parallel };w,\eta \right)+</math><math>\sum_e{\int }_{{\Omega }_e^{{'}}}P(w)\quad {\tau }_S\quad R(\eta )\quad d{\Omega }_e^{{'}}=</math><math>\left(w,s\right)\quad en\quad {\Omega }_{SL}^{{'}};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 53)
|}

con el tiempo intrínseco o parámetro de estabilización ''τ''<sub>''S''</sub> , que se describe más adelante en esta sección, e indicando con <math display="inline">{\Omega }_e^{{'}}</math>  que la integración del término se realiza sobre cada uno de los elementos. El residuo de la ecuación diferencial es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R(\eta )=\frac{\partial \eta }{\partial t}+\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\nabla }\eta -</math><math>s=L(\eta )-s;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 54)
|}

y el operador <math display="inline">P(w)</math>  es definido para una estabilización con SUPG en la forma:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>P(w)=\boldsymbol{v}_{\parallel }\cdot \overset{\mbox{ˆ}}{\nabla }w.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 55)
|}

El problema discretizado toma finalmente la forma de: ''hallar  ''<math display="inline">{\eta }^h\in S^h</math>''tal que''

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a\left(w^h,{\eta }^h\right)+c\left(\boldsymbol{v}_{\parallel };w^h,{\eta }^h\right)+</math><math>\sum_e{\int }_{{\Omega }_e^{{'}}}P(w^h)\quad {\tau }_S\quad R({\eta }^h)\quad d{\Omega }_e^{{'}}=</math><math>\left(w^h,s\right);</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 56)
|}

''para todo  ''<math display="inline">w^h\in W^h</math> , con <math display="inline">{\eta }^h\in S^h</math> , siendo <math display="inline">W^h</math>  y <math display="inline">S^h</math>  subespacios de dimensión finita de <math display="inline">W</math>  y <math display="inline">S</math> , respectivamente.      

El tiempo intrínseco ''τ''<sub>''S''</sub>  es tomado como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\tau }_S=\frac{h}{2\parallel \boldsymbol{a}\parallel };</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 57)
|}

en la cual ''h''  es un tamaño típico de la discretización en la SL y '''''a'''''  es una velocidad media, ambos en cada elemento, y considerando que en general las velocidades <math display="inline">\boldsymbol{v}_{\parallel }</math>  pueden ser distintas en cada nodo. Cuando la velocidad media '''''a'''''  es baja, no es necesario realizar una estabilización numérica, de manera que se adopta ''τ''<sub>S</sub>  = 0 cuando '''''a'''''  se encuentra por debajo de una velocidad umbral, que puede determinarse en función de dicha velocidad, del paso de tiempo y del tamaño medio de los elementos de la malla.      

===2.4. Movimiento de la malla===

En un paradigma ALE, en el que el dominio sufre deformaciones al evolucionar el problema en el tiempo, es preciso contar con una metodología para actualizar el teselado, generalmente en todos los pasos de tiempo, siendo posible optar, entre otras, por alguna de las alternativas siguientes [[#bib0025|[5]]] :
* '''Actualización algebraica''' : la discretización es modificada mediante expresiones algebraicas explícitas para el desplazamiento de los nodos interiores como función de los desplazamientos de la SL, conservando la definición topológica de los elementos, tal como en los métodos de ''spines''[[#bib0170|[34]]]  sobre mallas estructuradas.                          
* '''Relocalización de nodos interiores''' : consiste en la reubicación de los nodos interiores del dominio manteniendo la topología, pero a través de procedimientos auxiliares tales como la resolución de un problema pseudoelástico [[#bib0015|[3]]] , [[#bib0025|[5]]] , [[#bib0070|[14]]] , [[#bib0105|[21]]]  and [[#bib0160|[32]]]  o la minimización de un indicador de la distorsión de los elementos [[#bib0110|[22]]] .                          
* '''Remallado''' : la malla es generada nuevamente cada vez que el dominio es deformado, lo cual exige la interpolación/extrapolación de los valores de los campos en los nodos.                  

La alternativa de relocalización de los nodos de la discretización es la adoptada en este trabajo, en la variante de resolución del problema pseudoelástico. El procedimiento consiste en calcular las nuevas posiciones de los nodos resolviendo un problema elástico artificial sobre el dominio inicial ''Ω''<sub>0</sub> , en el cual las condiciones de contorno son siempre de tipo Dirichlet. Este problema pseudoelástico se formula como uno elástico convencional, esto es,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\sigma }_{ij,j}=0;</math>
|-
|<math>{\sigma }_{ij}=2\tilde{\mu }{\epsilon }_{ij}+\tilde{\lambda }{\delta }_{ij}{\epsilon }_{kk};</math>
|-
|<math>{\epsilon }_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i});</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 58)
|}

donde <math display="inline">\tilde{\mu }</math>  y <math display="inline">\tilde{\lambda }</math>  son las constantes elásticas de Lamé adoptadas arbitrariamente para el material, ''δ''<sub>''ij''</sub>  es el tensor de Kronecker y los desplazamientos nodales

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}_j=\boldsymbol{x}_j^{n+1}-\boldsymbol{x}_j^0;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 59)
|}

corresponden a las condiciones de contorno sobre la SL. Para contornos sólidos, las condiciones de contorno se imponen de la forma '''''u'''''  = '''''0'''''  para los sectores no deslizables y '''''u'''''  · '''''n'''''  = 0 para los de deslizamiento perfecto.      

Las propiedades artificiales del material pueden expresarse en términos del coeficiente de Poisson <math display="inline">\tilde{\nu }</math>  y del módulo de elasticidad longitudinal <math display="inline">\tilde{E}</math> , que son los parámetros tomados como referencia a la hora de ingresar los datos para un cálculo dado, aunque debido al tipo de condiciones de contorno, solo el primero de ellos tiene relevancia, recordando que el objetivo es calcular posiciones de nodos y no tensiones elásticas. Normalmente, se adopta <math display="inline">\tilde{\nu }=</math><math>0,3</math>  teniendo en cuenta que para <math display="inline">\tilde{\nu }\rightarrow 0,5</math> , es decir, tendiendo a incompresibilidad, el problema numérico está mal condicionado.      

Una vez que se ha determinado el campo de velocidades sobre la SL, los desplazamientos nodales sobre esta son calculados y pasan a ser datos para resolver el problema de actualización de la malla.

En lo relativo a condiciones de contorno, por ejemplo, en un tanque rectangular con una SL simple consisten en: (i) nodos con desplazamientos nulos en el fondo del tanque; (ii) nodos fijos en dirección horizontal pero libres de desplazarse verticalmente en las paredes laterales, y (iii) los desplazamientos de los nodos de la SL, determinados según lo indicado en la sección anterior. Obviamente, la selección de las condiciones de contorno debe ser acorde al problema considerado.

Esta alternativa de relocalización de nodos permite abordar problemas en los cuales las deformaciones son relativamente grandes, aunque puede proporcionar mallas distorsionadas de forma tal que los resultados numéricos se ven deteriorados y, eventualmente, puede llegar a fallar la actualización debido a la generación de elementos con jacobiano negativo o nulo.

==3. Ejemplos de aplicación==

===3.1. Agitación de pequeña amplitud en un dominio bidimensional===

El método propuesto es validado mediante un problema bidimensional de agitación de un fluido viscoso con solución analítica, empleado por Rabier y Medale [[#bib0160|[32]]]  y Battaglia et al. [[#bib0015|[3]]] . Este test muestra cómo el esquema numérico propuesto predice, además de la frecuencia, la tasa de amortiguamiento dominada por la viscosidad del fluido.      

El ejemplo consiste en resolver el problema de valores iniciales del movimiento de pequeña amplitud de la SL de un fluido viscoso en un tanque rectangular esquematizado en la [[#fig0015|figura 3]]  cuya SL tiene una posición inicial dada por

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h\left(x\right)=1,5\mbox{ }m\mbox{ }+\mbox{ }a_0sin\left[\pi \left(1/2-\right. \right. </math><math>\left. \left. x\right)\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 60)
|}

donde ''a''<sub>0</sub>  es la amplitud de la perturbación sinusoidal inicial del movimiento. El fluido se encuentra sometido a la aceleración de la gravedad g, ''y''  las fuerzas viscosas son responsables del amortiguamiento del movimiento.

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr3.jpg|center|251px|Dimensiones en metros y posición inicial de la SL para el problema de agitación ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Dimensiones en metros y posición inicial de la SL para el problema de agitación con solución analítica.

</span>
|}

La solución analítica del caso linealizado está dada por Prosperetti [[#bib0155|[31]]]  en la forma

<span id='eq0340'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>a(t)=\frac{4{\nu }^2k^4}{8{\nu }^2k^4+{{\omega }_0}^2}\quad a_0\quad erfc(\nu k^2t)^{1/2}+</math><math>\sum_{i=1}^4\frac{z_i}{Z_i}\left(\frac{{\omega }_0^2a_0}{z_i^2-\nu k^2}\right)exp[(z_i^2-</math><math>\nu k^2)t]\quad erfc(z_it^{1/2});</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 61)
|}

en la cual ''ν''  es la viscosidad cinemática del fluido, ''k  ''  es el número de onda, <math display="inline">{\omega }_0^2=</math><math>gk</math>  es la frecuencia natural para el caso no viscoso, y cada ''z''<sub>''i''</sub>  es una raíz de la siguiente ecuación algebraica:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>z^4+k^2\nu z^2+4(k^2\nu )^{3/2}z+{\nu }^2k^4+{\omega }_0^2=</math><math>0;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 62)
|}

donde ''Z''<sub>1</sub>  = (''z''<sub>2</sub>  − ''z''<sub>1</sub> )(''z''<sub>3</sub>  − ''z''<sub>1</sub> )(''z''<sub>4</sub>  − ''z''<sub>1</sub> ), ''Z''<sub>2</sub> , ''Z''<sub>3</sub> , ''Z''<sub>4</sub>  se obtienen por permutación circular de los índices y erfc (…) es la función error para variable compleja. Notar que la ecuación [[#eq0340|(61)]]  transcripta en [[#bib0160|[32]]]  a partir de [[#bib0155|[31]]]  tiene un error de tipeo.      

El ejemplo fue resuelto para una amplitud inicial de la deformación ''a''<sub>0</sub>  = 0, 01 m, aceleración gravitatoria unitaria ''g''  = − 1, 0 m/''s''<sup>2</sup>  y una malla de 40 × 60 elementos cuadrangulares para un tanque representado en la [[#fig0015|figura 3]]  de ''h''  = 1, 5 m de altura y ''d''  = 1, 0 m de ancho, tal que ''k''  = 2''d'' , con un paso de tiempo Δ''t''  = 2, 12 × 10<sup>−2</sup>  s. Las propiedades físicas del fluido fueron adoptadas para un primer análisis de la resolución de este mismo ejemplo en la literatura  [[#bib0030|[6]]] , [[#bib0160|[32]]]  and [[#bib0165|[33]]] , en el cual la viscosidad cinemática y la densidad son ''ν''  = 10<sup>−3</sup>  m<sup>2</sup> /''s''  y ''ρ''  = 1 kg/m<sup>3</sup> , respectivamente.      

La [[#fig0020|figura 4]]  muestra la curva de posición vertical del nodo superior izquierdo versus el tiempo superpuesta a la curva calculada analíticamente, donde se observa que tanto la frecuencia como el amortiguamiento viscoso del sistema son bien capturados por la solución numérica, siendo el error relativo menor al 3% en relación con la amplitud inicial ''a''<sub>0</sub> . Las condiciones de contorno elegidas son de deslizamiento perfecto sobre los contornos ''Γ''<sub>D</sub> , lo cual implica que la velocidad en dirección normal a estos y las tensiones tangenciales son nulas, tanto en los límites laterales como en el fondo del tanque. La condición de deslizamiento perfecto sobre las paredes responde a la hipótesis de dominio lateral e inferior infinitos para la solución analítica [[#bib0155|[31]]] . Sobre la SL, ''P''<sub>atm</sub>  = 0 y las tensiones tangenciales son '''''T'''''  · '''''n'''''  = '''''0''''' .

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr4.jpg|center|346px|Curva de solución analítica y resultados numéricos para el problema de agitación ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

Curva de solución analítica y resultados numéricos para el problema de agitación con viscosidad cinemática ''ν''  = 10<sup>−3</sup>  m<sup>2</sup> /''s'' .                  

</span>
|}

Las velocidades que se registran sobre la SL son de muy baja intensidad y de dirección prácticamente vertical, por tanto no se producen inestabilidades numéricas debidas a efectos convectivos y no es necesario realizar una estabilización. Independientemente de ello, la inclusión de los términos correspondientes a la estabilización en la formulación de la ecuación [[#eq0285|(50)]]  no influirían sensiblemente en los resultados, pues las velocidades <math display="inline">\boldsymbol{v}_{\parallel }</math>  son menores que las velocidades que intervienen en el término fuente ''s'' , en un factor del orden de 10<sup>−3</sup> .      

A los fines de la comprobación del método, el problema fue resuelto con diferentes valores de viscosidad cinemática, manteniendo el resto de los parámetros sin cambios. Así, para ''ν''  = 10<sup>−4</sup>  m<sup>2</sup> /''s''  se obtuvieron curvas con menor amortiguamiento que en el primer caso analizado y un error relativo de menos del 4% con respecto a ''a''<sub>0</sub> , en tanto que en la resolución con ''ν''  = 10<sup>−2</sup>  m<sup>2</sup> /''s'' , representada en la  [[#fig0025|figura 5]] , el error relativo a la amplitud inicial ''a''<sub>0</sub>  es del 2%.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr5.jpg|center|345px|Curva de solución analítica y resultados numéricos para el problema de agitación ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 5.

Curva de solución analítica y resultados numéricos para el problema de agitación con viscosidad cinemática ''ν''  = 10<sup>−2</sup>  m<sup>2</sup> /''s'' .                  

</span>
|}

La resolución numérica de este problema se realizó en un ordenador de escritorio con procesador ''Pentium IV''  y 2 Gb de memoria RAM, totalizando 56 minutos de ejecución para 300 pasos de tiempo. Del tiempo total de ejecución, el 25% fue insumido por los pasos de actualización de la malla. En lo relativo a la convergencia del método, el módulo de NS requirió 3 iteraciones de Newton y el de MMV, de elasticidad lineal, 2 iteraciones, considerando en ambos casos una tolerancia absoluta de 1 × 10<sup>−8</sup>  para la norma del residuo.      

===3.2. Agitación en un dominio cilíndrico inclinado===

Otro ejemplo consiste en resolver la agitación que se produce en un tanque cilíndrico de radio ''R''  = 0, 50 m y altura de líquido ''H''<sub>g</sub>  = 0,60 m de eje inclinado en un ángulo ''γ''  = 10<sup>''o''</sup>  con respecto a la aceleración gravitatoria '''''g''''' , como puede apreciarse en la [[#fig0030|figura 6]] . El movimiento del fluido se genera a partir de liberar la SL desde una posición inicial diferente de la hidrostática.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr6.jpg|center|133px|Geometría para el problema de agitación en un dominio cilíndrico de eje ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

Geometría para el problema de agitación en un dominio cilíndrico de eje inclinado.

</span>
|}

Con el objetivo de verificar la aptitud del método para resolver este tipo de problemas, se establece una comparación de resultados con la solución analítica no viscosa dada por Moiseev y Petrov [[#bib0115|[23]]] . Para ello, la posición inicial de la SL se propone proporcional a la deformada del primer modo de oscilación propio del sistema, cuyo autovalor adimensional es <math display="inline">{\kappa }_s^n\approx 1.746</math>  para ''γ''  = 10<sup>''o''</sup> , modo de orden ''n''  = 1 y número de autovalor ''s''  = 1. Dado el autovalor, la frecuencia de oscilación ''σ''  se calcula como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sigma ={\left(\frac{g\quad {\kappa }_1^{\left(1\right)}}{R}\right)}^{0,5}\approx 5,853\quad s^{-1};</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 63)
|}

con ''g''  = 9, 81 m/s<sup>2</sup>  la magnitud de la aceleración gravitatoria, tal que el período de oscilación resulta:

<span id='eq0355'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>T=\frac{2\pi }{\sigma }\approx 1,0735\quad s;</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 64)
|}

La resolución numérica del ejemplo fue realizada sobre una malla de 28.800 elementos finitos hexahédricos y 31.000 nodos mostrada en la [[#fig0035|figura 7]] . La discretización fue propuesta tal que el eje del cilindro inclinado coincide con la dirección espacial ''x''<sub>3</sub> , en tanto que la aceleración gravitatoria se impuso con componentes ''g''<sub>1</sub>  = ''g''  sin (''γ'' ), ''g''<sub>2</sub>  = 0 y ''g''<sub>3</sub>  = ''g''  cos (''γ'' ) en las direcciones ''x''<sub>1</sub> , ''x''<sub>2</sub>  y ''x''<sub>3</sub> , respectivamente. El fluido contenido en el recipiente es glicerina, cuyas propiedades físicas son: viscosidad cinemática ''ν''  = 6, 5<sup>−4</sup>  m<sup>2</sup> /''s''  y densidad ''ρ''  = 1.260 kg/''m''<sup>3</sup> . El paso de tiempo adoptado para el análisis es Δ''t''  = 0, 025 s, hasta un tiempo final de 5 s.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr7.jpg|center|226px|Malla de elementos finitos para el problema del cilindro de eje inclinado.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

Malla de elementos finitos para el problema del cilindro de eje inclinado.

</span>
|}

Las condiciones de borde para el problema de NS sobre la SL vienen dadas por un valor de presión ''P''<sub>atm</sub>  = 0 y tensiones tangenciales '''''T'''''  · '''''n'''''  = '''''0''''' . En la base y en el 70% inferior de las paredes del cilindro se propuso una condición de no deslizamiento <math display="inline">\boldsymbol{v}=</math><math>\boldsymbol{0}</math> , mientras que en el 30% superior de la superficie cilíndrica se adoptó una condición de deslizamiento perfecto a los fines de disminuir los efectos de la viscosidad en las adyacencias de la SL, dado el carácter no viscoso de la solución de referencia. En lo que respecta al problema de movimiento de malla, en la base y en el 70% inferior de la superficie cilíndrica los nodos se encuentran impedidos de desplazarse en cualquier dirección, en tanto que el 30% superior del contorno cilíndrico puede desplazarse libremente sobre la superficie lateral. Para la SL, los desplazamientos nodales son determinados en función de las velocidades obtenidas en NS. La máxima amplitud inicial de la deformación de la SL por sobre el plano hidrostático se adoptó de valor ''a''<sub>0</sub>  = 0, 012 m, considerando al fluido en reposo al comienzo del análisis.      

Una vez resuelto numéricamente el problema, se determinó el valor del período de oscilación calculado, en este caso en función de una pseudoamplitud, que consiste en la suma de los productos del desplazamiento por la posición ''x''<sub>1</sub>  de cada nodo de la superficie libre, normalizado con el radio ''R'' , representada en la  [[#fig0040|figura 8]] . El período de oscilación calculado con la pseudoamplitud fue de ''T''  = 1.0768 s, esto es, con un error relativo de 0,30% con respecto al período de oscilación para el problema no viscoso estimado mediante la ecuación [[#eq0355|(64)]] .

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr8.jpg|center|265px|Pseudoamplitud de desplazamientos en el ejemplo del tanque cilíndrico de eje ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

Pseudoamplitud de desplazamientos en el ejemplo del tanque cilíndrico de eje inclinado.

</span>
|}

La [[#fig0045|figura 9]]  muestra los desplazamientos de los nodos extremos en dirección ''x''<sub>1</sub>  sobre la SL, en la cual se indica con 1 el nodo en ''x''<sub>1</sub>  positivo, es decir, donde la deformación de la SL es mayor medida en dirección ''x''<sub>3</sub> , en tanto que 2 corresponde al nodo restante, de menor amplitud inicial. En dicha figura puede apreciarse la disminución progresiva de la amplitud de las deformaciones debida a los efectos viscosos, así como también la diferencia de amplitudes entre desplazamientos positivos y negativos debida a la acumulación de fluido que se produce cuando la interfase asciende sobre los laterales.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr9.jpg|center|265px|Desplazamientos de los nodos extremos en x1 sobre la SL del tanque cilíndrico de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

Desplazamientos de los nodos extremos en ''x''<sub>1</sub>  sobre la SL del tanque cilíndrico de eje inclinado.                  

</span>
|}

El problema fue resuelto numéricamente en 10 nodos del ''cluster''  Aquiles  [[#bib0005|[1]]] , empleando para ello 64 min para la resolución de 200 pasos de tiempo. La cantidad de iteraciones del lazo de Newton en cada paso de tiempo fue de 4, tanto para el problema de fluido como para el de actualización de la malla, con una tolerancia de 1 × 10<sup>−8</sup> . La proporción de tiempo empleada por el módulo de movimiento de malla fue del 30 al 35% del total de cada paso de tiempo.      

===3.3. Agitación producida por aceleración horizontal periódica===

La metodología descripta es aplicada a la agitación producida en un contenedor rectangular de ''W''  = 0, 80 m de ancho, con una altura de agua de ''D''  = 0, 30 m, como puede verse en la [[#fig0050|figura 10]] , y es modelado como un problema bidimensional. Las condiciones de contorno para el problema de NS en las paredes sólidas se propusieron como de deslizamiento perfecto debido a la baja viscosidad del fluido, ''ν''  = 10<sup>−6</sup>  m<sup>2</sup> /s, tal como en la bibliografía consultada  [[#bib0090|[18]]] , [[#bib0195|[39]]]  and [[#bib0200|[40]]] . La aceleración que actúa sobre el dominio es ''G''  = [''g''<sub>1</sub>  ; ''g''<sub>2</sub> ]<sup>''T''</sup> , de componente vertical de magnitud ''g''<sub>2</sub>  = − ''g''  = − 9, 81 m/s<sup>2</sup> , y componente horizontal dependiente del tiempo, propuesta de variación sinusoidal. El valor de densidad adoptado es de ''ρ''  = 1.000 kg/''m''<sup>3</sup> .

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr10.jpg|center|263px|Datos geométricos y condiciones de contorno para el tanque sometido a ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 10.

Datos geométricos y condiciones de contorno para el tanque sometido a aceleraciones horizontales, tanto periódicas como de origen sísmico.

</span>
|}

Este ejemplo fue tomado de Huerta y Liu [[#bib0090|[18]]] , y es reproducido por Soulaïmani et al. [[#bib0195|[39]]] , entre otros. Aunque la formulación considerada en el presente trabajo está enfocada a desplazamientos pequeños de la SL, este ejemplo permite mostrar la robustez del método al momento de resolver desplazamientos mayores y velocidades horizontales con intensidad mayor que las que se producen en los problemas antes presentados. La aceleración horizontal se propone como ''g''<sub>1</sub>  = ''A''  ''g''  sin ''ωt'' , en la cual el coeficiente de amplitud es ''A''  = 0, 01, ''t''  es el tiempo y ''ω''  la frecuencia circular. Esta frecuencia ha sido calculada de manera tal que excita el primer modo de agitación, teniendo en cuenta que la longitud de onda para este es ''λ''  = 2''W'' , por lo tanto,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>f=\sqrt{\frac{g}{4\pi W}\quad tanh\frac{\pi D}{W}}\approx 0,89825\quad 1/s.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 65)
|}

La frecuencia angular a imponer es entonces ''ω''  = 2''πf''  ≈ 5, 64 rad/s.      

La [[#fig0055|figura 11]]  muestra el desplazamiento del nodo superior izquierdo comparado con las envolventes de los resultados de la bibliografía [[#bib0090|[18]]]  and [[#bib0195|[39]]] , para el problema resuelto sobre una malla de 39 elementos en profundidad y 104 en ancho, con Δ''t''  = 0,009 s, haciendo uso de un esquema de integración temporal Crank-Nicolson, i.e. ''α''  = 0, 5, a lo largo de 2.000 pasos de tiempo. Teniendo en cuenta que el tanque es excitado con la frecuencia asociada al primer modo de ''sloshing'' , se espera un incremento en la amplitud del movimiento hasta que o bien se equilibren la potencia aplicada y la disipación viscosa, o bien se produzca la rotura de la SL, invalidando el método.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr11.jpg|center|350px|Desplazamiento vertical del nodo izquierdo sobre la SL, relativo a la ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 11.

Desplazamiento vertical del nodo izquierdo sobre la SL, relativo a la profundidad del tanque con integración temporal tipo Crank-Nicolson.

</span>
|}

En la [[#fig0055|figura 11]]  puede apreciarse que los desplazamientos en subida son mayores que las amplitudes en bajada, en coincidencia con resultados de modelos físicos [[#bib0090|[18]]] . Esto se debe a que el fluido en la zona de descenso muestra curvas más suaves que en ascenso, donde el líquido adquiere una forma empinada, como puede verse en la [[#fig0060|figura 12]] .

<span id='fig0060'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr12.jpg|center|310px|Magnitud de las velocidades y deformación del dominio en el problema del tanque ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 12.

Magnitud de las velocidades y deformación del dominio en el problema del tanque agitado horizontalmente en ''t''  = 16, 11 s.                  

</span>
|}

Se efectuó un estudio de convergencia en malla del presente método numérico, para el cual se realizaron análisis con las diversas discretizaciones que se detallan en la [[#tbl0005|tabla 1]] . Además, se calcularon y graficaron los valores de error absoluto ''E''  en función del tamaño de los elementos ''h'' , mostrados en la  [[#fig0065|figura 13]]  y calculados para los primeros 4 seg de la simulación. Se tomó como solución de referencia para la determinación del error el resultado del análisis obtenido con la malla indicada con 1 en la [[#tbl0005|tabla 1]] , que corresponde al mayor grado de refinamiento. Aproximando los puntos por mínimos cuadrados se obtiene una tasa de convergencia aproximada ''E''  ∝ ''h''<sup>1,67</sup> .

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Datos para el estudio de convergencia en el problema con aceleración horizontal periódica.

|-

! 
! colspan="2" | Cantidad elementos
! colspan="2" | Parámetros de la malla
! Error absoluto
|-

! Malla
! horizontal
! vertical
! ''h''  [m]                                                    
! Δ''t''  [s]                                                    
! [m]
|-

| 1 (referencia)
| 208
| 78
| 0,0038
| 0,0045
| –
|-

| 2 (resultados)
| 104
| 39
| 0,0077
| 0,0090
| 0,0050
|-

| 3
| 78
| 29
| 0,0103
| 0,0135
| 0,0106
|-

| 4
| 52
| 20
| 0,0154
| 0,0180
| 0,0163
|}

<span id='fig0065'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr13.jpg|center|248px|Convergencia en malla para el problema del tanque sometido a agitación ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 13.

Convergencia en malla para el problema del tanque sometido a agitación horizontal periódica, con escala logarítmica en abscisas y ordenadas. En los ejes derecho y superior se indican las magnitudes correspondientes de ''h''  y de ''E'' .                  

</span>
|}

La resolución numérica sobre la malla más refinada, de 208 × 78 elementos y para 900 pasos de tiempo se realizó en 6 hs : 25 min en un ordenador de escritorio con procesador ''Pentium IV''  y 2 Gb de memoria RAM, en tanto que en el mismo equipo se requirieron 2 hs : 31 min para resolver 500 pasos de tiempo en la malla de 104 × 39 elementos. Las características de convergencia para las dos instancias de ejecución, NS y MMV, coinciden con el ejemplo de agitación de pequeña amplitud en el dominio bidimensional antes presentado.      

===3.4. Agitación producida por aceleraciones sísmicas===

El ejemplo consiste en resolver un problema sobre la misma geometría del ejemplo de aceleración horizontal periódica, [[#fig0050|figura 10]] , imponiendo las aceleraciones horizontales registradas durante un sismo. Tanto la discretización espacial como los parámetros del fluido son los mismos que en el ejemplo anterior, esto es, 104 × 39 elementos en ancho y en profundidad, respectivamente, viscosidad cinemática ''ν''  = 10<sup>−6</sup>  m<sup>2</sup> /s y densidad ''ρ''  = 1.000 kg/''m''<sup>3</sup> . El paso de tiempo adoptado coincide con la frecuencia de registro de aceleraciones, es decir Δ''t''  = 0, 005 s.      

El acelerograma empleado corresponde al terremoto del día 3 de marzo de 1985 en Valparaíso (Chile), a las 19:47 hora local, medido en la estación La Ligua [[#bib0040|[8]]] . Se trata de un movimiento sísmico de algo más de 48 s de duración con registros cada 0, 005 s y amplitudes de aceleración horizontal mostradas en la [[#fig0070|figura 14]]  (arriba). En dicha figura se aprecia además la respuesta en desplazamientos relativos a la altura de agua ''D''  = 0, 30 m de los nodos superior izquierdo, [[#fig0070|figura 14]]  (centro), y superior derecho, [[#fig0070|figura 14]]  (abajo), mostrando amplitudes relativas menores que las del caso anterior, lo cual se debe a que las frecuencias predominantes en el acelerograma no son cercanas a modos propios de baja frecuencia del recipiente.

<span id='fig0070'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_524553912-1-s2.0-S0213131512000119-gr14.jpg|center|528px|Acelerograma (arriba) y desplazamiento relativo de los nodos extremos del tanque ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 14.

Acelerograma (arriba) y desplazamiento relativo de los nodos extremos del tanque sometido a aceleraciones sísmicas (centro y abajo).

</span>
|}

El problema insumió 20 hs : 20 min de tiempo para la ejecución de 9.740 pasos del registro de aceleraciones horizontales del movimiento sísmico, empleando para ello un nodo del ''cluster''  Aquiles  [[#bib0005|[1]]] . Al igual que en el ejemplo anterior, las características de la convergencia tanto del NS como del MMV coinciden con el caso de agitación de pequeña amplitud en un dominio bidimensional.      

==4. Conclusiones==

Se ha mostrado el desempeño de una estrategia lagrangiana-euleriana arbitraria en la resolución de problemas de agitación en tanques de almacenamiento o transporte de líquidos en dos y tres dimensiones espaciales, para desplazamientos pequeños o moderados de la SL, resolviendo mediante cálculo distribuido las instancias involucradas en cada paso de tiempo, en un esquema de acoplamiento débil. La validez del método ha sido mostrada a través de la resolución de ejemplos con solución analítica, así como también la capacidad de resolver problemas de agitación producidos tanto por deformaciones iniciales de la SL como por aceleraciones impuestas a lo largo del análisis, como en el caso del tanque excitado horizontalmente. Se verificó además la convergencia en malla, a través de la comparación de los errores absolutos medidos en análisis realizados con distintas discretizaciones. La incorporación de la estabilización numérica en la resolución de la ecuación de transporte de la SL introduce una mejora en relación con la metodología presentada con anterioridad [[#bib0015|[3]]] , que no resultaba apropiada para resolver problemas en los cuales la interfase se ve sometida a velocidades tangenciales de importancia o a desplazamientos moderados.

==Agradecimientos==

Los autores desean agradecer la financiación del trabajo mediante proyectos del Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) de Argentina, PIP–5271/05, de la Universidad Nacional del Litoral (UNL), CAI+D 2009–III–4–2 y CAI+D 2009 65/33 y de la Agencia Nacional de Promoción Científica y Tecnológica (ANPCyT) de Argentina, PICT 1506/06, PICT 1141/07.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] Cluster Aquiles en CIMEC. http://www.cimec.org.ar/aquiles , 2011.                                    </li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] S. Balay, K. Buschelman, V. Eijkhout, W. Gropp, D. Kaushik, M. Knepley, L.C. McInnes, B. Smith, H. Zhang; PETSc Users Manual. ANL 95/11 - Revision 3.0.0; Argonne National Laboratory (2008)</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] L. Battaglia, J. D’Elía, M.A. Storti, N.M. Nigro; Numerical simulation of transient free surface flows using a moving mesh technique; ASME J. Appl. Mech., 73 (6) (2006), pp. 1017–1025</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] L. Battaglia, J. D’Elía, M.A. Storti, N.M. Nigro; Parallel implementations of free surface flows; G. Buscaglia, E. Dari, O. Zamonsky (Eds.), Mec. Comp., volumen XXIII, San Carlos de Bariloche, Noviembre 08-11 (2004), pp. 3119–3132</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] M. Behr, F. Abraham; Free surface flow simulations in the presence of inclined walls; Comp. Meth. Appl. Mech. Eng., 191 (47-48) (2002), pp. 5467–5483</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] H. Braess, P. Wriggers; Arbitrary lagrangian eulerian finite element analysis of free surface flow; Comp. Meth. Appl. Mech. Eng., 190 (1-2) (2000), pp. 95–109</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] A.N. Brooks, T.J.R. Hughes; Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations; Comp. Meth. Appl. Mech. Eng., 32 (1-3) (1982), pp. 199–259</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] Departamento de Geofísica y Geodesia de Chile. Acelerogramas del terremoto de Valparaíso (Chile) en 1985. COSMOS Virtual Data Center, http://db.cosmos-eq.org , 2009.                                    </li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] J. D’Elía, M. Storti, S. Idelsohn; Un método de paneles para el cálculo de la resistencia de ola en barcos; Rev. int. métodos numér. cálc. diseño ing., 13 (4) (1997), pp. 515–530</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] J. Donea, A. Huerta; Finite Element Methods for Flow Problems; John Wiley & Sons (2003)</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] O.M. Faltinsen, O.F. Rognebakke, A.N. Timokha; Classification of three-dimensional non-linear sloshing in a square-base tank with finite depth; J. Fluids Struct., 20 (2005), pp. 81–103</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] A. Folch Duran. A numerical formulation to solve the ALE Navier-Stokes equations applied to the withdrawal of magma chambers. Tesis doctoral, Universitat Politécnica de Catalunya, Abril 2000.</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] J. García Espinosa. Un Método de Elementos Finitos para Análisis Hidrodinámico de Estructuras Navales. Tesis doctoral, Universitat Politécnica de Catalunya, 1999.</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] I. Güler, M. Behr, T.E. Tezduyar; Parallel finite element computation of free-surface flows; Comp. Mech., 23 (2) (1999), pp. 117–123</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] F.H. Harlow, J.E. Welch; Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface; Phys. Fluids, 8 (12) (1965), pp. 2182–2189</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] H. Hernández-Barrios, E. Heredia-Zavoni, A.A. Aldama-Rodríguez; Nonlinear sloshing response of cylindrical tanks subjected to earthquake ground motion; Eng. Struct., 29 (12) (2007), pp. 3364–3376</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] C.W. Hirt, B.D. Nichols; Volume of fluid (VOF) method for the dynamics of free boundaries; J. Comput. Phys., 39 (1) (1981), pp. 201–225</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] A. Huerta, W.K. Liu; Viscous flow with large free surface motion; Comp. Meth. Appl. Mech. Eng., 69 (1988), pp. 277–324</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] T.J.R. Hughes, W.K. Liu, T.K. Zimmermann; Lagrangian-Eulerian finite element formulation for incompressible viscous flows; Comp. Meth. Appl. Mech. Eng., 29 (3) (1981), pp. 329–349</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] S.R. Idelsohn, E. Oñate, F. Del Pin; The Particle Finite Element Method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves; Int. J. for Num. Meth. Eng., 61 (7) (2004), pp. 964–984</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] A.A. Johnson, T.E. Tezduyar; Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces; Comp. Meth. Appl. Mech. Eng., 119 (1-2) (1994), pp. 73–94</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] E.J. López, N.M. Nigro, M.A. Storti; Simultaneous untangling and smoothing of moving grids; Int. J. for Num. Meth. Eng., 76 (7) (2008), pp. 994–1019</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] N.N. Moiseev, A.A. Petrov; The calculation of free oscillations of a liquid in a motionless container;  Dryden,  Von Kármán (Eds.), Advances in Applied Mechanics, volumen 9, Academic Press (1966), pp. 91–155</li>
<li><span id='bib0120'></span>
[[#bib0120|[24]]] MPI. Message Passing Interface. http://www.mpi-forum.org , 2008.                                    </li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] A.D. Myshkis, V.G. Babskii, N.D. Kopachevskii, L.A. Slobozhanin, A.D. Tyuptsov; Low-Gravity Fluid Mechanics; Springer-Verlag (1987)</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] C.M. Oishi, M.F. Tome, J.A. Cuminato, S. McKee; An implicit technique for solving 3D low Reynolds number moving free surface flows; J. Comput. Phys., 227 (16) (2008), pp. 7446–7468</li>
<li><span id='bib0135'></span>
[[#bib0135|[27]]] S. Osher, J.A. Sethian; Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations; J. Comput. Phys., 79 (1) (1988), pp. 12–49</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] S. Papaspyrou, D. Valougeorgis, S.A. Karamanos; Sloshing effects in half-full horizontal cylindrical vessels under longitudinal excitation; ASME: J. Appl. Mech., 71 (2) (2004), pp. 255–265</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] J. Pedlosky; Waves in the ocean and atmosphere. Introduction to wave dynamics; Springer (2003)</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] PETSc-FEM. A general purpose, parallel, multi-physics FEM program, 2011. GNU General Public License (GPL), http://www.cimec.org.ar/petscfem .                                    </li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] A. Prosperetti; Motion of two superposed viscous fluids; Phys. Fluids, 24 (7) (1981), pp. 1217–1223</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] S. Rabier, M. Medale; Computation of free surface flows with a projection FEM in a moving mesh framework; Comp. Meth. Appl. Mech. Eng., 192 (41-42) (2003), pp. 4703–4721</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] B. Ramaswamy; Numerical simulation of unsteady viscous free surface flow; J. Comput. Phys., 90 (2) (1990), pp. 396–430</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] H. Saito, L.E. Scriven; Study of coating flow by the finite element method; J. Comput. Phys., 42 (1) (1981), pp. 53–76</li>
<li><span id='bib0175'></span>
[[#bib0175|[35]]] R. Scardovelli, S. Zaleski; Direct numerical simulation of free-surface and interfacial flow; Annu. Rev. Fluid Mech., 31 (1999), pp. 567–603</li>
<li><span id='bib0180'></span>
[[#bib0180|[36]]] J.A. Sethian; A Fast Marching Level Set Method for Monotonically Advancing Fronts; National Academy of Sciences (1995)</li>
<li><span id='bib0185'></span>
[[#bib0185|[37]]] W. Shyy, H.S. Udaykumar, M.M. Rao, R.W. Smith; Computational Fluid Dynamics with Moving Boundaries; Taylor and Francis (1996)</li>
<li><span id='bib0190'></span>
[[#bib0190|[38]]] V.E. Sonzogni, A.M. Yommi, N.M. Nigro, M.A. Storti; A parallel finite element program on a Beowulf Cluster; Adv. Eng. Software, 33 (7-10) (2002), pp. 427–443</li>
<li><span id='bib0195'></span>
[[#bib0195|[39]]] A. Soulaïmani, M. Fortin, G. Dhatt, Y. Ouellet; Finite element simulation of two- and three-dimensional free surface flows; Comp. Meth. Appl. Mech. Eng., 86 (3) (1991), pp. 265–296</li>
<li><span id='bib0200'></span>
[[#bib0200|[40]]] M. Souli, J.P. Zolesio; Arbitrary Lagrangian-Eulerian and free surface methods in fluid mechanics; Comp. Meth. Appl. Mech. Eng., 191 (3-5) (2001), pp. 451–466</li>
<li><span id='bib0205'></span>
[[#bib0205|[41]]] J.J. Stoker; Water waves. The mathematical theory with applications, volumen IV of Pure and Applied Mathematics; Interscience (1957)</li>
<li><span id='bib0210'></span>
[[#bib0210|[42]]] M. Storti, J. D’Elía, S. Idelsohn; Forma simétrica de la condición de contorno absorbente dnl para el problema de la resistencia de ola en barcos; Rev. int. métodos numér. cálc. diseño ing., 14 (3) (1998), pp. 317–338</li>
<li><span id='bib0215'></span>
[[#bib0215|[43]]] P.R.F. Teixeira, C.J. Fortes; Aplicação de um modelo numérico não-linear tridimensional na análise de ondas sobre quebra-mares trapezoidais submersos; Rev. int. métodos numér. cálc. diseño ing., 25 (4) (2009), pp. 313–336</li>
<li><span id='bib0220'></span>
[[#bib0220|[44]]] T.E. Tezduyar; Stabilized finite element formulations for incompressible flow computations; Adv. Appl. Mech., 28 (1991), pp. 1–44</li>
<li><span id='bib0225'></span>
[[#bib0225|[45]]] T.E. Tezduyar, S. Mittal, S.E. Ray, R. Shih; Incompressible flow computations with stabilized bilinear and linear-equal-order interpolation velocity-pressure elements; Comp. Meth. Appl. Mech. Eng., 95 (2) (1992), pp. 221–242</li>
<li><span id='bib0230'></span>
[[#bib0230|[46]]] T.E. Tezduyar, Y. Osawa; Finite element stabilization parameters computed from element matrices and vectors; Comp. Meth. Appl. Mech. Eng., 190 (3-4) (Octubre 2000), pp. 411–430</li>
</ol>