Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Solución transporte conservativo (1)
 Discretización
especial (2D)
x
j+1
Diferencias finitas
y
j
j-1
i-1
c

x
c
i 12
c
i 12

x
i, j
c
i 1
c
i 1
c
 c
2
2x
x
c
i
x

2
i 1
c c
i

i
i+1
i 1
x

x
i, j
c
i 1
 2c  c
i
x
i 1
2
Ecuación de transporte para nudo i,j

c
t
 qx
 qx
c
i 1, j
c
x
c
x
 qy
c
y
i 1 , j
 qy
 c
2
 D xx
c
i , j 1
x
c
y
2
 c
2
 D yy
y
i 1 , j
 D xx
c
2

i 1
 2c  c
i
x
i 1
  D yy
c
i , j 1
 2c
i, j
y
c
i , j 1
Solución transporte conservativo (2)
Elementos Finitos
Se basa en interpolación entre
nudos
Se especifica coordinadas de
nudos y los nudos de cada
elementos
Se pueden mezclar tipos de
elementos
Segmento
Triángulo
Valor (c)
Nudo
Elemento
Segmento 1D
Triángulo
Cuadrilátero
Cuadrilátero
Malla 2D de EF
Tetraedro
Prisma
Solución transporte conservativo (3)
 Para
todos los nudos y en notación de matriz:
Vector de concentraciones para cada nudo
Matriz para
adv./dif./disp.
E a c  Fa
c
t
 ga
Vector de fuentes/sumideros
Matriz diagonal de almacenamiento
 Discretización
temporal por DF
E a  c k 1  1   c k   F a
c k 1  c k
 ga
t
Factor de ponderación temporal (entre 0 y 1)
 Resolver
ck+1 mediante:
Fa 
Fa 





E

c

g



1
E


 k 1

c k
a
a
a
t 
t 


c
t

c k 1  c k
t
Resolver sistema lineal
 Métodos
Ax  b
directos
Descomposición en LU
Dos partes
Descomposición: ‘prepara’ matriz A
Solución: calcula x
La descomposición es más costosa que la solución
Si cambia b pero no cambia A, se puede aprovechar la
descomposición
Es un método robusto
Resolver sistema lineal (2)
 Métodos
iterativos
Ejemplos: Gradientes Conjugados, GMRES
Empiezan con una solución inicial, que se mejora
cada iteración
Pueden tener problemas de convergencia
Requieren una solución inicial y criterios de
convergencia, p.e. i  1
i
i 1
x
  Ax
x
max

b
max
Son mejores para mallas de 2D y 3D con muchos
nudos
Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Solución de las ec. de transporte reactivo
 Problema
Muchas ecuaciones/incógnitas: número de nudos
 número de componentes
En general, ecuaciones muy no lineales
 Para
problemas no lineales hay dos métodos
de resolución
Picard (en transporte reactivo también llamado:
SIA (Sequential Iteration Approach), Two-step)
Newton Raphson (en transporte reactivo también
llamado: DSA (Direct Substitution Approach), Onestep, Global Implicit)
Picard, principio
ecuación como A ( x ) x  b ( x )
 Resolvemos iterativamente xi+1 un sistema
lineal hasta que converja:
 Escribimos
i
A (x )x
 Para
i 1
 b(x )
i
1 incógnita
i = número de iteración
b(x)
x1
A(x1)x
A(x2)x
x2 x3
A(x)x
Criterios de convergencia
 Error
x
i 1
x
absoluto de la incógnita
i
max
 Error
x
i 1
x
x
relativo de la incógnita
i
  x , rel
i 1
max
 Error
Ax
i 1
  x , abs
de la ecuación
b
i 1
max
  ec
Picard, ejemplo
 Dos
incógnitas (x1 x2) y dos ecuaciones
 2 x1  x 22  x1 x 2  1
 3
2
x

x
 3 x2  2
1
 1
2

 2
 x1  x1
x 2   x1   1  x1 x 2 
    

3  x2  
2

Iter.
1
2
..
8
x
0.00
0.50
..
0.22
0.00
0.67
..
0.65
1.00
0.67
..
0.86
1.00
2.00
..
2.00
b
A
2.00
0.00
2.00
0.67
..
2.00
0.65
0.00
3.00
0.75
3.00
..
0.27
3.00
Ax-b -1.00
0.78
..
0.00
-2.00
0.38
..
0.00
Picard, divergencia
 En
lugar de converger (cada vez más cerca de
la solución), puede divergir (cada vez más
lejos de la solución)
 Se puede
solucionar
b(x)
A(x)x
refinando la
descretización
(disminuyendo x
o sobre todo t)
x3
x1
x2
Newton-Raphson, principio
Escribimos ecuación como
f (x)  0
 Resolvemos iterativamente xi+1 mediante un sistema
lineal hasta que converja:

J x
i
i 1
x
i
  f
Jacabiano

i
Residuo
Para 1 incógnita
f
J0
J1
x3 x2
x1
x0
  f1
 i
  x1
 f 2
f
i
 x i
J 

i
x
 1
 
 f n
 x i
 1
 f1
x
f 2
i
2
x

f n
i
2
x
i
2




 f1 

i
xn 
f 2 
i
xn 

 
f n 
i
 x n 
Newton Raphson, ejemplo
 2 x1  x 22  x1 x 2  1
 3
2
x

x
 3 x2  2
1
 1
 2 x1  x 22  x1 x 2  1 
 2  x2
 J  2
f  3
 3x  2x
 x  x2  3x  2 
1
 1
1
2
 1

Iter.
1
2
..
5
x
0.00
0.50
..
0.22
0.00
0.67
..
0.65
1.00
-0.78
..
0.00
2.00
-0.38
..
0.00
-f
J
xi+1-xi
2.00
0.00
2.67
1.83
..
2.65
1.51
0.00
3.00
1.75
3.00
..
0.58
3.00
0.50
-0.34
..
0.00
0.67
0.08
..
0.00
2 x 2  x1 

3

Newton-Raphson, divergencia
 También
Newton-Raphson puede divergir
f
x0
x1
 También
x2
se puede solucionar refinando
descretización
Pseudo Newton-Raphson
 No
(siempre) se actualiza el jacobiano
 Requiere más iteraciones, pero se ahorra
tiempo de cálculo en el ensamblaje del
jacobiano
f
J0
x3 x2 x1
x0
Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Si hay expresión explícita para química

Nc ecuaciones y Nc incógnitas (c1)
T
 c1 
*
f  Ua 
  u a  c1   S a  c a 2 (c1 )  u a  0
 c a 2 (c1 ) 
1: Se sabe la concentración total (ua)
f  c1  c1
2: Se sabe la concentración (cprescrita)
0
prescrita

f  c a 2 (c1 )  c 2
prescrita
or
Newton-Raphson
J
i
c
i 1
1
1: J 
i
c
f
c
2: J i  I
i
1
i
1
  f
 I  S
or
*
a
J 
i
0
log c a 2  S a log c 1  log k a
*
*
i

T
c a 2
 c1
i
c a 2
 c1
i
ca 2, j
 c1, k

c a 2 , j  ln c a 2 , j
 ln c a 2
c1, i  ln c1, k
 ln c 1
i
 Sa
*
Ejemplo
Primarias: H+, CO3-2
 Secundarias: HCO3-, CO2, OH Carbono total = 10-3, pH =7

 [C O 3-2 ]  [H C O 3- ]  [C O 2 ]  10  3   0 
f 
 
+
7
[H ]  10

 0
log c a 2
1

1

0

c a 2

[H C O 3 ] [C O 2 ]

1 
22J 
[C
O
]
[C
O
]
3
3


0

-
[H C O 3 ]
+
[H ]
[C O 2 ] 
2

+
[H ] 

1

 c1
 log[H C O 3- ] 


  log[C O 2 ]  
 log[O H - ] 


1 
 log k 1
2  log[C O 3 ]  
2 
 log k 2
 log[H + ]  

 
 1 
 log k 2
 [H C O 3- ]

2[C
O
]
3

 [C O ]
2

2 [C O 3 ]


0


[H C O 3 ] 

+
[H ] 
[C O 2 ] 
2

+
[H ] 
- 
[O H ]


+

[H ] 
-





Si no hay expresión explícita para química

Ecuaciones químicas
log c a 2  log γ a 2  S a log c 1  γ a 1   log k a
*
*
2
log  i  

Az i
1  Ba i
I
 bI
I 
1
2

2
ci z i
i
I
¿Cómo calcular ca2 y ca2/c1?
Se itera
fc
c a 2
c
j 1
a2
Después de la última iteración
c
j
a2

 fc
fc c a 2
 c a 2  c1

fc
dfc
 c1
d c1
Regla de la cadena
0
fc c 2
 c 2  c1

fc
 c1
O, montar Newton-Raphson con Ns incógnitas (c) y
ecuaciones (componentes y químicas)
 O, calcular ca2 y ca2/c1 iterativamente (i+1 = (ci))

Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Picard, aplicación a transporte reactivo
 Escribimos
u a
t
ec. de transporte reactivo como:
 L (u a )   U s
Términos lineales
 Primer
c s
t
 Um
c m
t
 US k rk ( c )  b
t
Términos no lineales
paso: transporte
Se puede resolver para cada componente por
i 1
separado
u a
i 1
i
t
 L (u a )  b
Ecuación discretizada
F a  c ,i 1
Fa  c


c
c ,i



E

u

g



1
E

u

b




a
a , k 1
a
a ,k

t

t




c
u a = vector de conc. acuosas de componente c para todos los nudos
Picard, aplicación a transporte reactivo (2)
 Segundo
paso: química
Cálculo de bi+1 a partir de uai+1 para cada nudo por
separado
Resolver c1 (Nc) y cm (Nm) aplicando Newton-Raphson
(pequeño)
u a ( c 1 )  u d ( c 1 )  U m c m  u a  u d  u m   0
i 1
i
i
log a m  0  S m log a 1  log k m
*
*
Nm
Calcular bi+1 a partir de c1 y cm
b
i 1
 U s
 c s (c1 )
t
 Um
c m
t
 US k rk ( c ( c 1 ))
t
Nc
Detalles ecuaciones químicas
u a ( c 1 )  u d ( c 1 )  U m c m  u
i 1
a
u u
i
d
i
m
 0
2
log  i  
Az i
I
 bI
1  Ba i
I 
I
1
2
cz
i
 c1 

u a  U a 
 c a 2 (c1 ) 
log c a 2  log γ a 2  S a log c 1  γ a 1   log k a
u d  U d c d (c1 )
log c d  S d log c 1  γ a 1   log k d
i
*
*
log a m  0  S m log a 1  log k m
*
*
Ec. de transporte tratada
como ec. química
Para componentes inmóviles, p.e.
[XNa] + 2[X2Ca] = CEC
U s c s  cte
*
*
2
i
Ejemplo

Reacciones en equilibrio
R1: CO32- = HCO3- - H+
R2: X2Ca = 2XNa + Ca2+ - 2Na+
R3: H2O = H+ + OHR4: CaCO3(s) = Ca2+ + CO32-

Componentes:
Ca2+, HCO3-, Na+, XNa (CEC), H+, OH-
Tratar como ecuación 'química'
Ejemplo, paso de transporte
 Resolver
las concentraciones acuosas
F a  Ca , i  1
F a  Ca


Ca
Ca , i



E

u

g



1
E

u

b

 a , k 1

 a ,k
a
a
t 
t 


F a  HCO 3 , i  1
F a  HCO 3


HCO 3
HCO 3 , i
g
    1 E a 
E a 
 u a , k 1
u a ,k  b
t 
t 


Fa

E a 
t

Fa
 Na , i  1

Na
 u a , k  1  g     1 E a 
t


 Na
Na , i
u

b
 a ,k

Fa  H ,i 1
Fa  H


H
H ,i
E a 
 u a , k  1  g     1 E a 
u a ,k  b
t 
t 


Fa

E a 
t

Fa
 OH , i  1

OH
    1 E a 
 u a , k 1  g
t


 OH
OH , i
u

b
 a ,k

Ejemplo, paso químico

Resolver Ca2+, HCO3-, Na+, H+, OH-, XNa, CaCO3
2
[ Ca
2
]
[ XNa ] [ Ca
2

]
2
[ Na ] K 2
-

[ HCO 3 ]
[ HCO ] 
3

2
2
[ Na ]  2
[ XNa ] [ Ca

]

[H ] 

[ H ][ OH ]


[ OH ] 
[ HCO 3 ]

[H ]K 1

[ H ][ OH ]
K1
[ XNa ] [ Ca
2

HCO , i
HCO , i
 0

OH , i  1
Na , i


H ,i 1

 [ CaCO 3 ]  u a , k 1  u d , k 1  u m , k 1  0
H ,i
H ,i

 u a , k 1  u d , k 1  u m , k 1  0
2
[ XNa ]  2
Ca , i
 u d , k 13  u m , k 13
Na , i
-
K3

Na , i  1
Ca , i
 u a , k 1  u d , k 1  u m , k 1  0
2
[ Na ] K 2

HCO , i  1
 [ CaCO 3 ]  u a , k 13
[H ]K 1

Ca , i  1
 [ CaCO 3 ]  u a , k 1  u d , k 1  u m , k 1  0
[ Na ] K 2
2
]
OH , i
 CEC  0
OH , i
[ Ca
2
-
][ HCO 3 ]

[H ]
 K4  0
Picard, variantes
 Utilizar
u en lugar de ua como variable
Paso de transporte
u
i 1
t
 L (u
i 1
)   U d L ( c d )  U m L ( c m )  US k rk ( c )  b
i
i
t
i
i
Paso químico
u a (c1 )  u d (c1 )  U m c m  u
i 1
0
log a m  0  S m log a 1  log k m
*
 SNIA (Sequential
*
Non Iterative Approach)
Lo mismo que SIA pero sin iterar
Valores iniciales
 Paso
de transporte: utilizar el valor del tiempo
anterior
k  1, i 1
ua
 ua
k
 Paso
químico: utilizar el valor de la iteración
de transporte anterior
c
k  1 , i , j 1
1
c
k  1 , i 1
1
i = iteración de transporte
j = iteración química
Picard y eliminación de minerales
 La
presencia de minerales puede depender
del espacio  E y EU dependen del espacio
 Se pierde la ventaja de calcular el paso de
transporte para cada componente por
separado
 La eliminación de minerales sólo se puede
incorporar a Picard, si la presencia de
minerales no depende del espacio.
Newton-Raphson, aplicación a trans. react.
 Escribimos
f t  EU
 Hace
ecuación de transp. react. como
 c o ( c 1 .1 )
o
t
 EU o L ( c o ( c 1 .1 ))  EUS k rk ( c ( c 1 .1 ))  0
t
falta
Discretizar
Jacobiano
 Si
las conc. secundarias se escriben
explícitamente en función de las primarias,
p.e. log c  S log c  log k
se pueden escribir las ec. de transporte en
función de conc. primarias
*
2
*
1
Newton-Raphson, sustitución ec. químicas
 Si
ecuaciones químicas no son explícitas
 f t (c1 , c 2 )  0

f c (c1 , c 2 )  0
(Nc-Np) ecuaciones de transporte
(Nr) ecuaciones químicas
Aplicar Newton-Raphson a ecuación de transporte
 ft
 f t  c 2  i 1

 c 1  c 1i    f t

 c


c

c
2
1 
 1
Newton-Raphson pequeño para química
Se itera
f c
c 2
c
j 1
2
Después de la última iteración

 c 2  fc
j
f c c 2
c 2 c1

f c
 c1
Valores iniciales
 Para
el valor inicial utilizar el valor del tiempo
anterior
k  1, i 1
c 1 .1
 c 1 .1
k
Contenido
 Resolver
ecuación de transporte conservativo
y sistemas lineales
 Resolver sistemas no lineales
Teoría
Aplicación de Newton-Raphson a especiación
Aplicación de Picard y Newton-Raphson a
transporte reactivo
Comparación entre dos métodos
Comparación teórica
Picard, SIA
Resuelve transporte y
química por separado
Requiere menos
memoria de ordenador
Más fácil de programar
Se usa mucho
Convergencia lineal
Más rígido
Newton-Raphson, DSA
Resuelve transporte y
química simultáneamente
Requiere más memoria
de ordenador
Más difícil de programar
Se usa poco
Convergencia cuadrática
Más robusto
Comparación mediante ejemplos
 Metodología
Calcular ejemplos
mediante los dos
métodos
Usar gestión
automática de t
Comparar número de
incrementos de
tiempo, número de
iteraciones y tiempo
de CPU
ca lcu la tio n s
fa ile d to
co n verg e
su cce ssful
co n verg e n ce
n ite r < th rm in
n ite r
n ite r > th rm a x
th rm in < n ite r < th rm a x
d t = d t/fd
d t = d t*fi
ye s
d t > d tm a x
d t = d t/fd
no
d t = d tm a x
n e xt
tim es te p
Ejemplo CAL
Disolución de CALcita
C A L -E
C A L -4
10
p H p ro file s a fte r 5 y rs
9
C A L -3
C A L -2
8
C A L -1
pH

7
C A L -0
6
5
 = 0 .1
 = 10 m
q = 2 m /y r
0
B a s ic g rid : 2 0 1 D e le m e n ts
o f e q u a l le n g th
50
100
D is ta n c e (m )

volumen
3 componentes
7 acuosas
1 mineral
21 nudos
1.0 volumen poros lavados
de agua entrado/sa lido
volumen
de poros

tiempo simulado
tiempo de tránsito
Ejemplo WAD
1E+000
q = 1 5 m m /y r
1 E -0 0 1
CO3
1 E -0 0 2
Na
B re a k th ro u g h c u rv e s a t
2 m fo r c a s e w ith o u t
c a lc ite (W A D -0 )
CO3
1 E -0 0 3
Mg
1 E -0 0 4
Ca
1 E -0 0 5
0
Na
1 E -0 0 6
1 E -0 0 1
CO3
Na
1 E -0 0 2
 = 0 .3
C E C = 6 1 m e q /l
 = 0 .1 m
Cl
1 E -0 0 7
B re a k th ro u g h c u rv e s a t
2 m fo r c a s e w ith c a lc ite
in e q u ilib riu m (W A D -E )
CO3
1 E -0 0 3
Mg
1 E -0 0 4
Ca
1 E -0 0 5
Na
1 E -0 0 6
1
Cl
2
1 E -0 0 7
0
400
800
T im e (y e a rs )
1200
1600
D e p th (m )
1E+000
B a sic g rid : 2 0 1 D e le m e n ts
o f e q u a l le n g th
T o ta l a q u e o u s c o n c . (m o l/l)
Intercambio iónico en el ‘WADdenzee’
T o ta l a q u e o u s c o n c . (m o l/l)

6 componentes
9 acuosas
3 adsorbidas
1 mineral
21 nudos
37.5 vol.por.lav.
Ejemplo DEDO

DEDOlimitación cerca de una fractura
Y ( m)
0 .0 1 0
V o lu m e fra c tio n o f d o lo m ite
a fte r 2 0 0 y e a rs fo r k in e tic
c a s e (D E D O -K )
0 .0 0 5
0 .9
V o lu m e fra c tio n o f d o lo m ite
a fte r 2 0 0 y e a rs fo r e q u ilib riu m
c a s e (D E D O -E )
0 .0
V o lu m e fra c tio n o f c a lc ite
a fte r 2 0 0 y e a rs fo r e q u ilib riu m
c a s e (D E D O -E )
7 componentes
15 acuosas
2 mineral
225 nudos
22 705 vol.por.lav.
0 .0
0 .0 0 0
Y (m )
0 .0 1 0
V o lu m e fra c tio n o f c a lc ite
a fte r 2 0 0 y e a rs fo r k in e tic
c a s e (D E D O -K )
0 .0 0 5
1 .0
0 .0 0 0
Y (m )
0 .0 1 0
B a s ic g rid :
15x15 = 225 nodes
3 9 2 e le m e n ts
0 .0 0 5
q = 4 .5 m /y r
M a trix :
 = 0 .1
 = 0 .1 m
D (m o l. d if.) = 0 .0 1 m 2 /y r
q = 4 5 0 m /y r
0 .0 0 0
0 .0
0 .2
0 .4
0 .6
X (m )
0 .8
1 .0
0 .0
0 .2
0 .4
F ra c tu re :
w id th = 0 .1 m m
 = 0 .9
 = 0 .1 m
D (m o l.d if.) = 0 .0 1 m 2 /y r
0 .6
X (m )
0 .8
1 .0
Ejemplo OSA
 Meteorización
en una mina de Uranio de
OSAmu Utsumi (Poços de Caldas, Brasil)
13 componentes
42 acuosas
8 mineral
101 nudos
40 000 vol.por.lav.
D
D
E
E
W
W
W
W
W
W
E
4
3
2
1
0
A
-E
-K
-E
-4
-3
-2
-1
-0
S
O
O
D
D
D
D
D
D
D
D
L-
L-
L-
L-
L-
L-
O
A
A
A
A
A
A
A
A
C
C
A
A
A
A
8
6
4
2
0
D
D
E
E
W
W
W
W
W
W
E
4
3
2
1
0
A
-E
-K
-E
-4
-3
-2
-1
-0
S
O
O
D
D
D
D
D
D
D
D
L-
L-
L-
L-
L-
O
A
A
A
A
A
A
A
A
A
A
C
C
C
C
L-
12
A
A
E
E
A
-E
-K
-E
-4
-3
-2
-1
-0
S
O
O
D
D
D
D
D
D
4
3
2
1
0
E
L-
L-
L-
L-
L-
L-
O
D
D
A
A
A
A
A
A
A
A
A
A
A
A
0
C
C
10
L o g n u m b e r o f ite ra tio n s
D
D
W
W
W
W
W
W
C
C
C
C
C
C
L o g n u m b e r o f tim e s te p s
8
C
C
C
C
L o g C P U tim e (s )
Resultados
10
S IA
DSA
6
4
2
12
10
8
6
4
2
0
Influencia malla
 ¿Qué
pasa con mallas finas y de más
dimensiones?
Cambiar las mallas 1D en 2D
Variar el número de nudos
G rid w ith sm a lle r
n u m b e r of n o d es
G rid w ith la rg er
n u m b e r of n o d es
lo g tie m p o d e c á lc u lo (s )
Resultados
12
C A L -E
W A D -E
D E D O -E
C A L -3
W A D -3
D E D O -K
OSA
8
4
0
lo g tie m p o d e c á lc u lo (s )
1
12
2
3
4
S IA c a lc u la d o
8
S IA m e d id o
D S A c a lc u la d o
D S A m e d id o
4
0
1
2
3
4
lo g n ú m e ro d e n u d o s
5
1
2
3
4
lo g n ú m e ro d e n u d o s
5
1
2
3
4
lo g n ú m e ro d e n u d o s
5
5
Descargar

Solución de las ecuaciones matemáticas