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

E. Moreno[[#aff0005|<sup>a</sup>]]<sup>, </sup> ,                             M. Cervera[[#aff0010|<sup>b</sup>]]<sup>, </sup><sup>, </sup>

<ul><li><span id='aff0005'></span>
<sup>a</sup>Departamento de Ordenación de Cuencas, Ingeniería Forestal, Universidad de los Andes, ULA, Vía Chorros de Milla, 5001, Mérida, Venezuela</li>
<li><span id='aff0010'></span>
<sup>b</sup>Centro 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</li>
</ul>

Under a Creative Commons [http://creativecommons.org/licenses/by-nc-nd/4.0/ license]

<span id='ppvPlaceHolder'></span>

<span id='refersToAndreferredToBy'></span>
<span id='referredToBy'></span>

<span id='authorabs00101'></span>
==Resumen==

<span id='spar0005'></span>
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.

<span id='spar0010'></span>
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.

<span id='spar0015'></span>
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.

<span id='spar0020'></span>
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).          

<span id='spar0025'></span>
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.

<span id='authorabs00052'></span>
==Abstract==

<span id='spar0030'></span>
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.

<span id='spar0035'></span>
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.

<span id='spar0040'></span>
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.

<span id='spar0045'></span>
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.

<span id='spar0050'></span>
The methodology proposed in this work provides a computational tool to study confined viscoplastic flows, common in industry.

<span id='kwd_1'></span>
==Palabras clave==

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

<span id='kwd_2'></span>
==Keywords==

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

<span id='embeddedPDF'></span>

<span id='SD_BA1P'></span>

<span id='sec0005'></span>
==1. Introducción==

<span id='par0005'></span>
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.

<span id='par0010'></span>
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.

<span id='par0015'></span>
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.

<span id='par0020'></span>
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).      

<span id='par0025'></span>
El modelo que dio inicio al estudio de los materiales viscoplásticos fue el modelo plástico de Bingham [[#bib0005|[1]]] , formulado por Eugene C. Bingham para describir el comportamiento de las pinturas. El modelo de Herschel-Bulkley [[#bib0190|[2]]]  se considera como un modelo generalizado de Bingham, aunque ha sido menos estudiado que el modelo de Bingham.      

<span id='par0030'></span>
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 [[#bib0015|[3]]] , Tanner y Milthorpe [[#bib0020|[4]]]  y Beris et al. [[#bib0025|[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 [[#bib0030|[6]]]  en las zonas de no fluencia y el modelo ideal de Bingham en la zona de fluencia. En 1987, Papanastasiou [[#bib0035|[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) [[#bib0040|[8]]]  han propuesto recientemente una modificación del modelo de Papanastasiou.      

<span id='par0035'></span>
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.

<span id='par0040'></span>
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) [[#bib0035|[7]]] , [[#bib0045|[9]]] , [[#bib0050|[10]]]  and [[#bib0055|[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.      

<span id='par0045'></span>
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 [[#bib0060|[12]]] . La formulación de igual interpolación lineal usada en este trabajo no cumple con la condición Babuška-Brezzi [[#bib0065|[13]]] .      

<span id='par0050'></span>
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])  [[#bib0070|[14]]]  para una ecuación escalar de difusión-reacción. Codina [[#bib0075|[15]]]  amplió esta aproximación algebraica aplicándola a sistemas escalares multidimensionales.      

<span id='par0055'></span>
Posteriormente, Codina [[#bib0075|[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  [[#bib0075|[15]]]  and [[#bib0080|[16]]] . La estabilización OSS ha sido reformulada en una nueva versión del método llamada estabilización ''split'' -OSS  [[#bib0085|[17]]] , computacionalmente más ventajosa. Actualmente se usan en problemas muy variados, tanto de mecánica de fluidos  [[#bib0075|[15]]] , [[#bib0085|[17]]] , [[#bib0090|[18]]] , [[#bib0095|[19]]] , [[#bib0100|[20]]]  and [[#bib0105|[21]]]  como de mecánica de sólidos  [[#bib0095|[19]]] , [[#bib0110|[22]]] , [[#bib0115|[23]]] , [[#bib0120|[24]]] , [[#bib0125|[25]]] , [[#bib0130|[26]]] , [[#bib0135|[27]]] , [[#bib0140|[28]]] , [[#bib0145|[29]]] , [[#bib0150|[30]]]  and [[#bib0155|[31]]] .      

<span id='par0060'></span>
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.      

<span id='sec0010'></span>
==2. Modelo continuo para el problema del flujo confinado==

<span id='par0065'></span>
El problema continuo de dinámica de fluidos incompresibles e isotérmicos puede resolverse completamente considerando las ecuaciones de Navier-Stokes.

<span id='par0070'></span>
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.

<span id='par0075'></span>
Considérese Ω, un dominio abierto y acotado dimensional de <math display="inline">\mathbb{R}^d</math> , 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) Γ<sub>''d''</sub>  = ∂Ω<sub>''d''</sub>  y el contorno con condiciones de Neumann (tracciones impuestas) Γ<sub>''n''</sub>  = ∂Ω<sub>''n''</sub>  de forma que Γ = Γ<sub>''d''</sub>  ∪ Γ<sub>''n''</sub> , [0, ''T'' ] es el intervalo de tiempo de análisis.      

<span id='par0080'></span>
El problema de Navier-Stokes consiste en encontrar una velocidad '''u'''  y una presión ''p''  tal que:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho (\partial _t\mbox{u}+\mbox{u}\cdot \nabla \mbox{u})-\nabla \cdot \sigma =\mbox{f}</math>
|}
| style="text-align: right;" | ( 1)
|}

con

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\nabla \cdot \sigma =-\nabla p+\nabla \cdot \tau </math>
|}
| style="text-align: right;" | ( 2)
|}

y

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\nabla \cdot \mbox{u}=0</math>
|}
| style="text-align: right;" | ( 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.      

<span id='par0085'></span>
A esas ecuaciones deben añadirse condiciones iniciales de la forma '''u'''  = '''u'''<sub>0</sub>  en Ω, ''t''<sub>0</sub>  = 0 y condiciones de contorno:      

<span id='par0090'></span>
Condiciones de Dirichlet:

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{u}=\bar{\mbox{u}}\mbox{   }\mbox{en}\mbox{   }\Gamma _d\times [0,t]</math>
|}
| style="text-align: right;" | ( 4)
|}

<span id='par0095'></span>
Condiciones de Neumann:

<span id='eq0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{n}\cdot \sigma =\mbox{t}\mbox{   }\mbox{en}\mbox{   }\Gamma _n\times [0,t]</math>
|}
| style="text-align: right;" | ( 5)
|}

donde '''n'''  es el vector unitario normal al contorno ''∂'' Ω. Por simplicidad se tomará en el contorno Γ<sub>''d''</sub>  la velocidad <math display="inline">\bar{\mbox{u}}=0,t\in [0,T]</math> . El vector '''t'''  es el vector de tracción sobre el contorno con condiciones de Neumann.      

<span id='sec0015'></span>
==3. Líneas de corriente==

<span id='par0100'></span>
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.

<span id='par0105'></span>
Las líneas de corriente para un flujo bidimensional con un campo de velocidades '''u'''  = (''u''<sub>''x''</sub>'',u''<sub>''y''</sub> ) coinciden con las líneas de nivel de la función ''ϕ'' , solución de la ecuación laplaciana:

<span id='eq0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\nabla ^2\phi (x,y)=\frac{\partial u_x}{\partial x}-\frac{\partial u_y}{\partial y}</math>
|}
| style="text-align: right;" | ( 6)
|}

con la condición de contorno ϕ = 0.      

<span id='sec0020'></span>
==4. Ecuación constitutiva para fluidos viscoplásticos==

<span id='par0110'></span>
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.

<span id='par0115'></span>
El tensor de tensiones '''σ'''  se descompone en su parte volumétrica y desviadora como:

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sigma =-p\mbox{I}+\tau </math>
|}
| style="text-align: right;" | ( 7)
|}

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

<span id='par0120'></span>
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:

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau =2\mu \epsilon (\mbox{u})</math>
|}
| style="text-align: right;" | ( 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:

<span id='eq0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\epsilon (\mbox{u})=\frac{1}{2}(\nabla \mbox{u}+\nabla \mbox{u}^t)=\nabla ^s\mbox{u}</math>
|}
| style="text-align: right;" | ( 9)
|}

donde ▿'''u'''  es el gradiente de la velocidad y (▿'''u'''<sup>t</sup> ) es la transpuesta del mismo.      

<span id='par0125'></span>
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:

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\dot{\gamma }=\sqrt{2\epsilon (\mbox{u}):\epsilon (\mbox{u})}</math>
|}
| style="text-align: right;" | ( 10)
|}

<span id='par0130'></span>
La magnitud del tensor desviador o tensión efectiva, ''τ'' , se toma como la raíz del segundo invariante del tensor de las tensiones desviadoras:

<span id='eq0055'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau =\sqrt{\frac{1}{2}J_2}=\sqrt{\frac{1}{2}(\tau :\tau )}</math>
|}
| style="text-align: right;" | ( 11)
|}

<span id='par0135'></span>
De acuerdo con lo anterior, la ecuación (7) se puede escribir como:

<span id='eq0060'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\sigma =-p\mbox{I}+2\mu \epsilon (\mbox{u})</math>
|}
| style="text-align: right;" | ( 12)
|}

<span id='par0140'></span>
Dependiendo de los valores de la viscosidad en función de la velocidad de deformación, <math display="inline">\mu =\mu \left(\dot{\gamma }\right)</math> , 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. [[#bib0160|[32]]] . Se estudian en este trabajo los modelos para fluidos viscoplásticos, en particular para el modelo de Bingham y de Herschel-Bulkley.      

<span id='sec0025'></span>
===4.1. Modelo ideal de Bingham===

<span id='par0145'></span>
Eugene C. Bingham describió las pinturas con este modelo en 1919, publicado en su libro ''Fluidity and Plasticity''[[#bib0005|[1]]] . El modelo fue analizado por Oldroyd [[#bib0165|[33]]] , Reiner [[#bib0170|[34]]]  y Prager [[#bib0175|[35]]] . Los plásticos de Bingham requieren de una tensión de corte mínima, ''τ''<sub>y</sub> , a partir de la cual comienzan a moverse ([[#fig0005|fig. 1]] ).      

<span id='figure_fig0005'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr1.jpg|center|320px|Curvas reológicas. Modelo bilineal.]]


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

Figura 1.

<span id='spar0055'></span>
Curvas reológicas. Modelo bilineal.

</span>
|}

<span id='par0150'></span>
En el modelo de Bingham la viscosidad está dada por:

<span id='eq0065'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 13)
|}

donde ''μ''<sub>0</sub>  es la viscosidad plástica, y la viscosidad aparente <math display="inline">\mu (\dot{\gamma })</math>  disminuye con el incremento en la magnitud de la velocidad de deformación <math display="inline">\dot{\gamma }</math> ; ''τ''  es la magnitud del tensor de tensiones desviadoras.      

<span id='par0155'></span>
En consecuencia, el tensor de tensiones desviadoras es:

<span id='eq0070'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 14)
|}

<span id='par0160'></span>
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, ''τ''<sub>''y''</sub> . 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.      

<span id='sec0030'></span>
====4.1.1. Magnitudes adimensionales====

4.1.1.1. Número de Bingham

<span id='par0165'></span>
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.  [[#bib0160|[32]]] , se define como:

<span id='eq0075'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Bn=\frac{\tau _yH}{\mu V}</math>
|}
| style="text-align: right;" | ( 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.      

<span id='par0170'></span>
El número de Bingham relaciona la tensión de fluencia, ''τ''<sub>''y''</sub> , con la tensión ocasionada por una velocidad de deformación característica <math display="inline">\left(\dot{\gamma }_0=\frac{V}{H}\right)</math> .      

<span id='par0175'></span>
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''  → ∝ [[#bib0050|[10]]] .              

4.1.1.2. Tensión adimensional de fluencia

<span id='par0180'></span>
En flujos viscoplásticos es conveniente mostrar los resultados en función de una tensión de fluencia adimensional, <math display="inline">\tau _y^\ast </math> , definida por Papanastasiou [[#bib0035|[7]]]  como:

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau _y^\ast =\frac{\tau _yH}{\mu V_N}</math>
|}
| style="text-align: right;" | ( 16)
|}

donde ''V''<sub>''N''</sub>  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

<span id='par0185'></span>
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.  [[#bib0180|[36]]]  y Mitsoulis y Huilgol [[#bib0185|[37]]] , entre otros, usan en trabajos previos para flujo newtoniano. Para el flujo laminar de Bingham:

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Re=\rho \frac{V_BH}{\mu }</math>
|}
| style="text-align: right;" | ( 17)
|}

donde ''ρ''  es la densidad del fluido, ''V''<sub>''B''</sub>  es la velocidad promedio del flujo, ''H''  es una longitud característica y ''μ''  es la viscosidad para el fluido de Bingham.      

<span id='sec0050'></span>
===4.2. Modelo ideal de Herschel-Bulkley===

<span id='par0190'></span>
En el modelo plástico de Herschel-Bulkley [[#bib0190|[2]]]  se combinan la tensión de fluencia y la ley potencial. En el modelo de Herschel-Bulkley la viscosidad aparente está dada por:

<span id='eq0090'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 18)
|}

<span id='par0195'></span>
Al igual que para el modelo de Bingham, los materiales de Herschel-Bulkley requieren de una tensión de corte mínima, ''τ''<sub>''y''</sub> , 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)'' .      

<span id='par0200'></span>
El tensor desviador resulta:

<span id='eq0095'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 19)
|}

<span id='par0205'></span>
Si ''n''  ''=''  1, se tiene el fluido de Bingham como caso particular [[#bib0175|[35]]]  y el índice de consistencia es igual a la viscosidad plástica del material ''k''  = ''μ''<sub>0</sub> . Si la tensión de fluencia es nula, ''τ''<sub>''y''</sub>  = 0, se recupera la ley potencial.      

<span id='sec0055'></span>
====4.2.1. Magnitudes adimensionales====

4.2.1.1. Número generalizado de Bingham

<span id='par0210'></span>
Para materiales que obedecen el modelo de Herschel-Bulkley se define el número generalizado de Bingham, ''Bn'' * o número de Oldroyd  [[#bib0165|[33]]] , ''Od'' , como:

<span id='eq0100'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Bn\ast =Od=\frac{\tau _y}{k\left(\frac{V}{H}\right)^n}=\frac{\tau _y}{k}\left(\frac{H}{V}\right)^n</math>
|}
| style="text-align: right;" | ( 20)
|}

donde ''τ''<sub>''y''</sub>  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

<span id='par0215'></span>
El número de Reynolds usado en flujos de Herschel-Bulkley viene dado por la ley potencial [[#bib0020|[4]]]  and [[#bib0010|[38]]] :

<span id='eq0105'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Re=\rho \frac{VH}{k\left(\frac{V}{H}\right)^{n-1}}=\rho \frac{V^2}{k}\left(\frac{H}{V}\right)^n</math>
|}
| style="text-align: right;" | ( 21)
|}

<span id='sec0070'></span>
==5. Modelos viscoplásticos regularizados. Fluido de Bingham==

<span id='par0220'></span>
Los modelos reológicos ideales con tensión de fluencia presentan 2 problemas:<span id='list_lis0005'></span>
* 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 <math display="inline">\left(lim\mbox{ }\underset{\dot{\gamma }\rightarrow 0}{\mu }\rightarrow \infty \right)</math> .                  

<span id='par0235'></span>
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 [[#bib0035|[7]]] , [[#bib0195|[39]]]  and [[#bib0200|[40]]] .      

<span id='par0240'></span>
Para evitar estas dificultades y lograr una conveniente formulación computacional, se han propuestos diferentes modelos regularizados.

<span id='par0245'></span>
Un modelo muy utilizado es el modelo de doble viscosidad, inicialmente propuesto por Tanner y Milthorpe [[#bib0020|[4]]] . Otro de los modelos más usados en nuestros días es el modelo de Papanastasiou [[#bib0035|[7]]] . Souza Mendes y Dutra (SMD) [[#bib0040|[8]]]  han propuesto recientemente una modificación del modelo de Papanastasiou.      

<span id='sec0075'></span>
===5.1. Modelo de Tanner y Milthorpe===

<span id='par0250'></span>
El modelo introducido originalmente por Tanner y Milthorpe [[#bib0020|[4]]]  es un modelo con doble viscosidad lineal que regulariza el fluido de Bingham.      

<span id='par0255'></span>
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:

<span id='eq0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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.</math>
|}
| style="text-align: right;" | ( 22)
|}

donde ''μ''<sub>''r''</sub>  es la viscosidad crítica y <math display="inline">\dot{\gamma }_c</math>  es la velocidad de deformación crítica. Para esta velocidad de deformación crítica la tensión es:

<span id='eq0115'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu _r\dot{\gamma }_c=\mu _0\dot{\gamma }_c+\tau _y</math>
|}
| style="text-align: right;" | ( 23)
|}

<span id='par0260'></span>
Por tanto, la velocidad de deformación crítica es:

<span id='eq0120'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\dot{\gamma }_c=\frac{\tau _y}{\mu _r-\mu _0}</math>
|}
| style="text-align: right;" | ( 24)
|}

<span id='par0265'></span>
El valor de ''μ''<sub>''r''</sub>  debe ser grande para aproximar el modelo ideal. Una buena aproximación recomendada por Beverly y Tanner es tomar <math display="inline">300\leq \frac{\mu _r}{\mu _0}\leq 1000</math> .      

<span id='par0270'></span>
El inconveniente de este modelo es que para la viscosidad crítica se tiene una tensión crítica ''τ''<sub>''c''</sub>  algo mayor que la tensión de fluencia, como se muestra en la [[#fig0005|figura 1]] . Esta tensión crítica es:

<span id='eq0125'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau _c=\tau _y\frac{\mu _r}{\mu _r-\mu _o}>\tau _y</math>
|}
| style="text-align: right;" | ( 25)
|}

<span id='sec0080'></span>
===5.2. Modelo de Papanastasiou===

<span id='par0275'></span>
Uno de los intentos por solventar la limitación debida a la singularidad de la viscosidad para <math display="inline">\dot{\gamma }\rightarrow 0</math>  se debe a Papanastasiou [[#bib0035|[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.      

<span id='par0280'></span>
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 <math display="inline">\mu (\dot{\gamma })</math>  del modelo ideal de la manera:

<span id='eq0130'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu \left(\dot{\gamma }\right)=\mu _o+\frac{\tau _y}{\dot{\gamma }}\left(1-exp(-m\dot{\gamma })\right)</math>
|}
| style="text-align: right;" | ( 26)
|}

<span id='par0285'></span>
En la [[#fig0010|figura 2]]  se puede apreciar la influencia del parámetro de regularización ''m'' .      

<span id='figure_fig0010'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr2.jpg|center|301px|Modelo regularizado de Papanastasiou para el fluido de Bingham con diferentes ...]]


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

Figura 2.

<span id='spar0060'></span>
Modelo regularizado de Papanastasiou para el fluido de Bingham con diferentes valores del parámetro de regularización, m.

</span>
|}

<span id='par0290'></span>
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:

<span id='eq0135'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu _{max}=\underset{\dot{\gamma }\rightarrow 0}{lim}\mu (\dot{\gamma })=\mu _0+m\tau _y</math>
|}
| style="text-align: right;" | ( 27)
|}

<span id='par0295'></span>
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 ''μ''<sub>''t''</sub>  < ''μ''<sub>max</sub>  para velocidades de deformación muy bajas <math display="inline">\dot{\gamma }<\dot{\gamma }_c</math> .      

<span id='sec0085'></span>
===5.3. Modelo de Souza Mendes y Dutra (SMD)===

<span id='par0300'></span>
El modelo regularizado SMD de Souza Mendes y Dutra [[#bib0040|[8]]]  es similar al modelo de Papanastasiou, pero la regularización exponencial afecta a todos los términos de la viscosidad:

<span id='eq0140'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 28)
|}

<span id='par0305'></span>
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 η<sub>''o''</sub>  y la tensión de fluencia ''τ''<sub>''y''</sub> , ''m''  = ''η''<sub>0</sub> /''τ''<sub>''y''</sub> .      

<span id='par0310'></span>
El límite para la viscosidad cuando la velocidad de deformación es cero es:

<span id='eq0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu _{max}=\underset{\dot{\gamma }\rightarrow 0}{lim}\mu (\dot{\gamma })=m\tau _y</math>
|}
| style="text-align: right;" | ( 29)
|}

<span id='sec0090'></span>
==6. Modelos viscoplásticos regularizados. Fluido de Herschel-Bulkley==

<span id='par0315'></span>
Para el fluido de Herschel-Bulkley se proponen modelos regularizados análogos a los del fluido de Bingham.

<span id='sec0095'></span>
===6.1. Modelo de Tanner y Milthorpe===

<span id='par0320'></span>
El modelo de Tanner y Milthorpe [[#bib0020|[4]]]  para el fluido de Herschel-Bulkley es:

<span id='eq0150'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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.</math>
|}
| style="text-align: right;" | ( 30)
|}

donde ''μ''<sub>''r''</sub>  es la viscosidad crítica y <math display="inline">\dot{\gamma }_c</math>  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:

<span id='eq0155'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu _r\dot{\gamma }_c=k\dot{\gamma }_c^n+\tau _y</math>
|}
| style="text-align: right;" | ( 31)
|}

<span id='sec0100'></span>
===6.2. Modelo de Papanastasiou===

<span id='par0325'></span>
La regularización propuesta por Papanastasiou ([[#fig0015|fig. 3]] ) es también aplicable al modelo de Herschel-Bulkley. La viscosidad aparente queda definida como:

<span id='eq0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu (\dot{\gamma })=k\dot{\gamma }^{n-1}+\frac{\tau _y}{\dot{\gamma }}\left(1-exp(-m\dot{\gamma })\right)</math>
|}
| style="text-align: right;" | ( 32)
|}

<span id='figure_fig0015'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr3.jpg|center|292px|Modelo regularizado de Papanastasiou para el fluido de Herschel-Bulkley con ...]]


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

Figura 3.

<span id='spar0065'></span>
Modelo regularizado de Papanastasiou para el fluido de Herschel-Bulkley con diferentes valores del parámetro de regularización m.

</span>
|}

<span id='par0330'></span>
La influencia del parámetro ''m''  en el fluido de Herschel-Bulkley-Papanastasiou puede verse en la  [[#fig0015|figura 3]] .      

<span id='par0335'></span>
Al igual que para el modelo de Bingham, el modelo de Herschel-Bulkley requiere en la implementación numérica un valor de truncamiento ''μ''<sub>''t''</sub> .      

<span id='par0340'></span>
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  [[#tbl0005|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.      

<span id='table_tbl0005'></span>

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

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

Tabla 1.

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

|-

! Exponente
! colspan="2" | Términos de la viscosidad
! Viscosidad
|-

| ''n''
| <math display="inline">\underset{\dot{\gamma }\rightarrow 0}{lim}k\dot{\gamma }^{n-1}</math>
| <math display="inline">\underset{\dot{\gamma }\rightarrow 0}{lim}\frac{\tau _y}{\dot{\gamma }}[1-exp(-m\dot{\gamma })]</math>
| <math display="inline">\underset{\dot{y}\rightarrow 0}{lim}\mu (\dot{\gamma })</math>
|-

| Si ''n''  > 1                                                    
| 0
| ''mτ''<sub>''y''</sub>
| ''mτ''<sub>''y''</sub>
|-

| Si ''n''  = 1                                                    
| ''k = μ''
| ''mτ''<sub>''y''</sub>
| ''μ''  + ''mτ''<sub>''y''</sub>
|-

| Si ''n''  < 1                                                    
| ∞
| ''mτ''<sub>''y''</sub>
| ∞
|}

<span id='sec0105'></span>
===6.3. Modelo Souza-Mendez-Dutra===

<span id='par0345'></span>
El modelo SMD (Souza-Mendez-Dutra) regulariza el modelo ideal de Herschel-Bulkley en la forma:

<span id='eq0165'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mu \left(\dot{\gamma }\right)=\left(k\dot{\gamma }^{n-1}+\frac{\tau _y}{\dot{\gamma }}\right)\left(1-exp(-m\dot{\gamma })\right)</math>
|}
| style="text-align: right;" | ( 33)
|}

con ''m''  = ''η''<sub>0</sub> /''τ''<sub>''y''</sub>

<span id='par0355'></span>
El valor de la viscosidad cuando la velocidad de deformación tiende a cero se muestra en la [[#tbl0010|tabla 2]] . La viscosidad está acotada para cualquier valor de ''n'' , lo cual es una ventaja en la implementación numérica.      

<span id='table_tbl0010'></span>

<span id='tbl0010'></span>

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

Tabla 2.

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

|-

! Exponente
! colspan="2" | Términos de la viscosidad
! Viscosidad
|-

| ''n''
| <math display="inline">\underset{\dot{\gamma }\rightarrow 0}{lim}k\dot{\gamma }^{n-1}[1-exp(-m\dot{\gamma })]</math>
| <math display="inline">\underset{\dot{\gamma }\rightarrow 0}{lim}\frac{\tau _y}{\dot{\gamma }}[1-exp(-m\dot{\gamma })]</math>
| <math display="inline">\underset{\dot{\gamma }\rightarrow 0}{lim}\mu (\dot{\gamma })]</math>
|-

| Si ''n''  > 1                                                    
| 0
| ''mτ''<sub>''y''</sub>
| ''mτ''<sub>''y''</sub>
|-

| Si ''n''  = 1                                                    
| 0
| ''mτ''<sub>''y''</sub>
| ''mτ''<sub>''y''</sub>
|-

| Si ''n''  < 1                                                    
| 0
| ''mτ''<sub>''y''</sub>
| ''mτ''<sub>''y''</sub>
|}

<span id='par0360'></span>
Estos modelos presentados han formado parte de los estudios de diferentes problemas resueltos por Papanastasiou [[#bib0035|[7]]] , Kelessidis et al. [[#bib0205|[41]]] , Westerberg et al. [[#bib0210|[42]]]  y Dall’Onder dos Santos et al. [[#bib0040|[8]]] , entre otros.      

<span id='sec0110'></span>
==7. Modelos viscoplásticos regularizados propuestos==

<span id='par0365'></span>
Se proponen a continuación sendos modelos viscoplásticos regularizados para los fluidos de Bingham y de Herschel-Bulkley.

<span id='par0370'></span>
Ambos son modelos de doble viscosidad, basados en los modelos descritos anteriormente.

<span id='par0375'></span>
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).

<span id='sec0115'></span>
===7.1. Modelo regularizado de Bingham de doble viscosidad===

<span id='par0380'></span>
Este modelo es idéntico al modelo bilineal, pero la viscosidad crítica, ''μ''<sub>''r''</sub> , se toma igual al valor límite regularizado de Papanastasiou y SMD correspondiente a <math display="inline">\dot{y}=0</math> , esto es, ''μ''<sub>''r''</sub>  = ''mτ''<sub>''y''</sub> , en función del parámetro de regularización ''m''  ( [[#fig0020|fig. 4]] ). Por tanto:

<span id='eq0170'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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.</math>
|}
| style="text-align: right;" | ( 34)
|}

<span id='figure_fig0020'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr4.jpg|center|320px|Modelo ideal de Bingham y modelo regularizado Bingham-DV.]]


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

Figura 4.

<span id='spar0070'></span>
Modelo ideal de Bingham y modelo regularizado Bingham-DV.

</span>
|}

<span id='par0385'></span>
En este caso, el valor de la velocidad de deformación crítica es:

<span id='eq0175'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\dot{\gamma }_c=\frac{\tau _y}{m\tau _y-\mu _o}</math>
|}
| style="text-align: right;" | ( 35)
|}

<span id='par0390'></span>
En la [[#fig0025|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 ''τ''<sub>''y''</sub>  = 10 ''Pa''  ''y''  ''μ''<sub>0</sub>  = 0.2 ''Pas'' .      

<span id='figure_fig0025'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr5.jpg|center|346px|Comparación entre el modelo bilineal propuesto para el fluido de Bingham (DV) ...]]


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

Figura 5.

<span id='spar0075'></span>
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.

</span>
|}

<span id='sec0120'></span>
===7.2. Modelo regularizado de Herschel-Bulkley de doble viscosidad===

<span id='par0395'></span>
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, ''μ''<sub>''r''</sub>  = ''mτ''<sub>''y''</sub> , en función del parámetro de regularización ''m''  ( [[#fig0030|fig. 6]] ). Esto es:

<span id='eq0180'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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.</math>
|}
| style="text-align: right;" | ( 36)
|}

<span id='figure_fig0030'></span>

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

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


[[Image:draft_Content_671798508-1-s2.0-S0213131515000243-gr6.jpg|center|300px|Modelo ideal de Herschel-Bulkley y modelos regularizados de ...]]


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

Figura 6.

<span id='spar0080'></span>
Modelo ideal de Herschel-Bulkley y modelos regularizados de Herschel-Bulkley-Papanastasiou y Hershel-Bulckley-DV, ''n'' >1.                  

</span>
|}

<span id='par0400'></span>
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 <math display="inline">\dot{\gamma }=\dot{\gamma }_c</math> ):

<span id='eq0185'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>m\tau _y\dot{\gamma }_c=k\dot{\gamma }_c^n+\tau _y</math>
|}
| style="text-align: right;" | ( 37)
|}

<span id='sec0125'></span>
==8. Modelo discreto. Formulación de elementos finitos==

<span id='par0405'></span>
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.

<span id='par0410'></span>
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.

<span id='par0415'></span>
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 <math display="inline">\mu =\mu (\dot{\gamma })</math> . 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.      

<span id='par0420'></span>
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.      

<span id='par0425'></span>
Para la discretización espacial mediante el método de elementos finitos es necesario construir subespacios discretos '''V'''<sub>''h''</sub>  ⊂ '''V''' , ''y''  ''Q''<sub>''h''</sub>  ⊂ ''Q''  que aproximen los espacios continuos.      

<span id='par0430'></span>
Sean '''V'''<sub>''h''</sub>  ''y''  ''Q''<sub>''h''</sub>  los espacios de elementos finitos para interpolar las funciones vectoriales (velocidad) y escalares (presión), respectivamente, y sea Ω una partición de elementos finitos ''Ω''  = ∪ ''Ω''<sup>''e''</sup> , ''e''  = 1, …, ''n''<sub>''ele''</sub> , donde ''n''<sub>''ele''</sub>  es el número de elementos.      

<span id='par0435'></span>
En la formulación estándar de Galerkin se toman las funciones de ''test''  iguales a las funciones de forma, así que '''v'''<sub>''h''</sub>  ∈ '''V'''<sub>''h''</sub>  ''y''  ''q''<sub>''h''</sub>  ∈ ''Q''<sub>''h''</sub> .      

<span id='par0440'></span>
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 [[#bib0065|[13]]]  ofrece una descripción comparativa de varios de los métodos de estabilización propuestos en las últimas décadas.      

<span id='par0445'></span>
Los métodos de estabilización más usados en la actualidad están basados en los métodos de subescalas [[#bib0215|[43]]]  and [[#bib0220|[44]]]''.''  Estos métodos consisten en descomponer la solución, por ejemplo, la velocidad '''u'''  en 2 componentes <math display="inline">\mbox{u}=\mbox{u}_h+\bar{\mbox{u}}</math> ; una componente '''u'''<sub>''h''</sub> , resuelta en la escala de la malla de elementos finitos considerada y una subescala <math display="inline">\bar{\mbox{u}}</math> , 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.      

<span id='par0450'></span>
La solución <math display="inline">\bar{\mbox{u}}</math>  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 <math display="inline">\bar{\mbox{u}}=0</math>  en el contorno de los elementos.      

<span id='par0455'></span>
El residuo de la ecuación de momento en la escala grande, '''R'''<sub>''h''</sub> , resulta en:

<span id='eq0190'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 38)
|}

<span id='par0460'></span>
La subescala de velocidad, <math display="inline">\bar{\mbox{u}}</math> , se aproxima de distinta forma en cada método de estabilización.      

<span id='par0465'></span>
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.      

<span id='par0470'></span>
En el método de las subescalas algebraicas (ASGS), <math display="inline">\bar{\mbox{u}}</math>  se toma proporcional al residuo '''R'''<sub>''h''</sub>

<span id='eq0195'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\bar{\mbox{u}}=-\tau _1\mbox{R}_h</math>
|}
| style="text-align: right;" | ( 39)
|}

donde ''τ''<sub>1</sub>  es un parámetro numérico.      

<span id='par0475'></span>
En el método de las subescalas ortogonales (OSS) la subescala <math display="inline">\bar{\mbox{u}}</math>  se toma proporcional a la proyección ortogonal de dicho residuo '''R'''<sub>''h''</sub> :

<span id='eq0200'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\bar{\mbox{u}}=-\tau _1P_h^{\perp }\mbox{R}_h=-\tau _1(\mbox{R}_h-P_h\mbox{R}_h)</math>
|}
| style="text-align: right;" | ( 40)
|}

donde ''P''<sub>''h''</sub>  es la proyección sobre el espacio de los elementos finitos y <math display="inline">P_h^{\perp }=\mbox{I}-P_n</math>  es la proyección ortogonal.      

<span id='par0485'></span>
Comparando ambos métodos se observa que la diferencia radica en sustituir el término '''R'''<sub>''h''</sub>  de la versión ASGS por <math display="inline">P_h^{\perp }\mbox{R}_h</math>  en la versión OSS.      

<span id='sec0130'></span>
===8.1. Método ASGS===

<span id='par0490'></span>
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:

<span id='par0495'></span>
Hallar <math display="inline">\mbox{u}_h^{n+1}\mbox{ }y\mbox{ }p_h^{n+1}</math>  tales que

<span id='eq0205'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 41)
|}

<span id='eq0210'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}\displaystyle \int _{\Omega }[q_h\nabla .\mbox{u }_h^{n+1,i}]d\Omega +\sum _{e=1}^{n_el}\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}</math>
|}
| style="text-align: right;" | ( 42)
|}

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

<span id='sec0135'></span>
===8.2. Método OSS===

<span id='par0510'></span>
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 [[#bib0225|[45]]] :

<span id='eq0215'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>P_{h\rho }\left(\mbox{R}_h^{n+1}\right)=\rho P_h\left(\frac{\mbox{R}_h^{n+1}}{\rho }\right)</math>
|}
| style="text-align: right;" | ( 43)
|}

<span id='par0515'></span>
Sustituyendo <math display="inline">\mbox{R}_h^{n+1}</math>  y la proyección modificada para el método OSS queda que:      

<span id='par0520'></span>
El problema discretizado con linealización de Picard, integración temporal BDF1 y estabilización OSS consiste en hallar <math display="inline">\mbox{u}_h^{n+1}\mbox{ }y\mbox{ }p_h^{n+1}</math>  tales que:

<span id='eq0220'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{c}\int _{\Omega }[\frac{\rho }{\delta t}(\mbox{u }_h^{n+1}-\mbox{u }_h^n)\cdot \mbox{v }_hd\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}</math>
|}
| style="text-align: right;" | ( 44)
|}

<span id='eq0225'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 45)
|}

<span id='sec0140'></span>
===8.3. Método ''split-'' OSS                        ===

<span id='par0525'></span>
El método ''split'' -OSS  [[#bib0075|[15]]]  and [[#bib0225|[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:

<span id='eq0230'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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 }_hd\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}</math>
|}
| style="text-align: right;" | ( 46)
|}

<span id='eq0235'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 47)
|}

para las iteraciones ''i''  = 1,2… hasta la convergencia, es decir, hasta que '''u'''<sup>''n'' +1,''i'' −1                            </sup>  ≈ '''u'''<sup>''n'' +1,''i''</sup>  y <math display="inline">p_h^{n+1,i}\approx p_h^{n+1,i-1}</math>  en la norma elegida.      

<span id='sec0145'></span>
===8.4. Parámetros de estabilización===

<span id='par0530'></span>
El parámetro ''τ''<sub>1</sub>  de las ecuaciones (39) y (40) se elige con el fin de obtener esquemas numéricos estables y velocidades de convergencia óptimas (ver [[#bib0230|[46]]] ). Este parámetro se calcula para cada elemento ''Ω''<sub>''e''</sub> . Para ''τ''<sub>1</sub>  se toma:

<span id='eq0240'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\tau _1=\left[c_1\frac{\mu }{(h)^2}+c_2\frac{\rho \left|\mbox{u}^e\right|}{h}\right]^{-1}</math>
|}
| style="text-align: right;" | ( 48)
|}

donde ''h  ''  es la longitud del elemento <math display="inline">e\mbox{   }\mbox{y}\mbox{   }\left|\mbox{u}^e\right|</math>  es la norma de la velocidad en el elemento ''e'' , ''c''<sub>1</sub>  y ''c''<sub>2</sub>  son coeficientes a elegir, μ y ''ρ''  son la viscosidad dinámica y densidad del fluido, respectivamente.      

<span id='par0535'></span>
Se recomiendan los valores de ''c''<sub>1</sub>  = 1 y ''c''<sub>2</sub>  = 2 [[#bib0230|[46]]] .      

<span id='sec0150'></span>
==9. Formulación matricial del problema==

<span id='par0540'></span>
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.      

<span id='sec0155'></span>
===9.1. Método ASGS===

<span id='par0545'></span>
En la versión matricial para el método ASGS la proyección <math display="inline">P_h\left(\mbox{u}_h^{n+1}\cdot \nabla \mbox{u}_h^{n+1}+\frac{1}{\rho }\nabla p_h^{n+1}\right)</math>  se trata con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

<span id='eq0245'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 49)
|}

<span id='par0550'></span>
En lo que sigue se usa la notación compacta <math display="inline">(a,b)=\int_{\Omega }a\cdot b\mbox{ }\mbox{d}\Omega </math> . En esta notación, la proyección de la ecuación (49) es la solución de:

<span id='eq0250'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 50)
|}

para todo <math display="inline">\mbox{v}_h^\ast \in \mbox{V}_h^\ast </math> , donde <math display="inline">\mbox{V}_h^\ast </math>  es el espacio '''V'''<sub>''h''</sub>  ampliado con los vectores de funciones continuas asociados a los nodos del contorno.      

<span id='par0555'></span>
El sistema algebraico resultante es:

<span id='eq0255'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 51)
|}

<span id='eq0260'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{D}\mbox{U}^{n+1}+\mbox{S}_p\left(\tau _1\right)\mbox{ }\mbox{P}^{n+1}=0</math>
|}
| style="text-align: right;" | ( 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.      

<span id='par0560'></span>
La ecuación (51) corresponde a la ecuación (41). La ecuación (52) corresponde a la ecuación (42).

<span id='par0565'></span>
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''<sup>a</sup>  y la función de forma de los nodos ''b''  por ''N''<sup>b</sup> , entonces las matrices de las ecuaciones anteriores, las cuales son válidas para los métodos restantes, son:

<span id='eq0265'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{M}_{ij}^{ab}=(N^a,\rho N^b)\delta _{ij}</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0270'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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 ^sN^b)\delta _{ij}\end{array}</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0275'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{G}_i^{ab}=(N^a,\partial _iN^b)</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0280'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0285'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{D}_j^{ab}=(N^a,\partial _jN^b)</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0290'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{S}_p\left(\tau _1\right)^{ab}=\left(\tau _1\mbox{u}_h^{n+1}\cdot \nabla N^a,\nabla N^b\right)</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0295'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{F}_i^a=\left(N^a,\mbox{f}_i\right)</math>
|}
| style="text-align: right;" | ( 53)
|}

donde ''δ''<sub>''ij''</sub>  es la delta de Kronecker''.''

<span id='sec0160'></span>
===9.2. Método OSS===

<span id='par0570'></span>
En la versión matricial para el método OSS la proyección <math display="inline">P_k\left(\mbox{u}_k^{n+1}\cdot \nabla \mbox{u}_k^{n+1}+\frac{1}{\rho }\nabla p_h^{n+1}\right)</math>  también se trata con un ciclo iterativo, al igual que para la linealización del término convectivo. Se definen:

<span id='eq0300'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 54)
|}

<span id='par0575'></span>
En notación compacta, la proyección de la ecuación (54) es la solución de:

<span id='eq0305'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 55)
|}

para todo <math display="inline">\mbox{v}_h^\ast \in \mbox{V}_h^\ast </math> , donde <math display="inline">\mbox{V}_h^\ast </math>  es el espacio '''V'''<sub>''h''</sub>  ampliado con los vectores de funciones continuas asociados a los nodos del contorno.      

<span id='par0580'></span>
El sistema algebraico resultante es:

<span id='eq0310'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 56)
|}

<span id='eq0315'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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</math>
|}
| style="text-align: right;" | ( 57)
|}

<span id='eq0320'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{C}(\mbox{U}^{n+1})\mbox{U}^{n+1}+\mbox{G}_{\pi }\mbox{P}^{n+1}=0</math>
|}
| style="text-align: right;" | ( 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.      

<span id='par0585'></span>
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'' .      

<span id='par0590'></span>
Las matrices de las ecuaciones (56) a (58) no definidas anteriormente son:

<span id='eq0325'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{G}_{\pi i}^{ab}=(N^a,\partial _iN^b/\rho )</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0330'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0335'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{S}_z\left(\tau _1\right)_j^{ab}=(\tau _1\partial _jN^a,\rho N^b)</math>
|}
| style="text-align: right;" | 
|}

<span id='eq0340'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{C}(\mbox{U}^{n+1})_{ij}^{ab}=(N^a,\mbox{u}_h^{n+1}\cdot \nabla N^b)\mbox{ }\delta _{ij}</math>
|}
| style="text-align: right;" | ( 59)
|}

<span id='sec0165'></span>
===9.3. Método ''split'' -OSS                        ===

<span id='par0595'></span>
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 <math display="inline">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)</math>  también se tratan con un ciclo iterativo al igual que para la linealización del término convectivo. Se definen:

<span id='eq0345'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{y}_h^{n+1}=P_h\left(\mbox{u}_h^{n+1}\cdot \nabla \mbox{u}_h^{n+1}\right)</math>
|}
| style="text-align: right;" | ( 60)
|}

y

<span id='eq0350'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{z}_h^{n+1}=P_h\left(\frac{1}{\rho }(\nabla p_h^{n+1}-\mbox{f}^{n+1})\right)</math>
|}
| style="text-align: right;" | ( 61)
|}

<span id='par0600'></span>
En la notación compacta, las proyecciones de las ecuaciones (60) y (61) son la solución de:

<span id='eq0355'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 62)
|}

<span id='eq0360'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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)</math>
|}
| style="text-align: right;" | ( 63)
|}

para todo <math display="inline">\mbox{v}_h^\ast \in \mbox{V}_h^\ast </math> , donde <math display="inline">\mbox{V}_h^\ast </math>  es el espacio '''V'''<sub>''h''</sub>  ampliado con los vectores de funciones continuas asociados a los nodos del contorno.      

<span id='par0605'></span>
El sistema algebraico resultante es:

<span id='eq0365'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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}</math>
|}
| style="text-align: right;" | ( 64)
|}

<span id='eq0370'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\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</math>
|}
| style="text-align: right;" | ( 65)
|}

<span id='eq0375'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{M}_{\pi }\mbox{Y}^{n+1}-\mbox{C}(\mbox{U}^{n+1})\mbox{U}^{n+1}=0</math>
|}
| style="text-align: right;" | ( 66)
|}

<span id='eq0380'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{M}_{\pi }\mbox{Z}^{n+1}-\mbox{G}_{\pi }\mbox{P}^{n+1}=0</math>
|}
| style="text-align: right;" | ( 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.      

<span id='par0610'></span>
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.      

<span id='par0615'></span>
La matriz de las ecuaciones (64) a (67) no definida anteriormente es:

<span id='eq0385'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\mbox{M}_{\pi ij}^{ab}=(N^a,N^b)\delta _{ij}</math>
|}
| style="text-align: right;" | ( 68)
|}

<span id='sec0170'></span>
==10. Conclusiones==

<span id='par0620'></span>
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.

<span id='par0625'></span>
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.

<span id='par0630'></span>
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.

<span id='ack001'></span>
==Agradecimientos==

<span id='par0635'></span>
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.

<span id='bibl0005'></span>
==Bibliografía==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] E. Bingham; Fluidity and Plasticity; McGraw-Hill, New York (1922)</li>
<li><span id='bib0190'></span>
[[#bib0190|[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.</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] M. Bercovier, M. Engelman; A finite-element method for incompressible non-Newtonian flows; J. Comput. Phys., 36 (1980), pp. 313–326</li>
<li><span id='bib0020'></span>
[[#bib0020|[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.</li>
<li><span id='bib0025'></span>
[[#bib0025|[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</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] R. Von Mises; Mechanik der festen Korper im plastisch deformablen Zustand; Gottinger Nachr, math-phys Kl (1913), pp. 582–592</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] T. Papanastasiou; Flow of material with yield; J. Rheol., 36 (1987), pp. 389–407</li>
<li><span id='bib0040'></span>
[[#bib0040|[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</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] S. Abdali, E. Mitsoulis; Entry and exit flows of Bingham fluids; J. Rheol., 36 (2) (1992)</li>
<li><span id='bib0050'></span>
[[#bib0050|[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</li>
<li><span id='bib0055'></span>
[[#bib0055|[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</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] F. Brezzi, M. Fortin; Mixed and Hibrid Finite Element Methods; Springer Series in Computational Mathematics 15, Springer, New York (1991)</li>
<li><span id='bib0065'></span>
[[#bib0065|[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</li>
<li><span id='bib0070'></span>
[[#bib0070|[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</li>
<li><span id='bib0075'></span>
[[#bib0075|[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</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] Principe J. (2008). «Subgrid scale stabilizad finite elements for low speed flows». Ph. D thesis, Universidad Politécnica de Cataluña.</li>
<li><span id='bib0085'></span>
[[#bib0085|[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</li>
<li><span id='bib0090'></span>
[[#bib0090|[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.</li>
<li><span id='bib0095'></span>
[[#bib0095|[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.</li>
<li><span id='bib0100'></span>
[[#bib0100|[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</li>
<li><span id='bib0105'></span>
[[#bib0105|[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</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] M. Cervera, M. Chiumenti; Size effect and localization in J2 plasticity; Int. J. Solids Struct., 46 (2009), pp. 3301–3312</li>
<li><span id='bib0115'></span>
[[#bib0115|[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</li>
<li><span id='bib0120'></span>
[[#bib0120|[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</li>
<li><span id='bib0125'></span>
[[#bib0125|[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</li>
<li><span id='bib0130'></span>
[[#bib0130|[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</li>
<li><span id='bib0135'></span>
[[#bib0135|[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</li>
<li><span id='bib0140'></span>
[[#bib0140|[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</li>
<li><span id='bib0145'></span>
[[#bib0145|[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</li>
<li><span id='bib0150'></span>
[[#bib0150|[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</li>
<li><span id='bib0155'></span>
[[#bib0155|[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</li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] R. Bird, G. Dai, B. Yarusso; Rheology and flow of viscoplastic materials; Rev. Chem. Eng., 1 (1983), pp. 1–70</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] J. Oldroyd; A rational formulation of the equations of plastic flow for a Bingham solid; Proc. Camb. Philos. Soc., 43 (1947), p. 100</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] E. Reiner; Handbuch der Phisik; Springer, Berlin (1958)</li>
<li><span id='bib0175'></span>
[[#bib0175|[35]]] W. Prager; Introduction to Mechanics of Continua; Ginn and Co, Boston, MA (1961)</li>
<li><span id='bib0180'></span>
[[#bib0180|[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</li>
<li><span id='bib0185'></span>
[[#bib0185|[37]]] E. Mitsoulis, L.R. Huilgol; Entry flows of Bingham plastic in expansions; J. Non-Newtonian Fluid Mech., 122 (2003), pp. 45–54</li>
<li><span id='bib0010'></span>
[[#bib0010|[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</li>
<li><span id='bib0195'></span>
[[#bib0195|[39]]] R. Bird, R. Armstrong, O. Hassager; Dynamics of Polymeric Liquids: I. Fluid Mechanics;  (2 nd ed.)Jhon Wiley and Sons, New York (1987)</li>
<li><span id='bib0200'></span>
[[#bib0200|[40]]] D. Peric, S. Slijecpcevic; Computational modelling of viscoplastic fluids based on a stabilized finite element method; Engineering Computations, 18 (2001), pp. 577–591</li>
<li><span id='bib0205'></span>
[[#bib0205|[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</li>
<li><span id='bib0210'></span>
[[#bib0210|[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.</li>
<li><span id='bib0215'></span>
[[#bib0215|[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</li>
<li><span id='bib0220'></span>
[[#bib0220|[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</li>
<li><span id='bib0225'></span>
[[#bib0225|[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.</li>
<li><span id='bib0230'></span>
[[#bib0230|[46]]] Codina R. (2000c). «Stabilized finite element approximation of transient incompresible flows using orthogonal subscale». Cimne, No. 197.</li>
</ol>

<span id='footer-area-dummy'></span>