MAT-Serie A, 5 (2001), 11-26 11

Método de Diferencias Finitas para un Problema de Bingham Unidimensional

Germán A. Torres ^{∗} y Cristina Turner ^{†} Facultad de Matemática, Astronomía y Física

Universidad Nacional de Córdoba, 5000 Córdoba, Argentina

Abstract

In this work, a numerical method of finite differences is proposed to solve an unidimensional Bingham problem, and some theoretical properties of the solution are proved, and corroborate numerical results. A Bingham fluid is a non-newtonian fluid, whose viscous behaviour makes his layers move only if shear stress is greater than a threshold value τ_{0} . There are some previous theoretical results that allows us to affirm existence and uniqueness of the solution under certain initial conditions. The proposed method is a finite difference scheme with spatial variable step, that is, while we move a fixed time step, we adjust the grid in such a way that the spatial step represents the free boundary advance. Using an internal calculus of fixed point is possible to reduce the number of flops, making feasible the implementation of the method in a computer. Besides that, a theorem of existence and uniqueness is proved for more general cases, and also a result about asymptotic behaviour is demonstrated.

Resumen

En este trabajo se propone un método numérico de diferencias finitas para resolver un problema de Bingham unidimensional y se demuestran algunas propiedades teóricas de la solución que sirven para corroborar los resultados numéricos. Un fluido de Bingham es un fluido no newtoniano cuyo comportamiento viscoso hace que sus capas se desplacen entre sí sólo si la tensión tangencial supera una cierta tensión umbral τ_{0}. Hay varios resultados teóricos previos que nos permiten afirmar la existencia y unicidad de la solución bajo ciertas condiciones iniciales. El método consiste en el uso de diferencias finitas, más precisamente en la discretización de la variable temporal con un paso fijo, mientras que en la variable espacial se elige un paso variable, de modo tal que la grilla se vaya adaptando a la frontera libre a medida que el tiempo transcurre. Usando una iteración interna de punto fijo es posible reducir notablemente el costo computacional, haciendo factible la implementación práctica del método. Además se prueba un resultado de existencia y unicidad para problemas más generales y se demuestra cuál es el comportamiento asintótico.

Key words: Bingham fluid, finite differences, non-newtonian fluids Palabras Claves: Fluido de Bingham, diferencias finitas, fluidos no newtonianos AMS Subject Classification: 35A40, 35B40, 35R35, 39A12, 65M06

∗torres@mate.uncor.edu †turner@mate.uncor.edu

1 Introducción

Consideramos el fluido encerrado entre dos placas paralelas, y mediante el uso de la ecuación de Navier-Stokes para la región viscosa y la ecuación de movimiento de Newton para la parte rígida es posible modelar el comportamiento del sistema. La frontera que separa las dos regiones previamente mencionadas (la viscosa y la rígida) es desconocida y se encuentra en movimiento, y es una de las incógnitas más importantes del problema.

Existe una formulación débil, en forma variacional, para problemas de frontera libre del tipo del fluido de Bingham. Se puede consultar en [4],[5], [8], [9] y [10]. Además en [11] hay detallada una extensa bibliografía acerca del tema.

Según sea la relación entre el tensión de corte y la tensión tangencial, es posible realizar una clasificación de los fluidos. Si la tensión tangencial es un múltiplo de la tensión de corte (la tensión de corte es la velocidad de cambio de la velocidad de las capas del fluido), donde el coeficiente multiplicativo es la viscosidad, entonces el fluido es denominado newtoniano, caso contrario se llamará no newtoniano. Ejemplos de fluidos newtonianos son: agua, benceno, alcohol, glicerina. Ejemplos de fluidos no newtonianos son: plasma de la sangre, chocolate, mayonesa, salsa de tomate, pasta dental, algunas grasas, aguas cloacales, y en general líquidos con partículas en suspensión.

Entre los fluidos no newtonianos tenemos los fluidos de Bingham, que son aquéllos cuya relación entre la tensión tangencial τ y la tensión de corte σ es lineal. Es decir,

τ = τ_{0} + ησ.  (1.1)

donde η es la viscosidad y τ_{0} es la tensión umbral.

Suponemos que el fluido es incompresible, laminar, y con densidad ρ constante. Fijado el sistema de coordenadas, donde la dirección y es la vertical, la x es la dirección de movimiento, y z la restante, se hacen las siguientes hipótesis:

1. El gradiente de presión se aplica en un solo sentido, es decir, \frac{∂p}{∂y} = \frac{∂p}{∂z} = 0.

2. El fluido es laminar, es decir, v = w = 0.

3. La componente no nula de la velocidad es sólo función del tiempo y de la componente

vertical de la posición, es decir, \frac{∂u}{∂z} = \frac{∂u}{∂x} = 0.

4. No hay intercambio de material entre la región viscosa y la región rígida, es decir, u_{y}(s(t), t) = 0.

5. El fluido tiene la condición de adherencia en las paredes de las placas, es decir, que la velocidad del fluido allí es cero.

Usando todas estas hipótesis simplificativas se llega a un sistema de ecuaciones diferenciales parciales, que llamaremos Problema (P). Haciendo un simple cambio de variables, podemos independizarnos de algunas constantes quedando:

u_{t} − u_{yy} = f(t),  s(t) < y < 1, t > 0,  (1.2)
u(1, t) = 0,  t > 0,  (1.3)

u(y, 0) = u_{0}(y),  s(0) = s_{0}, 0 < s_{0} < y < 1,  (1.4)
u_{y}(s(t), t) = 0,  t > 0,  (1.5)
u_{t}(s(t), t) = f(t) − \frac{τ_{0}}{s(t)},  t > 0.  (1.6)

Como vemos, en el sistema ha quedado un problema del tipo de conducción del calor, donde el término fuente es el gradiente de presión f que, según las hipótesis, depende sólo de la variable temporal. Este sistema es llamado un problema de frontera libre puesto que la función s(t), que es la frontera que separa dos regiones, forma parte de las incógnitas. Para que haya desplazamientos de las capas del fluido, debemos suponer que el gradiente de presión supera la tensión umbral τ_{0}, es decir,

f(t) > τ_{0}, ∀ t > 0.  (1.7)

También se puede colocar una condición más general en lugar de tener una condición de contorno fija (como en (1.3)), simulando la situación de dos fluidos de distinta naturaleza en contacto, quedando una ecuación alternativa:

u(1, t) = g(t),  t > 0.  (1.8)

Se puede observar que tal función g no acarrea mayores problemas puesto que se puede reescribir el sistema considerando una nueva función para el gradiente de presión.

Un sencillo cálculo de derivación nos permite obtener el problema derivado, que nos será de utilidad más adelante.

Si w = u_{y} se tiene el problema (P_{y})

w_{t} − w_{yy} = 0,  s(t) < y < 1, t > 0,  (1.9)
w_{y}(1, t) = − f(t),  t > 0,  (1.10)
w(y, 0) = u0^{′}(y),  s(0) = s_{0}, 0 < s_{0} < y < 1,  (1.11)
w(s(t), t) = 0,  t > 0,  (1.12)
w_{y}(s(t), t) = − \frac{τ_{0}}{s(t)},  t > 0.  (1.13)

2 Resultados previos

Si suponemos que el gradiente de presión tiende a un valor estacionario, es decir, lim_{t→∞} f(t) = f_{0}, y suponemos u independiente del tiempo en (1.2)-(1.6), se puede obtener la solución estacionaria.

u_{∞}(y) = −\frac{f_{0}}{2} (y^{2} − 1) + f_{0}s_{∞} (y − 1) ,  s_{∞} = \frac{τ_{0}}{f_{0}} .  (2.1)

Elena Comparini probó en [1] la existencia y unicidad del problema (P) para todo tiempo, dependencia continua respecto de los datos, un resultado de monotonicidad y la convergencia asintótica, cuando el gradiente de presión tiende a un valor constante y cuando las condiciones iniciales cumplen ciertas hipótesis. Llamemos

u_{0}(1) = 0, u0^{′}(y) ≤ 0, u0^{′′}(y) ≤ 0, (A1)

u0^{′}(y) ≥ −f_{0} (y − \frac{τ_{0}}{f_{0}} ) , u0^{′′}(y) ≥ −f_{0}. (A2)

Proposición 2.1 Si s_{0} > s_{∞}, (A1), (A2), f′ ≥ 0 y u0^{′′′}(y) ≤ 0 se tiene existencia y unicidad del problema (P) para todo tiempo y

s(t) > s_{∞},  t > 0,  (2.2)
0 ≤ u(y, t) ≤ u_{∞}(y),  s(t) < y < 1, t > 0,  (2.3)
u_{∞y}(y) ≤ u_{y}(y, t) ≤ 0,  s(t) < y < 1, t > 0,  (2.4)
f(t) − f_{0} ≤ u_{t}(y, t) ≤ f(t) ≤ f_{0},  s(t) < y < 1, t > 0,  (2.5)
s′(t) < 0,  t > 0.  (2.6)

Se observa además que la siguiente integral está acotada.

∫_{0}^{∞} (\frac{s_{∞}}{s(t)} − \frac{f(t)}{f_{0}}) dt.  (2.7)

Esto nos dice que el integrando converge a cero cuando t tiende a infinito, concluyendo que la frontera libre es, asintóticamente, inversamente proporcional al término fuente multiplicada por la tensión umbral.

s(t) ≈ \frac{τ_{0}}{f(t)},  t → ∞.  (2.8)

3 Generalizaciones

En esta sección se generalizarán los resultados expuestos en [1]. Llamemos (B1) a:

f(t) ≥ 0 ∀ t ∈ [0, T],  u_{0}(1) = 0,  u0^{′}(y) ≤ 0 ∀ y ∈ [s_{0}, 1].  (B1)

Notar que (B1) dice que u_{0}(y) ≥ 0.

Proposición 3.1 Se puede ver que:

1. Si la hipótesis (B1) se cumple, entonces

u(y, t) ≥ 0,  (3.1)
u_{y}(y, t) ≤ 0.  (3.2)

2. Si se cumple (B1), u0^{′′′}(y) ≤ 0 ∀ y ∈ [s_{0}, 1] y f′(t) ≥ 0 ∀ t > 0, entonces s′(t) < 0 para los t donde la solución exista y para los t tales que s(t) > 0.

3. Si u es solución de (1.2)-(1.6) se cumple la siguiente relación

∫_{0}^{t} [f(τ) − \frac{τ_{0}}{s(τ)}] dτ = u(s(t), t) − u_{0}(s_{0}).  (3.3)

para los t donde la solución exista.

Demostración: (1) y (2) son una consecuencia obvia del Principio del Máximo y el Lema de Hopf y (3) resulta de las ecuaciones de Green.

□

Observación 1 Se observa que la acotación de la integral depende de la velocidad del fluido en la frontera libre. Pero ahora no podemos asegurar que el integrando tienda a cero, puesto que los experimentos numéricos del método que se propondrá a continuación revelaron que la velocidad del fluido en la frontera libre tiende a infinito.

Proposición 3.2 Supongamos que se cumple la hipótesis (B1), que u0^{′′′}(y) ≤ 0, ∀ y ∈ [s_{0}, 1], y que f′(t) ≥ 0, ∀ t > 0. Entonces

w_{yy}(y, t) ≤ 0,  s(t) < y < 1, t > 0,  (3.4)
w_{y}(y, t) ≤ 0,  s(t) < y < 1, t > 0,  (3.5)
w_{t}(y, t) ≤ 0,  s(t) < y < 1, t > 0,  (3.6)
s(t) > \frac{τ_{0}}{f(t)},  t > 0.  (3.7)

siempre que s(t) > 0.

Demostración: La función v ≐ w_{yy} resuelve el siguiente sistema de ecuaciones,

v_{t} − v_{yy} = 0,  s(t) < y < 1, t > 0,  (3.8)
v_{y}(1, t) = −f′(t) ≤ 0,  t > 0,  (3.9)
v(y, 0) = u0^{′′′}(y) ≤ 0,  0 < s_{0} < y < 1,  (3.10)
v(s(t), t) = \frac{τ_{0}s′(t)}{s(t)} ≤ 0,  t > 0.  (3.11)

Por el lema de Hopf la solución no tiene máximo en y = 1. Luego es fácil ver por el Principio del Máximo que v ≤ 0, lo cual nos dice que w_{yy} ≤ 0. Como w_{yy} = w_{t}, entonces w_{t} ≤ 0. También, como w_{yy} ≤ 0 resulta que w_{y}(y, t) es decreciente respecto de y para cada t. Luego el máximo de w_{y} está en la frontera libre. Como w_{y}(s(t), t) = −\frac{τ0}{s(t)} ≤ 0 entonces w_{y} ≤ 0. Otra consecuencia de ser w_{y} decreciente es que w_{y}(s(t), t) ≥ w_{y}(1, t), lo cual, reemplazando sus valores respectivos da que s(t) > \frac{τ_{0}}{f(t)}. Esto concluye la demostración.

□

Teorema 3.3 Sea f una función diferenciable que cumple f(t) > τ_{0} ∀ t > 0, f′(t) ≥ 0 ∀ t > 0, lim_{t⟶∞} f(t) = ∞. Asumimos también que se cumplen la hipótesis (C1) conformada por estas cinco ecuaciones:

\{ \begin{matrix}u_{0}(1) = 0, & u0^{′}(y) ≤ 0 ∀ y ∈ [s_{0}, 1], & u0^{′′}(y) ≤ 0 ∀ y ∈ [s_{0}, 1], \\ s_{0} > 0, & u0^{′′′}(y) ≤ 0 ∀ y ∈ [s_{0}, 1] & . \end{matrix} (C1)

Luego existe única solución del problema (P_{y}) para todo tiempo y además

lim_{t⟶∞} s(t) = 0,  s(t) > \frac{τ_{0}}{f(t)} ∀ t > 0.  (3.12)

Demostración: Llamamos problema (P_{y}^{n}) al sistema (P_{y}) salvo que el término f(t) es reemplazado por f_{n}(t), donde, para n ∈ N y ϵ > 0 arbitrario

f_{n}(t) = \{ \begin{matrix}f(t), & t ≤ n \\ polinomio que empalme a f en forma diferenciable, & n < t < n + ϵ \\ f(n + ϵ), & t ≥ n + ϵ \end{matrix} (3.13)

Como la función f tiene límite infinito, a partir de un N en adelante se cumplen las hipótesis de existencia y unicidad de la Proposición (2.1), lo cual nos dice que cada problema (P_{y}^{n}) con n ≥ N tiene solución única para t arbitrariamente grande, cumpliéndose también que

lim_{t⟶∞} s^{n}(t) = \frac{τ_{0}}{f(n + ϵ)}, s^{n}(t) > \frac{τ_{0}}{f(n + ϵ)}.  (3.14)

A la solución de este problema la llamamos {w^{n}, s^{n}}.

Sea T un número positivo arbitrario. Tomo n tal que n ≥ N y n ≥ T. Luego el problema (P_{y}) coincide con el problema (P_{y}^{n}) en el intervalo [0, T]. Luego, de acuerdo al párrafo anterior, el problema (P_{y}) tiene única solución en el intervalo [0, T] para todo T arbitrario. Esto me dice que el problema (P_{y}) admite solución para t arbitrariamente grande.

A partir de ahora consideraremos que n ≥ N para que se cumplan las hipótesis de existencia y unicidad. Llamamos {w^{n}, s^{n}} solución de (P_{y}^{n}) para n ≥ N y llamaremos {w, s} solución de (P_{y}). Probaremos que esa sucesión de soluciones converge a la solución general.

Por unicidad de la solución, y debido a que f_{n}(t) = f(t) para t ≤ n, se tiene que

\{ \begin{matrix}s^{n}(t) = s(t) & t ≤ n \\ w^{n}(y, t) = w(y, t) & s(t) ≤ y ≤ 1, t ≤ n \end{matrix}  (3.15)

Esto dice que lim_{n⟶∞} s^{n}(t) = s(t). Por un teorema de monotonía para problemas de este tipo, se sabe que si f_{n} < f_{n+1} entonces s^{n} > s^{n+1}. Ver [1]. También sabemos que lim_{t⟶∞} s^{n}(t) = \frac{τ_{0}}{f(n+ϵ)} para todo n.

Sea ϵ_{1} > 0. Luego existe N_{1} ≥ N tal que si n ≥ N_{1} entonces \frac{τ_{0}}{f(n+ϵ)} < ϵ_{1}. En particular \frac{τ_{0}}{f(N_{1}+ϵ)} < ϵ_{1}. Luego existe M, un número real, tal que si t ≥ M, entonces s^{N_{1}}(t) < ϵ_{1}. Por monotonía, s^{n}(t) < s^{N_{1}}(t) < ϵ_{1} para n ≥ N_{1} y t ≥ M. Dado t ≥ M, tomo n ≥ t, luego s^{n}(t) = s(t), que da como resultado que s(t) < ϵ_{1}. Notar que s(t) > 0, ∀ t > 0. Haciendo uso de la Proposición (3.2) se obtiene por (3.7) el otro resultado que buscábamos. Esto concluye la demostración.

□

4 Método de Diferencias Finitas con Paso Espacial Variable

Por simplicidad trabajaremos con el sistema (1.9)-(1.13). El método consiste en discretizar la ecuación principal. Fijado un paso de tiempo Δt, se construye una grilla de tal forma que la frontera libre pase por los nodos de ésta. En cada paso de tiempo se calcula un paso espacial Δy adecuado para tal fin. El algoritmo se ha construido para el caso cuando la frontera libre es decreciente. Para métodos similares ver [3], [6] y [7]. Definimos entonces Δt arbitrario y

y_{i+1} = 1 − (Δy_{1} + ... + Δy_{i}) ,  i ≥ 1, y_{1} = 1,  (4.1)
t_{n} = (n − 1) Δt,  n ≥ 1.  (4.2)

Dividimos el intervalo [s_{0}, 1] en m−1 partes, donde m es un número natural arbitrario mayor o igual que tres. Definimos

Δy_{i} = \frac{1 − s_{0}}{m − 1},  i = 1, ..., m − 1.  (4.3)

Con esta elección los puntos y_{1},...,y_{m} son una partición de [s_{0}, 1] igualmente distribuida. Se ve claramente que y_{m} = s_{0}. De aquí en adelante construiremos Δy_{m}, Δy_{m+1}, ..., de manera tal que (y_{n+m−1}, t_{n}) sea un punto de la frontera libre. Denotamos w_{in} ≐ w(y_{i}, t_{n}). Realizamos aproximaciones mediante diferencias finitas para las derivadas y reemplazando en (1.9)-(1.13) queda lo que llamamos problema (P d_{y}):

w_{1,n+1} − w_{2,n+1} = −Δy_{i}f_{n+1},  n ≥ 1,  (4.4)

\frac{Δt}{Δy_{i−1}} w_{i−1,n+1} − (\frac{Δt}{Δy_{i−1}} + \frac{Δt}{Δy_{i}} + Δy_{i}) w_{i,n+1} + \frac{Δt}{Δy_{i}} w_{i+1,n+1} =
= −Δy_{i}w_{i,n},  i = 2...n + m − 1, n ≥ 1,  (4.5)

w_{n+m,n+1} = 0,  n ≥ 1,  (4.6)
w_{i1} = u0^{′}(y_{i}),  i = 1, ..., m,  (4.7)

Δy_{n+m−1} = (\frac{w_{n+m−1,n+1}}{w_{n+m−1,n+1} − τ_{0}} ) y_{n+m−1},  n ≥ 1.  (4.8)

La idea del proceso es la siguiente:

1. Elegidos y_{1},...,y_{m}, se puede calcular {w_{i1}}_{i=1}^{m} mediante la ecuación (4.7). De este modo el primer paso quedaría resuelto.

2. Hasta el momento sólo disponemos de Δy_{1},...,Δy_{m−1}. Si pudiéramos conocer el valor de Δy_{m}, entonces utilizando el sistema (4.4)-(4.6) para n = 1, entonces sería posible calcular {w_{i2}}_{i=1}^{m+1}. Luego nos dedicaremos a calcular Δy_{m}.

3. Esto nos permite realizar un paso inductivo. Supongamos conocidos los valores {w_{in}}_{i=1}^{n+m−1} y Δy_{1},..., Δy_{n+m−2}. Si supiéramos cuánto vale Δy_{n+m−1}, entonces podríamos resolver el sistema lineal (4.4)-(4.6) y obtener {w_{i,n+1}}_{i=1}^{m+n}.

Como vemos, lo único que nos queda es tratar de encontrar el valor de Δy_{n+m−1}, de tal forma que sea consistente con la ecuación (4.8). Lo que proponemos es lo siguiente:

1. Dar un valor inicial Δy_{n+m−1}^{(1)} tal que 0 < Δy_{n+m−1}^{(1)} < y_{n+m−1}.

2. Calcular una iteración interna en r determinada por la ecuación (4.8), es decir,

Δy_{n+m−1}^{(r+1)} = (\frac{w_{n+m−1,n+1}^{(r)}}{w_{n+m−1,n+1}^{(r)} − τ_{0}} ) y_{n+m−1},  n ≥ 1.  (4.9)

donde w_{n+m−1,n+1}^{(r)} se calcula resolviendo (4.4)-(4.6) usando Δy_{1},...,Δy_{n+m−2}, Δy_{n+m−1}^{(r)}.

Si la sucesión {Δy_{n+m−1}^{(r)}} fuera convergente a un número Δy_{n+m−1} > 0, resolveríamos (4.4)-(4.6) usando Δy_{1},...,Δy_{n+m−2}, Δy_{n+m−1}, lo cual nos daría {w_{i,n+1}}_{i=1}^{m+n}, y de esta manera quedaría resuelto el problema para el nivel n + 1.

Proposición 4.1 Sea A una matriz cuadrada tridiagonal, donde los elementos de la diagonal principal y las dos adyacentes son todos distintos de cero. Entonces A es una matriz irreducible, es decir, el sistema no puede desacoplarse.

Demostración: Ver [12].

□

Proposición 4.2 Si Δt > 0 y Δy_{1},...,Δy_{n+m−1} son números positivos, entonces la matriz del sistema lineal formado por las ecuaciones (4.4)-(4.6) es no singular.

Demostración: Claramente la matriz del sistema es diagonalmente dominante en sentido fuerte. Usando la proposición anterior se observa que esta matriz es también irreducible. Como un corolario del Segundo Teorema de Gershgorin resulta que la matriz es invertible.

□

Ahora despejaremos explícitamente la solución del sistema lineal (4.4)-(4.6) con el objeto de obtener propiedades de la solución.

De la ecuación (4.4) se deduce que:

w_{1,n+1} = w_{2,n+1} − Δy_{1}f_{n+1}.  (4.10)

\underline{Caso i = 2:} Reemplazamos por el valor de w_{1,n+1} y despejamos w_{2,n+1}, quedando:

w_{2,n+1} = \frac{−\frac{Δt}{Δy_{2}} f_{n+1} + w_{2n} + \frac{Δt}{Δy_{2}^{2}} w_{3,n+1}}{1 + \frac{Δt}{Δy_{2}^{2}}} .  (4.11)

Notemos que

w_{2,n+1} = a_{2,n+1} + b_{2,n+1}w_{3,n+1},  (4.12)

donde

a_{2,n+1} = \frac{\frac{−Δt}{Δy_{2}}f_{n+1} + w_{2n}}{1 + \frac{Δt}{Δy_{2}^{2}}} ,  b_{2,n+1} = \frac{\frac{Δy}{Δy_{2}^{2}}}{1 + \frac{Δt}{Δy_{2}^{2}}} ,  (4.13)

donde

a_{2,n+1} = a_{2,n+1} (Δt, Δy_{2}, f_{n+1}, w_{2n}) , b_{2,n+1} = b_{2,n+1} (Δt, Δy_{2}) , 0 < b_{2,n+1} < 1. (4.14)

\underline{Caso i = k, donde 3 ≤ k ≤ n + m − 2:} Reemplazamos por el valor de w_{k−1,n+1} y despejamos w_{k,n+1}, quedando:

w_{k,n+1} = \frac{\frac{Δt}{Δy_{k−1}Δy_{k}} a_{k−1,n+1} + w_{kn} + \frac{Δt}{Δy_{k}^{2}} w_{k+1,n+1}}{1 + \frac{Δt}{Δy_{k−1}Δy_{k}} (1 − b_{k−1,n+1}) + \frac{Δt}{Δy_{k}^{2}}} .  (4.15)

Notemos que:

w_{k,n+1} = a_{k,n+1} + b_{k,n+1}w_{k+1,n+1},  (4.16)

donde

a_{k,n+1} = \frac{\frac{Δt}{Δy_{k−1}Δy_{k}} a_{k−1,n+1} + w_{kn}}{1 + \frac{Δt}{Δy_{k−1}Δy_{k}} (1 − b_{k−1,n+1}) + \frac{Δt}{Δy_{k}^{2}}} ,  (4.17)

b_{k,n+1} = \frac{\frac{Δt}{Δy_{k}^{2}}}{1 + \frac{Δt}{Δy_{k−1}Δy_{k}} (1 − b_{k−1,n+1}) + \frac{Δt}{Δy_{k}^{2}}} ,  (4.18)

donde

a_{k,n+1} = a_{k,n+1} (Δt, Δy_{2}, ..., Δy_{k}, f_{n+1}, w_{2n}, ..., w_{kn}) ,  (4.19)
b_{k,n+1} = b_{k,n+1} (Δt, Δy_{2}, ..., Δy_{k}) ,  (4.20)
0 < b_{k,n+1} < 1.  (4.21)

\underline{Caso i = n + m − 1:} Siguiendo el mismo procedimiento de antes, queda:

w_{n+m−1,n+1} = \frac{\frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}} a_{n+m−2,n+1} + w_{n+m−1,n}}{1 + \frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}} (1 − b_{n+m−2,n+1}) + \frac{Δt}{Δy_{n+m−1}^{2}}} ,  (4.22)

donde

w_{n+m−1,n+1} = a_{n+m−1,n+1},  (4.23)

donde

a_{n+m−1,n+1} = \frac{\frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}} a_{n+m−2,n+1} + w_{n+m−1,n}}{1 + \frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}} (1 − b_{n+m−2,n+1}) + \frac{Δt}{Δy_{n+m−1}^{2}}} ,  (4.24)

donde

a_{n+m−1,n+1} = a_{n+m−1,n+1} (Δt, Δy_{2}, ..., Δy_{n+m−1}, f_{n+1}, w_{2n}, ..., w_{n+m−1,n}) .  (4.25)

Proposición 4.3 Los {b_{i,n+1}}_{i=2}^{n+m−2} no dependen de n. Demostración: Es una simple observación de (4.13), (4.14), (4.18), (4.20) para 3 ≤ k ≤ n + m − 2.

□

Usando la proposición anterior, y resumiendo todas las fórmulas, queda:

w_{1,n+1} = w_{2,n+1} − Δy_{1}f_{n+1},  (4.26)
w_{i,n+1} = a_{i,n+1} + b_{i}w_{i+1,n+1}, i = 2, ..., n + m − 2, (4.27)
w_{n+m−1,n+1} = a_{n+m−1,n+1}.  (4.28)

Donde los coeficientes están dados por:

a_{2,n+1} = \frac{\frac{−Δt}{Δy_{2}}f_{n+1} + w_{2n}}{1 + \frac{Δt}{Δy_{2}^{2}}} ,  (4.29)
a_{i,n+1} = \frac{\frac{Δt}{Δy_{i−1}Δy_{i}} a_{i−1,n+1} + w_{in}}{1 + \frac{Δt}{Δy_{i−1}Δy_{i}} (1 − b_{i−1}) + \frac{Δt}{Δy_{i}^{2}}} ,  i = 3, ..., n + m − 1, (4.30)
b_{2} = \frac{\frac{Δt}{Δy_{2}^{2}}}{1 + \frac{Δt}{Δy_{2}^{2}}} ,  (4.31)
b_{i} = \frac{\frac{Δt}{Δy_{i}^{2}}}{1 + \frac{Δt}{Δy_{i−1}Δy_{i}} (1 − b_{i−1}) + \frac{Δt}{Δy_{i}^{2}}} ,  i = 3, ..., n + m − 2.  (4.32)

Y los coeficientes son función de las siguientes variables:

a_{i,n+1} = a_{i,n+1} (Δt, Δy_{2}, ..., Δy_{i}, f_{n+1}, w_{2n}, ..., w_{in}) , i = 2, ..., n + m − 1, (4.33)
b_{i} = b_{i} (Δt, Δy_{2}, ..., Δy_{i}) , i = 2, ..., n + m − 2,  (4.34)
0 < b_{i} < 1, i = 2, ..., n + m − 2.  (4.35)

Proposición 4.4 Supongamos que el gradiente de presión es positivo, es decir, f(t) > 0 ∀ t, y que Δt, Δy_{1}..., Δy_{n+m−1} son números positivos, es decir, que el método está bien definido hasta el nivel n. Asumiendo que w_{in} ≤ 0, i = 1, ..., n + m − 1, entonces a_{i,n+1} < 0, i = 2, ..., n + m − 1.

Demostración: Lo que dice esta proposición es que si en el nivel n las incógnitas son no negativas, entonces los coeficientes resultan todos negativos, lo cual, como se verá más adelante, nos será útil para averiguar el signo de la solución en el nivel n + 1.

La demostración consiste en un sencillo paso inductivo. Si observamos la ecuación (4.29) es claro que a_{2,n+1} < 0. Si suponemos ahora que a_{i−1,n+1} < 0, con la ayuda de (4.30) y (4.35) se obtiene que a_{i,n+1} < 0, donde i = 3, ..., n + m − 1.

□

Proposición 4.5 Supongamos que el gradiente de presión es positivo, es decir, f(t) > 0 ∀ t, y que Δt, Δy_{1}..., Δy_{n+m−1} son números positivos, es decir, que el método está bien definido hasta el nivel n. Asumiendo que w_{in} ≤ 0, i = 1, ..., n + m − 1, entonces w_{i,n+1} < 0, i = 1, ..., n + m − 1.

Demostración: Se realiza un sencillo paso inductivo hacia atrás. Por la fórmula (4.28) y la Proposición (4.4) se observa que w_{n+m−1,n+1} < 0. Supongamos que w_{i+1,n+1} < 0, entonces, usando la ecuación (4.27) se tiene que

w_{i,n+1} = a_{i,n+1} + b_{i}w_{i+1,n+1} < 0.  (4.36)

Observación 2 Si comenzamos con la condición inicial sobre u0^{′} determinada por (B1), resulta que los w_{i,n+1} son todos números negativos, coincidiendo con el resultado teórico expresado en (3.2).

Proposición 4.6 Supongamos que el gradiente de presión es positivo, es decir, que f(t) > 0 ∀ t, que Δt, Δy_{1}..., Δy_{n+m−1} son números positivos, es decir, que el método está bien definido hasta el nivel n,que w_{in} ≤ 0, i = 1, ..., n + m − 1, y que y_{n+m−1} = 1 − (Δy_{1} + ... + Δy_{n+m−2}) es un número en el intervalo (0, 1). Entonces si comenzamos la iteración interna propuesta en (4.9) con un valor inicial Δy_{n+m−1}^{(1)} que cumpla 0 < Δy_{n+m−1}^{(1)} < y_{n+m−1}, se genera una sucesión en r, bien definida.

Demostración: Mediante el uso de la Proposición (4.5), y mediante la ecuación que define la iteración, es decir, la fórmula (4.9), resulta que 0 < Δy_{n+m−1}^{(r)} < y_{n+m−1} ∀ r.

□

Definición 4.7 Definimos los siguientes números:

A_{in} = \frac{w_{i,n} − w_{i+1,n}}{Δy_{i}} ,  i = 1, ..., n + m − 2, n ≥ 1,  (4.37)
B_{in} = \frac{w_{in} − w_{i,n−1}}{Δt} ,  i = 1, ..., n + m − 2, n ≥ 1.  (4.38)

Proposición 4.8 Supongamos que el gradiente de presión es positivo, es decir, que f(t) > 0 ∀ t, que Δt, Δy_{1}..., Δy_{n+m−1} son números positivos, es decir, que el método está bien definido hasta el nivel n, que w_{in} ≤ 0, i = 1, ..., n + m − 1, que y_{n+m−1} = 1 − (Δy_{1} + ... + Δy_{n+m−2}) es un número en el intervalo (0, 1). Asumiendo que A_{in} ≤ 0, i = 1, ..., n + m − 2, entonces A_{i,n+1} < 0, i = 1, ..., n + m − 1.

Demostración: Recordando que las w_{in} cumplían las ecuaciones (4.4)-(4.6) se observa que {A_{i,n+1}}_{i=1}^{n+m−1} satisfacen el siguiente sistema lineal:

A_{1,n+1} = −f_{n+1},  (4.39)

\frac{Δt}{Δy_{i}^{2}} A_{i−1,n+1} − (1 + \frac{Δt}{Δy_{i}^{2}} + \frac{Δt}{Δy_{i}Δy_{i+1}} ) A_{i,n+1} + \frac{Δt}{Δy_{i}Δy_{i+1}} A_{i+1,n+1} =
= −A_{in},  i = 2, ..., n + m − 2,  (4.40)

A_{n+m−1,n+1} = \frac{w_{n+m−1,n+1}}{Δy_{n+m−1}} .  (4.41)

Por hipótesis, A_{1,n+1} < 0. Además, por la Proposición (4.5) se ve que A_{n+m−1,n+1} < 0. Despejando explícitamente los {A_{i,n+1}}_{i=2}^{n+m−2}, se obtiene que:

A_{i,n+1} = Q_{i,n+1} + R_{i,n+1}A_{i+1,n+1},  i = 2, ..., n + m − 2,  (4.42)

donde

Q_{2,n+1} = \frac{\frac{−Δt}{Δy_{2}}f_{n+1} + A_{2n}}{1 + \frac{Δt}{Δy_{2}^{2}} + \frac{Δt}{Δy_{2}Δy_{3}}} , R_{2,n+1} = \frac{\frac{Δt}{Δy_{2}Δy_{3}}}{1 + \frac{Δt}{Δy_{2}^{2}} + \frac{Δt}{Δy_{2}Δy_{3}}} ,  (4.43)

Q_{i,n+1} = \frac{\frac{Δt}{Δy_{i}^{2}} Q_{i−1,n+1} + A_{in}}{1 + \frac{Δt}{Δy_{i}^{2}} (1 − R_{i−1,n+1}) + \frac{Δt}{Δy_{i}Δy_{i+1}}} ,  i = 3, ..., n + m − 2,  (4.44)

R_{i,n+1} = \frac{\frac{Δt}{Δy_{i}Δy_{i+1}}}{1 + \frac{Δt}{Δy_{i}^{2}} (1 − R_{i−1,n+1}) + \frac{Δt}{Δy_{i}Δy_{i+1}}} ,  i = 3, ..., n + m − 2.  (4.45)

Notar que Q_{2,n+1} < 0 y 0 < R_{2,n+1} < 1. Suponiendo que Q_{i−1,n+1} < 0 y 0 < R_{i−1,n+1} < 1, y mirando (4.44) y (4.45) es fácil ver que Q_{i,n+1} < 0 y 0 < R_{i,n+1} < 1. Luego, como estamos haciendo un paso inductivo, resulta Q_{i,n+1} < 0 y 0 < R_{i,n+1} < 1, i = 2, ..., n + m − 2.

Suponiendo ahora que A_{i+1,n+1} < 0, observando (4.42) y debido a lo dicho en el párrafo anterior, resulta que A_{i,n+1} < 0. Luego, como estamos frente a un paso inductivo, tenemos que A_{i,n+1} < 0, i = 2, ..., n + m − 2. Con esto queda demostrada la proposición.

□

Proposición 4.9 Supongamos que el gradiente de presión es positivo, es decir, que f(t) > 0 ∀ t, que f′(t) > 0 ∀ t, que Δt, Δy_{1}..., Δy_{n+m−1} son números positivos, es decir, que el método está bien definido hasta el nivel n, que w_{in} ≤ 0, i = 1, ..., n + m − 1, que y_{n+m−1} = 1 − (Δy_{1} + ... + Δy_{n+m−2}) es un número en el intervalo (0, 1). Asumiendo que B_{in} ≤ 0, i = 1, ..., n + m − 2, entonces B_{i,n+1} < 0, i = 1, ..., n + m − 1.

Demostración: Análoga a la demostración de la Proposición (4.8).

□

Observación 3 Hay que destacar que los números {A_{i,n+1}}_{i=1}^{n+m−1} son aproximaciones numéricas de w_{y}. Luego, el resultado de la Proposición (4.8) confirma los resultados teóricos, expresados en la Proposición (3.2), ecuación (3.5). Además, los números {B_{i,n+1}}_{i=1}^{n+m−1} son aproximaciones numéricas de w_{t}. Luego, el resultado de la Proposición (4.9) confirma los resultados teóricos, expresados en la Proposición (3.2), ecuación (3.6).

Hasta este punto, la iteración interna propuesta es muy cara computacionalmente, puesto que cada iteración interna involucra la resolución de un sistema lineal de tamaño considerable a medida que el tiempo transcurre. Lo que haremos es simplificar los cálculos para que eso se reduzca notablemente.

Supongamos que el algoritmo está definido hasta el tiempo t_{n}. Sabemos que:

Δy_{n+m−1}^{(r+1)} = (\frac{w_{n+m−1}^{(r)}}{w_{n+m−1}^{(r)} − τ_{0}} ) y_{n+m−1},  (4.46)

w_{n+m−1,n+1}^{(r+1)} = \frac{\frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}^{(r)}} a_{n+m−2,n+1} + w_{n+m−1,n}}{1 + \frac{Δt}{Δy_{n+m−2}Δy_{n+m−1}^{(r)}} (1 − b_{n+m−2}) + \frac{Δt}{Δy_{n+m−1}^{(r)}2}} .  (4.47)

Es claro que w_{n+m−1,n} = 0 ya que así se requirió en el nivel anterior o bien porque las condiciones iniciales son compatibles entre sí. Por simplicidad, renombramos algunas cosas y analizamos sus signos, quedando,

a = \frac{Δt a_{n+m−2,n+1}}{Δy_{n+m−2}} < 0,  c = \frac{Δt(1 − b_{n+m−2})}{Δy_{n+m−2}} > 0.  (4.48)

Δy^{(r)} = Δy_{n+m−1}^{(r)} > 0,  w^{(r)} = w_{n+m−1,n+1}^{(r)} < 0,  y = y_{n+m−1}.  (4.49)

Reemplazando (4.47) en (4.46) y haciendo uso de los símbolos escritos en (4.48) y (4.49) se llega a la siguiente expresión:

Δy^{(r+1)} = \frac{ay}{b − τ_{0}Δy^{(r)} − \frac{Δt}{Δy^{(r)}} τ_{0}} ,  b = a − cτ_{0}.  (4.50)

Estas últimas fórmulas nos dicen que el Δy es un punto fijo de la función

F (x) = \frac{ax}{b − τ_{0}x − \frac{Δt}{x}τ_{0}} ,  (4.51)

Torres-Turner, Bingham Unidimensional, MAT-Serie A, 5 (2001), 11-26 24

cuyo gráfico nos muestra que tiene a lo sumo dos puntos fijos, el x = 0 y/o un punto fijo positivo. Luego, el algoritmo se reformula, de tal forma que mediante el cálculo recursivo de algunos coeficientes se pueda trabajar directamente sobre la función F sin tener que resolver sistemas lineales. El número de iteraciones se reduce notablemente. Si ese punto fijo positivo es contractivo, es claro que satisface (4.8), con lo cual el método queda bien definido. Hay que destacar el hecho de que en numerosos experimentos numéricos el punto fijo resulta contractivo.

5 Resultados numéricos

El algoritmo fue programado en Matlab en su totalidad. En todos los ensayos que se realizaron se comprobó que los resultados numéricos coincidían con los resultados teóricos. Como muestra la Proposición (3.3), si graficáramos la frontera libre se observaría que ésta converge a cero de manera monótona cuando el tiempo tiende a infinito.

• Caso 1: s_{0} = 0.8, τ_{0} = 1, f(t) = 2 − \frac{1}{(t+1)^{3}} , v0^{′}(y) = 0.

• Caso 2: s_{0} = 0.8, τ_{0} = 1, f(t) = 2 + t, v0^{′}(y) = 0.

Referencias

[1] Comparini E., A One Dimensional Bingham Flow, Journal of Mathematical Analysis and its Applications, 1992, 127-139.

[2] Crank J., Free and Moving Boundary Problems, Claredon Press, Oxford (1984).

[3] Douglas J., Jr. & Gallie T. M., Jr., On the Numerical Integration of a Parabolic Equation Subject to a Moving Boundary Condition, Duke Math. J., 22 (1955), pp. 557-571.

[4] Duvaut G. & Lions J. L., Inequalities in Mechanics and Physics, vol. 219, Springer Verlag, 1976.

[5] Glowinky-Lions Tremolieres, Analyse Numerique des Inequalities Variationales, vol. 1,2, Dumomd, 1976.

[6] Gupta R. S. & Kumar D., Variable Time Step Methods for One-Dimensional Stefan Problem with Mixed Boundary Condition, Int. J. Heat Mass Transfer., vol 24, pp. 251-259, 1981.

[7] Murray W. D. & Landis F., Numerical and Machine Solutions of Transient Heat-Conduction Problems Involving Melting or Freezing, J. Heat Transfer. 81C, pp. 106-112, 1959.

[8] Primicerio M., Problemi di Diffusione a Frontiera Libera, Bolletino U.M.I. (5) 18-A (1981), 11-68.

[9] Rubinstein L.I., The Stefan Problem, Trans. Math. Monographs-vol. 27, Amer Math. Soc. ,Providence 1971.

[10] Tarzia D. A., Introducción a las Inecuaciones Variacionales Elípticas y sus Aplicaciones a los Problemas de Frontera Libre, CLAMI, 1981.

Torres-Turner, Bingham Unidimensional, MAT-Serie A, 5 (2001), 11-26 25

[11] Tarzia D., A Bibliography on Moving-Free Boundary Problems for the Heat-Diffusion Equation. The Stefan and Related Problems, MAT - Serie A, 2 (2000).

[12] Trefethen L. N. & Bau D., Numerical Linear Algebra, Society for Industrial and Applied Mathematics (SIAM), 1997, pp. 227-228.

Figura 1: Caso 1.

Torres-Turner, Bingham Unidimensional, MAT-Serie A, 5 (2001), 11-26 26 Figura 2: Caso 2.