PRÁCTICA 1, CIV-334

SOLUCIÓN NUMÉRICA
DE PERFIL DE SUPERFICIE DE FLUJO EN CANAL ABIERTO CON
ENTRADA LATERAL Y APLICACIÓN DE NEWTON RAPHSON

1. OBJETIVOS
Establecer la solución numérica a la ecuación diferencial ordinaria de flujo permanente
y representar el perfil de flujo para el caso analizado.
2. MARCO TEÓRICO
La solución numérica ha recobrado mucha atención para aplicar y hallar soluciones a
los problemas hidrodinámicos. Es un adelanto y una alternativa a los estudios y
soluciones analíticas, o conocidas como soluciones meramente aritméticas y donde los
resultados son singulares o puntual; en cambio la solución numérica puede abarcar a
una escala de espacio.
Las bases para la solución de éste tipo de casos, mejor dicho la solución de ecuaciones
diferenciales, se fundamenta en la teoría de los métodos numéricos, como los más
simples, los otros sobre los conceptos de diferencias finitas y finalmente sobre los
elementos finitos. En esta oportunidad se aplicará el método numérico más popular:
Newton Raphson.
Para ello, iniciaremos con la descripción de la ecuación diferencial fundamental que
describe el movimiento de flujo y el perfil de flujo, conocida por la ley de Saint Venant
principal para flujo.
dQ
q (1)
dx

dE d   Q2  Q
 h 2 
 S  2 q (2)
dx dx  2 gA  gA

Donde,

h - nivel de agua sobre el nivel base.

Q - descarga.

A - área transversal de flujo.

x - longitud, distancia horizontal.
q - entrada lateral de flujo.

 - factor de corrección energético.

S - factor de fricción energético:

Estas ecuaciones son aplicables para el flujo en todo tipo de canales. La ecuación
diferencial ordinaria clásica que describe el perfil de la superficie del agua para un flujo
gradualmente variado derivada, para un tramo recto, la ecuación de Saint Venant tiene
la forma:

dH sS
 (3)
dx 1  Fr 2

Donde
H - profundidad de flujo.
s - pendiente del canal.
Fr - número de Froude.

El factor de fricción energético S esta descrito por la siguiente ecuación

n 
2
Q2
S M (4)
R 4/3 A2

Donde

nM - coeficiente de Manning.

R - radio hidráulico.

El número de Froude se calcula mediante la siguiente ecuación.

U
Fr  (5)
gH

Este número de Froude establece el valor crítico denotado por

H C cuando adquiere el
valor de Fr  1 .

La ecuación (4.60) es una ecuación diferencial de H con respecto a x , que necesita

valores iniciales para ser resuelta. Estos valores de H inicial pueden ser cualquiera,
H C . Puesto que en estos sistemas
siempre y cuando estos sean mayores al valor crítico
también se tienen parámetros como la altura normal, se puede tomar a este como otro
valor inicial.

Altura normal.

La altura normal está dada por la siguiente ecuación.

1 2/3 1/2
Q R s  A (6)
nM

Tomando en cuenta las ecuaciones anteriores se obtiene la siguiente expresión.

Q
1

 bH n  mH n 
2 3 1
 s2  0 (7)

2

 b  2H
nM
n 1 m 2 3

H
Esta es una ecuación algebraica no lineal. Se resuelve esta ecuación para n de forma
numérica. Entre los métodos se podría usar bisección, Newton-Raphson, punto fijo, o
cualquier otro método numérico existente.

Altura critica.

La altura crítica puede ser calculada por la siguiente ecuación.

  Q 2 A3
 (8)
g B

Tomando en cuenta las ecuaciones anteriores, se obtiene la siguiente ecuación.

g
2 
bH c  mH c 2   0
3
 b  2mH c   (9)
Q

Otra vez se tiene una ecuación algebraica no lineal. Esta, al igual que para la altura
Hc .
normal se puede utilizar cualquier método numérico para obtener el valor de

3. PROCEDIMIENTO DE ESTUDIO DEL DERFIL DE FLUJO

3.1 Perfil de superficie de flujo en canal abierto con aporte lateral

La ecuación que describe el movimiento, también derivada de Saint Vennat para el caso
de presente flujo está dado:

U U h
U   g  gSo  gS f (10)
t x x
Donde
U - velocidad promedio.

h - profundidad de agua o tirante del flujo.

SO
- pendiente de la base del canal.

Sf
- pendiente de la línea de energía.
g - aceleración debido a la gravedad.

x - coordenada a lo largo del canal.

t - tiempo.

La pendiente de la línea energética puede ser determinada de la ecuación de Manning:

n 2U 2
Sf  4 (11)
h3

Para condiciones de flujo permanente si tiene la siguiente situación:

U
0 (12)
t

La ecuación (10) puede reescribirse como:

U h
U   g  gSo  gS f (13)
x x

Esta ecuación al depender de una sola variable, puede escribirse como una ecuación
diferencial total. Además, el primer término de la derecha se la puede cambiar por su
forma alterna.

d U2  dh
   g  gSo  gS f (14)
dx  2  dx

Reordenando,

d  U 2  dh
   gSo  gS f (15)
dx  2 g  dx

En esta última ecuación, tanto U como h están en función de solamente x .

Ecuación de continuidad.
Ahora es necesario hacer un análisis de la ecuación de continuidad. Que para un canal
con aporte lateral de flujo está dada por la siguiente ecuación:

h 
  Uh   p (16)
t x

Dónde:

p - entrada lateral de flujo.

Esta última ecuación para un flujo permanente si la variación de h en relación de la

variación de t es cero; entonces la ecuación (16) se reduce:


 Uh   p (17)
x

Y de la misma forma, ya que depende de una sola variable, esta pasa a ser una ecuación
diferencial total.

d
 Uh   p (18)
dx

Análisis numérico.

Las ecuaciones (15) y (18) deben ser resueltas simultáneamente para obtener los dos
valores desconocidos U y h . En general esto debe ser realizado numéricamente ya que
la solución analítica no es posible.

Decimos sistema de ecuaciones de diferencias totales:

d  U 2  dh
   So  S f (19)
dx  2 g  dx

d
 Uh   p (20)
dx

En una aproximación numérica las derivadas son expresadas como diferencias.

Suponiendo que el canal está dividido en N celdas, cada una con la misma longitud x .

Esquema: Flujo con aporte lateral, P

La ecuación (20) puede ser escrita de la siguiente forma, o conocidas por su forma
discretizada:

U i 1hi 1  U i hi
p (21)
x

Donde el índice i se refiere al límite aguas arriba de la celda, y el subíndice i  1 se

refiere al límite aguas debajo de la celda (ver en Figura 1).

Figura 1: Celda en la posición (i) e (i+1)

La ecuación (21) puede ser reescrita como:

 U i 1hi 1  px Ap
Ui   (22)
hi hi

Se supone que se tiene un flujo subcrítico; por tanto, los cálculos procederán aguas
arriba. Suponiendo que los datos en i  1 son conocidos, además de que el flujo de
aporte lateral canal es
p y el paso x también son conocidos, es posible introducir el
coeficiente A .

U i 1hi 1  px  Ap (23)

La ecuación (19) puede ser discretizada como,

1  U i21  U i2  hi 1  hi 1
   So   S f ,i  S f ,i 1  (24)
2 g  x  x 2

Donde la pendiente de la línea de energía se toma como el promedio entre las

condiciones en i e i  1 . Multiplicando por x y agrupando todos los términos con
subíndice i al lado izquierdo se tiene:

U i2 1 U2 1
hi   xS f ,i  hi 1  i 1  xSo  xS f ,i 1 (25)
2g 2 2g 2

En esta última ecuación se puede observar que todos los términos de la derecha son
conocidos; por tanto, también si se puede introducir el coeficiente B de la siguiente
forma:

U i21 1
hi 1   xSo  xS f ,i 1  B (26)
2g 2

Ya que es una igualdad este coeficiente también puede utilizarse con el lado izquierdo.

U i2 1
hi   xS f ,i  B (27)
2g 2

La pendiente de la línea de energía S puede expresarse mediante la ecuación de

Manning. Reemplazando en la ecuación (27) se tiene,

U i2 1 n 2U i2
hi   x 4/3  B (28)
2g 2 hi

De la ecuación (22),
 Ap
Ui  (29)
hi

Utilizando la ecuación (29) en la ecuación (28) puede reescribirse como:

2 2
1  Ap  1 n 2  Ap 
hi     x   B (30)
2 g  hi  2 hi4/3  hi 

Reescribiendo la ecuación (30) se tiene,

Ap 2 1 C
hi  2
 10/3  B (31)
2 g hi hi

C es otro coeficiente introducido. Su expresión agrupa los términos x , n , y A .

1
C  xn 2 Ap 2 (32)
2

3.2 Solución del modelo (31) mediante el método numérico

El modelo base cálculo y mediante un método numérico está dado para la solución de
f(h):

Ap 2 1 C
f (h)  hi  2
 10/3  B
2 g hi hi
(33)

De este modelo se presentan las ecuaciones que comprenden el modelo numérico.

A
Ui  (44)
hi

A  U i 1hi 1  px (45)

U i21 1
B  hi 1   xSo  xS f ,i 1 (46)
2g 2

1
C  xn 2 A2 (47)
2

Ap 2 1 C
hi  2
 10/3  B (48)
2 g hi hi
3.3 Datos de entrada:

P  0.01[m3 / m 2 / s]

n  0.02
so  0.003

L  100[m]

g  9.81 m
s2

El paso para realizar los cálculos es de

x  10[m]

3.4 Algoritmo
Se deben conocer los valores para el vertical o punto más bajo ((N+1)

Para el punto más bajo el carácter de flujo (Fr) es 1:

U N 1
Fr  1
g * hN 1
(a)

De la continuidad, el flujo en ese punto es,

qN 1  pL
(b)

Dónde: L - es la longitud del canal.

Ya que

qN 1  U N 1hN 1
(c)

La velocidad de flujo en punto más bajo:

U N 1   PLg  , m / s
1/3

(d)

Derivando el tirante en la ecuación de Froude (Fr), tenemos.

  U N 1 
2

hN 1  ,  m
g
(e)

Ahora se procede a calcular los coeficientes Ap, C y B:

Ap  U i 1 * hi 1  Px
(f)

h11  hN 1

1
C x * n 2 Ap2
2 (g)

Considerando las ecuaciones 28 y 29, B será:

1 Ap 2 1 n 2 Ap
B  hi  ( )  x 4/3 ( ) 2
2 g hi 2 hi hi (h)

Donde Ui es:

Ap
ui  ( )
hi (i)aa01

Finalmente computar la ecuación principal de perfil de flujo aplicando la ecuación de Newton-

Raphson:
Ap 2 1 C
f (h)  hi  2
 10/3  B
2 g hi hi
(j)
Se recomienda para la iteración iniciar hi=hi+1, m que es resultado de la ecuación (h)
Ecuación de Newton-Raphson:
f ( h)
hi 1  hi 
f ' ( h) (k)

para iniciar hi=h i+1 despues 0.40, 0.5

 f (h  )  f (h )
f'
 (l)

  0, 001

Aproxactual  Aproxanterior
a  100%
Aproxactual (m)
De esa manera el siguiente paso es calcular el resto de los tirantes del flujo hasta que x=100, con
intervalos de 10m.

Tabla de control para el cómputo manual:

- Al calcular la ec. (e) para el cual x=0, usted ha obtenido las primeras coordenadas para
la gráfica
x=0 m, h1= …… m ?
1ra iteración
hn+1 Δx, m Ap C B hi f(h) f(h+Δ) f’ hi+1 Error, Ec
(m)
Ec. (e) es el 10
valor de h1

 0,00
Con la 1ra iteración usted ha hallado la segunda coordenada del perfil de flujo
x2=10, h2
2da iteración
hn+1 Δx, m Ap C B hi f(h) f(h+Δ) f’ hi+1 Error
Ec. 10
(h)

 0,00
Con la 2ra iteración usted ha hallado la segunda coordenada del perfil de flujo
x3=20, h3

.
.
.
Ultima iteración
hi+1 Δx, m Ap C B hi f(h) f(h+Δ) f’ hi+1 Error
Ec. 10
(e)

 0,00
Con la 1ra iteración usted ha hallado la segunda coordenada del perfil de flujo
x10=100, h10

3.5 Solución: representación gráfica del perfil de flujo

Siguiendo el procedimiento establecido, debe obtener un perfil de flujo de agua como el
siguiente, con los valores de x, h:

1.2

0.8
Profundidad, [m]

0.6

0.4

0.2

0
10 20 30 40 50 60 70 80 90 100
Distancia, [m]

Literatura
Steven Chapra, Métodos numéricos, MC GRAW HILL, 2010
(Repasar capítulo 6)