Capítulo 7
Solución Numérica de la Ecuación de Flujo
Subterráneo
Teoría de Flujo Subterráneo
Semestre 2008-1
Alberto Rosas Medina
1
Índice
 Polinomios
de Lagrange
 Diferencias
Finitas en una Dimensión
 Diferencias
Finitas en dos Dimensiones
2
Teoría de Aproximación Polinomial
Motivación
El punto de partida es la definición de un polinomio.
Un polinomio es un expresión matemática que consiste de
una suma de potencias en una o más variables cada una
multiplicada por un coeficiente.
Ejemplo, un polinomio de una variable.
Polinomio de dos variables
3
Interpolación polinomial
Suponga que se dan n+1 puntos como
donde
son las absicas de los puntos (de la malla)
con separación arbitraría. Entonces el polinomio de orden
n que pasa por los n+1 puntos puede escribirse como una
serie de potencias
4

El ajuste de serie de potencias a los n+1 puntos da un
sistema de ecuaciones

Aunque los coeficientes ai pueden determinarse
resolviendo el sistema de ecuaciones. No es factible
debido a que primero se necesita un programa para
resolver ecuaciones lineales y segundo por que el orden
de los polinomios puede ser tan grande que induzca
errores de redondeo.
5
Polinomios de Lagrange



Los polinomios de Lagrange de mayor interés son los de
primero, segundo y tercero orden.
Ejemplo de Polinomio lineal
Supóngase que se tiene solo dos puntos en la malla,
entonces la manera más sencilla de ajustar un polinomio
es utilizando polinomios lineales, es decir,
6
Datos
propuesta de funciones lineales
Entonces se necesita plantear un polinomio lineal que será el
que ajuste los dos puntos utilizando
y
7
Polinomio lineal
Entonces el polinomio lineal es
Equivalente a la ecuación de la recta dados 2 puntos.
8
Polinomio de Orden n

Para ilustrar la aproximación con polinomios de orden n.
se considera n=2. Se tienen 3 puntos x0, x1, x2.

Se necesita tres polinomios de orden 2.
9

Entonces el polinomio que se ajusta a los tres datos es
10

Entonces de manera general se tiene

De manera abreviada es

Por lo tanto el polinomio de ajuste es
11

Un punto importante a considerar es el error de
aproximación, es decir, se tiene que

Donde
es el polinomio de aproximación
12
Método de Diferencias Finitas

El método de diferencias fintas representa derivadas
continuas para una ecuación diferencial parcial con
expresiones envolviendo la evaluación de la función
desconocida en puntos discretos.

La consecuencia natural de esta acción es la
transformación de un problema envolviendo derivadas
clásicas a uno envolviendo ecuaciones algebraicas.
13



El primer paso es representar f(x) usando
derivando esta expresión para aproximar la derivada
Se deriva el polinomio de aproximación
y
La derivada es solo en el polinomio, n es el grado del
polinomio de Lagrange. El término f(xj) es el valor
especifico de f(x)en el punto xj.
14

Si se considera un polinomio de orden n=2 y j=o para
esta expresión

Se obtiene
15

Al hacer la derivada se obtiene

Si se considera el punto xj , se obtiene

Usando similarmente para j=1,2 se tiene
16

El siguiente paso es saber donde se quiere la derivada
evaluada, sea x=x0 , entonces se tiene

Entonces la derivada de la función f en x0 está dada por

Donde el error es proporcional al cuadrado de
.
17

El procedimiento anterior se hace de manera semejante
para x1 y x2 . Lo importante es notar el resultado para x1.
Entonces evaluando x1 en

Se obtiene
18

Simplificando la ecuación anterior se obtiene

Nótese que en esta aproximación la información la
información en el nodo en el cual se está
aproximando no aparece en la fórmula.
19

Para encontrar la derivada de segundo orden

Si se considera j=0, y teniendo en cuenta que n=2
20

Desarrollando para j=1,2

Un punto importante a notar es que la segunda derivada
en el punto x1 tiene un error de aproximación cuadrático.
Esto es derivando 2 veces y evaluando x1 en la
siguiente expresión se obtiene el orden del error. Cosa
que no pasa con x0 y x2
21

Si se consideran polinomios lineales esto es

Se obtiene la primera derivada
22

Si se evalúa en x=x0 se obtiene

Nótese que debido a que se usaron sólo dos nodos, el
error de truncamiento es
y es mayor que
.
23

Para finalizar se plantea una manera que facilita las
expresiones de ecuaciones en diferencias finitas.
Ejemplo
24
Representación de la Ecuación de Flujo
Subterráneo con Diferencias Finitas

Ecuación de flujo en un medio poroso saturado

Usando la notación de la tabla anterior se tiene

Donde
25

Expandiendo y se define hik =h(xi , tk ) se tiene

Asociado a esta fórmula podemos resolver para h(xi , tk+1)
Véase la siguiente figura
26

Esquema computacional de Diferencias finitas. Para
hallar el valor un valor desconocido al tiempo k+1
27

Comparación de las ecuaciones

El término A corresponde al espacio en el nivel de tiempo
k, y el término B es la derivada en tiempo y está en
diferencia hacía adelante (forward), está localizado en i,k
y es proyectadoal nivel de tiempo k+1
28

Si se considera un tiempo inicial k=0, entonces en la
ecuación se modifican los subíndices y se tiene

Es decir, el sistema se resolverá para k+1
29
Problema bien planteado


Para resolver el sistema de ecuaciones se necesita que
el problema este bien planteado, es decir, que tenga
condiciones de iniciales y de frontera y la solución
dependa continuamente de estas y sea única.
Condiciones de Frontera Dirichlet
h(z,t)=f1(t),
z =0,L
En términos de diferencias finitas
h(zi,t)=hi i=0 y L
30

Condiciones Neumann: condiciones de flujo

Teniendo presente que

Entonces se tiene que

La representación de derivadas en los nodos finales
puede puede emplearse para diferencias centradas
forward y backward
31

Condiciones Robin

Simplificado se tiene

La ecuación anterior describe la salida o escape en un
acuífero cuando k es considerado como el coeficiente
de salida, h0 (z) es el valor de la carga en el exterior, y
h(z,t) es la carga en el acuífero.
32
Sistema de ecuaciones de las diferencias
finitas

Del sistema discretizado

Se puede escribir de manera abreviada como

Donde f son las condiciones de frontera, estas modifican
la matriz K en el primer y último renglón.
33

De manera general K está dada por

Y la matriz S por
34

Se tiene que resolver el siguiente sistema de ecuaciones

Se factoriza el término hk+1 y se tiene

Sea

Y se resuelve el sistema por un método efectivo excepto
hallando c-1.
entonces se tiene
35
Método de Diferencias Finitas en dos
Dimensiones

Se extiende el problema a dos dimensiones, para ello
considérese la ecuación

Tal que
Reescribiendo la ecuación del operador se tiene

36

La aproximación del operador en diferencias finitas es

Donde
Considere que el coeficiente de trasmisividad es menor
en xy entonces se anula el tercer sumando de la
ecuación

37

El esquema de discretización es
38
39
Ejemplo de la ecuación de flujo en dos
dimensiones
Para el caso de un acuífero isotrópico y homogenero:
Es decir:
Considerando el estado estacionario:
Lo que es igual a:
40
El método de diferencias finitas
En este método, el valor de la función desconocida en cada nodo es
aproximada por su desarrollo en series de Taylor
Esta ecuación se aplica en cada nodo
incógnita
41

Obtenemos un sistema Ax=b, de n ecuaciones con n incógnitas x y
una matriz A como la siguiente:
Donde el vector b contiene a las condiciones de frontera y no todas sus
entradas son cero.
Resolviendo el sistema obtenemos los valores de la función buscada
que es solución de la ecuación diferencial parcial
42
Experimento 1
Se tiene una porción de un acuífero confinado, homogéneo e
isotrópico, en el estado estacionario. Se quiere conocer la
distribución de la carga hidráulica h en toda la región. Esta se
encuentra limitada por un rectángulo como el siguiente:
lo que nos da las condiciones de frontera
condiciones de Dirichlet
43
Gráfica de la solución del experimento 1
Longitud de los lados:
XL=YL=100L.
La malla tiene 35x35
rectángulos y
36x36 nodos.
Solución analítica:
44
Experimento 2
Tenemos un acuífero como el del caso anterior. Las condiciones de
frontera en el problema se modifican de la siguiente manera:
Y en uno de los lados tenemos una condición de Neumann
45
Solución experimento 2
Longitud de los lados:
XL=4L YL=3L.
La malla tiene 35x35
rectángulos y 36x36
nodos.
Solución analítica es
igual que en el caso
anterior.
46
Experimento 3
Condiciones de frontera:
47
Experimento 4
Tenemos el cuadrado unitario en cuyo centro está ubicado el punto
(1/2,1/2). Esto representa una isla con un pozo en su centro.
El movimiento del agua que circula por debajo de esta región está
gobernado por la misma ecuación
48

Sabemos que:

Sustituyendo la velocidad de Darcy, y el area de la superficie del
cilindro

Despejando

E integrando

Esta ultima igualdad es solución de la ecuación diferencial parcial
que gobierna el fenómeno, así:
49

Nuestra condición de frontera es la siguiente:

Así

Donde:

Esta solución esta compuesta por dos funciones, una de ellas
no está determinada, y es precisamente la que encontraremos con
el algoritmo desarrollado

50
Solución
Esta es la gráfica de:
Aquí vemos la gráfica
encontrada:
51
Y finalmente la suma de las dos anteriores es la solución al problema
planteado
52
Aproximación en series de Taylor
53
54
Diferencias progresivas
55

Diferencias progresivas
56

Diferencias progresivas
57
Diferencias Regresivas
58
Diferencias Centrales
59
Diferencias Centrales
60
61
Resumen
62



Bibliografía
Referencia Principal, George F. Pinder y Michael A. Celia,
Subsurface Hydrology, Jhon wiley & Sons, Inc., Hoboken New
Jersey, 2006
“Desarrollo de un algoritmo para evaluar el suministro de agua
subterranea” Nora Isabel Pérez Quezadas. Tesis de Lic. en
Matemáticas, Facultad de Matemáticas. UV. México 2005.
http://mmc2.igeofcu.unam.mx/norapeq/
63
Descargar

Problemas en Flujo Bidimensional