Embed Size (px)
Citation preview
Modelación numérica en microcomputadora de un acuífero confinado
Leopoldo Rodarte Ramón
Hidrolegro
En este trabajo se informa sobre los resultados de la primera etapa del estudio que emprendió el Instituto Mexicano de Tecnología del Agua para analizar la explotación eficiente de los recursos de agua subterráneos. Se presenta tanto el estado del arte de la aplicación de las microcomputadoras al estudio de las aguas subterráneas como el diseño de u n programa de cómputo para la modelación matemática en microcompu- tadora de un acuífero confinado.
México es un país con una distribución de los recursos hidráulicos muy irregular, por lo que el estudio del agua para su aprovechamiento racio- nal adquiere especial importancia.
De acuerdo con la distribución espacial de la lluvia y la temperatura, el 31% del territorio nacio- nal se clasifica como desértico y árido, el 36% es semiárido y el 33% restante es semihúmedo y hú- medo. Un 42% del territorio nacional, ubicado En nuestro país el aprovechamiento de las aguas principalmente en el norte del país, tiene precipi- subterráneas mediante la perforación de pozos taciones medias anuales inferiores a los 500 milí- cobró importancia en los años 20 del presente metros. En estas regiones áridas y semiáridas los siglo, a fin de suministrar agua potable a la ciudad escurrimientos superficiales son escasos y en de México. Debido al desarrollo económico y al muchas ocasiones torrenciales, lo que dificulta su crecimiento de la población, la demanda de agua aprovechamiento económico y hace necesario re- subterránea tanto para el consumo urbano- currir a la explotación de las aguas subterráneas industrial como para la producción agrícola, ha como Única o principal fuente de abastecimiento. crecido con mucha rapidez; esto ha creado la ne-
Las peculiaridades en el estudio de las aguas cesidad de llevar a cabo estudios más profundos subterráneas han hecho que su explotación en el acerca de las características de las formaciones país sea menos importante que la de las aguas acuíferas. superficiales. Muchas veces se tiene un conoci- A partir de los sesenta se introdujeron en Méxi- miento preliminar sobre los acuíferos y por des- co las primeras computadoras electrónicas en las conocer su potencial, se ha ocasionado su so- tareas de la ingeniería civil y, entre ellas, en el breexplotación o degradación, o bien su subex- estudio de las aguas subterráneas. plotación. La modelación numérica de acuíferos es una
Con el advenimiento de la era electrónica, las herramienta que puede ser de gran utilidad tanto computadoras han adquirido gran importancia en las etapas previas del conocimiento de un para resolver problemas prácticos en ingeniería. acuífero como en las finales, en las cuales se pre- En especial, las microcomputadoras han permiti- tende predecir el comportamiento de una estruc- do calcular problemas complejos en el nivel de tura sujeta a distintas políticas de explotación.
oficina sin necesidad de recurrir a computadoras centrales, cuyo empleo impide la interacción en- tre el especialista y este poderoso equipo de cál- culo, así como la utilización eficiente del personal técnico disponible en cada dependencia.
Estado del arte
Diferentes investigadores (Carlos Cruickshank y Rubén Chávez Guillén, 1968; Darío de Hoyos, 1968 y Leopoldo Rodarte, 1969) trabajaron con los programas elaborados bajo el patrocinio de la Dirección de Aguas Subterráneas y diseñados pa- ra su procesamiento en una computadora de se- gunda generación CDC-3100, con una reducida capacidad disponible de memoria (de 32k palabras).
A fines de esa década aparecieron también en la literatura técnica norteamericana los primeros modelos numéricos de importancia. Una descrip- ción amplia sobre las técnicas de modelación pa- ra evaluar los recursos hidráulicos se encuentra en la obra de Prickett (1979). A juicio de este autor, destacan los modelos numéricos en los que se emplea el método de las diferencias finitas, ela- borado por Pinder y Bredehoeft (1968), y Prickett y Lonnquist (1971), así como el de los elementos finitos aplicados en los trabajos de Witherspoon et. al., (1968). En todos estos casos, los progra- mas de modelación se elaboraron para computa- doras de tamaño mediano y grande.
Las prime ras microcomputadoras se introduje- ron también a fines de los sesenta y tienen gran- des ventajas en relación con las minicomputado- ras y las computadoras medianas y grandes. Destaca el hecho de que pueden ser instaladas en cualquier lugar, funcionan de manera satisfacto- ria en un amplio rango de temperaturas y se evita el uso de acondicionadores de aire.
En 1973, Prickett y Lonnquist publicaron una versión de su programa procesable en computa- doras de tamaño pequeño que funciona en las microcomputadoras actuales.
Hoy en día existen variados programas para el análisis de problemas de aguas subterráneas. Gran parte de este software puede ser adquirido a través del International Ground Water Modeling Center del Holcomb Research Institute, Butler University, Indianapolis, EUA, a precios diferen- tes, según la complejidad del modelo. Entre los programas disponibles para microcomputadoras IBM-PC o compatibles se pueden ci tar los siguientes:
A modular Three- Dimensional Finite -Differen- ce Ground Water Flow Model, Computer Model of Two-Dimensional Solute. Transport and Dispertion i n Ground Water (MOC). Numerical Method of Pumping Test Analysis Using Microcomputers. El South Australia Department of Mines and
Energy ofrece un paquete de 25 programas para
el análisis en microcomputadoras de problemas de aguas subterráneas. El International Ground Water Modeling Center, como resultado del curso Practical Analysis of Well Hydraulic and Acuifer Pol- lution, ofrece también un paquete de 35 progra- mas con soluciones analíticas y numéricas. Por su parte, en Hidrolegro existen tres programas sobre este mismo tema. En dos de ellos aparece el modelo de Prickett y Longquist. El tercero es el que acompaña al libro Modeling Ground Water Flow and Pollution de Bear y Verruijt (1987). Estos programas permiten la modelación numérica de diferentes variantes de problemas de acuíferos tanto estacionarios como transitorios.
Por otra parte, en la actualidad se cuenta con una gran variedad de modelos digitales y progra- mas para la resolución de problemas prácticos de geohidrología.
Modelo digital para un acuífero confinado
La parte medular de este trabajo la constituye el diseño o adaptación de un programa para la mo- delación numérica de un acuífero confinado para una computadora IBM-XT o compatible, en el que es necesario tomar en cuenta varios aspectos.
Características de los acuíferos reales
La ecuación diferencial que representa a un acuí- fero confinado, anisótropo, no homogéneo, en el que actúa un sistema de pozos de extracción o de recarga, es:
donde: Xi= X xj = Y Ti¡= Componentes del tensor de
trasmisibilidad en un sistema no principal de referencia (L2/T). Gastos de extracción o recarga, por unidad de área producidos por un sistema de pozos (L/T).
h = Carga hidráulica en el acuífero
S = Coeficiente de almacenaje del
Q =
( L ) .
acuífero (Lo).
La ecuación anterior se puede simplificar me- diante la elección de un sistema principal de refe- rencia en el cual los términos Txy y Tyx sean igual
a cero. En este caso, la ecuación anterior se escribirá:
Esta ecuación ha sido resuelta en forma analíti- ca para algunos casos de condiciones iniciales y de frontera simples. Si se considera que el acuífe- ro es isótropo, homogéneo, de espesor constante, extensión radial infinita y que en él actúa un pozo del cual se extrae un gasto constante Q, la solu- ción analítica a este problema está dada por la expresión (Theis, 1935):
sólo tiene solución analítica para condiciones ini- ciales y de frontera relativamente simples, por lo que para resolver problemas prácticos en que se requiera conocer el comportamiento a largo plazo de un acuífero ante la acción de un sistema de pozos, casi siempre se recurre a la modelación matemática en computadoras digitales. Estas son las que han logrado una mayor difusión en nues- tro país.
AI elaborar un modelo matemático, las deriva- das parciales con respecto a las coordenadas y al tiempo se representan mediante operadores aproximados que transforman la ecuación dife- rencial en un sistema de ecuaciones algebraicas, sistema, que una vez resuelto, permite conocer la magnitud de la carga hidráulica en puntos discre- tos de una malla de integración para diferentes
De acuerdo con ésta, los abatimientos que pro- Los métodos numéricos aplicados con mayor duce el pozo de bombeo son directamente pro- éxito al análisis de acuíferos confinados pueden porcionales al gasto y crecen con el transcurso clasificarse en dos grandes grupos: los de dife- del tiempo; alcanzan su mayor magnitud en las rencias finitas y los variacionales. En ambos la cercanías del pozo y decrecen rápidamente a me- calidad de una aproximación depende en alto dida que el punto analizado se encuentra a mayor grado del error de truncamiento de la misma, de distancia de éste. su convergencia y su estabilidad.
A fin de determinar las políticas de explotación Error de truncamiento. Supóngase que se quie- por seguir en la resolución de problemas de acuí- re saber el valor de una función f ( x ) en un punto x feros, éstos se presentan en la naturaleza como cuando' se conoce el valor que toma en x . anisótropos, no homogéneos, de forma irregular Para ello se recurre a su desarrollo en serie de en el plano horizontal, y con condiciones de fron- Taylor. tera complejas, además de que esxplotados
mediante un sistema de pozos, con recargas y pérdidas producidas por riego, lluvias, evapora- ción, transpiración, etc. Esto origina la necesidad de estudiarlos mediante modelos, ya que la ob- tención de una solución analítica exacta no sería posible para condiciones tan complejas. En nues- tro país, los modelos más eficientes y populares son los digitales, ya que se carece de la posibili- dad de utilizar tipos distintos (analógicos, eléctri- cos, Hele-Shaw, etc.).
Conviene recalcar que un buen modelo digital para un acuífero confinado debe permitir un aná- lisis general de la ecuación (1), tomando en cuen- ta las peculiaridades que produce su explotación mediante un sistema de pozos, que en forma sim- plificada están dadas por la solución de Theis (ecuación 3) y el principio de superposición de causas y efectos.
Aspectos matemáticos del diseño
La ecuación (1) es diferencial lineal, en derivadas parciales, de tipo parabólico. Ya se mencionó que
valores valores de tiempo. tiempo.
valores de tiempo.
Por ejemplo, para la primera derivada:
El error de truncamiento es aquél que se come- te al cortar un desarrollo en serie. Normalmente cuando se habla de este tipo de error siempre se establece el orden de la aproximación (se dice que un término A es de orden (Ax)", lo cual se es- cribe 0 ( ( A x ) " ) si existe una constante positiva C, independiente de A x , tal que A<C Ob- sérvese que al definir una aproximación, cuanto mayor sea ésta, será más precisa. Así, por ejemplo, el operador de la primera derivada por el
método de diferencias finitas presenta las si- guientes variantes con su correspondiente error de truncamiento.
diferencias hacia adelante:
diferencias hacia atrás:
diferencias centrales:
La aproximación más exacta es la última por ser la que tiene un mayor error de truncamiento.
Convergencia y estabilidad. Para que un méto- do numérico sea aplicable conviene que la solu- ción obtenida con él sea muy parecida a la solu- ción analítica del problema. Por ejemplo, sea h ( x , y , t ) , la solución analítica exacta de un pro- blema y g ( i A x , j A y , k A t ) , la solución numérica del mismo obtenida por un método cualquiera. La diferencia entre ambas, o sea el error que se co- mete al aplicar un método numérico en particular, es Ei , j , k = h ( x , y ,
Convergencia. Se dice que un método numéri- co es convergente si cuando A x , A y y A t tien- den a cero, el error Ei , j ,
Error de truncamiento y convergencia son con- ceptos íntimamente ligados; cuanto mayor sea el primero, mejor será la convergencia de una aproximación.
Un aspecto muy importante del comportamien- to de una solución numérica es aquel que ocurre cuando el tiempo tiende a infinito mientras Ax, A y y At permanecen fijos.
Estabilidad. Se dice que un método numérico es estable si cuando el tiempo tiende a infinito, con A x , A y y A t fijos, los errores contenidos en su aplicación permanecen constantes o decrecen.
A fin de normar criterios de cómo se aplican estos conceptos a la ecuación diferencial (1) con- viene comparar dos métodos de diferencias fin¡- tas; el implícito y el de Crank-Nicolson. Una expli- cación más amplia sobre los conceptos de convergencia y estabilidad y acerca de los opera- dores que definen estos métodos puede encon- trarse en Remson et a / (1971).
Por lo anterior, el método de Crank-Nicolson permite obtener un modelo de mejor calidad que el implícito ai aplicarse a la ecuación 1, ya que es más convergente en el tiempo.
Además de los factores arriba señalados, en la eficiencia de un modelo digital influye el método que se emplee para dar solución al sistema de ecuaciones algebraicas. Estos métodos se divi- den en dos grandes grupos: los directos y los iterativos.
Los primeros son aquellos que mediante un nu- mero finito de operaciones permiten obtener la solución exacta de un sistema de ecuaciones al- gebraicas, como el método de suma y resta, el de sustitución, el de Gauss, etcetera.
Los iterativos son aquellos en los que se propo- ne inicialmente un valor numérico para las incóg- nitas y después, mediante iteraciones, se obtiene un valor aproximado satisfaciendo un error máxi- mo determinado con anterioridad.
El tiempo de cálculo requerido para resolver un sistema de ecuaciones algebraicas depende prin- cipalmente del número de multiplicaciones y divi- siones que es necesario efectuar para llegar al resultado final. Así por ejemplo, para resolver un sistema de n ecuaciones con n incógnitas es ne- cesario, según el método de Gauss, efectuar n3/3 multiplicaciones o divisiones, mientras que el al- goritmo que permite resolver matrices tridiagona- les (en las que como máximo hay tres coeficien- tes diferentes de cero en cada renglón), efectúa únicamente (5n-4) multiplicaciones o divisiones. Como es evidente el tiempo requerido para resol- ver un problema es de gran importancia, ya que si es muy grande puede bloquear el uso de la com-
putadora. Para evitar esto, en función de la veloci- dad de cálculo de la computadora, normalmente se restringe el número máximo de puntos (incóg- nitas), que pueden plantearse para un problema determinado.
Parametros de microcomputadora que influyen en la eficiencia de un modelo
Hasta hace poco, la mayor parte de los programas de modelación de trabajos de hidrología y geohi- drología se elaboraba para ser procesada en com- putadoras medianas o grandes que por lo común efectúan operaciones matemáticas, a alta veloci- dad; cuentan además con un procesador central, de por lo menos un megabyte, y una memoria ex- terna en disco duro Prácticamente ilimitada. En el diseño de programas para la modelación de acuí- feros en estas computadoras la restricción mayor la impone el hecho de que, aun cuando cuentan con un procesador central (CPU) de gran capaci- dad, al captar simultáneamente varios programas para su procesamiento (multiprogramación), si uno de ellos utiliza una gran parte de la memoria, se reducen las posibilidades de la multiprograma- ción. En estos casos, en la modelación de acuífe- ros se recurre normalmente al almacenamiento parcial de la información, tanto de datos como de variables, en unidades de disco duro, que tienen un tiempo de acceso mucho menor que las cintas magnéticas cuando se trata de recuperar informa- ción. Es claro que este almacenamiento externo incrementa el tiempo de procesamiento en rela- ción con la velocidad que se lograría si toda la información estuviera almacenada directamente en el procesador central.
Las matrices que representan al sistema de ecuaciones algebraicas en que se transforma la ecuación diferencial de un acuífero confinado presentan dos ventajas para su análisis: e Son simétricas con respecto a la diagonal
principal. e Si se elige de manera adecuada la numeración
de los nudos de la malla de integración, la ma- yor parte de los coeficientes diferentes de cero se agrupa en torno de la diagonal principal. Por tanto, en computadoras medianas y gran-
des puede utilizarse algún método para resolver sistemas de ecuaciones que aproveche la simetría de la matriz para obtener una solución más rápi- da. En cuanto a la rapidez de resolución del siste- ma de ecuaciones algebraicas para el análisis de
un acuífero confinado, merece atención especial el algoritmo de matrices tridiagonales.
AI resolver la ecuación diferencial (1) por el método implícito de dirección alternante (ADI), la matriz de coeficientes que se obtiene del sistema de ecuaciones algebraicas es tridiagonal. El algo- ritmo de matrices tridiagonales es muy eficiente para resolver este tipo de matrices. El método ADI fue utilizado por Prickett y Lonnquist (1971), para diseñar su modelo de acuíferos confinados. Una versión de este programa, traducido al español, especialmente adaptado a microcomputadoras, forma parte del informe entregado al Instituto Me- xicano de Tecnología del Agua. La segunda ven- taja radica en que la mayor parte de los coeficien- tes de la matriz del sistema de ecuaciones son ceros, en tanto que los diferentes de cero se agru- pan en torno de la diagonal principal.
Con objeto de ahorrar almacenamiento, estas matrices pueden ser transformadas en matrices rectangulares con el mismo número de renglones pero con un número de columnas tal que sólo se analicen los coeficientes diferentes de cero. Una mayor información sobre este algoritmo se puede encontrar en Thurnau (1963).
Otra posibilidad para resolver sistemas de ecuaciones grandes es el análisis de matrices me- diante particiones, lo que también permite reducir el almacenamiento en el procesador central (CPU). Las particiones que no se trabajan se pue- den guardar en el almacenamiento externo; para nuestro caso, en el disco duro. Un ejemplo de este tipo de programas puede verse en Rodarte (1974). Conviene aclarar que en computadoras medianas y grandes, en general, la velocidad de cálculo del procesador central no es una limitante para el análisis de acuíferos confinados, ya que la mayor parte de los problemas prácticos de la modelación de éstos puede resolverse en unos cuantos minutos de CPU.
Las microcomputadoras presentan diferencias fundamentales con respecto a las computadoras arriba mencionadas. Hasta la generación de mi- crocomputadoras IBM-XT sólo es posible cargar y procesar un solo programa, es decir, no aceptan multiprogramación. Por ello, en el programa que se diseñe se puede utilizar toda la capacidad del procesador central.
La configuración que se tomó como base para el diseño o adaptación del programa de modela- ción de acuíferos tiene las siguientes característi- cas: e Capacidad de memoria del procesador central:
512K.
Capacidad del disco duro: 10 MB. CPU INTEL 8086 con coprocesador matemático 8087. Si bien la capacidad de memoria en apariencia
no es una limitante para el diseño o adaptación del programa de análisis de acuíferos confinados, la velocidad del procesador central sí establece un limite máximo del tamaño (número de nudos o incógnitas), del acuífero que se estudie.
Método del elemento finito en la modelación de acuíferos
Los problemas de flujo transitorio en medios po- rosos sólo se pueden estudiar mediante métodos analíticos en aquellos casos en que las caracterís- ticas del medio y las condiciones de frontera son suficientemente simples. En los problemas prácti- cos de geohidrología es frecuente encontrar con- diciones de frontera complejas, tales como un medio homogéneo y anisótropo, forma irregular de la región de estudio, etcétera.
El método del elemento finito es un método numérico con grandes ventajas para el análisis de problemas geohidrológicos, entre las que sobre- salen las siguientes:
No requiere el desarrollo de expresiones espe- ciales para definir las condiciones de frontera, como ocurre, por ejemplo, cuando se utiliza el método de diferencias finitas. Las dimensiones de los elementos pueden ser fijadas en forma arbitraria. En aquellas zonas del acuífero donde la función varía rápidamen- te, por ejemplo en la vecindad de un pozo, pue- den elegirse elementos pequeños, en tanto que en aquéllas donde lo haga con lentitud, se pue- de optar por elementos grandes. AI utilizarse elementos de forma triangular se puede obtener una mayor coincidencia con las fronteras geométricas reales. El método del elemento finito se puede plan-
tear de diferentes maneras; en este trabajo se uti- lizará una formulación variacional del problema matemático original.
Transformación de la ecuación diferencial
La ecuación diferencial (1) para el caso de un acuífero confinado bidimensional se escribe en la forma:
Para simplificar el planteamiento variacional del problema se integrará la ecuación (5) de t=o a t=t , aplicando para ello el operador:
luego,
de manera semejante:
Donde h' es el valor de la carga hidráulica al final del intervalo y ho, el valor al principio del mismo. AI substituir en (5) las expresiones (7) a (10) se obtiene finalmente la ecuación:
donde h y Q deben ser considerados como valo- res promedio en el intervalo de tiempo. Para sim- plificar la notación, se ha omitido en estas Últimas variables el símbolo que indica valor promedio. Si además se considera que el valor promedio de la carga durante el intervalo de tiempo t puede ser expresado en función de los valores que toma al principio y final del mismo mediante la relación:
i es una constante de interpolación O < i < 1 AI substituir esta última expresión en (11), pero
escritaen la forma: h'-ho= se obtiene final- mente: I
Esta ecuación puede ser empleada para el aná- lisis numérico de un acuífero confinado anisótro- po y heterogéneo; es similar a la de un acuífero semiconfinado estacionario (Rodarte, 1978), y servirá de base para el análisis numérico de un
acuífero confinado mediante el método del ele- mento finito.
Aplicación del método del elemento finito
Considérese la funcional con:
donde:
Se puede demostrar (Kantorovich y Krilov, 1949), que la funcional (13) es equivalente a la ecuación diferencial (12) para aquellos casos en que las fronteras del acuífero sean impermeables o de carga constante.
Para el análisis, la región R definida por el acuí- fero se divide en elementos triangulares tal como se muestra en la ecuación (5), (véase ilustración 1).
Supóngase, además, que los nodos de la re- gión R de integración se numeran consecutiva- mente de 1 a n’. Si los elementos son lo suficiente- mente pequeños, la carga hidráulica en cada uno de ellos puede expresarse mediante la función lineal:
Para un triángulo cualquiera con nodos j , k, I, tomados en un orden previamente definido (véase ilustración 1), la carga hidráulica en cada nodo se calcula mediante las expresiones:
AI resolver este sistema de ecuaciones, con respecto a las incógnitas a, b y c es posible eva- luar h en cada elemento substituyendo (15) en (14):
donde:
mediante las expresiones (17) a (19), la ecuación (16) puede ser escrita en la forma siguiente:
y a partir de (20), se pueden definir las ecuaciones:
y de estas últimas expresiones:
al substituir (20) y (21) en (13), para cada elemen- to se define la funcional:
si además se define
U, V y W son matrices de 3 x 3 elementos y el superíndice T significa en (26) la transpuesta de la matriz. Para definir las matrices conviene desa- rrollar el primer renglón:
AI diferenciar con respecto a hj y agrupar los términos comunes (23) se transforma en:
AI substituir estas expresiones en (23) se obtiene:
La ecuación anterior puede ser escrita en un sistema centroidal de referencia (el origen del sis- tema de referencia coincide con el centroide de cada elemento), con el f in de simplificar su plan- teamiento. Para ello se define: De manera semejante se puede calcular
a fin de determinar:
donde:
La ecuación de minimización de la funcional se expresa en la forma:
donde la suma se efectúa sobre todos los elemen- tos. Con base en (30), la ecuación (29) se escribe de la siguiente manera:
donde las sumas se efectúan sobre todos los ele- mentos y sobre los nudos de la malla de integra- ción. La ecuación anterior puede ser descrita con la siguiente notación:
El término F que corresponde al sistema de ecuaciones algebraicas (32) representa el efecto que produce en las cargas hidráulicas la acción de un sistema de pozos de bombeo o de recarga. En la ecuación (32) puede observarse que la va- riación de la carga hidráulica con el tiempo puede ser producida tanto por la acción de un sistema de pozos (F) como por la de condiciones iniciales y de frontera (término ho).
Es conveniente recalcar que los intervalos de tiempo t que se utilizan en la resolución de la ecuación anterior, no pueden ser elegidos arbitra- riamente, ya que es fundamental respetar las hi- pótesis planteadas con el fin de que sea posible emplear en lugar de la ecuación (5) la (12) que es mucho más simple de resolver y que hace más sencilla la aplicación del método del elemento fi- nito. Una solución más rigurosa de la ecuación diferencial (5) por medio de este método se en- cuentra en Rodarte (1974, 1978).
Uso Uso de elementos cuadrangulares y matrices donde P y R son matrices cuadradas de n x n elementos (n es el número de nudos de la malla) y F es un vector de n elementos.
Programación en Basic de la ecuación diferencial
Como se ha visto, la ecuación diferencial para el análisis de un acuífero confinado transitorio bidi- mensional es (12) y, al aplicarle el método del elemento finito se transforma en:
rectangulares
AI obtener la ecuación (32) por el método del ele- mento finito se ha dividido la región R definida por el acuífero en elementos triangulares tal como se muestran en la ilustración 1.
La utilización de cuadriláteros presenta venta- jas en el análisis numérico, ya que aquéllos per- miten reducir el número de incógnitas que se plantean sin que por ello se pierda precisión en los cálculos. Para emplear cuadriláteros por el método del elemento finito se procede de la si- guiente manera:
donde P y R son matrices cuadradas de n x n elementos (n es el número de nudos de la malla) y F es un vector de n elementos; P representa la contribución de los dos primeros términos de la ecuación (12) y es también directamente propor- cional a los valores de la transmisibilidad Txx y
R representa la contribución del tercer término de (12) y es también directamente proporcional al valor del coeficiente de almacenaje. Para el caso en que S = O en todos los puntos del acuífero, la ecuación (32) se simplifica transformándose en:
Tyy.
ecuación que es válida para el análisis de un acuí- fero estacionario.
Se divide el cuadrilátero en cuatro triángulos (I, II, Ill, IV) que tienen como vértice común el punto E, cuyas coordenadas se colocan por comodidad en el centroide de la figura. La contribución del cuadrilátero se analiza sumando las de cada uno de los triángulos para después eliminar (despe- jando) el nudo ficticio E. De esta manera la preci- sión obtenida en los cálculos es mayor o igual
que la que se obtendría al dividir el cuadrilátero en cuatro o en dos triángulos, reduciendo al mis- mo tiempo (para el caso de 4 triángulos) el núme- ro final de incógnitas que es necesario plantear para resolver el problema.
En el programa planteado la solución se obtie- ne a partir de elementos cuadrangulares. Cuando en lugar de un cuadrilátero sea necesario utilizar un triángulo, basta con repetir uno de los nudos para que aquél se convierta en triángulo.
En la ecuación (32) P y R son matrices cuadra- das de n x n elementos donde n es el número de nudos de la malla. En estas matrices, la mayor parte de los elementos son ceros y si la numera- ción de la malla se elige de manera adecuada, los elementos diferentes de cero quedan agrupa- das cerca de la diagonal principal. Por ejemplo, para el elemento mostrado en la ilustración 2 los elementos que constituyen a la matriz P son los siguientes:
De 121 coeficientes, para el problema plantea- do, sólo 41 son diferentes de cero lo que repre- sentaría un gran desperdicio de memoria si para resolver el sistema de ecuaciones se utilizara un método que sólo pueda resolver matrices cuadra- das.
Obsérvese, además, que en cada renglón sólo se necesita conservar cinco posiciones para al- macenar los coeficientes, y que el semiancho de la banda ( B ) (distancia del coeficiente de la dia- gonal principal al coeficiente más alejado diferen- te de cero) es de 2, por lo tanto el número máximo de posiciones de almacenamiento por renglón es de 2B + 1 = 5, para este caso. Mucho más lógico sería almacenar la matriz anterior de la manera siguiente:
Esta matriz rectangular de 11 X 5 sólo almace- na los coeficientes que se encuentran a una dis- tancia menor o igual a B que, en nuestro caso, son todos los coeficientes diferentes de cero. Los coeficientes de la diagonal principal se encuen- tran en la posición B + 1, a la izquierda y derecha los coeficientes diferentes de cero, conservando su posición original con respecto a la diagonal principal.
Si se utiliza un método apropiado para resolver este tipo de matrices (en este ejemplo se emplea el del gradiente conjugado) el ahorro en tiempo de cálculo es significativo. Una explicación más amplia sobre este método se encuentra en Gam- bolati, G. et al, (1984).
Por tanto, en síntesis, la solución numérica que se propone para el análisis bidimensional de flujo transitorio en acuíferos confinados se plantea mediante el método del elemento finito, usando elementos cuadrangulares o triangulares donde las matrices de coeficientes que se utilizan son rectangulares y el sistema de ecuaciones resul- tante se resuelve mediante el método del gradien- te conjugado. Por último, conviene mencionar que la capacidad de almacenamiento de la computadora para la cual se diseñó el programa permite como máximo utilizar 150 nudos y un se- miancho de banda de 11 (23 coeficientes diferen- tes de cero por renglón).
Un ejemplo de aplicación del programa
U n programa de la complejidad del elabora- do por Hidrolegro, para el análisis de un acuífero confinado, implica comparar la solución numéri- ca con una analítica (exacta) para un mismo pro- blema, a fin de comprobar que el programa fun- ciona correctamente y con la aproximación suficiente para ser aplicado a la resolución de ta- reas prácticas. A este proceso se le llama calibra- ción del modelo. En el caso que nos ocupa, para la calibración se hará uso de un problema tomado de Bear J. y Verruijt (1987). En el acuífero mostrado en la ilustración 3 (1000 m de longitud y 100 m de espesor) actúa una infil- tración neta uniforme de 0.001 m/día. El estrato
que separa al acuífero estudiado del acuífero infe- rior se considera impermeable, al igual que su frontera izquierda ( X = O). Inicialmente, la carga hidráulica en el acuífero es igual a cero y para el tiempo t= 0 se inicia la infiltración I. La transmisi- bilidad del acuífero se considera como T = 1 0 0 0 m2/día y el coeficiente de almacenaje S = 0.40; para resolver el problema se usa una malla de integración con un paso de 100 m. Los datos con- centrados de la malla se anexan a continuación. Los cálculos se efectuaran para los 13 diferentes intervalos de tiempo indicados. Los resultados obtenidos mediante el programa que definen el valor de la carga hidráulica ( F en metros) para diferentes valores del tiempo se enlistan a conti- n nuación. n.
La solución numérica obtenida por Hidrolegro es prácticamente igual a la del programa de Bear y Verruijt (1987). De acuerdo con estos autores, la solución numérica a este problema se aproxima a la solución analítica con errores menores del 1%. En el libro antes citado se comparan gráficamente las soluciones analíticas y numéricas.
Bibliografía
Bear, Jacob y Arnold Verruijt. Modeling Ground Water Flow and Pollution. Reidel Publishing Company, Dordrecht, 1987.
Bittinger, M.W., Duke H.R. y Longenbaugh, R.A. Math- ematical Simulations for Better Aquifer Manage- ment, International Association of Science and Hydrology Symposium Artificial Recharge Manage of Aquifers, Haifa, pp. 509-519, 1967.
Métodos de ordenador para evaluación de recursos hi- dráulicos subterráneos, Boletin núm. 41, Ministerio de Obras Públicas, Dirección General de Obras Hi- dráulicas, España. 1976.
Butler, S.S. Engineering Hydrology, Prentice-Hall, En- glewood Cliffs, New Jersey, 1957.
Cruickshank Carlos y Rubén Chávez Guillén, “Modelo DAS para el estudio de aguas subterráneas“, Inge- niería Hidráulica en México, 1968.
Gambolati, G. y A. Perdon. “The Conjugate Gradients in Subsurface Flow and Subsidence Modeling”, Fun- damentals of Transport Phenomena in Porus Media, J. Bear y M.Y. Corapciuglu (eds.) pp. 953-984, Martí- nez Nijnoff, Dordrecht, 1984.
Hoyos, Darío de. “Estudios de simulación de funciona- miento de acuíferos mediante modelos matemáti- cos”. Informe presentado por IPESA constructores S.A. a la Dirección de Aguas subterráneas, SRH. 1968.
Javadel I., y P.A. Witherspoon. Analysis of Transient Fluid Flow in Multilayered System, Contribution R4, Water Resources Center University of California, EUA.
Kantorovich, L.B. y V.I. Krilov. Pniblizhennye Metody By S Shego Analiza, Giztexizdot, 1946.
Karplus, W.J. Analog Simulation, McGraw-Hill, Nueva York, 1958.
Lin, C.L. “Digital simulation of an Outwash Aquifer”, Ground Water, 11(2), 38-43, 1973.
Peaceman, P.W. y Rachford, H.H., Jr. The Numerical Solution of Parabolic and Elliptic Differential Equa- tions. Journal o f Soc. Ind. Appl. Math., 3(1), 28-41, 1955.
Pinder G.F. y Bredehoeff J.D. “Application of the Digi- tal Computer for Aquifer Evaluation”, Water Resour- ces Research, 4(5), 1069-1093, 1968.
Prickett, Thomas A. “Modeling Techniques for Ground- water Evaluation, Advances in Hydroscience”, 10, 1- 143, 1979.
Prickett, T.A. y Lonnquist, C.G. Selected Digital Com- puter Techniques for Groundwater Resources Eva- luation, Boletín 55, Illinois State Water Surv., Urbana, 1971.
Prickett, T.A. y Lonnquist C.G. “Aquifer Simulation Mo- del for Use on Disk Supported Small Computer System”, Circular 114 Illinois State Water Surv., Ur- bana.
Remson, I., G.H. Horneberg y F.J. Molz. Numerical Meth- ods in Subsurface Hidrology with Introduction to the Finite Element Method, Wiley-lnterscience, Nue- va York, 1971.
Rodarte, Leopoldo. Tesis doctoral, Instituto de Ingenie- ros Civiles V.V. Kuibuishey, Moscú, URSS, 1974.
Rodarte, Leopoldo. Application o f Finite Elements. Method to the Numerical Analysis of a Leaky Aqui- fer, 1.3-1.18, 1976: Finite Element in Water Resour- ces, Pentech Press, Londres, 1978.
Rodarte R. Leopoldo. “Un modelo matemático para el estudio de acuíferos semiconfinados múltiples”. In- forme presentado por el Instituto de Geofísica a la Dirección de Aguas Subterráneas, SRH, 1969.
Theis, C.V. “The Relationship Between the Lowering of the Piezometric Surface and the Rate and Duration of the Discharge of a Well Using Groundwater Stora- ge”, Transl. Am. Geophys, Union, 16, pp. 519-524, 1935.
Thurnau, D.F. “Algorithm 195, BANDSOLVE”, Commu- nications of the ACM, 6, núm. 8, p. 441, 1963.
Tyson, H.N. Jr. y Weber, E.M. “Groundwater Manage- ment for the Nation’s Future, Computer Simulation of Grounwater Basins”, Journal of Hydraulics, Div. Proc. Amer. Soc. Civil Eng. 90 (HY 4), 59-77, 1964.
Witherspoon, P.A., Jayvel I. y Newman S.P. “Use of Finite Elements Method in Solving Transient Flow Problems in Aquifer Systems”, Inst. Ass. Asi. Hydrol. Simp, Publ 81. 687-698, 1968.
