'
P ,La Ecuación de Poisson-Boltzmann en Dos Dimensiones para
Geometría Cilíndrica.
Tesis que presenta
/Jesus Enrique Díaz Herrera
Para la obtención del grado de
/Maestro en Física
Mayo 1 9 8 6 ~
I
. Universidad Autónoma M e t r o p o l i t a n a - I z t a p a l a p a ADivisidn de Ciencias e Ingenieria
A Cractela fuente de vita/ estimulo
A mi madre y hermanos
Agradecimientos:
Deseo agradecer al M. en F. J. Sanchez Sanches y al Dr. M Lazada Cassou
por su valiosa dirección en la elaboración de este trabajo.
A la Universidad Autónoma Metropolitana-Isbapalapa el apoyo en la realisación
de este trabajo.
Agradeaco tambien al Fis. E. Gonsiles Tovar sus comentarios y la importante
ayuda para llevar a cabo las wmparciones realidas, al Dr. F. Uribe por ms valiosos
comentarios que permitieron aclarar algunos puntos y al Dr. MAntonio Duran por
BU apoyo para la solución de los sistemae de ecuaciones
INDICE
Introducción .".i
CAPITULO I
A).- La Ecuación de Poisson-Boltamann
B)- Distribución de Iones Alrededor de un Polielectrolito
C).- Diatribución de Carga Alrededor de un Cilindro Cargado
CAPITULO I1 A)- El Método de Elemento Finito
B)- El Cilculo Variacionaf
C).- El Método de Rayleigh-Rib
DI- El Método de Galerkin
E)- Funciones de hterpolación y Elementos Finitos
il- Elementos Triangulares
F).- Solucidn a la Ecuacidn de Poisson-Boltxmann Lineal
G).- Solución a la Ecuación de Poisson-Boltzrnann
CAPITULO III
Resultados
Discusión de Resultados y Conclusiones
Referencias
Bibliografía
."4
.... 7
""10
.I 13
,.59
""21
"-23
.." "-27
.... 31
""34
.,.52
"59
.... 61
INTRODUCCION
Existen varios métodos de la Mminica Estadietica para estudiar I a n propiedades
de siatemas en estado liquido, uno de ellos es el método de funciones de distri-
bucidn, que consiste en construir un sistema de ecuaciones integrales en términos
de dichas funciones de distribucion radia1,y que dan resubados aproximados.
A saber: Hypernetted Chain Aproximation (HNC), Percus-Yevick (PY), Mean Spherical
Model [MSM), Born-Creen-Yvon (BGY), Kirkwood-Poirier (KP), etc.
Aplicar estas teorías para evaluar las propiedades termodinimicas de soluciones
electroliticas y polielectroliticas, que son de gran interés en varias ramas como:
Bioquirnica, Electroquirnica, Biofisica y Química de Coloides, ha sido una trabajo
reciente que sólo fue posible hasta que se desarrollaron métodos numéricos que
permitieron resolver dichas ecuaciones integrales.
Sin embargo la primera teoria que permitió estudiar las propiedades termodinámicas
de una solución electrolitica es la bien conocida Teoría de Debye-Hückel (DH)
(1920), que consiste en resolver la ecuación linealizada de Poisson-Boltmann (PBL).
Es bien conocido que, la ecuación de Poisson-Boltzmann no linealiaada (PB) es
inconsistente con las bases de la M&nica Estadistica, pero se ha mostrado
ampliamente su validez comparando con resultados de otras teorias como HNC tanto
como datos de experimentos en DNA y cilculos de Monte Carlo, que para bajm
concentraciones de electrchito resulta ser una muy buena aproximación
La ecuación de PB es una ecuación diferencial no lineal con condiciones a la frontera,
la no linealidad de esta ecuacidn en general no permite solucidn analhica.
Sin embargo, para problemas con geometría simple (una dimensión y problemas con
simetría esférica o cilíndrica) dicha ecuación es posible reoolverla aproxi-
1
madamente utilisando métodos numéricos como: el método de iteración de Picard (en
donde la ecuación diferencial se resuelve en su forma integral) y el método clásico
de Runge-Kutta.
Sin embargo para geometrías mas complicadas, por ejemplo: cilindros curvos ó bicilindros
que son bidimensionales, los métodos numéricos mencionados no parecen ser los más
adecuados. Hasta donde sabemos la ec. de PB no ha sido resuelta para este tipo de
geometrias debido a la dificultad de obtener su solución.
Nuestra intención en esta tesis es explorar la posibilidad de resolver este tipo
de sistemas bidimensionales, tratando un problema intrinsecamente unidimensional
[un cilindro) con las ecs. de PB y PBL y resolviendolas con el método de elemento
finito (méhdo de Rayleigh-Rits-Galerkin) en dos dimensiones.
En el capítulo I, presentaremos la ecuación de Poisson-Boltamann y su aplicación
al problema de polielectrolitos cilíndricos ampliamente estudiados por varios
autores.
En el capítulo I1 se explicará la técnica de Elemento Finito para resolver ecuaciones
diferenciales con condiciones a la frontera en do6 dimensiones utilizando los métodos
de Rayleigh-Rits-Galerkin con una aproximación lineal
En el capitulo 111 damos los resultados de aplicar Elemento Finito en la
solución a la ecuación de PB y PBL para determinar la distribución de iones
alrrededor de un cilindro infinito .uniformemente cargado y compararemos con los
resultados de E.Gondea, M . L o d a y DHenderson.
En la comparación mostraremos, que la solución en dos dimensionea para la función de distribucidn radial, es tan precisa como la obtenida con el mQtodo de iteracidon
de Picard en una dimensión, pero que, debido a la bidimensionalidad del problema el
número de ecuaoiones algebráicas que es necesario resolver para obtener esa precisión,
es grande.
Esto representa la principal dificultad encontrada y para superarla, eo necesario
que en aplicaciones futuras se incluyan aigunas modificaciones que d i v u f' tremos
en esa sección.
3
CAPITULO I
A) La Ecuación de Poisson-Boltxmann.
Para estudiar la estructura microscópica de una solución electrolitica,
Debye y Huckel' desarrollaron una teoria que consiste en calcular E
en 1923 !i potencia
e\ectrost$t,ico promedio alrededor de un i6n fijo. En esta teoría hay 3 ingredientes
principales:
1). El Modelo de la Solución Electrolitica.
En este modelo el ión central es considerado como una esfera dura con carga eléctrica
en su centro y con una constante dielktrica igual a la del solvente, de la solución
que es considerada como un medio dieléctrico uniforme, y el resto de los ionea como
cargas puntuales. Eete modelo puede ser extendido a electrodos metilicos o a membranas
biológicas pero considerando al elecrodo o a la membrana como pared rigida con carga
eléctrica en su superficie y los iones en la solución (donde el sistema esti inmerso)
como cargas puntuales. En el caso de una proteina esta puede pensarse como una
partícula grande cargada y rodeada cargas puntuales.
2). La ecuación de Poisson-Boltsmann*
La ecuación de Poisson-Boltmnann, como su nombre lo indica, es una ecuación que
resulta de consideraciones electrostáticas y mecánico-estadísticas y permite describir
la distribución de carga (cargas puntuales) de un electrolito en presencia de un campo
externo, cuya fuente es la carga de una partícula fija o la carga depositada en la
superficie de una particula grande de referencia.
En la teoría de Debye-Hücltel primero se considera el potencial en un punto r como:
""( 1) 1
donde la suma es sobre todos l o s iones en el sistema.
Este potencial es promediado canónicamente conservando a una partícula fija en el
origen y puede mostrarse que dicho potencial promedio satisface la relación:
v'( <+Ir)>) = -(4~/c)<p(r)>
4
donde <p(r)> es la densidad de carga promedio, que puede emxibirse como:
""( 2)
Donde P k es la densidad numérica de los iones del tipo k y gk es la función
de distribución radial de los iones del tipo k alrrededor del ión fijo en el origen.
La función de distribución radial en (2) puede ahora expresarse en kérminos del potencial
de la fuena promedio para obtener: n
V*(<*(r)>) = (-4n/&) pk q k EXP(-pwk(r) 1 ....( 3) k=l
La primera aproximación en la teroria de Debye-Huckel es considerar el potencial
de la fueraa promedio ~2 igual al potencial electrostático promedio
w(r) = cpk(r) q k con pk E <*k(r)> para r>a ....( 4)
y sustituir (4) en (3) para obtener la ecuación de Poisson-Boltsmann:
Vz$dr) = 4 4 ~ / & 1 2 /Jk q k EXP(-hkdrl 1 para r)a k:l
si: q k = &e donde & es la Valencia de los iones del tipo k y e la
carga elemental, tenemos finalmente: n
V2 q(r) = (-4~/c) EXP(-Zke/? (PCr) 1 para r>a ....( 5) k= I
Que es una ecuación diferencial no lineal para el potencial promedio, con condiciones
a la frontera en cp o en su derivada.
b t a ecuación fue estudiada por Fowler3, Onager' y Kirkwood6 quienes mostraron
que dicha ecuación es inconsidente con las bases de la Mecánica Estadistica
3). Linealisación de la ecuqción de Poisson-Boltsmann.
La segunda aproximación para completar la teoría de Debye-Hückel consisbe en, linealirar
el lado derecho de (5) desarrollando la exponencial de la siguiente manera:
k=t k*l kr\
donde el primer sumando en la segunda igualdad se anula por electroneutralidad
y se desprecian términos de orden superior. Por lo tanto obtenemos:
Vz cpfr) = K;' cplr) ....( 6)
5
donde:
Esta ecuación fue propuesta por Kirkwood-Poirier', y por mucho tiempo Fue así,
una ley vBlida y límite de teorías correctas que describieran las propiedades
termodinimicas de una solución electrolítica, cuando la concentración tendiera
a cero. Sin embargo recientemente' se ha mostrado que para el modelo primitivo
de electrolito, la ecuación de Poismn-Boltsmann no lineal e ~ l mejor que la ec. 6 ) para bajas concentraciones, al compararse estas ecuaciones con teorías como HNC y
simulaciones de Monte Carlo.
B) Distribucidn de Iones Alrrededor de un Polielectrolito
En Bioquímica, Biofisica o Quimica de Coloides un buen número de sistemas pueden
pensarse compuestos de una partícula grande (partícula coloidal,molécula, eleckrodos,
etc,.) con alguna geometría (cilíndrica,esférica,etc.) inmerm en una solución de ionerr
pequeños. Estos sistemas pueden modelarse de una Forma simple suponiendo que los iones
son cargas puntuales, el solvente como un medio dieléctrico uniforme y a la particula
grande como una fase sólida y rígida hecha de un material con la misma constante
dieléctrica que la del mlvente.
Aunque otros modelos pueden s e r usados, como por ejemplo: bmar en cuenta el tamaño
de loa iones'l, introducir efectos debido al solvente', tomar contribuciones de,
potencial de corto alcance que describan fuersas tipo Van der Waals, fuerzas de
hidratación', etc., el m i s utiliaado es el modelo descrito arriba al que FR le agrega
un aspecto adicional, consistenk en considerar una distancia de máximo acercamiento
de l o s ionea a la pared cargada del polielectrolito. &to a equivalente a considerar
que 10s iones interaccionan con el polielectrolito como esferas duras cargadas
y entre ellos como iones puntuales. Este modelo es conocido como el modelo de Stern. Este modelo se ha utilizsdo para estudiar estos sistemas empleando la ecuación de
Poisson-Boltzrnann no lineal que result3 ser una muy buena aproximación para
concentraciones bajas de electrolito. Sin embargo e8 debido a la no linealidad de la
ecuacidn de Poisson-Boltmann lo que hace dificil encontrar soluciones analitim y sólo
soluciones numéricas para geometrías sencillas han permitido su e~aluación.'~~''
C) Distribución de Carga Alrrededor de un Cilindro Cargado.
Aplicaciones de la ecuación de Poineon-Boltsmannl' pueden encontrarse en un buen
número de sistemas, tales como: proteinas fibrosas, microelectrodos, moléculas
polielectroliticas como el DNA" y recientes modelos de contracciones musculares.
Estos siskmas usualmente mn modelados como un cilindro infinito cargado inmerso
en una solución iónica Fig( I ) y han sido eehdiados por varios autores, entre
ellos Brenner y M~Quarrie'~ quienes en 1972 dan la solución analltica siguiente a
la ecuación de Poisson-Boltsrnann lineal.
cP(r) = - ( 4 7 ~ R / c ) LN(r/R) R I r 5 (R+a/2)
cp(r) = A KO (m) rZ R+a/2
Fig (1) Representación esquemática del modelo.
7
donde:
yq, = potencial sobre el cilindro
= densidad de carga en la superficie del cilindro
R = radio del cilindro
a/2 = radio de los iones o distancia de máximo acercamiento
KO = función de Bessel modificada de orden cero
Las constantes A y U pueden ser calculadas igualando ( 8 ) J ( 9 1 para obtener:
Q = c p * / 4 7 m l { ( K o ( a ) / aKt(a) - LN(1+a/2Rll
donde:
a = ~;(R+a/2)
K, = función de Bessel modificada de orden uno.
En cuanto a la solución a la ecuación de PB para eske sistema es p o r primera
vez evaluada en 1975 por Stigter", mostrando en 19781b medianbe comparaciones con
datos experimentales en DNA que, a bajas concentraciones la ecuación de
Poisson-Boltzmann lineal lleva a substanciales errores y que la ec. de PB es mejor.
En 1983 M.Losada'6 empleando la aproximación HNC a este sistema, muestra que en
el limite del radio de IOEI iones tendiendo a cero la aproximación HNC recupera la
ecuación de PB y posteriormente EGonsiles et alJ7 demuestran que PB en el límite
de bajas concentraciones resulta tan aproximada como los ciilculos hechoa con HNC y que sólo en el caso de 'concentraciones y densidades de carga sobre el cilindro
muy bajas PI3 y PBL son aproximadamente iguales.
Comparaciones posteriores utilisando simulaciones de Monte Carlo" han mostrando
nuevamente la valides de la ecuación de Poisson-Boltzmnnn para bajas concentraciones.
8
El problema ¿e polielectrolitos esférico^'^-^' y la doble c a p s plana, al igual que loa
polielectrolitos cilíndricos y las soluciones electrolíticas han sido tambien
ampliamente estudiados y Poisson-Boltzmann no Lineal, el modelo mas simple que los
describe con cierta precisión, en el límite de las bajas concentraciones.
Todos estos sistemas en gran medida han podido eetudiarse gracias a la simplicidad
de su geometría y al desarrollo de las técnicas nurnéricaa, pero es claro que para
modelos mas realistas, ver por ejemplo Fig(2), la geometría y las ecuaciones resultantes
son mas complicadas. En general, a menudo se requiere de la solución de ecusciones
en dos y tres dimensiones y en este caso métodos numéricos mas eficientes que los
tradicionalmente usados para este tipo de problemas, se hacen necesarios. El m6bodo
de Elemento Finito ha probado ser de gran utilidad en otras areas como hidrodinámice
y electrodinámica. Mostrar la posibilidad de utilixar dicho método en dos dimensiones
para calcular la distribución de carga alrededor de un cilindro es el objetivo
de la^ siguientes secciones.
22-26
26-28
Fig ( 2 ) El modelo del DNA
9
CAPITULO I1
A) El Método de Elemento Finito
El método de elemento finito es un método aproximado de resolver muaciones
diferenciales con condiciones a la frontera y/o iniciales.
En este método un continuo es dividido en elemenbs muy pequenos de formas
convenientes (trihgulos, cuadrilsteros, etc ... ).Escogiendo puntos llamados "nodos",
dentro de los elementos, la variable en la ecuación diferencial se escribe como
una combinación lineal de funciones de interpoiación construidas con el valor en
los nodos de las variables o el de sus distintas derivadas. Usando principios
variacionales o métodos de residuos pesados, las ecuaciones diferenciales gobernantes
mn transformadas en "ecuaciones de elemento finito" que gobiernan todos los elementos
aislados. Estos elementos locales son finalmente recopilados para Formar un sistema
global de ecuaciones diferenciales o algebriicaa con condiciones a la frontera y/o
iniciales impuestas apropiadamente. Finalmente los ralores nodales de las variables
son determinados al resolver ese sistema de ecuaciones.
El método de elemento finito fue originalmente desarrollado por ingenieros alrrededor de los años cincuentas, ya que los problemas en ingenierla frecuentemente presentan
dificultades que no se deben a su complejidad matemática, sino simplemente a causa
del número de componentes individuales que es th presentes (Fig I).
10
\ \ Elementos-,
Algunos Problemas Discretos en Ingeniería y sus elementos.
Fig( 1 1 Para el ingeniero Fue muy natural pensar en un "continuo" en términot3 de un
eneamble discreto de pequeñas componentes. En 1956 Tomer,Clough,Marbin y Topp,
presentaron el primer articulo donde el método de elemento finito ea utiliaado
para analisar grandes sistemas de elementos estructurales en aeronaves.
Aplicaciones de elemento finito a problemas no estructurales, tales como, el flujo
en fluidos y electromagnetismo, empezaron con Zienkiewks y Chung en 1956, asi
como aplicaciones a una amplia clase de problemas de interés en mecánica no lineal
fueron hechas por Wen en 1972.
El método variacional de Rayleigh-Ritdl877-1909) y el rnétdo de Calerkin(l915)
establecieron, principalmente, el método de elemento finito como una importante rama
de aproximación para resolver ecuaciones diferenciales con condiciones a la frontera
y/o iniciales.
11
Cabe aquí mencionar que es &lo recientemente (iW@-i972), que se reconsideran
las ventajas del uso del procedimiento de residuos pesados y en particular el método
de Galerkin en la formulación de problemas mas generales, asi como el desarrollo de
su teoría matemtitica
En este capítulo recordaremos el m W o variacional de Rayleigh-Rits y el método de
Galerkin y los utililisaremos para resolver las ecuaciones de Poisson-Boltsmann lineal
y Poisson-Boltzmann no lineal respectivamente.
12
B) El Cálculo Variacional.
El método de Rayleigh-Ritz está basado en el cálculo variacional que, en su forma más
simple, consiete en determinar una función desconocida y=y(x) para la cual la integral
entre los puntos Po(xo, yo) y Pl(xl , y*) es un mlnimo. Fyx, y, y ) es una
función que tiene derivadas parciales de Begundo orden continurrs y suponemos que
hay una funcibn y=y(x) con derivada continua que minimiga a I, ver pig( 2 ).
1
PI _""""
PO - I 1
Fig( 2 1 Escogemos una familia de curvas para comparar, de la siguiente manera:
q='p7(x) cualquier función t a l que tenga derivada segunda continua y además
~(xo)=(rl(xl)=O. Si a es un parámetro pequeño entonces:
13
representa una familia de curvas que pasan por Po y PI ver Fig( 3 )
I I I
I I I
Sustituyendo ( 2 ) en ( I ) obbenemos una integral que ahora es función de a
Si a = O de ( 2 ) entonces y(x)=y,(x) y para que y(x) minimice la integral ( 1 )
para a = O, la integral I(a) debe tener un minimo. La condición
necesaria para que I(a) tenga un extremo en a = O, como sabe,
BE la anulación de 8u derivada para a4. Entonces bI(a)=(aI/aa)da=O
y desarrollando tenemos que:
Y .
14
como 61(0)=0
Mediante integración por partes obtenemos que:
x 0 *o
Para llegar al resultado anterior se usó el hecho de que q(xo)=rl(xl)=O,
entoncea
o finalmente:
x1
entonces, para cualquier a x ) ,
(aF /ay) - (d/dx)(aF/ay)) = O,
que es la bien conocida ecuac@n de Euler-Lagrange.
.”( 3 )
El principal interés lo encontramos al generalisar lo anterior, y aplicarlo
al problema de minimizar la doble integral; ver Pig( 4 )
""( 4 ) n
// -u =u C X , Y \ ,
r l /@ I l l
/ I -;cr t
Y
r Pig( 4 )
La ecuación de Euler-Lagrange correspondiente puede ser derivada por un pro-
cedimienbo similar al caso anterior-
En este caso la función se considera derivable tres veces, la superficie U=U(x,y)
en la cual se realiza el extremo, se supondrá derivable dos veces.
Consideremos nuevamente la familia monoparamétrica de superficies:
U = U(x,y,a) = U(x,y) + a8U, siendo 6 U = U,[x,y) - U(x,y) que
contiene, cuando a = O , a la superficie U=U(x,y), en la cual se realisa el extremo
y, para a=i, cierka superficie admisible U=U,(x,y).
En términos de l a s funciones de la familia U[x,y,a), la funcional I se trasforma
en una función de a , la .cual debe tener un extremo cuando a = O , por lo tanto
=lJ{ (&,F)(acrU) + (a,F#aaP) + (aQF)(aaQ) }dx dy n
16
donde: P=(& U) y Q=(a,U)
como;
Sustituyendo en ( 5 ) de ( 6 ) y ( 7 ) (&F)bP, (a,F)¿%Q entonces:
61 = { (duFNU + a,(&F)bU + d,(d~P) 8 U }dx dy n
n
.""( 5 )
"-( 6 1
.....I 7 )
La integral segunda de ( 8 ) utilisando el teorema de Green se puede trasformar
en una integral de contorno, que resulk s e r cero, puedo que, en el contorno
la variación 6 U 4 ya que las superficiert admisibles pasan por el mismo mntorno.
Por lo tanto :
61 = SS( (&F) + dl(apF) + a,(a,F) }8U dx dy
De la condición de extremo la integral anterior es cero y como &I es
arbitraria entonces :
(auF) + a,(apF) + d,(aqF) = O ."( 9 )
~8 s. ( 9 ) la extensión en dos dimensiones de la ecuación de Euler-Lagrange,
ec. ( 3 ) para una dimensión.
17
Si suponemos por ejemplo en [ 4 )
donde f[x,y) es una función desconocida, la sustitución de esta función en ( 9 )
nos dará.:
O
-V2U = ~(xJ) ,
que es la ecuación de Poisson.
Si f[x,y) es la distribución de carga canónica,i.e.
i=l
donde Hx,y) es el potencial electrostático promedio, entonces la ecuación de
Poisson se convierte en la ecuación de Poisson-Boltamann: y\
V2p = -(~T/E) pie& EXp(-IpeZ,sp)
En la siguiente sección presentamos el método de Rayleigh-Rib que posteriormenb
utilixaremos para resolver la ecuación de Poisson-Bolbmann lineal, y en Is seccion
D presentamos el método de Galerkin con el que resolveremos la correspondiente ecuación no IineaL
i= I
18
C) El Método de Rayleigh-Rits.
Cualquier medio continuo consiste teóricamente de un número infinito de puntos
en los cuales todas las variables pertinentes estan definidas.
El m6todo de Rayleigh-Ribs es un procedírníenh aproximado por medio del cual tales sistemas continuos son reducidos a sistemas con un nGmero finito de variables y
es una aplicación directa de los principios variacionales discutidos anteriormente-
Consideremos el problema general de encontrar el mínimo de la doble integral
.""( 10 ]
con las condiciones de frontera:
u = ca(s) sobre el contorno r ""( 11 )
donde I' es el contorno que encierra al dominio 0
Consideremos ahora una familia de funciones que dependen de varios par6metros
u = M X , Y, at, ."-an1 "( 12 ) tal que la condición de frontera dada en la ecuación [ 11 ] se satisfaga para
todos los valores de l o s pargmetros.
Sustituyendo [ 12 ) en ( 10 ) convertimos la integral I(U) en una función de n
variables al, +, -"-, 8,.
4U) = 1~~1,a*v".,anl Ya que buscamos el mínimo de estas funciones, los números ak deben satisfacer
el sistema de ecuaciones siguiente:
(&,I) = O para (k=l,"..--n)
19
Resolviendo para estos parámetros desconocidos y sustituyendo en ( 12 ) obtenemos
la solución aproximada requerida.
&,Y) = cp(x, Y, at,.----,an) "-( 13 )
Si consideramos que la solución ( 13 ) es solamente una aproximación a la verdadera
solución U(x,y), psra la cual I(U)=M,donde M es el valor del mínimo, entonces puede
asegurarse que la secuencia de integrales I(*), I(%),... tiende al verdadero
minimo, esto es.
limn- ( I ( 4 ) = I[U) = M
20
D) E\ Método de Galerkin.
El m6todo de Galerkin es parte de los metodos conocidos como mdtodos de residuos
pesados. Donde la idea es obtener una solución aproximada a ecuaciones diferenciales
de la forma:
-Au-f = O en Q -I 14 I con A un operador diferencial, f una función cualquiera y sujeta a condicionen a la frontera:
u[x,yj=u0 sobre el contorno I' La solución aproximada a [ 14 ] es construida introduciendo un conjunto de
funciones de la forma.
.....( 15 )
donde Ci son constantes y qi(x,y) son funciones linealmente independientes que
mtisfacen l a s condiciones a la frontera del problema.
Y a que ( 15 ) es una solucidn aproximada, al sustituirla en ( 14 ) tenemos que no
satisface la ecuación diferencial y entonces.
- b f = e
donde e es un error residual Ahora introducimos un conjunto de funciones de peso Wi (i=l,".m) y construimos
el producto interno (e,Wi). Además, pedimos que dicho producto sea cero, esto es.
(e,Wi) = O
Lo cual ei equivalente a forzar que el error en promedio de la ecuación diferencial
aproximada sea cero. Existen varias formas de escoger las funciones de peso Wi; llevando esto a
distintos métodos: 1) Método de Calerkin, 2) Método de Mínimos Cuadrados
3) Método de Momentos, 4) Método de Colocación.
En el m M o de Galerkin las Funciones de peso Wi s o n el conjunto de Funciones
con las que se construye la solución aproximada u.
Entonces:
letpi) = 0
Que podemos entender como una proyeccidn ortogonal de l o s residuos al conjunto
de funciones linealmente independientes pi-
Hacer el producto interno igual a cero es hacer cero la siguiente integral
16 )
La expresidn ( 16 ) servir5 para determinar las constantes Ci y al sustituirlas
en la expresión para u(x,y) obtendremos la solución aproximada que requerimos.
Antes de proceder a aplicar los dos métodos antes descritos a nuestro problema,
en la siguiente sección describiremos la función de interpelación y l o s elementos
finitos que Fueron utilizados.
22
E) Funciones de Interpolación y Elementos Finitos.
El siguiente paso en la técnica de Elemento Finito es escoger adecuadamente las
funciones de interpolación del problema planteado.
Dichas interpolaciones de elemento finito están caracterizadas por la Forma del
dominio local para el elemento finito (dicho dominio local tambien es llamado
elemento Finito) y el orden de aproximación. En general la elección del elemento
finito depende de la geometría del dominio global, el grado de precisión deseado,
la facilidad de integración sobre el dominio, etc. En la Fig( 5 ) mostrarnos un dominio
en dos dimensiones que es discretkdo por elementos triangulares, donde como ya
mencionamos el dominio global est6 construido por muchos subdominios llamados
elementos finitos. Los dominios pueden ser uni, bi, ó tri-dimensionales y las
correspondientes geometrías de elementos finitos son mosbradas en la Fig ( 6 1 Cada esquina del elemento recibe el nombre de nodos locales. En general, l a s funciones
de interpolación son polinomios de diferentes grados, sin excluir, en algunos casos
productos de polinomios con funciones triionométricas o exponenciales.
En seguida describiremos la Función de inbrpolación, consistente en un polinomio
de grado uno en x y y, definidos sobre elementos triangulares ya que fueron los
utiliaados en es& trabajo.
i) Elementos Triangulares.
Las propiedades del elemento son determinadas en tkrminos de las coordenadas
oarteoianaia looslea oon BU centro en el origen del centroide del trisngufo Fig ( 7 ).
Y '
o x
Fig( 7 I Considerando como función de interpolación
s a + bx + cy .....( 17 )
Que representa una variación lineal de w tanto en la dirección x como en y dentro
del elemento triangular.
Para evaluar a, b, c construimos 3 ecuaciones en términos de valores conocidos
de (I en cada uno de los tres nodos locales.
Sean estos valores q, ccp, (18 entonces:
= a + bxl + cy1
% = a + bx2 + cy2
u3 = a + bxS + cy3
24
Resolviendo para 18s 3 constantes tenemos:
cl = ( l/det)(x3 - x*)
c2 = (l/det)(x, - xg)
cs = (ljdetj(x2 - x,)
donde:
1 x1 Y1
det = 1 x2 ya
1 x3 Y3
=2 A ; A= &rea del triangulo
EB importante notar que la numeración (1,2,3) de los nodos en el triSngulo son
asignadas en el sentido cpntrario a las manecillas del reloj como se muestra
en la Fig( 7 ). Si esta numeración es asignada en el sentido de las manecillas
del reloj el determinate resulta igual a - !?(área del triángulo).
Uno de los requisitos Fundamentales de la función de interpolación es:
25
Equivalente a que los parámetm oWescan las restricciones; al + a2 + a3 = 1
bl + b2 + b3 = O
c1 + c2 + c3 = o
El siguiente paso en la técnica de elemento Finito consiste en construir el dominio
del problema con estos elementos finitos y numerar globalmente los vértices de
todos los elementos Finitos, ver por ejemplo Fig( 8 ), esta numeración corresponderá
a la numeración de los nodos globales donde se encontrar& la solución y de la forma
que se haga dicha numeración dependerá la estructura del sistema de ecuaciones
algebriiicas que será necesario resolver,
Fig( 8 1 Con esto concluimos la presentsción de la técnica de elemento finito que ahora
aplicaremos en las siguientes secciones para nuestro problema.
26
F) Solución a la Ecuación de Poisson-Boltsmann Lineal
Esta ecuación para geometría cilíndrica de la cual m conoce su solución analítica
aeri resuelta numéricamente utilitando el método de Rayleigh-Rits en el plano
cartesiano. Esta ecuación se escribe de la siguiente manerrr:
V'$ = K2*
donde k es el número de especies en el electrolib.
Si consideramos un electrolito simétrico (electrolib 1:i)
p=pi=p* ,& =-z* = 1
entonces. K;' = [4ne'P/&) (plZ: + p&', = [8.1re2plp/&)
p = densidad numérica
e = carga elemental
= constante Boltxmann
c = constante dieléctrica
Antes de continuar construimos la ecuacion dimensional llamando:
1L* = eS1L
X* = x/[a/2)
Y* = ~ / ( a / 2 )
donde [a/2) = radio de los iones.
De esta manera obtenemos la eouacion adimensional :
(9) 7k = & I 1cI 2 * * 2 * * donde [K?* = (dXa/Z)"
Para simplificar la notacih omitiremos el símbolo * quedando establecido
que el problema 6 sin dimesiones.
37
Para resolver V'I) - n2q = o .....( 18 )
con el fin de comparar la solución numérica con la analitica, ya que esta puede conocerse
y además usarla como solución de entrada pars resolver la ecuación no lineal
utilizaremos el método variacional de Rayleigh-Rits donde la funcional ser6 ;
1 = J J (m{ (ax$)' + (a,W' + K 2 V 1 dxdr ..-I 19a ) n.
que satisface ( 18 ) y las siguienk condiciones a la fronbra sobre el contorno I' que
donde KO y Kt son funciones de B e s s e l de orden O y 1 respectivamente.
Ahora proponemos como solución aproximada la iguiente:
Con c p i (x,y) polinomios lineales como los discutidos en la sección anterior
20
Sustituyendo( 20 ) en I y en virtud del método variacional aI/an = O tenemos:
."..I 21 )
De todas las y, S610 son desconocidas desde i=l, .. .m, y de i=n+i,..m
son dadae como condiciones de frontera. La relación ( 21 ) es un sistema de n ecuaciones
lineales con n i n & n i t , a s que escribimos de la forma.
& Y = / ? .....( 22 )
donde los elementos de la matris .d como los del recbr son:
Como el dominio Q fue discretido en elementos f i n i b triangulares, l a s
integrales anteriores deben realimse sobre los triángulos que camparkan los
nodos i, j, con i=l,-.-n j=l,"n que es donde queremos resolver, para los
elementos . a,k dichas integrales deberh realizarse sobre aquellm
triángulos que compartan los nodos i, k con i=l,,n k=n+¶,.,,m donde los
nodos k estin sobre el contornor que encierra al dominio Q y en ellos se
está imponiendo las condiciones de fronfera
29
Por otro lado cp,, lp,, v k corresponderán a una de l a s tres funciones de interpolación definidas en el triingulo, sobre el que se está
realizando la integración y que coincide con el nodo global i-esimo.
Por lo tanto en el entendido de que, al construir a,, O Qt,k debemos Sumar
sobre todas las contribuciones de aquellos triángulos que compartan los nodos i,j
o i,k
Sustituyendo en forma explícita la función de interpolación en IOS elemenbs
de matris ati i , para el triángulo I-ésimo,tenemoa que:
O(ii = JJ bib, dxdy + Sf cici dxdy + SSK'vivj dxdy 7,. 9 7;,
= bib, A + C ~ C , A + K'JJ Vivi dxdx 7 9
donde A= área del triángulo 1-ésimo.
Encontrando todos los elementos del tistema de ecuaciones [ 22 ), resolvemos este
para calcular el vector 7 que es directamente la solución en l o s nodos.
30
G) Solución a la Ecuación de Poisson-Boltsmann Resolveremos numericamente esta ecuacidn utilizando el metodo de Galerkin.
Sea:
"".( 28 )
donde j es el número de especies iónicas en el electrolib.
De igual manera al caso lineal el electrolib se considerará simétrico y la ecuación
se adimensionalim, obteniendo: z
Vzl(/ =-E CF, EXP(-Z#] i: 1
donde:
E = (4ne2b/&)(a/2)*
Proponemos como solución aproximada: VA
U = yi~p, con ~ p , = a,+biX+ciy t"1
donde i=número de nodos globales, y de n+l h a s t a m corresponden a l a s condiciones a la frontera.
Utilisando el m6tudo de Galerkin tenemos que:
-{ 1 J { V 2 ~ } ~ i dxdy + $$ EzP,EXP(-Zi UFpi dxdy = O } 2
n n j:t ".( 24 1
con i=l,"-,m.
Como las funciones de inkrpolación lineales que usaremos no tienen segundas derivadas, utilisamos el teorema .de Green-Gauss sobre el primer integrando para
pasar a primeras derivadas, entonces:
31
Ya que sobre los conbrnos no existe y,, puesto que sobre ellos damos
las condiciones de frontera, la integral de conbrno es cero.
Entonces:
Y sustituyendo w tenemos;
El segundo integrando ai r e a l i z a r s e sobre cada elemento triangulsr, de la suma sobre el
indice k sólo contribuyen aquellos correspondientes a los 3 vértices del triángulo
sobre el que ge está llevando a cabo la inhgración.
Entonces escribimos para el triángulo 1-ésirno:
La relación ankrior corresponde a un conjunto de ecuaciones algeb6icas no
lineales que podemos escribir de la siguiente manera:
32
7 que es directamente la solución en los nodos.
33
CAPITULO 111
RESULTADOS
En esta sección presentamos loa resultadoe obtenido8 al resolver las ecuaciones
de Poisson-Boltamann no lineal y Poisson-Boltxmann lineal utilizando el método
de Elemento Finito descrito en el capitulo anterior para un cilindro inmerso en
una solución electrolitica 1:l. Los parámetros del sistema son IOS siguientes:
Radio del cilindro 5 w
Temperatura 298 K O
Constante dieléctrica 78.5
Diámetro de los iones 425 A
En la TABLA 1 mostramos los parhetros de l a a corridw de l a s cuales preaentamos
los resultados.
TABLA 1
CORRIDA DENSIDAD MOLAR POTENCIAL DEL CELINDRO(mr)
A 1
B 1 C 1
D - .o1
40
80
120
m
La malla de elemento Finito fue construida con 704 nodos y 1344 triángulos como
mostramos en la PIC( 1 1
34
. .-
FIG( 1 1 Los 704 nodos fueron distribuidos en 21 círculos (equipotenciales) concéntricos
numerados del O al 21 y sobre cada círculo se colocaron 32 nodos uniformemente
espaciados. En el círculo O se da la condición de frontera para el potencial H, que para el caso de Poisson-Boltrsmann lineal puede calcularse utilizando la relación
(19b) del capitulo I1 y para el caso de Poisson-Boltmann fue tomado de lw resultados
obtenidos con el método de iteración de Picard (E.GonsSles y MLoxada)
De el círculo 1 al 20 hay 640 nodos donde Be encontrará la eoluoi6n eproximada
de elemento finito. En el círculo 21 se da la condición de frontera para el potencial
cp(r -+ m) = O.
En la Fig( 2 ) se muestra un quema de lo anterior;
. ""
-
35
El criterio utilisado para
DENSIDAD MOLAR
1
.o1
hacer el potencial cero ea mosrtrado en la TABLA 2
TABLA 2
RADIO DE CORTE ( A ) 1pp~1* (PPB*
42.5 O(E-06) O@-07)
212.5 qE-03) qE-03)
La distribución de los 20 circulos no fue uniforme dependiendo del radio de corte
y la8 tablas 3 y 4 muestran esfa distribución de los radios reducidos en el radio
del ión para l a s densidades 1M y .OlM respectivamenbe.
TABLA 3
r* 1.5 2.0 25 3.1 3.7 43 5.0 5.7 6.5 7.3
8.2 9.1 10.1 111 12.2 13.3 145 15.8 17.1 18.5
TABLA 4
r* 15 2.1 2.8 3.7 4.7 6.0 7.6 9.4 116 14.3
17.4 21.0 25.4 30.5 36.4 43.4 51.5 61.1 72.1 85.0
En las gráficas, 1 a 4, presentamos el error relativo porcentual de la solucih
a la ecuación de Poisson Boltsmann lineal con elemento finito, respecto a la
solución analitica.
En las gráficas, 5 a 8, presentamos la solución con elemento finito para la
ecuación de Poisson-Boltamann no lineal comparandola con la obtenida con Picard, y su error relativo porcentual en las grhficas 9 a 12
36
Finalmente en la TABLA 7 presentamos el cálculo de la densidad de carga en la
superficie del cilindro calculada con las Funciones de distribución radial obtenida
con el potencial calculado con elemento finito para la ecuación de Poisson-Boltzmann
no lineal y comparada con los resultados obtenidos usando el método de Picard.
Este cálculo de la densidad puede obtenerse de la siguiente manera:
Como sabemos, la carga sobre la superficie del cilindro debe ser igual a la carga
en el volumen que lo rodea
Entonces:
J ads = - J Pelectrolib dv
de manera que:
entonces: f-¿
Q = (-Zpe/R) (g+ - g-) rdr 9 W Z
utilizando las unidade; reducidas r* = [2r/a)obtenemos finalrnenk
cr = (-Zpe/R)(a/2)* (g+ - g-) r*dr* b
Q donde :
a = R * + 1 y b = R * + %
donde: .% es el radio de corte dado en la TABLA 2
La integral (1) fué evaluada utilizando integración de Simpson.
37
TABLA 7
CORRIDA v ELEMENTO FINITO(scoul/cm2) Q PICARD(scoul/cm2) A .2380411305 -2367674305 B -4893716305 .4861100305
C "1633972E05 D -6181699305
-7564360305 -6156324305
-.I
Ir,
Error relativo‘% para la ecuación de PBL en g ( r ) obtenida con elemento finito respecto a g ( r ) analítica. Para un
. cilindro de radio 5 A y potencial sobre su superficie de $3 mv, electrolito 1:l y densidad 1M.
Lo mismo que en la gráfica 1 .Para 80 mv
4 0
Lo mismo que eh l a gráfica 1. Para 120 mv
4 1
CWlPIcA 4 COMIM D
. .
Lo mismo que en la gráfica 1. Para 200 mv y densidad .01M
42
P
A i - PICRRD I
!
!
43
4.
3.
t.
I_ Pi CARD
-0 FINITO
1.
O O 1 il t d lb ir S lb 1
” Lo mismo que en la gráfica 5. Para 80 mv,densidad 1M
i
4 4
ELECLENTO FINITO
4r b 8 da lh 14 16 18 ; w
Lo mismo que e; la gráfica 5. Pero 120 mv, densidad 1M
4 5
1
i Fw
Lo mismo que en la gráfica 5. Para g+ ( r ) , 200 mv y densidad .01M
4 6
4 7
* Error relativo-% para la ecuación de PB en g(r) obtenida con elemento finito respecto a g(r) obtenida con Picard. Para los parámetros de la gráfica 1.
4 8
M Lo mismo que en la gráfica 9. Para 80 mv densidad 1M
4 9
r i I
s.,.
b .
Lo mismo que en la gráfica 9 . Para 120 mv densidad 1M
5 0
.
.
I i I I I I I I I I 19
Pf Lo mismo que.en la gráfica 9. Para 200mv densidad .01M
5 1
DISCUSION DE RESULTADOS Y CONCLUSIONES
Para obtener la precisión que mostramos en los resultados anteriores en la Función
de distribución radial, fue necesario construir una malla con 704 nodos. Esto significó resolver un sistema de ecuaciones algebriiicas (no lineales en el caso
de Poisson-Boltsmann y lineales en el caso de Poisson-Boltsmann Lineal) de
640 x 640. Dicho sistema de ecuaciones fue resuelto utilizando las subrutinas
MA28A y MANA de la bibliobka WARWELL, que son subrutinas muy eficientes para
resolver sistemas de ecuaciones ralas, a pesar de esto, el sistema de 640x640
fue el mas grande que se pudo resolver en una CYBER-173 que es el sistema de
cómputo donde estos cilculos se realizaron.
Sabiendo eska limitación, se procedió a encontrar la distribución mas adecuada
de esos 704 nodos. Para ello utilisamos la solución a la ecuación de Poisson-
Boltsmann como calibrador ya que conocemos su solución analitica.
La primera distribución que intentamos consistió en mantener u n 3 distribución
de equipotenciales igualmente espaciadas y variar el número de estas así como
el número de nodos en cada equipotencial Esta manera de distribuir los nodos
sólo recuperaba, en todos los casos, la primera cifra de la soiución analítica.
Fue claro que la poca precisón es debida al hecho de que estas mallas resultan
poco refinadas para cubrir el dominio con 8610 704 nodos.
Una de las ventajas del método de Elemento Finito frente al de Diferencias Finitas, es que puede variarse el grano de la malla sin el compromiso de hacerlo en todo
el dominio de integración, de manera que esto permite refinar la malla donde se quiera mejorar la precisón.
La siguiente distribución tomó en cuenta lo anterior y se decidió sacrificar la
52
precisión lejos del cilindro para mejorarla cerca. Vrtriando de nuevo el nilmcro
de equipotenciales y el número de nodos sobre ellas pero distribuyendo los radios
de manera exponencial para llegar rapidamente al radio de corte.
Nuestros resultados finalmente pudieron mejorarse hasta recuperarse la tercera cifra,
separando los radios de las equipotenciales de acuerdo a como se muestra en las
gráficas 13 y 14 con 21 equipotenciales y 32 nodos en cada una.
Con esta distribución de nodos, como mostramos en los resultados, pudimos resolver
e¡ problema en 2 dimensiones con una muy buena precisibn comparada con Is solución
obtenida en 1 dimensión. En este momento es necesario mencionar que, la solución
contra la que estarnos comparando, es una solución obtenida dividiendo el dominio
unidimensional en 131 partes y si comparamw esto con las 21 equipotenciales
utilixadas, nueertra partición es bastanke menor y m obtienen muy buenos resultados.
Otra ventaja al utiligar elemento finita es que podernos obtener la solución, no
sólo en los nodos sino también en todo el continuo mediante los polinomios de
interpolación característicos del método de elemento finito. Y pudo comprobarse
en el cálculo de la densidad de carga, (TABLA 7) que la inberpolación con polinomioa
lineales, pars este problema, resulta una buena aproximación.
Ahora bien, si recordamos, la solución a la ecuación de Poisson-Boltsmann m obtiene resolviendo un sistema de ecuaciones algebráicas no lineales. En este trabajo
utilimrnoa, para resolver dicho sistema, el método de Newton-Raphson y como función
de prueba se utilizó la respectiva solución obtenida con elemento Finito a la ecuación
de Poisson-Boltsamann linegl. kta solución resulb una muy buena aproximación
a, la solución no lineal, ya que en l a s cuatro corridas que realixamos en promedio
bastaron 3 iteraciones para que la diferencia entre la solución de entrada y la de
salida fuera menor a IE-O4 , criterio que utilissmoe para la convergencia a la solución.
1
Hasta este momento hemos hablado exclusivamente de la precisih de l o s resultados
53
obbnidos con elemento finito en 2 dimenaionea, pero ea deseabíe de toda técnica
numérica que el tiempo de CPU que requiera sea el menor posible. Decimos esto porque mientras una corrida utilizando el método de iteración de Picard en 1
dimensión el tiempo de CPU es del orden de 10 minutos(CYBER-l73), utilizando
elemento finito en dos dimensiones con 704 nodos una corrida ocupa aproximadamente
2 horas, que en el caso de Poisson-Boltmnann ese tiempo debe multiplicarse por
el número de iteraciones. En estos momentos, tal y como tenemos los programas, de
esa8 2 horas el 80% ea ocupado en wndruir el sistema de ecuaciones dgebriiccrs.
Es .entonces la dimensionalidad del sistema de ecuaciones la principal dificultad
encontrada en la tknica de elemenb finito en 2 dimensiones. Es claro entonces
que la conotrucción de la malla en elemento finito debe hacerse utilizando al miximo
la simetria del problema. Esto reduce el niímero de nodos y por lo tanto el tiempo
de CPU.
El error numérico sobre la equipotencial está determinado por el número de nodos
sobre ella. Por otra parte, en nuestro m a debido a l a s limitaciones de memoria
el error numérico introducia una asimetría ficticia en l o s cálculos. Por este motivo,
para hacer mas eficiente la convergencia del programa se tomó sobre los nodos de
la equipotencial el valor promedio de cp-
Otra manera de disminuir el número de nodos es utilizar la simetría del problema.
Y resolver s610 un cuadranb del plano cartesiano Fig( 3 1 Sin embargo esto será
objeto de trabajo posterior- Para esto es claro que hay que modificar las condiciones a la frontera del problema, ya que es: necesario agregar condiciones a la frontera
sobre los nuevos contornos que forman :los ejes Y y X, estas condiciones pueden ser
sobre la derivada normal ya que p o r la simetría sabemos que la derivada normal en
e m 8 fronteras debe ser cero.
54
Y
Fig (33
Si utilizamos lo anterior es necesario cambiar las Funciones de interpolación
ya que las funciones l ineales que aquí utilizamos no satisfacen dichas condiciones
de frontera.
Resumiendo, es factible aplicar el método de elemento finito utilizando polinómios
de interpolación lineales y elementos triangulares para resolver las ecuaciones de
Poisson-Boltzamann y Poisson-Boltamann Lineal en dos dimensiones con una muy buena
precisión. Debido a la gran dimensionalidald de los sistemas algebráicos que es necesario
resolver, se requiere de importantes recursos de cómputo.
En futuras aplicaciones del método, es imporbante bomar en cuenta que para h e r
buena precisión, es necexário graduar la malla de manera adecuada al problema que
se quiera resolver. Hemos airendido también, que para la precisión es muy importante
que la densidad de nodos en la dirección angular sea del mismo orden que en la
dirección radiaL
Para hacer más eficiente la aplicación de elemento finito ser& necesario buscar
en el problema a resolver algun plano de simetría, que haga disminuir el dominio
de integración, recordando que entonces es necesario utilizar otros polinomios de interpolación diferentes a lo que aquí empleamos.
55
Finalmente es importante remarcar que utilislamos el m6todo de Elemento Finito
en su forma m& eencilla y creerno0 necesario explorar Gcnicaa mas refinad= para
hacer eficiente su aplicacidn a estos proble!mas.
. .
. .
5 7
iae !
.
58
REFERENCIAS.
1.- P-Debye and E.Wuckel, PhysikZ. !14,185 [1923)
2.- D.A.McQuarrie, Statistical Mechanics (Harper,Row) Cap 15.
3.- R.H.Fowler, Statistical Mechanics (Cambridge University Press.) Cap 8.
4- LOnsager, Chem Rev. 13,73(1933)
5.- J.C.Kirkwood, J. Chem. Phys 2,76711934)
6.- J.G.Kirkwood and J.C.Poiriere, J. phys. Chem. 58,591(1954)
7.- D.Henderson, J.A. Barker and M.Lozada-Cassou, J. Chem. Soc.
Faraday Trans.:! 81,45711985).
8.- J.Frahm and S-Diekmann, J. Colloid Interface Sci. 70,440(1979)
9.- Y.Gor, 1.Ravina and A.J.Badchin, J. Colloid Interface Sei. 64,33311977)
10.- P.T.Wal1 and J.Berkowibs, J. Chenn. Phys. 26,11q1957)
11.- G-V-Ramanathan, J. Chem. Phys. 78,3223[1983)
12- J.A.Schellman and D.Stigter, Biopolymers 16,1415(1977)
13.- S.L.Brenner and D-A-McQuarrie, ,J. Theor. Biol. 39,343(1973)
14.- D.Stigter, J. Colloid Interface Sci. 53,296(1975)
15.- DStigter, J. Phys. Chem 82,1603(1~978)
16.- M.Losada-Cassou, J. Phys. Chem. 87,3729(1983)
17.- EGonsáles, MLomda-Cassou and D-Henderson, J. Chem Phys. 83,36111985)
18.- P.Milla, Ch.F.Anderson and M.T.:Record Jr.,J. Phye. Chern. 89,3984(1985)
19.- Verwey, EJ.W. and Overbeek, J. Th C., Theory of the Stability of Lyophobic Colloids", Elsevier Putdishing Company, Inc., N.Y., (1948),
Caps. I1 Y 111.
20.- G.N.Patey, J. Chem. Phys. 72,5763(1980)
21.- M.Losada-Cassou and RSaavedra "Resulbados Preliminares Sobre un Modelo
para una Proteina", Reporte Interno UAMI(1980)
59
22.- M-Lozada-Cassou, R.Saavedra and D-Henderson, J. Chem. Phys. 77,515011982)
23.- M-Losada-Cassou and D. Henderson, J. Phys. Chem. 87,2821(1982)
24.- Ch.W.Outhwaite and L-B-Bhuiyan, J. Chem. Soc. Faraday Trans.:!
79,707(1983)
25.-
26.-
27.-
28.-
T.L.Croxton and D.A.McQuarrie, IMol. Phys. 42,141(1981)
D-Henderson and M-Lozada-Cassou, J. Chem. Phys. 79,3055(1983)
CRassaiah and LFriedmann, J. CIhem Phys. 48,2742f193)
T.L.Croxton and D.A.McQuarrie, ,J. Phys. Chem. 83,184of1979)
60
BIBLIOGRAFIA
T.J.CHUNC, "Finite Elements Analysis In Fluid Dynamics"
Ed. Mcgraw-Hill [ 1978)
BRUCE A. FINLAYSON, 'The Metshod of Weighted Residuals and Variational
Principles"
Ed. Academic Press (1972)
BRUCE A. FINLAYSON, "Nonlinea,r Anaiyais in Chemioal Engineering"
E" McGraw-Hill 11980)
J-R-WHITEMAN, The Mathematics of Finite Elemenbs and Applications"
M. Academic Press (1973)
L-ELSGOLTZ, "Ecuaciones Diferenciarles y Cálculo Variacional"
Ed. Ediciones de Cultura Popular (19'75)
"LOZADA-CASSOU "Satistical Mechanics of Model Charge Systems"
Reporte Interno UAM (1980)
61
Top Related