Re is a In e nacional de Mé odos Numé icos pa a Cálculo
y
Diseño en Ingenie ía. Vol.
9,1, 15-34( 1993)
OPTIMIZACION DE LA SOLUCION NUMERICA
POR ELEMENTOS FINITOS
DE SONDEOS GEBELECTIRICOS
O.A. SOTO
Y
C. MBLANO
Dp o. de Ingenie ía Ci il,
Uni e sidad de los Andes A.A.4976,
Bogo á, Colombia.
RESUMEN
En el a ículo se p esen a una nue a me odología de solución pa a esol e p oblemas
de lujo y po encial modelados con elemen os ini os. El sis ema de ecuaciones esul an es
de dicha solución conocido como CGEIS (G adien es Conjugados con P econdicionamien o
y
Escalonamien o, u ilizando la Descomposición Incomple a de Cholesky). Es e mé odo puede
se u ilizado sólo pa a el caso de ma ices simé icas y posi i amen e de inidas, las cuales
son encon adas en p oblemas de po encial ales como sondeos geoeléc icos, lujo de aguas
sub e áneas, hid aúlica, con aminación, pe óleos, e c.. El mé odo ue implemen ado con
el in de simula sondeos geoeléc icos y lujo de aguas sub e áneas. Su mon aje se ealizó
en un compu ado VAX-2, y se u ilizó el FORTRAN 77, como lenguaje de p og amación.
Adicionalmen e el mé odo op imiza los ecu sos compu acionales ales como iempo de ejecución
y memo ia p incipal.
SUMMARY
A new solu ion me hod is p esen ed o sol e po encial and low p oblems, when a ini e
elemen ep esen a ion is used. The esul ing equa ions can be a anged in a ma ix o m and
apply a new solu ion me hod called CGEIS (Conjuga ed G adien s Escaled and P econdi ioned,
using he Incomple Cholesky Descomposi ion). This me hod can be used o syme ic and
posi i ely de ined ma ices, which a e ound in many po en ials p oblems such as in ese oi
enginee ing, geoelec ics, hid aulics, oil, e c.. The me hod was implemen ed o sol e complex
ma ices ob ained in s anda and con inuous geolec ics soundings, and in g oundwa e low.
A FORTRAN
77
code was de eloped in a VAX-2 compu e . Also, he me hod op imizes he
compu e esou ces (Time and memo y).
Recibido: Junio
1991
@Uni e si a Poli kcnica de Ca aiunya (España)
ISSN
0213-1315
4
O.A. SOTO
Y
C.
MOLANO
INTRODUCCION
El obje i o p incipal del p esen e a ículo, es p esen a las bondades del mé odo de
los G adien es Conjugados con P econdicionamien o y Escalamien o, pa a da solución
al sis ema de ecuaciones lineales esul an es de la modelación po elemen os ini os de
sondeos geoléc icos, y de p oblemas de lujo y po encial en gene al (Flujo de aguas
sub e áneas, lujo de calo , con aminación, e c
...),
Es uc u as (Flujo de es ue zos
po elemen os ini os, mé odos ma iciales de cálculo es uc u al, e c
...),
Geo ecnia,
Ingenie ía de Pe óleos, e c.
MARCO TEORICO
La ecuación di e encial que ige la dis ibución de po encial en cualquie medio
puede se ob enida a pa i del p incipio de con inuidad y de las leyes ísicas que
modelen el p oblema especí ico, p.ej. la ley de Ohm en el caso del lujo de co ien e
eléc ica y la ley de Da cy en el caso de lujo de agua en medios po osos. Es a es la
siguien e:
Ki
=
P opiedad del medio la cual da la asa de lujo en la di ección
i
en cada
pun o (x, y, z) del dominio (Q), po unidad de caída de po encial en la
di ección i-ésima en dicho pun o.
@(x, y, z)
=
Función de po encial.
W
=
Densidad de lujo po unidad de olumen en cada pun o (x,
y,
z) del
dominio de lujo (O).
Rep esen ación po Elemen os Fini os
Pa a la simulación numé ica po elemen os ini os de la ecuación an e io
gene almen e se a a el p oblema en o ma bidimensional. Es o es álido cuando
la espues a del medio es básicamen e en dos di ecciones como en el caso de sondeos
geoeléc icos, o debido a azones de sime ía como en el caso de aguas sub e áneas. La
ecuación
1
puede se en onces educida a su o ma bidimensional, la cual se p esen a a
con inuación:
en es e caso
W
ep esen a la densidad del lujo po unidad de á ea, y las demás a iables
son las mismas de la ecuación (1).
La simulación de la ecuación
(2)
po elemen os ini os ipo Gale kin, consis e
p imo dialmen e en di idi el dominio de lujo en una se ie de subdominios homogéneos
SOLUCION NUMERICA POR MEF DE SONDEOS GEOELECTRICOS
5
e iso ópicos (K,
=
K
-
K) llamados elemen os, y eemplaza la solución exac a del
"7
po encial po una solucion ap oximada de la siguien e o ma:
a;
=
Valo del po encial en cada uno de los nodos de la ed. Los nodos se localiza on
en los é ices de los elemen os, los cuales se oma on iangula es pa a es e
caso.
4;
=
Función base de in e polación del nodo
i.
Es a se omó lineal eniendo el
alo de uno en el nodo
i,
y
de ce o en los nodos adyacen es. Rep esen a el
po cen aje de in luencia o ponde ación del alo del po encial en el nodo
i
sob e
cada pun o (x,
y)
del dominio de lujo
(O).
n
=
Núme o de nodos de la ed.
Reemplazando la unción del po encial es imado ecuación
(3)
en la ecuación
(2),
se ob iene un esidual en cada pun o del dominio. El mé odo de Gale kin o za el alo
p omedio del esidual sob e odo el dominio de ' lujo a ce o, y u iliza como ac o es
de ponde ación pa a calcula dicho p omedio las unciones bases de in e polación
p esen adas en la ecuación
(3).
Haciendo es o se llega a
un
sis ema de ecuaciones
lineales simul áneas las cuales pueden se ep esen adas ma icialmen e como:
[A]
=
Ma iz de
n
x
n.
{a)
=
Vec o de po enciales de
n
x
1;
en su posición i-ésima gua da la a.~iable
a;,
las cuales son las incógni as del p oblema.
{b)
=
Vec o de in ensidad o de inducción. En su posición i-ésima gua da el alo de
~
la in ensidad del lujo inducido en el nodo
i.
Pa a las condiciones de on e a se secciona la ed de elemen os ini os en es
pa es: Zona de in e és ( ed más ina), semi-in ini os la e ales y semi-in ini o con la
p o undidad.
A
los nodos en las on e as ex e nas de la ed puede ijá seles el po encial
en un alo dado (F on e a ipo Di ichle ), o deja los lib es (F on e as ipo Newman).
También puede oma se como esul ado el alo p omedio de los po enciales ob enidos
de la simulación con las dos condiciones de on e a an e io es.
l
Mé odo de los G adien es Conjugados (Kaaschie e ,
1985)
Con el algo i mo del mé odo CG (G adien es Conjugados), básicamen e se
compu an ap oximaciones al ec o de incógni as, pa iendo de un ec o solución
inicial, y ajus ándolo de una mane a óp ima po medio del siguien e algo i mo:
O.A. SOTO
Y
C.
MOLANO
I
Algo i mo No.1
(Kaaschie e , 1985):
I
i
{ ;}
=
O
hen s op
T.
b-1
:=
({~i}
{~s})/({~;-l}T{ ;-l})(~-l
:=
o)
{Pi}
:=
i ;}
+
Pi-1{p;-1}
T.
ai
:=
{~*}/({~i}~[A]{pi})
{@;+1)
:=
{@;)
+
a;{p;}
{~i+i}
:=
{T;}
-
a;[A]{p;}
End do
Pa a ga an iza la con e gencia del mé odo la ma iz del sis ema de ecuaciones [A],
debe se simé ica
y
posi i amen e de inida; es deci :
[Al
=
[AlT
Y
{xlT[~l{x}
>
0 V{x)
#
0
El ec o { } calculado en el algo i mo an e io ep esen a el esidual; es o es:
Se hab á llegado a la espues a exac a del sis ema de ecuaciones cuando el ec o { }
sea igual a ce o, como puede se obse ado en el algo i mo No.1. Sin emba go, desde
el pun o de is a p ác ico es a condición es di ícil de log a , po lo cual el c i e io de
con e gencia o e minación del algo i mo debe se modi icado.
C i e io
de
Te minación
El c i e io de e minación pa a es a clase de mé odos i e a i os puede se de dos
ipos: Rela i o y absolu o. El algo i mo del CGEIS u iliza una combinación de es os
dos ipos como se mues a a con inuación.
El c i e io de e o absolu o dice básicamen e, que la ap oximación ob enida en una
i e ación de e minada del CGEIS es lo su icien emen e buena pa a se omada como
la solución del sis ema, si la no ma del e o (di e encia en e la solución exac a
y
su
i-ésima ap oximación) es meno o igual a un alo posi i o muy ce cano a ce o. En
la implemen ación del algo i mo se omó como c i e io de e o absolu o la siguien e
exp esión (Kaaschie e , 1985):
de donde se sigue di ec amen e,
ll{ ;}ll
5
u?)
*
e1
SOLUCION NUMERICA POR MEF DE SONDEOS GEOELECTRICOS
7
u?)
=
Ap oximación i-ésima al meno alo p opio de la ma iz [A].
el
=
Valo posi i o ce cano a ce o.
La ecuación
(7)
se implemen ó como c i e io de e o absolu o.
u ilizando un c i e io de e o ela i o se exige que la no ma del e o sob e la no ma
de la solución exac a del sis ema, sea meno o igual a un alo posi i o muy ce cano a
ce o. Se omó como c i e io de e o ela i o lo siguien e (Kaaschie e , 1985):
de donde se sigue di ec amen e,
e2
=
Valo posi i o ce cano a ce o.
Si se exigen los dos ipos de e o es, se asegu a que la solución encon ada po el
mé odo sea muy ce cana a la solución exac a (e o absolu o pequeño); y que además,
no alga la pena ajus a más la ap oximación ob enida debido a que la magni ud de
las co ecciones en una nue a i e ación se á muy pequeña (e o ela i o pequeño).
Pa a es ima el meno alo p opio de la ma iz [A], se u iliza el mé odo de la
bisección.
P econdicionamien o
y
Escalonamien o de la Ma iz del Sis ema de
Ecuaciones
Pa a acele a la con e gencia del mé odo, se p econdiciona la ma iz [A] al que:
[Á]
=
[c]-'[A][c]~
de donde [A]
=
[c][A][c]~
(10)
Si se eemplaza lo an e io en el sis ema o iginal,
y
se p emul iplica po la in e sa de
la ma iz [C], se llega a un sis ema de ecuaciones p econdicionado de la o ma:
[C]-'[A][C]-~[C]~{@}
=
[C]-'{b); o lo que es igual
[A]{&)
=
(6)
(11)
donde{&)
=
[CIT{Q) y
{b}
=
[~]-l{b)
(12)
A
la ma iz[M]
=
[C][CIT
se le llama comunmen e ma iz p econdicionado a del
sis ema.
Pa a implemen a en o ma e icien e el p econdiciona nien o del sis ema de
ecuaciones, se ealiza un o denamien o ma emá ico a pa i del algo i mo No.1, y se
llega que pa a el sis ema p econdicionado se debe inclui en las in e aciones un ec o
{z}
dado po :
{.i)
=
{Ti}
(13)

8
O.A.
SOTO
Y
C.
MOLANO
1
y se deben calcula los é minos be a y al a del algo i mo como:
Con es as modi icaciones, los ajus es del ec o de incógni as y del ec o esidual
con inuan igual
a
lo p esen ado.en el algo i mo No.1.
Pa a la ob ención de la ma iz p econdicionado a puede se u ilizada la
Descomposición Incomple a de Cholesky, la cual consis e en ob ene una ma iz
[C]
simé ica al que:
c(i, i)
=
1
pa a i,
j
=
1,
...,
n
(14)
;
c(i, j)
=
-
c(j7 8;
j=l
Pa a que es e ipo de ma iz p econdicionado a pueda se u ilizada, la ma iz
[A]
debe
se al que los é minos c(i, i) no esul en imagina ios
al
calcula las aices cuad adas
(Ve ecuación 14). El p econdicionamien o iene como obje i o asemeja la ma iz del
sis ema de ecuaciones a la iden idad, y de esa o ma acele a la con e gencia.
Con el in de aumen a la e iciencia del mé odo, se de inen dos ma ices de la
siguien e o ma:
d(i, i)
=
c(i,
i)
V
i; d(i, j)
=
O
V
i
#
j;
[E]
=
[DI-'[c]
Reemplazando es a úl ima exp esión en el sis ema p econdicionado (ll), se llega al
siguien e sis ema escalado:
donde
[A]
=
[D]-'[A][D]-~
=
[DI-'[A][D]-';
{&)
=
[DIT{@)
=
[D]{@)Y{~}
=
[DI-'{b)
(16)
La ecuación
(15)
es exac amen e igual a p econdiciona el sis ema de ecuaciones:
[A]{&)
=
{
6)
con la ma iz p econdicionado a
[M]
=
[E][E]~.
Se puede e i ica ácilmen e que los é minos de la ma iz iangula in e io de
[E]
son exac amen e iguales a los de la ma iz del sis ema ep esen ado en la ecuación (17),
y
que los é minos de su diagonal son iguales a uno; es as ca a e ís icas pe mi en aho a
almacenamien o y simpli ican el cálculo del ec o
{z)
(Ve ecuación 14), disminuyendo
conside ablemen e el iempo de ejecución en compu ado .
SOLUCION NUMERICA POR MEF DE SONDEOS GEOELECTRICOS
9
Almacenamien o Compac o
La ma iz
[A]
puede se almacenada de mane a compac a con el in de aho a
memo ia de compu ado , y de agiliza la implemen ación de las mul iplicaciones del
ipo
{w)
=
[A]{p) que deben se ejecu adas en el algo i mo.
El almacenamien o es ealizado de la siguien e mane a: se de ine el ec o
{a)
de
n
x
1,
en el cual son almacenados los é minos de la diagonal de la ma iz [A]; se de ine
un ec o {m) de
n
x
1,
el cual con iene en su i-ésima posición el núme o de elemen os
di e en es de ce o que iene la ma iz iangula in e io de [A] en sus p ime as i- ilas;
se de ine el ec o {a), en el cual son almacenados los é minos di e en es de ce o de la
ma iz in e io de [A] ila po ila, po lo que la dimensión de es e ec o se á m(n)
x
1;
inalmen e se de ine un ec o {m), el cual en su i-ésima posición con end á el alo
de la columna en el que se encuen a el é mino a(i) en la ma iz iangula in e io
o iginal. Con los ec o es {a), {a), {m)
y
{m) queda de inida comple amen e la ma iz
[Al
PRESENTACION DE RESULTADOS
Compa ación de Resul ados con Soluciones Analí icas
Con el in de e alua la bondad del mé odo de. simulación numé ica u ilizando
el algo i mo del CGEIS, se compa a on cu as de esis i idad (Va iable
1/K
de
p opo cionalidad pa a lujo de co ien e eléc ica, e ecuación (1) de modelos eó icos
con a simulaciones numé icas de las mismas, u ilizando como algo i mo de solución al
sis ema de ecuaciones, el mé odo desc i o en el p esen e a ículo.
Se simula on modelos eó icos de dos, es y cua o capas de di e en e esis i idad,
compa ándose las cu as ob enidas analí icamen e con las ob enidas po el mé odo.
Las cu as de esis i idad se calculan a pa i de los alo es de po encial en e dos
nodos de la ed; lo cual quie e deci que si se ob ienen buenas ap oximaciones a la
solución del sis ema de ecuaciones las cu as de esis i idad simuladas se asemeja án
muy bien a las calculadas eó icamen e. En las Figu as la), b)y c) se p esen an las
cu as eó icas (cu a con inua) y simuladas (cu a pun eada) pa a los modelos de dos
es y cua o capas espec i amen e. Las lineas a azos ec os que se supe ponen sob e
dichas g á icas, ep esen an el modelo eó ico simulado.
Análisis de Sensibilidad
Fue on ealizados análisis de sensibilidad con los modelos eó icos con el in de
de e mina el g ado de sucep ibilidad de los esul ados de la simulación espec o al
c i e io de e o . Se encon ó que colocando an o las cons an es del e o absolu o como
las del ela i o en un alo igual a 0.01, se ob ienen ap oximaciones muy buenas desde
el pun o de is a p ác ico. Modi ica las a alo es meno es no mejo a ap eciablemen e
la calidad de los esul ados, y po el con a io, sí aumen a los iempos de ejecución del
algo i mo.
Po o a pa e se hicie on ensayos de p ueba cambiando el c i e io de e o mix o
po un c i e io de e o absolu o de la siguien e o ma:
l{~ilI
<
e
(18)
1
0
O.A. SOTO
Y
C.
MOLANO
Resis i idad
Apa en e en
a
{Oh n-m)
2
Modelo
2
Capas
P o
(m)-Resis.(Oh n-m1
0-10
10
Sepa ación de Elec odos (m).
Figu a la. Cu as de Resis i idad Teó ica
y
Numé ica Modelo de
2
Capas.
Modelo
3
Ca as
P o
(m)-~esis.(~hm-m)
0-13 20
Resis i idad
Apa en e en
(Oh n-m)
Sepa ación de Elec odos (m).
Figu a lb. Cu as de Resis i idad Teó ica
y
Numé ica Modelo de
3
Capas.
Resis i idad
Apa en e en
(Ohm-m)
Sepa ación de Elec odos(m).
Modelo
4
Capas
P o
(m)-Resis.íOhm-m)
0-3
20
Figu a lc. Cu as de Resis i idad Teó ica
y
Numé ica Modelo de
2
Capas.
-
SOLUCION NUMERICA POR MEF
DE
SONDEOS GEOELECTRICOS
11
La espues a siguiendo es a me odología ue de igual calidad que con el c i e io mix o, y
el núme o de ap oximaciones del algo i mo se educía conside ablemen e, disminuyendo
el iempo de ejecución del p og ama. Lo an e io demues a que el c i e io de e o
mix o p esen ado es muy es ic o pa a ines p ác icos, y que és e puede se modi icado
po un c i e io de e o absolu o sin pe de calidad en los esul ados.
Simulación de una Sección T ans e sal Real
Con el algo i mo implemen ado, se e ec ua on simulaciones numé icas sob e una
sección ans e sal eal de ca ac e ís icas más o menos conocidas. Básicamen e con la
solución numé ica se ob ienen una se ie de cu as de iso esis i idad que son dibujadas
en el pe il del e eno (pseudosecciones), las cuales deben se simila es a las ob enidas
po medio de la ejecución de sondeos geoeléc icos eales. El mé odo u ilizado pa a
cons ui dichas cu as ue el WTP (Wenne T ipo encial); po medio del cual se
gene án dos ipos de seudosecciones: Una ipo Al a, o a ipo Be a sob e Gama. Se
abajó sob e una sección suminis ada po el TNO (Ins i u e o Applied Geoscience
de Holanda), el cual enía in e és en e i ica median e simulaciones numé icas las
condiciones es a ig á icas de la misma. Pa a su es udio el ins i u o mencionado
e aluó pseudosecciones u ilizando el mé odo WTP, y cons uyó pseudosecciones ipo
Al a y Be a sob e Gama a pa i de los esul ados de campo, sin ob ene esul ados
sa is ac o ios al compa a las con las ob enidas a pa i de pe o aciones manuales y
es udios de suelos. En la Uni e sidad de los Andes se ealiza on simulaciones numé icas
con la es a ig a ía de la sección en iada po el TNO, u ilizando el algo i mo p esen ado
en el a ículo.
En la Figu a 2a) puede se obse ada la es a ig a ía ap oximada de la zona de
es udio. En la Figu a 2b) puede se obse ada la pseudosección ipo Al a ob enida
a pa i de las mediciones de campo. En las Figu as 2c) y d) se mues an las
pseudosecciones ob enidas a pa i de la simulación numé ica del e eno en cues ión.
La pseudosección ipo Al a ob enida numé icamen e es muy pa ecida an o en los
alo es como en la o ma de las cu as de iso esis i idad, a la ob enida a pa i de
las mediciones de campo. Cabe ano a sin emba go, que la pseudosección Be a sob e
Gama ob enida en el e eno, no es muy pa ecida a la ob enida con la simulación
numé ica; y que ampoco sigue la es a ig a ía que se supone iene el pe il del suelo.
Es o hace pensa que en ealidad el e eno posee una con igu ación es a ig á ica un
poco di e en e a la modelada numé icamen e, p esen ando algunos len es de ma e ial o
calidades de agua aún no iden i icados, que a ec an la pseudosección Be a sob e Gama
en g an medida, pe o muy poco a la pseudosección ipo Al a. Lo an e io es posible si
se iene en cuen a que las seudosecciones ipo Al a mues an una esis i idad apa en e
p omedio del e eno, po lo cual len es de ma e ial de espeso pequeño pueden no
llega a se egis ados; y po el con a io, las pseudosecciones ipo Be a sob e Gama
ienden a segui la con igu ación es a ig á ica del e eno, po lo que cualquie len e
de ma e ial in luye no ablemen e en la o ma de la pseudosección.
--