Source: http://studylib.es/doc/1261904/cd-6167.pdf
Timestamp: 2017-12-15 10:11:39+00:00

Document:
CD-6167.pdf
A Sergio González, por la oportunidad de trabajar en su proyecto, por la paciencia, el incentivo diario y la confianza entregada. Por el privilegio de ser parte de su
A Juan Carlos De los Reyes, por su motivación y permanente ayuda durante el
desarrollo de este proyecto. Por la oportunidad de trabajar en el Centro de Modelización Matemática.
A Pedro Merino y Luis Miguel Torres, quienes han contribuido enormemente en
A Miguel Yangari, por su ayuda con la búsqueda de textos bibliográficos.
A Enrique, Bibiana y Daniela. En especial, a mi hermano David y a quien se convirtió en el
2.1. Diferenciabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2. Operadores en coordenadas polares . . . . . . . . . . . . . . . . . . .
3. Estudio del flujo de un fluido de Bingham en una tubería
3.1. Leyes constitutivas de un fluido visco-plástico . . . . . . . . . . . . .
3.2. Flujo en la sección transversal de una tubería . . . . . . . . . . . . . .
3.3. Regularización . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4. Algoritmo de Newton semi suave
4.1. Discretización en coordenadas polares . . . . . . . . . . . . . . . . . .
4.2. Aproximación por diferencias finitas . . . . . . . . . . . . . . . . . . .
4.2.1. Operador Laplaciano . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2. Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3. Divergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3. Método de Newton semi suave (Semismooth Newton Method) . . .
4.4. Aproximación del paso de Newton . . . . . . . . . . . . . . . . . . . .
4.4.1. Resultados numéricos del problema de Bingham en coordenadas polares . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5. Método de Schwarz aplicado al algoritmo de Newton semi suave
5.1. Métodos de descomposición de dominio . . . . . . . . . . . . . . . . .
5.2. Método de Schwarz . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1. Forma multiplicativa del método de Schwarz . . . . . . . . . .
5.2.2. Forma aditiva del método de Schwarz . . . . . . . . . . . . . .
5.3. Método de Schwarz aplicado al problema de Bingham . . . . . . . .
5.4. Generalización a varios subdominios . . . . . . . . . . . . . . . . . . .
5.4.1. Detalles de implementación . . . . . . . . . . . . . . . . . . . .
5.4.2. Resultados numéricos de la paralelización del fluido de Bingham en varios subdominios . . . . . . . . . . . . . . . . . . . .
5.5. Resultados de convergencia . . . . . . . . . . . . . . . . . . . . . . . .
El estudio de los fluidos visco-plásticos ha cobrado gran importancia debido a la
gran cantidad de materiales que presentan este comportamiento y su presencia en
distintos procesos industriales.
En el campo de la mécanica de fluidos han surgido problemas en los que los algoritmos numéricos y las técnicas computacionales tradicionales resultan en un alto
costo computacional. Existe una verdadera necesidad y oportunidad en el análisis
numérico para plantear propuestas que permitan mejorar la simulación numérica
de problemas que involucren fluidos visco-plásticos. La industria, a través de sus
procesos de producción, ha generado y propuestos nuevos retos y problemas. (Ver
[24], [14] y [3]).
La investigación generada en la industria involucra simulaciones a gran escala.
Las técnicas de descomposición por multidominios constituyen una importante herramienta en la aproximación de estos problemas. Por este motivo, en este proyecto
de titulación se plantea el uso del método de descomposición de dominios de Schwarz para plantear un esquema paralelizado del algoritmo de Newton semi-suave
para la resolución del flujo de un fluido de Bingham.
En este trabajo de titulación revisaremos algunos conceptos y teoremas de mecánica de fluidos para la mejor comprensión del problema a tratar. Presentaremos las
leyes constitutivas de Bingham y el modelo matemático que describe el movimiento de este fluido. La geometría en la que se trabaja es una tubería cilíndrica. Por lo
tanto, para una precisa simulación numérica, se plantea el modelo en coordenadas
Revisaremos también el método de Newton semi-suave, (Semismooth Newton
Method) que se propone en [8] para simular el flujo laminar y estacionario de un
fluido de Bingham en la sección transversal de una tubería.
Este proyecto de titulación consiste en el planteamiento de un esquema parale1
lizado del algoritmo de Newton semi-suave para el flujo laminar y estacionario de
un fluido de Bingham. En concreto, se propone dividir el problema de Bingham en
una familia de subproblemas que se resuelven simultáneamente utilizando el algoritmo de Newton semi-suave. Para abordar el tema de la paralelización se propone
el método de descomposición de dominios de Schwarz [25]. Esta técnica permite
elaborar algoritmos paralelizados que se potencian gracias al uso de computadoras
con multiprocesadores.
Finalmente, presentaremos los resultados numéricos de la aplicación del método de Schwarz para el problema de Bingham. En estos experimentos numéricos
detallaremos el tiempo de ejecución del programa y compararemos el esquema paralelizado en el que se usa la técnica de descomposición de dominio con el esquema
sin paralelizar, es decir, si resolvemos el problema en un dominio general. Además,
presentaremos experimentos numéricos para analizar la convergencia del algoritmo.
2.1. Diferenciabilidad
D EFINICIÓN 2.1. (Espacios de Banach). Se dice espacio de Banach a ( X, k · k) un es-
pacio normado completo.
D EFINICIÓN 2.2. (Lipschitz continuidad). Sean X e Y dos espacios normados y F :
X → Y se dice Lipschitz continua si existe una constante positiva L tal que
k F( x1 ) − F( x2 )kY ≤ Lk x1 − x2 kX
para todo x1 , x2 ∈ X .
D EFINICIÓN 2.3. (Derivada de Gâteaux). Sean X un espacio normado, U un abierto
no vacío de X y F : U ⊂ X → R una función. Fijado un vector no nulo h ∈ X y
u0 ∈ U ; si
F(u0 + th) − F(u0 )
DF(u0 ) = lı́m
existe y es lineal, se dice que F es Gâteaux diferenciable en la dirección h y en el
punto u0 ∈ U .
D EFINICIÓN 2.4. (Derivada de Fréchet). Sean X y Y dos espacios normados, A ⊂ X
un conjunto abierto, F : A ⊂ X → V una función. Se dice que F es diferenciable en
u0 en el sentido de Fréchet, si existe una aplicación lineal y continua L : X → R tal
F ( u0 + h ) + F ( u0 ) − L ( h )
P ROPOSICIÓN 2.1. Sea F : U ⊆ X → R una función diferenciable en el sentido de
Fréchet en u0 . Entonces, F es diferenciable en el sentido de Gâteaux en ese punto.
Demostración. Ver [1].
D EFINICIÓN 2.5. (Derivada de Newton) Sean X y Y dos espacios de Banach. Adi-
cionalmente, sea D ⊂ X un conjunto abierto no vacío y f : D → Y una función
dada. f se denomina Newton diferenciable en D si existe una familia de funciones
G : D → L(X, Y ) tal que
k f ( x + h) − f ( x ) − G f ( x + h)hkY = 0,
h →0 k h k X
para todo x ∈ D. Llamamos a G f derivada de Newton de f en D
A continuación introduciremos el concepto de función semismooth. Esta discusión
es esencial para el análisis del método de Newton semisuave (Semismooth Newton
Method) que se tratará en los capítulos posteriores.
Sea D ⊂ R n un conjunto abierto no vacío y f : D ⊂ R n → R m . Se dice que f es
una función semismooth si
| f ( x + h) − f ( x ) − Vh| = o(khk),
para x ∈ D. Aquí, V sirve como una derivada de Newton en el sentido de la definición 2.5.
La noción de Newton diferenciabilidad implica un concepto más general de función semismooth.
Por lo tanto, la noción de función semismooth introducida por la derivada de
Newton implica las siguientes caracterizaciones:
P ROPOSICIÓN 2.2. Sea F : V → R m definida en el conjunto abierto V ⊂ R n . Luego,
para x ∈ V las siguientes afirmaciones son equivalentes:
• f es semismooth en x
• f es Lipschitz continua en una vecindad de x, f ′ ( x, .) existe, y
M∈∂ f ( x +s)
k Ms − ( x, s)k = o(ksk) cuando s → 0
k f ( x + s) − f ( x ) − Msk = o(ksk) as s → 0
Demostración. Ver [27] Prop. 2.7
E JEMPLO 1. Sea e ∈ R n un vector cuyas componentes son unos y g > 0 una cons-
tante. La función máx : R n −→ R n , ~y 7→ máx( g~e, ~y) es Newton diferenciable (ver
[17, Lem.3]). La derivada de Newton de esta función está dada por la matriz diagonal Gm definida de la siguiente manera:
( Gm ( −
y )ii ) =
si yi ≥ 0,
para 1 ≤ i ≤ n,
si yi < 0.
Tenemos que Gm ∈ L(R n ) y | Gm −
y | está acotado para todo y ∈ R n , luego, es posible
| máx(0, −
y + h ) − máx(0, −
y ) − Gm ( −
y + h ) h | = 0 si khk∞ < mı́n |yi | : yi 6= 0,
lo que implica que Gm es la derivada de Newton de −
y −→ máx(0, −
2.2. Operadores en coordenadas polares
En este trabajo se plantea un esquema paralelizado del algoritmo de Newton
semi-suave para el flujo laminar y estacionario de un fluido de Bingham en la sección transversal de una tubería circular. Por lo tanto, debido a la geometría escogida,
en el Capítulo 4 trataremos con distintos operadores que requieren ser transformados a coordenadas polares. En esta sección presentamos el desarrollo para transformar los operadores involucrados a coordenadas polares.
Sean Ω ⊂ R2 un conjunto abierto y u : Ω → R tal que u ∈ C2 (Ω). Si tenemos que
( x, y) 7→ u( x, y) y hacemos uso de la transformación clásica a coordenadas polares
x = r cos θ, y = r sin θ,
r ∈ [0, ∞),
θ ∈ [0, 2π ).
u( x, y) → u( x (r, θ ), y(r, θ )).
En los capítulos posteriores trabajaremos con los operadores: Laplaciano (∆), Gra5
diente (∇) y Divergencia (div) en coordenadas polares. A continuación mostramos
el desarrollo para transformar los operadores de coordenadas rectangulares a polares.
En coordenadas rectangulares tenemos que
∂2 u( x, y) ∂2 u( x, y)
Para transformar ∆ a coordenadas polares, debemos expresar las segundas deriva∂2 u ( x,y)
∂2 u ( x,y)
das parciales ∂x2 y ∂y2 en función de r y θ. En principio, como x e y son
funciones de r y θ, de (2.8) tenemos
∂x (r, θ )
∂y(r, θ )
∂2 x (r, θ )
∂ x (r, θ )
∂ y(r, θ )
∂2 y(r, θ )
= −r sen θ,
= −r cos θ,
= −r sin θ.
Ahora, derivamos parcialmente u( x (r, θ ), y(r, θ )) respecto de r y, del desarrollo
previo, obtenemos:
∂u( x (r, θ ), y(r, θ ))
La segunda derivada parcial de u( x (r, θ ), y(r, θ )) respecto de r corresponde a:
∂ ∂u ∂y ∂x
∂ ∂u ∂x ∂y ∂u ∂2 x ∂u ∂2 y
∂2 u ∂x
∂ u ∂y
∂x ∂y ∂r ∂r
∂y ∂x ∂r ∂r
∂x ∂r2
∂y ∂r2
∂x2 ∂r
∂y2 ∂r
cos θ sin θ + 2 sin2 θ.
= 2 cos2 θ + 2
Análogamente, calculamos la primera y segunda derivada parcial de u( x (r, θ ), y(r, θ ))
respecto de θ:
∂x ∂θ
∂y ∂θ
∂x ∂θ 2
∂y ∂θ 2
∂x2 ∂θ ∂x ∂y ∂θ ∂θ
∂y2 ∂θ ∂y ∂x ∂θ ∂θ
= 2 r2 sin2 θ − 2
r sin θ.
cos θ sin θ + 2 r2 cos2 θ − r cos θ −
Si reescribimos (2.12) y (2.13) tenemos:
∂2 u 2
r cos θ + r sin θ.
Para determinar las derivadas parciales respecto de r y θ, sumamos (2.14) y (2.15) y
obtenemos la siguiente igualdad
cos θ sin θ +
∂2 u 1 ∂u
Por la regla de candena, podemos escribir (2.16) de la siguiente manera
1 ∂ ∂u
(r ) + 2 2 .
Las igualdades (2.16) y (2.17) se conocen como las versiones no coservativa y conservativa, respectivamente, del operador Laplaciano.
Para los operadores Divergencia y Gradiente se procede de manera análoga y se
Forma no conservativa:
divu =
∂(u1 ) ∂u2
Forma conservativa:
1 ∂(ru1 ) 1 ∂u2
Estudio del flujo de un fluido de
Bingham en una tubería
Un fluido de Bingham es un fluido visco-plástico que se caracteriza por poseer
un umbral de plasticidad o yield stress: si la fuerza que se impone sobre el fluido
excede este umbral o punto crítico, entonces el material inicia su movimiento con
deformación continua. Por lo tanto, si la fuerza o stress es menor al yield stress el
material se comporta como un sólido rígido. El estudio de fluidos de tipo viscoplástico ha cobrado gran importancia debido a su uso en innumerables procesos
industriales tales como producción de comida (salsas y pastas), encementación de
pozos petrolíferos (ver [24]), diseño de sistemas de suspensión para automóviles,
lubricación, producción de cerámicas, entre otros (ver [3], [14]).
En este capítulo se analizan las leyes constitutivas de un fluido visco-plástico
que es tratado desde la perspectiva de medio continuo. La idea de medio continuo
es una abstracción, no se requiere estudiar la materia desde una perspectiva molecular. Se considera una porción de materia formada por entidades pequeñas pero, de
dimensión finita y que interactúan con sus semejantes de acuerdo a ciertas leyes. En
este contexto, un fluido visco-plástico obedece las tres leyes básicas de conservación
de la mécanica de un medio continuo y, además, leyes constitutivas particulares que
se mencionan en la siguiente sección.
3.1. Leyes constitutivas de un fluido visco-plástico
Un fluido visco-plástico obedece a las tres leyes de conservación que constituyen los principios básicos de la mécanica clásica de un medio continuo, éstas son:
conservación de masa, conservación de momentum y conservación de la energía.
Sea D una región en un espacio n-dimensional, n = 2, 3 ocupada por un fluido.
Sea x ∈ D un punto, consideremos que una partícula del fluido pasa por x en el
tiempo t, con x = ( x, y, z). Sea y(x, t) la velocidad de la partícula del fluido. Para
cada tiempo fijo y es un campo vectorial en D. El fluido tiene densidad de masa
ρ(x, t), es decir, la masa por unidad de volumen del material en la posición x en el
D EFINICIÓN 3.1. (Conservación de masa.) El principio de conservación de masa en
forma diferencial establece que:
+ div (ρy) = 0
esta ecuación es conocida como ecuación de continuidad. Ver [4]
D EFINICIÓN 3.2. (Conservación de momentum.) En un fluido en movimiento, se pro-
duce un gradiente de velocidad y una transferencia de momento por velocidad de
área y su forma diferencial es (ver [4]):
= −∇ p + ρ f
Tensor de deformación o strain tensor E . El tensor de deformación mide la deformación del material cuando está en movimiento. Si el tensor de deformación es
igual a cero, significa que le movimiento del material es sin deformación, es decir,
rígido. Por el contrario, si el movimiento del material es con deformación, la tasa de
deformación está dada por la magnitud de E . (Ver [8]).
Tensor de estrés σ. El tensor de estrés describe las fuerzas internas en un material elástico debido a su deformación. Consideremos una región V del fluido con
superficie S. Si asumimos que sobre un elemento de superficie dS, el fluido (visto
desde la perspectiva de un medio continuo) fuera de V transmite una fuerza sobre
el fluido contenido en V, esta fuerza se denomina fuerza de superficie y se denota
por σ n dS. Las componentes σ n representan la tracción transmitida desde fuera hacia el interior de V. Estos vectores de tracción son conocidos como los vectores de
estrés. Sea σij la j-ésima componente de σi y σi n la i-ésima componente de σ n . Del
principio de equilibrio local de las fuerzas de estrés (stress) podemos escribir:
σi n = σij −
donde σij son las componentes del tensor de segundo orden σTot denominado tensor
de estrés de Cauchy. Por lo tanto, podemos reescribir la ecuación (3.3) de la siguiente
σ n = σTot · −
Ecuación de Cauchy del movimiento. (Ver [4]). La ecuación de Cauchy relaciona la
definición de balance de momento en un volumen regular con el concepto de tensor
= ρf + ∇ · σ
Incomprensibilidad. (Ver [20]). Del principio de conservación de la energía sabemos
que la tasa de incremento de masa en D es igual a la tasa a la cual la masa cruza la
frontera ∂D.
ρdV = −
ρy · −
n dS,
donde y es el campo de velocidad del fluido y −
n es el vector normal a D. Integrando
por partes, obtenemos la siguiente expresión.
+ div(ρy) = 0.
Si ρ es constante en el tiempo, la ecuación de la continuidad (3.7) implica que
div y = 0.
3.2. Flujo en la sección transversal de una tubería
El problema del flujo de un fluido de Bingham estacionario y laminar en la sección transversal de una tubería debido a una diferencia de presión ha sido tratado
en trabajos como [23] y [20]. En este último trabajo, se analiza el flujo de un fluido visco-plástico desarrollando el tensor de estrés de Cauchy, σTot , y el tensor de
deformación E como se muestra a continuación.
σTot = − p · I + σ
σij = g q
+ 2µEij
donde g corresponde al umbral de plasticidad o yield stress y µ es la viscosidad del
fluido. Si E = 0, tenemos que el material se mueve como un sólido rígido pues,
no existe deformación. Sin embargo, cuando E = 0, el tensor de estrés de Cauchy
(3.9) está indeterminado. A pesar de esto, podemos obtener información de (3.10)
analizando si E 6= 0.
Si E 6= 0, obtenemos:
k σ k2 =
2 kσ k
g + 2µ
De (3.10) tenemos que:
Por otro lado, si kE k = 0, concluimos que 12 kσk < g. Para resumir, las leyes
constitutivas de un fluido de Bingham se pueden escribir de la siguiente manera
(ver [20, Ch.VI] y [8]):
 12 kσk < g
 12 kσk ≥ g
ssi E = 0
ssi E =
1 − √1
2 kσk
Es así que, si el estrés aplicado al material, cuya magnitud está dada por 12 kσk
no excede el umbral de estrés g, el fluido se mueve como sólido rígido. Si el estrés
aplicado es mayor a g, el material se comporta como un fluido.
Ahora, las ecuaciones que describen el movimiento de fluidos incompresibles
= Div σTot + f ,
div y = 0,
En este trabajo nos interesa el flujo estacionario de un fluido. Por lo tanto, podemos
reescribir las ecuaciones de la siguiente manera:
Div σTot + f = 0.
Finalmente, el sistema que modela el flujo de un fluido de Bingham en una geometría en R d sigue de la siguiente manera: hallar y : Ω ⊂ R d → R d y p : Ω ⊂ R d → R
σ Tot = − p · I + σ,
σ = 2µε(y ) + 2g ε(y) ,
kε(y)k
kσk ≤ g,
Div σ − ∇ p + f = 0,
si ε(y ) 6= 0,
si ε(y ) = 0,
Consideremos ahora una tubería cuyo rayo generador es paralelo al eje x3 en el
sistema de ejes x1 x2 x3 . Sea Ω ⊂ R2 el dominio que representa la sección transversal
de una tubería como se muestra a continuación.
րΩ
Estudiaremos el flujo entre x3 = 0 y x3 = L, donde L corresponde a una longitud
dada. En las secciones transversales x3 = 0 y x3 = L imponemos lo siguiente:
p( x3 )| x3 =0 = 0,
p( x3 )| x3 = L = −cL,
donde c es un escalar positivo que representa la caída de presión por longitud.
Estamos interesados en los campos de velocidades y y de presiones p que satisfacen: la condición de incompresibilidad (3.14), las leyes constitutivas de un fluido
de Bingham (3.13), la condición de adhesión dada por:
y = 0 sobre
y las condiciones (3.18). Si asumimos que el flujo es laminar el problema tiene solución única (ver [20]). Consideramos que el fluido se mueve debido a un efecto
de caída de presión en la tubería y no actúan fuerzas de cuerpo sobre el material,
por lo tanto, f = 0 y el campo de velocidad se reduce a (0, 0, y), De la condición de
incompresibilidad tenemos que el tensor de velocidades es:
Luego, la ecuación (3.14) toma la forma:
∂σ31 ∂σ32
p = −cx3 .
Por lo tanto, buscamos un campo escalar y = y( x1 , x2 ) en Ω tal que satisfaga:
= −c,
ε3i
σ3i = g
1 + 2µ ε 3i , para
(ε13 +ε223 ) 2
i = 1, 2 en
El sistema anterior corresponde a la forma diferencial de las ecuaciones del flujo del
fluido de Bingham. Si reemplazamos (3.19) en (3.22) tenemos lo siguiente:
2 ∂x 1
2 2 ! 12
= g ∂x 1 + µ
|∇y|
Procedemos de la misma manera con σ32 y obtenemos:
= g ∂x 2 + µ
De (3.23) y (3.24) podemos escribir la primera ecuación de (3.22) de la siguente manera:
-div(σ31 , σ32 ) − c
∂y 
= −div  g ∂x 1 + µ
, g ∂x 2 + µ
−c
∂x 1 |∇y|
= -div
∇ y + µ∇ y − c
= −µ∆y − div
∇y − c
Finalmente, para el flujo laminar y estacionario de un fluido de Bingham en la
sección transversal de una tubería, el tensor de velocidades E corresponde a ∇y y
tensor de estrés está dado por q. El modelo se escribe como se muestra a continuación.
−µ∆y − div q = f , en Ω,
q = g ∇y ,
si ∇y 6= 0
|q| ≤ g
si ∇y = 0
y = 0,
en Γ.
Para aproximar numéricamente el campo vectorial de la velocidad y del flujo, es
necesario determinar las regiones en las que el material se mueve como sólido rígido
(∇y = 0) y las regiones en las que el movimiento del material presenta deformación
continua (∇y 6= 0). Las zonas donde el material se comporta como un sólido y como
un fluido se denominan zonas inactivas y activas, respectivamente, y se definen de
A := { x ∈ Ω : |∇y| 6= 0}
I := Ω\A.
Determinar las zonas activas (∇y 6= 0) e inactivas (∇y = 0) no es posible de
forma directa en (3.26). En efecto, el tensor de estrés q, cuando ∇y = 0, no está
determinado. Este hecho sugiere un tratamiento especial para el término no diferenciable |∇y|. Una forma de analizar este problema es a través de un proceso de
regularización. En [8] se propone la regularización de Huber. Este proceso consiste
en una regularización local, es decir, se aproxima la función |∇y| por una función
continuamente diferenciable y se realiza un suavizamiento local en la vecindad donde existen problemas. Así, se supera el problema de la no diferenciabilidad de una
función al reemplazarla por una versión continuamente diferenciable de la misma.
Para ejemplificar lo mencionado anteriormente revisemos la regularización de
Huber para el valor absoluto. Consideremos la función diferenciable ψγ (w) definida
 g|w| − g2 , si |w| ≥ g
ψγ ( w ) : = γ
 | w |2 ,
para γ >> 1.
Figura 3.1: Regularización de Huber para |.|
En la figura 3.1 observamos que en x = 0 resolvemos el problema de la no
diferenciabilidad del valor absoluto con una aproximación diferenciable.
Revisemos el proceso de regularización local para la norma euclidiana k · k. Consideremos la siguiente función continuamente diferenciable ψ : R2 → R definida
 gkwk − g2 , si kwk ≥ g
 k w k2 ,
Esta función es una regularización local de la norma euclidiana. En la siguiente
figura observamos como ψ suaviza la curva de k · k. Si tomamos γ lo suficientemente
grande, ψ aproxima k · k.
Figura 3.2: Regularización de Huber de la norma euclidiana.
Como se mencionó previamente, nuestro objetivo consiste en usar una regularización local para determinar las zonas activas e inactivas y determinar el tensor de
estrés q cuando ∇y = 0. Podemos revisar el proceso de regularización de Huber
para el modelo de Bingham de la siguiente manera.
Si consideramos ∇y 6= 0, podemos verificar que
máx (0, |∇y|) q = |∇y| g
= g∇y.
Sin embargo, las zonas inactivas se mantienen indeterminadas, es decir, nuestra mayor dificultad recae cuando ∇y = 0. En este caso, para determinar q, consideremos
una versión regularizada de (3.31).
, |∇yγ | qγ = g∇y
donde γ >> 1 corresponde al parámetro de regularización.
Es así, que proponemos tomar un parámetro de regularización γ >> 1 de tal
forma que el tensor de estrés quede definido de la siguiente manera:
|∇yγ |
γ∇y ,
|∇yγ | ≥ γ ,
|∇yγ | < γ .
Esta técnica nos conduce a la obtención del tensor de estrés regularizado. Gracias a (3.33) podemos determinar las regiones activas e inactivas del material. De
esta manera tenemos una forma explícita de q para poder obtener la solución numérica del flujo. Podemos entender (3.33), desde la perspectiva de las zonas activas e
inactivas como:
c.t.p en
Aγ ,
Aγ = { x ∈ Ω : |∇yγ | ≥
Iγ = Ω\Aγ .
De (3.33) y (3.32) podemos reescribir el tensor de estrés regularizado como sigue
qγ = g
γ∇yγ
máx( g, γ|∇yγ |)
Luego, la ecuación previa es equivalente a:
máx( g, γ|∇yγ ( x )|)qγ ( x ) − gγ∇yγ ( x ) = 0
Una vez obtenido qγ , de (3.37), podemos escribir el sistema (3.26) en su versión
regularizada de la siguiente manera:
−µ∆yγ − div qγ = f ,
máx( g, γ|∇yγ ( x )|)qγ ( x ) − gγ∇yγ ( x ) = 0,
c.t.p en Ω,
en ∂Ω.
La metodología propuesta en esta sección, que describe el proceso y la técnica de
regularización local para obtener la versión regularizada (3.33) se presenta en [8].
Algoritmo de Newton semi suave
En este capítulo, plantearemos el modelo del flujo de Bingham en coordenadas
polares para obtener una simulación más precisa en la geometría deseada, es decir, en la sección transversal de una tubería cilíndrica. Para nuestros propósitos de
simulación, proponemos una discretización en diferencias finitas (en coordenadas
polares) para el dominio circular. Además, revisaremos el método de Newton semi
suave (Semismooth Newton method, SSN-method) para la resolución del problema de
Bingham en una tubería. El objetivo principal de este trabajo es desarrollar un esquema paralelizado del algoritmo SSN y su implementación. Por lo tanto, revisaremos
nociones y conceptos básicos sobre el método de Newton, así como el algoritmo
planteado en [8] para la resolución numérica del fluido de Bingham.
4.1. Discretización en coordenadas polares
Como vimos en capítulos anteriores, el problema regularizado del flujo de un
fluido de Bingham está dado por:
−µ∆yγ ( x ) − div qγ ( x ) = fγ ,
a.e. Ω,
El dominio en el que se desea resolver el sistema (4.1) corresponde a la sección
transversal de una tubería cilíndrica. Sin embargo, (4.1) está dado en el sistema ortonormal de ejes 0x1 x2 x3 (coordenadas rectangulares). En esta sección se propone
plantear el sistema de ecuaciones en coordenadas polares para así obtener la solución en la tubería cilíndrica. Además, para la discretización espacial, se propone un
esquema en diferencias finitas (para coordenas polares). Por lo tanto, es necesario
transformar el sistema de optimalidad a coordenadas polares, esto incluye los operadores y las funciones involucrados.
−µ∆yγ − div qγ = f,
máx( g, γ|∇yγ |)qγ − gγ∇yγ = 0,
r ∈ [0, ∞), θ ∈ [0, 2π );
los operadores que se requieren transformar son los siguientes:
(r ) + 2 2
r ∂r r2 ∂θ 2
forma conservativa,
forma no conservativa.
1 ∂(ru) 1 ∂v
r∂u 1 ∂v
Para verificar estas igualdades hacemos referencia al Capítulo 2 y a [15].
Una vez transformados los operadores, es necesario expresar la dependencia que
tienen el campo de velocidades y, el tensor de estrés q y f de r y θ, pues, con el
cambio de coordenadas, x e y son ahora funciones de r y θ.
De (4.3) podemos escribir:
yγ ( x, y) → yγ ( x (r, θ ), y(r, θ )) → uγ (r, θ )
qγ ( x, y) → qγ ( x (r, θ ), y(r, θ )) → wγ (r, θ ).
fγ ( x, y) → fγ ( x (r, θ ), y(r, θ )) → hγ (r, θ ).
Luego, el sistema de optimalidad tranformado es el siguiente:
∂u γ (r,θ )
−µ r ∂r (r ∂r ) +
1 ∂ u γ (r,θ )
1 ∂(rw)
+ 1r ∂w
∂θ = hγ ,
máx( g, γ|∇uγ (r, θ )|)wγ (r, θ ) − gγ∇uγ (r, θ ) = 0,
4.2. Aproximación por diferencias finitas
4.2.1. Operador Laplaciano
Consideremos el dominio Ω ⊂ R2 . Buscaremos una solución para el problema
(4.12) utilizando un esquema en diferencias finitas.
Consideremos el caso [0, R] × Ω = [0, 2π ).
En la Figura previa podemos ver una malla de nodos en el plano r − θ. Como
nuestro objetivo es resolver el sistema (4.12) en coordenadas polares, la configuración de los nodos cambia de la siguiente manera: los nodos ubicados en las filas de
la malla r − θ se disponen a lo largo de los rayos de la malla circular o, las columnas
de nodos de la malla r − θ se ubican en los anillos de la malla circular. Luego, los
nodos en el rayo θ = 0 coinciden con los nodos en el rayo θ = 2π. Este cambio de
coordenadas advierte una complicación: el origen se encuentra dentro del dominio
de discretización.
Como el origen está dentro del dominio de interés y 0 ≤ θ < 2π se tiene condiciones de frontera periódicas en la dirección de θ. Sin embargo, en la dirección radial
r el origen requiere un trato especial. En [22] se da un tratamiento numérico de las
coordenadas polares evitando la colocación de un nodo de la cuadrícula en el polo.
Para ello se utiliza una grilla escalonada de la siguiente manera:
ri = (i − )∆r,
i = 1, 2, · · · , m.
, rm = R.
Es así que se evita tener un nodo en r = 0. Luego,
θ = j∆θ,
j = 1, 2, · · · , n.
Para la discretización del operador Laplaciano se tienen dos esquemas en forma
conservativa y no conservativa. Matemáticamente, las derivadas originales (4.4) y
(4.5) son las mismas. Sin embargo, en su forma discreta existen diferencias.
Para un mayor entendimiento sobre la diferencia de los dos esquemas, en general, revisemos la siguiente derivada:
Al discretizarla obtenemos:
(ρu)i − (ρu)i −1
Ahora, en su forma no conservativa, la derivada corresponde a:
Usando la misma aproximación numérica tenemos:
u − ui − 1
ρ − ρi −1
+ u ≈ ρi i
+ ui i
Por lo tanto, la forma no conservativa es diferente a la conservativa y su nombre se debe a la siguiente justificación: si sumamos los términos correspondientes
a cada uno de los nodos sobre la grilla, sólo los términos de los nodos extremos
deben permanecer para balancear la ecuación. Tomemos una grilla unidimensional
de 4 puntos, por lo tanto, 0 ≤ i ≤ 3. Si sumamos todos los términos en su forma
conservativa tendríamos:
(ρu)1 − (ρu)0 (ρu)2 − (ρu)1 (ρu)3 − (ρu)2
es decir, se mantienen sólo los términos correspondientes a los nodos i = 0 e i = 3.
Los otros términos se cancelan. Por otro lado, al realizar la suma con los términos en
la forma no conservativa obtenemos
ρ − ρ2
+ u1 1
+ ρ2 2
+ ρ3 3
+ u3 3
En este caso, ningún término se cancela, es decir, cada vez que se añade un término
la suma aumenta. Desde el punto de vista de las ecuaciones de la conservación que
se mencionan en el Capítulo 3, la forma conservativa de los operadores balancea la
ecuación, es decir, el flujo que sale es igual al que ingresa, por lo tanto, al sumar los
términos, se cancelan todos excepto los correspondientes a los nodos en la frontera.
Una vez mencionada la diferencia surge la pregunta sobre qué esquema usar.
En [19, Ch.12] se menciona que para los problemas que involucren choques (shocks)
en la solución se prefiere la forma conservativa. Choque, en un fluido, corresponde
a cambios discontinuos en las propiedades del fluido así como: velocidad, presión
o temperatura. En [18] se presentan varios trabajos en los que para problemas de
flujos que involucran discontinuidades en sus dominios espaciales es mejor usar
un esquema numérico en su forma conservativa. Sin embargo, hay que tomar en
cuenta otras consideraciones. En algunas aplicaciones en ingeniería se prefiere el
esquema no conservativo, aún si se resuelven problemas con choques o shocks. En [5]
se trabajó con las ecuaciones de flujo transónico y se obtuvieron mejores resultados
con el esquema no conservativo.
Tomando en cuenta las consideraciones mencionadas previamente, en el desarrollo de este proyecto de titulación se llevaron a cabo varios experimentos numéricos. Los operadores Laplaciando y Divergencia admiten ambas formulaciones, conservativa y no conservativa. Por lo tanto, para estos operadores se probaron, ambos
esquemas y sus combinaciones. Sin embargo, la mejor aproximación numérica, al
compararla con la solución exacta, se alcanzó si se elige, para el operador Laplaciano, la formulación no conservativa, en i = 1 y la conservativa para los otros
nodos de la grilla. En el caso de la Divergencia, la mejor aproximación se alcanzó
usando el esquema no conservativo.
Se conoce que la solución exacta u para el flujo de un fluido de Bingham en la
sección transversal de una tubería está dada por (ver c.f. [13]) :
ū( x ) =
 √ 2 2 h q
R − x1 + x2
2 + x2 − 2g
2g h 
 R− f
si 0 ≤
x12 + x22 ≤ R
x12 + x22 ≤
La sección transversal de la tubería está dada por el disco definido por:
Ω := { x = ( x1 , x2 ) ∈ R2 : x12 + x22 < R2 , R > 0.}
El valor máximo de la solución exacta es: máx u = 1,6.
En la siguiente tabla mostramos el error de aproximación en la norma H01 usando
los diferentes esquemas.
Cuadro 4.1: Esquemas
ku − −
u k H1
3.9 ×10− 4
13.6503
máx −
Cuadro 4.2: Error de aproximación
De los experimentos numéricos mostrados previamente, vemos que los esquemas 3 y 4, cuando se usa la formulación no conservativa del Laplaciano, el error de
aproximación en la norma H01 es de aproximadamente 13,7. Por lo tanto, descartamos este esquema de discretización para el operador Laplaciano. Sin embargo, si
combinamos los esquemas de discretización, para el Laplaciano, de la siguiente manera: la formulación no conservativa, en i = 1 y la conservativa para los otros nodos
de la grilla, obtuvimos:
no conservativa en i = 1, conservativa para
los otros nodos de la grilla
Cuadro 4.3: Esquema para el Laplaciano
conservativa en i = 1, no conservativa para
Cuadro 4.4: Esquema para la Divergencia
4.97 ×10− 4
Cuadro 4.5: Error de aproximación
De estos experimentos encontramos que las mejores aproximaciones son las entregadas por los esquemas 5 − 7 y 5 − 8. Finalmente, el esquema 5 − 8 es el que mejor
aproxima la velocidad máxima, máx −
u respecto de máx u = 1,6.
En resumen, para el Laplaciano, la forma conservativa del esquema en diferencias finitas es la que se usa para los nodos que corresponden a los índices: i =
2, · · · , m − 1. Para i = 1, se usa el esquema que corresponde a la forma no conservativa. Por lo tanto, para la aproximación del Laplaciano ∆u se utiliza la estrella
de 5 puntos y el esquema final es el siguiente:
u0j − 2u1j + u2j
1 u1,j−1 − 2u1j + u1,j+1
1 u2j − u0j
r1 2∆r
(∆r2 )
(∆θ )2
para i = 1, j = 1, · · · , n,
1 ri − 21 ui −1,j − (ri − 12 + ri + 21 )ui,j + ri + 21 ui +1,j
1 ui,j−1 − 2ui,j + ui,j+1
para i = 2, · · · , m − 1, j = 1, · · · , n.
Los coeficientes u0,j corresponden a la aproximación de la solución en el nodo
fantasma r = 0, pues, con la grilla escalonada, evitamos tener un nodo en ese polo. Por lo tanto, estos coeficientes son igual a cero y, el esquema simplificado es el
−2u1j + u2j
1 u2j
(∆r )
r1 2∆r r1
para i = 1, j = 1, · · · , n.
Al fijar los coeficientes u0,j = 0 se asegura que la matriz resultante sea diagonal
Tomaremos n = m como pasos de discretización tanto para ∆R y ∆θ para que la
grilla escalonada tenga el mismo número de nodos en los rayos y los anillos de la
malla circular.
El operador laplaciano se aproxima por la siguiente matriz por bloques:
−µ∆ ≈ Ah
 L1 D2 R2
 . ...
Lk−1 Dk
 0 · · · · · · ...
0 Lm−2 Dm−1
con k = 1, · · · , m − 1. Usaremos el índice k para determinar la posición de las submatrices en Ah . El índice i se reserva para la discretización nodal.
D1 R1 0 · · · · · · 0
corresponde al esquema (4.23). Cada bloque tiene la siguiente forma:
D1 = 
−2 1 +
 0 .. · · ·
 . · · · ..
Es decir, R1 = 2I.
Para k = 2, · · · , m − 1 usamos el esquema (4.24)
(∆θ )
 −2 
Dk = 
−2 1 + 
Rk = 
(i − 21 )
4.2.2. Gradiente
El esquema de discretización para el gradiente (4.6), usando diferencias centradas, es el siguiente:
ui +1,j − ui −1,j
para i = 1, · · · , m − 1, j = 1, · · · , n.
1 ui,j+1 − ui,j−1
El gradiente está aproximado por la siguiente matriz por bloques:
∇h :=
∂1h
∈ R m× m
donde denotaremos a las aproximaciones de ∂1h por ∇h,1 y 1r ∂2h por ∇h,2 . De esta
0m × m
− I 0m × m
donde I corresponde a la matriz identidad y 0m×m corresponde a la matriz nula.
∆r∆θ
 −1 0
1 ··· 0 
H= 0
0 · · · −1 0
4.2.3. Divergencia
De igual manera que para el operador Laplaciano, para la Divergencia existen
dos esquemas, el conservativo y el no conservativo. Una vez más nos encontramos
con la pregunta acerca de qué esquema usar. Se realizaron varios experimentos numéricos usando ambos esquemas y sus combinaciones, numéricamente, los experimentos que mejor aproximaron la solución fueron aquellos realizados usando una
combinación de los dos esquemas. La aproximación por diferencias finitas para la
divergencia, div w, se realizó utilizando el esquema conservativo (4.7) en los índices
que corresponden a i = 1 y j = 1, · · · , n. Para los índices i = 2, · · · , m, j = 1, · · · , n
se usó el esquema no conservativo (4.8),
1 ri + 21 wi + 21 ,j − ri − 12 wi − 21 ,j ≈
para i = 1, j = 1, · · · , n
wij wi +1,j − wi −1,j
r∂w w ∂w
para i = 2, · · · , m, j = 1, · · · , n
1 wi,j+1 − wi,j−1
para i = 1, · · · , m, j = 1, · · · , n
El operador divergencia está aproximado por la siguiente matriz por bloques:
divw ≈ ((∂1h ), (∂2h )) = Qh
donde denotaremos a las aproximaciones de ∂1h por Qh,1 y ∂2h por Qh,2 .
Qh,1
 B C
 −C B
∆r 
0 E 0 ···
corresponde al esquema (4.28) , donde:
e I ∈ R m×m es la matriz identidad.
 (2i − 1)
De (4.6), (4.7) y (4.8) podemos concluir que Qh,2 corresponde a la matriz que aproxih,2
ma 1r ∂u
∂θ , es decir, ∇ .
Una vez que tenemos las versiones discretas de los operadores en coordenadas polares podemos escribir la versión discreta de la primera ecuación del sistema
Ah~u + Qh w
~ − ~h = 0, .
Para la segunda ecuación del sistema tenemos:
máx( g~e, γξ (∇h ~u)) ⋆ w
~ − gγ∇h w
~ = 0 para γ > 0,
Se denota por ⋆ al producto de Hadamard de vectores, i.e.
~v ⋆~z = (v1 z1 , · · · , vm zm )⊤
y ~e ∈ R2m corresponde al vector de componentes iguales a 1. Para calcular |∇uγ | en
términos de ~u, definimos la siguiente función ([8]).
ξ : R2m → R2m
(ξ ( p))k = (ξ ( p))k+m = |( pk , pk+m )⊤ | para p ∈ R2m , k = 1, · · · , m,
donde (ξ (∇h~y ))k representa la aproximación del valor de |∇uγ | en cada nodo de la
Como vimos anteriormente, la versión discreta del gradiente es:
∈ R2m×m .
Finalmente, obtenemos el sistema de optimalidad (4.12) discretizado (ver [8]):
h ~ − ~h = 0,
 A ~u + Q w
máx( g~e, γξ (∇h~u)) ⋆ w
for γ > 0.
En [8] se obtuvo el sistema de optimalidad para el flujo estacionario de un fluido de Bingham. Sin embargo, en este trabajo el dominio corresponde al cuadrado
unitario y, para la discretización espacial, se usa el método de elementos finitos.
4.3. Método de Newton semi suave (Semismooth Newton Method)
Sea F : R n −→ R n . Estamos interesados en un sistema de la forma:
F( x ) = 0.
Si suponemos que F es Fréchet diferenciable, el sistema F( x ) = 0 es un sistema
smooth (ver [16]). Por lo tanto, un método tipo Newton puede ser utilizado para
Sin embargo, si tenemos una función F : R n → R n , que no es diferenciable
en el sentido de Fréchet y estamos interesados en resolver un sistema de la forma
F( x ) = 0, es posible definir una noción más débil de diferenciabilidad de F de tal
manera que sea posible plantear un esquema iterativo tipo Newton. (Ver [6]).
En el Capítulo 2 revisamos la definición de Newton diferenciabilidad. Gracias a
esta noción más débil de diferenciabilidad podemos aplicar un método tipo Newton
para resolver sistemas de la forma (4.36). Es decir, cuando F no es diferenciable en
el sentido de Fréchet (por lo tanto, F es un sistema nonsmooth). (Ver [6] y [16]).
Podemos plantear (4.35) como un sistema de la forma (4.36). La función máx :
−→ R n , definida en R n , es Newton diferenciable en R n , (ver [17, Lem.3]). La
derivada de Newton de esta función está dada por la matriz diagonal Gm definida
(Gmax (−
En el Capítulo 2 se muestra que Gm es, en efecto, la derivada de Newton de la
función máx.
La noción de Newton diferenciabilidad permite plantear un esquema tipo Newton como el método de Newton semi-suave (Semismooth Newton Method).
T EOREMA 4.1. (ver [6]) Sea x una solución de (4.36), con F Newton diferenciable en
una vecindad abierta V que contiene a x. Si
kG( x )−1 kL(Z,X) ≤ C,
donde G es la derivada de Newton de F, para alguna constante C > 0 y para todo
x ∈ V , entonces la iteración de Newton
x k+1 = x k − G ( x k )−1 F ( x k )
converge superlinealmente a x.
Demostración. Ver [6] Th. 5.4
Algoritmo SSN
A continuación mostramos el algoritmo de Newton semi-suave Semismooth Newton method:
Algorithm 1 Algoritmo de Newton semi suave
Elegir un punto inicial x0 y fijar k = 0. Si f ( xk ) = 0, parar.
Calcular G ( xk )
Calcular sk tal que:
G ( x k ) s k = − f ( x k ).
Actualizar xk+1 = xk + sk , y k = k + 1.
Como se mencionó previamente, la función máx es Newton diferenciable si está
definida de R n en R n para cualquier n ∈ N. Por lo tanto, se puede obtener la derivada de Newton del sistema discreto (4.35) . Así, se propone el algoritmo tipo Newton
(método de Newton semi-suave) para su resolución.
4.4. Aproximación del paso de Newton
El problema a resolverse es el siguiente:
F(~u , w
~ ) :=
~ − ~h
~ − gγ∇h~u
De los resultados de c.f. [27] sabemos que el operador máx y la función ξ son semismooth de R m → R m y de los resultados de c.f [27] (Th. 2.10 y Th. 2.11) sabemos
que la composición de funciones semismooth es semismooth. Ahora, como estamos
trabajando con el sistema de optimalidad discretizado, tenemos la función:
R m −→ R m
~u 7−→ máx( g~e, ~u)
que es Newton diferenciable en R m (ver [17, Lem.3]), ~e ∈ R m es un vector cuyas
componentes son unos y g > 0 es una constante . La derivada de Newton para esta
función está definida por la matriz diagonal Gm ∈ R m×m :
(Gm (~u )ii ) :=
si ui > g,
si ui ≤ g.
Por lo tanto, la matriz Jacobiana de F es la siguiente R3m×3m (ver [8])
G (~u, w
γ(χA D (~
w ) Ph (∇h~u) − gI )∇h D (mh )
donde mh := max ( g~e, γξ (∇h~u)) ∈ R2n , y χA = D (th ) ∈ R2m×2m con
( t h )i : =
si ξ (∇h~u)i ≥ γ ,
Ph ∈ R2m×2m es la matriz Jacobiana de ξ, luego, tenemos que:
Ph (~p )
∂ξ i+l
∂p j+l
Consecuentemente, el paso de Newton de (4.40) es (ver [8]):
uk , w
= − F(~
~ k ),
wk ) Ph (∇h~uk ) − gI )∇h D (mhk )
− Ah u~k − Qh w
~ k + fh
− D (mhk )~uk + γg∇h~uk
con mhk := máx( g~e, γξ (∇h~uk )) ∈ R2m y χAk+1 = D (thk ) ∈ R2m×2m .
La obtención del paso de Newton que se muestra en esta sección corresponde al
trabajo realizado en [8] para la simulación numérica de un fluido de Bingham.
Algoritmo de Newton semi-suave para el modelo del flujo de un fluido de Bingham
A continuación presentamos el algoritmo de Newton semi-suave (SSN) para resolver el sistema (4.40).
Algorithm 2 Algoritmo de Newton semi suave (ver [8])
Inicializar (~
u0 , w
Fijar k = 0
Determinar χAk+1 .
Calcular G (~
~ k ).
Calcular F(~
uk, w
~ k ). δu~k
6: Calcularδ =
δw~ k
7: Actualizar u~
~k + δu y w~k+1 := w
~ k + δw
k+1 : = u
8: parar, o actualizar k := k + 1 e ir al paso 2.
4.4.1. Resultados numéricos del problema de Bingham en coordenadas polares
En este ejemplo se considera el flujo de un fluido de Bingham en la sección transversal de una tubería dada por el disco definido por:
Se conoce que la solución exacta para este problema está dada por (ver c.f. [13]) :
Para resolver el sistema previo se utilizó una malla circular como la que se muestra
Figura 4.1: Malla circular
se consideraron los siguientes parámetros: ∆r = 0,0050 y ∆θ = 0,0314.
A continuación presentamos la solución exacta y la solución dada por el algoritmo SSN con los siguientes parámetros: R = 1, µ = 1, g = 1, ∆r = 0,0050, ∆θ =
0,0314, f = 10.
Figura 4.2: Solución exacta del flujo de un fluido de Bingham en la sección transversal de una tubería, max u=1.6000
Figura 4.3: Solución del flujo de un fluido de Bingham en la sección transversal de
una tubería con el algoritmo SSN, max u=1.5962
4.9718e-04
2.7616e-04
1.7531e-04
Cuadro 4.6: Comparación del error
Método de Schwarz aplicado al
En este capítulo introduciremos la teoría sobre los métodos de descomposición
de dominio. En particular, nos centraremos en el método de Schwarz. Esta técnica,
en combinación con el método de Newton Semi-suave, será utilizada para resolver
 Ah~u + Qh w
~ − ~h = 0,
máx( g~e, γξ (∇h~u)) ⋆ w
5.1. Métodos de descomposición de dominio
Las técnicas de decomposición de dominio para ecuaciones diferenciales parciales se basan en una idea general, que consiste en reformular la ecuación diferencial
parcial, o una aproximación de la misma, en problemas acoplados más pequeños sobre subdomnios que forman una partición del dominio original. Dado un dominio
computacional Ω, se lo divide en subdominios Ωl , l = 1, · · · , K, que pueden o no estar solapados, es decir, cada subdominio puede o no cubrir parte de los subdominios
vecinos. Con la reformulación del problema se obtiene una familia de subproblemas
de tamaño reducido que están acoplados unos a otros gracias a los valores que toma
la solución en las fronteras de los subdominios. Por lo tanto, en lugar de resolver el
problema original se resuelven pequeños problemas sobre los nuevos subdominios
un determinado número de veces.
En esta sección introduciremos los fundamentos matemáticos de los métodos
de descomposición de dominio. En general, se particiona el dominio computacional en subdominios y se introduce una ecuación para la interfaz entre subdominios.
El problema de la interfaz se define en términos del operador de Steklov-Poincaré.
En [25, Ch.1 Sec. 1.1] se encuentran los detalles y las referencias para un completo
desarrollo matemático de este operador. Para el entendimiento de los métodos de
descomposición de dominio y, en especial, para la interfaz introducida entre subdominios, exhibimos el concepto de transmisión de condiciones en las fronteras de los
subdominios (Ver [25]). Nos valdremos del problema de Poisson para mostrarlo:
−∆u = f en Ω,
u = 0 sobre ∂Ω.
Aquí, Ω es un dominio d-dimensional (d = 2, 3) y con frontera ∂Ω Lipschitz continua y f es una función de L2 (Ω). Asumamos que Ω es particionado en dos subdominios no solapados Ω1 y Ω2 , donde la frontera entre subdominios Γ := Ω1 ∩ Ω2 .
El problema de Poisson (5.2) puede ser reformulado para la partición de Ω en dos
subdominios de la siguiente manera:
−∆u1 = f en Ω1 ,
u1 = 0 sobre ∂Ω1 ∩ Ω2 ,
u1 = u2 sobre Γ,
= 1 sobre Γ,
u2 = 0 sobre ∂Ω2 ∩ ∂Ω,
−∆u = f en Ω .
Las igualdades 3 y 4 del sistema (5.3) corresponden a las ecuaciones de transmisión de condiciones para u1 y u2 sobre la frontera Γ. El problema (5.3) tiene sentido
mientras la solución del problema (5.2) sea C(Ω ).
Consideremos los dos problemas de Dirichlet siguientes para escribir el problema (5.3).
Sea λ el valor desconocido de u en la frontera Γ y
−∆wi = f
 wi = λ
en Ωi ,
sobre ∂Ωi ∩ ∂Ω
para i = 1, 2. Podemos formular la solución de (5.4) de la siguiente manera
wi = u0i + ui∗
donde u0i y ui∗ son las soluciones de los problemas:
−∆ui = 0
 u0
= λ sobre Γ,
−∆ui = f
ui∗ = 0 sobre ∂Ωi ∩ ∂Ω
u∗ = 0 sobre Γ.
Para cada i = 1, 2. Para que los problemas (5.3) y (5.4) sean iguales se debe cumplir
wi = ui para i = 1, 2. si y solo si
Por lo tanto, para que el problema original se pueda particionar en subproblemas,
se debe cumplir la condición (5.8). Para su demostración se introduce la ecuación de
interfaz de Steklov-Poincaré. Para más detalles ver [25, Ch.1 Sec. 1.1].
5.2. Método de Schwarz
Uno de los primeros métodos sobre descomposición de dominios para la resolución de ecuaciones diferenciales parciales es el método alternante de Schwarz que
fue introducido en 1869. En su forma clásica, el método consiste en resolver problemas elípticos con valores en la frontera en un dominio que resulta de la unión de
dos dominios que no necesariamente tienen la misma forma. Se resuelve el mismo
problema alternativamente en cada subdominio (ver [2]).
Consideremos el problema de Poisson (5.2) y el dominio original propuesto por
←− Γ2
−→ Γ1
Dominio original propuesto por Schwarz.
El método de Schwarz se inicializa con una primera solución u02 para el dominio
Ω2 . De esta solución se requieren los valores en la frontera Γ1 . Luego, iterativamente,
para k = 1, 2, 3 · · · , se resuelve el siguiente problema:
−∆u1 = f en Ω1 ,
u1k+1 = u2k sobre Γ1 ,
uk+1 = 0 sobre ∂Ω1 \ Γ1 .
Una vez que se obtiene la solución u1k+1 se resuelve:
−∆u2 = f en Ω2 ,
u2k+1 = u1k+1 sobre Γ2 ,
uk+1 = 0 sobre ∂Ω2 \ Γ2 .
El método alternante de Schwarz se usa en general en dominios que, al ser discretizados, sus mallas no coinciden en la región solapada como sucede en el dominio
planteado por Schwarz y cuya discretización se muestra en la siguiente figura:
Figura 5.1: Dominio de Schwarz discretizado [2]
Consideremos una discretización por diferencias finitas para los problemas (5.9)
y (5.10).
 A1 u1 = f 1 en Ω1 ,
= u2k sobre Γ1 ,
uΓk+
 uk+1
= 0 sobre ∂Ω1 \ Γ1 ,
∂Ω1 \Γ1
A uk+1 = f en Ω2 ,
= u1k+1 sobre Γ2 ,
= 0 sobre ∂Ω2 \ Γ2 .
∂Ω2 \Γ2
Los coeficientes ui asociados a la versión discreta de la solución se los denota de
uΩi
ui =  u∂Ωi \Γ 
uΓ i
Notemos que los coeficientes u∂Ωi \Γ son conocidos, pues, están dados por las
condiciones de Dirichlet homogéneas. La forma discreta del operador −∆ está dada
por la matriz Ai restringida a cada subdominio i = 1, 2. Esta matriz está compuesta
por tres submatrices:
Ai = AΩi A∂Ωi \Γ AΓi .
Ai es una matriz cuadrada. La submatriz AΩi corresponde a la matriz de los nodos al interior del dominio i. A∂Ωi \Γ representa el acoplamiento entre los nodos al
interior del subdominio i y los nodos en la frontera ∂Ωi . AΓi corresponde al acoplamiento de los nodos al interior del dominio y los nodos en la frontera artificial
Γi . Si el dominio corresponde al cuadrado unitario y se usa una aproximación por
diferencias finitas centradas, la matriz Ai posee una columna por cada nodo en el
dominio. Si un nodo se encuentra al interior de Ωi , su aproximación por diferencias
finitas está en la matriz Ai . Si el nodo pertenece a la interfaz, su aproximación está
en la submatriz AΓi . Este esquema de aproximación por diferencias finitas está dado
por el stencil o estrella de 5 puntas.
Usando esta descomposición podemos reescribir los sistemas (5.11) y (5.12) de la
siguiente manera (ver [2]):
 Ω1
 +1
AΩ1 A∂Ω1 \Γ AΓ1  uk∂Ω
1 \Γ
uΓ1
 k+1
∂Ω2 \Γ
 ∂Ω2 \Γ
∂Ω \Γ
 = f 1 en Ω1 ,
 = f 2 en Ω2 ,
El método alternante de Schwarz propone interpolar los valores al interior del
dominio Ω j (j = 1, 2) para hallar los valores en la frontera Γi (i = 1, 2), es decir, para
determinar los coeficiente uΓi (i = 1, 2). La interpolación es un proceso para estimar
valores que se encuentran entre puntos conocidos. En este caso, los puntos conocidos corresponden a los valores resultantes (en los nodos del primer subdominio) al
término de la solución del primer subproblema (5.11). Los valores obtenidos de la
interpolación constituyen una solución inicial para los nodos en la frontera artificial.
En esta sección introduciremos un operador de interpolación lineal, que constituye
el caso más general de interpolación. Cuando se tienen dos puntos conocidos, éstos
se pueden unir mediante una línea recta. Es decir, se utiliza la ecuación de la recta
para estimar diferentes puntos entre los puntos conocidos. Con el proceso de interpolación se logra aproximar el valor de una función en un valor desconocido de su
dominio. Para la interpolación lineal se utiliza un polinomio de primer grado.
Sea IΩ j −→Γi el operador discreto que interpola los valores de los nodos desde
interior de Ω j en la frontera Γi . El algoritmo de Schwarz se escribe de la siguiente
Algorithm 3 Algoritmo alternante de Schwarz
Elegir w01 ←− 0.
for k = 0, · · · do
Resolver para u1k+1 :
A1 u1k+1 = f 1 en Ω1 ,
= w1k sobre Γ1 ,
= 0 sobre ∂Ω1 \ Γ1 .
u∂Ω
w2k+1 ←− IΩ1 −→Γ2 uΩ1 .
Resolver para u2k+1
A2 u2k+1 = f 2 en Ω2 ,
= w2k+1 sobre Γ2 ,
w1k+1 ←− IΩ2 −→Γ1 uΩ2 .
if kw1k+1 − w1k k ≤ tolΓ1 y kw2k+1 − w2k k ≤ tolΓ2 then
if ku1k+1 − u1k k ≤ tolΓ1 y ku2k+1 − u2k k ≤ tolΓ2 then
Una vez que el algoritmo termina, el vector u1k+1 contiene la aproximación discreta de la solución para Ω1 . A su vez, u2k+1 contiene la solución para el subdominio
Ω2 . Para la región solapada se pueden usar cualquiera de las dos soluciones, pues,
ambas convergen a la misma solución en la región solapada si la malla es lo suficientemente fina.
En [2] se presenta la convergencia numérica del método alternante de Schwarz.
Se resuelven los subproblemas (5.11) y (5.12) usando los métodos de Gauss-Seidel
y Jacobi. Además, se muestra el comportamiento del algoritmo cuando la malla es
refinada pero, la zona solapada se mantiene constante.
5.2.1. Forma multiplicativa del método de Schwarz
En la sección previa mostramos el método alternante tal como Schwarz lo introdujo por primera vez en 1869 y en el que se considera el caso general, cuando,
al discretizar los dos subdominios de Ω, la discretización no coincide en la región
solapada. (Ver Figura 5.1).
En general, el método de Schwarz consiste en descomponer el domino Ω en dos
subdominios solapados Ω1 y Ω2 tales que Ω = Ω1 ∪ Ω2 como se muestra en la
Figura 5.2: Dominio solapado.
donde Γ1 = ∂Ω1 ∩ Ω2 , Γ2 = ∂Ω2 ∩ Ω1 y Ω1,2 := Ω1 ∩ Ω2 .
El método de Schwarz, en su forma multiplicativa para el problema de Poisson,
consiste en resolver sucesivamente los siguientes problemas:
Sea u0 una función de inicialización definida en Ω. Para facilitar el estudio del método multiplicativo podemos fijar en cero los valores de u0 en ∂Ω. Sea û02 := u0 |Ω2 ,
para k ≥ 0, índice de iteraciones del método de Schwarz, definimos dos sucesiones
û1k+1 y û2k+1 resolviendo respectivamente:
−∆û1 = f en Ω1 ,
û1k+1 = û2k sobre Γ1 ,
ûk+1 = 0 sobre ∂Ω1 ∩ ∂Ω
−∆û2k+1 = f en Ω2 ,
ûk+1 = ûk+1 sobre Γ ,
û2 = 0 sobre ∂Ω2 ∩ ∂Ω.
En la región solapada, Ω1,2 , los valores de los nodos están dados tanto por la solución de û1k+1 como de û2k+1 . Al término de la resolución del método, ambas soluciones convergen a los mismos valores.
A continuación, presentamos el algoritmo del método de Schwarz en su forma
Algorithm 4 Algoritmo multiplicativo de Schwarz
Resolver (5.2) en Ω0
û02 ←− IΩ0 −→Γ2 uΩ0 .
Elegir w01 ←− û02
A1 û1k+1 = f 1 en Ω1 ,
ûkΓ+
ûk∂Ω
w2k+1 ←− û1k+1
A2 û2k+1 = f 2 en Ω2 ,
û∂Ω
w1k+1 ←− û2k+1
if kû1k+1 − û1k k ≤ tolΓ1 y kû2k+1 − û2k k ≤ tolΓ2 then
En este contexto, el operador discreto de interpolación, interpola los valores obtenidos de una solución inicial, en una malla gruesa, y nos entrega los valores en
los nodos de frontera de la malla refinada, donde vamos a resolver el problema. A
diferencia del caso alternante, en donde el proceso de interpolación se lleva a cabo
una vez resuelto el problema en el primer subdominio y se comienza con valores en
la frontera iguales a cero.
5.2.2. Forma aditiva del método de Schwarz
En la sección anterior revisamos una discretización general de Ω. Sin embargo,
si el dominio de interés es particionado en dos subdominios de tal manera que la
malla en la zona solapada coincide tanto en el primer como en el segundo dominio,
entonces, se puede optimizar el algoritmo de Schwarz para aprovechar la información en los nodos del dominio compartido y resolver los subproblemas (5.13) y (5.14)
paralelamente. Es decir, si tomamos una solución inicial del problema restringida a
los dos subdominios, esta define dos sucesiones para cada subproblema. Como los
nodos coinciden en la zona solapada, se toman los valores de las sucesiones en los
nodos que corresponden a las fronteras Γ1 y Γ2 para resolver cada subproblema al
Si tomamos Û10 := u0 |Ω1 y Û20 := u0 |Ω2 se pueden resolver los dos subproblemas
independientemente de la siguiente manera:
−∆Û1 = f en Ω1
Û1k+1 = Û2k sobre Γ1
Û k+1 = 0 sobre ∂Ω1 ∩ ∂Ω
−∆Û2 = f en Ω2
Û2k+1 = Û1k sobre Γ2
La aproximación (5.15) y (5.16) corresponde a la forma aditiva del método de
Schwarz. Estas dos nuevas formulaciones fueron planteadas en [9] y [21].
Una vez que el algoritmo termina, Û1k+1 contiene la aproximación discreta de la
solución para Ω1 y Û2k+1 contiene la solución para el subdominio Ω2 . En este caso se
obtienen las soluciones en los nodos de Ω1,2 en ambos vectores Û1k+1 y Û2k+1 , pues,
se resuelven los dos subproblemas de manera simultánea. Ambas soluciones convergen a los mismos valores en los nodos de Ω1,2 . Para determinar la solución final,
una alternativa consiste particionar Ω1,2 por la mitad. La solución de los nodos pertenecientes a la primera mitad de Ω1,2 estará dada por los componentes del vector
Û1k+1 en aquellos nodos. El valor de los nodos en la segunda mitad de Ω1,2 estará
dado por la solución de Û2k+1 en aquellos puntos.
A continuación presentamos el algoritmo de Schwarz en su forma aditiva.
Algorithm 5 Algoritmo aditivo de Schwarz
Resolver en (5.2) en Ω0 ,
Û20 ←− IΩ0 −→Γ2 uΩ0 .
Û10 ←− IΩ0 −→Γ1 uΩ0 .
Elegir w0i ←− Ûi0 para i = 1, 2
Resolver para Uik+1 :
Ai Ûik+1 = f i en Ωi ,
ÛΓk+
= wik sobre Γi ,
= 0 sobre ∂Ωi \ Γi .
wik+1 ←− Ûik+1
if kwik+1 − wik k ≤ tolΓi then
if kÛik+1 − Ûik k ≤ tolΓi then
5.3. Método de Schwarz aplicado al problema de Bingham
Los métodos de descomposición de dominio se introdujeron como técnicas para
la resolución de ecuaciones diferenciales parciales en base a una descomposición del
dominio espacial del problema en varios subdominios [21]. En las secciones previas
hemos introducido el método de Schwarz para un operador lineal y elíptico a través del problema de Poisson. Sin embargo, nuestro objetivo es aplicar esta técnica
al problema de Bingham que se discutió en el Capítulo 3, que constituye un sistema
quasilineal. Se han propuesto distintos algoritmos de descomposición de dominio
para la resolución de problemas quasilineales. En [26] se considera una función en
un espacio de Banach reflexivo que se supone que es Gateaux diferenciable y se demuestra la convergencia para un método de descomposición de dominio aplicado a
este problema. Nos valdremos de los resultados presentados en [11] y [10], en donde se adapta la teoría de descomposición de dominio a funcionales non-smooth y con
restricciones de variación total para extender la aplicación del método de Schwarz
al problema de Bingham.
El objetivo principal de este trabajo es desarrollar un esquema paralelizado para
la resolución del problema de Bingham. Si analizamos la familia de subproblemas de
la forma aditiva y multiplicativa del método podemos concluir que el método multiplicativo no se presta para la resolución del problema en paralelo. Supongamos
que estamos trabajando en un dominio como el de la Figura 5.2, en cada iteración
del método se requiere resolver el problema en uno de los subdominios y utilizar
esta solución para fijar las condiciones en la frontera artificial del segundo subdominio. Por lo tanto, se resuelve el problema secuencialmente y no paralelamente. Es
por esta razón que hemos elegido el método de Schwarz en su forma aditiva para
construir un esquema paralelizado que nos permita utilizar el algoritmo de Newton
semi suave y resolver el problema de Bingham.
Recordemos que el problema (4.35) está planteado en coordenadas polares de tal
manera que podamos resolverlo en la sección transversal de una tubería cilíndrica.
En este caso, si particionamos el dominio en dos subdominios solapados como se
muestra en la Figura 5.2, vemos que la malla en la zona solapada coincide para los
dos subdominios.
Figura 5.2: Zona solapada Ω1,2 .
Por lo tanto, para nuestros propósitos, la forma alternante del método no es aplicable pues, resultaría innecesario interpolar los valores de cada subdominio en cada
Mostraremos la aplicación del método resolviendo el problema en dos subdominios (Figura (5.4) ).
Figura 5.3: dominio original
Figura 5.4: subdominios Ω1 y Ω2 .
Sin embargo, como hemos visto en las secciones previas, el problema se resuelve
en dos subdominios solapados como se muestra en la Figura 5.5
Figura 5.5: Subdominios solapados
Sea N la cantidad de nodos en el dominio completo (Figura 5.5) y K = 2 el número de subdominios en los que se particiona Ω. Para cada subdominio Ωl , con
l = 1, 2, denotaremos por M a la cantidad de nodos por subdominio solapado, es
decir, cada subdominio, sin ser solpadado, tendría N
K nodos. En nuestro caso vamos
a aplicar el método de descomposición de dominio con solapamiento. Por lo tanto,
los subdominios Ω1 y Ω2 comparten una determinada cantidad de nodos que se
encuentran en la región Ω1,2 . Es así que M = N
K + t, donde t es el número de nodos que el primer subdominio comparte con el segundo subdominio. Finalmente, el
dominio compartido Ω1,2 tiene 2t nodos.
En la aplicación que vamos a desarrollar a continuación la formulación del problema corresponde a:
−µ∆Û1 − div W1 = h
máx( g, γ|∇Û1k+1 |)W1k+1 − gγ ∇Û1k+1 = 0
Û k+1 = Û k
−µ∆Û2k+1 − div W2k+1 = h
máx( g, γ|∇Û k+1 |)W k+1 − g ∇Û k+1 = 0
Û2 = Û1
Û2 = 0
en Ω1 ,
c.t.p Ω1 ,
en Γ1 ,
en Ω2 ,
c.t.p Ω2 ,
en Γ2 ,
en ∂Ω2 \ Γ2 .
El problema (5.17), que corresponde al primer subdominio, sólo posee condiciones en la interfaz Γ1 , a diferencia del segundo problema (5.18) que posee condiciones
en la frontera exterior ∂Ω y en la frontera Γ2 .
Como vimos en el capítulo anterior, la versión discreta del sistema que vamos a
resolver corresponde a (4.35). Por lo tanto, si aplicamos la forma aditiva del método
de Schwarz al sistema discreto obtendríamos los siguientes subproblemas:
→ k+1
+ Q1h W1 k+1 = h 1k+1
 A1 U1
máx( g e , γξ (∇1h U1 k+1 )) ⋆ W1 k+1 − gγ∇1h U1 k+1 = 0 paraγ > 0
→ k+1 −
= U2 k
A2h U2 k+1 + Q2h W2 k+1 = h 2k+1
máx( g e , γξ (∇h U k+1 )) ⋆ W k+1 − gγ∇h U k+1 = 0 paraγ > 0
U2 k+1 = U1 k
Con el fin de determinar las condiciones de frontera en Γ1 y Γ2 , es necesario
interpolar estos valores desde un dominio inicial Ω0 (Figura 5.6) y a partir de una
solución inicial uΩh0 (Figura 5.7).
Figura 5.6: Dominio inicial Ω0 .
Figura 5.7: Solución inicial uΩh0 en el dominio Ω0 .
Ω el operador discreto que interpola los valores desde los nodos
Sea IΩ0 →Ω := IΩ
en el interior del dominio Ω0 a los nodos en el dominio Ω. En el proceso de interpolación se estiman los valores que se encuentran entre los puntos conocidos de la
malla del dominio Ω0 . El valor interpolado en cada nodo de la malla más refinada,
es decir Ω, proviene de una interpolación lineal de los valores en los nodos vecinos
de la grilla en cada dirección.
La forma más sencilla de interpolación es la lineal, este método consiste en unir
dos puntos adyacentes conocidos mediante un polinomio de primer grado. Con la
ecuación de la recta que pasa por los dos puntos conocidos se puede calcular el
valor de y para cualquier x. Como nuestro ejemplo corresponde a una malla en dos
dimensiones, hemos realizado una interpolación bilineal, la interpolación bilineal es
una extensión de la interpolación lineal para interpolar funciones de dos variables
en una malla de dos dimensiones. La idea principal es realizar una interpolación
lineal en una dirección, y después en la otra.
Para calcular el valor de la variable independiente y, que se encuentra entre dos
valores conocidos x1 y x2 , usando interpolación lineal tenemos que el valor buscado
( y − y1 ) + y1 .
( x2 − x1 ) 2
Si lo hacemos en la otra dirección, es decir, si calculamos el valor de la variable x,
( y − y1 )
( x − x1 ) + x1 .
( y2 − y1 ) 2
MATLAB puede realizar interpolaciones mediante varias técnicas programadas. Hemos utilizado un método programado de MATLAB para realizar el proceso de interpolación a través de la función interp2. Este comando interpola los valores entre
puntos de datos conocidos. Encuentra los valores de una función en dos dimensiones entre datos de puntos intermedios.
Zi = interp2 (X, Y, Z, Xi , Yi ) devuelve la matriz Zi que contiene los valores interpolados correspondientes a la malla determinada por Xi y Yi . interp2 determina la
interpolación de la función de dos dimensiones especificada por las matrices X, Y y
La solución uΩh con la que se inicializa el método de Schwarz proviene de la
siguiente interpolación:
(u0Ωh0 ).
uΩh = IΩ
La interpolación de la solución uΩh0 corresponde a la siguiente figura:
Figura 5.8: Solución uΩh proveniente de la interpolación.
En la primera iteración del método, los valores en las fronteras Γ1 y Γ2 correspon0
den a los valores que la solución interpolada uΩh toma en aquellos nodos, es decir,
U 1 := u0h |Ω1 y U 02 := u0h |Ω2 .
Una vez que se obtienen los valores en los nodos de las frontera Γ1 y Γ2 , el siguiente paso consiste en determinar las matrices Ah , Qh y ∇h restringidas a cada
Usamos el comando spy de MATLAB para trazar el patrón de la matriz Ah . En
un solo dominio, la matriz Ah se ve de la siguiente manera:
Figura 5.9: matriz para un solo dominio
Sin embargo, a los subdominios solapados Ω1 = Ω1 ∪ Γ1 y Ω2 = Ω2 ∪ Γ2 ∪ ∂Ω les
corresponden las matrices A1h y A2h respectivamente. Estas matrices corresponden
a submatrices de la matriz Ah . En la Figura 5.10 se visualizan las porciones de la
matriz Ah que corresponden a las submatrices A1h y A2h .
Figura 5.10: matrices para dominios solpados
Las matrices A1h y A2h tienen entradas comunes en los nodos que corresponden a
la región solapada Ω1,2 .
AhΩ1
A1h =  · · · 
AhΓ1
donde AhΩ1 es la matriz con nodos pertenecientes al subdominio Ω1 solapado
sin la frontera, mientras que AhΓ1 es la matriz de los nodos en la interfaz Γ1 . Para la
matriz A2h del segundo dominio se tiene:
AhΓ2
 ··· 
A2h = 
 Ω2 
 ··· 
Ah∂Ω
donde los coeficientes de AhΓ2 corresponden a los nodos de la frontera Γ2 . AhΩ2 es
la matriz con los coeficientes de los nodos que se encuentran interior del dominio
Ω2 y Ah∂Ω es la matriz de los nodos en la frontera externa del dominio original. De
la Figura 5.10 podemos ver que las matrices A1h y A2h son tridiagonales, representaremos cada diagonal por una línea punteada. Gráficamente, si introducimos las
condiciones de frontera, estas matrices se visualizarían de la siguiente manera:
 .. ... ...
.. .. .. 
donde I ∈ R m×m es la matriz identidad con entradas correspondientes a los m
nodos que se encuentran en las fronteras Γ1 , Γ2 y ∂Ω. Como el primer subdominio
posee una sola frontera con el segundo subdominio, a A1h se le añade únicamente la
matriz I ∈ R m×m en la parte inferior. Por el contrario, como el segundo subdominio
posee dos fronteras, la interfaz con Ω1 y la frontera del dominio general, a A2h se le
añaden dos submatrices I ∈ R m×m en la parte superior e inferior.
Como se indicó previamente, en el primer subdominio existen condiciones de
frontera en la interfaz Γ1 entre subdominios, estas condiciones de frontera se añaden
al lado derecho de cada uno de los subsistemas:
 1s 
= k 
 Û2
s +1 
Û2kM
Û1k1
 Û1km
 h2
 2s − m
 0s − m +1
donde s = M − m es la cantidad de nodos por subdominio solapado menos los
nodos de frontera. Como las fronteras Γ1 y Γ2 corresponden a un anillo de la malla
circular entonces poseen m nodos cada una. Cabe notar que el lado derecho del
segundo sistema h 2k+1 contiene un vector (m × 1) de ceros que corresponde a la
solución en la frontera ∂Ω del dominio general.
Para determinar las particiones que corresponden a cada subdominio en las versiones discretas del gradiente ∇ih y la divergencia Qih , con i = 1, 2, el tratamiento
que se da a las matrices es similar al propuesto para el laplaciano Ahµ .
En el dominio original, las particiones de la divergencia, Qih , para i = 1, 2 correspondientes a cada subdominio se visualizan de la siguiente manera:
Figura 5.11: Particiones de Qh .
Q1h =
Q1h,1 Q1h,2
Q1h = 
. · 
Q2h,1 Q2h,2
Q2h = 
donde 0 ∈ R m×m .
En estas matrices se añaden m filas de ceros (en las partes inferior, Ω1 , y superior,
Ω2 , respectivamente ) que corresponden a los nodos en la frontera. En este caso, a
diferencia del tratamiento de las matrices Aihµ , para i = 1, 2, no existen condiciones
en la frontera del tensor de estrés.
Las particiones del gradiente, ∇ih , para i = 1, 2 correspondientes a cada subdominio, se visualizan de la siguiente manera:
Figura 5.12: Particiones de D h .
D1h =
D1h,1
D1h,2
. · 
D2h =
D2h,1
D2h,2
 ·
Análogamente al caso de la divergencia, se añaden m filas de ceros en las partes
inferior y superior de las matrices que corresponden a los nodos en las fronteras de
los subdominios Ω1 y Ω2 respectivamente.
Una vez que las matrices correspondientes a las versiones discretas del laplaciano, la divergencia y el gradiente se hayan particionado para cada subdominio, se
resuelven los subproblemas (5.19) y (5.20) con el algoritmo de Newton semi-suave
(SSN) (Algoritmo 2).
La solución entregada por el algoritmo, en cada subproblema, corresponde a los
nodos en cada subdominio solapado, es decir, en los M = N
K + t nodos. Como cada
subdominio comparte t nodos con el subdominio adyacente, y estamos resolviendo
los subproblemas paralelamente, esto significa que encontramos la solución, en los
t nodos compartidos, dos veces, al resolver el primer subproblema y el segundo
subproblema independientemente. Para obtener la solución en el dominio general,
se eliminan las soluciones encontradas en los t nodos del dominio Ω1,2 , es decir, se
eliminan las soluciones encontradas dos veces al resolver cada subproblema. Por lo
tanto, al eliminar las t soluciones, en cada subdominio, obtenemos los valores en los
K nodos.
La forma aditiva del método de Schwarz es la que se utilizará para resolver el
problema (4.35) en el dominio de la Figura (5.3) con N = 3540 nodos que se distribuyen en una malla circular discretizada en m − 1 = 59 anillos con m = 60 nodos
por anillo más un nodo en el origen. Es decir, la malla posee 59 anillos con diferentes
radios, t = 360.
Finalmente, cada solución se visualiza de la siguiente manera:
Figura 5.13: Primer subdominio.
Figura 5.14: Solución del problema
Figura 5.15: Segundo subdominio.
Figura 5.16: Solución del problema
La solución en el dominio original se presenta en la Figura 5.17.
Figura 5.17: Solución en el dominio completo (Figura 5.4).
Al inicio de la sección determinamos que la forma aditiva del método es la más
conveniente para nuestros propósitos, ya que nos permite resolver el problema en
paralelo. Sin embargo, es conveniente comentar sobre la convergencia del método
pues, la forma aditiva del método, en términos de convergencia, es más lenta que
la forma multiplicativa. A continuación mostramos una tabla para ejemplificar lo
mencionado anteriormente. El ejemplo se realizó en una malla con 3541 nodos de
discretización y variando la cantidad de nodos en la zona solapada o dominio compartido (Ω1,2 ) entre los dos subdominios. La convergencia fue determinada cuando
el error inicial, en la norma infinito, alcanzó una tolerancia igual a 10−5 .
No. de nodos en Ω1,2
% del dominio total Ω
Cuadro 5.1: Iteraciones del método de Schwarz.
Cuadro 5.2: Tiempo en segundos.
A pesar de que el método multiplicativo converge en un número menor de iteraciones, el método aditivo es más rápido cuando comparamos los tiempos de ejecución.
5.4. Generalización a varios subdominios
En el ejemplo anterior se resolvió el problema de Bingham en dos subdominios,
sin embargo, podemos generalizar el problema (4.35) a múltiples subdominios. Sea
Ω el dominio original del problema. La generalización consiste en dividirlo en subdominios solapados Ωl , l = 1, · · · , K tales que
donde K es el número total de subdominios. En la Figura 5.18 podemos ver un ejemplo del dominio original particionado en más de dos subdominios.
Figura 5.18: Dominio general particionado en 4 subdominios.
A continuación mostramos el algoritmo de Schwarz aplicado al problema de
Bingham en varios subdominios.
Algorithm 6 Algoritmo de Schwarz aplicado al problema de Bingham
Calcular los valores iniciales (u~0 , w
Interpolar la solución inicial.
Fijar s=1, (iteraciones de Schwarz.)
Extracción de los valores de frontera u~sk para cada l-subdominio.
for l = 1 : D do,
donde D es el número de subdominios.
fijar k = 0 (SSN-iterations).
Inicializar (u~sl,k , w~sl,k )
Calcular G (u~s , w
Calcular F(u~sl,k , w
Calcular δ.
~ := u~s + δu y ws~ := w~s + δw
Actualizar usl,k
l,k+1
parar, o fijar k := k + 1 y volver al paso 8.
Reunir la solución de los subdominios.
fijar s = s + 1 e ir al paso 4.
En este caso se formula la siguiente familia de subproblemas:
Para l = 1 el subproblema a resolverse sigue de la siguiente manera:
máx( g−
e , γξ (∇1h U1 k+1 )) ⋆ W1 k+1 − gγ∇h U1 k+1 = 0
paraγ > 0,
en Γ1 .
En este caso, el subdominio Ω1 posee únicamente condiciones en la frontera Γ1 =
∂Ω1 ∩ Ω2 (Figura 5.19).
Figura 5.19: Ω1 solapado.
Para 1 < l < K los subproblemas a resolverse son:
en Ωl ,
Ahl Ul k+1 + Qhl Wl k+1 = h lk+1
= U l −1
en Γl,1 ,
= U kl+1
en Γl,2 ,
Los subproblemas formulados anteriormente poseen dos fronteras, pues, en ca68
da subdominio Ωl , con 1 < l < K, existe una interfaz entre dicho dominio y su
dominio consecutivo y antecesor. En el ejemplo de la Figura 5.18 los dominios Ω2 y
Ω3 poseen dos fronteras cada uno, el dominio Ω2 presenta una interfaz con el dominio Ω1 denominada Γ2,1 y una interfaz Γ2,2 con el dominio Ω3 . (Ver Figuras 5.20
y 5.21).
Para cualquier dominio Ωl , 1 < l < K, las interfaces están dadas por Γl,1 =
∂Ωl ∩ Ωl −1 y Γl,2 = ∂Ωl ∩ Ωl +1.
Figura 5.20: Ω2 solapado.
Figura 5.21: Ω3 solapado.
Finalmente, para l = K tenemos:
AhK UK k+1 + QhK WK k+1 = h K
en ΩK ,
−→k+1
máx( g−
= 0 paraγ > 0
− gγ∇ UK
e , γξ (∇K UK )) ⋆ WK
UK k + 1 = UK − 1 k
en ΓK ,
Este problema es ligeramente diferente, existe la interfaz ΓK entre ΩK y ΩK −1 y
la frontera externa del dominio general ∂Ω.
Figura 5.22: Ω4 solapado.
Para resolver la familia de subproblemas (5.24), (5.25) y (5.26) el procedimiento
de particionamiento de matrices es análogo al realizado para dos subdominios.
Para la matriz Ahl se toman las filas correspondientes a los nodos dentro de los
subdominios Ωl \ Γl . En las entradas correspondientes a los nodos de las fronteras
Γl o Γl,i , con i = 1, 2 se fija el valor de 1 en cada entrada.
En las matrices Qhl y ∇hl se procede de manera similar al escoger las filas que corresponden a los nodos dentro del subdominio l sin la frontera, sin embargo, en los
nodos fronterizos se fija el valor de 0, pues, las condiciones de frontera se imponen
en la versión discreta del laplaciano, Ahl , como se indicó previamente.
Si resolvemos la familia de subproblemas (5.24), (5.25) y (5.26) para K = 4 en el
dominio de la Figura 5.18) con el algoritmo 6 obtenemos:
Figura 5.23: Primer subdominio.
Figura 5.24: Segundo subdominio.
Figura 5.25: Tercer subdominio.
Figura 5.26: Cuarto subdominio.
y la solución en el dominio general se visualiza de la siguiente manera:
Figura 5.27: Solución en el dominio completo (Figura 5.18).
5.4.1. Detalles de implementación
Hemos reformulado el problema general en base a una descomposición del dominio espacial en subdominios. Resolvemos cada subproblema, restringido a cada
subdominio, en pararelo. Para ejecutar el código en paralelo usamos el comando
parfor de MATLAB.
parfor nos permite ejecutar un lazo de iteraciones o loop, en paralelo. Este comando paraleliza una declaración, o una parte de un código que se escribe en paralelo,
ejecutándolo en varios procesadores y/o workers. Revisemos la sintaxis del lazo parfor:
parfor variable = valor inicial: valor final
2: declaraciones · · ·
El toolbox para computación en paralelo de MATLAB, que incluye el coman72
do parfor, ejecuta códigos, programas y aplicaciones en distintos workers (motores
computacionales de MATLAB) que se ejecutan localmente. Es decir, permite al usuario ejecutar un trabajo en paralelo usando diferentes labs o workers con copias de
MATLAB adicionales que asisten a la copia original o cliente. Si la computadora de
cálculo tiene múltiples procesadores, los labs los activarán. Por lo tanto, sin cambiar
el código, se puede ejecutar una misma aplicación en un cluster de computadora.
Esta herramienta computacional de MATLAB permite resolver problemas usando procesadores multi-núcleo o clusters de computadoras. Además, este toolbox paraleliza aplicaciones de MATLAB sin necesidad de programación MPI (Message
Passing Interface). MPI es una especificación para programación de paso de mensajes, proporciona librerías que son empleadas en programas para comunicar datos
entre procesos. Permite desarrollar códigos y programas que pueden ser migrados
a diferentes computadores en paralelo.
MATLAB Distributed Computing Server, DCS, permite computación en paralela
en clusters, así, las tareas se distribuyen a distintos labs. DCS constituye una especificación similar al paradigma de memoria distribuida MPI. MATLAB usa una versión de MPI denominada MPICH2. Sin embargo, el usuario no necesita programar
explícitamente el paso de mensajes mediante MPI.
Para iniciar la computación en paralelo de MATLAB se requiere ejecutar una
sesión en paralelo. Para ello, al iniciar matlabpool a través de la declaración open
matlabpool se permite trabajar en sesiones interactivas con una configuración dada
por el usuario. Una vez iniciado el ambiente en paralelo, es posible usar el lazo
parfor. Este lazo es la versión para computación en paralelo del comando tradicional
for. Por lo tanto, dentro del toolbox en paralelo, simplemente se usa parfor en lugar
de for. Es importante notar que la variable del lazo puede ser únicamente un vector
con enteros positivos consecutivos.
parfor i = 1 : 100
2: código · · ·
En este caso, la variable i corresponde a un vector de enteros consecutivos del 1 al
100. Además, el código dentro del lazo parfor es la porción de código que se ejecutará
en paralelo. Se debe tomar en cuenta que, al usar parfor, todas las iteraciones del lazo
deben ser independientes unas de otras. Esto se debe a que el lazo no ejecuta las
iteraciones consecutivamente sino, en otro orden decidido por MATLAB. De esta
manera, no se puede trabajar con variables que se modifiquen dentro del proceso
pues, no se pueden modificar dentro de todos los procesadores y al mismo tiempo.
Es importante notar también que el código dentro del lazo no puede depender de la
variable de conteo.
El lazo parfor puede ser usado únicamente en un nivel, es decir, no puede existir
un parfor dentro de otro lazo parfor. Sin embargo, el lazo tradicional for sí puede estar
incluido en la porción del código. Revisaremos un ejemplo sencilllo del uso del lazo
parfor.
x=zeros(100,5);
y=zeros(1,5);
y(j)= i;
x(i,:)=y;
A continuación revisemos un ejemplo del uso incorrecto del lazo:
parfor i = 1 : 50
En este caso, el código no se ejecutará, pues, estamos usando fuera del lazo parfor la
variable x, que se modifica dentro del lazo. Incluso si es una constante, el código no
Si revisamos nuevamente nuestro problema (4.35), al paralelizarlo usando parfor,
el pseudocódigo seguiría de la siguiente manera:
~0 , w
1: Calcular los valores iniciales (u
parfor l = 1 : D, ( D es el número de subdominios y además, el número de
parar, o fijar k := k + 1 y volver al paso 7.
El lazo parfor nos permite ejecutar las iteraciones en paralelo. Por lo tanto, paraleliza la porción de código que incluye las declaraciones para resolver el problema
de Bingham en cada subdominio simultáneamente. Es decir, ejecutándolo en varios
procesadores o workers. De esta manera, cada uno de los subproblemas, restringido
a cada subdominio, se resuelve, simultáneamente, en un procesador diferente.
5.4.2. Resultados numéricos de la paralelización del fluido de Bingham en varios subdominios
A continuación se presentan tres experimentos numéricos que muestran el comportamiento del algoritmo paralelizado. Se muestra el porcentaje de reducción en el
tiempo respecto al tiempo de resolución del problema sin descomposición de dominio.
E JEMPLO 2. Este ejemplo se realizó con 40400 nodos y el tiempo de ejecución del
algoritmo, para un solo dominio, fue de 58,49 segundos.
No. nodos por
subdominio solapado
% de reducción respecto
Cuadro 5.3: Comportamiento del algoritmo paralelizado para el ejemplo (2)
E JEMPLO 3. Este ejemplo se realizó con 160800 nodos y el tiempo de ejecución del
algoritmo, para un solo dominio, fue de 260,16 segundos.
Cuadro 5.4: Comportamiento del algoritmo paralelizado para el ejemplo (3)
E JEMPLO 4. Este ejemplo se realizó con 361200 nodos y el tiempo de ejecución del
algoritmo, para un solo dominio, fue de 1814,35 segundos.
Cuadro 5.5: Comportamiento del algoritmo paralelizado para el ejemplo (4)
Figura 5.28: Problema resuelto con 40400 nodos en 50 subdominios.
Figura 5.29: Problema resuelto con 40400 nodos en un dominio.
Figura 5.30: Problema resuelto con 160800 nodos en 100 subdominios.
Figura 5.31: Problema resuelto con 160800 nodos en un dominio.
Figura 5.32: Problema resuelto con 361200 nodos en 100 subdominios.
Figura 5.33: Problema resuelto con 361200 nodos en un dominio.
Velocidad máxima -solución
Cuadro 5.6: Comparación con la solución exacta
5.5. Resultados de convergencia
En esta sección discutiremos el comportamiento del algoritmo en términos de
convergencia y comentaremos resultados generales del método de Schwarz. Si tenemos un dominio como el que se muestra en la Figura 5.34, para operadores lineales
y elípticos se espera: (Ver [2])
Figura 5.34: Tamaño del dominio, solapamiento y elementos. Ver [2]
• el número de iteraciones crezca a medida de 1/H,
• si δ es proporcional a H, el número de iteraciones está acotado y no depende
de h, es decir, si se fijan los subdominios, el número de iteraciones requerido
es independiente del grado de refinamiento de la discretización,
• la tasa de convergencia no es independiente de H, es decir, depende de las
dimensiones de la región solapada, si δ crece la tasa de convergencia mejora.
Cabe recalcar que la familia de sistemas (5.24), (5.25) y (5.26) que estamos resolviendo es no lineal y los resultados de convergencia presentados anteriormente no
necesariamente se cumplen, sin embargo, a continuación mostramos varios experimentos numéricos sobre el comportamiento del método aplicado al problema de
Bingham y en concreto al ejemplo con N = 3540. En la siguiente tabla se muestran
dos escenarios, se fija el número de subdominios variando la cantidad de nodos en
las zonas solapadas y en el segundo escenario se mantiene constante la cantidad de
nodos en las zonas solapadas y se varía el número de subdominios.
No. de subdominios
4 (885 nodos x
5 (708 nodos x
6 (590 nodos x
No. de nodos en las
zonas solapadas
Cuadro 5.7: Iteraciones de Schwarz
Si el número de nodos en las regiones solapadas aumenta, el número de iteraciones de Schwarz se reduce manteniendo fija la cantidad de subdominios, sin embargo, si se mantiene la cantidad de nodos fija en las zonas solapadas pero se varía el
número de dominios, las iteraciones de Schwarz aumentan.
A continuación presentamos el siguiente experimento, se particiona el dominio
del Ejemplo con N = 3540 en 4, 5 y 6 subdominios, en lugar de mantener fija la
cantidad de nodos en la zona solapada, esta se mantiene proporcional a la cantidad
de nodos por subdominio.
% de nodos en las zonas
Cuadro 5.8: Iteraciones de Schwarz
En el siguiente experimento se muestra el comportamiento del método variando
el número de nodos en el dominio general, es decir, se trabajará con tres discretizaciones diferentes y el mismo número de subdominios, en este caso, 10 subdominios.
El dominio del ejemplo con N = 3540 se discretizará en un número fijo de subdominios y se mantendrá, proporcional al número de nodos por subdominio, la cantidad
de nodos en la región solapada.
% solapamiento
It. Schwarz
Cuadro 5.9: Iteraciones de Schwarz
Si el paso de discretización varía, pero, se mantiene constante la proporción de
nodos en la zona solapada y el número de subdominios, vemos que el número de
iteraciones no mejora significativamente.
El objetivo principal de este trabajo consistió en plantear un esquema paralelizado para simular el flujo de un fluido de Bingham en la sección de una tubería
cilíndrica. Para lograr este objetivo, a su vez, se trataron dos temas: la geometría
en la que se resolvió el problema (sección transversal de una tubería) y el planteamiento de la técnica de descomposición de dominio. El sistema a resolver se había
planteado en trabajos como [8] en el sistema ortonormal de ejes 0x1 x2 x3 . Por lo tanto,
para obtener la simulación en una tubería, se transformó el sistema de ecuaciones
que modelan el flujo de Bingham bidimensional a coordenadas polares. Una vez
transformado el sistema, para la discretización espacial, se propuso un esquema en
difrencias finitas centradas para coordenadas polares.
El segundo reto de este trabajo consistió en plantear la técnica de descomposición
de dominio de Schwarz para el problema de Bingham en coordenadas polares. A
través de esta técnica, se descompuso el dominio espacial del problema y se generó
una familia de subproblemas restringidos a los nuevos subdominios. Con esto, se
logró replantear el problema, de un tamaño determinado, en varios subproblemas
cuyo tamaño fue menor que el problema original.
En la literatura, ([2] y [25]), se encuentra la técnica de Schwarz en sus formas aditiva y multiplicativa. Se plantearon las dos versiones para analizar su aplicación al
problema de Bingham. Sin embargo, la forma multiplicativa del método no permite
la construcción de un esquema paralelizado, pues, se resuelven los subproblemas
generados de manera secuencial. La forma aditiva del método resultó ideal para
resolver los subproblemas de manera independiente y, por lo tanto, en paralelo.
Se realizaron experimentos numéricos para analizar la convergencia del algoritmo de Schwarz. En los experimentos numéricos, se encontró que si el número de
nodos en las regiones solapadas aumenta, el número de iteraciones de Schwarz se
reduce manteniendo fija la cantidad de subdominios, sin embargo, si se mantiene la
cantidad de nodos fija en las zonas solapadas pero, se varía el número de dominios,
las iteraciones de Schwarz aumentan.Por lo tanto, concluimos experimentalmente
que la convergencia del método depende del tamaño de la región solapada, es decir,
de la proporción de nodos que los subdominios comparten.
En los experimentos en los que se fijó inicialmente el número de iteraciones de
Schwarz, se comparó el tiempo del algoritmo paralelizado con el tiempo de resolución del problema sin paralelizar, es decir, en un solo dominio. Fijando el número
de iteraciones de Schwarz se obtuvo aproximadamente el 90 % de reducción en el
tiempo de ejecución en la resolución del problema de Bingham.
[1] Ben Andrews y Christopher Hopper. The ricci flow in riemannian geometry.
Lect. Notes Math, 2011.
[2] Petter Bjorstad y William Gropp. Domain decomposition: parallel multilevel methods for elliptic partial differential equations. Cambridge University Press, 2004.
[3] Daniel Bonn y Morton M Denn. Yield stress fluids slowly yield to analysis.
Science, 324(5933):1401–1402, 2009.
[4] Alexandre Joel Chorin, Jerrold E Marsden, y Jerrold E Marsden. A mathematical
introduction to fluid mechanics, volume 3. Springer, 1990.
[5] JD Cole y EM Murman. Calculation of plane steady transonic flows. AIAA
Journal, 9(1):114–121, 1971.
[6] Juan Carlos De los Reyes. Numerical PDE-Constrained Optimization. Springer,
[7] Juan Carlos De Los Reyes y Sergio González-Andrade. Numerical simulation
of thermally convective viscoplastic fluids by semismooth second order type
methods. Journal of Non-Newtonian Fluid Mechanics, 193:43–48, 2013.
[8] Juan Carlos De Los Reyes y Sergio González. Path following methods for
steady laminar bingham flow in cylindrical pipes. ESAIM: Mathematical Modelling and Numerical Analysis, 43:81–117, 1 2009.
[9] Maksymilian Dryja. An additive schwarz algorithm for two-and threedimensional finite element elliptic problems. Tony Chan, Roland Glowinski, Jacques P eriaux, and Olof Widlund, editors, Domain Decomposition Methods, Philadelphia, PA, 1989.
[10] Massimo Fornasier, Andreas Langer, y Carola-Bibiane Schönlieb. A convergent
overlapping domain decomposition method for total variation minimization.
Numerische Mathematik, 116(4):645–685, 2010.
[11] Massimo Fornasier y Carola-Bibiane Schönlieb. Subspace correction methods
for total variation and l1-minimization. SIAM Journal on Numerical Analysis,
47(5):3397–3428, 2009.
[12] Martin Fuchs y Gregory Seregin. Regularity results for the quasi–static bingham variational inequality in dimensions two and three. Mathematische Zeitschrift, 227(3):525–541, 1998.
[13] Roland Glowinski. Numerical methods for nonlinear variational problems, volume 4. Springer, 1984.
[14] Faramarz Gordaninejad. Proceedings of the 10th International Conference on Electrorheological Fluids and Magnetorheological Suspensions: Lake Tahoe, USA, 18-22
June, 2006. World Scientific Publishing Company Incorporated, 2007.
[15] Bertil Gustafsson. Fundamentals of Scientific Computing, volume 8. Springer,
[16] Michael Hintermüller.
Semismooth newton methods and applications.
Oberwolfach-Seminar, 2010.
[17] Michael Hintermüller, Kazufumi Ito, y Karl Kunisch. The primal-dual active
set strategy as a semismooth newton method. SIAM Journal on Optimization,
13(3):865–888, 2002.
[18] Richard W Johnson. Handbook of fluid dynamics. Crc Press, 1998.
[19] Randall J LeVeque y Randall J Le Veque. Numerical methods for conservation laws,
volume 132. Springer, 1992.
[20] G Duvaut JL Lions. Inequalities in mechanics and physics. Springer, 1976.
[21] Pierre-Louis Lions. On the schwarz alternating method. i. In First international
symposium on domain decomposition methods for partial differential equations, págs.
1–42. Paris, France, 1988.
[22] Kamran Mohseni y Tim Colonius. Numerical treatment of polar coordinate
singularities. Journal of Computational Physics, 157(2):787–795, 2000.
[23] PP Mosolov y VP Miasnikov. Variational methods in the theory of the fluidity of a viscous-plastic medium. Journal of Applied Mathematics and Mechanics,
29(3):545–577, 1965.
[24] MA Moyers-Gonzalez y IA Frigaard. Numerical solution of duct flows of multiple visco-plastic fluids. Journal of non-newtonian fluid mechanics, 122(1):227–
241, 2004.
[25] Alfio Quarteroni y Alberto Valli. Domain decomposition methods for partial differential equations. Number CMCS-BOOK-2009-019. Oxford University Press,
[26] Xue-Cheng Tai y Jinchao Xu. Global and uniform convergence of subspace
correction methods for some convex optimization problems. Mathematics of
Computation, 71(237):105–124, 2002.
[27] Michael Ulbrich. Nonsmooth Newton-like methods for variational inequalities and
constrained optimization problems in function spaces. PhD thesis, Habilitation thesis, Fakultät für Mathematik, Technische Universität München, 2002.
[28] Michael Ulbrich. Semismooth Newton methods for variational inequalities and constrained optimization problems in function spaces. SIAM, 2011.
Incertidumbre: Los Axiomas de la probabilidad
Transitividad • Redes de Decisión • Teoría de la utilidad • Sustitución • Ciclos Dirigidos • Redes de Creencia • Dominio Incierto • Probabilidad Condicional
Bimestrales - Colegio del Valle
Red multipunto - dgeti quintana roo
Ciudades Duales • Nodos metrapolitanos o glocales • Telecomunicaciones
Receptor • Transmisor • Simplex • Capa de Enlace de Datos • Niveles de comunicación
COLEGIO DE SEÑORITAS “EL SAGRADO CORAZÓN” EL NARANJO
COLEGIO SAN JOSÉ DE CHICUREO Dirección Académica
LEYES DE NEWTON ISAAC NEWTON:
Acción y reacción • Mecánica • Dinámica clásica • Principio de inercia

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución