FACULTAD DE INGENIERÍA, CIENCIAS FÍSICAS Y MATEMÁTICA
ESCUELA DE INGENIERÍA CIVIL
PAQUETES COMPUTACIONALES EN HIDRÁULICA
DR. ING. JOAO LEAL
CARRERA ANDRÉS
CHACASAGUAY MÓNICA
CUASPA FREDDY
DURÁN ANITA
NOVENO
TERCERO
11 – 06 06 – 2013 2013
MARZO 2013 AGOSTO 2013 –
QUITO ECUADOR –
1. CÁLCULO HIDRÁULICO BÁSICO ENUNCIADO DEL PROBLEMA 1.1 Desarrolle una función (“function”) en MATLAB que calcule la curva de descarga (altura
del flujo vs. caudal) para las siguientes secciones transversales, admitiendo que el cauce es compuesto por arena gruesa ( d 50 50 = 4 mm): a) Sección rectangular con 25 m de ancho; b) Sección trapezoidal con 25 m de ancho en el fondo fondo y pendientes de las las orillas de 1:2 (vertical: horizontal). c) Sección de un río con los siguientes datos medidos: Tabla 1 Datos del cauce de un río
Distancia PUNTO Transversal (m)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
0.00 2.98 6.74 9.33 12.57 14.77 17.36 19.96 22.55 24.75 27.34 30.06 33.30 36.28 40.78 42.20
Elevación del cauce (m)
184.00 182.00 180.00 178.00 176.00 174.00 172.00 170.00 168.00 166.00 164.00 162.00 160.00 158.00 156.00 154.00
PUNTO
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
Distancia Transversal (m)
65.53 68.64 74.21 78.75 82.89 86.78 90.15 93.39 97.02 101.04 107.00 109.07 111.27 114.12 117.36 120.34
Elevación del cauce (m)
154.00 156.00 158.00 160.00 162.00 164.00 166.00 168.00 170.00 172.00 174.00 176.00 178.00 180.00 182.00 184.00
(NOTA: la función debe ser lo más general posible, el área debe ser llamada como una función y la pendiente deberá ser un parámetro de entrada)
2
SOLUCIÓN AL PROBLEMA 1.1 1.1-1 Función Descarga:
Primero creamos una función que en este caso se llama descarga y el código de programación es el que se muestra a continuación:
1 2 3 4 5 6
function Q=descarga(h,area,perimetro,pendiente) function d50=4; d50=d50/1000;%en d50=d50/1000;%en milimetros n=0.047*(d50^(1/6)); Q=(1/n).*(((area).^(5/3))./((perimetro).^(2/3))).*sqrt(pendiente); end
En la función (“función”) anterior se puede observar que los datos de entrada necesarios son: ( h,area,perimetro,pendiente h,area,perimetro,pendiente ), de los cuales el área y el perímetro serán llamados como una función mientras que h y pendiente serán parámetros de entrada. En este caso se utiliza la ecuación de Manning en la línea 5 y es la siguiente:
Dónde: Q= Caudal, m³/s n= Coeficiente de rugosidad de Manning A= Área de la sección transversal P= Perímetro mojado de la sección transversal
El dato del problema para la rugosidad del cauce fue una arena gruesa ( d 50 50 = 4 mm), y en todo el desarrollo del problema nosotros utilizaremos el coeficiente n de Manning mismo que se obtiene mediante la siguiente expresión: 1.1-1.1 Cálculo de coeficiente n de Manning:
Expresión de: Garde & Raju, 1978; Subramanya, 1982:
Donde en d50 está en metros
3
Como se puede observar en el código dentro de la función llamada “descarga” ya se incluye el d50=4mm y por ende también se calcula el n de Manning . Una vez solucionado el coeficiente n de Mannig procedemos a crear las funciones de área y perímetro para cada una de las secciones dadas: 1.1-2 Cálculo del área para la sección rectangular:
La función creada para resolver el área de la sección rectangular fue llamada como area_rec y el código de programación es el que se muestra a continuación: function S=area_rec(h) S=area_rec(h) 1 function 2 %PROGRAMA PARA CALCULAR EL ÁREA DE UN CANAL 3 %RECTANGULAR S=B*h 4 %DONDE: 5 %S=Área de la sección transversal del canal en m² 6 %B=ancho del canal en m 7 %h=altura del espejo de agua en m 8 B=25; 9 S=(B.*h)'; 10 End
Descripción del contenido de la función area_rec
En el código anterior se observa que el único parámetro de entrada es h ya que el ancho B=25m ya viene incluido dentro de esta función y se lo observa claramente en la línea 8 del código anterior, la sección y dimensiones adoptadas son las que se muestra en la siguiente figura. El área es calculada con la formula
Superficie libre del agua
h
B=25m
Fig. 1 Sección Rectangular
1.1-3 Cálculo del perímetro para la sección rectangular:
De la misma forma y criterio que para el cálculo cálculo del área se ha desarrollado la función función per_rec para el cálculo del perímetro de una sección rectangular y el código de programación es el siguiente:
4
1 function function P=per_rec(h) P=per_rec(h) 2 %PROGRAMA PARA CALCULAR EL perímetro DE UN CANAL 3 %RECTANGULAR P=B+2*h 4 %DONDE: 5 %P=Perímetro del canal en m 6 %B=ancho del canal en m 7 %h=altura del espejo de agua en m 8 B=25; 9 P=B+2.*h; 10 P=P'; 11 end La fórmula utilizada en la línea 9 para el cálculo del perímetro es:
1.1-4 Cálculo del área para la sección trapezoidal:
Superficie libre del agua
1
1 m=2
h
m=2
B=25m
Fig. 2 Sección Trapezoidal
La función utilizada para el cálculo del área de una sección trapezoidal fue llamada con el nombre de area_tra y el código de programación es el siguiente: function S=area_tra(h) S=area_tra(h) 1 function 2 %PROGRAMA PARA CALCULAR EL ÁREA DE UN CANAL 3 %RECTANGULAR S=h*(B+m*h) 4 %DONDE: 5 %S=Área de la sección transversal del canal en m² 6 %B=ancho del canal en m 7 %h=altura del espejo de agua en m 8 %m=inversa de la pendiente 9 B=25; 10 m=2; 11 S=(h.*(B+m.*h))'; 12 End
5
Descripción del contenido de la función area_tra
El área que consta en en la línea línea 11 se calculó con la siguiente expresión:
Dónde:
S= área de la sección transversal, m² h= altura del espejo de agua, m m= relación horizontal de pendiente o pendiente inversa. B= ancho de la base de la sección. En ancho B=25 fue ingresado dentro de la función y consta en la línea 9 del código de programación de la función area_tra, así como también m=2 consta en la línea 10 del mismo código. 1.1-5 Cálculo del perímetro para la sección trapezoidal:
La función utilizada para el cálculo del perímetro de una sección trapezoidal fue llamada con el nombre de per_tra y el código de p rogramación es el siguiente:
1 function function P=per_tra(h) P=per_tra(h) 2 %PROGRAMA PARA CALCULAR EL perímetro DE UN CANAL 3 %Trapezoidal P=B+2*h*(1+m^2)^(1/2) 4 %DONDE: 5 %P=Perímetro del canal en m 6 %B=ancho del canal en m 7 %h=altura del espejo de agua en m 8 %m= inversa de la pendiente de las paredes del canal. 9 B=25; 10 m=2; 11 P=B+2.*h.*sqrt(1+m^2); 12 P=P' 13 end En esta función per_tra ya se incluye el ancho B y la inversa de la pend iente m en las líneas 9 y 10 respectivamente. En la línea 11 se incluye el cálculo del perímetro utilizando la siguiente formula:
√ 6
1.1– 6 Cálculo del área para la sección del río:
Cota máx=184.00
Superficie libre del agua
Pto 32
Pto 1 hmáx=30 Cota mín=154.00
Pto 16
23,33
Pto 17
Fig. 3 Sección del río
La función utilizada para el cálculo del área de una sección d el río fue llamada con el nombre de area_rio y el código de programación es el siguiente:
1 function function S=area_rio(h) S=area_rio(h) 2 clc 3 CONTADOR=0; 4 G=h; seccion=xlsread('DATOS_CAUCE.xlsx', ,'hoja1' 'hoja1', ,'A2:B33' 'A2:B33'); ); 5 seccion=xlsread('DATOS_CAUCE.xlsx' 6 cotamax=max(seccion(:,2)); 7 cotamin=min(seccion(:,2)); 8 conta=length(seccion); for k=1:length(G) k=1:length(G) %para hacer valor por valor 9 for 10 vec=G(1,k); %cada valor del vector. 11 cota=vec+cotamin; 12 if if cota cota > cotamax 13 error('COTA error('COTA MAYOR QUE LA MÁXIMA PERMITIDA') PERMITIDA') 14 end if cota cota < cotamin 15 if error('COTA error('COTA MENOR QUE LA MÍNIMA PERMITIDA') PERMITIDA') 16 17 end 18 area=0; 19 ii=0; while 1 1 20 while i=conta/2+1+ii; %i va de 17 a 31 y j de 16 a 2 21 ii=ii+1; 22 j=conta+1-i; 23 B=seccion(i,1)-seccion(j,1); 24 25 hh=seccion(i+1,2)-seccion(i,2); 26 mizq=(seccion(j,1)-seccion(j-1,1))./hh;
7
27 mder=(seccion(i+1,1)-seccion(i,1))./hh; 28 h1=cota-seccion(i,2); if h1<=hh if h1<=hh 29 area2=h1.*(B+(mizq.*h1+mder.*h1)./2); 30 CONTADOR=CONTADOR+1; 31 32 S1=area+area2; 33 S(CONTADOR)=S1; 34 break 35 end 36 area=hh.*(B+(mizq.*hh+mder.*hh)./2)+area; area=hh.*(B+(mizq.*hh+mder.*hh)./2)+area;%area %area cada tramo completo 37 end 38 end 39 S=S'; 40 End Descripción del contenido de la función area_río
Los datos del cauce del río mostrados en la tabla1, fueron almacenados en una hoja de Excel llamada “DATOS_CAUCE” para luego ser llamados a una matriz dentro de la función area_rio por medio del comando de importación de datos de Excel tal como se muestra en la línea 5 del código asi: seccion=xlsread('DATOS_CAUCE.xlsx' seccion=xlsread('DATOS_CAUCE.xlsx', ,'hoja1' 'hoja1', ,'A2:B33' 'A2:B33'); );
NOMBRE DEL ARCHIVO DE EXCEL CON EXTENSIÓN .xlsx
NOMBRE DE LA HOJA DONDE SE ENCUENTRAN LOS DATOS
RANGO DE CELDAS A IMPORTAR
En las líneas 6 y 7 se calculan las cotas máxima y mínima del cauce, es decir los valores máximos y mínimos de la segunda columna de la matriz llamada “sección” ya que la primera columna de esta matriz corresponde a las distancias transversales mientras que la segunda columna corresponde a las cotas o elevaciones de cada uno de los puntos que conforman la sección del cauce. En la línea 8 se determina el número de puntos que tiene la sección del río con el comando length de MATLAB y este valor es almacenado en la variable llamada “conta”. La sentencia for que consta en la línea 9 va desde un valor de 1 hasta el número de valores que tenga nuestro vector h, esto nos va a servir para poder calcular cada una de las áreas que correspondan a cada uno de los elementos del vector de alturas llamado h que en este caso caso ya fue reemplazado por G como consta en la línea 4.
8
La variable de la línea 10 llamada “vec” se encargara de proporcionarn proporcionarn os cada uno de los valores de h, mismos que consta n en la columna 1 de la matriz “G”
En la línea 11 se calcula las cotas que corresponden a cada uno de los valores ingresados y esto se hace sumando la cota más baja del cauce a cada una de las alturas del vector h o G que en este caso es el mismo. Las sentencias if que se observan en las líneas 12 y 15 detendrán el cálculo en caso de que algún dato de altura se encuentre sobre la cota máxima o bajo la cota mínima y arrojara los correspondientes mensajes de error. La sentencia while que consta en la línea 20 se encargara de sumar las áreas desde la cota más baja hasta la cota que el usuario requiera y se detendrá cuando ya tenga calculada el área a la altura dada. Los contadores que se muestran en las líneas 21 y 23 van desde los puntos más bajos del cauce hasta los puntos más altos del mismo, en este caso el contador i va desde 17 hasta 31 y el contador j va desde 16 hasta 2, el motivo de no llegar hasta los puntos 32 y 1 se debe a que estos puntos son el espejo de agua más alto y solo necesitamos llegar hasta la última base de las secciones trapezoidales parciales que se presentan a lo largo de toda la altura del cauce. En la línea 24 se calcula cada uno de los anchos B restando cada una de las distancias transversales de los puntos i menos cada una de las distancias transversales de los puntos j en forma ascendente desde la cota más baja hasta la altura deseada. De la misma forma se calculan cada una de las alturas parciales entre punto y punto con la diferencia que en este caso se restan cada una de las abscisas de los puntos i+1 menos las abscisas de los puntos i de esta forma se resta el punto más alto menos el punto más bajo para evitar cualquier problema con los signos. En las líneas 26 y 27 se calculan cada una de las pendientes inversas de las paredes laterales del rio y se hace una variable para el lado izquierdo y otra variable para el lado derecho debido a que las dos pendientes laterales son diferentes si se tiene algún caso en que las dos pendientes son iguales no habrá ningún problema puesto que en este caso la variable mizq y mder serán las mismas. Para entender mejor las variables utilizadas y el proceso de cálculo se muestra a continuación un gráfico con los contadores i, j, pendientes,alturas, anchos que se utilizaron para cada uno de los trapecios parciales que forman la sección completa del río además de esto también se incluye un contador ii que se encarga de contar para posteriormente sumar cada una de las áreas calculadas:
9
hh(ii)
(seccion(i+1,1),seccion(i+1,2)) (seccion(i+1,1),seccion(i+1,2))
(seccion(j+1,1),seccion(j+1,2))
Superficie libre del agua
1
1 hh(ii) mizq
mder
coordenada x del punto i
(seccion(j,1),seccion(j,2)) (seccion(i,1),seccion(i,2)) coordenada x del punto j
B
coordenada x del punto j coordenada x del punto i Fig. 4 Sección parcial del río
El h1 que consta en la línea 28 representa cada una de las cotas dato del rio mismo que nos ayudara a detener el while cuando esta altura o cota sea menos o igual a la cota hh que es la suma de cada altura de trapecio parcial con contador i de 17 a 31 El if que consta en la línea 29 nos ayudara a detener todo el proceso con la ayuda del break que fue insertado en la línea 34, cuando se cumpla con la condición de h1<=hh, es decir, que detendrá todo el proceso cuando la altura que ingresamos se encuentra en uno de los trapecios parciales, una vez realizado esto se procede a calcular el área para una altura h1 que es menor o igual a la altura del trapecio parcial en análisis. Las dos fórmulas utilizadas para el cálculo del área parcial del trapecio como para el área cuya altura sea menor o igual a la del trapecio parcial son: Para un área menor o igual a la de un trapecio parcial es:
( ( )
Para un área completa de trapecio parcial sumado las áreas parciales de cada trapecio bajo el trapecio en análisis:
( ( )
Si la altura ingresada es menor o igual a la altura del primer trapecio la suma se realizara correctamente ya que la primera condición de “area”=0 ya fue ingresada
10
en la línea 18, por lo tanto no se suma nada en el caso del primer trapecio, pero si la altura dada es mayor que la la del primer trapecio de abajo hacia hacia arriba entonces como se observa el área se va almacenando den la memoria con variable denominada “area”
La fórmula de la variable “a rea” fue ingresada estratégicamente en la línea 36 ya que está ubicada luego del if y por ende no será contabilizada si la altura ingresada es menor que la del primer trapecio, si no se hacía esto la función tendría un error ya que en el primer trapecio parcial se sumaría un trapecio inicial que no existe. En la línea 32 se suma todos los trapecios parciales que se encuentren bajo el trapecio en análisis y el área cuya altura es menor o igual que la altura del trapecio en análisis y se esta forma se obtiene el área buscada con la suma de las áreas parciales. En la línea 33 se almacena todas las áreas S1 en un vector llamado “S” , y esta es
la salida final de las áreas correspondientes a cada una de las alturas ingresadas en la función area_rio. Finalmente en la línea 39 se obtiene la transpuesta del vector de áreas (S) para visualizar de una mejor manera los resultados. 1.1– 7 Cálculo del perímetro para la sección del río:
El perímetro de la sección de río fue calculado con un criterio similar al utilizado en el cálculo del área, la función creada para este fin tiene el nombre de p er_rio y el código de programación es el que se muestra a continuación:
1 function function [P]=per_rio(h) [P]=per_rio(h) 2 CONTADOR=0; 3 G=h; seccion=xlsread('DATOS_CAUCE.xlsx', ,'hoja1' 'hoja1', ,'A2:B33' 'A2:B33'); ); 4 seccion=xlsread('DATOS_CAUCE.xlsx' 5 cotamax=max(seccion(:,2)); 6 cotamin=min(seccion(:,2)); 7 for for k=1:length(G) k=1:length(G) %para hacer valor por valor 8 vec=G(1,k); %cada valor del vector. 9 cota=vec+min(seccion(:,2)); if cota cota > cotamax 10 if error('COTA error('COTA MAYOR QUE LA MÁXIMA PERMITIDA') PERMITIDA') 11 12 end 13 if if cota cota < cotamin 14 error('COTA error('COTA MENOR QUE LA MÍNIMA PERMITIDA') PERMITIDA') 15 end 16 conta=length(seccion); 17 per=0; 18 ii=0; 19 B2=seccion(17,1)-seccion(16,1);
11
20 while while 1 1 21 i=conta/2+1+ii; %i va de 17 a 31 y j de 16 a 2 ii=ii+1; 22 j=conta+1-i; 23 B=seccion(i,1)-seccion(j,1); 24 25 hh=seccion(i+1,2)-seccion(i,2); 26 mizq=(seccion(j,1)-seccion(j-1,1))./hh; 27 mder=(seccion(i+1,1)-seccion(i,1))./hh; 28 h1=cota-seccion(i,2); 29 if if h1<=hh h1<=hh 30 per2=sqrt(((h1)^2)+((mizq.*h1)^2))+sqrt(((h1)^2)+((mder.*h1)^2)) 31 32 33 34 35 36 37 38 39 40 41 42 43 44
; CONTADOR=CONTADOR+1; P1=per+per2; P(CONTADOR)=P1; break end per=sqrt((hh^2)+((mizq.*hh)^2))+sqrt((hh^2)+((mder.*hh)^2))+per; %area cada tramo completo end end ku1=zeros(1,length(P)); for ku=1:length(P) for ku=1:length(P) ku1(1,ku)=P(1,ku)+B2; end P=(ku1)'; end
La principal diferencia con la función del área “area_rio” es la fórmula del
perímetro, que en este caso se calcula cada uno de los perímetros pero sin tomar en cuenta la base B ya que esta afecta la secuencia de cálculo y es por eso que cuando ya se tiene el perímetro de las paredes calculado y almacenado en memoria, se procede al sumar el ancho B únicamente al final del cálculo a cada uno de los elementos del vector de perímetros y con esto ya se tiene el perímetro final, este paso adicional que se ha implementado en este caso se hace constar en las líneas de código 39 hasta la 42. Finalmente todos los perímetros correspondientes a cada una de las alturas ingresadas se almacenan en la variable P que es de salida de resultados. Las fórmulas utilizadas para los perímetros tanto parcial completo como parcial incompletos son las las siguientes: Para perímetro parcial incompleto (solo paredes de altura incompleta). Entiéndase por alturas incompletas a aquellas que son menores o iguales a las alturas de un trapecio parcial, la fórmula utilizada es la siguiente:
12
La fórmula anterior es la que consta en la línea 30 del código de programación. Ahora para perímetros parciales completos completos se tiene:
Esta fórmula es la línea 36 y al a l parámetro per es el perímetro del paso anterior a nterior que se va sumando al siguiente paso hasta que el proceso se detenga en función del valor de h1 Con esto quedan terminadas las funciones de área y perímetro de las tres secciones y ahora la función descarga ya está lista para utilizarse, y lo haremos uso de esta función función el numeral 1.2 que se muestra a continuación. continuación.
ENUNCIADO DEL PROBLEMA 1.2 Desarrolle un programa (“script”) en MATLAB que llame la función de la curva de
descarga definida anteriormente y realice una figura con 3 gráficos presentados en la vertical, correspondientes a cada una de las alineas de 1.1 con pendientes longitudinales de 0,1%, 1% y 2%.
SOLUCIÓN AL PROBLEMA 1.2 1.2-1 Programa Script – curvas de descarga
El programa script que se desarrolló para obtener las curvas de descarga de cada una de las secciones con tres pendientes diferentes tiene un código de 59 líneas y es el que se muestra a continuación:
1 %PROGRAMA SCRIPT CURVAS DE DESCARGA CON PENDIENTE DADA: 2 clc 3 %SECCIÓN RECTANGULAR: 4 ni1=0; 5 deltha1=1; 6 nf1=30; 7 h_rec=ni1:deltha1:nf1; 8 H1=descarga(h_rec,area_rec(h_rec),per_rec(h_rec),0.1/100); 9 H2=descarga(h_rec,area_rec(h_rec),per_rec(h_rec),1/100); 10 H3=descarga(h_rec,area_rec(h_rec),per_rec(h_rec),2/100); 11 %SECCIÓN TRAPEZOIDAL: 12 ni2=0; 13 deltha2=1;
13
14 nf2=30; 15 h_tra=ni2:deltha2:nf2; 16 H4=descarga(h_tra,area_tra(h_tra),per_tra(h_tra),0.1/100); 17 H5=descarga(h_tra,area_tra(h_tra),per_tra(h_tra),1/100); 18 H6=descarga(h_tra,area_tra(h_tra),per_tra(h_tra),2/100); 19 %SECCIÓN DEL RÍO: 20 ni3=0; 21 deltha3=1; 22 nf3=30; 23 h_rio=ni3:deltha3:nf3; 24 H7=descarga(h_rio,area_rio(h_rio),per_rio(h_rio),0.1/100); 25 H8=descarga(h_rio,area_rio(h_rio),per_rio(h_rio),1/100); 26 H9=descarga(h_rio,area_rio(h_rio),per_rio(h_rio),2/100); 27 h1=ni1:deltha1:nf1; 28 h2=ni2:deltha2:nf2; 29 h3=ni3:deltha3:nf3; 30 Q1=H1; 31 Q1(:,2)=H2; 32 Q1(:,3)=H3; 33 Q2=H4; 34 Q2(:,2)=H5; 35 Q2(:,3)=H6; 36 Q3=H7; 37 Q3(:,2)=H8; 38 Q3(:,3)=H9; 39 %subplot===================== 40 subplot(3,1,1) 41 plot(h1,Q1), grid on legend('i=0.1%', ,'i=1%' 'i=1%', ,'i=2%' 'i=2%') ) 42 legend('i=0.1%' title('SECCIÓN RECTANGULAR') RECTANGULAR') 43 title('SECCIÓN 44 xlabel('h(m)' xlabel('h(m)') ) 45 ylabel('Q(m³/seg.)' ylabel('Q(m³/seg.)') ) 46 %subplot===================== 47 subplot(3,1,2) 48 plot(h2,Q2), grid on legend('i=0.1%', ,'i=1%' 'i=1%', ,'i=2%' 'i=2%') ) 49 legend('i=0.1%' title('SECCIÓN TRAPEZOIDAL') TRAPEZOIDAL') 50 title('SECCIÓN xlabel('h(m)') ) 51 xlabel('h(m)' ylabel('Q(m³/seg.)') ) 52 ylabel('Q(m³/seg.)' 53 %subplot===================== 54 subplot(3,1,3) 55 plot(h3,Q3), grid on 56 legend('i=0.1%' legend('i=0.1%', ,'i=1%' 'i=1%', ,'i=2%' 'i=2%') ) title('SECCIÓN DEL RÍO') RÍO') 57 title('SECCIÓN xlabel('h(m)') ) 58 xlabel('h(m)' 59 ylabel('Q(m³/seg.)' ylabel('Q(m³/seg.)') )
14
Descripción Descripció n del contenido del programa (“SCRIPT”) para curvas de descarga:
Este programa básicamente consiste en ingresar la función ya desarrollada (“descarga”) en el numeral anterior con la diferencia que para cada curva de descarga se le asigna un vector de alturas por ejemplo h=1:30 y adicionalmente se ingresa las funciones de área y perímetro para cada una de las secciones. En las líneas 4,5 y 6 se define el vector de alturas para la curva de descarga de la sección rectangular. En las líneas 12,13 y 14 se define el vector de alturas para la curva de descarga de la sección trapezoidal. En las líneas 20,21 y 22 se define el vector de alturas para la curva de descarga de la sección del río. En las líneas 8, 9 y 10 se calculan los caudales correspondientes a la sección rectangular. En las líneas 16, 17 y 18 se calculan los caudales correspondientes a la sección trapezoidal. En las líneas 24, 25 y 26 se calculan los caudales correspondientes a la sección del río. En la línea 27 se define el vector de alturas de la sección rectangular. En la línea 28 se define el vector de alturas de la sección trapezoidal En la línea 29 se define el vector de alturas de la sección del río. Desde la línea 30 a la línea 38 se intercambian las letras H por Q para los caudales correspondientes a cada una de las secciones y pendientes por separado. De la línea 39 a la línea 59 se procesa la información para ser graficada con los comandos plot y subplot que permitirán p ermitirán visualizar todos los gráficos en una sola ventana así como también las leyendas títulos y nom bres de los ejes a utilizar u tilizar encada caso. Las curvas de descarga que se obtienen para los datos ingresados son las que se muestra en la siguiente figura:
15
1.2-2 Gráficos de las curvas de descarga para las tres secciones dadas
3
x 10
4
SECCIÓN RECTANGULAR i=0.1% i=1% i=2%
) . g 2 e s / ³ m ( Q 1
0 0 15
5
x 10
10
4
15 h(m) SECCIÓN TRAPEZOIDAL
15
0
5
x 10
10
4
15 h(m) SECCIÓN DEL RÍO
30
20
25
30
25
30
i=0.1% i=1% i=2%
) . g 10 e s / ³ m ( Q 5
0
25
i=0.1% i=1% i=2%
) . g 10 e s / ³ m ( Q 5
0
20
0
5
10
15 h(m)
20
Fig. 5 Curvas de descarga para las tres secciones
Como se puede observar en las curvas de descarga para una misma altura dada, cuanto mayor sea la pendiente mayor será la velocidad del fluido y por lo tanto tendrá mayor capacidad de transporte de caudal y eso es justamente lo que se observa en las curvas de descarga cuanto mayor pendiente tengan mayor será el caudal que pase por esa misma sección. Además de todo esto se observa que por ejemplo e jemplo para una altura de 30 m comparando las tres secciones, el mayor caudal pasa por la sección trapezoidal con aproximadamente 12x10^4 m³/s mientras que por la sección del río pasan un caudal muy cercano con 10x10^4 m³/s y un caudal mucho más bajo por la sección rectangular para la misma altura solo pasan aproximadamente 2.5 x10^4 m³/s, este análisis se lo realiza para la pendiente del fondo de cauce de 2% y es por esa pendiente tan fuerte y altura de 30 metros que los caudales que se observan son muy grandes.
16
ENUNCIADO DEL PROBLEMA 1.3 Desarrolle un programa (“script”) en MATLAB que calcule y presente la curva de la
altura crítica en función del caudal para cada una de las alineas de 1.1.
SOLUCIÓN AL PROBLEMA 1.3 1.3-1 Programa Script – curvas de altura critica.
El programa (“script”) que se ha desarrollado para graficar la altu ra crítica en función del caudal se llama “ECUACION_CRITICO” y es el que se muestra a continuación:
1 %altura critica 2 clc 3 B=25; 4 hi=0; 5 paso=1; 6 hf=30; 7 h=(hi:paso:hf); 8 n=length(h); 9 g=9.81; 10 QQ1=ones(n,1); 11 QQ2=ones(n,1); 12 QQ3=ones(n,1); 13 for for i=1:n i=1:n 14 Q1=sqrt((g.*(area_rec(h(i))).^3)./(B)); Q1=sqrt((g.*(area_rec(h(i))).^3)./(B));%RECTANGULAR %RECTANGULAR 15 QQ1(i,1)=Q1; Q2=sqrt((g.*(area_tra(h(i))).^3)./(Darea_tra(h(i))));%TRAPEZOIDAL %TRAPEZOIDAL 16 Q2=sqrt((g.*(area_tra(h(i))).^3)./(Darea_tra(h(i)))); 17 QQ2(i,1)=Q2; Q3=sqrt((g.*(area_rio(h(i))).^3)./(Darea_rio(h(i))));%RÍO %RÍO 18 Q3=sqrt((g.*(area_rio(h(i))).^3)./(Darea_rio(h(i)))); 19 QQ3(i,1)=Q3; 20 end 21 QQ1 22 QQ2 23 QQ3 24 HH=h' 25 subplot(1,2,1);%GRAFICA subplot(1,2,1);%GRAFICA 1 26 plot(QQ1,HH,QQ2,HH,QQ3,HH), grid on on; ; 27 axis([0 10000 0 30]) legend('rectangular', ,'trapezoidal' 'trapezoidal', ,'río' 'río'); ); 28 legend('rectangular' title('hc=f(Q)'); ); 29 title('hc=f(Q)' 30 xlabel('Q(m³/s)' xlabel('Q(m³/s)'); ); ylabel('hc(m)'); ); 31 ylabel('hc(m)' subplot(1,2,2);%GRAFICA 2 32 subplot(1,2,2);%GRAFICA 33 plot(QQ1,HH,QQ2,HH,QQ3,HH), grid on on; ; 34 axis([0 1000 0 7]) legend('rectangular', ,'trapezoidal' 'trapezoidal', ,'río' 'río'); ); 35 legend('rectangular'
17
36 title('AMPLIACIÓN title('AMPLIACIÓN TRAMO INICIAL'); INICIAL'); 37 xlabel('Q(m³/s)' xlabel('Q(m³/s)'); ); ylabel('hc(m)'); ); 38 ylabel('hc(m)' El programa se basa fundamentalmente en obtener el caudal despejado de la ecuación de energía específica derivada e igualada a cero debido al criterio de mínimos, como se puede observar en esta ecuación se tiene la derivada del área en función de la altura y es por este motivo que se han creado funciones necesarias de las derivadas de la ecuación del área correspondiente a cada una de las secciones. El proceso es el siguiente:
Como sabemos la ecuación de energía es:
Derivando esta ecuación tenemos:
Y la ecuación a utilizar será la siguiente:
Como se puede observar en esta ecuación necesitamos las funciones de área las cuales ya tenemos creadas pero adicional a esto también necesitamos las funciones de las derivadas del área para pa ra la sección trapezoidal y sección de río, no se menciona la derivada de la sección rectangular debido a que es muy elemental y en este caso es una constante que la podemos ingresar directamente en la ecuación de caudal ya deducida. Los caudales correspondientes a las alturas críticas están constando en las líneas 14, 16 y 18 para secciones rectangular, trapezoidal y de rio respectivamente.
18
Como ya se mencionó en este programa se necesita las funciones que calculen la derivada del área para la sección trapezoidal y la sección del río y estas funciones son las siguientes: 1.3-2 función para obtener la derivada del área trapezoidal Esta función ha sido llamada llamada Darea_tra y es la que se muestra a continuación: continuación:
1 2 3 4 5 6
function DS=Darea_tra(h) function DS=Darea_tra(h) %PROGRAMA PARA CALCULAR LA DERIVADA DEL ÁREA B=25; m=2; DS=(B+2.*m.*h); end
Como se puede observar esta función permitirá calcular la derivada para cualquier altura dada y la ecuación que gobierna esto es la que se deduce a continuación: Como sabemos el área de una sección trapezoidal es:
Y su derivada es:
Y esta es la ecuación que se puede observar en la línea 5 del código anterior. Al igual que se hizo con el área aquí también se introducen introduce n los valores constantes de B y m dentro de la función mismos que se los puede observar en las líneas 3 y 4. 1.3-3 función para obtener la derivada del área del río Esta función ha sido llamada llamada Darea_rio y es la que se muestra a continuación: continuación:
1 2 3 4 5
function DS=Darea_rio(h) function DS=Darea_rio(h) %calcula la derivada del area del río %en base a la ecuacion ajustada del área S=c*h^n: DS=16.801733803014294.*1.413750007042020.*(h.^(1.4137500070420201)); end
En este código se puede observar una ecuación aproximada que gobierna la sección del rio misma que ha sido obtenida mediante en ajuste por mínimos cuadrados y posteriormente derivada, esta ecuación se jus tifica a continuación:
19
Con la ayuda de la función area_rio obtenemos datos de áreas para varios puntos, y una vez que tenemos estos puntos y con la ayuda de logaritmos linealizamos estos datos para posteriormente realizar el ajuste por mínimos cuadrados Los datos obtenidos son los siguientes: Tabla 2 Datos para ajuste de la ecuación de área del río
h(m) 30 28 26 24 22 20 18 16 14 12 10 8 6 4 2
S(m²) 2209.51 1974.79 1753.03 1543.71 1345.27 1156.54 980.63 819.89 671.99 535.75 410.91 298.64 200.36 116.98 51.19 SUMA=
Ln(h) 3.40119738 3.33220451 3 .33220451 3.25809654 3.17805383 3 .17805383 3.09104245 3 .09104245 2.99573227 2 .99573227 2.89037176 2.77258872 2.63905733 2.48490665 2.30258509 2.07944154 2 .07944154 1.79175947 1 .79175947 1.38629436 0.69314718 38.2964791 SX
Ln(S) 7.70052605 7.58821734 7.469101 7.34194389 7.20435002 7.05318807 6.88819522 6.70917018 6.51024346 6.28366763 6.01837421 5.69923883 5.30011575 4.76200298 3.9355442 96.4638788 SY
Ln(h)*Ln(S) 26.191009 25.2854921 24.3350521 23.3330929 22.2689517 21.1294631 19.9094449 18.6017696 17.1809057 15.6143275 13.8578187 11.851234 9.49653258 6.60153788 2.72791137 258.384543 SXY
(Ln(h))² 11.5681436 11.1035869 10.6151931 10.1000261 9.55454345 8.97441185 8.3542489 7.68724822 6.96462359 6.17476106 5.30189811 4.32407713 3.210402 1.92181206 0.48045301 106.335429 SXX
La ecuación de la recta de ajuste es de la forma:
La ecuación del área del rio es de la forma:
Por lo tanto:
y
b
ax
20
Dónde
El ajuste que se obtiene es el siguiente:
:
8 7.5
y = 1.4138x + 2.8215 R² = 0.9974
7 6.5 6 5.5 5 4.5 4 3.5
0
0.5
1
1.5
2
2.5
3
3.5
4
Fig. 6 Ajuste del área de río
Por lo tanto los valores de n y C son los siguientes: n= 1.413750007042020 c= 16.801733803014294 La ecuación que resulta de esto es la que necesitamos para la sección del río.
Y realizando el cálculo la derivada quedaría de la siguiente manera: dS/dh=16.801733803014294.*1.413750007042020.*(h.^(1.413750007042020-1))
Y esta es la fórmula que consta en la línea numero 4
21
También se comprueba este ajuste utilizando el programa que lo hemos llamado cal_area Y el código de programación es el siguiente: 1.3-4 Código para el ajuste del área de rio por mínimos cuadrados.
1 clear 2 clc seccion=xlsread('DATOS_CAUCE.xlsx', ,'hoja1' 'hoja1', ,'A2:B33' 'A2:B33'); ); 3 seccion=xlsread('DATOS_CAUCE.xlsx' 4 conta=length(seccion)/2-1; 5 hsum=0; 6 Ssum=0; for i=conta:-1:1 i=conta:-1:1 7 for 8 j=32-i; 9 B=seccion(j,1)-seccion(i+1,1); 10 hh=seccion(i,2)-seccion(i+1,2); 11 mizq=(seccion(i+1,1)-seccion(i,1))/hh; 12 mder=(seccion(j+1,1)-seccion(j,1))/hh; 13 area=hh.*(B+(mizq.*hh+mder.*hh)/2); 14 hsum=hh+hsum; 15 Ssum=Ssum+area; 16 hvec(i)=hsum; 17 Svec(i)=Ssum; 18 end 19 tabla=[hvec' Svec']; 20 tabla(:,3)=log(tabla(:,1)); 21 tabla(:,4)=log(tabla(:,2)); 22 tabla(:,5)=tabla(:,3).*tabla(:,4); 23 tabla(:,6)=(tabla(:,3)).^2; 24 suma=sum(tabla); 25 nn=length(tabla); 26 Sx=suma(1,3); 27 Sy=suma(1,4); 28 Sxy=suma(1,5); 29 Sxx=suma(1,6); 30 a=(nn*Sxy-Sx*Sy)/(nn*Sxx-Sx*Sx); 31 b=(Sxx*Sy-Sx*Sxy)/(nn*Sxx-Sx*Sx); 32 format long 33 n=a 34 c=exp(b) En este código básicamente lo que se hace es calcular los valores de a y b que sirven para calcular los valores de C y n que necesitamos para nuestra ecuación Y ejecutando este código lo que se obtiene es exactamente lo mismo que en la fig.6 a= 1.413750007042020 b= 2.821482083644570
22
n= 1.413750007042020 c= 16.801733803014294 Ahora que ya disponemos de todos los datos necesarios necesarios podemos ejecutar el programa ECUACION_CRITICO. Para una altura que va de 0 hasta 30 y se obtiene la siguiente gráfica: hc=f(Q)
30
AMPLIACIÓN TRAMO INICIAL
7 rectangular trapezoidal río
rectangular trapezoidal río 6
25
5 20 4 ) m ( 15 c h
) m ( c h
3 10 2
5
0
1
0
2000
4000 6000 Q(m³/s)
8000
10000
0
0
200
400
600
800
1000
Q(m³/s)
Fig. 7 Altura crítica en función de Caudal
En las gráficas se observa que la altura crítica es mayor para la sección rectangular con caudales grandes pero al realizar un zoom del tramo inicial se observa que para caudales menores a 600 m³/s la altura critica ya no es mayor para el canal rectangular sino que ahora es mayor para el canal de sección rio.
23
ENUNCIADO DEL PROBLEMA 1.4 Desarrolle un programa (“script”) en MATLAB que presente el gráfico tridimensional de la energía específica, H , en función de la altura del flujo, h, y del caudal, Q, para
la sección transversal de la alinea c) de 1.1. (NOTA: presente dos gráficos con vistas/ángulos de visualización diferentes)
SOLUCIÓN AL PROBLEMA 1.4 1.4-1 Código para grafica de energía especifica tridimensional.
El programa desarrollado para este fin se llama energía_3d y su código de programación es el siguiente:
1 %ENERGIA 3D 2 clc 3 g=9.81; 4 Hi=0.01; 5 paso=0.5; 6 Hf=40; 7 H=Hi:paso:Hf; 8 n=length(H); 9 HH=ones(n,1); 10 hh=ones(n,1); 11 QQ=ones(n,1); 12 for for ii=1:n ii=1:n 13 deltha=H(ii)/20; 14 h=0:deltha:H(ii); k=length(h); 15 for jj=1:k jj=1:k 16 for 17 QQ(ii,jj)=sqrt((H(ii)-h(jj)).*2.*g.*((area_rec(h(jj))).^2)); 18 HH(ii,jj)=H(ii); 19 hh(ii,jj)=h(jj); 20 end 21 end 22 HH; 23 hh; 24 QQ; subplot(2,2,1)%GRAFICO Nº 1 25 subplot(2,2,1)%GRAFICO 26 cs1=contour3(hh,HH,QQ,8); 27 clabel(cs1) 28 view([0,90]) xlabel('h(m³/s)'); ); 29 xlabel('h(m³/s)' ylabel('H(m)'); ); 30 ylabel('H(m)' zlabel('Q(m)'); ); 31 zlabel('Q(m)' 32 % xlabel('x'); 33 % ylabel('y'); 34 % zlabel('z');
24
35 %camorbit(37.5,55) 36 subplot(2,2,2)%GRAFICO subplot(2,2,2)%GRAFICO Nº 2 37 cs2=contour3(QQ,hh,HH,8); 38 clabel(cs2) 39 view([0,90]) 40 xlabel('Q(m³/s)' xlabel('Q(m³/s)'); ); 41 ylabel('h(m)' ylabel('h(m)'); ); 42 zlabel('H(m)' zlabel('H(m)'); ); 43 %camorbit(37.5,65) 44 subplot(2,2,3)%GRAFICO subplot(2,2,3)%GRAFICO Nº 3 45 surf(QQ,hh,HH) xlabel('Q(m³/s)'); ); 46 xlabel('Q(m³/s)' ylabel('h(m)'); ); 47 ylabel('h(m)' zlabel('H(m)'); ); 48 zlabel('H(m)' 49 subplot(2,2,4)%GRAFICO subplot(2,2,4)%GRAFICO Nº 4 50 mesh(QQ,HH,hh) xlabel('Q(m³/s)'); ); 51 xlabel('Q(m³/s)' ylabel('H(m)'); ); 52 ylabel('H(m)' zlabel('h(m)'); ); 53 zlabel('h(m)' La solución para crear el gráfico tridimensional de energía específica Al ver el comportamiento de la curva bidimensional de caudal en función de altura para una energía específica constante se puede observar claramente como la altura máxima que se puede alcanzar es un valor igual al de la energía específica que en este caso fue constante, entonces adoptando este criterio como criterio básico para la creación de la curva tridimensional se procedió de la siguiente manera: El primer paso fue crear una matriz de que contenga diversos valores de energía específica cuyos límites y paso de cálculo son los que se muestran en las líneas de programación 4, 5 y 6 para posteriormente crear crear el vector H en la línea línea 7. Luego para realizar las matrices correspondientes a caudal, altura y energía específica se procede a crear dos contadores uno llamado ii que representa las filas de cada matriz y otro llamado jj que representa las columnas de cada matriz. Luego de esto para cada una de las energías especificas se creara una fila de valores que van desde cero hasta un valor máximo de h que en este caso es la energía específica para esto nos ayudamos de un paso de calculo que lo llamamos deltha y en lugar de darle un valor constante le pedimos que nos divida cada altura especifica en un numero dado de espacios que en este caso adoptamos 20 los mismos que constan claramente en la línea l ínea número 13, este número de divisiones lo podemos reducir, pero este cambio se verá reflejado en la calidad de curva que se dibuje ya que a mayor número de tramos el paso de cálculo será menor y por lo tanto la exactitud de la curva será mucho mejor visualmente hablando.
25
Como ya se explicó el procedimiento adoptado solo queda por mostrar la ecuación resultante obtenida a partir de la energía específica y es esta misma ecuación la que se hace constar en la línea 17 y calculara cada uno de los caudales correspondientes correspondientes a las energías específicas específicas y rango de valores de h creado exclusivamente para cada valor de energía especifica. Ecuación utilizada para graficar la curva tridimensional de energía especifica:
De donde despejamos el caudal de la siguiente manera:
Y finalmente de las líneas 25 a 53 consta la representación gráfica de esta curva con sus s us diferentes vistas quedando finalmente de la siguiente forma: 9.4e+003 9.4e+003 8.23e+003 7.05e+003 5.88e+003 5.88e+003
35 30 25
35 30
) m20 ( H
5
10
13.2
10
17.6 5 4.4 8.79 0 0 2000
5 0
22
15
1.18e+003
10
30.7
) m20 ( h
2.35e+003
15
35.1
25
4.7e+003 3.53e+003
15
20 h(m³/s)
25
30
35
40
26.3
4000
6000 Q(m³/s)
8000
10000
40
30
30
) m20 ( H
) m20 ( h
10
10
0 40
0 40
20 h(m)
0
0
5000 Q(m³/s)
10000
15000
20 H(m)
0
0
5000 Q(m³/s)
10000
15000
Fig. 8 Energía especifica tridimensional
26
ENUNCIADO DEL PROBLEMA 1.5 Desarrolle un programa (“script”) en MATLAB que presente lo gráfico Q=Q(h), para
un momento total constante de 100 MN y para la sección transversal de la alinea c) de 1.1.
SOLUCIÓN AL PROBLEMA 1.5 1.5-1 Código utilizado para representar Q=Q(h) con momento constante: El programa desarrollado para este caso lo hemos h emos llamado momento_cte y su código de programación es el que se muestra a continuación:
1 %Momento constante: 2 clc 3 MI=100;%MN MI=100;%MN 4 MI=MI*1e6;%A MI=MI*1e6;%A NEWTON 5 gama=9800;%N/m³ gama=9800;%N/m³ 6 %gama=rho/g 7 %g=gravedad=9.81m/s² rho=1000;%Kg/m³ 8 rho=1000;%Kg/m³ 9 alfha=1; 10 h1=(0:21); 11 Q=(sqrt(((area_rio(h1))'))).*((MI12 13 14 15 16 17 18
gama.*((area_rio(h1))').*((center(h1))'))./(rho)); Q=Q'; h1=h1'; plot(Q,h1), grid on ylabel('h(m)' ylabel('h(m)'); ); xlabel('Q(m³/s)' xlabel('Q(m³/s)'); ); title('Q=Q(h)' title('Q=Q(h)'); ); legend('Caudal legend('Caudal para: MI=cte=100MN') MI=cte=100MN')
El código utilizado es muy sencillo donde básicamente ingresamos el valor del momento en la línea 3, ya que el momento está presentado en MN, este valor es convertido a Newton en la línea 4, la fórmula utilizada en este caso fue:
De donde:
̅
27
̅ Y despejando el caudal la ecuación es la siguiente:
̅ Pero como se puede observar en esta ecuación se tiene un problema y que aún no se ha resuelto mismo que tiene que ver con el centro de gravedad a cualquier altura de la sección transversal en análisis que en este caso es la sección del río. Para solucionar el problema antes expuesto es que se ha creado la función llamada centroide_rio misma que es capaz de obtener el centro de gravedad de la sección dada únicamente ingresando el valor de la altura y el código de programación que se ha desarrollado para esta función es el siguiente: function ccgg=centroide_rio(h) ccgg=centroide_rio(h) 1 function 2 clc 3 G=h; 4 seccion=xlsread('DATOS_CAUCE.xlsx' seccion=xlsread('DATOS_CAUCE.xlsx', ,'hoja1' 'hoja1', ,'A2:B33' 'A2:B33'); ); 5 cotamax=max(seccion(:,2)); 6 cotamin=min(seccion(:,2)); 7 conta=length(seccion); 8 ii=0; 9 pro_suma=0; 10 hh_suma=0; 11 CONTADOR=0; 12 CONTADOR2=0; 13 n_con=length(G); 14 vec=ones(n_con,1) for k=1:n_con k=1:n_con %para hacer valor por valor 15 for 16 vec(k,1)=G(1,k); %cada valor del vector. 17 cota=vec+cotamin; 18 if if cota cota > cotamax 19 error('COTA error('COTA MAYOR QUE LA MÁXIMA PERMITIDA') PERMITIDA') 20 end if cota cota < cotamin 21 if error('COTA error('COTA MENOR QUE LA MÍNIMA PERMITIDA') PERMITIDA') 22 23 end 24 while while 1 1 25 i=conta/2+1+ii; %i va de 17 a 31 y j de 16 a 2 26 ii=ii+1; j=conta+1-i; 27 B=seccion(i,1)-seccion(j,1); 28
28
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
hh=seccion(i+1,2)-seccion(i,2); mizq=(seccion(j,1)-seccion(j-1,1))./hh; mder=(seccion(i+1,1)-seccion(i,1))./hh; h1=cota-seccion(i,2); cg1=(2/3).*hh;%centroide cg1=(2/3).*hh;%centroide de los triángulos. cg2=hh./2;%centroide cg2=hh./2;%centroide del rectángulo central. aa1=(mizq.*(hh.^2))./2+(mder.*(hh.^2))./2; aa2=B.*hh; cg_parcial1=((cg1.*aa1+cg2.*aa2)./(aa1+aa2))+hh_suma;%con cg_parcial1=((cg1.*aa1+cg2.*aa2)./(aa1+aa2))+hh_suma; %con función área. area_parcial1=aa1+aa2; producto1=cg_parcial1.*area_parcial1; if h1==0 if h1==0 h_total=0; centros_g=0; centros_g2=0; break end if h1<=hh if h1<=hh cg11=(2/3).*h1;%centroide cg11=(2/3).*h1;%centroide de los triángulos. cg22=h1./2;%centroide cg22=h1./2;%centroide del rectángulo central. aa11=(mizq.*(h1.^2))./2+(mder.*(h1.^2))./2; aa22=B.*h1; cg_parcial2=((cg11.*aa11+cg22.*aa22)./(aa11+aa22))+hh_suma; area_parcial2=aa11+aa22; producto2=cg_parcial2.*area_parcial2; h_total=hh_suma+h1; centro_final1=(pro_suma+producto2)./(area_rio(h_total)); CONTADOR=CONTADOR+1; centros_g(CONTADOR)=centro_final1; CONTADOR2=CONTADOR2+1; centros_g2(CONTADOR2)=h_total-centros_g; break end pro_suma=pro_suma+producto1 hh_suma=hh_suma+hh; end ccgg=centros_g2;%centro ccgg=centros_g2;%centro de gravedad desde el espejo de agua. end end
En resumen el código antes mostrado no es más que un desarrollo adicional al ya creado y explicado código correspondiente para el cálculo del área de la sección de río, y es que en esta función las diferencias más importantes son por ejemplo el razonamiento y procedimiento de cálculo de centros de gravedad que consiste fundamentalmente en lo siguiente: Primero consideramos el cálculo de cada una de las áreas parciales en sub áreas
29
(dos triangulares extremas y un a rectangular central) y las multiplicamos por sus respectivos centros de gravedad tomando en cuenta como origen vertical el fondo del cauce, luego de esto procedemos con el mismo criterio de áreas parciales y vamos sumando el producto de cada una de las sub secciones hasta alcanzar la altura deseada, una vez sumados todos los productos de área por centros de gravedad procedemos a dividir para el área total que para este caso se puede utilizar la función creada para para este fin y de esta forma ya se dispone del centro de gravedad de la sección analizada a la altura deseada pero con el centro de gravedad desde el fondo y para corregir esto es que este centro de gravedad es restado de la altura total para obtener el complemento que en este caso será el centro de gravedad pero medido desde el espejo de aguas hacia abajo y que es justamente lo que estábamos buscando. buscando. Debido a que la función “centroide_rio” calcula el centroide para u na altura dada es que se ha desarrollado también la función “center” , misma que calcula el centroide para un conjunto de valores de h, o lo que es lo mismo un vector de alturas pero con la ayuda de la función anteriormente desarrollada(“centroide_rio”) , y es en esta nueva función se introduce la sentencia for para calcular los centroides uno a uno en función de cada altura y como resultado nos da un vector donde ya están calculados los centroides correspondientes a cada un a de las alturas ingresadas a esta función , el código de programación de la función “center” es el que se muestra a continuación:
1 function function cen_gra=center(h) cen_gra=center(h) 2 n=length(h); 3 g=h; for i=1:n i=1:n 4 for 5 hh=g(1,i); 6 cen1=centroide_rio(hh); 7 conta=i; 8 vector(conta)=cen1; 9 end 10 cen_gra=vector'; Como se puede observar los cálculos valor a valor que proporciona esta función son fácilmente calculados en la línea 6 del código reciente y el la sentencia for que hace esto paso a paso o valor por valor está en la línea 4 y finalmente presenta el vector de respuesta llamado cen_gra. Ahora con todos los parámetros de entrada entrada listos ya podemos pensar en graficar graficar la curva solicitada en este numeral y dicha grafica es la siguiente:
30
Q=Q(h)
25
Caudal para: MI=cte=100MN 20
15
10
5
0 -2
0
2
4
6
8 10 Q(m³/s)
12
14
16
18 x 10
5
Fig. 9 Gráfica Caudal para momento Cte=100MN
En la gráfica se puede observar de forma muy clara como, ya no es posible tener una altura mayor a 20.5 m ya que esta es la altura máxima que se puede alcanzar con el momento de 100MN, se observa que a partir de esa altura en adelante se tienen caudales negativos, lo cual no es posible sino que ya no se tiene suficiente fuerza para una altura mayor.
2. CÁLCULO DE CURVAS DE REMANSO CON CAUDAL CONSTANTE ENUNCIADO DEL PROBLEMA 2.1 Desarrolle un programa (“script”) en MATLAB que permita dibujar la curva de
remanso M2 para la sección transversal de la alinea a) de 1.1 y la pendiente longitudinal de 0.1%. Para ese fin solucione la siguiente EDO usando: dh d s
i
1
J 2
Fr
31
a) el método clásico de Runge-Kutta de 4ta orden (desarrolle el programa necesario); b) la función función ode45 del MATLAB; c) la función función ode23s del MATLAB. Presente los resultados en el mismo gráfico. Presente también un “zoom” de la región donde h está próximo de la altura crítica y comente los resultados.
SOLUCIÓN AL PROBLEMA 2.1 2.1-1 procedimientos para curvas de remanso M2: La curva de remanso tipo M2 fue grafica utilizando los tres métodos antes solicitados en un programa al cuál se lo ha llamado “curvas_remanso” y el código de programación para lograr todo esto es el siguiente:
1 clc hu=5;%m 2 hu=5;%m 3 n=0.047*(0.004)^(1/6); 4 g=9.81;%m/s² g=9.81;%m/s² i=0.001;%0.1% 5 i=0.001;%0.1% B=25;%m 6 B=25;%m Q=(1/n)*((B*hu)^(5/3))*sqrt(i)/((B+2*hu)^(2/3));%m/s³ %m/s³ 7 Q=(1/n)*((B*hu)^(5/3))*sqrt(i)/((B+2*hu)^(2/3)); hc=((Q^2)/(g*(B^2)))^(1/3)+0.001;%dato %dato de partida altura critica. 8 hc=((Q^2)/(g*(B^2)))^(1/3)+0.001; 9 %======================================== 10 %Método Runge Kutta 4 ORDEN s2=-1;%paso de cálculo 11 s2=-1;%paso 12 n2=2900; 13 hh=zeros(1,1); 14 ss=zeros(1,1); 15 suma=0; 16 hpartida=hc; for i2=1:n2 i2=1:n2 17 for k1=dh_ds(hpartida,i,n,g,B,hu,Q); 18 k2=dh_ds(hpartida+k1.*(s2/2),i,n,g,B,hu,Q); 19 k3=dh_ds(hpartida+k2.*(s2/2),i,n,g,B,hu,Q); 20 k4=dh_ds(hpartida+k3.*s2,i,n,g,B,hu,Q); 21 kk=k1+2*k2+2*k3+k4; 22 23 hpartida=hpartida+(1/6)*kk*s2; 24 hh(1,1)=hc; 25 hh(i2+1,1)=hpartida; suma=suma+s2; 26 ss(i2+1,1)=s2+suma; 27 28 end 29 s3=ss; 30 h3=hh; 31 %======================================== 32 dhds=@(s,h) (i-(Q^2)*(n^2)*((B+2.*h).^(4/3))/((B.*h).^(10/3)))/(133 34 35
(Q^2)/((B.*h).^2)/(g.*h)); [s1 h1]=ode45(dhds,[0:-1:-2900],hc);%PASO h1]=ode45(dhds,[0:-1:-2900],hc);%PASO DE CALCULO ES NEGATIVO [s2 h2]=ode23s(dhds,[0:-1:-2900],hc); subplot(2,1,1)
32
36 plot(s1,h1,'b--' plot(s1,h1,'b--'), ), grid on on, , hold on 37 xlabel('ABSCISAS xlabel('ABSCISAS EN METROS') METROS') ylabel('ALTURAS EN METROS') METROS') 38 ylabel('ALTURAS title('LONGITUD TOTAL DE LA CURVA M2') M2') 39 title('LONGITUD plot(s2,h2,'r'), ), grid on on, , hold on 40 plot(s2,h2,'r' 41 plot(s3,h3,'c--' plot(s3,h3,'c--'), ), grid on on, , hold on %Runge Kutta 4 42 legend('ode45' legend('ode45', ,'ode23s' 'ode23s', ,'R-4' 'R-4') ) 43 subplot(2,2,3)%ZOOM subplot(2,2,3)%ZOOM CERCA A LA ALTURA CRÍTICA PASO 44 [s1 h1]=ode45(dhds,[0:-1:-10],hc);%PASO h1]=ode45(dhds,[0:-1:-10],hc);%PASO DE CALCULO 45 [s2 h2]=ode23s(dhds,[0:-1:-10],hc); plot(s1,h1,'b--'), ), grid on on, , hold on 46 plot(s1,h1,'b--' xlabel('ABSCISAS EN METROS') METROS') 47 xlabel('ABSCISAS ylabel('ALTURAS EN METROS') METROS') 48 ylabel('ALTURAS plot(s2,h2,'r'), ), grid on on, , 49 plot(s2,h2,'r' 50 plot(s3,h3,'c' plot(s3,h3,'c'), ), grid on on, , 51 axis([-10 0 3.4 3.70]) title('ZOOM HASTA ABSCISA= -10 metros') metros') 52 title('ZOOM legend('ode45', ,'ode23s' 'ode23s', ,'R-K4' 'R-K4') ) 53 legend('ode45' subplot(2,2,4)%ZOOM CERCA A LA ALTURA CRÍTICA PASO 54 subplot(2,2,4)%ZOOM h1]=ode45(dhds,[0:-1:-10],hc);%PASO DE CALCULO 55 [s1 h1]=ode45(dhds,[0:-1:-10],hc);%PASO 56 [s2 h2]=ode23s(dhds,[0:-1:-10],hc); 57 plot(s1,h1,'b--' plot(s1,h1,'b--'), ), grid on on, , hold on xlabel('ABSCISAS EN METROS') METROS') 58 xlabel('ABSCISAS ylabel('ALTURAS EN METROS') METROS') 59 ylabel('ALTURAS plot(s2,h2,'r'), ), grid on on, , 60 plot(s2,h2,'r' 61 plot(s3,h3,'c' plot(s3,h3,'c'), ), grid on on, , 62 axis([-2 0 3.40 3.6]) 63 title('ZOOM title('ZOOM HASTA ABSCISA= -2 metros') metros') legend('ode45', ,'ode23s' 'ode23s', ,'R-K4' 'R-K4') ) 64 legend('ode45'
A ES NEGATIVO
B ES NEGATIVO
Descripción del contenido del programa (“SCRIPT”) para curvas de descarga:
En el código se muestra los tres métodos solicitados para la obtención de la curva de remanso mismos que son: Método de Runge Kutta 4 orden, función de matlab ode45 y también la función de matlab ode23s En las líneas líneas 1 a 8 se ingresan fundamentalmente fundamentalmente los los datos de partida en este caso como en el enunciado no se dispone de ningún tipo de información como caudal o altura sino que únicamente tenemos disponible la pendiente que hemos de utilizar y es justamente por eso que para efectos de obtener resultados visibles se adopta una altura de 5 metros que consta en la línea 2 y con el resto de datos se procede a calcular el caudal mismo que consta en la línea 7 y la fórmula utilizada es la siguiente:
33
Una vez realizado este cálculo se procede a calcular la altura crítica misma que es calculada con la siguiente formula:
Esta fórmula consta en la línea 8 pero como se puede observar en la línea 8 al resultado obtenido con esta fórmula se le ha sumado 0.001 m (1mm) estratégicamente debido a que si no se lo hace de esta manera cualquier método de solución de ecuaciones diferenciales que adoptemos no pasara del primer paso ya que al utilizar la altura critica como dato de partida el número de Froude tomará un valor de 1 ya que la la ecuación del número de Froude correspondería a flujo subcrítico y en este caso solo se anularía el denominador de la ecuación diferencial de flujo gradualmente variado (FGV),este fenómeno lo podemos observar claramente en la siguiente ecuación:
Como se puede observar en la ecuación anterior si el número de Froude vale 1 entonces el denominador se anula, y es por esta razón en la línea 8 a la altura crítica se la ha sumado 0.001 y con este pequeño p equeño cambio el denominador ya no se anula pero se obtiene una pendiente muy grande en el primer tramo debido a que solo estamos a 0.001m de distancia horizontal de espacio del punto crítico dónde la pendiente seria infinita, pero desde ahí a hí ya podemos empezar a calcular y dibujar la curva de remanso con los métodos ya mencionados . Desde la línea 11 hasta la línea 30 se tiene programado el método de Runge Kutta de 4º orden, en el cual básicamente se ha utilizado las siguientes formulas: Dónde:
Y posteriormente en las línea 32 consta la ecuación diferencial a resolverse por los métodos de MATLAB ode45 Y ode23s, la solución por ode45 consta en la línea 33 y la solución por ode23s consta en la línea 34.
34
Finalmente de las las líneas 35 a 64 consta los los códigos básicos para obtener las las gráficas de los resultados obtenidos por los lo s diferentes métodos. Realizando y ordenando los gráficos de curva de remanso M2 se obtiene la siguiente gráfica:
LONGITUD TOTAL DE LA CURVA M2
5 S O4.5 R T E M N E 4 S A R U T L 3.5 A
ode45 ode23s R-4
3 -3000
3.7
-2500
-2000
-1500 ABSCISAS EN E N METROS
ZOOM HASTA ABSCISA= -10 metros
-8
-6 -4 ABSCISAS EN METROS METROS
-2
0
ode45 ode23s R-K4
S O3.55 R T E M N E 3.5 S A R U T 3.45 L A
3.45 3.4 -10
-500
ZOOM HASTA ABSCISA= -2 metros
ode45 ode23s R-K4
3.65
S O R T 3.6 E M N E 3.55 S A R 3.5 U T L A
-1000
0
3.4 -2
-1. 5
-1 -0. 5 ABSCISAS EN METROS METROS
0
Fig. 10 Curva de Remanso Tipo M2
En la Fig.10 se puede observar claramente lo que ya se mencionó con anterioridad respecto a la pendiente exagerada que se obtendrá cuando el paso de cálculo este muy cercano al punto crítico está pendiente exagerada se muestra en al grafica completa donde no se ha realizado el zoom y más específicamente en la línea obtenida con el metodo de Runge Kutta de 4º orden que consta en color cian, en la ampliación del tramo inicial se muestra como los métodos son mejores o perores según el error que se obtenga y es en el siguiente orden: Primero el método de MATLAB ode23s es el que menos error tiene y la pendiente baja hace que se ajuste mas a la curva real. Segundo el método de MATLAB ode45 tiene un error mayor que el método ode23s pero mucho menor que el método de runge Kutta de 4º orden
35
Finalmente el tercero que es el método de Runge Kutta de 4º orden el error en el tramo inicial es muy grande y se aleja casi 1metro de la realidad en el primer punto lográndose disipar este error a lo largo de la curva conforme se aleja del punto crítico el error es menor pero aun así es muy grande el error que se cometería utilizando este método y sobre todo con pasos de cálculo muy mu y grandes a modo de ejemplo ahora se mostrara en el siguiente grafico lo que pasa si en lugar de utilizar un paso de cálculo de 1 m utilizamos un paso de cálculo de 0.10 m LONGITUD TOTAL DE LA CURVA M2
5 S O 4.5 R T E M N E 4 S A R U T L 3.5 A
ode45 ode23s R-4
3 -3000
-2500
-2000
-1500 ABSCISAS EN METROS
ZOOM HASTA ABSCISA= -10 metros
3.7 3.65
-8
-6 -4 ABSCISAS EN METROS
-2
0
ode45 ode23s R-K4
S O 3.55 R T E M N E 3.5 S A R U T 3.45 L A
3.45
3.4 -10
-500
ZOOM HASTA ABSCISA= -2 metros ode45 ode23s R-K4
S O R T 3.6 E M N E 3.55 S A R U 3.5 T L A
-1000
0
3.4 -2
-1. 5
-1 -0. 5 ABSCISAS EN METROS
0
Fig. 11 Curva de Remanso con paso de cálculo=0.10m
Ahora ya se puede observar como al realizar este e ste cambio en e n el paso pa so de cálculo el error más alto disminuye considerablemente, en este caso disminuye de 1 m a 0.05 m , entonces se puede concluir que los mejores métodos en estos caso son la ode45 y ode23s ya que el error será muy pequeño, aunque el método de Runge Kutta de 4º orden se lo puede mejorar considerando un tramo inicial de 1000 pasos de cálculo de 0.001m pero esto es mucho calculo comparado con la ode 45 o con la ode23s pero de todas maneras se muestra el resultado que se obtendría en la curva proveniente del método de Runge Kutta de 4º orden sobre todo observando el zoom del tramo inicial que quedaría de la siguiente manera:
36
LONGITUD TOTAL DE LA CURVA M2
5
ode45 ode23s R-4
S O 4.5 R T E M N E 4 S A R U T L 3.5 A
3 -3000
-2500
-2000
-1500 ABSCISAS EN METROS METROS
ZOOM HASTA ABSCISA= -10 metros
3.7
3.55 3.5 3.45 -8
-6 -4 ABSCISAS EN ME TROS TROS
0
-2
ode45 ode23s R-K4
S O 3.55 R T E M N E 3.5 S A R U T 3.45 L A
3.6
3.4 -10
-500
ZOOM HASTA HASTA AB SCISA= -2 metros ode45 ode23s R-K4
3.65
-1000
3.4 -2
0
-1.5
-1 -0.5 ABSCISAS EN ME TROS TROS
0
Fig. 12 Curva de remanso con paso=0.001m en RK4
Ahora ya se observa cómo se ha eliminado la mayor parte de er ror y en este caso el error no superara de 1 o 2 cm, lo cual en la práctica ya es muy bueno.
ENUNCIADO DEL PROBLEMA 2.2 Desarrolle un programa (“script”) en MATLAB que permita presentar el andamiento
de la superficie libre en el siguiente canal
A i AB =
0,006
B i BD =
0,0003
D
El canal AD tiene sección transversal rectangular, con 4 m de ancho, y es de hormigón liso ( K = = 80 m 1/3/s). El caudal que pasa en el canal de A para D es de 11,2 m3/s y los tramos AB y BD son largos lo suficiente (por ejemplo 5 km) para que el régimen uniforme sea atingido. Considere que la cota del fondo en A es igual a
37
200,00. SOLUCIÓN AL PROBLEMA 2.2 2.2-1 Planteamiento del problema:
En este caso no se puede trabajar con energías específicas ya que se trata de un fenómeno local en el cual es mucho mejor trabajar con los momentos y en este caso calcularemos un momento aguas arriba y otro momento aguas abajo para verificar si la fuerza con la que baja el agua es suficiente para empujar el resalto aguas abajo de la la sección de cambio o no lo es, en este caso caso el resalto resalto se formara aguas arriba, para esto utilizamos la ya conocida ecuación de momento:
̅
Como se trata de un canal rectangular la ecuación queda de la siguiente forma:
Esta es la fórmula que nos sirve para comprobar que nuestro cálculo es correcto y consta claramente el cálculo de los momentos aguas arriba y aguas debajo de la sección de cambio en las líneas 21 y 22 del código llamado “andamiento” mismo que se lo muestra a continuación:
1 %DATOS PARA EL ANDAMIENTO: 2 clc 3 B=4; 4 g=9.81; 5 K=80; 6 n=1/K; i1=0.006;%pendiente del primer tramo 7 i1=0.006;%pendiente 8 Q=11.2; 9 fun=@(h) Q*n/sqrt(i1)-(((B*h)^(5/3))/((B+2*h)^(2/3))); 10 PROX=5; 11 h_u1 = fzero(fun,PROX); 12 h_c1=((Q.^2)./(g.*(B.^2))).^(1/3); 13 i2=0.0003;%pendiente i2=0.0003;%pendiente del segundo tramo 14 fun=@(h) Q*n/sqrt(i2)-(((B*h)^(5/3))/((B+2*h)^(2/3))); 15 h_u2 = fzero(fun,PROX); 16 h_c2=((Q.^2)./(g.*(B.^2))).^(1/3); 17 %=========================== 18 %CÁLCULO DE MOMENTOS: GAMA=9800;%N/m³ 19 GAMA=9800;%N/m³ 20 RHO=1000;%Kg/m³ RHO=1000;%Kg/m³ 21 MI1=GAMA*B*((h_u1)^2)/2+RHO*(Q^2)/(B*h_u1) 22 MI2=GAMA*B*((h_u2)^2)/2+RHO*(Q^2)/(B*h_u2)
38
23 %========================== 24 %CURVA DE REMNASO S1 i=i1;%0.1% 25 i=i1;%0.1% 26 hu=h_c1+0.0001; 27 dhds=@(s,h) (i-(Q^2)*(n^2)*((B+2.*h).^(4/3))/((B.*h).^(10/3)))/(128 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
(Q^2)/((B.*h).^2)/(g.*h)); cont=0; while 1 while 1 cont=cont+1; [s2 h2]=ode23s(dhds,[0:0.01:cont],hu); maximo=max(h2); if maximo>=h_u2 if maximo>=h_u2 break end end cont2=0; while 1 while 1 cont2=cont2+1; if h2(cont2)>=h_u2 if h2(cont2)>=h_u2 limitemax=s2(cont2); break end end [s2 h2]=ode23s(dhds,[0:limitemax],hu); s1=[0;100]; h1=[h_u1;h_u1]; incremento1=s1(2); s2=s2+incremento1; s3=[s2(length(s2));s2(length(s2))+100]; h3=[h_u2;h_u2]; s4=[s2(length(s2));s2(length(s2))]; h4=[0.5;2.5]; h1=200+h1; h2=200+h2; h3=200+h3; h4=200+h4; subplot(3,1,1) plot(s2,h2,'b' plot(s2,h2,'b'), ), grid on on, , hold on plot(s4,h4,'r' plot(s4,h4,'r'),grid ),grid on on, , hold on plot(s1,h1,'b' plot(s1,h1,'b'),grid ),grid on on, , hold on plot(s3,h3,'b' plot(s3,h3,'b'),grid ),grid on on, , hold on xlabel('ABSCISAS xlabel('ABSCISAS EN METROS') METROS') ylabel('ALTURAS ylabel('ALTURAS EN METROS') METROS') title('SUPERFICIE title('SUPERFICIE LIBRE CON CURVA TIPO S1') S1') legend('Espejo legend('Espejo de agua', agua','Sección de Cambio') Cambio')
39
Descripción Descripció n del contenido del programa (“SCRIPT”) para andamiento de superficie
libre:
En el código anterior se observa en las primeras líneas únicamente los datos del problema a resolver, luego con estos datos se calcula las alturas uniformes aguas arriba que conocemos como h_u1 y la altura uniforme aguas abajo que se le conoce como h_u2, mismas que constan en las líneas 11 y 15 respectivamente, todo esto fue posible únicamente evaluando para cero las funciones de las líneas 9 y 14. El las líneas 19 a 22 se observa claramente el cálculo de los momentos de estos dos tramos en análisis y los resultados resu ltados obtenidos son los siguientes: MI1 =
5.4402e+004
MI2 = 9.5182e+004 Como se puede observar en los resultados la fuerza aguas abajo de la sección de cambio es mayor mayor que la fuerza aguas arriba, razón por la cual la la altura en la sección de cambio será igual a la altura uniforme en suscritico es decir aguas abajo, y a partir de la sección de cambio hacia arriba se tiene la curva “S 1” misma que es graficada utilizando la función ode23s que ya fue detallada en el numeral 2.1 y fue calculada a partir de la altura crítica hacia abajo como consta en las línea 45. El cálculo que se realizó en las líneas 29 a 44 fue únicamente para encontrar la longitud con la que se logra alcanzar la altura uniforme del tramo 2. Las líneas 46 y 47 definen la superficie uniforme aguas arriba mientras que en las líneas 50 y 51 se define la superficie libre aguas abajo. Desde las líneas 54 a 57 se incremente la cota del fondo del cauce de 200m, y se asume igual para todas las alturas ya que la pendiente del fondo es muy baja y los resultados no varían mayormente al hacer este a rtificio. Finalmente de las líneas 59 a 66 se unen todos los tramos de grafico en uno solo para visualizar todo el andamiento, con excepción del tramo inicial de la curva ya que parte desde la altura crítica y el tramo discontinua que se observa es justamente donde se produce la mayor turbulencia debido al resalto hidráulico generado en este tramo, y la gráfica obtenida es la que se muestra a continuación. S 202.5 O R T 202 E M N E 201.5 S A R 201 U T L A 200.5
SUPERFICIE LIBRE CON CURVA TIPO S1 Espejo de agua Sección de Cambio
0
50
100
150 200 ABSCISAS E N METROS
Fig. 13 Andamiento de la superficie libre
250
300
350
40
Contenido ESCUELA DE INGENIERÍA CIVIL ................................. ................................................. ................................. ................................. ................................. .........................1 ........1 1. CÁLCULO HIDRÁULICO BÁSICO ............................ ............................................. .................................. .................................. ................................. ........................2 ........2 ENUNCIADO DEL PROBLEMA 1.1 .................................. .................................................. ................................. .................................. .................................2 ................2 SOLUCIÓN AL PROBLEMA 1.1 ....................................... ....................................................... ................................. .................................. .................................3 ................3 1.1-1 Función Descarga: ................................ ................................................ ................................. .................................. .................................. ...........................3 ..........3 1.1-1.1 Cálculo de coeficiente n de Manning: M anning: ......................................... ......................................................... ................................. ...................3 ..3 1.1-2 Cálculo del área para la sección rectangular: ...................................... ....................................................... ..............................4 .............4 Descripción del contenido de la función area_rec .................................. .................................................. ................................. ...................4 ..4 1.1-3 Cálculo del perímetro para la sección rectangular: ............................ ............................................. ..............................4 .............4 1.1-4 Cálculo del área para la sección trapezoidal: .............................................. ............................................................... ......................5 .....5 Descripción del contenido de la función area_tra ..................................... ...................................................... .................................6 ................6 1.1-5 Cálculo del perímetro para la sección trapezoidal:........................................... trapezoidal:...........................................................6 ................6 1.1 – 6 Cálculo del área para la sección s ección del río: ....................... ....................................... ................................. .................................. ...................7 ..7 Descripción del contenido de la función area_río ............................... ................................................ ................................... ......................8 ....8 1.1 – 7 Cálculo del perímetro para la sección s ección del río: ................... ..................................... .................................. ...........................11 ...........11 ENUNCIADO DEL PROBLEMA 1.2 .................................. .................................................. ................................. .................................. ...............................13 ..............13 SOLUCIÓN AL PROBLEMA 1.2 ....................................... ....................................................... ................................. .................................. ...............................13 ..............13 1.2-1 Programa Script – curvas de descarga ..................... ...................................... .................................. ................................. ......................13 ......13 Descripción del contenido del programa (“SCRIPT”) para curvas de descarga: ........................15 ........................15
1.2-2 Gráficos de las curvas de descarga para las tres secciones dadas .......................... ...................................16 .........16 ENUNCIADO DEL PROBLEMA 1.3 .................................. .................................................. ................................. .................................. ...............................17 ..............17 SOLUCIÓN AL PROBLEMA 1.3 ....................................... ....................................................... ................................. .................................. ...............................17 ..............17 1.3-1 Programa Script – curvas de altura critica. .................... .................................... ................................. .................................. .................17 17 1.3-2 función para obtener la derivada del área trapezoidal ........................................... ...................................................19 ........19 1.3-3 función para obtener la derivada del área del río ............................... ................................................ ............................19 ...........19 1.3-4 Código para el ajuste del área de rio por mínimos cuadrados. ....................................... .......................................22 22 ENUNCIADO DEL PROBLEMA 1.4 .................................. .................................................. ................................. .................................. ...............................24 ..............24 SOLUCIÓN AL PROBLEMA 1.4 ....................................... ....................................................... ................................. .................................. ...............................24 ..............24 1.4-1 Código para grafica de energía especifica tridimensional. ...................................... ..............................................24 ........24 La solución para crear el gráfico tridimensional de energía específica.......................... específica.....................................25 ...........25 ENUNCIADO DEL PROBLEMA 1.5 .................................. .................................................. ................................. .................................. ...............................27 ..............27
41
SOLUCIÓN AL PROBLEMA 1.5 ....................................... ....................................................... ................................. .................................. ...............................27 ..............27 1.5-1 Código utilizado para representar Q=Q(h) con momento constante: .......................... .............................27 ...27 2. CÁLCULO DE CURVAS DE REMANSO CON CAUDAL CONSTANTE ................................. .................................................. .................31 31 ENUNCIADO DEL PROBLEMA 2.1 ................................. ................................................. .................................. .................................. ...........................31 ...........31 SOLUCIÓN AL PROBLEMA 2.1 ........................................ ......................................................... .................................. ................................. .........................32 .........32 2.1-1 procedimientos procedimientos para curvas de remanso M2: .................... ...................................... .................................. ...........................32 ...........32 ........................33 Descripción del contenido del programa (“SCRIPT”) para curvas de descarga: ........................33 ENUNCIADO DEL PROBLEMA 2.2 .................................. .................................................. ................................. .................................. ...............................37 ..............37 SOLUCIÓN AL PROBLEMA 2.2 ........................................ ......................................................... .................................. ................................. .........................38 .........38 2.2-1 Planteamiento Planteamiento del problema: ....................... ........................................ ................................. ................................. .................................. .................38 38 Descripción del contenido del programa (“SCRIPT”) para andamiento de superficie libre: .....40
Figuras Fig. 1 Sección Rectangular .................................. ................................................... .................................. ................................. ................................. ..............................4 .............4 Fig. 2 Sección Trapezoidal ............................... ................................................ .................................. ................................. ................................. .................................5 ................5 Fig. 3 Sección del río ............................... ................................................ ................................... .................................. ................................. ................................. ........................7 ........7 Fig. 4 Sección parcial del río ............................... ................................................ .................................. ................................. ................................. ............................10 ...........10 Fig. 5 Curvas de descarga para las tres secciones............................. secciones.............................................. .................................. ...............................16 ..............16 Fig. 6 Ajuste del área de río .......................... ........................................... .................................. ................................. ................................. .................................. .................21 21 Fig. 7 Altura crítica en función de Caudal ........................... ............................................. .................................. ................................. ............................23 ...........23 Fig. 8 Energía especifica tridimensional ................................. ................................................. ................................. .................................. .........................26 ........26 Fig. 9 Gráfica Caudal para momento Cte=100MN ..................................... ...................................................... ................................. ......................31 ......31 Fig. 10 Curva de Remanso Tipo M2 ............................... ................................................ .................................. ................................. ................................. .................35 35 Fig. 11 Curva de Remanso con paso de cálculo=0.10m .................. .................................. ................................. .................................. .................36 36 Fig. 12 Curva de remanso con paso=0.001m en RK4 ................................. .................................................. ................................. ......................37 ......37 Fig. 13 Andamiento de la superficie libre ....................................... ........................................................ ................................. ................................. .................40 40
Tablas Tabla 1 Datos del d el cauce de un río ................................. .................................................. ................................. ................................. .................................. ...................2 ..2 Tabla 2 Datos para ajuste de la l a ecuación de área ár ea del río .............................................. .............................................................. ...................20 ...20
42
