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 [14] u otras estructuras civiles [43] , la hidrodinámica naval [9]  and [42] y la agitación de líquidos contenidos en tanques [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 [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 [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 [37] . En el primer grupo, los métodos más difundidos hoy día son el de Volume of Fluid (VOF) [17]  and [35] y el de pseudoconcentración, o de Level Set (LS) [27]  and [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) [15]  and [26] , los métodos de partículas, combinaciones de elementos finitos con partículas [20] , y la estrategia denominada lagrangiana euleriana arbitraria (ALE, por Arbitrary Lagrangian-Eulerian ) [19]  and [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 [30] , que usa herramientas de Portable Extensible Toolkit for Scientific Computations (PETSc) [2] y de Message Passing Interface (MPI) [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 Ω_X 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 Ω_x se compone de puntos correspondientes a las coordenadas espaciales x . Por último, la configuración de referencia Ω_χ , 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:

${\boldsymbol {v}}\left({\boldsymbol {X}},t\right)={\frac {\partial {\boldsymbol {x}}}{\partial t}}_{\boldsymbol {X}};$

( 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 (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 [10] o en Folch Durán [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 [10]  and [12] :

${\overset {\mbox{ˆ}}{\boldsymbol {v}}}\left({\boldsymbol {\chi }},t\right)=$${\frac {\partial {\boldsymbol {x}}}{\partial t}}_{\boldsymbol {\chi }}.$

( 2)

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

${\boldsymbol {w}}\left({\boldsymbol {X}},t\right)={\frac {\partial {\boldsymbol {\chi }}}{\partial t}}_{\boldsymbol {X}}.$

( 3)

Por último, la velocidad convectiva c se define de la siguiente manera [10]  and [12] :

${\boldsymbol {c}}:={\boldsymbol {v}}-{\overset {\mbox{ˆ}}{\boldsymbol {v}}}=$${\frac {\partial {\boldsymbol {x}}}{\partial {\boldsymbol {\chi }}}}{\boldsymbol {w}};$

( 4)

y corresponde a la velocidad relativa entre las configuraciones material y de referencia. Recordando que ${\textstyle {\boldsymbol {w}}}$ es la velocidad relativa de la partícula con respecto a la configuración de referencia, la ecuación (4) implica que la velocidad convectiva coincide con ${\textstyle {\boldsymbol {w}}}$ , es decir ${\textstyle {\boldsymbol {c}}=}$${\boldsymbol {w}}$ , 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»[10] , que para una cantidad física escalar f es:

${\frac {\partial f}{\partial t}}_{\boldsymbol {X}}={\frac {\partial f}{\partial t}}_{\boldsymbol {\chi }}+$${\boldsymbol {c}}\cdot \nabla f.$

( 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:

$\rho \left({\partial }_{t}{\boldsymbol {v}}+{\boldsymbol {v}}\cdot \nabla {\boldsymbol {v}}-\right.$$\left.{\boldsymbol {f}}\right)-\nabla \cdot {\boldsymbol {\sigma }}=$$0;$

( 6)

sobre el dominio del flujo Ω para el tiempo t  ∈ [0, T   ], en la cual ${\textstyle {\boldsymbol {v}}}$ 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, ∂_t indica derivación parcial con respecto al tiempo y ∇ = ∇ _x es el operador de derivación con respecto a las coordenadas espaciales x del marco euleriano. El tensor ${\textstyle {\boldsymbol {\sigma }}=}$${\boldsymbol {\sigma }}({\boldsymbol {v}},p)$ 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:

${\boldsymbol {\sigma }}=-p{\boldsymbol {I}}+{\boldsymbol {T}};$

( 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 ${\textstyle {\boldsymbol {\epsilon }}=}$${\boldsymbol {\epsilon }}({\boldsymbol {v}})$ como sigue,

${\boldsymbol {T}}=2\mu {\boldsymbol {\epsilon }};$

( 9)

en la cual

${\boldsymbol {\epsilon }}={\frac {1}{2}}\left[\nabla {\boldsymbol {v}}+\right.$$\left.{\left(\nabla {\boldsymbol {v}}\right)}^{T}\right];$

( 10)

con μ la viscosidad dinámica del fluido y donde (. . .)^T 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 Ω_x , y las que relacionan estas con un dominio de referencia arbitrario Ω_χ . Para ello, se utilizan las formas ALE de las ecuaciones de conservación de cantidad de movimiento y de continuidad [10] , dadas por:

$\rho \left({\partial }_{t}{\boldsymbol {v}}+{\boldsymbol {c}}\cdot \nabla {\boldsymbol {v}}-\right.$$\left.{\boldsymbol {f}}\right)-\nabla \cdot {\boldsymbol {\sigma }}=$$0;$

( 11)

Nótese que solo la ecuación de cantidad de movimiento se ve afectada por la velocidad del dominio de referencia ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {v}}}}$ , que interviene a través de la velocidad convectiva c dada por la ecuación (4) .

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

${\boldsymbol {v}}={\boldsymbol {v}}_{D}\quad en\quad {\Gamma }_{D};$

( 13)

donde las condiciones Dirichlet están dadas sobre Γ_D , que corresponde a los contornos sólidos cuya velocidad ${\textstyle {\boldsymbol {v}}_{D}}$ es conocida, y Γ_t 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 Γ  = Γ_D  ∪ Γ_t y que Γ_D ∩ Γ_t  = ∅.

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, Γ_SL  ⊆ Γ_t , 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, σ_l y σ_g , respectivamente,

${\boldsymbol {\sigma }}_{l}\cdot {\boldsymbol {n}}={\boldsymbol {\sigma }}_{g}\cdot {\boldsymbol {n}}\quad \quad en\quad {\Gamma }_{SL}.$

( 15)

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

${\boldsymbol {T}}\cdot {\boldsymbol {n}}={\boldsymbol {0}}\quad en\quad {\Gamma }_{SL};$

( 16)

en las cuales P_atm 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 Γ_SL valen t  = − P_atm  n , o bien:

${\boldsymbol {\sigma }}_{g}\cdot {\boldsymbol {n}}=-P_{atm}\quad {\boldsymbol {n}}\quad sobre\quad {\Gamma }_{SL}.$

( 18)

En el caso de adoptarse P_atm  = 0, el término de contorno sobre Γ_SL 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 (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:

${\boldsymbol {v}}({\boldsymbol {x}},0)={\boldsymbol {v}}_{0}({\boldsymbol {x}});$

( 19)

esto es, ${\textstyle {\boldsymbol {v}}_{0}}$ es el campo de velocidades inicial definido sobre todo el dominio en tiempo t  = 0.

La discretización se realiza sobre las ecuaciones  and [38] . Para ello, el dominio Ω se divide en n_el elementos finitos Ω^e , siendo ${\textstyle E}$ el conjunto de todos los elementos y H^1h el espacio de dimensión finita definido como sigue,

$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 ,$

( 20)

donde P¹ es el espacio de polinomios de primer orden. Los espacios de funciones de peso e interpolación son

$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 ;$

( 21)

siendo n_d el número de dimensiones espaciales. El método estabilizado con streamline upwind/Petrov-Galerkin (SUPG) [7] y pressure stabilizing/Petrov-Galerkin (PSPG) [44]  and [45] se escribe de la siguiente forma: hallar  ${\textstyle {\boldsymbol {v}}^{h}\in S_{v}^{h}}$y  ${\textstyle p^{h}\in S_{p}^{h}}$tal que:

${\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$ $+{\int }_{\Omega }{\boldsymbol {\epsilon }}({\boldsymbol {w}}^{h}):{\boldsymbol {\sigma }}^{h}d\Omega$ $+\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}$ $+\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}$ $+{\int }_{\Omega }q^{h}\quad \nabla \cdot {\boldsymbol {v}}\quad d\Omega$ $+\sum _{e=1}^{n_{el}}{\int }_{{\Omega }^{e}}{\boldsymbol {\kappa }}^{h}\quad \nabla \cdot {\boldsymbol {v}}^{h}d{\Omega }^{e}$ $={\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};$

( 24)

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

${\boldsymbol {\delta }}^{h}=\tau _{SUPG}\quad ({\boldsymbol {v}}^{h}\cdot \nabla ){\boldsymbol {w}}^{h};$

( 25)

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

$\tau _{SUPG}=\left({\frac {1}{\tau _{S1}^{2}}}+{\frac {1}{\tau _{S2}^{2}}}+\right.$$\left.{\frac {1}{\tau _{S3}^{2}}}\right)^{-{\frac {1}{2}}};$

( 28)

siendo

$\tau _{S1}={\frac {h}{2\parallel {\boldsymbol {v}}\parallel }};\quad \tau _{S2}=$${\frac {\Delta t}{2}};\quad \tau _{S3}={\frac {h^{2}}{4\nu }}.$

( 29)

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

$\tau _{PSPG}=\left({\frac {1}{\tau _{P1}^{2}}}+{\frac {1}{\tau _{P2}^{2}}}+\right.$$\left.{\frac {1}{\tau _{P3}^{2}}}\right)^{-{\frac {1}{2}}};$

( 30)

con

$\tau _{P1}={\frac {h^{\#}}{2\parallel {\boldsymbol {v}}\parallel }};\quad \tau _{P2}=$${\frac {\Delta t}{2}};\quad {\tau }_{P3}={\frac {(h^{\#})^{2}}{4\nu }}.$

( 31)

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

${\nu }_{LSIC}={\frac {\parallel {\boldsymbol {v}}\parallel h}{2}}\quad z(Pe);$

( 32)

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

$Pe={\frac {\parallel {\boldsymbol {v}}\parallel h}{2\nu }}.$

( 33)

Nótese que si ${\textstyle \parallel {\boldsymbol {v}}\parallel =}$$0$ , algunas de las expresiones en las ecuaciones (29) o (31) son singulares; sin embargo, las restantes no lo son, y los parámetros de estabilización de las ecuaciones (28) y (30) resultan siempre regulares.

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

$h=2{\left(\sum _{a=1}^{n_{en}}\vert {\boldsymbol {s}}\cdot \nabla w_{a}\vert \right)}^{-1};$

( 34)

donde las ${\textstyle w_{a}}$ son las funciones asociadas al nodo a , n_en 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^# para las expresiones en la ecuación (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:

$z(Pe)={\begin{array}{ll}Pe/3,&0\leq Pe<3;\\1,&3\leq Pe.\end{array}}$

( 35)

Una vez realizada la discretización espacial de la ecuación (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 θ[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 Γ_D , esto es, sobre la velocidad ${\textstyle {\boldsymbol {v}}_{D}}$ , 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

${\boldsymbol {v}}_{D}={\boldsymbol {0}},$

( 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 (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»[14] :

$({\boldsymbol {I}}-{\boldsymbol {n}}{\boldsymbol {n}})\cdot ({\boldsymbol {n}}\cdot {\boldsymbol {\sigma }})=$$-{\frac {1}{\beta }}({\boldsymbol {I}}-{\boldsymbol {n}}{\boldsymbol {n}})\cdot ({\boldsymbol {v}}_{S}-$${\boldsymbol {v}}_{D});$

( 37)

en la cual ${\textstyle {\boldsymbol {v}}_{S}}$ 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⁻⁶  m² /s. En estos casos, puede emplearse la condición de deslizamiento perfecto ${\textstyle {\boldsymbol {v}}_{D}\cdot {\boldsymbol {n}}=}$$0$[18] , [39]  and [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, ${\textstyle {\boldsymbol {v}}_{D}\cdot {\boldsymbol {n}}=}$$0$ .

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 ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {s}}}_{j}}$ , fija, de manera que la nueva posición del nodo j sobre la SL en el tiempo t es [4] :

${\boldsymbol {x}}_{j}(t)={\boldsymbol {x}}_{0,j}+{\eta }_{j}(t)\quad {\overset {\mbox{ˆ}}{\boldsymbol {s}}}_{j};$

( 38)

en la cual ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {s}}}_{j}}$ da la dirección fija o «espina»(spine) , x_0,j es la posición inicial del nodo j y η_j la coordenada escalar, todos ellos indicados en la figura 1 . Usualmente, ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {s}}}_{j}}$ es normal a la SL en reposo para el problema a resolver.

Figura 1. Desplazamiento de los nodos de la SL.

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 [13] , [25] , [29]  and [41] , y puede expresarse como:

${\boldsymbol {v}}\cdot {\boldsymbol {n}}={\partial }_{t}{\boldsymbol {\eta }}\cdot {\boldsymbol {n}}\quad sobre\quad {\Gamma }_{SL};$

( 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 ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {s}}}}$ , es decir ${\textstyle {\boldsymbol {\eta }}=}$$\eta \quad {\overset {\mbox{ˆ}}{\boldsymbol {s}}}$ , con η representando el desplazamiento escalar, esquematizado en la figura 1 . Remplazando en la ecuación (39) ,

${\partial }_{t}\eta ={\frac {{\boldsymbol {v}}\cdot {\boldsymbol {n}}}{{\overset {\mbox{ˆ}}{\boldsymbol {s}}}\cdot {\boldsymbol {n}}}}.$

( 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 (38) y que el campo escalar en 3 dimensiones es η  = η (x₁ , x₂ , t ), la SL puede expresarse de manera implícita como:

$F({\boldsymbol {x}},t)=x_{3}-\eta =0;$

( 41)

asumiendo que la dirección de las espinas se adopta vertical, i.e. ${\textstyle {\overset {\mbox{ˆ}}{\boldsymbol {s}}}\parallel x_{3}}$ , la normal a la SL puede calcularse como

${\boldsymbol {n}}=\nabla F({\boldsymbol {x}})=\left[-{\frac {\partial \eta }{\partial x_{1}}};\quad -\right.\left.{\frac {\partial \eta }{\partial x_{2}}};\quad +\right.\left.1\right]^{T}.$

( 42)

Remplazando en la ecuación (40) con (42) , e introduciendo el parámetro de proyección

$H={\overset {\mbox{ˆ}}{\boldsymbol {s}}}\cdot {\boldsymbol {n}};$

( 43)

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

${\partial }_{t}\eta ={\frac {1}{H}}{\boldsymbol {v}}\cdot {\boldsymbol {n}}=$${\frac {1}{H}}\left[-v_{1}{\frac {\partial \eta }{\partial x_{1}}}-\right.$$\left.v_{2}{\frac {\partial \eta }{\partial x_{2}}}+v_{3}\right];$

( 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 (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,

${\partial }_{t}\eta +{\boldsymbol {v}}_{\parallel }\cdot {\overset {\mbox{ˆ}}{\nabla }}\eta =$$s\quad en\quad {\Omega }_{SL}^{'}$

( 45)

en la que la velocidad normal a las espinas es:

${\boldsymbol {v}}_{\parallel }={\frac {1}{H}}\quad {\left[v_{1};\quad v_{2}\right]}^{T};$

( 47)

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

${\overset {\mbox{ˆ}}{\nabla }}\eta ={\left[{\frac {\partial \eta }{\partial x_{1}}};\quad {\frac {\partial \eta }{\partial x_{2}}}\right]}^{T};$

( 48)

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

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

La formulación explícita del desplazamiento de la SL dada por la ecuación (45) , al igual que la dada en la ecuación (40) , resulta numéricamente inestable para ondas de superficie, dado su carácter advectivo [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) [7] aplicado a la SL por Soulaïmani et al. [39] y Güler et al. [14] , y el denominado Galerkin/Least-Squares (GLS), elegido por Behr y Abraham [5] .

El método de estabilización SUPG, introducido por Brooks y Hughes [7] , es empleado en este trabajo para la resolución del problema de advección de la ecuación (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 [10] , y se basa en los espacios funcionales ${\textstyle S_{t}:=}$$\lbrace \eta \vert \eta (\cdot ,t)\in H^{1}({\Omega }_{SL}^{'}),t\in [0,T]\vert \eta =$${\eta }_{D}\quad sobre\quad {\Gamma }_{D}^{'}\rbrace$ y ${\textstyle W:=}$$\lbrace w\in H^{1}({\Omega }_{SL}^{'})\vert w=0\quad sobre\quad {\Gamma }_{D}^{'}\rbrace$ . El problema se propone de la siguiente forma: para cualquier tiempo t  ∈ [0, T ], hallar  ${\textstyle \eta \in S}$tal que

${\int }_{{\Omega }_{SL}^{'}}w\quad {\frac {\partial \eta }{\partial t}}d{\Omega }_{SL}^{'}+$${\int }_{{\Omega }_{SL}^{'}}w\left({\boldsymbol {v}}_{\parallel }\cdot {\overset {\mbox{ˆ}}{\nabla }}\eta \right)d{\Omega }_{SL}^{'}=$${\int }_{{\Omega }_{SL}^{'}}w\quad s\quad d{\Omega }_{SL}^{'};$

( 49)

o, expresado en forma compacta,

$a\left(w,\eta \right)+c\left({\boldsymbol {v}}_{\parallel };w,\eta \right)=$$\left(w,s\right)\quad en\quad {\Omega }_{SL}^{'};$

( 50)

con condición inicial

$\left(w,\eta (\cdot ,0)\right)=\left(w,{\eta }_{0}\right).$

( 51)

En la ecuación (50) , los términos se identifican de la siguiente manera,

$a\left(w,\eta \right)={\int }_{{\Omega }_{SL}^{'}}w\quad {\frac {\partial \eta }{\partial t}}\quad d{\Omega }_{SL}^{'};$ $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}^{'};$ $\left(w,s\right)={\int }_{{\Omega }_{SL}^{'}}w\quad s\quad d{\Omega }_{SL}^{'}.$

( 52)

Una estabilización consistente para la formulación de la ecuación (50) puede escribirse como:

$a\left(w,\eta \right)+c\left({\boldsymbol {v}}_{\parallel };w,\eta \right)+$$\sum _{e}{\int }_{{\Omega }_{e}^{'}}P(w)\quad {\tau }_{S}\quad R(\eta )\quad d{\Omega }_{e}^{'}=$$\left(w,s\right)\quad en\quad {\Omega }_{SL}^{'};$

( 53)

con el tiempo intrínseco o parámetro de estabilización τ_S , que se describe más adelante en esta sección, e indicando con ${\textstyle {\Omega }_{e}^{'}}$ que la integración del término se realiza sobre cada uno de los elementos. El residuo de la ecuación diferencial es:

$R(\eta )={\frac {\partial \eta }{\partial t}}+{\boldsymbol {v}}_{\parallel }\cdot {\overset {\mbox{ˆ}}{\nabla }}\eta -$$s=L(\eta )-s;$

( 54)

y el operador ${\textstyle P(w)}$ es definido para una estabilización con SUPG en la forma:

$P(w)={\boldsymbol {v}}_{\parallel }\cdot {\overset {\mbox{ˆ}}{\nabla }}w.$

( 55)

El problema discretizado toma finalmente la forma de: hallar  ${\textstyle {\eta }^{h}\in S^{h}}$tal que

$a\left(w^{h},{\eta }^{h}\right)+c\left({\boldsymbol {v}}_{\parallel };w^{h},{\eta }^{h}\right)+$$\sum _{e}{\int }_{{\Omega }_{e}^{'}}P(w^{h})\quad {\tau }_{S}\quad R({\eta }^{h})\quad d{\Omega }_{e}^{'}=$$\left(w^{h},s\right);$

( 56)

para todo  ${\textstyle w^{h}\in W^{h}}$ , con ${\textstyle {\eta }^{h}\in S^{h}}$ , siendo ${\textstyle W^{h}}$ y ${\textstyle S^{h}}$ subespacios de dimensión finita de ${\textstyle W}$ y ${\textstyle S}$ , respectivamente.

El tiempo intrínseco τ_S es tomado como:

${\tau }_{S}={\frac {h}{2\parallel {\boldsymbol {a}}\parallel }};$

( 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 ${\textstyle {\boldsymbol {v}}_{\parallel }}$ 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 τ_S  = 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 [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[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 [3] , [5] , [14] , [21]  and [32] o la minimización de un indicador de la distorsión de los elementos [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 Ω₀ , 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,

${\sigma }_{ij,j}=0;$ ${\sigma }_{ij}=2{\tilde {\mu }}{\epsilon }_{ij}+{\tilde {\lambda }}{\delta }_{ij}{\epsilon }_{kk};$ ${\epsilon }_{ij}={\frac {1}{2}}(u_{i,j}+u_{j,i});$

( 58)

donde ${\textstyle {\tilde {\mu }}}$ y ${\textstyle {\tilde {\lambda }}}$ son las constantes elásticas de Lamé adoptadas arbitrariamente para el material, δ_ij es el tensor de Kronecker y los desplazamientos nodales

${\boldsymbol {u}}_{j}={\boldsymbol {x}}_{j}^{n+1}-{\boldsymbol {x}}_{j}^{0};$

( 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 ${\textstyle {\tilde {\nu }}}$ y del módulo de elasticidad longitudinal ${\textstyle {\tilde {E}}}$ , 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 ${\textstyle {\tilde {\nu }}=}$$0,3$ teniendo en cuenta que para ${\textstyle {\tilde {\nu }}\rightarrow 0,5}$ , 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 [32] y Battaglia et al. [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 figura 3 cuya SL tiene una posición inicial dada por

$h\left(x\right)=1,5{\mbox{ }}m{\mbox{ }}+{\mbox{ }}a_{0}sin\left[\pi \left(1/2-\right.\right.$$\left.\left.x\right)\right]$

( 60)

donde a₀ 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.

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

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

$a(t)={\frac {4{\nu }^{2}k^{4}}{8{\nu }^{2}k^{4}+{{\omega }_{0}}^{2}}}\quad a_{0}\quad erfc(\nu k^{2}t)^{1/2}+$$\sum _{i=1}^{4}{\frac {z_{i}}{Z_{i}}}\left({\frac {{\omega }_{0}^{2}a_{0}}{z_{i}^{2}-\nu k^{2}}}\right)exp[(z_{i}^{2}-$$\nu k^{2})t]\quad erfc(z_{i}t^{1/2});$

( 61)

en la cual ν es la viscosidad cinemática del fluido, k   es el número de onda, ${\textstyle {\omega }_{0}^{2}=}$$gk$ es la frecuencia natural para el caso no viscoso, y cada z_i es una raíz de la siguiente ecuación algebraica:

$z^{4}+k^{2}\nu z^{2}+4(k^{2}\nu )^{3/2}z+{\nu }^{2}k^{4}+{\omega }_{0}^{2}=$$0;$

( 62)

donde Z₁  = (z₂  − z₁ )(z₃  − z₁ )(z₄  − z₁ ), Z₂ , Z₃ , Z₄ se obtienen por permutación circular de los índices y erfc (…) es la función error para variable compleja. Notar que la ecuación (61) transcripta en [32] a partir de [31] tiene un error de tipeo.

El ejemplo fue resuelto para una amplitud inicial de la deformación a₀  = 0, 01 m, aceleración gravitatoria unitaria g  = − 1, 0 m/s² y una malla de 40 × 60 elementos cuadrangulares para un tanque representado en la figura 3 de h  = 1, 5 m de altura y d  = 1, 0 m de ancho, tal que k  = 2d , con un paso de tiempo Δt  = 2, 12 × 10⁻²  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 [6] , [32]  and [33] , en el cual la viscosidad cinemática y la densidad son ν  = 10⁻³  m² /s y ρ  = 1 kg/m³ , respectivamente.

La 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₀ . Las condiciones de contorno elegidas son de deslizamiento perfecto sobre los contornos Γ_D , 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 [31] . Sobre la SL, P_atm  = 0 y las tensiones tangenciales son T  · n  = 0 .

Figura 4. Curva de solución analítica y resultados numéricos para el problema de agitación con viscosidad cinemática ν  = 10−3  m2 /s .

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 (50) no influirían sensiblemente en los resultados, pues las velocidades ${\textstyle {\boldsymbol {v}}_{\parallel }}$ son menores que las velocidades que intervienen en el término fuente s , en un factor del orden de 10⁻³ .

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⁻⁴  m² /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₀ , en tanto que en la resolución con ν  = 10⁻²  m² /s , representada en la figura 5 , el error relativo a la amplitud inicial a₀ es del 2%.

Figura 5. Curva de solución analítica y resultados numéricos para el problema de agitación con viscosidad cinemática ν  = 10−2  m2 /s .

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⁻⁸ 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_g  = 0,60 m de eje inclinado en un ángulo γ  = 10^o con respecto a la aceleración gravitatoria g , como puede apreciarse en la figura 6 . El movimiento del fluido se genera a partir de liberar la SL desde una posición inicial diferente de la hidrostática.

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

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 [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 ${\textstyle {\kappa }_{s}^{n}\approx 1.746}$ para γ  = 10^o , modo de orden n  = 1 y número de autovalor s  = 1. Dado el autovalor, la frecuencia de oscilación σ se calcula como:

$\sigma ={\left({\frac {g\quad {\kappa }_{1}^{\left(1\right)}}{R}}\right)}^{0,5}\approx 5,853\quad s^{-1};$

( 63)

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

$T={\frac {2\pi }{\sigma }}\approx 1,0735\quad s;$

( 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 figura 7 . La discretización fue propuesta tal que el eje del cilindro inclinado coincide con la dirección espacial x₃ , en tanto que la aceleración gravitatoria se impuso con componentes g₁  = g  sin (γ ), g₂  = 0 y g₃  = g  cos (γ ) en las direcciones x₁ , x₂ y x₃ , respectivamente. El fluido contenido en el recipiente es glicerina, cuyas propiedades físicas son: viscosidad cinemática ν  = 6, 5⁻⁴  m² /s y densidad ρ  = 1.260 kg/m³ . El paso de tiempo adoptado para el análisis es Δt  = 0, 025 s, hasta un tiempo final de 5 s.

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

Las condiciones de borde para el problema de NS sobre la SL vienen dadas por un valor de presión P_atm  = 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 ${\textstyle {\boldsymbol {v}}=}$${\boldsymbol {0}}$ , 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₀  = 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₁ de cada nodo de la superficie libre, normalizado con el radio R , representada en la 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 (64) .

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

La figura 9 muestra los desplazamientos de los nodos extremos en dirección x₁ sobre la SL, en la cual se indica con 1 el nodo en x₁ positivo, es decir, donde la deformación de la SL es mayor medida en dirección x₃ , 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.

Figura 9. Desplazamientos de los nodos extremos en x1 sobre la SL del tanque cilíndrico de eje inclinado.

El problema fue resuelto numéricamente en 10 nodos del cluster Aquiles [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⁻⁸ . 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 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⁻⁶  m² /s, tal como en la bibliografía consultada [18] , [39]  and [40] . La aceleración que actúa sobre el dominio es G  = [g₁  ; g₂ ]^T , de componente vertical de magnitud g₂  = − g  = − 9, 81 m/s² , y componente horizontal dependiente del tiempo, propuesta de variación sinusoidal. El valor de densidad adoptado es de ρ  = 1.000 kg/m³ .

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

Este ejemplo fue tomado de Huerta y Liu [18] , y es reproducido por Soulaïmani et al. [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₁  = 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 λ  = 2W , por lo tanto,

$f={\sqrt {{\frac {g}{4\pi W}}\quad tanh{\frac {\pi D}{W}}}}\approx 0,89825\quad 1/s.$

( 65)

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

La figura 11 muestra el desplazamiento del nodo superior izquierdo comparado con las envolventes de los resultados de la bibliografía [18]  and [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.

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

En la figura 11 puede apreciarse que los desplazamientos en subida son mayores que las amplitudes en bajada, en coincidencia con resultados de modelos físicos [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 figura 12 .

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

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 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 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 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^1,67 .

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

	Cantidad elementos		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

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 .

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, 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⁻⁶  m² /s y densidad ρ  = 1.000 kg/m³ . 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 [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 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, figura 14 (centro), y superior derecho, 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.

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

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 [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 [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

1. [1] Cluster Aquiles en CIMEC. http://www.cimec.org.ar/aquiles , 2011.
2. [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)
3. [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
4. [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
5. [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
6. [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
7. [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
8. [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.
9. [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
10. [10] J. Donea, A. Huerta; Finite Element Methods for Flow Problems; John Wiley & Sons (2003)
11. [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
12. [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.
13. [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.
14. [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
15. [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
16. [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
17. [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
18. [18] A. Huerta, W.K. Liu; Viscous flow with large free surface motion; Comp. Meth. Appl. Mech. Eng., 69 (1988), pp. 277–324
19. [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
20. [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
21. [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
22. [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
23. [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
24. [24] MPI. Message Passing Interface. http://www.mpi-forum.org , 2008.
25. [25] A.D. Myshkis, V.G. Babskii, N.D. Kopachevskii, L.A. Slobozhanin, A.D. Tyuptsov; Low-Gravity Fluid Mechanics; Springer-Verlag (1987)
26. [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
27. [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
28. [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
29. [29] J. Pedlosky; Waves in the ocean and atmosphere. Introduction to wave dynamics; Springer (2003)
30. [30] PETSc-FEM. A general purpose, parallel, multi-physics FEM program, 2011. GNU General Public License (GPL), http://www.cimec.org.ar/petscfem .
31. [31] A. Prosperetti; Motion of two superposed viscous fluids; Phys. Fluids, 24 (7) (1981), pp. 1217–1223
32. [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
33. [33] B. Ramaswamy; Numerical simulation of unsteady viscous free surface flow; J. Comput. Phys., 90 (2) (1990), pp. 396–430
34. [34] H. Saito, L.E. Scriven; Study of coating flow by the finite element method; J. Comput. Phys., 42 (1) (1981), pp. 53–76
35. [35] R. Scardovelli, S. Zaleski; Direct numerical simulation of free-surface and interfacial flow; Annu. Rev. Fluid Mech., 31 (1999), pp. 567–603
36. [36] J.A. Sethian; A Fast Marching Level Set Method for Monotonically Advancing Fronts; National Academy of Sciences (1995)
37. [37] W. Shyy, H.S. Udaykumar, M.M. Rao, R.W. Smith; Computational Fluid Dynamics with Moving Boundaries; Taylor and Francis (1996)
38. [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
39. [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
40. [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
41. [41] J.J. Stoker; Water waves. The mathematical theory with applications, volumen IV of Pure and Applied Mathematics; Interscience (1957)
42. [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
43. [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
44. [44] T.E. Tezduyar; Stabilized finite element formulations for incompressible flow computations; Adv. Appl. Mech., 28 (1991), pp. 1–44
45. [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
46. [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