CAPyAH 2013 – San Juan (Argentina) Clasificación: (Proyecto de Presas)
MODELAMIENTO NM!"ICO DE IN#ILT"ACI$N EN P"ESAS CON EL M!TODO DE %OL&MENES #INITOS' CASO P"ESA CCHOESE"A A*ACCHO
Autor: Empresa/Organismo: Email:
Isaías Antonio Vilca Tueros Proyecto Especial ierra Centro ur ! "inisterio #e Agricultura $ni%ersi#a# &acional an Cristó'al #e uamanga Per* ci%ilca+gmail,com
Resumen El presente trabajo de investigación comprende el modelamiento numérico del fenómeno de infiltración en medios porosos, en particular en problemas de infiltración en presas, gobernada por la ecuación diferencial de flujo potencial (Ec. Laplace), en los límites de validez de la Ley de Darcy, donde la variable principal son las alturas piezométricas y el gradiente hidráulico, que son factores fundamentales para la ocurrencia del fenómeno. La técnica numérica principal en esta modelación es el método de volúmenes finitos (MVF); sin embargo se aborda también la modelación con el método de elementos finitos (MEF); la discretización geométrica utilizada es la malla triangular no estructurada, generadas por la triangulación de Delaunay. En síntesis el trabajo comprende: la formulación matemática, modelación numérica, programación computacional del modelo con MVF y utiliza un Software comercial para la modelación con MEF . La validación de estos modelos se hace con los datos históricos de medición en piezómetros de la Presa Cuchoquesera en Ayacucho - Perú , la cual en su sección de máxima altura alcanza 42m. y está diseñada para almacenar un volumen de 80MMC.
Palabras claves : Infiltración en Presas, Modelamiento de Infiltración, Método de Volúmenes
Finitos. El objetivo
El objetivo de esta investigación es hacer una formulación matemática del fenómeno de infiltración en medios porosos y particularmente en presas, esto significa determinar las alturas piezométricas en el dominio físico como variable principal y como secundarias las alturas de y caudal de filtración . Establecer un modelo presión , gradiente hidráulico , velocidades velocidades de filtración y que permitan el análisis del movimiento del agua a través de las presas (medios porosos) con el método de volúmenes finito , hacer su programación con la finalidad de validar el modelo utilizando los datos piezométricos registrados históricamente en la Presa Cuchoquesera. Planteamiento del problema Definición. El fenómeno de infiltración a través de las presas de tierra , y la fundación por debajo de la cimentación en presas de gravedad o similares , compromete las condiciones de estabilidad
frente al deslizamiento en masa y la estabilidad superficial frente al fenómeno de arrastre y tubificación, ocasionando los problemas como el gasto excesivo del agua en el embalse y la rotura de la presa , en casos más críticos. La ampliación del campo de aplicabilidad de los métodos numéricos como elementos finitos (MEF) y volúmenes finitos (MVF), ya sea en el campo de mecánica de estructuras, mecánica de medios sólidos, mecánica de medios continuos y otros , es precariamente desarrollado en nuestro
medio; sin embargo cualquier oficina técnica, por pequeña que sea, las usa; pero hoy en día resulta relativamente frecuente que se lleven a cabo cálculos por personal que desconoce casi absolutamente los fundamentos de los métodos que utiliza el Software y sus limitaciones, por tanto, es incapaz de evaluar la bondad de los resultados que está obteniendo. Planteamiento. Siendo el suelo, una materia heterogénea conformada por una fase sólida y los espacios vacíos (poros ), entonces el fenómeno ocurre netamente en un medio poroso o discontinuo . El agua que fluye a través del él desarrolla un flujo laminar , por tanto aplicable la Ley de Darcy , dependiente solo del gradiente hidráulico y propiedades del suelo y agua. En tal sentido, al fenómeno de infiltración gobierna la ecuación de flujo potencial (ecuación diferencial de segundo orden en derivadas parciales ), tanto en estado transitorio como estacionari o, en este último reduciéndose a la ecuación de Laplace ; la solución de esta ecuación permiten determinar la alturas piezométricas como variable principal y la secundarias la velocidades de filtración , caudales de filtración , las alturas de presión , etc.
Por otra parte la ecuación planteada no tiene una solución analítica, por tanto, nuestra propuesta es resolver con el método de volúmenes finitos (MVF), una de las técnicas numéricas que ha sido concebidas precisamente para resolver ecuaciones diferenciales en forma aproximada; en tanto que su precisión de sus resultados es dependen de la modelación adecuada del fenómeno. Sumando a esto el avance de la tecnología computacional y el desarrollo de las técnicas de discretización del dominio, esta metodología resultan poderosa técnica de cálculo en el campo de la ingeniería. La discretización del dominio será con una malla triangular construido con el método de triangulación de Delaunay , a partir del cual se construirá además los volúmenes finitos o volúmenes de control con él método de las medianas para la integración numérica utilizando el Método de la Divergencia de Gauss . Formulación matemática del problema.
El suelo, como un medio heterogéneo, contiene entre si, aire, agua y la matriz sólida , para un volumen elemental representativo (VER), sea Vol w volumen del agua y Vol T volumen total o VER. Para el fin ,y , entonces ; en el fenómeno es importante
Volw Vol β Vol βVol β m → dβ m dμ →μ hzγ
considerar la variación del contenido volumétrico del agua ( ), respecto a la presión de poros de agua en el suelo µw, entonces
; por la Ley de Bernuolli
, sea =µwγ. Finalmente dβ= d(h-z). siendo µw pendiente la curva presión de poros y el contenido volumétrico de agua en el suelo.
Figura 1: Balance en el volumen elemental representativo (VER) del suelo.
Considerando la figura 1 un diferencial de volumen (VER) del suelo, a través de cual fluye agua desarrollando un flujo laminar , entonces la ecuación gobernante del fenómeno se deduce a partir de la ley de balance de masas (ecuación de continuidad ) y la Ley de Darcy (figura 1), Según nuestro propósito, la ecuación reducido a su forma bidimensional (2D) es (% (1) ()
!!" #$" !!"%& !'! #$' !!'%& * (%
+,-./ *0 . $1 $0 / ! ! 0+,- ! !
y en su forma vectorial () !" !' !" !'.Siendo x,y,z posición, t tiempo, h altura piezométrica, kx, kz coeficiente de permeabilidad del suelo, Q condiciones de borde, =m w γ . El completo desarrollo deductivo se encuentra en la tesis del cual es síntesis este artículo. También las contribuciones en trabajos similares podemos señalar a LAMBE [1], TERZAGHI [2] MIGUEL TOLEDO [3], KUMAR [4], de cuyos investigaciones es deudor el presente en muchos aspectos. Modelamiento con volúmenes finitos En esencia para formulación del método de volúmenes finitos (MVF) es el principio de
conservación de una determinada cantidad física expresada por la ecuación gobernante (1) sobre cualquier volumen finito o volumen de control , [5, 6, 7, 3, 8, 9, 10, 11].
Figura 2: Construcción de los volúmenes finitos Una forma de obtener los volúmenes de control con una malla triangular es a través del método de las medianas , la cual consiste en unir los baricentros de los triángulos con los puntos medios de los lados generando un polígono al rededor de un vértice común de los triángulos adyacentes como se muestra en la figura 2, resultando el volumen de control Vc 1 delimitado por la frontera C 1. Modelo general de infiltración con volúmenes finitos
En primer término, la solución de la ecuación (1) expresado en la forma vectorial aplicando el MVF consiste en hacer la integración sobre un intervalo de tiempo $\Delta t$ y el volumen de control Vc , resultando la ecuación (2)
23 ++4 567 +,-./+3 23*
8
Modelamiento para flujo en estado estacionario, lo que nos interesa en este trabajo es
Solución para flujo estacionario
Para flujo estacionario t=∞, entonces (% () , por tal, la ecuación (1) se reduce a (2). Aplicando el Teorema de divergencia de Gauss , con la cual la integral volumétrica de la ecuación (3) es igual a !% y N es el vector la integral de superficie en el volumen de control Vc (4), siendo, " !% ' !" !' normal al superficie dirigido hacia a fuera en el volumen de control.
9
: $ $
;67 +,-./+3 23* 9 ;67 +,-./+3 23* ;>7 := +?
< @
Figura 3: Interpolación en el elemento triangular} Ahora bien, asumiendo valores del potencial h en los vértices del elemento, se hace necesario una función de interpolación para los puntos internos; por tal, nuestra función de interpolación para h sea lineal h=ax+bz+c
(5)
Calculamos a, b y c a partir de la interpolación de h 1, h 2 y h 3 en los vértices 1, 2 y 3 . siendo h 1=ax 1+bz 1+c para el vértice 1, h 2=ax 2+bz 2+c para el vértice 2 y h 3=ax 3 +bz 3+c para el vértice 3 ; ordenado en forma matricial y resolviendo por la matriz inversa (6) obtenemos los valores de a, b y c , ecuación (7), siendo det=x 1y 2+x 2 y 3+x 3 y 1-y 1x 2 -y 2x 3-y 3x 1
BACE F11GI GI HHKLG FGGK D zI 1JMJhJG HzJ zGGhI zG zI B P P h P dNOP h P P h C P z PJ z hI G P zG dNOPJz hI PI z G PJ z h Q D I J J I G J G dNOG J I G I !% I G !%J : $" !"!% $' !' ! % ! % 3+ $, $# &0 : $ $ " ' !" !R"ST'S7 !" !' !%!" !R"ST'S7 !% BU C0 !" !' !' : $" W$' X" $" X' W$' : X" X' 1 Y "ZS"I [ 0Y 'ZS'I [ 1\ "ZS"I ] 0\ 'ZS'I ] 1^ "[S"I ] 0^ '[S'I ]
Derivada del potencial h . De la ecuación (4)
velocidad de filtración en medios porosos es
(6)
. Según la Ley de Darcy la
entonces
.
De las ecuación (5) determinamos y por la Ley de Darcy ; , , donde . Integración. De la figura 3, siendo m,n,ñ puntos medios de los lados del triángulo, par 12, para 13, para 23. Desarrollando la integración nodo por nodo en el elemento triangular i formado por los puntos 1,2,3 , figura 3: Integración en el nodo 1.
Figura 4: Representación del vector F en las caras de la integración Según la figura 4, multiplicando F por N en los puntos r y t ; -f x(z m-z o )+f z(x m-x o ) es el flujo que atraviesa la cara mo , f x(z n -z o )-f z(x n-x o ) es el flujo que atraviesa la cara on . Ahora de la ecuación
(4),
;>7 := +? ;Y\ := +? ;Y_ := +? ;_\;:=\ :=+?+? X X 1 1 "\ Y 'Y \ Y \ 5Y := +? H8 `B$"I \JC$'1J 1Ia b ;Y := +? cGGHcGI8cGJ<
sub volumen de control , figura 4, resulta
, desarrollando la integral sobre cada cara del , remplazando los
términos f x, f z, x m, z m, x n y z n .
Con los valores de a y b en (8) obtenemos , integración en el nodo 1 con aportaciones de los nodos 2 y 3 del elemento i , y formando el sistema de ecuaciones para el elemento triangular 123 resulta la ecuación (9 ), en forma matricial la ecuación (10) de la forma Dh , siendo
gcGG 8+4HH `ijzI zJzI zJ ikPJ PIPJ PIa cGI 8+4 `ijzJ zGzI zJ ikPG PJPJ PIa H `ijzG zIzI zJ ikPI PGPJ PIa cGJ 8+4 H `ijzI zJzJ zG ikPJ PIPG PJa cIG 8+4 H `ijzJ zGzJ zG ikPG PJPG PJa cII 8+4 H `ijzG z8)(zJ − zG)+ik(PI − PG)(PG − PJ)a cIJ 8+4 H `i
cJG = 8h4 j(zI − zJ)(zG − zI)+ik(PJ − PI)(PI − PG)a cJI = 8h4H `ij(zJ − zG)(zG − zI)+ik(PG − PJ)(PI − PG)a cJJ = 8h4H `ij(zG − zI)(zG − zI)+ik(PI − PG)(PI − PG)a \ l5Y : ?mG = cGGℎG + cGIℎI + cGJℎJ Y l5^ : ?mI = cIGℎG + cIIℎI + cIJℎJ (n) ^ l5\ : ?mJ = cJGℎG + cJIℎI + cJJℎJ cFcGGIG ccGIII ccGJIJKFℎℎGGK cJG cJI cJJ ℎG cℎ = 2Vo ;67 (,-(./ℎ))3 −23* = 9
Solución. Retomando la ecuación (3),
(10)
, entonces
y la
solución es la ecuación (11); Donde, D matriz de propiedades (geometría y el coeficiente de permeabilidad ), h vector de alturas piezométricas (incógnita ), Q vector de condiciones de borde.
cℎ = o→ℎ = cLGo (HH)
De este modo llegamos a la solución de flujo de agua en medios porosos en su forma estacionaria. Para la programación utilizaremos la ecuación (11) como un algoritmo final . Cálculo de variables secundarios. Una vez resuelto el problema, el potencial h es conocido en todos los nodos del dominio de filtración . Entonces es posible determinar altura de presión
!%!" = −B 0 !%!' = C
(p = qrs(ℎ −)) 3' = X' = C$'
3" = X" = −B$" 0
, gradiente hidráulico ( . ), velocidad de filtración ) , definir líneas equipotenciales y líneas de corriente }, calcular caudal de filtración y volumen filtrado a través de la sección de la presa. Programación del MVF en MATLAB
Para este fin se ha utilizado el interfaz gráfica de usuario GUIDE de MATLAB. El programa elaborado se denomina CHUYMI.MVF v2009 (una terminología QUECHUA que en español significa infiltración ); programado bajo sistema de organización mostrada en el diagrama de flujo (figura 5), donde también se presenta el interfaz principal.
Figura 5: Diagrama de flujo de la programación estacionaria Validación del modelo MVF Condiciones de contorno
Estas condiciones resumidamente son de tipo Dirichlet (impone el valor del potencial h en la frontera del dominio), Neumann ( impone el valor del caudal que atraviesa la frontera del dominio), Cauchy (impone de forma combinada el potencial y el caudal), mixto (aplicado por partes, Dirichlet y Neumann ), [3].
Figura 6: Condiciones de borde en presas de tierra (a) y presa derivadora (b) Modelado de ejemplos genéricos Presa de tierra homogénea. Las figuras 7 y 8, representan la modelación de una presa de tierra -6 (k x=k m/s ), con base=52m, $altura=12m. y NAE=11m.$, discretizada en 543 nodos y 976 z=1x10
elementos.
Figura 7: Presa de tierra - Enmallado, líneas equipotenciales y Alturas piezométricas
Figura 8: Presa de tierra - Líneas de corriente, altura de presiones y velocidad de filtración Bocatoma con tablaestaca. Las figuras 9 y 10, es la modelación la fundación debajo de la -5 bocatoma, (k x=k m/s ) y NA=1.40m, discretizada en 1147 nodos y 3228 elementos. z=1x10
Figura 9: Bocatoma con tablaestaca - Enmallado, líneas equipotenciales y alturas piezométricas
Figura 10: Bocatoma con tablaestaca - Líneas de corriente, altura de presiones y velocidad de filtración Modelado de la presa Cuchoquesera Condición: Año 2009, Volumen en la presa 80MMC. En la fecha 20/04/09 la presa alcanza
llenar a su máxima capacidad operativa llegando a la cota 3741.70msnm y un volumen de 80MMC en el embalse. Las figuras 13, 14 y 15, es la modelación del dique principal de la PRESA CUCHOQUESERA , según condiciones mostradas en la figura 12 y el nivel de agua en el pozo de alivio PA-3=3702.48msnm como se muestra en la figura 11, la discretizado 1145 nodos y 2105 elementos.
Figura 11: Fotografía de la Presa Cuchoquesera y Pozo de alivio PA-3 desfogando
Figura 12: Modelo y Condiciones de borde- Dique Principal sección 0+420
Figura 13: Alturas Piezométricas y enmallado, Dique principal Sección 0+420: NAE=3741.70, PA- 3=3702.48 (20/04/09)
Figura 14:Alturas de presión y velocidad del flujo (líneas de corriente), Dique principal Sección 0+420: NAE=3741.70, PA-3=3702.48 (20/04/09)
Figura 15: Alturas de presión y velocidad del flujo (líneas de corriente), toma acercada
Conclusiones
La presente investigación cumpliendo con los objetivos trazados en su planteamiento, tiene un carácter científico y aplicativo en el campo de la ingeniería civil; primero por el modelamiento con una técnica avanzada y actual como lo es el MVF, segundo por la aplicación a un caso real como es la Presa Cuchoquesera como un elemento de validación del modelo propuesto. La ecuación de flujo potencial (1) tiene un rango de aplicabilidad sólo dentro de la Ley de Darcy (válido sólo para un flujo laminar ). MVF resulta formulación numéricas de altísima potencia para la solución de ecuaciones diferenciales, cuando el modelamiento es adecuado. Para el modelamiento del fenómeno en estado estacionario, el tiempo de modelación es infinito , por tanto, las líneas de corriente siempre alcanzará el nivel de interpolación de gradiente hidráulico , dependiendo de la velocidad de filtración , este último, directamente proporcional la coeficiente de permeabilidad . Para problemas de campo escalar resulta adecuada la formulación del MVF con mallas triangulares no estructuradas , y además el método de las medianas tiene mayores facilidades en cuanto a la sistematización de datos, aunque tiene mayor complejidad en la formulación, respecto a otros métodos como el diagrama de Voronoi . En cuanto a la validación del modelo propuesto con el modelamiento de la Presa Cuchoquesera , los resultados de modelación de la sección 0+420 del dique principal concuerda con los datos reales de los niveles piezométricos registrados históricamente, por tanto el modelo es válido en esta modelación. Más allá de estos resultados ningún modelo comercial ha resultado mejor que la nuestra en cuanto al modelamiento y convergencia de los resultados. Por tanto nuestra conclusión es nuevamente afirmar que el modelo es correcto y los resultados de algún modelamiento de un problema caen en la responsabilidad del modelador. •
•
•
•
•
•
Referencias
[1] LAMBE T. W. & WHITMAN R. V., Mecánica de suelos , LIMUSA, Noriega Editorese, Ingenieros de caminos, canales y puertos, Madrid - España. [2] TERZAGHI, K.; PECK, R. B. & MESRI, G. (1996), Soil mechanics in engineering practice }, A Wiley interscience Publication (Trird Edition)- New York. [3] TOLEDO MUNICIO, M. A. (1997), Tesis Doctoral: Presas de escollera sometidas a sobrevertido. Estudio del movimiento del agua a través de la escollera y de la estabilidad frente al deslizamiento en masa" , Universidad politécnica de Madrid - ETS. Ingenieros de caminos,
canales y puertos. [4] KUMAR, C. P. (2002), Introduction to groundwater modelling , Scientist E1 of National Institute of Hydrology, Roorkee – 247667, Uttaranchal -India [5] MALISKA, C. R. (1995), Transferencia de calor e mecánica dos fluidos computacional , Universidad federal de Santa Catarina. [6] Seepage Modeling with SEEP/W 2007 , GEO-SLOPE International Ltd. (Third Edition, March 2008). [7] ZIENKIEWICZ, O. C. & TAYLOR, R. L.(1994), El método de los elementos finitos - Volumen II (Mecánica de sólidos y fluidos, dinámica y no linealidad ), Mc Graw-Hill, cuarta edición. [8] PATANKAR, S. V. (1998), Numerical heat transfer and fluid flow , Hemisphere publishing Corporation - New York. [9] VERSTEEG, H. K. & MALALASEKERA, W. (1995), An itroducction to computational fluid dynamics - The finite volume method , Longman Scietific & Technical, New York.
[10]
GOIS, J. P.; ESTÁCIO, K. C. & OISCHI, C. M.(2007), Aplicaçao de volumes finitos na simulación numérica de contaminaçao em lençóis freáticos , Instituto de ciencias matemáticas e de de computaçao - Brasil. [11] EYMARD, R & GALLOUET, T. \& HERBIN, R (2006), Finite volume methods , Universite de Provence, Marseille.
