Elementos finitos mixtos estabilizados para flujos confinados de Bingham y de Herschel-Bulkley. ParteI: Formulación

Resumen

En este trabajo se presenta una metodología para la resolución de las ecuaciones de Navier-Stokes para los fluidos viscoplásticos de Bingham y de Herschel-Bulkley mediante el método de los elementos finitos mixtos estabilizados velocidad/presión. Se desarrolla una formulación teórica y se realiza la implementación computacional.

Los fluidos viscoplásticos se caracterizan por presentar una tensión de corte mínima, denominada tensión de fluencia. Por encima de esta tensión de corte mínima el fluido comienza a moverse. En caso de no superarse esta tensión de fluencia, el fluido se comporta como un cuerpo rígido o cuasirrígido, con velocidad de deformación nula.

Se presentan inicialmente las ecuaciones de Navier-Stokes para un fluido incompresible. Se incluye una revisión de los modelos reológicos viscoplásticos. Se hace una descripción detallada de los mismos. Se describen los modelos viscoplásticos regularizados de Papanastasiou. Se proponen modelos regularizados de doble viscosidad como alternativa a los comúnmente usados.

Se formula el modelo discreto, así como la formulación estabilizada con los métodos de subescalas algebraica (Algebraic SubGrid Scale [ASGS]), de subescalas ortogonales (Orthogonal Subgrid Scale [OSS]) y de subescalas ortogonales desacopladas (split-OSS).

La metodología descrita en este trabajo proporciona la base para el desarrollo de una herramienta computacional para estudiar flujos viscoplásticos confinados, muy comunes en la industria.

Abstract

This work presents a methodology for the solution of the Navier-Stokes equations for Bingham and Herschel-Bulkley viscoplastic fluids using stabilized mixed velocity/pressure finite elements. The theoretical formulation is developed and implemented in a computer code.

Viscoplastic fluids are characterized by a minimum shear stress called yield stress. Above this yield stress, the fluid is able to flow. Below this yield stress, the fluid behaves as a quasi-rigid body, with zero strain-rate.

First, the Navier-Stokes equations for incompressible fluid are presented. A review of the viscoplastic rheological models is included, with a detailed description of these models. The regularized viscoplastic models due to Papanastasiou are described. Double viscosity regularized models are proposed as an alternative to the models commonly used.

The discrete model is developed, and the Algebraic SubGrid Scale (ASGS) stabilization method, the Orthogonal Subgrid Scale (OSS) method and the split orthogonal subscales method are introduced.

The methodology proposed in this work provides a computational tool to study confined viscoplastic flows, common in industry.

1

Introducción

En el presente trabajo se presenta la formulación continua y su correspondiente versión discreta para modelos de elementos finitos mixtos velocidad/presión de flujos confinados viscoplásticos de Bingham y de Herschel-Bulkley.

Los fluidos viscoplásticos de Bingham y de Herschel-Bulkley son fluidos no-newtonianos que se caracterizan por presentar una tensión de corte mínima, denominada «tensión de fluencia». Por encima de esta tensión de corte mínima el fluido comienza a moverse. En caso de no superar esta tensión de fluencia, el fluido se comporta como un cuerpo rígido o cuasirrígido, con velocidad de deformación nula.

En la industria, los fluidos de Bingham pueden modelar el comportamiento de las pinturas, de los plásticos, de productos alimenticios como la mayonesa y el kétchup, entre otros. Los fluidos de Herschel-Bulkley incluyen, por ejemplo, el comportamiento de las pastas, algunos geles y los fluidos de perforación. En el medio ambiente, estos fluidos pueden modelar flujos de detritos, entre otros.

El movimiento de los fluidos isotérmicos se describe mediante las ecuaciones de conservación de masa y momentum, representado por las ecuaciones de Navier-Stokes. Numerosos ensayos experimentales han demostrado que las ecuaciones de Navier-Stokes bajo condiciones isotérmicas describen exactamente el flujo incompresible de los fluidos. Las ecuaciones de Navier-Stokes requieren de una ecuación constitutiva para caracterizar el tipo de fluido. Esta ecuación define el valor de las tensiones en función de la dinámica del flujo y está asociada con la viscosidad del fluido (modelo reológico).

El modelo que dio inicio al estudio de los materiales viscoplásticos fue el modelo plástico de Bingham [1], formulado por Eugene C. Bingham para describir el comportamiento de las pinturas. El modelo de Herschel-Bulkley [2] se considera como un modelo generalizado de Bingham, aunque ha sido menos estudiado que el modelo de Bingham.

Ambos fluidos exhiben una fuerte discontinuidad en su comportamiento reológico debido a la existencia de la tensión de fluencia que es difícil de tratar numéricamente. Para solventar este problema, autores como Bercovier y Engelman [3], Tanner y Milthorpe [4] y Beris et al. [5], entre otros, han propuesto diferentes formulaciones regularizadas. Tanner y Milthorpe fueron los primeros que simularon el problema utilizando un modelo de doble viscosidad aplicable a ambos fluidos. Beris y sus colegas centraron sus estudios en el fluido de Bingham, utilizando el criterio de Von Mises [6] en las zonas de no fluencia y el modelo ideal de Bingham en la zona de fluencia. En 1987, Papanastasiou [7] propuso un modelo regularizado aplicable tanto en las zonas de no fluencia como en las zonas de fluencia para estos 2 fluidos. Souza Mendes y Dutra (SMD) [8] han propuesto recientemente una modificación del modelo de Papanastasiou.

En el presente trabajo se proponen nuevos modelos regularizados para el fluido de Bingham y el fluido de Herschel-Bulkley como alternativa a los modelos regularizados comúnmente usados.

En el caso de los materiales viscoplásticos, el método numérico más utilizado es el método de los elementos finitos (MEF) [7,9–11]. Para abordar el problema de flujo incompresible mediante el MEF, se emplea la formulación mixta de velocidad/presión (u/p). La formulación estándar de Galerkin presenta 2 fuentes de inestabilidades.

La primera es la presencia del término convectivo en las ecuaciones de gobierno que puede resultar en oscilaciones numéricas en el campo de la velocidad. La segunda fuente de inestabilidad es la combinación inapropiada de espacios de interpolación para los campos de velocidad y presión. Esta falta de estabilidad produce oscilaciones numéricas en el campo de las presiones. Para que el problema discreto sea estable, los espacios de interpolación usados para la velocidad y la presión deben satisfacer la condición inf-sup de compatibilidad o condición de Babuška-Brezzi [12]. La formulación de igual interpolación lineal usada en este trabajo no cumple con la condición Babuška-Brezzi [13].

En ambos casos el problema necesita estabilizarse para poder probar convergencia a la solución del problema. Los métodos de estabilización más usados en la actualidad están basados en los métodos de subescalas. Hughes fue el pionero en estos métodos de subescalas (SubGrid Scale [SGS]), proponiendo el método de estabilización de subescalas algebraicas (Algebraic SubGrid Scale stabilization method [ASGS]) [14] para una ecuación escalar de difusión-reacción. Codina [15] amplió esta aproximación algebraica aplicándola a sistemas escalares multidimensionales.

Posteriormente, Codina [15] propuso adoptar un espacio de subescalas ortogonales al espacio de los elementos finitos, fundamentando así el método de estabilización de subescalas ortogonales (Orthogonal Subscale Stabilization method [OSS]). El método OSS se ha aplicado al problema de Stokes, al problema de convección-difusión-reacción y a las ecuaciones de Navier-Stokes, entre otros [15,16]. La estabilización OSS ha sido reformulada en una nueva versión del método llamada estabilización split-OSS [17], computacionalmente más ventajosa. Actualmente se usan en problemas muy variados, tanto de mecánica de fluidos [15,17–21] como de mecánica de sólidos [19,22–31].

En la parte i de este trabajo se presenta el problema del flujo confinado con un amplio desarrollo de los modelos constitutivos para flujos viscoplásticos de Bingham y de Herschel-Bulkley, así como los modelos viscoplásticos regularizados propuestos en este trabajo. Finalmente, se presenta el modelo discreto incorporando los modelos de Bingham y de Herschel-Bulkley. Como métodos de estabilización se discuten los métodos ASGS, OSS y split-OSS.

2

Modelo continuo para el problema del flujo confinado

El problema continuo de dinámica de fluidos incompresibles e isotérmicos puede resolverse completamente considerando las ecuaciones de Navier-Stokes.

Las ecuaciones de Navier-Stokes, planteadas inicialmente para fluidos newtonianos, pueden usarse conjuntamente con los modelos reológicos viscoplásticos de Bingham y Herschel-Bulkley en la ecuación constitutiva.

Considérese Ω, un dominio abierto y acotado dimensional de ℝd, donde d = 2 o 3 es el número de dimensiones del espacio, Γ = ∂Ω es su contorno, que puede ser dividido en el contorno con condiciones de Dirichlet (velocidad impuesta) Γ_d = ∂Ω_d y el contorno con condiciones de Neumann (tracciones impuestas) Γ_n = ∂Ω_n de forma que Γ = Γ_d ∪ Γ_n, [0, T] es el intervalo de tiempo de análisis.

El problema de Navier-Stokes consiste en encontrar una velocidad u y una presión p tal que:

con

y

en Ω, t ∈ [0,t], donde σ es el tensor de tensiones, τ es el tensor de tensiones desviadoras, f es el vector de fuerzas de volumen y ρ es la densidad del fluido.

A esas ecuaciones deben añadirse condiciones iniciales de la forma u = u₀ en Ω, t₀ = 0 y condiciones de contorno:

Condiciones de Dirichlet:

Condiciones de Neumann:

donde n es el vector unitario normal al contorno ∂Ω. Por simplicidad se tomará en el contorno Γ_d la velocidad u¯=0,t∈[0,T]. El vector t es el vector de tracción sobre el contorno con condiciones de Neumann.

3

Líneas de corriente

Las líneas de corriente son curvas tangentes en cada punto al campo de velocidades. En un flujo estacionario, las líneas de corriente no varían con el tiempo, mientras que en flujo transitorio sí lo hacen.

Las líneas de corriente para un flujo bidimensional con un campo de velocidades u = (u_x,u_y) coinciden con las líneas de nivel de la función ϕ, solución de la ecuación laplaciana:

con la condición de contorno ϕ = 0.

4

Ecuación constitutiva para fluidos viscoplásticos

La ecuación constitutiva relaciona las tensiones con la presión y la velocidad de deformación. En el caso de los fluidos, esta relación se denomina también modelo reológico.

El tensor de tensiones σ se descompone en su parte volumétrica y desviadora como:

donde p es la presión, I es el tensor de identidad de segundo orden y τ es el tensor de las tensiones desviadoras.

Para un fluido newtoniano, usando la hipótesis de Stokes y usando la ecuación de incompresibilidad, el tensor desviador de tensiones se expresa como:

donde u es el vector de velocidades, μ es la viscosidad dinámica (constante en caso de fluido newtoniano) y ¿(·) es el gradiente simétrico de la velocidad:

donde ▿u es el gradiente de la velocidad y (▿u^t) es la transpuesta del mismo.

El valor de la magnitud del tensor de la velocidad de deformación, γ, se toma como la raíz del segundo invariante del tensor simétrico:

La magnitud del tensor desviador o tensión efectiva, τ, se toma como la raíz del segundo invariante del tensor de las tensiones desviadoras:

De acuerdo con lo anterior, la ecuación (7) se puede escribir como:

Dependiendo de los valores de la viscosidad en función de la velocidad de deformación, μ=μγ˙, pueden distinguirse diferentes ecuaciones constitutivas que representan diferentes modelos reológicos no-newtonianos. Un amplio número de estos materiales pueden verse en Bird et al. [32]. Se estudian en este trabajo los modelos para fluidos viscoplásticos, en particular para el modelo de Bingham y de Herschel-Bulkley.

4.1

Modelo ideal de Bingham

Eugene C. Bingham describió las pinturas con este modelo en 1919, publicado en su libro Fluidity and Plasticity[1]. El modelo fue analizado por Oldroyd [33], Reiner [34] y Prager [35]. Los plásticos de Bingham requieren de una tensión de corte mínima, τ_y, a partir de la cual comienzan a moverse (fig. 1).

Figura 1.

Curvas reológicas. Modelo bilineal.

En el modelo de Bingham la viscosidad está dada por:

donde μ₀ es la viscosidad plástica, y la viscosidad aparente μ(γ˙) disminuye con el incremento en la magnitud de la velocidad de deformación γ˙; τ es la magnitud del tensor de tensiones desviadoras.

En consecuencia, el tensor de tensiones desviadoras es:

Para definir si una partícula del fluido se mueve o no, es decir, si está en fluencia o no, se comprueba si la magnitud del tensor de tensiones desviadoras, τ, excede o no el valor de la tensión de fluencia, τ_y. Cuando la magnitud del tensor de tensiones del fluido, τ, supera la tensión de fluencia, el comportamiento es similar al de un fluido newtoniano; en caso contrario, el fluido no presenta deformaciones por corte.

4.1.1

Magnitudes adimensionales

4.1.1.1

Número de Bingham

Para estudiar flujos de Bingham se define un número adimensional denominado número de Bingham, Bn. El número de Bingham, sugerido por Bird et al. [32], se define como:

donde V es una velocidad característica del flujo viscoplástico, H es una longitud característica y μ es la viscosidad del fluido de Bingham.

El número de Bingham relaciona la tensión de fluencia, τ_y, con la tensión ocasionada por una velocidad de deformación característica γ˙0=VH.

En el caso de un fluido newtoniano el valor Bn es nulo, Bn = 0; en el límite opuesto, para fluidos en no fluencia (sólido) el número de Bingham puede tener valores muy altos, Bn → ∝ [10].

4.1.1.2

Tensión adimensional de fluencia

En flujos viscoplásticos es conveniente mostrar los resultados en función de una tensión de fluencia adimensional, τy*, definida por Papanastasiou [7] como:

donde V_N es una velocidad característica tomada como la velocidad promedio del líquido newtoniano de la misma viscosidad del modelo viscoplástico.

4.1.1.3

Número de Reynolds

El número de Reynolds, Rn, define si el régimen del flujo es laminar o turbulento. El empleado en este trabajo es el que Scott et al. [36] y Mitsoulis y Huilgol [37], entre otros, usan en trabajos previos para flujo newtoniano. Para el flujo laminar de Bingham:

donde ρ es la densidad del fluido, V_B es la velocidad promedio del flujo, H es una longitud característica y μ es la viscosidad para el fluido de Bingham.

4.2

Modelo ideal de Herschel-Bulkley

En el modelo plástico de Herschel-Bulkley [2] se combinan la tensión de fluencia y la ley potencial. En el modelo de Herschel-Bulkley la viscosidad aparente está dada por:

Al igual que para el modelo de Bingham, los materiales de Herschel-Bulkley requieren de una tensión de corte mínima, τ_y, para que el material fluya. Para niveles de tensión por encima de la tensión de fluencia, el material fluye con una relación no lineal tensión-velocidad de deformación como un fluido pseudoplástico (n < 1) o dilatante (n > 1) determinado por el exponente de la ley de potencia (n).

El tensor desviador resulta:

Si n = 1, se tiene el fluido de Bingham como caso particular [35] y el índice de consistencia es igual a la viscosidad plástica del material k = μ₀. Si la tensión de fluencia es nula, τ_y = 0, se recupera la ley potencial.

4.2.1

Magnitudes adimensionales

4.2.1.1

Número generalizado de Bingham

Para materiales que obedecen el modelo de Herschel-Bulkley se define el número generalizado de Bingham, Bn* o número de Oldroyd [33], Od, como:

donde τ_y es la tensión de fluencia, H es una longitud característica, k es el índice de consistencia y n es el índice potencial. La velocidad V es una velocidad característica.

4.2.1.2

Número de Reynolds

El número de Reynolds usado en flujos de Herschel-Bulkley viene dado por la ley potencial [4,38]:

5

Modelos viscoplásticos regularizados. Fluido de Bingham

Los modelos reológicos ideales con tensión de fluencia presentan 2 problemas:

• •

  Existe una singularidad para la viscosidad cuando la velocidad de deformación es nula.

• •

  Además, en algunos casos no está acotada la función de la viscosidad cuando la velocidad de deformación tiende a cero lim μγ˙→0→∞.

Estos problemas no constituyen una limitación en soluciones analíticas para problemas simples, pero sí constituyen un serio inconveniente de cara a la solución numérica [7,39,40].

Para evitar estas dificultades y lograr una conveniente formulación computacional, se han propuestos diferentes modelos regularizados.

Un modelo muy utilizado es el modelo de doble viscosidad, inicialmente propuesto por Tanner y Milthorpe [4]. Otro de los modelos más usados en nuestros días es el modelo de Papanastasiou [7]. Souza Mendes y Dutra (SMD) [8] han propuesto recientemente una modificación del modelo de Papanastasiou.

5.1

Modelo de Tanner y Milthorpe

El modelo introducido originalmente por Tanner y Milthorpe [4] es un modelo con doble viscosidad lineal que regulariza el fluido de Bingham.

Este modelo de doble viscosidad sustituye el comportamiento rígido del modelo ideal para valores de tensiones por debajo de la tensión de fluencia por una dependencia lineal entre la tensión y la velocidad de deformación. El modelo para el fluido de Bingham es:

donde μ_r es la viscosidad crítica y γ˙c es la velocidad de deformación crítica. Para esta velocidad de deformación crítica la tensión es:

Por tanto, la velocidad de deformación crítica es:

El valor de μ_r debe ser grande para aproximar el modelo ideal. Una buena aproximación recomendada por Beverly y Tanner es tomar 300≤μrμ0≤1000.

El inconveniente de este modelo es que para la viscosidad crítica se tiene una tensión crítica τ_c algo mayor que la tensión de fluencia, como se muestra en la figura 1. Esta tensión crítica es:

5.2

Modelo de Papanastasiou

Uno de los intentos por solventar la limitación debida a la singularidad de la viscosidad para γ˙→0 se debe a Papanastasiou [7], que propuso una regularización exponencial para el término de la tensión de fluencia del modelo de Bingham. La misma idea se ha usado posteriormente con el modelo de Herschel-Bulkley.

La ventaja que presenta el modelo es que describe con una sola ecuación tanto las zonas de fluencia como las de no fluencia, mediante una función suavizada de la viscosidad que depende de la velocidad de deformación y de un parámetro de regularización (m), modificando la viscosidad aparente μ(γ˙) del modelo ideal de la manera:

En la figura 2 se puede apreciar la influencia del parámetro de regularización m.

Figura 2.

Modelo regularizado de Papanastasiou para el fluido de Bingham con diferentes valores del parámetro de regularización, m.

La viscosidad en la ecuación (26) está acotada cuando el gradiente de la velocidad de deformación tiende a cero. Desarrollando en serie de Taylor y despreciando los términos de más de segundo orden, se tiene que:

Para valores muy altos del parámetro de regularización, m, este valor límite de la viscosidad puede causar problemas numéricos. En este caso es aconsejable definir un valor de truncamiento μ_t < μ_max para velocidades de deformación muy bajas γ˙<γ˙c.

5.3

Modelo de Souza Mendes y Dutra (SMD)

El modelo regularizado SMD de Souza Mendes y Dutra [8] es similar al modelo de Papanastasiou, pero la regularización exponencial afecta a todos los términos de la viscosidad:

Además, el parámetro de regularización m se sustituye por un parámetro reológico que depende de la viscosidad del fluido a cero cizallamiento η_o y la tensión de fluencia τ_y, m = η₀/τ_y.

El límite para la viscosidad cuando la velocidad de deformación es cero es:

6

Modelos viscoplásticos regularizados. Fluido de Herschel-Bulkley

Para el fluido de Herschel-Bulkley se proponen modelos regularizados análogos a los del fluido de Bingham.

6.1

Modelo de Tanner y Milthorpe

El modelo de Tanner y Milthorpe [4] para el fluido de Herschel-Bulkley es:

donde μ_r es la viscosidad crítica y γ˙c es la velocidad de deformación crítica. La velocidad de deformación crítica se obtiene resolviendo la siguiente ecuación no lineal implícita:

6.2

Modelo de Papanastasiou

La regularización propuesta por Papanastasiou (fig. 3) es también aplicable al modelo de Herschel-Bulkley. La viscosidad aparente queda definida como:

Figura 3.

Modelo regularizado de Papanastasiou para el fluido de Herschel-Bulkley con diferentes valores del parámetro de regularización m.

La influencia del parámetro m en el fluido de Herschel-Bulkley-Papanastasiou puede verse en la figura 3.

Al igual que para el modelo de Bingham, el modelo de Herschel-Bulkley requiere en la implementación numérica un valor de truncamiento μ_t.

El valor límite de la viscosidad cuando la velocidad de deformación tiende a cero varía de acuerdo con el valor n. El límite para cada término de la ecuación (32) y el límite resultante para la viscosidad se muestran en la tabla 1. Se observa que para fluidos pseudoplásticos (n < 1) la viscosidad no está acotada. En estos casos es imprescindible la aplicación del procedimiento de truncamiento.

Tabla 1.

Modelo regularizado de Papanastasiou. Valores límites para la viscosidad cuando la velocidad de deformación tiende a cero

Exponente	Términos de la viscosidad		Viscosidad
n	limγ˙→0kγ˙n−1	limγ˙→0τyγ˙[1−exp(−mγ˙)]	limy˙→0μ(γ˙)
Si n > 1	0	mτy	mτy
Si n = 1	k = μ	mτy	μ + mτy
Si n < 1	∞	mτy	∞

6.3

Modelo Souza-Mendez-Dutra

El modelo SMD (Souza-Mendez-Dutra) regulariza el modelo ideal de Herschel-Bulkley en la forma:

con m = η₀/τ_y

El valor de la viscosidad cuando la velocidad de deformación tiende a cero se muestra en la tabla 2. La viscosidad está acotada para cualquier valor de n, lo cual es una ventaja en la implementación numérica.

Tabla 2.

Modelo regularizado SMD. Valores límites para la viscosidad cuando la velocidad de deformación tiende a cero

Exponente	Términos de la viscosidad		Viscosidad
n	limγ˙→0kγ˙n−1[1−exp(−mγ˙)]	limγ˙→0τyγ˙[1−exp(−mγ˙)]	limγ˙→0μ(γ˙)]
Si n > 1	0	mτy	mτy
Si n = 1	0	mτy	mτy
Si n < 1	0	mτy	mτy

Estos modelos presentados han formado parte de los estudios de diferentes problemas resueltos por Papanastasiou [7], Kelessidis et al. [41], Westerberg et al. [42] y Dall’Onder dos Santos et al. [8], entre otros.

7

Modelos viscoplásticos regularizados propuestos

Se proponen a continuación sendos modelos viscoplásticos regularizados para los fluidos de Bingham y de Herschel-Bulkley.

Ambos son modelos de doble viscosidad, basados en los modelos descritos anteriormente.

Presentan, por tanto, viscosidad constante, pero de diferente magnitud en las zonas de fluencia (por encima de la velocidad de deformación crítica) y en las de no fluencia (por debajo de la velocidad de deformación crítica).

7.1

Modelo regularizado de Bingham de doble viscosidad

Este modelo es idéntico al modelo bilineal, pero la viscosidad crítica, μ_r, se toma igual al valor límite regularizado de Papanastasiou y SMD correspondiente a y˙=0, esto es, μ_r = mτ_y, en función del parámetro de regularización m (fig. 4). Por tanto:

Figura 4.

Modelo ideal de Bingham y modelo regularizado Bingham-DV.

En este caso, el valor de la velocidad de deformación crítica es:

En la figura 5 se compara el modelo bilineal propuesto con los modelos de Papanastasiou y SMD para m = 10 s, 100 s con tensión de fluencia τ_y = 10 Pa y μ₀ = 0.2 Pas.

Figura 5.

Comparación entre el modelo bilineal propuesto para el fluido de Bingham (DV) con el modelo regularizado de Papanastasiou y el modelo SMD para m = 10, 100 s.

7.2

Modelo regularizado de Herschel-Bulkley de doble viscosidad

De forma análoga, se propone un modelo regularizado de doble viscosidad para el fluido de Herschel-Bulkley, en el que la viscosidad crítica se toma como el valor límite de la viscosidad del modelo de Papanastasiou y SMD, μ_r = mτ_y, en función del parámetro de regularización m (fig. 6). Esto es:

Figura 6.

Modelo ideal de Herschel-Bulkley y modelos regularizados de Herschel-Bulkley-Papanastasiou y Hershel-Bulckley-DV, n>1.

El valor crítico de la velocidad de deformación cuando n ≠ 1 ha de determinarse de forma iterativa imponiendo continuidad en las tensiones de corte (para γ˙=γ˙c):

8

Modelo discreto. Formulación de elementos finitos

En esta sección se describe el modelo discreto de elementos finitos para las ecuaciones de Navier-Stokes, correspondiente al modelo continuo que se describe en la sección anterior.

La estrategia de discretización adoptada en este trabajo consiste en 2 pasos. El primero es discretizar las ecuaciones en el tiempo usando un esquema de integración de diferencias finitas, y luego, la aproximación de elementos finitos se desarrolla en el espacio. Este procedimiento desacopla errores que vienen de la discretización temporal y espacial. Es de destacar que la estrategia más común es al revés: primero se discretiza en el espacio y luego en el tiempo.

El segundo y tercer término de la ecuación (1) es el término convectivo y el de la ecuación constitutiva para fluidos no-newtonianos, en particular para los modelos viscoplástico de Bingham y de Herschel-Bulkley, con μ=μ(γ˙). Ambos términos son no lineales y, por tanto, es necesaria la linealización del mismo. En este trabajo se linealizan el término convectivo y el término viscoso con el método de Picard por su robustez.

La presencia de una derivada temporal en la ecuación (1) precisa de un algoritmo de integración en el tiempo. La discretización temporal puede hacerse por diferencias finitas usando la regla trapezoidal generalizada. Este es el método de diferencias finitas más simple y de un solo paso. Incluye como caso particular el método de diferenciación hacia atrás de Euler Backward Differentiation Formula (BDF1), entre otras posibilidades.

Para la discretización espacial mediante el método de elementos finitos es necesario construir subespacios discretos V_h ⊂ V, y Q_h ⊂ Q que aproximen los espacios continuos.

Sean V_h y Q_h los espacios de elementos finitos para interpolar las funciones vectoriales (velocidad) y escalares (presión), respectivamente, y sea Ω una partición de elementos finitos Ω = ∪ Ω^e, e = 1, …, n_ele, donde n_ele es el número de elementos.

En la formulación estándar de Galerkin se toman las funciones de test iguales a las funciones de forma, así que v_h ∈ V_h y q_h ∈ Q_h.

El cálculo con elementos finitos de flujos incompresibles con la formulación estándar de Galerkin deben estabilizarse debido a que la formulación de igual interpolación lineal tanto para la velocidad como para la presión usada en este trabajo no cumple con la condición Babuška-Brezzi. La referencia [13] ofrece una descripción comparativa de varios de los métodos de estabilización propuestos en las últimas décadas.

Los métodos de estabilización más usados en la actualidad están basados en los métodos de subescalas [43,44]. Estos métodos consisten en descomponer la solución, por ejemplo, la velocidad u en 2 componentes u=uh+u¯; una componente u_h, resuelta en la escala de la malla de elementos finitos considerada y una subescala u¯, que no puede ser capturada por la partición de elementos finitos y que se resuelve analíticamente. La aproximación particular usada para la escala submalla define el modelo numérico.

La solución u¯ en la escala fina se obtiene a partir del residuo de la solución en la escala gruesa. La solución de esta escala fina se localiza en el interior de cada elemento finito y se supone u¯=0 en el contorno de los elementos.

El residuo de la ecuación de momento en la escala grande, R_h, resulta en:

La subescala de velocidad, u¯, se aproxima de distinta forma en cada método de estabilización.

En este trabajo se utilizarán el método ASGS y el método OSS en los problemas para flujo confinado. El método split-OSS se utiliza para problemas grandes, como en flujos con superficie libre.

En el método de las subescalas algebraicas (ASGS), u¯ se toma proporcional al residuo R_h

donde τ₁ es un parámetro numérico.

En el método de las subescalas ortogonales (OSS) la subescala u¯ se toma proporcional a la proyección ortogonal de dicho residuo R_h:

donde P_h es la proyección sobre el espacio de los elementos finitos y Ph⊥=I−Pn es la proyección ortogonal.

Comparando ambos métodos se observa que la diferencia radica en sustituir el término R_h de la versión ASGS por Ph⊥Rh en la versión OSS.

8.1

Método ASGS

De acuerdo con la formulación anterior, el problema discreto con linealización de Picard, integración temporal BDF1, la estabilización ASGS consiste en:

Hallar uhn+1 y phn+1 tales que

con los términos con segundas derivadas de las funciones de elementos finitos nulos para elementos lineales.

8.2

Método OSS

En el método OSS, en el caso de considerar distintas densidades, como por ejemplo en problemas con interfase entre fluidos inmiscibles, el residuo en puntos de integración de lados opuestos a la interfase varía fuertemente y en proporción a la densidad. En tales casos puede adoptarse una proyección modificada [45]:

Sustituyendo Rhn+1 y la proyección modificada para el método OSS queda que:

El problema discretizado con linealización de Picard, integración temporal BDF1 y estabilización OSS consiste en hallar uhn+1 y phn+1 tales que:

8.3

Método split-OSS

El método split-OSS [15,45] separa la proyección del término convectivo y de la presión en 2 proyecciones, lo que permite una convergencia a la solución más rápida que en los métodos ASGS y OSS. Esto lo hace ventajoso para la resolución de problemas en 3 D. Este método resulta en:

para las iteraciones i = 1,2… hasta la convergencia, es decir, hasta que u^n+1,i−1 ≈ u^n+1,i y phn+1,i≈phn+1,i−1 en la norma elegida.

8.4

Parámetros de estabilización

El parámetro τ₁ de las ecuaciones (39) y (40) se elige con el fin de obtener esquemas numéricos estables y velocidades de convergencia óptimas (ver [46]). Este parámetro se calcula para cada elemento Ω_e. Para τ₁ se toma:

donde h es la longitud del elemento e   y   ue es la norma de la velocidad en el elemento e, c₁ y c₂ son coeficientes a elegir, μ y ρ son la viscosidad dinámica y densidad del fluido, respectivamente.

Se recomiendan los valores de c₁ = 1 y c₂ = 2 [46].

9

Formulación matricial del problema

Se presenta a continuación la versión matricial para los métodos ASGS, OSS y split-OSS. Se utiliza linealización de Picard y discretización en tiempo BDF1.

9.1

Método ASGS

En la versión matricial para el método ASGS la proyección Phuhn+1⋅∇uhn+1+1ρ∇phn+1 se trata con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

En lo que sigue se usa la notación compacta (a,b)=∫Ωa⋅b dΩ. En esta notación, la proyección de la ecuación (49) es la solución de:

para todo vh*∈Vh*, donde Vh* es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

donde U y P son los vectores de las incógnitas nodales para la velocidad u y la presión p, respectivamente. F es el vector de fuerzas nodales.

La ecuación (51) corresponde a la ecuación (41). La ecuación (52) corresponde a la ecuación (42).

Si se denotan los índices nodales a, b, los índices espaciales con i, j, la función de forma de los nodos a por N^a y la función de forma de los nodos b por N^b, entonces las matrices de las ecuaciones anteriores, las cuales son válidas para los métodos restantes, son:

donde δ_ij es la delta de Kronecker.

9.2

Método OSS

En la versión matricial para el método OSS la proyección Pkukn+1⋅∇ukn+1+1ρ∇phn+1 también se trata con un ciclo iterativo, al igual que para la linealización del término convectivo. Se definen:

En notación compacta, la proyección de la ecuación (54) es la solución de:

para todo vh*∈Vh*, donde Vh* es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

donde U, P y Y son los vectores de las incógnitas nodales para la velocidad u, la presión p y la proyección y, respectivamente. F es el vector de fuerzas nodales. Solo se consideran las fuerzas nodales del residuo en la proyección de los términos respectivos.

La ecuación (56) corresponde a la ecuación (44). La ecuación (57) corresponde a la ecuación (45). La ecuación (58) es la proyección del residuo en la ecuación de conservación de momentum.

Las matrices de las ecuaciones (56) a (58) no definidas anteriormente son:

9.3

Método split-OSS

Se presenta a continuación la versión matricial para el método split-OSS, algo más compleja que los métodos ASGS y OSS. Las proyecciones Phuhn+1⋅∇uhn+1 y Ph1ρ∇phn+1 también se tratan con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

y

En la notación compacta, las proyecciones de las ecuaciones (60) y (61) son la solución de:

para todo vh*∈Vh*, donde Vh* es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

donde U, P, Y y Z son los vectores de las incógnitas nodales para la velocidad u, la presión p y las proyecciones y y z, respectivamente. F es el vector de fuerzas nodales.

La ecuación (64) corresponde a la ecuación (46). La ecuación (65) corresponde a la ecuación (47). Las ecuaciones (66) y (67) son las proyecciones de los residuos de la ecuación de conservación de momentum y la ecuación de incompresibilidad, respectivamente.

La matriz de las ecuaciones (64) a (67) no definida anteriormente es:

10

Conclusiones

En este trabajo se presentan el modelo continuo y la correspondiente formulación discreta para la resolución de las ecuaciones de Navier-Stokes para los flujos viscoplásticos de Bingham y de Herschel-Bulkley usando elementos finitos mixtos estabilizados con interpolación lineal tanto para la velocidad como para la presión.

Se tratan en detalle los fluidos viscoplásticos ideales y regularizados de Bingham y de Herschel-Bulkley. Se proponen, asimismo, nuevos modelos viscoplásticos regularizados para el fluido de Bingham y de Herschel-Bulkley.

Posteriormente, se presenta el modelo discreto de elementos finitos estabilizados para flujos confinados. Los métodos de estabilización usados son el método de estabilización de subescalas algebraicas (ASGS), el método de subescalas ortogonales (OSS) y el método de subescalas ortogonales desacopladas (split-OSS). El modelo discreto se ha extendido a los modelos viscoplásticos de Bingham y de Herschel-Bulkley.

Agradecimientos

Los autores agradecen el apoyo del Prof. R. Codina, mostrado a través de múltiples sugerencias y discusiones fructíferas, y a la Universidad de los Andes, en Venezuela, por su apoyo económico en el desarrollo de esta investigación, dentro de su programa de becarios en el exterior.