A LA SIMULACION

INTRODUCCION

NUMERICA
DE PROBLEMAS

HIDRAULICOS

Angel
N. Men
endez
Jefe del Departamento de Modelos Matematicos y Estudios Especiales

INSTITUTO NACIONAL DE CIENCIA Y TECNICA

HIDRICAS
Laboratorio de Hidraulica Aplicada

Informe LHA - 064-003-87

Ezeiza, Septiembre de 1987

Resumen
Este trabajo es un apunte, donde se introducen los conceptos y metodologas basicas
para la resolucion numerica de ecuaciones diferenciales en derivadas parciales. Se
describen los distintos tipos de formulaciones analticas que admiten los problemas
de la mecanica del continuo, y, a partir de ellas, los metodos numericos mas utilizados
en la actualidad.

Indice general
1. Introducci
on
1.1. Simulacion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Formulacion continua y discretizacion . . . . . . . . . . . . . . . . . .
1.3. Programas computacionales . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5

2. Formulaci
on continua
2.1. Formulacion diferencial . . . . . . . . . . . . . . . . .
2.1.1. Ecuaciones diferenciales ordinarias . . . . . .
2.1.2. Ecuaciones diferenciales en derivadas parciales
2.2. Formulacion debil . . . . . . . . . . . . . . . . . . . .
2.3. Formulacion variacional . . . . . . . . . . . . . . . .
2.4. Formulacion integral . . . . . . . . . . . . . . . . . .
2.4.1. Formulacion directa . . . . . . . . . . . . . . .
2.4.2. Funcion de Green . . . . . . . . . . . . . . . .
2.4.3. Formulacion indirecta . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

7
7
7
9
11
13
15
15
18
19

3. Discretizaci
on
3.1. Diferencias finitas . . .
3.2. Residuos Ponderados .
3.3. Elementos finitos . . .
3.4. Elementos de contorno

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

21
21
26
28
32

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

Captulo 1
Introducci
on
1.1.

Simulaci
on

Por simulacion de un fenomeno fsico, se entiende su reproduccion bajo condiciones de control. Los objetivos de la simulacion son el estudio del fenomeno y/o la
prediccion de su evolucion en situaciones particulares de interes.
A partir de la observacion de un fenomeno fsico se pasa a un primer nivel de
abstraccion, que se denominara modelo teorico fsico. Este se basa en la conceptualizacion del fenomeno (identificacion del sistema, fuerzas interiores significativas,
accion externa), que surge a partir de los conocimientos brindados por la Fsica y el
conocimiento emprico del problema. Este modelo queda plasmado en un cuerpo de
hipotesis. A partir de el, es posible proceder a la simulacion fsica del fenomeno, es
decir, su reproduccion a escala (reducida o aumentada, dependiendo del problema).
En la jerga de la ingeniera hidraulica, a esto se lo denomina modelo fsico.
Una vez establecido el modelo teorico fsico, es posible avanzar a un segundo
nivel de abstraccion: el modelo teorico matematico. Este no es mas que la representacion simbolica del problema mediante las herramientas provistas por el Analisis
Matematico (ecuaciones diferenciales, integrales, etc.). Para dar este paso, es necesario recurrir a las leyes fsicas generales (por ejemplo, teoremas de conservacion)
que gobiernan el fenomeno. A partir del modelo teorico matematico, es posible proceder a la simulacion analogica, que consiste en estudiar un fenomeno distinto (por
ejemplo, un fenomeno electrico en el caso de utilizarse una computadora analogica),
pero regido por las mismas ecuaciones. Tambien es posible, en algunas situaciones
ya sea muy simplificadas o correspondientes a geometras muy simples), buscar la
solucion del problema matematico recurriendo a tecnicas analticas de resolucion.
Este procedimiento se denomina simulacion analtica.
En muchos problemas practicos la solucion por metodos analticos no es actualmente posible. En consecuencia, puede procederse a discretizar el modelo teorico
matematico, para transformar el problema analtico en uno algebraico, y as poder
recurrir al uso de la computadora digital. Esto se denomina simulacion numerica.
En la jerga de la ingeniera, a estos se los identifica como modelos matematicos. En
la practica, Cambien se incluye en esta denominacion a cualquier procedimiento de
calculo implementado en una computadora.
En el presente apunte, se desarrollaran con caracter introductorio los principales
3

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

aspectos de la simulacion numerica.

1.2.

Formulaci
on continua y discretizaci
on

La mayora de los problemas de la Hidraulica se formulan sobre la hipotesis de que

el sistema es continuo, es decir, sus propiedades (densidad, presion, velocidad, etc.) se
especifican a traves de funciones continuas. La aplicacion de las leyes de la Fsica (en
general, en la forma de teoremas de conservacion) a un tal modelo teorico requiere,
en general, la utilizacion de procesos de paso al lmite, que conducen a expresiones
diferenciales y/o integrales (que constituyen el modelo teorico matematico).
En primer lugar, se considera que el medio (agua, aire, etc.) es continuo. Esto,
denominado hipotesis del continuo, significa que se hace abstraccion de su estructura
microscopica, ya que las escalas de movimiento que interesan son mucho mayores que
las moleculares. Sin embargo, los movimientos microscopicos no estan desacoplados
de los macroscopicos, sino que los primeros adquieren su energa de los segundas, en
un proceso que se conoce como disipacion de la energa mecanica. Para expresar
matematicamente esta disipacion, es necesaria recurrir a leyes empricas. Las ecuaciones de Navier-Stokes surgen, precisamente, de incorporar relaciones empricas
entre el tensor de las tensiones (que engloba la influencia de las escalas moleculares) y la velocidad de deformacion del medio (que es una cantidad macroscopica),
introduciendo los coeficientes de viscosidad (de Newton y de volumen) [1].
El tema de la existencia de movimientos de escala muy diferente entre s, y el efecto de las escalas no resueltas sobre aquellas resueltas, no se agota en la division entre
movimientos microscopicos y macroscopicos. A
un dentro de lo macroscopico existen
escalas de movimiento dismiles, la cual conduce a que, en muchas situaciones, sea
aconsejable, e incluso necesario, proceder a nuevos filtrados . Por ejemplo, en el
caso de flujos turbulentos es com
un, desde Reynolds, descomponer a las variables en
valores medios y fluctuaciones. Los primeros son promedios estadsticos. Las fluctuaciones turbulentas estan alimentadas por el movimiento medio (en realidad, existe
una cascada de energa desde las escalas mayores a las menores). Si se aplica este
tratamiento a las ecuaciones de Navier-Stokes, se obtiene, en primer lugar, las ecuaciones de Reynolds para el movimiento medio [7]. En ellas, la influencia de las escalas
no resueltas (es decir, las escalas de las fluctuaciones turbulentas) se manifiesta a
traves de las tensiones de Reynolds. Para expresar a estas tensiones en terminos
de cantidades medias, se recurre a los denominados modelos de turbulencia , de
variados grados de complejidad [3].
Cuando los movimientos de escalas dismiles estan alimentados externamente, es
decir, no hay un trasvase de energa del uno al otra, cada uno de ellos puede estudiarse independientemente. Tal es lo que sucede, por ejemplo, entre los movimientos
generados por las olas y aquellos producidos por las mareas (en una primera aproximacion).
Dada la gran variedad de problemas a los que se enfrenta la Hidraulica, el modelo teorico matematico cambia, en general, de un estudio a otro. En efecto, si bien
la base camon esta constituida, en principio, por las ecuaciones de Navier-Stokes,
los necesarios filtrados y simplificaciones conducen a expresiones matematicas que,
en muchos casos, solo conservan un lejano parecido con aquellas. Esto significa que
4

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

las tecnicas de resolucion a las que debe apelarse son tambien variadas. De ah la
necesidad de adquirir, al menos, una vision panoramica de los diferentes metodos
disponibles, con sus alcances y limitaciones, que es lo que pretende el presente
apunte.
Tal cual se ha dicho en la seccion anterior, la resolucion de las ecuaciones del
modelo teorico matematico por tecnicas analticas solo es posible en muy contados
casos. Generalmente, se hace necesario recurrir a metodos numericos. Esto significa
que debe procederse a discretizar el problema, hasta aqu planteado en terminos
analticos. La discretizacion, que se realiza de acuerdo a tecnicas que se discutiran
en captulos siguientes, transforma el problema analtico en uno algebraico, el cual
puede resolverse utilizando la computadora.
El hecho es que la discretizacion introduce un error de truncamiento, denominado
error de discretizacion. Este error puede manifestarse en la solucion numerica como
si se tratara de fenomenos de difusion y/o dispersion. En algunos problemas, estos
fenomenos no estaban presentes en la formulacion continua original, lo cual da lugar
a imprecisiones. En otros, se produce una contradifusion que provoca la inestabilidad
de la solucion numerica. Precisamente, una de las claves de una adecuada simulacion
numerica esta en poder controlar los efectos del error de discretizacion.
La conviccion de que el proceso de discretizacion puede llevar, en algunos casos,
a variar sustancialmente el problema que se esta tratando de resolver, condujo a
desarrollar una metodologa especifica de discretizacion basada en consideraciones
fsicas. Se trata de interpretar el proceso de discretizacion como una fase mas del filtrado de escalas de movimiento. La idea es integrar las ecuaciones originales sobre las
celdas, elementos, segmentos, etc., que componen la grilla (la cual determina un sistema discretizado), y definir adecuadamente las variables medias. Los terminos que
no admiten una expresion directa en funcion de las variables medias, constituyen
los terminos de interaccion entre las escalas resueltas y las no resueltas (denominadas escalas sub-grilla). Para cerrar el planteo, se necesitan leyes para expresar
esos terminos de interaccion, lo cual no es un problema sencillo. En definitiva, los
resultados de las formulaciones basadas en este punto de vista, en poco difieren de
los obtenidos directamente a partir de la formulacion continua.

1.3.

Programas computacionales

El proceso de discretizacion convierte al problema analtico en uno algebraico,

apto para ser resuelto por una computadora digital. La resolucion requiere dise
nar
un procedimiento de calculo; este se denomina algoritmo. La descripcion detallada
de los pasos de calculo que constituyen un algoritmo es un programa, que puede,
o no, involucrar el uso de una computadora. Un programa, para ser implementado
en una computadora, debe ser escrito en uno de los lenguajes que la maquina es
capaz de interpretar. Generalmente, se programa en lenguajes de alto nivel (BASIC,
FORTRAN, PASCAL, etc.), que son aquellos adaptados a la forma expresiva del
usuario (palabras, smbolos matematicos, etc.). La computadora misma lo traduce
(mediante interpretes o compiladores) a lenguaje de maquina, y los convierte en
instrucciones.
Una computadora no vale solo por lo que es como equipo, es decir, por su hard5

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

ware , sino tambien por la calidad y cantidad de programas de que dispone, es decir,
por su software. Este incluye no solo el propio sistema operativo de la maquina,
absolutamente necesario para su operacion, sino tambien las libreras de funciones
matematicas, de graficacion, etc.
El programador (que, en la actualidad, puede ser practicamente cualquier profesional o tecnico), utilizando la base de hardware y software provista por la maquina,
desarrolla programas especficos para una determinada disciplina. Un paquete de
programas interconectados, dirigido a resolver una determinada familia de problemas, constituye un sistema computacional. Estos, en general, estan preparados para
ser amigables con el usuario. Es decir, la forma de alimentarlo de datos es clara y
simple, y la presentacion de resultados (tablas, graficos) concisa y significativa. En
este caso, los sistemas son productos comerciales y, a veces, altamente cotizados.
Un programa de calculo requiere, en general, que los datos se entren y almacenen
de una forma muy particular. Para lograrlo, partiendo de la informacion virgen de
que se dispone (mediciones, esquemas, etc.), es necesario un proceso de elaboracion
(depuracion, tabulacion, analisis estadstico, etc.). Este proceso es conveniente automatizarlo, dando lugar a los denominados programas de preprocesamiento. Estos,
en algunos casos, pueden involucrar el uso de perifericos especiales, tales como un
digitalizador, o de interfases, tales como un conversor analogico-digital.
Todo proceso de calculo requiere la eleccion de ciertos parametros numericos: intervalos de discretizacion, intervalos de tabulacion, etc. El criterio de eleccion de esos
parametros depende de consideraciones particulares al problema a resolver; basicamente, esta asociado a la escala de resolucion deseada y/o posible. Sin embargo,
ellos se hallan, en general, interrelacionados, es decir, la eleccion de alguno de ellos
determina el rango de eleccion de los restantes, de modo de conseguir la precision
requerida y, en algunos casos, de mantener la estabilidad del proceso de calculo.
Los resultados directos del calculo no son, en general, los mas significativos para
el usuario, sino que estos requieren a
un cierta elaboracion para proveer salidas u
tiles.
Por ejemplo, en muchos casos se obtienen valores nodales de una o mas variables; si se
trata de un dominio bidimensional, estos resultados adquieren mayor significacion, y
producen mayor provecho, si se representan vistas tridimensionales, curvas de nivel,
perfiles a lo largo de cortes especificados, etc. Estos procedimientos se implementan
en los denominados programas de posprocesamiento.

Captulo 2
Formulaci
on continua
2.1.

Formulaci
on diferencial

La forma mas com

un de expresar matematicamente las leyes fsicas que gobiernan un fenomeno, es a traves de ecuaciones diferenciales, es decir, expresiones que
involucran derivadas de las variables.

2.1.1.

Ecuaciones diferenciales ordinarias

Las ecuaciones diferenciales se denominan ordinarias cuando solo existe una variable independiente. Por ejemplo, si la funcion incognita es u, y esta depende solo de
la variable t, una ecuacion diferencial ordinaria de orden n para u(t) es una expresion
de la forma


nu
u 2 u
F t, u, , 2 , ..., n = 0
(2.1)
t t
t
Para determinar completamente a la funcion u(t), deben agregarse n condiciones
complementarias, en la forma de condiciones iniciales y/o de contorno (las cuales
pueden involucrar derivadas de hasta orden n1). Si el problema esta bien planteado, la solucion existe y es u
nica. Si bien hay teoremas de existencia y unicidad para
formas especificas de ecuaciones diferenciales, es, en general, la fsica del problema la mejor gua para determinar correctamente cuales deben ser las condiciones
complementarias a especificar.
Cuando las condiciones complementarias se especifican en un solo punto (condiciones iniciales), se esta frente a un problema de valores iniciales. En caso contrario,
se trata de un problema de valores de contorno. Los primeros se hallan asociados
a procesos evolutivos, mientras que los segundos corresponden a situaciones estacionarias (o, en particular, de equilibrio).
Un tema importante en relacion a las soluciones de ecuaciones diferenciales, es
el de la estabilidad. Esta esta asociada a la sensibilidad de la solucion frente a
cambios en las condiciones iniciales o de borde. Cuando cambios peque
nos en los
valores iniciales y/o de borde producen cambios tambien peque
nos en la solucion,
esta se dice estable, y viceversa. En problemas de valores iniciales, la estabilidad
requiere que la diferencia entre la solucion perturbada y la no perturbada tienda a
cero a medida que estas evolucionan. Si la diferencia tiende a crecer, la solucion es
7

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

inestable. En el caso intermedio en el que la diferencia permanece finita, se habla

de una solucion neutralmente estable. La inestabilidad de una solucion muestra la
existencia de un proceso fsico inestable, que, en general, evoluciona hacia un estado
que ya no puede ser descrito por la ecuacion diferencial resuelta (es decir, debe
redefinirse el modelo teorico).
A ttulo de ejemplo, sea el siguiente problema de valores iniciales:
du
+ au = 0
dt

u(t = 0) = u0

(2.2)

para t > 0, donde a y u, son valores constantes. La Ec. (2.2) puede transformarse
en
Z u
Z t
du
a dt
(2.3)
=
u0 u
0
Procediendo a la integracion, se obtiene
u(t) = u0 eat

(2.4)

La solucion (2.4) es decreciente o creciente, dependiendo de si a es positivo o

negativo, respectivamente. Para establecer la estabilidad de la solucion, supongase
que se perturba el valor inicial en la cantidad u. La solucion perturbada esta dada,
entonces, por una expresion identica a la (2.4), solo que debe reemplazarse u0 por
u0 + u. La diferencia entre ambas soluciones vale
u u = u eat

(2.5)

La Ec. (2.5) muestra que la solucion u es estable si a > 0, y viceversa (si a = 0

es neutralmente estable, pero este es un caso trivial).
El ejemplo anterior desarrollo una tecnica de resolucion directa de una ecuacion
diferencial ordinaria. Sin embargo, el tipo de solucion obtenida, de forma exponencial, es caracterstica de cualquier ecuacion diferencial lineal homogenea. Sea, por
ejemplo, el siguiente problema de valores de contorno:
du
d2 u
+u+x=0
;
u(0) = 0,
(1) = q
(2.6)
2
dx
dx
para 0 < x < 1. Esta ecuacion es lineal, pero inhomogenea. Su solucion general
es suma de una solucion particular y de la solucion general de la homogenea. Una
solucion particular de la Ec. (2.6) es
u(x) = x

(2.7)

(notese que esta no verifica las condiciones de contorno). La forma de la solucion de

la ecuacion homogenea es, de acuerdo a lo discutido,
u(x) = A ex

(2.8)

Introduciendo la Ec. (2.8) en la ecuacion homogenea se obtiene la siguiente

ecuacion caracterstica
2 + 1 = 0
(2.9)
8

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

que provee las dos soluciones 1 = +i y 2 = i. La solucion general de la homogenea

es una superposicion de las dos formas de solucion resultantes de introducir estos
dos valores de . En consecuencia, la solucion general de la inhomogenea es de la
forma
u(x) = x + a cos(x) + b sen(x)

(2.10)

donde a y b son constantes a ser determinadas en base a las condiciones de contorno.

De esta forma se obtiene la solucion final:
u(x) = (1 + q)

sen(x)
x
cos(1)

(2.11)

Es facil demostrar, siguiendo la metodologa anterior, que la solucion (2.11) es

neutralmente estable frente a variaciones del valor q.
Existen una variedad de tecnicas y procedimientos para obtener soluciones exactas de ecuaciones diferenciales ordinarias mas generales que las mostradas en los
ejemplos anteriores. Sin embargo, con ellos solo se pueden resolver una cantidad
limitada de problemas. El espectro se amplia si se recurre a metodos analticos para
obtener soluciones aproximadas. En este sentido, existen tecnicas de perturbacion y
analisis asintoticos que, en muchos casos, proveen resultados mas u
tiles que la propia
solucion exacta, al poner de manifiesto la importancia relativa de los mecanismos
fsicos que intervienen. Todos estos metodos pueden generalizarse, en principio, para
resolver sistemas de ecuaciones diferenciales ordinarias [2].

2.1.2.

Ecuaciones diferenciales en derivadas parciales

Cuando existe mas de una variable independiente (por ejemplo, el tiempo t y la

coordenada espacial x), las ecuaciones diferenciales involucran derivadas parciales
de las variables dependientes. Los problemas resultantes adquieren una complejidad,
en general, bastante mayor que los asociados a ecuaciones diferenciales ordinarias.
Las ecuaciones diferenciales en derivadas parciales se clasifican en tres tipos:
hiperbolicas, parabolicas y elpticas. Las dos primeras clases estan asociadas a problemas evolucionarios (es decir, a problemas puros de valores iniciales o mixtos de
valores iniciales y de contorno). Las ecuaciones hiperbolicas corresponden a problemas de propagacion de ondas, es decir, de transmision de informacion a velocidad
finita. Las ecuaciones parabolicas, en cambio, corresponden a un proceso de difusion,
es decir, transmision de informacion con velocidad infinita. Las ecuaciones elpticas,
por su parte, estan asociadas a problemas estacionarios (es decir, a problemas puros
de valores de contorno).
La determinacion del caracter de una ecuacion, se lleva a cabo utilizando la
metodologa de hallar las curvas caractersticas asociadas a dicha ecuacion. Desde el
punto de vista matematico, estas son las curvas en el plano x t (o, en general, las
hipersuperficies en el espacio de las variables independientes) a lo largo de las cuales
pueden ser discontinuas las derivadas de mayor orden de la ecuacion diferencial.
Fsicamente, las curvas pueden interpretarse como la trayectoria de propagacion de
la informacion. En rigor, ellas tienen existencia real solo en problemas evolucionarios.
En problemas elpticos, las curvas caractersticas resultan complejas, lo cual es el
9

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

criterio para determinar este caracter. En problemas parabolicos, algunas curvas son
reales y otras complejas, o son reales pero perpendiculares al eje de los tiempos (con
este criterio se diferencian los problemas parabolicos de los hiperbolicos, ya que,
entonces, la velocidad de propagacion resulta infinita).
Para ilustrar la aplicacion del metodo de obtencion de las curvas caractersticas,
se lo desarrollara para la siguiente ecuacion diferencial lineal de segundo orden (a la
cual pueden reducirse muchas ecuaciones de mayor orden):
a

2u
u
2u
u
2u
+
c
+
e
+ fu = g
+
b
+
d
t2
xt
x2
t
x

(2.12)

donde a, b, c, d, e, f y g son constantes. Si se define

v=

u
t

(2.13)

u
(2.14)
x
la Ec. 2.12 puede ser transformada en una ecuacion diferencial de primer orden
w=

v
w
v
+b
+c
+ dv + ew + f u = g
t
x
x

(2.15)

La variacion de las cantidades v y w (es decir, de las derivadas de u ) a lo largo

de un segmento infinitesimal en el plano x t esta dada por
dv =

v
v
dt +
dx
t
x

(2.16)

w
w
dt +
dx
(2.17)
t
x
Si se complementan las Ecs. 2.15 a 2.17 con la siguiente relacion entre v y w
dw =

v
w
=
x
t

(2.18)

queda determinado el siguiente sistema de ecuaciones diferenciales lineales de primer

orden

a b
0
c
v/t
dt dx 0

0 v/x

0 0 dt dx w/t
0 1 1 0
w/x

dv ew f u + g

dv

dw
0

(2.19)

Notese que las incognitas de las Ecs. 2.19 son las derivadas de mayor orden de
la ecuacion diferencial. Para que estas puedan ser discontinuas sobre ciertas trayectorias del plano x, t, dicho sistema debe de estar indeterminado a lo largo de ellas.
En particular, el determinante de la matriz de coeficientes debe anularse, lo cual
conduce a la condicion
 2
 
dx
dx
a
b
+c=0
(2.20)
dt
dt
10

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Las races de la Ec. 2.20, que definen la ecuacion para las curvas caractersticas
(en este caso particular son rectas), estan dadas por
p
b (b2 4ac)
dx
=
(2.21)
dt
2a
Si el discriminante es positivo, la Ec. 2.21 define dos familias de trayectorias
sobre el plano x t, que corresponden a sendas familias de curvas caractersticas, y
el problema es hiperbolico. En cambio, si el discriminante es negativo dichas trayectorias son complejas, lo cual muestra que el problema es elptico. En el caso de
transicion, es decir cuando el discriminante es nulo, las dos familias de trayectorias
coinciden (aunque apuntan en sentidos contrarios), y el problema es parabolico. En
este caso, si t es la coordenada temporal, consideraciones fsicas llevan a concluir
que las trayectorias deben ser normales al eje de los tiempos.
Dado que muchas de las formulaciones diferenciales de problemas fsicos parten
de leyes de conservacion, una de las formas basicas y mas comunes de ecuaciones
diferenciales en derivadas parciales (de tipo hiperbolico) es
u F (u)
+
=G
t
x

(2.22)

donde u puede ser un vector (que incluye masa, cantidad de movimiento, etc.), F (u)
es una densidad de flujo de la cantidad u, y G es un termino fuente (aporte de masa,
campo de fuerzas, etc.). Por ejemplo, la ecuacion de conservacion de la masa (o de
continuidad) para el escurrimiento de agua en un canal es
h (hu)
+
=q
t
x

(2.23)

donde h es el tirante, u la velocidad y q el aporte lateral por unidad de ancho.

Como en el caso de las ecuaciones diferenciales ordinarias, existen una variedad de
procedimientos para obtener soluciones exactas de ecuaciones en derivadas parciales,
muchos de ellos siendo extensiones de metodos desarrollados para las primeras. En
problemas hiperbolicos resulta de gran utilidad el metodo de las curvas caractersticas. Para las ecuaciones parabolicas y elpticas, una tecnica muy utilizada es separacion de variables, que transforma al problema en un conjunto de ecuaciones diferenciales ordinarias. Tambien pueden utilizarse transformaciones de similaridad o las
transformadas de Fourier, Laplace, etc. [11].

2.2.

Formulaci
on d
ebil

Existen formas alternativas de formular un problema en lugar de hacerlo en forma

diferencial. Estas formas alternativas permiten, como se vera, ampliar la categora
de soluciones posibles al problema en cuestion. Ademas, constituyen nuevos puntos
de partida para desarrollar metodos aproximados de resolucion.
Para ilustrar que se entiende por una formulacion debil, se partira del siguiente
problema diferencial elptico en el plano x y:
2 u + f (u, x, y) = 0
11

en

(2.24)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

u = u
en 1
(2.25)
u
= q
en 2
(2.26)
n
donde f es una funcion arbitraria, aunque conocida, u y q son datos, 1 + 2 = ,
siendo el contorno total del dominio y n es la normal exterior a . Notese que
la solucion u(x, y) a este problema debe ser continua hasta la segunda derivada, es
decir, u debe pertenecer a la clase C 2 .
Una forma de encarar la obtencion de una solucion aproximada al sistema 2.24 2.26, es requerir que esas ecuaciones se satisfagan en un sentido medio. Por ejemplo,
se puede pedir que

Z
Z 
Z
u
2
q w d = 0 (2.27)
( u + f (u, x, y)) w d +
(u u) w
d +
n

1
2
donde w, w y w
son funciones de peso (o de prueba) que distribuyen el error.
Esta claro que si la Ec. 2.27 debe satisfacerse para funciones de peso arbitrarias,
esta es equivalente a las Ecs. 2.24 - 2.26.
Si se integra por partes el primer termino de la Ec. 2.26, es decir, se aplica la
primera formula de Green:
Z
Z
Z

2
d
(2.28)
d = d +
n

donde y son funciones arbitrarias, se obtiene

Z
Z
Z
Z
u
u w d + f w d + w
d +
(u + u) w
d +
n

1

Z 
u
+ q w d = 0
n
2

(2.29)

Ahora pueden realizarse algunas simplificaciones: en primer lugar, se toma (sin

perder generalidad) w = w; en segundo lugar, se elige que w = 0 sobre 1 (lo
cual tampoco restringe la generalidad del planteo); finalmente, se impone que u = u
sobre 1 (aunque la solucion sea aproximada en el resto del dominio). Entonces, la
Ec. 2.29 se reduce a
Z
Z
Z
u w d + f w d +
q w d = 0
(2.30)

La Ec. 2.30, junto con la Ec. 2.25, es la formulacion debil del problema diferencial
2.24 - 2.26. Si la funcion de peso w es arbitraria (pero perteneciente a C 1 y haciendose
nula sobre 1 ) esta formulacion se reduce a la anterior. Sin embargo, ella es mas
amplia, en el sentido de que permite la existencia de soluciones con una restriccion
menor en cuanto a continuidad. En efecto, la Ec. 2.30 muestra que solo es necesario
que u pertenezca a C 1 es decir, puede tener derivada segunda discontinua. De existir,
estas se denominan soluciones debiles o generalizadas [15].
La consideracion de soluciones debiles es crucial en problemas hiperbolicos no
lineales. Sea, por ejemplo, un problema de valores iniciales definido por la Ec. 2.22,
con G = 0, es decir
u F (u)
+
=0
(2.31)
t
x
12

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

con condiciones iniciales dadas por

u(x, t = 0) = u0 (x)

(2.32)

Notese que este borde (t = 0) es del tipo 1 . Tal como esta formulado, la solucion
de este problema debe pertenecer a la clase C 0 (ya que las derivadas primeras pueden
ser discontinuas sobre las curvas caractersticas). Introduciendo, como antes, una
funcion de peso w, la Ec. 2.31 conduce a


Z Z
u F
dx
dt
+
w=0
(2.33)
t
x

0
Si se integra por partes la Ec. 2.33 se obtiene

 Z
Z Z
Z
w
w
t
dx u
dt
dt(F w)|x

dx(uv)|t=0 +
+F
+
x = 0 (2.34)
t
x

0
0

De acuerdo a las simplificaciones discutidas mas arriba, es conveniente tomar

w = 0 sobre el borde t = 0 (ya que es del tipo 1 ). Si, ademas, se admite que w
se anula fuera de un dominio finito (es decir, tiene soporte acotado), la Ec. 2.34 se
reduce a:


Z Z
w
w

dt
dx u
+F
=0
(2.35)
t
x
0

Si la Ec. 2.35 debe satisfacerse para cualquier funcion de prueba (perteneciente

a C 1 y sujeta a las restricciones ya explicitadas), constituye una formulacion alternativa a la Ec. 2.31 (en realidad, la Ec. 2.35 es la definicion de derivada en el
sentido de las distribuciones). Esta es la formulacion debil del problema que, como
en el caso anterior, es algo mas general que la formulacion diferencial, al permitir
soluciones discontinuas. Estas soluciones aparecen asociadas a fenomenos fsicos de
gran interes: ondas de choque en gases, ondas de frente abrupto en canales, etc. A
partir de la Ec. 2.35 es posible obtener explcitamente la relacion algebraica entre
los valores de la solucion inmediatamente hacia ambos lados de la discontinuidad.
Esta es
x[u]
= [F ]
(2.36)
donde x = dx/dt es la velocidad con que se desplaza la discontinuidad, y el
smbolo [F ] significa la diferencia entre los valores limites de la funcion F sobre
ambos lados de la discontinuidad. La Ec. 2.36 se denomina relacion de RankineHugoniot [15].

2.3.

Formulaci
on variacional

La formulacion variacional de un problema, como alternativa a la formulacion

diferencial, ya es clasica en la Mecanica. Sin embargo, ella no siempre existe.
El problema que plantea una formulacion variacional es encontrar una funcion
(o funciones) que hace estacionaria (es decir, maximiza o minimiza) una funcional
(o sistema de funcionales), y que este sujeta a ciertas condiciones de borde. Una
13

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

funcional es una funcion de funciones que, en la mayora de los problemas, representa

una cantidad fsica (por ejemplo, la energa).
El funcional, como en el caso de la formulacion debil, contiene derivadas de
un orden menor que el operador diferencial, por lo cual permite soluciones de una
clase mas amplia de funciones. En realidad, existe un parentesco entre la formulacion
debil y variacional. Por ejemplo, supongase que la funcion de peso w es una variacion
(debil) de u, y denotesela, entonces, como u. Si, ademas, se supone que en la Ec.
2.30 la funcion f (u, x, y) es de la forma f = u + g(x, y), donde g es una funcion
arbitraria, dicha ecuacion puede escribirse como
Z
Z
Z
u u d + (u + g) u d +
q u d = 0
(2.37)

que es equivalente a:

Z
Z 
1
1
2
2
(u) + (u ) + (gu) d +
(
q u) d = 0
2
2

2
Si se define el funcional

Z 
Z
(u)2 u2
I(u) =

gu d
(
q u) d = 0
2
2

(2.38)

(2.39)

la Ec. 2.38 es equivalente a requerir que

I = 0

(2.40)

es decir, que el funcional sea estacionario. Esta es la formulacion variacional del

problema en cuestion [5].
Recuerdese que, para que la Ec. 2.37 sea valida, se requiere que u = 0 sobre
1 . Esta se denomina condicion de borde geometrica, ya que debe ser impuesta
explcitamente. En cambio, la condicion de contorno sobre 2 se ha de cumplir
automaticamente, por lo cual se denomina condicion de borde natural.
En el caso particular de una dimension espacial, y tomando g(x) = x, x = 0 como
el contorno 1 , y x = 1 como el contorno 2 , se tiene el problema cuya formulacion
diferencial se presenta en la Ec. 2.6. De acuerdo a la Ec. 2.39, el funcional asociado
es
!
 2
Z 1
u2
1 du

xu dx q u(1)
(2.41)
I(u) =
2 dx
2
0
Para verificar que la formulacion variacional provee el mismo resultado que la
diferencial, se tomara una solucion de la forma
u(x) = a sen(x) + bx

(2.42)

donde las constantes a y b se ajustaran de manera de hacer estacionario el funcional.

Reemplazando la Ec. 2.42 en la Ec. 2.41 y operando, se obtiene


sen(2)
b2
1
I(u) =
a + + cos(1) ab + (cos(1) (1 + q)sen(1))a q +
b (2.43)
2
3
3
14

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Para extremar la funcional 2.45 se impone que I/a = I/b = 0, lo cual

conduce al siguiente sistema de ecuaciones:


sen(2)/2 cos(1)
cos(1)
2/3



a
b


=

cos(1) + (1 + q)sen(1)
q + 1/3


(2.44)

cuya solucion es
a = (1 + q)/cos(1)
b = 1

(2.45)

la cual, reemplazada en la 2.42, da la solucion 2.11, tal como se quera mostrar.

2.4.

Formulaci
on integral

La formulacion integral de un problema surge de avanzar un paso mas en la

integracion por partes que condujo a la formulacion debil. Aplicando, entonces,
nuevamente la primera formula de Green, Ec. 2.28 a la Ec. 2.29 se obtiene
Z
Z
Z
Z
w
u
2
d + w
d +
u w d + f w d u
n
n


Z
Z 
u
(u u) w
d +
q w d = 0
(2.46)
n
1
2
Si se elige, sin perdida de generalidad, que w = w y w
= u/n la Ec. 2.44 se
reduce a
Z
Z
Z
Z
u
2
w
u w d + f w d
d +
w
q d
n
1

2
Z
Z
w
w
d
u
d = 0
(2.47)
u
1 n
2 n
Notese que la restriccion sobre la categora de la solucion u es a
un mas suave que
antes: solo es necesario que sea continua (en cambio, w debe pertenecer a C 2 ).
A partir de la Ec. 2.45 se generan dos clases distintas de formulaciones integrales:
la directa y la indirecta. En lo que sigue se analizara el caso particular en que
f (u, x, y) = u + g(x, y)

(2.48)

donde u es una constante y g una funcion arbitraria.

2.4.1.

Formulaci
on directa

Cuando la funcion f esta dada por la Ec. 2.46, las integrales de area de la Ec.
2.45 pueden reescribirse como
Z
Z
2
u( w + w) d + gw d
(2.49)

15

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Sea, ahora, la siguiente ecuacion diferencial:

2 u + u + (x xp , y yp ) = 0

(2.50)

donde es la funcion generalizada delta/de/Dirac y xp , yp , las coordenadas, fijas

pero arbitrarias, de un punto P . La solucion u (x, y; xp , yp ) de la Ec. 2.48 sobre
un dominio infinito (con u 0 en el infinito) se denomina solucion fundamental
del operador diferencial ( 2 + ), y representa el efecto de una fuente (de intensidad unitaria) concentrada en el punto P . La solucion fundamental para el presente
problema es
1 (2)
H ( r)
(2.51)
u (x, y; xp , yp ) =
4i 0
donde r = [(x xp )2 + (y yp )2 ]1/2 .
Si se elige w = u , el primer termino de la Ec. 2.47 da

Z
u(xp , yp ) si P es interior a
2
u( w + w)d =
(2.52)
0
si P es exterior a

Si el punto P cae sobre la frontera es necesario desarrollar un tratamiento

especial para evaluar su contribucion, ya que la solucion fundamental es singular
para r = 0. El tratamiento, clasico en este contexto, consiste en evitar el punto
P , deformando la frontera en su derredor mediante una semicircunferencia e de
radio , que luego se hace tender a cero. Si el punto P permanece en el interior,
tal como se muestra en la Fig. 2.1, la integral de area da, de acuerdo a la Ec. 2.52,
una contribucion u(xp , yp ). Pero tambien hay que considerar, en la Ec. 2.47, los
terminos de la forma
Z
u
d
(2.53)
u
n
e
En efecto, para la solucion fundamental dada en la Ec. 2.51 el termino dominante
de su expansion es
du
1
u
=

(2.54)
n
dr
2r
por lo cual la contribucion 2.53 puede evaluarse como


Z
u
1
1
u
d u(xp , yp )
 = u(xp , yp )
(2.55)
n
2
2
e
que es exacta en el lmite  0. Es facil demostrar, siguiendo la misma metodologa,
que los terminos de la forma
Z
u
u d
(2.56)
n
e
dan una contribucion nula en ese limite.
Entonces, si se toma al punto P sobre el contorno , la Ec. 2.47, de acuerdo a
las Ecs. 2.49, 2.52 y 2.55, se convierte en
Z
Z
1
u
u
u(xp , yp ) +
u
d
u
d =
2
n
n
1
2
Z
Z
Z
u

gu d
u qd +
u
d
(2.57)
n

2
1
16

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Figura 2.1: Artificio para calcular la contribucion cuando el punto P cae sobre la
frontera
La Ec. 2.57 es una ecuacion integral, cuyas incognitas son los valores de u sobre
2 y los de u/n sobre 1 . Los terminos del miembro de la derecha no contienen
incognitas, por lo cual pueden evaluarse (en forma exacta o aproximada).
Una vez resuelta la Ec. 2.57, se conocen u y u/n completamente sobre todo
el contorno . Para obtener la solucion en el interior de , puede tomarse el punto
P en esa zona, con lo cual la Ec. 2.47 conduce a
Z
Z
Z
u

u
d u
d
(2.58)
u(xp , yp ) =
gu d + u
n
n

que es una formula explicita, ya que los terminos de la derecha no contienen incognitas. Se completa as el planteo de la formulacion integral directa [4].
A titulo de ejemplo, considerese el caso unidimensional ya discutido en la Seccion
2.4. La solucion fundamental, teniendo en cuenta que = 1, es
1
(2.59)
u(x; xp ) = sen|r|
2
donde r = x xp . La ecuacion correspondiente a la Ec. 2.57, cuando P coincide con
el borde x = 0, es



Z 1


1
u 
u 

u 
u(0) + u
u
=
xu dx u |x=1 q + u(0)u
(2.60)


2
n x=0
n x=1
n x=0
0

17

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Dado que

u(0) = 0

n x=0
x x=0




=

n
x x=1

x=1

u = sen(x)

u = 1 cos(x)
x
2
la Ec. 2.60 da
u(x = 1) =

sen(1)
(1 + q) 1
cos(1)

(2.61)

(2.62)

que, logicamente, coincide con lo calculado por medio de la Ec. (2.11).

En forma totalmente similar puede plantearse la Ec. 2.57 para el caso en que P
coincida con el borde x = 1, de donde surge el valor de u/x en x = 0. Finalmente,
con 0 < xp < 1, puede utilizarse la Ec. 2.58 para obtener la solucion completa.

2.4.2.

Funci
on de Green

Una variante de la formulacion integral directa es el metodo de la funcion de

Green. La funcion de Green, o de influencia, esta emparentada a la solucion fundamental, pero, a diferencia de esta, debe satisfacer ciertas condiciones de borde sobre
.
Supongase, para simplificar, que el borde es enteramente del tipo 1 (es decir,
se trata de un problema con condiciones de borde del tipo de Dirichlet). La funcion
de Green G(x, y; xp , yp ) asociada a este problema es la que satisface la Ec. 2.50 y la
condicion de borde
G=0
(2.63)
si (x, y) esta sobre Notese que si se conoce la funcion de Green para el problema
en cuestion, puede utilizarse directamente la Ec. 2.57 para calcular la solucion, es
decir
Z
Z
G
d
(2.64)
u(xp , yp ) =
gGd u

x
Obviamente, la dificultad se traslada a obtener la funcion de Green para un
dado problema, lo cual requiere, precisamente, resolver ecuaciones integrales del
tipo de la Ec. 2.57. La ventaja del metodo, entonces, solo se manifiesta cuando es
necesario resolver una familia de problemas con contornos fijos, ya que, en este caso,
la ecuacion integral debe resolverse una sola vez.

18

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

2.4.3.

Formulaci
on indirecta

La formulacion integral indirecta, como se vera, no resulta en una ventaja operativa. Sin embargo, conduce a una interpretacion del proceso de calculo conceptualmente mas rica.
La idea basica consiste en plantear en el exterior de (que se denominara 0 )
un problema simplificado, pero que produce, en el interior de el mismo efecto
que el problema fsico real. Especficamente, se considera que esa zona esta libre de
fuentes, es decir, se plantea que la solucion u0 en 0 satisface
2 u0 + u0 = 0

(2.65)

y que existe una distribucion de elementos (fuentes o dipolos) sobre que da cuenta
del efecto neto sobre la region .
Obviamente, es posible plantear el problema en 0 en forma integral. Si se elige
el punto P en el exterior de 0 , es decir, en el interior de ), la Ec. 2.47, con w = u ,
conduce a
Z
Z

0
0 u
u
d

u
d = 0
(2.66)
u
n0
n0

donde se ha supuesto que u0 0 en el infinito, y donde /n0 = /n. A partir de

aqu existen dos alternativas posibles. La primera, que se desarrollara en cierto detalle, consiste en suponer que u0 = u sobre , pero que u/n puede ser discontinuo,
lo cual equivale a admitir que hay una distribucion de fuentes sobre .
Si el pun to P esta sobre , puede utilizarse la metodologa ya descripta para
circunvalar este punto. Si se hace como en la Fig.2.1, la integral de area de la Ec.
2.47 a
un se anula, pero hay una contribucion extra dada por
Z
Z

w
1
0 u
de =
u0
de = u0 (xp , yp )
(2.67)
u
0
n
r
2
e
e
con lo cual la Ec. 2.47 da
1
u0 (xp , yp ) +
2

Z
u

n0

Z
d

0 u

n0

d = 0

Entonces, sumando las Ecs. 2.57 y 2.68 se obtiene

Z
Z

u(xp , yp ) + u d = gu d

(2.68)

(2.69)

donde = u/n + u0 /n0 = u/n u0 /n la intensidad, por unidad de area,

de la distribucion de fuentes. Diferenciando la Ec. 2.69 en la direccion de la normal
a se llega a
Z
Z
u(xp , yp )
u
1
u
+
d (xp , yp ) = g
d = 0
(2.70)
np
2
np
np
Las Ecs. 2.69 y 2.70 son ecuaciones integrales para la funcion incognita . La
primera se utiliza cuando el punto P esta sobre (ya que, entonces, u(xp , yp ) es
dato), y la segunda cuando P esta sobre 2 . Una vez conocida , la solucion en
19

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

se calcula de acuerdo a la formula que surge de sumar miembro a miembro las Ecs.
2.48 y 2.66, es decir
Z
Z

u(xp , yp ) =
gu d + u d
(2.71)

Queda as completada la formulacion indirecta con fuentes [4].

La segunda alternativa es la formulacion indirecta con dipolos. En este caso se
supone que u0 /n = u/n sobre , pero que u puede ser discontinuo. La ecuacion
equivalente a la Ec. 2.67 es
Z
Z
1
u
1
d = gu d
(2.72)
u(xp , yp ) +
2
2
n

donde = u0 u es la intensidad, por unidad de area, de la distribucion de dipolos.

Para obtener la ecuacion equivalente a la Ec. 2.68 hay que diferenciar 2.72. Esto
conduce a una expresion mas complicada desde el punto de vista operativo, ya que
aparece la derivada segunda de u . Es por eso que la formulacion con fuentes es la
mas utilizada en la practica.

20

Captulo 3
Discretizaci
on
3.1.

Diferencias finitas

El metodo de las diferencias finitas es una tecnica de resolucion aproximada que

parte de la formulacion diferencial del problema. La idea basica del metodo consiste
en aproximar las derivadas por cocientes en diferencias, que involucran valores de la
solucion en una serie de puntos del espacio de definicion del problema denominados
nodos ([6]; [14]; [9] y [10]).
Por ejemplo, sea un problema unidimensional donde debe aproximarse la derivada du/dx. La posibilidad mas simple consiste en expresarla de la siguiente manera:

uj+1 uj
du 
=
(3.1)

dx x= xj+1 xj
donde el subndice identifica el nodo. Ahora bien, la precision de la aproximacion
expresada en la Ec. 3.1 depende de cual es el punto g en el cual se quiere aproximar
la derivada, en relacion con los puntos xj y xj+1 . Si coincide con xj la Ec. (3.1)
es una aproximacion en adelanto para du/dx, y al cual se observa en la Fig. 3.1, es
menos precisa, en general, que en el caso en que y es interior al intervalo (xj , xj+1 ).
Lo mismo sucede cuando coincide con xj+1 que constituye una aproximacion en
atraso.
Para determinar el orden de precision de la aproximacion, se desarrolla la funcion
en serie de Taylor alrededor de f , es decir,


1 d2 u 
du 
(xj ) +
(xj )2 +
uj = u() +
dx x=
2 dx2 x=

1 d3 u 
+
(xj )3 + O(4)
(3.2)

3
6 dx x=

uj+1



du 
1 d2 u 
= u() +
(xj+1 ) +
(xj+1 )2 +


2
dx x=
2 dx x=

3 
1 d u
+
(xj+1 )3 + O(4)
6 dx3 
x=

21

(3.3)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Figura 3.1: Aproximacion de una derivada por diferencias finitas

donde O(4) significa terminos de orden cuatro en el respectivo incremento. Reemplazando las Ecs. 3.2 y 3.3 en la expresion en diferencias dada en la Ec. 3.1 se
obtiene

 2 

du 
d u 
uj+1 uj
xj+1 xj
=

+
+

xj+1 xj
dx x=
2
dx2 x=

 3 
1 x2j+1 + xj+1 xj + x2j
d u 
2
+ (xj+1 xj ) +
+
+ ...
(3.4)
2
3
dx3 
x=

Cuando = xj , la Ec. 3.1 se reduce a




uj+1 uj
du 
h d2 u 
h2 d3 u 
=
+
+
+ ...
xj+1 xj
dx x=xj 2 dx2 x=xj
6 dx3 x=xj

(3.5)

donde h = xj+1 xj . Es decir, la aproximacion en adelanto tiene un orden 1 de

precision en el incremento h, ya que ese es el mayor orden del error de truncamiento.
La aproximacion en atraso tiene el mismo orden de precision. En cambio, cuando
= (xj + xj+1 )/2 , que se denomina aproximacion centrada, se tiene que


uj+1 uj
du 
h3 d3 u 
=
+
+ ...
(3.6)
xj+1 xj
dx x= 24 dx3 x=
es decir, el orden de precision es 2.
Para obtener un orden de precision mayor que 2 en una derivada primera, es
necesario utilizar mas de dos nodos. Para aproximar una derivada segunda deben
22

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

utilizarse no menos de tres nodos. La formula que relaciona el orden de la derivada (OD), el orden de precision de la aproximacion (OP ) y el n
umero de nodos
involucrados (N N ) es
OP = N N OD
(3.7)
salvo en el caso de tratarse de una aproximacion centrada y que la paridad de N N
y OD sean distintas, en cuyo caso se gana un orden de precision (por ejemplo, una
aproximacion centrada con 5 nodos para una derivada segunda tiene un orden de
precision 4, en lugar de 3).
Si bien existen diversas maneras de construir una aproximacion en diferencias
para una derivada, la tecnica mas sistematica es el conocido metodo de los coeficientes indeterminados [9].
Lo dicho hasta ahora apunta a construir una aproximacion en diferencias para
una derivada. Pero este es solo el primer escalon en el establecimiento de una aproximacion en diferencias para una ecuacion diferencial. En efecto, no cualquier combinacion constituye un esquema numerico apto para el calculo.
Para ilustrar el procedimiento de discretizacion y calculo se tomara, en primer
lugar, el siguiente problema de valores iniciales:
du
+ f (u, t) = 0
dt

u(t = 0) = u0

(3.8)

Notese que se trata de una ecuacion diferencial ordinaria de primer orden, siendo f
una funcion arbitraria y u* una constante. La forma mas sencilla de discretizar la
Ec. 3.8 es mediante el metodo de Euler:
uj+1 uj
+ f (un , tn ) = 0
tn+1 tn

(3.9)

donde el supraindice n identifica el nodo en la dimension temporal t. Notese que,

respecto del punto de evaluacion de la funcion f , la aproximacion de la derivada es
en adelanto. La. Ec. 3.9 puede reescribirse como
un+1 = un + (tn+1 tn )f n

(3.10)

donde f n = f (tn , un ). La Ec. 3.10 permite evaluar en forma explcita la solucion en

el nivel temporal n + 1, en funcion de su valor en el nivel n. Por eso es que se dice
que el metodo de Euler pertenece a la categora de los metodos explcitos. La Ec.
3.10 es la formula operativa del metodo, es decir, partiendo de u0 = u(t = 0) = u0
dato, permite calcular, mediante un procedimiento de marcha, la solucion un , para
n 1.
Por ejemplo, en el caso particular en que f (u, t) = u y los nodos estan equiespaciados, es decir, tn+1 tn = k = cte. para todo valor de n, la Ec. 3.10 se reduce
a
un+1 = (1 k)un
(3.11)
Si u0 = 1, la solucion exacta del problema diferencial, Ec. 3.8, es, de acuerdo a
la Ec. 2.4,
U (t) = et
(3.12)
23

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

n
0
1
2
3
4
5
6
7
8
9
10

t
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0

U (tn )
1
0.9048
0.8187
0.7408
0.6703
0.6065
0.5488
0.4966
0.4493
0.4066
0.3679

un (Euler)
1
0.9000
0.8100
0.7290
0.6561
0.5905
0.5314
0.4783
0.4305
0.3874
0.3487

un (Crank N icolson)
1
0.9048
0.8186
0.7406
0.6701
0.6063
0.5485
0.4963
0.4490
0.4063
0.3676

Tabla 3.1: Resultados numericos para el problema de la Ec. 3.8

Tomando k = 0,1, la Tabla 3.1 muestra los valores de la solucion numerica
obtenidos mediante la Ec. 3.11, y los correspondientes valores de la solucion analtica,
dados por la Ec. 3.12, para los primeros pasos de calculo. All se observa que el
metodo de calculo es adecuado, aunque no muy preciso.
Una alternativa posible al metodo de Euler es el esquema fuertemente implcito,
que, aplicado a la Ec.3.8, da
uj+1 uj
+ f n+1 = 0
n+1
n
t
t

(3.13)

es decir, ahora la aproximacion de la derivada es en atraso. Notese que la incognita,

u es uno de los argumentos de f , por lo cual la Ec. 3.13 es una formula implcita
para un+1 , de ah el nombre del metodo. Salvo en el caso en que f sea lineal en
u, la obtencion de un+1 debe realizarse por medio de tecnicas iterativas (secante,
Newton-Raphson, etc.).
Otra alternativa, muy conocida, es el metodo de Crank-Nicolson, que da
uj+1 uj 1 n
+ (f + f n+1 ) = 0
n+1
n
t
t
2

(3.14)

Notese que se trata de un metodo implcito, y que la aproximacion esta centrada

en el punto intermedio entre tn y tn+1 . Para el caso particular del problema f = u,
u0 = 1, los valores de la solucion numerica tambien se muestran en la Tabla 3.1.
Observese que estos son mucho mas precisos que los obtenidos por el metodo de
Euler. Es que este metodo solo tiene un orden de precision 1 (es decir, el error de
truncamiento o discretizacion es de orden k), mientras que el Crank-Nicolson, al ser
centrado, tiene un orden de precision 2 (error del orden de k 2 ). Esto puede comprobarse reemplazando la solucion analtica en el esquema numerico, y efectuando
desarrollos en serie de Taylor alrededor del nivel temporal n. Para el metodo de
Euler se obtiene que el error de discretizacion es
E =

k d2 u k 2 d3 u
+
+ O(3)
2 dt2
6 dt3
24

(3.15)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

M etodo
Solucion Analtica Ec. 3.8
Diferencias Finitas Ec. 3.21
Residuos Ponderados, momentos Ec. 3.34
Residuos Ponderados, Galerkin Ec. 3.36
Elementos Finitos Ec. 3.57
Elementos de Contorno Ec. 3.80

u1 = u(1/3)
0.05550
0.05609
0.05433
0.05540
0.05494
0.05550

u2 = u(2/3)
0.06820
0.06891
0.05E88
0.06805
0.06751
0.06820

Tabla 3.2: Resultados numericos para el problema de la Ec. 3.17

mientras que para el Crank-Nicolson se tiene




k d2 u df
k 2 1 d3 u 1 d2 f
CN =
+
+
+
+ O(3)
2 dt2
dt
2 3 dt3
2 dt2

(3.16)

Pero el factor entre parentesis del primer termino del segundo miembro de la Ec.
3.16 se anula, ya que u satisface la Ec. 3.8, lo cual muestra que el metodo es de
segundo orden de precision. Sea ahora el siguiente problema de valores de contorno:
d2 u
+u+x=0
dx2

u(0) = u(1) = 0

(3.17)

cuya solucion analtica es

senx
x
(3.18)
sen1
A fin de ilustrar el procedimiento, divdase al intervalo (0, 1) en tres tramos
iguales de longitud h = 1/3. Ahora se trata de calcular la solucion en los dos nudos
interiores. La siguiente es una aproximacion centrada para la Ec. 3.17:
u(x) =

uj+1 2uj + uj1

+ uj + jh = 0
h2

j = 1, 2

(3.19)

que puede reescribirse como

uj+1 (2 h2 )uj + uj1 = jh3

(3.20)

Sabiendo que u0 = u3 = 0, la Ec. 3.20 conduce al siguiente sistema de ecuaciones

lineales: que puede reescribirse como


 

17/9
1
u1
1/27
=
(3.21)
1
17/9
u2
2/27
cuya solucion se presenta en la Tabla 3.2. Comparando con los valores que surgen
de la solucion analtica, Ec. 3.18, que tambien se muestran en esa tabla, se observa
que el metodo de calculo es adecuado a pesar del escaso numero de nodos utilizados.
Como u
ltimo ejemplo, supongase que se quiere resolver el siguiente problema
elptico:
2u 2u
+
+ f (u, x, y) = 0
en :
(3.22)
x2 y 2
25

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

donde es un dominio definido sobre el plano x y, con frontera sobre la cual se

especifican condiciones de borde del tipo de Dirichlet. La Ec. 3.22 se puede discretizar
con un esquema centrado como sigue:
ui+1j 2uij + ui1j uij+1 2uij + uij1
+
+ f (uij , x, y) = 0
x2
y 2

(3.23)

donde se ha supuesto que la grilla de discretizacion es rectangular, con pasos x y

y en las direcciones x e y, respectivamente, y donde los indices i y j identifican
los nodos a lo largo de los ejes x e y, respectivamente. Utilizando la Ec. 3.23 y las
condiciones de borde, puede construirse un sistema de ecuaciones lineales en el cual
todas las incognitas (los uij ) aparecen acopladas. La matriz de coeficientes es rala,
pero no presenta una estructura banda. La resolucion de sistemas de este tipo es
com
un encararla por medio de metodos iterativos (Jacobi, Gauss-Seidel, SOR).

3.2.

Residuos Ponderados

En la Seccion 2.3 se explico que una forma de encarar la obtencion de una solucion
aproximada a un problema diferencial, consiste en distribuir el error (o residuo) sobre
el dominio de caculo de acuerdo a ciertas funciones de peso. Por ejemplo, para el
problema definido por las Ecs. 2.24-2.26, se plantea la Ec. 2.27. Esta no es mas que
la formulacion del metodo de los residuos ponderados. En general, se impone que
la solucion aproximada satisfaga exactamente las condiciones de borde, Ecs. 2.252.26. La funcion de peso w se toma como una combinacion lineal de una familia de
funciones i , linealmente independientes, es decir,
w=

N
X

i i

(3.24)

i=1

donde los i son coeficientes arbitrarios. En consecuencia, la Ec. 2.27 conduce a

Z
(2 u + f )i d = 0
;
i = 1, 2, ..., N
(3.25)

Por su parte, la solucion u se expresa en terminos de funciones aproximantes i ,

es decir
N
X
u=
i i
(3.26)
i=1

Los coeficientes i , en la Ec. 3.26, se determinan de acuerdo a las Ecs. 3.25, lo cual
permite obtener la solucion aproximada.
Existe una gran variedad de metodos de residuos ponderados, determinados por
la particular eleccion de las funciones de peso i . En el metodo de los momentos, las
funciones de peso son tales que las Ecs. 3.25 son momentos del residuo. Por ejemplo,
en el caso unidimensional se tiene que
i = xi1
26

(3.27)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Supongase que se quiere resolver el problema definido en el Ec. 3.17 con este
metodo. Como funciones aproximantes se puede elegir
i = xi (1 x)

(3.28)

que satisfacen identicamente las condiciones de borde. Si se toma N = 2, la forma

de la solucion es, de acuerdo a la Ec. 3.26,
u = 1 x(1 x) + 2 x2 (1 x)

(3.29)

Entonces, el residuo  se expresa como

=

d2 u
+ u + x = (2 + x x2 )1 + (2 6x + x2 x3 )2 + x
dx2

Las ecuaciones equivalentes a las Ecs. 3.25 son, de acuerdo a la Ec. 3.27,
Z 1

1dx

0
Z

(3.30)

(3.31)

xdx

Reemplazando la Ec, 3.30 en las Ecs. 3.31 y operando se llega al siguiente sistema
 



1/2
11/6 11/12
1
=
(3.32)
2
1/3
1/2 19/20
de donde surge que

122

1 = 649

110

2 =
649
Es decir que la solucion, Ec. 3.29, es


122 110
u = x(1 x)
+
x
649 649

(3.33)

(3.34)

En la Tabla 3.2 se presentan valores calculados con la Ec. 3.34.

Otra variante del metodo de los residuos ponderados, denominada metodo de
Galerkin, consiste en tomar las funciones de peso identicas a las funciones aproximantes, es decir,
i = i
(3.35)
Aplicando este metodo al ejemplo antes considerado, se obtiene


71
7
u = x(1 x)
+ x
369 41
que da los valores mostrados en la Tabla 3.2.
27

(3.36)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Tambien existen los metodos de colocacion, de subdominio, etc. [12].

La diferencia conceptual entre el metodo de los residuos ponderados y el de
diferencias finitas, radica en la forma en que el error se distribuye sobre el dominio
de la solucion. Mientras que el segundo minimiza el error en un n
umero finito de
nodos. Es decir, es un metodo nodal, el primero lo distribuye sobre todo el dominio
de un modo determinado (es decir, es un metodo modal). Sin embargo, el metodo
de las diferencias finitas puede ser formalmente considerado como un caso particular
del metodo de los residuos ponderados, cuando las funciones de peso son deltas de
Dirac, es decir
i = (x xi )

(3.37)

en cuyo caso los coeficientes i devienen los valores nodales ui .

3.3.

Elementos finitos

El metodo de los elementos finitos parte de la formulacion debil o, en caso de

existir, la formulacion variacional del problema. Basicamente, el metodo consiste en
aproximar la funcion (o funciones) incognita por medio de una funcion continua a
trozos. Esta funcion aproximante es continua sobre cada uno de los elementos en
que se ha dividido al dominio de calculo, y all se expresa en terminos de los valores
nodales de la funcion incognita. Estos trozos de funciones continuas se denominan
funciones interpelantes o de forma. La forma y distribucion de los elementos depende,
en cierta medida, del problema particular ([16]; [13]; [8]).
La forma de elemento mas simple y, a la vez, la mas utilizada cuando el dominio
es bidimensional, es la triangular. Si los vertices de un elemento triangular generico
se identifican con los ndices i, j, k la funcion incognita en dicho elemento se expresa
como
(e)

(e)

(e)

u(e) (x, y) = Ni (x, y)ui + Nj (x, y)uj + Nk (x, y)uk

(3.38)

donde las funciones de forma Nl , con l = i, j o k, son bilineales, es decir, tienen la

forma
(e)
(e)
(e)
al + bl x + cl y
(e)
(3.39)
Nl (x, y) =
2(e)
Los coeficientes de la Ec. 3.39 valen, en terminos de las coordenadas de los
vertices,

ai = xj yk xk yj

bi = yj yk

ci = x k y j
(3.40)

 1 xi yi 

 = Area del elemento

1
x
y

j
j

 1 xk y k 

28

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

donde se ha suprimido, por simplicidad, el suprandice (e), y se ha supuesto que los

nodos i,j,k estan ordenados en sentido antihorario. Siendo as, las expresiones para
los coeficientes correspondientes a los otros nodos surgen de las Ecs. 3.40 mediante
permutaciones cclicas.
Si se supone que las funciones u(e) dadas por la Ec. 3.38, se anulan identicamente
fuera del elemento (e), la funcion incognita se puede expresar como
u(x, y) =

M
X

u(e) (x, y)

(3.41)

e=1

donde M es el n
umero total de elementos.
Ahora bien, si el problema admite una formulacion variacional, la tecnica consiste
en reemplazar en el funcional I(u) la expresion dada por la Ec. 3.41, cuyas incognitas son, en u
ltima instancia, los valores nodales un . En consecuencia, la condicion
expresada en la Ec. (2.42) conduce a
N
X
I
un = 0
I =
u
n
n=1

(3.42)

donde N es el n
umero total de nodos. Como las variaciones un son arbitrarias, la
Ec. 3.42 conduce a las siguientes ecuaciones:
I
=0
n

n = 1, 2, ..., N

(3.43)

La resolucion del sistema 3.43 permite hallar los valores nodales un , y, por medio de
la Ec.3.41, calcular la funcion en cualquier punto.
Considerese, a titulo de ejemplo, nuevamente el problema unidimensional definido
en la Ec. 3.17. El funcional asociado es el de la Ec. (2.43) si se elimina el u
ltimo
termino, es decir,
!
 2
Z 1
1 du
u2

I(u) =
xu dx
(3.44)
2 dx
2
0
Las ecuaciones equivalentes a las Ecs. 3.38 - 3.40 son

(e)
(e)

u(e) (x) = Ni (x)ui + Nj (x)uj

xj x
(e)
Ni (x) =
xj xi

xi x
(e)

Nj (x) =
xi xj

(3.45)

de donde surge que

(e)

(e)
dNj
dNi
uj ui
du(e)
=
ui +
uj =
dx
dx
dx
x j xi

(3.46)

Notese que, siendo, el problema unidimensional, los elementos son segmentos.

29

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

La parte del funcional correspondiente al elemento (e) es, de acuerdo a las Ecs.
3.44 - 3.46
2
Z xj 




1
1 uj ui
(e)
Nie ui + Nje uj x Nie ui + Nje uj dx (3.47)

I (ui , uj ) =
xj x i
2
xi 2
Las derivadas de I(e) respecto de las incognitas ui y uj , pueden calcularse directamente a partir de la Ec. 3.47:

  Z xj

 Z xj

uj ui
I(e)
2
=

Ni dx ui
Ni Nj dx uj
ui
xj xi
xi
xi
Z xj


Nj xdx
(3.48)
xi

I(e)
=
uj

uj ui
xj xi

Z

xj

Nj2 dx

Z
uj

xi

xj

Ni Nj dx ui
 Z xj


Ni xdx

xi

(3.49)

xi

Teniendo en cuenta que, de acuerdo a las Ecs. 3.45

Z xj
Z xj
xj xi

2
2

N
dx
=
N
dx
=

j
i

xi
xi

Z xj

xj xi

Ni Nj dx =

6
xi
Z xj

Ni xdx = (x2j + xi xj + 2x2i )

xi

Z xj

1 2

Nj xdx =
(xi + xi xj + 2x2j )

6
xi
las Ecs. 3.48 y 3.49 se transforman en




I(e)
1
1
xj xi
xj xi
ui
uj
=

+
ui
xj xi
3
xj xi
6
1
(x2j + xi xj + 2x2i )
6
I(e)
=
uj

1
xj xi

xj xi
3


uj


1
xj xi
+
ui
xj xi
6
1
(x2i + xi xj + 2x2j )
6

(3.50)

(3.51)

(3.52)

Obviamente, el funcional total es la superposicion de las partes correspondientes

a cada elemento, es decir
N
X
I=
Ie
(3.53)
e=1

30

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

Figura 3.2: Elementos finitos en un problema unidimensional

Para plantear las Ecs. 3.43 es necesario, entonces, proceder a ensamblar las contribuciones correspondientes a cada elemento que converge al nodo n. En un problema unidimensional, como el presente, a cada nodo convergen solo dos elementos;
para uno de ellos el nodo n juega el rol del nodo i del elemento generico considerado
mas arriba, mientras que para el restante el nodo n es como el nodo j (ver Fig. 3.2).
En consecuencia, de acuerdo a las Ecs. 3.51 y 3.52, se tiene que




1
xn+1 xn
xn+1 xn
1
I
=

+
un
un+1
un
xn+1 xn
3
xn+1 xn
6


1 2
xn xn1
1
2
(xn+1 + xn+1 xn + 2xn ) +

un
6
xn xn1
3


1
xn xn1

+
un1 +
xn xn1
6
1
(3.54)
+ (x2n1 + xn1 xn + 2x2n ) = 0
6
que puede reescribirse como



1
xn xn1

+
un1 +
xn xn1
6


1
1
xn+1 xn1
+
+

un
xn+1 xn xn xn1
3


1
xn+1 xn

+
un+1
xn+1 xn
6


xn+1 xn

(xn+1 + xn1 + xn ) = 0
6

(3.55)

Si, como cuando se resolvio este problema por diferencias finitas, se toma un
paso h = 1/3 uniforme, la Ec. 3.55 se reduce






h2
2h2
h2
1+
un+1 2
un + 1 +
un1 = nh3
(3.56)
2
3
6
La Ec. 3.56 es muy similar a la Ec. 3.20 obtenida por diferencias finitas. En
realidad, puede obtenerse directamente de un planteo como el de la Ec. 3.19, si se
reemplaza el segundo termino, uj por el promedio pesado (uj1 + 4uj + uj+1 )/6.
31

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

La formula generica 3.56, juntamente con las condiciones de contorno geometricas

en los nodos de borde, conduce al siguiente sistema


 

52/27 55/54
u1
1/27
=
(3.57)
55/54 52/27
u2
2/27
de donde surgen los valores mostrados en la Tabla 3.2.
Supongase ahora que se parte de una formulacion debil, como la Ec. 2.30, en
lugar de una variacional. Es necesario, entonces, definir la funcion de peso w. Si,
como en el metodo de los residuos ponderados, esta se expresa de acuerdo a la Ec.
3.24, lo que deben definirse son las funciones. La eleccion mas com
un, que da lugar
al metodo de Bubnov-Galerkin (o Galerkin a secas), consiste en tomar las funciones
de peso igual a las de forma, es decir
(e)

(e)

i = Ni

(3.58)

Si esta condicion no se cumple, se trata de un metodo de Petrov-Galerkin.

Si, una vez mas, se toma como ejemplo el problema definido por la Ec. 3.17, su
formulacion debil es

Z 1
du dw
+ (u + x)w dx = 0
(3.59)

dx dx
0
En un elemento generico, la funcion incognita se aproxima de acuerdo a las
(e)
(e)
Ecs. 3.45. Tomando como funcion de peso i = Ni , la parte de la Ec. 3.45 que
corresponde a ese elemento da lugar a


uj ui
xj xi

Z

xj

+
xi

Ni2 dx

Z

xj

Z

xj

Nj xdx

Ni Nj dx uj +

+
xi


(3.60)

xi

que es identica a la contribucion 3.48 de la formulacion variacional, Es decir, el

metodo de Bubnov-Galerkin, que se aplica a la formulacion debil de un problema, da
el mismo resultado que el metodo de los elementos finitos aplicado a la formulacion
variacional (cuando esta existe).

3.4.

Elementos de contorno

El punto de partida del metodo de los elementos de contorno es la formulacion

integral del problema. La esencia del metodo consiste en dividir al contorno del
dominio en segmentos o elementos, sobre 1os cuales se formulan aproximaciones
para evaluar las integrales que involucran a la funcion incognita [4].
Supongase que se parte de la formulacion directa dada en la Ec. 2.57. Si se toman
elementos constantes (tambien denominado metodo panel), es decir, se considera
que u y q u/n son constantes sobre cada segmento j , dicha ecuacion puede
aproximarse por
!
!
Z
Z
Z
N
N

X
X
1
u

ui +
u d qj +
d uj + gu d = 0
2
n
j
j

j=1
j=1
i = 1, 2, ..., N
32

(3.61)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

donde N es el n
umero total de elementos de contorno. Las integrales que aparecen en
las Ecs. 3.61 pueden calcularse, ya sea en forma exacta o aproximada. Se denominara
Z
Gij =
u d
(3.62)
j

ij =
H

u
d
n

Z
j

gu d

Bi =

(3.63)
(3.64)

Si, ademas, se define


Hij =

ij
H
ij +
H

si i 6= j
si i = j

1
2

(3.65)

las Ecs. 3.61 pueden reescribirse como N

N
X

Gij qj

N
X

j=1

Hij uj + Bi = 0

(3.66)

j=1

Si de los N elementos de la frontera , N1 estan asociados a la parte 1 y

N2 z(= N N1 ) a la parte 2 , las Ecs. 3.66 forman un sistema de N ecuaciones
algebraicas lineales con N1 incognitas qj y N2 incognitas uj , es decir, esta bien
determinado.
Una vez conocidos u y q sobre los elementos de contorno, los valores de la solucion
en puntos interiores a se obtienen dscretizando la Ec. 2.58, es decir
ui =

N
X

Gij qj

j=1

N
X

ij uj + Bi
H

(3.67)

j=1

En general, las integrales Gi j y Hij ,dadas por las Ecs. 3.62 y 3.63, respectivamente, se eval
uan numericamente utilizando la cuadratura de Gauss mas simple,
salvo las correspondientes a i = j , que requieren mayor precision. En el caso de elementos constantes estas ultimas pueden ser calculadas anallticamente; por ejemplo,
en problemas bidimensionales se tiene que

Hii = 0
(3.68)
1
ri

Gii = [ln + 1]
ri
donde ri es la semi-longitud del elemento i.
Por su parte, para calcular la integral Bi , dada por la Ec. 3.61, es com
un dividir el
dominio en peque
nos subdominios (por ejemplo, de forma triangular), y efectuar
aproximaciones locales (en general, cuadratura de Gauss), es decir
Z
X
Bi =
M(
gu d)
(3.69)
k

k=1

33

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

donde M es el n
umero total de subdominios o elementos internos.
Si se toma, nuevamente, como ejemplo de aplicacion al problema unidimensional definido por la Ec. 3.17, los elementos de contorno no son mas que los puntos
extremos del intervalo, x = 0 y x = 1, que se denominaran a y b, respectivamente.
La solucion fundamental esta dada en la Ec. 2.59, con lo cual se obtiene inmediatamente, de acuerdo a la Ec.3.62,

Gaa = Gbb = 0
(3.70)

Gab = Gba = 1 sen(1)

2
Las integrales dadas por las Ecs. 3.63 no necesitan ser evaluadas, ya que u se anula
en los bordes. Por su parte, las integrales correspondientes a las Ecs. 3.64 son

Z
1 1

xsen(x)dx
Ba =
2 Z0
(3.71)
1 1

Bb =
xsen(1 x)dx
2 0
Siguiendo la idea del metodo, las integrales de las Ecs. 3.71 se calcularan dividiendo
el dominio en tres intervalos iguales, y utilizando cuadratura de Gauss para cada
uno de ellos (a pesar de que esas integrales pueden evaluarse exactamente). Por
ejemplo, Ba puede escribirse como
Z 1
(f1 () + f2 () + f3 ())d
(3.72)
Ba =
1

donde

f1 () =

f2 () =

f () =
3

1 ( + 1)
g
6
6
1 ( + 3)
g
6
6

(3.73)

1 ( + 5)
g
6
6

siendo

1
g(x) = xsenx
(3.74)
2
Entonces, aplicando la formula de cuadratura de Gauss con dos puntos, la Ec.
3.72 da
Ba = (f1 + f2 + f3 )|=0,57735 + (f1 + f2 + f3 )|=0,57735 = 0,15059

(3.75)

De manera analoga, pero con

1
g(x) = xsen(1 x)
2
34

(3.76)

Introduccion a la Simulacion Numerica de Problemas Hdraulicos

en lugar de la Ec. 3.74, se obtiene que

Bb = 0,07926
Las Ecs. 3.66 dan, en este caso,

Gaa qa + Gab qb + Ba = 0
Gba qa + Gbb qb + Bb = 0

(3.77)

(3.78)

Reemplazando los valores dados en las Ecs. 3.70, 3.75 y 3.77 en las Ecs. 3.78, se
obtiene

qa = 0, 18838
(3.79)
qb = 0, 35792
(que tienen cuatro dgitos correctos, comparados con los valores exactos).
Ahora pueden calcularse los valores de la solucion en el interior, de acuerde con
la Ec. 3.67. Para comparar con los metodos discutidos anteriormente, se evaluara la
solucion en x = 1/3 y x = 2/3, que se identificaran como los puntos 1 y 2, respectivamente. Entonces, se tiene que

u1 = G1a qa + G1b qb + B1
(3.80)
u2 = G2a qa + G2b qb + B2
Pero

1
1

G1a = G2b = 2 sen 3

(3.81)

2
1

G1b = G2a = sen

2
3
Ademas, con g(x) = (x/2)sen(x 1/3) y con g(x) = (x/2)sen(x 2/3) se
obtienen, respectivamente,

B1 = 0,08598
(3.82)
B2 = 0,04860
Reemplazando los va/ores dados en las Ecs. 3.79, 3.81 y 3.82 en las Ecs. 3.80, se
obtienen los resultados mostrados en la Tabla 3.2, cuyos digitos significativos son
todos correctos.

35

Bibliografa
[1] Batchelor, G. K., An Introduction to Fluid Dynamics, Cambridge Univ. Press.,
(1967)
[2] Bender, C. M., Orszag, S. A., Advanced Mathematical Methods for Scientists
and Engineers, McGraw-Hill, (1978)
[3] Bradshaw, F. (Ed.), Turbulence, Springer-Verlag, (1976)
[4] Brebbia, C. A., Walker, S., Boundary Element Techniques in Engineering,
Newnes-Butterworths, (1980)
[5] Elsgoltz, L., Ecuaciones diferenciales y calculo variacional, Mir, (1969)
[6] Forsythe, G. E. Wasow, W., Finite Difference Methods for Partial Differential
Equations, Wiley, (1960)
[7] Hinze, J., Turbulence, McGraw-Hill, (1975)
[8] Huebner, K. H., Thornton, E. A., The Finite Element Method for Engineers,
Wiley, (1982)
[9] Marshall, G., Solucion numerica de ecuaciones diferenciales, Reverte (Tomo 1),
(1985)
[10] Marshall, G., Solucion numerica de ecuaciones diferenciales, Reverte (Tomo 2),
(1986)
[11] Morse, P. M., Feshback, H., Methods of Theoretical Physics, McGraw-Hill,
(1952)
[12] Nakamura, S., Computational Methods in Engneering and Science, Wiley,
(1977)
[13] Norrie, D. H., de Vries, G., An Introduction to Finite Element Analysis, Academic Press, (1978)
[14] Richtmyer, R. D., Morton, K. W., Difference Methods for Initial-Value Problems, Wiley, (1967)
[15] Richtmyer, R. D., Principles of Advanced Mathematical Physics, SpringerVerlag (Vol I), (1978)
[16] Zienkiewicz, U. C., The Finite Element Method, McGraw-Hill, (1977)
36