Revision as of 17:14, 24 May 2016 (view source) Scipediacontent (talk | contribs) (Created page with " <span id='centerContent'></span> <span id='centerInner'></span> <span id='tit0005'></span> =Elementos finitos mixtos estabilizados para flujos confinados de Bingham y de He...")	Revision as of 09:46, 25 May 2016 (view source) Scipediacontent (talk | contribs) m (Scipediacontent moved page Draft Content 671798508 to Moreno Cervera 2015) Newer edit →
(No difference)

---

Revision as of 09:46, 25 May 2016

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

Stabilized finite elements for Bingham and Herschel-Bulkley confined flows. Part I : Formulation

E. Moreno^a, , M. Cervera^{b, ,}

• ^aDepartamento de Ordenación de Cuencas, Ingeniería Forestal, Universidad de los Andes, ULA, Vía Chorros de Milla, 5001, Mérida, Venezuela
• ^bCentro Internacional de Métodos Numéricos en Ingeniería (CIMNE), Universidad Politécnica de Cataluña, UPC, Módulo C1, Campus Norte, Jordi Girona 1-3, 08034, Barcelona, España

Under a Creative Commons license

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.

Palabras clave

Método de los elementos finitos estabilizados ; Subescalas ortogonales OSS ; Incompresibilidad ; Fluidos viscoplásticos ; Modelo de Bingham ; Modelo de Herschel-Bulkley

Keywords

Stabilized finite elements ; Orthogonal sub-grid scales OSS ; Incompressibility ; Viscoplastic fluid ; Bingham model ; Herschel-Bulkley model

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] , [10]  and [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]  and [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] , [18] , [19] , [20]  and [21] como de mecánica de sólidos [19] , [22] , [23] , [24] , [25] , [26] , [27] , [28] , [29] , [30]  and [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 ${\textstyle \mathbb {R} ^{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:

${\displaystyle \rho (\partial _{t}{\mbox{u}}+{\mbox{u}}\cdot \nabla {\mbox{u}})-\nabla \cdot \sigma ={\mbox{f}}}$

( 1)

con

${\displaystyle \nabla \cdot \sigma =-\nabla p+\nabla \cdot \tau }$

( 2)

y

${\displaystyle \nabla \cdot {\mbox{u}}=0}$

( 3)

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:

${\displaystyle {\mbox{u}}={\bar {\mbox{u}}}{\mbox{   }}{\mbox{en}}{\mbox{   }}\Gamma _{d}\times [0,t]}$

( 4)

Condiciones de Neumann:

${\displaystyle {\mbox{n}}\cdot \sigma ={\mbox{t}}{\mbox{   }}{\mbox{en}}{\mbox{   }}\Gamma _{n}\times [0,t]}$

( 5)

donde n es el vector unitario normal al contorno ∂ Ω. Por simplicidad se tomará en el contorno Γ_d la velocidad ${\textstyle {\bar {\mbox{u}}}=0,t\in [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:

${\displaystyle \nabla ^{2}\phi (x,y)={\frac {\partial u_{x}}{\partial x}}-{\frac {\partial u_{y}}{\partial y}}}$

( 6)

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:

${\displaystyle \sigma =-p{\mbox{I}}+\tau }$

( 7)

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:

${\displaystyle \tau =2\mu \epsilon ({\mbox{u}})}$

( 8)

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:

${\displaystyle \epsilon ({\mbox{u}})={\frac {1}{2}}(\nabla {\mbox{u}}+\nabla {\mbox{u}}^{t})=\nabla ^{s}{\mbox{u}}}$

( 9)

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:

${\displaystyle {\dot {\gamma }}={\sqrt {2\epsilon ({\mbox{u}}):\epsilon ({\mbox{u}})}}}$

( 10)

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

${\displaystyle \tau ={\sqrt {{\frac {1}{2}}J_{2}}}={\sqrt {{\frac {1}{2}}(\tau :\tau )}}}$

( 11)

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

${\displaystyle \sigma =-p{\mbox{I}}+2\mu \epsilon ({\mbox{u}})}$

( 12)

Dependiendo de los valores de la viscosidad en función de la velocidad de deformación, ${\textstyle \mu =\mu \left({\dot {\gamma }}\right)}$ , 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:

${\displaystyle {\begin{array}{c}\mu \left({\dot {\gamma }}\right)=\mu _{o}+{\frac {\tau _{y}}{\dot {\gamma }}}{\mbox{ }}para{\mbox{ }}\tau >\tau _{y}\\{\dot {\gamma }}=0{\mbox{ }}para{\mbox{ }}\tau \leq \tau _{y}\end{array}}}$

( 13)

donde μ₀ es la viscosidad plástica, y la viscosidad aparente ${\textstyle \mu ({\dot {\gamma }})}$ disminuye con el incremento en la magnitud de la velocidad de deformación ${\textstyle {\dot {\gamma }}}$ ; τ es la magnitud del tensor de tensiones desviadoras.

En consecuencia, el tensor de tensiones desviadoras es:

${\displaystyle {\begin{array}{c}\tau =2\left(\mu _{0}+{\frac {\tau _{y}}{\dot {\gamma }}}\right)\epsilon ({\mbox{u}}){\mbox{ }}para{\mbox{ }}\tau >\tau _{y}\\{\dot {\gamma }}=0{\mbox{ }}para{\mbox{ }}\tau \leq \tau _{y}\end{array}}}$

( 14)

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:

${\displaystyle Bn={\frac {\tau _{y}H}{\mu V}}}$

( 15)

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 ${\textstyle \left({\dot {\gamma }}_{0}={\frac {V}{H}}\right)}$ .

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, ${\textstyle \tau _{y}^{\ast }}$ , definida por Papanastasiou [7] como:

${\displaystyle \tau _{y}^{\ast }={\frac {\tau _{y}H}{\mu V_{N}}}}$

( 16)

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:

${\displaystyle Re=\rho {\frac {V_{B}H}{\mu }}}$

( 17)

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:

${\displaystyle {\begin{array}{c}\mu \left({\dot {\gamma }}\right)=k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}{\mbox{ }}para{\mbox{   }}\tau >\tau _{y}\\{\mbox{ }}{\dot {\gamma }}=0{\mbox{ }}para{\mbox{   }}\tau \leq \tau _{y}\end{array}}}$

( 18)

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:

${\displaystyle {\begin{array}{c}\tau =2\left(k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}\right)\epsilon ({\mbox{u}}){\mbox{ }}para{\mbox{   }}\tau >\tau _{y}\\{\dot {\gamma }}=0{\mbox{ }}para{\mbox{   }}\tau \leq \tau _{y}\end{array}}}$

( 19)

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:

${\displaystyle Bn\ast =Od={\frac {\tau _{y}}{k\left({\frac {V}{H}}\right)^{n}}}={\frac {\tau _{y}}{k}}\left({\frac {H}{V}}\right)^{n}}$

( 20)

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]  and [38] :

${\displaystyle Re=\rho {\frac {VH}{k\left({\frac {V}{H}}\right)^{n-1}}}=\rho {\frac {V^{2}}{k}}\left({\frac {H}{V}}\right)^{n}}$

( 21)

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 ${\textstyle \left(lim{\mbox{ }}{\underset {{\dot {\gamma }}\rightarrow 0}{\mu }}\rightarrow \infty \right)}$ .

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]  and [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:

${\displaystyle \mu ({\dot {\gamma }})=\left\{{\begin{array}{c}\mu _{0}+{\frac {\tau _{\gamma }}{\dot {\gamma }}}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}>{\dot {\gamma }}_{c}\\\mu _{r}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}\leq {\dot {\gamma }}_{c}\end{array}}\right.}$

( 22)

donde μ_r es la viscosidad crítica y ${\textstyle {\dot {\gamma }}_{c}}$ es la velocidad de deformación crítica. Para esta velocidad de deformación crítica la tensión es:

${\displaystyle \mu _{r}{\dot {\gamma }}_{c}=\mu _{0}{\dot {\gamma }}_{c}+\tau _{y}}$

( 23)

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

${\displaystyle {\dot {\gamma }}_{c}={\frac {\tau _{y}}{\mu _{r}-\mu _{0}}}}$

( 24)

El valor de μ_r debe ser grande para aproximar el modelo ideal. Una buena aproximación recomendada por Beverly y Tanner es tomar ${\textstyle 300\leq {\frac {\mu _{r}}{\mu _{0}}}\leq 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:

${\displaystyle \tau _{c}=\tau _{y}{\frac {\mu _{r}}{\mu _{r}-\mu _{o}}}>\tau _{y}}$

( 25)

5.2. Modelo de Papanastasiou

Uno de los intentos por solventar la limitación debida a la singularidad de la viscosidad para ${\textstyle {\dot {\gamma }}\rightarrow 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 ${\textstyle \mu ({\dot {\gamma }})}$ del modelo ideal de la manera:

${\displaystyle \mu \left({\dot {\gamma }}\right)=\mu _{o}+{\frac {\tau _{y}}{\dot {\gamma }}}\left(1-exp(-m{\dot {\gamma }})\right)}$

( 26)

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:

${\displaystyle \mu _{max}={\underset {{\dot {\gamma }}\rightarrow 0}{lim}}\mu ({\dot {\gamma }})=\mu _{0}+m\tau _{y}}$

( 27)

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 ${\textstyle {\dot {\gamma }}<{\dot {\gamma }}_{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:

${\displaystyle \mu \left({\dot {\gamma }}\right)=\left(\mu _{0}+{\frac {\tau _{y}}{\dot {\gamma }}}\right)\left(1-exp(-{\frac {\eta _{0}}{\tau _{y}}}{\dot {\gamma }})\right)}$

( 28)

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:

${\displaystyle \mu _{max}={\underset {{\dot {\gamma }}\rightarrow 0}{lim}}\mu ({\dot {\gamma }})=m\tau _{y}}$

( 29)

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:

${\displaystyle \mu ({\dot {\gamma }})=\left\{{\begin{array}{c}k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}>{\dot {\gamma }}_{c}\\\mu _{r}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}\leq {\dot {\gamma }}_{c}\end{array}}\right.}$

( 30)

donde μ_r es la viscosidad crítica y ${\textstyle {\dot {\gamma }}_{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:

${\displaystyle \mu _{r}{\dot {\gamma }}_{c}=k{\dot {\gamma }}_{c}^{n}+\tau _{y}}$

( 31)

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:

${\displaystyle \mu ({\dot {\gamma }})=k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}\left(1-exp(-m{\dot {\gamma }})\right)}$

( 32)

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	${\textstyle {\underset {{\dot {\gamma }}\rightarrow 0}{lim}}k{\dot {\gamma }}^{n-1}}$	${\textstyle {\underset {{\dot {\gamma }}\rightarrow 0}{lim}}{\frac {\tau _{y}}{\dot {\gamma }}}[1-exp(-m{\dot {\gamma }})]}$	${\textstyle {\underset {{\dot {y}}\rightarrow 0}{lim}}\mu ({\dot {\gamma }})}$
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:

${\displaystyle \mu \left({\dot {\gamma }}\right)=\left(k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}\right)\left(1-exp(-m{\dot {\gamma }})\right)}$

( 33)

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	${\textstyle {\underset {{\dot {\gamma }}\rightarrow 0}{lim}}k{\dot {\gamma }}^{n-1}[1-exp(-m{\dot {\gamma }})]}$	${\textstyle {\underset {{\dot {\gamma }}\rightarrow 0}{lim}}{\frac {\tau _{y}}{\dot {\gamma }}}[1-exp(-m{\dot {\gamma }})]}$	${\textstyle {\underset {{\dot {\gamma }}\rightarrow 0}{lim}}\mu ({\dot {\gamma }})]}$
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 ${\textstyle {\dot {y}}=0}$ , esto es, μ_r  = mτ_y , en función del parámetro de regularización m ( fig. 4 ). Por tanto:

${\displaystyle \mu ({\dot {\gamma }})=\left\{{\begin{array}{c}\mu _{0}+{\frac {\tau _{y}}{\dot {\gamma }}}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}>{\dot {\gamma }}_{c}\\m\tau _{y}{\mbox{ }}para{\mbox{ }}{\dot {\gamma }}\leq {\dot {\gamma }}_{c}\end{array}}\right.}$

( 34)

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:

${\displaystyle {\dot {\gamma }}_{c}={\frac {\tau _{y}}{m\tau _{y}-\mu _{o}}}}$

( 35)

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:

${\displaystyle \mu ({\dot {\gamma }})=\left\{{\begin{array}{c}k{\dot {\gamma }}^{n-1}+{\frac {\tau _{y}}{\dot {\gamma }}}{\mbox{   }}para{\mbox{   }}{\dot {\gamma }}>{\dot {\gamma }}_{c}\\m\tau _{y}{\mbox{   }}para{\mbox{   }}{\dot {\gamma }}\leq {\dot {\gamma }}_{c}\end{array}}\right.}$

( 36)

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 ${\textstyle {\dot {\gamma }}={\dot {\gamma }}_{c}}$ ):

${\displaystyle m\tau _{y}{\dot {\gamma }}_{c}=k{\dot {\gamma }}_{c}^{n}+\tau _{y}}$

( 37)

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 ${\textstyle \mu =\mu ({\dot {\gamma }})}$ . 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]  and [44]. Estos métodos consisten en descomponer la solución, por ejemplo, la velocidad u en 2 componentes ${\textstyle {\mbox{u}}={\mbox{u}}_{h}+{\bar {\mbox{u}}}}$ ; una componente u_h , resuelta en la escala de la malla de elementos finitos considerada y una subescala ${\textstyle {\bar {\mbox{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 ${\textstyle {\bar {\mbox{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 ${\textstyle {\bar {\mbox{u}}}=0}$ en el contorno de los elementos.

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

${\displaystyle {\begin{array}{l}\displaystyle {\mbox{R }}_{h}=\int _{\Omega }[-\rho {\mbox{f }}+\rho \left(\partial _{t}{\mbox{u }}_{h}\right)+\rho \left({\mbox{u }}_{h}\cdot \nabla {\mbox{u }}_{h}\right)-2\mu \left({\dot {\gamma }}\right)\left(\nabla \cdot \nabla {\mbox{u }}_{h}\right)+\\\displaystyle +\nabla p_{h}]d\Omega \end{array}}}$

( 38)

La subescala de velocidad, ${\textstyle {\bar {\mbox{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), ${\textstyle {\bar {\mbox{u}}}}$ se toma proporcional al residuo R_h

${\displaystyle {\bar {\mbox{u}}}=-\tau _{1}{\mbox{R}}_{h}}$

( 39)

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

En el método de las subescalas ortogonales (OSS) la subescala ${\textstyle {\bar {\mbox{u}}}}$ se toma proporcional a la proyección ortogonal de dicho residuo R_h :

${\displaystyle {\bar {\mbox{u}}}=-\tau _{1}P_{h}^{\perp }{\mbox{R}}_{h}=-\tau _{1}({\mbox{R}}_{h}-P_{h}{\mbox{R}}_{h})}$

( 40)

donde P_h es la proyección sobre el espacio de los elementos finitos y ${\textstyle P_{h}^{\perp }={\mbox{I}}-P_{n}}$ 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 ${\textstyle P_{h}^{\perp }{\mbox{R}}_{h}}$ 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 ${\textstyle {\mbox{u}}_{h}^{n+1}{\mbox{ }}y{\mbox{ }}p_{h}^{n+1}}$ tales que

${\displaystyle {\begin{array}{l}\displaystyle \int _{\Omega }[{\frac {\rho }{\delta t}}\left({\mbox{u }}_{h}^{n+1}-{\mbox{u }}_{h}^{n}\right)\cdot {\mbox{v }}_{h}+\rho \left({\mbox{u }}_{h}^{n+1,i-1}\cdot \nabla {\mbox{u }}_{h}^{n+1,i}\right)\cdot {\mbox{v }}_{h}+\\\displaystyle +2\mu \left({\dot {\gamma }}\right)^{n+1,i}\nabla ^{s}{\mbox{u }}_{h}^{n+1,i}:\nabla ^{s}{\mbox{v }}_{h}-p_{h}^{n+1,i}\nabla \cdot {\mbox{v }}_{h}-{\mbox{f }}^{n+1}\cdot \\\displaystyle \cdot {\mbox{v }}_{h}]{\mbox{d }}\Omega +\sum _{e}\int _{\Omega ^{e}}\tau _{1}\left(\rho ({\mbox{u }}_{h}^{n+1,i-1}\cdot \nabla {\mbox{v }}_{h}\right)\cdot \left(-\rho {\mbox{f }}^{n+1,i}+\right.\\\left.\displaystyle +\rho \left({\mbox{u }}_{h}^{n+1,i}\cdot \nabla {\mbox{u }}_{h}^{n+1,i}\right)+\nabla p_{h}^{n+1,i}\right){\mbox{d }}\Omega =0\end{array}}}$

( 41)

${\displaystyle {\begin{array}{l}\displaystyle \int _{\Omega }[q_{h}\nabla .{\mbox{u }}_{h}^{n+1,i}]d\Omega +\sum _{e=1}^{n_{e}l}\int _{\Omega ^{e}}\tau _{1}^{n+1,i}\nabla q_{h}\cdot \\\displaystyle \cdot [\rho \left({\mbox{u }}_{h}^{n+1,i-1}\cdot \nabla \right){\mbox{ }}{\mbox{u }}_{h}^{n+1,i}+\nabla p_{h}^{n+1,i}-{\mbox{f }}]d\Omega =0\end{array}}}$

( 42)

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

${\displaystyle P_{h\rho }\left({\mbox{R}}_{h}^{n+1}\right)=\rho P_{h}\left({\frac {{\mbox{R}}_{h}^{n+1}}{\rho }}\right)}$

( 43)

Sustituyendo ${\textstyle {\mbox{R}}_{h}^{n+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 ${\textstyle {\mbox{u}}_{h}^{n+1}{\mbox{ }}y{\mbox{ }}p_{h}^{n+1}}$ tales que:

${\displaystyle {\begin{array}{c}\int _{\Omega }[{\frac {\rho }{\delta t}}({\mbox{u }}_{h}^{n+1}-{\mbox{u }}_{h}^{n})\cdot {\mbox{v }}_{h}d\Omega +2\mu \left({\dot {\gamma }}\right)^{n+1,i}\nabla ^{s}{\mbox{u }}_{h}^{n+1}:\nabla ^{s}{\mbox{v }}_{h}+\\\displaystyle +\rho ({\mbox{u }}_{h}^{n+1}.\nabla {\mbox{u }}_{h}^{n+1,i})\cdot {\mbox{v }}_{h}-p_{h}^{n+1}\nabla \cdot v_{h}-{\mbox{f }}^{n+1}\cdot v_{h}]d\Omega +\\\displaystyle +\sum _{e}\int _{\Omega ^{e}}\tau _{1}\rho \left({\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{v }}_{h}\right)\cdot \\[][(\rho {\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1}+\nabla p_{h}^{n+1}-{\mbox{f }}^{n+1})-\rho P_{h}({\mbox{u }}_{h}^{n+1}\cdot \\\displaystyle \cdot \nabla {\mbox{u }}_{h}^{n+1}+{\frac {\nabla p_{h}^{n+1}}{\rho }}-{\frac {{\mbox{f }}^{n+1}}{\rho }})]d\Omega =0\end{array}}}$

( 44)

${\displaystyle {\begin{array}{l}\displaystyle \int _{\Omega }[q_{h}\nabla \cdot {\mbox{u }}_{h}^{n+1}]d\Omega +\sum _{e}\int _{\Omega ^{e}}\tau _{1}\nabla q_{h}\cdot [(\rho {\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1}+\\\displaystyle +\nabla p_{h}^{n+1}-{\mbox{f }}^{n+1})-\rho P_{h}(({\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1})+{\frac {\nabla p_{h}^{n+1}}{\rho }}-\\\displaystyle -{\frac {{\mbox{f }}^{n+1}}{\rho }})]d\Omega =0\end{array}}}$

( 45)

8.3. Método split- OSS

El método split -OSS [15]  and [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:

${\displaystyle {\begin{array}{l}\displaystyle \int _{\Omega }\left[{\frac {\rho }{\delta t}}\left({\mbox{u }}_{h}^{n+1}-{\mbox{u }}_{h}^{n}\right)\cdot {\mbox{v }}_{h}d\Omega +2\mu ({\dot {\gamma }})^{n+1,i}\nabla ^{s}{\mbox{u }}_{h}^{n+1}:\nabla ^{s}{\mbox{v }}_{h}+\right.\\\left.\displaystyle +\rho \left({\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1,i}\right)\cdot {\mbox{v }}_{h}-p_{h}^{n+1}\nabla \cdot {\mbox{v }}_{h}-{\mbox{f }}^{n+1}\cdot {\mbox{v }}_{h}\right]d\Omega +\\\displaystyle +\sum _{e}\int _{\Omega ^{e}}\tau _{1}\rho \left({\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{v }}_{h}\right)\cdot \left[\left(\rho {\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1}\right)\right]-\\\displaystyle -\rho P_{h}\left({\mbox{u }}_{h}^{n+1}\cdot \nabla {\mbox{u }}_{h}^{n+1}\right)]d\Omega =0\end{array}}}$

( 46)

${\displaystyle {\begin{array}{l}\displaystyle \int _{\Omega }[q_{h}\nabla \cdot {\mbox{u }}_{h}^{n+1}]d\Omega +\sum _{e}\int _{\Omega ^{e}}\tau _{1}\nabla q_{h}\cdot [\nabla p_{h}^{n+1}-{\mbox{f }}^{n+1}]-\\\displaystyle -\rho P_{h}\left({\frac {\nabla p_{h}^{n+1}}{\rho }}-{\frac {{\mbox{f }}^{n+1}}{\rho }}\right)]d\Omega =0\end{array}}}$

( 47)

para las iteraciones i  = 1,2… hasta la convergencia, es decir, hasta que u^{n +1,i −1}  ≈ u^{n +1,i} y ${\textstyle p_{h}^{n+1,i}\approx p_{h}^{n+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:

${\displaystyle \tau _{1}=\left[c_{1}{\frac {\mu }{(h)^{2}}}+c_{2}{\frac {\rho \left|{\mbox{u}}^{e}\right|}{h}}\right]^{-1}}$

( 48)

donde h   es la longitud del elemento ${\textstyle e{\mbox{   }}{\mbox{y}}{\mbox{   }}\left|{\mbox{u}}^{e}\right|}$ 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 ${\textstyle P_{h}\left({\mbox{u}}_{h}^{n+1}\cdot \nabla {\mbox{u}}_{h}^{n+1}+{\frac {1}{\rho }}\nabla p_{h}^{n+1}\right)}$ se trata con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

${\displaystyle {\mbox{y}}_{h}^{n+1}=P_{h}\left({\mbox{u}}_{h}^{n+1}.\nabla {\mbox{u}}_{h}^{n+1}+{\frac {1}{\rho }}\left(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1}\right)\right)}$

( 49)

En lo que sigue se usa la notación compacta ${\textstyle (a,b)=\int _{\Omega }a\cdot b{\mbox{ }}{\mbox{d}}\Omega }$ . En esta notación, la proyección de la ecuación (49) es la solución de:

${\displaystyle \left({\mbox{y}}_{h}^{n+1},{\mbox{v}}_{h}^{\ast }\right)=\left({\mbox{u}}_{h}^{n+1}.\nabla {\mbox{u}}_{h}^{n+1}+{\frac {1}{\rho }}\left(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1}\right),{\mbox{v}}_{h}^{\ast }\right)}$

( 50)

para todo ${\textstyle {\mbox{v}}_{h}^{\ast }\in {\mbox{V}}_{h}^{\ast }}$ , donde ${\textstyle {\mbox{V}}_{h}^{\ast }}$ es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

${\displaystyle {\mbox{M}}{\frac {1}{\delta t}}{\mbox{U}}^{n+1}+{\mbox{K}}({\mbox{U}}^{n+1}){\mbox{U}}^{n+1}+{\mbox{G}}{\mbox{P}}^{n+1}+{\mbox{S}}_{u}(\tau _{1};{\mbox{U}}^{n+1}){\mbox{U}}^{n+1}={\mbox{F}}^{n+1}}$

( 51)

${\displaystyle {\mbox{D}}{\mbox{U}}^{n+1}+{\mbox{S}}_{p}\left(\tau _{1}\right){\mbox{ }}{\mbox{P}}^{n+1}=0}$

( 52)

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:

${\displaystyle {\mbox{M}}_{ij}^{ab}=(N^{a},\rho N^{b})\delta _{ij}}$

${\displaystyle {\begin{array}{l}\displaystyle {\mbox{K }}({\mbox{U }}^{n+1})_{ij}^{ab}=(N^{a},\rho {\mbox{u }}_{h}^{n+1}\cdot \nabla N^{b})\delta _{ij}+(\nabla N^{a}\\\displaystyle 2\mu \nabla ^{s}N^{b})\delta _{ij}\end{array}}}$

${\displaystyle {\mbox{G}}_{i}^{ab}=(N^{a},\partial _{i}N^{b})}$

${\displaystyle {\mbox{S}}_{u}(\tau _{1};{\mbox{ }}{\mbox{U}}^{n+1})_{ij}^{ab}=(\tau _{1}{\mbox{u}}_{h}^{n+1}\cdot \nabla N^{a},\rho {\mbox{u}}_{h}^{n+1}\cdot \nabla N^{b})\delta _{ij}}$

${\displaystyle {\mbox{D}}_{j}^{ab}=(N^{a},\partial _{j}N^{b})}$

${\displaystyle {\mbox{S}}_{p}\left(\tau _{1}\right)^{ab}=\left(\tau _{1}{\mbox{u}}_{h}^{n+1}\cdot \nabla N^{a},\nabla N^{b}\right)}$

${\displaystyle {\mbox{F}}_{i}^{a}=\left(N^{a},{\mbox{f}}_{i}\right)}$

( 53)

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 ${\textstyle P_{k}\left({\mbox{u}}_{k}^{n+1}\cdot \nabla {\mbox{u}}_{k}^{n+1}+{\frac {1}{\rho }}\nabla p_{h}^{n+1}\right)}$ también se trata con un ciclo iterativo, al igual que para la linealización del término convectivo. Se definen:

${\displaystyle {\mbox{y}}_{h}^{n+1}=P_{h}\left({\mbox{u}}_{h}^{n+1}.\nabla {\mbox{u}}_{h}^{n+1}+{\frac {1}{\rho }}\left(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1}\right)\right)}$

( 54)

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

${\displaystyle \left({\mbox{y}}_{h}^{n+1},{\mbox{v}}_{h}^{\ast }\right)=\left({\mbox{u}}_{h}^{n+1}.\nabla {\mbox{u}}_{h}^{n+1}+{\frac {1}{\rho }}\left(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1}\right),{\mbox{v}}_{h}^{\ast }\right)}$

( 55)

para todo ${\textstyle {\mbox{v}}_{h}^{\ast }\in {\mbox{V}}_{h}^{\ast }}$ , donde ${\textstyle {\mbox{V}}_{h}^{\ast }}$ es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

${\displaystyle {\begin{array}{l}\displaystyle {\mbox{M }}{\frac {1}{\delta t}}{\mbox{U }}^{n+1}+{\mbox{K }}({\mbox{U }}^{n+1}){\mbox{U }}^{n+1}+{\mbox{G }}{\mbox{P }}^{n+1}+{\mbox{S }}_{u}(\tau _{1}\\\displaystyle {\mbox{ }}{\mbox{U }}^{n+1}){\mbox{U }}^{n+1}-{\mbox{S }}_{y}(\tau _{1};{\mbox{ }}{\mbox{U }}^{n+1}){\mbox{Y }}^{n+1}={\mbox{F }}^{n+1}\end{array}}}$

( 56)

${\displaystyle {\mbox{D}}{\mbox{U}}^{n+1}+{\mbox{S}}_{p}\left(\tau _{1}\right){\mbox{ }}{\mbox{P}}^{n+1}-{\mbox{S}}_{z}\left({\mbox{U}}^{n+1}\right){\mbox{Y}}^{n+1}=0}$

( 57)

${\displaystyle {\mbox{C}}({\mbox{U}}^{n+1}){\mbox{U}}^{n+1}+{\mbox{G}}_{\pi }{\mbox{P}}^{n+1}=0}$

( 58)

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:

${\displaystyle {\mbox{G}}_{\pi i}^{ab}=(N^{a},\partial _{i}N^{b}/\rho )}$

${\displaystyle {\mbox{S}}_{y}\left(\tau _{1}{\mbox{u}}_{h}^{n+1}\right)_{ij}^{ab}=\left(\tau _{1}{\mbox{u}}_{h}^{n+1}\cdot \nabla N^{a},\rho N^{b}\right)\delta _{ij}}$

${\displaystyle {\mbox{S}}_{z}\left(\tau _{1}\right)_{j}^{ab}=(\tau _{1}\partial _{j}N^{a},\rho N^{b})}$

${\displaystyle {\mbox{C}}({\mbox{U}}^{n+1})_{ij}^{ab}=(N^{a},{\mbox{u}}_{h}^{n+1}\cdot \nabla N^{b}){\mbox{ }}\delta _{ij}}$

( 59)

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 ${\textstyle P_{h}\left({\mbox{u}}_{h}^{n+1}\cdot \nabla {\mbox{u}}_{h}^{n+1}\right){\mbox{ }}y{\mbox{ }}P_{h}\left({\frac {1}{\rho }}\nabla p_{h}^{n+1}\right)}$ también se tratan con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

${\displaystyle {\mbox{y}}_{h}^{n+1}=P_{h}\left({\mbox{u}}_{h}^{n+1}\cdot \nabla {\mbox{u}}_{h}^{n+1}\right)}$

( 60)

y

${\displaystyle {\mbox{z}}_{h}^{n+1}=P_{h}\left({\frac {1}{\rho }}(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1})\right)}$

( 61)

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

${\displaystyle \left({\mbox{y}}_{h}^{n+1},{\mbox{v}}_{h}^{\ast }\right)=\left({\mbox{u}}_{h}^{n+1}\cdot \nabla {\mbox{u}}_{h}^{n+1},{\mbox{v}}_{h}^{\ast }\right)}$

( 62)

${\displaystyle \left({\mbox{z}}_{h}^{n+1},{\mbox{v}}_{h}^{\ast }\right)=\left({\frac {1}{\rho }}\left(\nabla p_{h}^{n+1}-{\mbox{f}}^{n+1}\right),{\mbox{v}}_{h}^{\ast }\right)}$

( 63)

para todo ${\textstyle {\mbox{v}}_{h}^{\ast }\in {\mbox{V}}_{h}^{\ast }}$ , donde ${\textstyle {\mbox{V}}_{h}^{\ast }}$ es el espacio V_h ampliado con los vectores de funciones continuas asociados a los nodos del contorno.

El sistema algebraico resultante es:

${\displaystyle {\begin{array}{l}\displaystyle {\mbox{M }}{\frac {1}{\delta t}}{\mbox{U }}^{n+1}+{\mbox{K }}({\mbox{U }}^{n+1}){\mbox{U }}^{n+1}+{\mbox{G }}{\mbox{P }}^{n+1}+{\mbox{S }}_{u}(\tau _{1}\\\displaystyle {\mbox{ }}{\mbox{U }}^{n+1}){\mbox{U }}^{n+1}-{\mbox{S }}_{y}(\tau _{1};{\mbox{ }}{\mbox{U }}^{n+1}){\mbox{Y }}^{n+1}={\mbox{F }}^{n+1}\end{array}}}$

( 64)

${\displaystyle {\mbox{D}}{\mbox{U}}^{n+1}+{\mbox{S}}_{p}(\tau _{1}){\mbox{P}}^{n+1}-{\mbox{S}}_{z}(\tau _{1}){\mbox{Z}}^{n+1}=0}$

( 65)

${\displaystyle {\mbox{M}}_{\pi }{\mbox{Y}}^{n+1}-{\mbox{C}}({\mbox{U}}^{n+1}){\mbox{U}}^{n+1}=0}$

( 66)

${\displaystyle {\mbox{M}}_{\pi }{\mbox{Z}}^{n+1}-{\mbox{G}}_{\pi }{\mbox{P}}^{n+1}=0}$

( 67)

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:

${\displaystyle {\mbox{M}}_{\pi ij}^{ab}=(N^{a},N^{b})\delta _{ij}}$

( 68)

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.

Bibliografía

1. [1] E. Bingham; Fluidity and Plasticity; McGraw-Hill, New York (1922)
2. [2] Herschel W., & Bulkley R. (1926). «Measurement of consistency as applied to rubber-benzene solutions». Proceding of American Society of Testing Material, 26, part. II, 621-633.
3. [3] M. Bercovier, M. Engelman; A finite-element method for incompressible non-Newtonian flows; J. Comput. Phys., 36 (1980), pp. 313–326
4. [4] Tanner R. & Milthorpe. (1983). «Numerical simulation of flow fluids with yield stress». Num. meth. lam. turb. flow. En: Proc. 3 rd Int. Conf., Scattle, Swasea, UK.
5. [5] A. Beris, J. Tsamopoulos, R. Armstrong, R. Brown; Creeping motions of sphere through a Bingham Plastic; J. Fluid Mech., 158 (1985), pp. 219–244
6. [6] R. Von Mises; Mechanik der festen Korper im plastisch deformablen Zustand; Gottinger Nachr, math-phys Kl (1913), pp. 582–592
7. [7] T. Papanastasiou; Flow of material with yield; J. Rheol., 36 (1987), pp. 389–407
8. [8] D. Dall’Onder dos Santos, S. Frey, M. Naccache, P. Souza Mendes; Numerical approximations for flow of viscoplastic fluid in a lid-driven cavity; Non-Newtonian Fluid Mech., 166 (2011), pp. 667–679
9. [9] S. Abdali, E. Mitsoulis; Entry and exit flows of Bingham fluids; J. Rheol., 36 (2) (1992)
10. [10] E. Mitsoulis, T. Zisis; Flow of Bingham plastics in a lid-driven square cavity; J. Non-Newtonian Fluid Mech., 101 (2001), pp. 173–180
11. [11] O. Zienkiewicz, P. Jain, E. Oñate; Flow of solids during forming and extrusion: Some aspect of numerical solutions; Int. J. Solids Struct., 14 (1978), pp. 15–38
12. [12] F. Brezzi, M. Fortin; Mixed and Hibrid Finite Element Methods; Springer Series in Computational Mathematics 15, Springer, New York (1991)
13. [13] R. Codina; Comparison of some finite element methods for solving the diffusion-convection-reaction equation; Comput. Methods Appl. Mech. Eng., 156 (1998), pp. 185–210
14. [14] T. Hughes, G. FeijÓo, L. Mazzei, J. Quincy; The variational multiscale method-A paradigm for computational mechanics; Comput. Methods Appl. Mech. Eng., 166 (1998), pp. 3–24
15. [15] R. Codina; On stabilized finite element methods for linear systems of convection-diffusion-reaction equations; Comput. Methods Appl. Mech. Eng., 188 (2000), pp. 61–82
16. [16] Principe J. (2008). «Subgrid scale stabilizad finite elements for low speed flows». Ph. D thesis, Universidad Politécnica de Cataluña.
17. [17] R. Codina; Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods; Comput. Methods Appl. Mech. Eng., 190 (2000), pp. 1579–1599
18. [18] Coppola, H. & Codina R. (2005). «An improved level-set approach using finite elements wiyh discontinuous gradient pressure shape functios». En: Computational Methods in Marine Engineering, Oslo, Norway.
19. [19] Chiumenti, M., Cervera, M. & Codina, R. (2013). «A mixed three-field FE formulation for stress accurate analysis including the incompressible limit». Submitted to Comput. Methods Appl. Mech. Eng.
20. [20] O. Guash, R. Codina; An algebraic subgrid scale finite element method for the convected Helmholtz equation in two dimensions with applications in aeroacoustics; Comput. Methods Appl. Mech. Eng., 31 (2007), pp. 919–930
21. [21] R. Planas, S. BadÍa, R. Codina; Aproximation of the inductionless MHD problem using a stabilized finite element method; J. Comput. Phys., 230 (2011), pp. 2977–2996
22. [22] M. Cervera, M. Chiumenti; Size effect and localization in J2 plasticity; Int. J. Solids Struct., 46 (2009), pp. 3301–3312
23. [23] M. Cervera, M. Chiumenti, R. Codina; Mixed stabilized finite element methods in nonlinear solid mechanics. Part I: Formulation; Comput. Methods Appl. Mech. Eng., 199 (2010), pp. 2559–2570
24. [24] M. Cervera, M. Chiumenti, R. Codina; Mixed stabilized finite element methods in nonlinear solid mechanics. Part II: Strain localization; Comput. Methods Appl. Mech. Eng., 199 (2010), pp. 2571–2589
25. [25] M. Cervera, M. Chiumenti, R. Codina; Mesh objetive modeling of cracks using continuos linear strain and displacement interpolations; Int. J. Numer. Meth. Eng., 87 (10) (2011), pp. 962–987
26. [26] M. Cervera, M. Chiumenti, D. di Capua; Benchmarking on bifurcation and localization in J2 plasticity for plane strain conditions; Comput. Methods Appl. Mech. Eng. (2012) 241-244, 206-224
27. [27] M. Cervera, M. Chiumenti, Q. Valverde, A. de Saracibar; A mixed linear/linear simplicial elements for incompressible elasticity and plasticity; Comput. Methods Appl. Mech. Eng., 192 (2003), pp. 5249–5263
28. [28] M. Cervera, M. Chiumenti, A. Saracibar; Shear band localization via local J2 continuum damage mechanics; Comput. Methods Appl. Mech. Eng., 193 (2004), pp. 849–880
29. [29] M. Cervera, A. de Saracibar; Softening, localization and stabilization: Capture of discontinuous solutions in J2 plasticity; Int. J. Numer. Anal. Met., 28 (2004), pp. 373–393
30. [30] M. Chiumenti, Q. Valverde, A. de Saracibar, M. Cervera; A stabilized formulations for incompressible elasticity using linear displacement and pressure interpolations; Comput. Methods Appl. Mech. Eng., 191 (2002), pp. 5253–5264
31. [31] A. De Saracibar, M. Chiumenti, Q. Valverde, M. Cervera; On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity; Comput. Methods Appl. Mech. Eng., 195 (9-12) (2006), pp. 1224–1251
32. [32] R. Bird, G. Dai, B. Yarusso; Rheology and flow of viscoplastic materials; Rev. Chem. Eng., 1 (1983), pp. 1–70
33. [33] J. Oldroyd; A rational formulation of the equations of plastic flow for a Bingham solid; Proc. Camb. Philos. Soc., 43 (1947), p. 100
34. [34] E. Reiner; Handbuch der Phisik; Springer, Berlin (1958)
35. [35] W. Prager; Introduction to Mechanics of Continua; Ginn and Co, Boston, MA (1961)
36. [36] P. Scott, F. Mirza, J. Vlachopoulos; Finite-element simulation of laminar viscoplastic flows with regions of recirculations; J. Rheol., 32 (1988), pp. 387–400
37. [37] E. Mitsoulis, L.R. Huilgol; Entry flows of Bingham plastic in expansions; J. Non-Newtonian Fluid Mech., 122 (2003), pp. 45–54
38. [38] S. Panda, R. Chhabra; Laminar flow of power-law fluids past a rotating cylinder; Jl. Non-Newtonian Fluid Mech., 165 (2010), pp. 1442–1461
39. [39] R. Bird, R. Armstrong, O. Hassager; Dynamics of Polymeric Liquids: I. Fluid Mechanics; (2 nd ed.)Jhon Wiley and Sons, New York (1987)
40. [40] D. Peric, S. Slijecpcevic; Computational modelling of viscoplastic fluids based on a stabilized finite element method; Engineering Computations, 18 (2001), pp. 577–591
41. [41] V. Kelessidis, R. Maglione, C. Tsamantaki, Y. Aspirtakis; Optimal determination of rheological parameters for Herschel-Bulkley drilling fluids and impact on pressure drop, velocity profiles and penetration rates during drilling; J. Petrol. Sci. Eng., 53 (2006), pp. 203–224
42. [42] Westerberg L., Lundström T., Höglund E. & Lugt P. (2010). «Investigation of grease flow in a rectangular channel including wall slip effects using microparticle image velocimetry». Tribology Transaccions.
43. [43] A. Brooks, T. Hughes, A. Russo; Streamlines upwind/Petrov-Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier-Stokes equation; Comput. Methods Appl. Mech. Eng., 32 (1982), pp. 199–259
44. [44] T. Hughes, L. Franca, G. Hulbert; A new finite element formulations for computational fluid dynamics: VIII. The Galerkin/least-square method for advective-diffusive equations; Comput. Methods Appl. Mech. Eng., 73 (1989), pp. 173–189
45. [45] Coppola H. (2009). «A finite element model for free surface and two fluid flows on fixed meshes». Ph. D thesis, Universidad Politécnica de Cataluña.
46. [46] Codina R. (2000c). «Stabilized finite element approximation of transient incompresible flows using orthogonal subscale». Cimne, No. 197.