Esquema de descripción geométrica del crecimiento de grietas usando una aproximación geométrica basada en Level Sets y Fast

Marching Method

Geometric description scheme for crack propagation by a geometric description based on Level Sets and Fast Marching Method

Vicente Francisco González- Albuixech ¹, Eugenio Giner², José Tarancón³

¹ Departamento de Ingenieria Mecánica. Universidad Carlos III de Madrid. Legánes. España. Email: vicgonza@ing.uc3m.es

²Centro de Investigación en Ingeniería Mecánica. Departamento de Ingeniería Mecánica y Materiales. Universidad Politecnica de Valencia. Valencia. España. España. Email: eginerm@mcm.upv.es

³ 2Centro de Investigación en Ingeniería Mecánica. Departamento de Ingeniería Mecánica y Materiales. Universidad Politecnica

de Valencia.. Valencia. España. Email: jetaranc@mcm.upv.es

RESUMEN

El level set method es una técnica que se utiliza habitualmente para describir matemáticamente la propagación de superficies mediante la aproximación y resolución de ecuaciones diferenciales. Esta técnica ha sido ampliamente utilizada en combinación con el XFEM, especialmente en el análisis del crecimento de grietas. Sin embargo, existen otras opciones que se basan en consideraciones geometricas como el fast marching method. En este trabajo se implementa una técnica geométrica basada en el fast marching method y level sets para la descripción de propagación de grietas en XFEM y, por tanto, compatible para su uso directo con esta técnica en cálculos de mecánica de la fractura.  

PALABRAS CLAVE: LSM, FMM, XFEM, Mecánica de la fractura.

ABSTRACT

The Level Set Method is used to describe mathematically the propagation of surfaces, by the resolution and approximation of a set of differential equations. This technique has been widely used on the XFEM framework, specially for crack propagation. However, there are other techniques based on geometrical considerations, as the Fast Marching Method. In this study, it is proposed a geometrical related technique based on the Fast Marching Method and Level Set Method for crack growth description in XFEM. Therefore, the results could be used directly on fracture mechanics analysis using XFEM.

KEYWORDS: LSM, FMM, XFEM, Fracture mechanics.

1. INTRODUCCIÓN

Para simular el crecimiento de la fisura se usan técnicas numéricas de propagación de superficies, como son el Level Set Method –LSM– y el Fast Marching Method – FMM– [1-5]. La descripción básica de la aplicación de dichas técnicas a la descripción del crecimiento de grietas ha sido introducida en [6-10].

Sin embargo, la aplicación de la técnica del FMM a la descripción de la grieta es relativamente escasa [11-13]. En [11] y [12] se introduce la metodología FMM con elementos finitos extendidos –XFEM– pero se aplica a la Este artículo puede compartirse bajo la licencia CC BY-ND 4.0 y se referencia usando el siguiente formato V.F. González-Albuixech , E. Giner ,

J.E. Tarancón, “Esquema de descripción geométrica del crecimiento de grietas usando una aproximación geométrica basada en level sets y fast marching method),” UIS Ingenierías, vol. 16, no. 1, pp. 47-58, Enero-Junio 2017.            

descripción de curvas cerradas o inclusiones en vez de curvas abiertas, y no usan para ello dos funciones distancia, con lo que no se corresponde con la descripción que se realiza habitualmente en el estudio de mecánica de la fractura mediante XFEM. En [13] la implementación es ya más parecida a la usual en crecimiento de grietas, pero aun así existen diferencias significativas.En dichos trabajos, en la resolución del FMM se usa como espacio de aproximación y cálculo el método de las diferencias finitas. Sin embargo, se proporciona información sobre la aplicación de dicha técnica en mallas no estructuradas, como la que se puede encontrar en [14-17], que permite su extensión a otras discretizaciones como las usadas en el método de elementos finitos. 

Una de las principales ventajas del uso del FMM, en comparación con el LSM, es su mayor rapidez y que permite de forma natural una descripción tipo narrow band (descripción local alrededor de la grieta).

El objetivo de este trabajo es establecer un algoritmo que permita usar FMM en la descripción del crecimiento de la grieta. Para ello se utiliza el algoritmo descrito para LSM en [13, 17], aprovechando la implementación en elementos finitos que se puede encontrar en [5]. 

La novedad del trabajo es incorporar el FMM en esquemas ya disponibles en LSM utilizando un espacio de elementos finitos, de forma que la descripción narrow band sea natural De esta manera no hay que resolver la ecuación del LSM, optimizando así la descripción geométrica de crecimiento de la grieta para realizar cálculos de fractura en XFEM.

2. DESCRIPCIÓN DE LA GRIETA MEDIANTE FUNCIONES DISTANCIA

Para describir la posición de cualquier punto del espacio respecto al frente de grieta –e identificar, simultáneamente, la posición y geometría de la grieta en el dominio– se pueden usar dos funciones distancia con signo (level set o LS) [6-8].

Estas dos funciones distancia se suelen designar mediante  y . Sea(x, )t el LS que mide la distancia a la superficie de grieta, proporcionando (x, )t = 0 la posición de la superficie de grieta, y sea (x, )t el LS que mide la distancia a una superficie que contiene el frente de grieta y es ortogonal a la misma. La intersección de las superficies (x, )t = 0 y (x, )t = 0 da la posición del frente de la grieta en el instante t, Figura1.

Figura 1. Funciones distancia en la grieta, [9]. Fuente.

Elaboración modificada de [9].

En resumen, los LS (x, )t y (x, )t , distancias con signo, definen la posición y geometría de la grieta según las siguientes relaciones:

(x, )t      = 0,(x, )t < 0  Caras de grieta

(x, )t      = 0,(x, )t = 0  Frente de grieta (1)

(x, )t > 0           Zona sin agrietar

Además, las dos funciones (x, )t y (x, )t se construyen de forma que cumplan la siguiente relación de ortogonalidad:

(x, )t (x, )t  0       t (2)

Esta ortogonalidad de los dos LS, provee de un camino para definir un sistema ortogonal de coordenadas curvilíneas propio a la grieta, Figura 2. Los vectores que componen la base de dicho sistema de coordenadas se definen a continuación.

Figura 2. Vectores de la base propia a la grieta, [9]. Fuente.

Elaboración modificada de [9].

Un vector normal al frente de grieta:

n_ ^ (3)



Un vector normal a la superficie de grieta:

n_ ^ (4)



Un vector tangente al frente de grieta, que se construye mediante:

e n n_t    (5)

Usando esta definición para la base del sistema de coordenadas curvilíneas propio de la grieta, y por comparación con los vectores de la base euclídea, se pueden definir dos funciones equivalentes a las coordenadas polares r y θ en función de  y , según:

r   2  2

^tan₁  _ (6)  



Estas coordenadas son las que se necesitan para definir las funciones de enriquecimiento asociadas a la grieta en XFEM [6-10].

Además, si se tiene en cuenta la velocidad de crecimiento de la grieta, ésta se puede dividir en dos componentes perpendiculares a las superficies descritas por las funciones distancia:

vv   n vn      (7)

FAST MARCHING METHOD

Una de las técnicas numéricas que se usan para la descripción del movimiento de superficies es el FMM [14]. Esta técnica resuelve la evolución de un sistema en movimiento descrito mediante una función distancia. Para una función distancia d, esta aproximación se traduce en las siguientes expresiones:

d x y( , ) F x( )      en , F x y( , )  0

(8) d  g x y( , )     en 

donde F(x,y) es la velocidad de propagación de la función d(x,y) que evaluamos en el dominio Ω. De d se conocen los valores iniciales g(x,y) en un contorno o región Γ del dominio.

La técnica intuitivamente se basa en suponer que la velocidad de modificación no cambia de signo, y así lo único que hay que hacer es actualizar la distancia según una cantidad que dependerá la velocidad en el punto; una vez que se ha actualizado se verifica si la cantidad es correcta o se prueba otra distancia. Es una técnica sencilla y potente, pero con limitaciones evidentes para muchas aplicaciones.

Detalles de la implementación de la técnica general en dominios triangulados se pueden consultar en [14-16].

3. IMPLEMENTACIÓN DEL FAST MARCHING METHOD

USANDO ELEMENTOS FINITOS PARA CRECIMIENTO

DE GRIETAS

4.1. Descripción

Si se usa un espacio de elementos finitos para la descripción de los LS, entonces:

(x t, ) N_I (x)_I ( )t

I

(9)

(x t, ) N_I (x)_I ( )t

I

donde N x_I ( ) son las funciones de forma de elementos finitos. Igualmente, para las velocidades de los LS:

v x t_(       , ) N_I (x v) _I ( )t

I (10)

v_(x t, ) N_I (x v) _I ( )t

I  

Para los gradientes se obtiene:

(x t, )   N_I (x)_I ( )t

I

(11)

(x t, )   N_I (x)_I ( )t

I

4.2.     Reconstrucción de una función distancia, reinicialización

La idea es resolver la siguiente ecuación eikonal en el elemento mediante relaciones algebraicas, tal como se describe en [17]:

d x y( , )1 (12)

Esta relación describe una función distancia euclídea y coincide con el proceso de reinicialización, es decir, sirve para reconstruir la función distancia.

Si los elementos finitos utilizados corresponden a triángulos o tetraedros lineales, las derivadas son constantes en el elemento y, por tanto, es sencillo encontrar un valor nodal de la función distancia a partir de los otros valores nodales.

Por simplicidad, se va a desarrollar para el caso del triángulo lineal:

d x y( , )_x² d x y( , )_y² 1²          (13)

Si los valores de la función distancia en los nodos son t₁, t₂ y t₃, se tiene:

d x y( , )x  N t1,x 1  N2,xt2  N t3,x 3

(14)

d x y( , )y  N t1,y 1  N2,yt2  N3,yt3

Si t₁ y t₂ son conocidos, se pueden agrupar sus contribuciones en:

tx  N t1,x 1  N2,xt2

(15) ty  N t1,y 1  N2,yt2

Sustituyendo:

tx N t_{3,x 3 }²  ty N t_{3,y 3} ² 1 (16)

Ecuación que se desarrolla hasta:

N3,x2  N3,y2 t32 2tx N      3,x  2 ty N     3,yt3 tx2 ty2 1 0 

(17)  

Agrupando términos:

a  N3,x2  N3,y2  b  2tx N _3,x  2ty N _3,y  (18) c  tx² ty² 1

Se consigue, por tanto, la siguiente ecuación de segundo grado:

a t ₃²  b t₃ c 0 (19)

La distancia en el tercer nodo puede tomarse como la máxima de las dos raíces reales obtenidas, es decir:

^{ }b        b²  4ac  b    b²  4ac

t₃  max(   ,       )

2a     2a

(20)  

Hasta aquí se ha considerado que las funciones distancia son siempre positivas. Para considerar funciones distancia con signo, se trabaja con el valor absoluto de los valores nodales y se introduce el signo al final del cálculo. De esta forma, si el signo que tiene el valor nodal t₃ antes del proceso se denota por signo(t’₃), entonces el t₃ final es:

t₃  signo t( ' ) max₃             (^{ }b b²  4ac ,  b  b²  4ac )

2a     2a

(21)

4.3. Resolución de las ecuaciones de ortogonalidad y de extensión ortogonal

La misma idea de que las derivadas de las funciones de forma en un elemento triangular lineal o un tetraedro lineal son constantes, puede aplicarse al resolver otras ecuaciones que son necesarias para modelar la propagación de grieta. Por ejemplo, conocidos los valores del LS  en todo el dominio, se puede partir de una región donde se conoce el valor del LS   y obtener su valor en todo el espacio.

Para resolver la siguiente ecuación de ortogonalidad:

    0 (22)

se desarrolla, para un elemento triangular lineal, del siguiente modo:

3

(  i Ni x, ) 1N1,x  2N2,x 3N3,x 

i1

3

(  i Ni y, ) 1N1,y  2N2,y 3N3,y 0

i1

(23)

Suponiendo que la información nodal de  es conocida, y suponiendo que el nodo incógnito es el tercero, podemos agrupar términos y queda:

3

4.4.     Extensión de las velocidades de crecimiento de los LS

Aún no se ha abordado la construcción de las velocidades de actualización de los LS. Se tiene que v se puede descomponer en dos componentes según las bases definidas por los LS.

Esta velocidad de crecimiento se calcula a partir de alguna ley física de propagación y crecimiento de grietas, como la ley de Paris y el estado tensional en los elementos o relaciones entre los diversos factores de intensidad de tensiones. Una vez obtenida la velocidad de crecimiento del extremo de grieta se tiene que calcular su valor en todos los nodos del dominio.

Para construir los campos de velocidades v y v en el dominio de los LS, se define primero su valor en los nodos de los elementos que contienen el frente de grieta. Entonces se extienden los campos de velocidades, de forma que estas sean normales a los LS  y , [18]:

Para v:

𝛻𝛻v = 0

(  iNi x, ) 1N1,x  2N2,x  𝛻𝛻v = 0

i1

3      3      3

(  iNi y, ) 1N1,y  2N2,y (iN Ni x, ) 3,x (iN Ni y, ) 3,y 3 0   

(27)

i1    i1 i1   

(24)

Designando

3

ex(  i Ni x, ) 1N1,x  2N2,x 

i1

3

ey(  i Ni y, ) 1N1,y  2N2,y  (25)

i1

 3   3       h(i N Ni x, ) 3,x (i N Ni y, ) 3,y 

 i1       i1   

La solución obtenida es:

 (ex ey)

₃  (26) h

Para v:

𝛻𝛻v = 0

(28)

𝛻𝛻v = 0

Se puede observar que estas ecuaciones son formalmente iguales a la ecuación de ortogonalidad entre LS. Por tanto, es posible usar el mismo algoritmo para su resolución. 

4.5. Selección de elementos y nodos

El algoritmo de selección de qué elemento y nodo se ha de resolver en cada momento, es el punto más crítico del proceso. La técnica usada se basa en construir un vector con la información del número de nodos calculados por cada elemento. Este vector permite separar los elementos en dos tipos, los que ya tienen toda la información y los que aún tienen los datos de algún nodo por obtener. El siguiente paso es seleccionar entre los elementos que aún tienen algún nodo con información incompleta, aquellos elementos a los que solo les faltan los datos de un nodo. Se aplica el proceso en estos elementos y se

Figura 3. Grieta inicial del caso de referencia de Duflot, [10].

Fuente. Elaboración modificada de [10].

Los LS iniciales, para el caso de referencia, de tipo puramente geométrico, se obtienen a partir de las siguientes expresiones:

x y 0.1 ₀  2

(29) x y 0.1

₀ 

2

Que se pueden observar en la siguiente figura, donde la línea roja indica el LS  y la azul el LS .

A continuación, definimos las velocidades con las que el extremo de grieta se va a desplazar. Estas velocidades vienen determinadas solo por razones matemáticas y no se corresponden con un caso real. Permiten estudiar la actualización en la descripción geométrica, que es el objetivo de este trabajo, de forma genérica.

v_ 0

(31) v_ 0.1

El caso propuesto no se corresponde con ningún modelo físico real, siendo solo un ejemplo matemático de descripción de la actualización de los LS. En un caso real, la descripción vendría determinada por la grieta existente y las velocidades de actualización se deberían obtener de relaciones como la ley de Paris junto con los estados tensionales en el elemento o a partir de distintas relaciones entre los factores de intensidad de tensiones.

6. ACTUALIZACIÓN DE LOS LEVEL SETS  

El método geométrico utilizado se basa en el que Duflot introduce como smooth level set method [10]. La idea es actualizar la función  mediante una expresión que origine unos gradientes suaves mediante razonamientos geométricos.

Una vez que se tiene actualizado , se actualiza  mediante una ecuación no diferencial, de forma que se cumple exactamente el criterio del nivel cero.

En primer lugar, hay que extender ortogonalmente los campos de velocidades por todo el dominio, mediante las expresiones (27) y (28). Hay dos casos para los vectores de velocidad, en el primer paso de propagación la velocidad viene dada por la ecuación (30), y en los pasos siguientes viene dada por la ecuación (31). Las figuras 5 y 6 ilustran las velocidades para el paso inicial y los siguientes, donde las flechas indican la velocidad en cada nodo.

Figura 5. Vectores de velocidades iniciales, donde el eje x la velocidad en dirección  y el eje y es la velocidad en dirección

. 

Figura 6. Vectores de velocidades en el resto de los pasos, donde el eje x la velocidad en dirección  y el eje y es la velocidad en dirección.

Una vez extendidos ortogonalmente los campos de velocidad a todo el dominio, hay que actualizar . Para realizar esta actualización, se divide el espacio de los LS en seis áreas dependiendo de la pendiente de variación de la grieta y de los valores de los LS  y .

En el sistema de coordenadas de la grieta, es decir, usando los LS  y  como sistema de coordenadas, la pendiente de crecimiento de la grieta se obtiene de la relación entre velocidades de avance de los LS o de crecimiento de la grieta:

v^

m_ (32)

v

El área 1 corresponde a una zona donde la grieta ya se ha propagado en la parte positiva de . Matemáticamente viene definida por las condiciones  m₀  0 y  ₀ m ₀ . El valor actualizado de  es el valor inicial, ₀ .

El área 2 corresponde a una zona donde la grieta ya se ha propagado en la parte negativa de . Matemáticamente se define mediante las condiciones  m₀  0 y ₀ 0 El  valor actualizado de  es el valor inicial, ₀ .

El área 3 viene dada por m₀  0 y 0   ₀ m ₀ .

El valor actualizado de  es:

 sign(  ₀)      ₀²  ₀² (33)

El área 4 se define por m₀  0 y  ₀ m ₀ . Se actualiza con:

 ^{0 }m ⁰ (34)

1m²

El área 5 corresponde a la zona de propagación de la grieta y se define por m₀  0 y ₀ 0. El nuevo  es:

 ^{0 }m ⁰ (35)

1m²

El      área 6      es     determinada       por m₀  0      y

m ₀  ₀ 0. Aquí se obtiene  de forma que mantenga la continuidad en el dominio. Concretamente, se construyen arcos de circunferencia que unen los valores de  en el área 1 con los valores en el área 5.

Con esta construcción del campo, éste ya no representa una función distancia signada, por lo que se renombra como ' .Este campo ' actualizado se calcula con la expresión:

'  (2cos2() 1) sign() b2 ac (36) b a

donde 

Figura 7. Localización de las áreas en el caso de referencia [10]. Fuente. [10].

Para obtener el otro LS , primero se actualiza el valor del LS  a un valor intermedio ' , que obtenemos mediante una rotación del LS  inicial. Esta transformación viene expresada en el sistema de coordenadas de los LS por:

 '  ⁰v^{ 0}v^ (38)

Ahora, si ₀ 0 y  '  0 , entonces  '  ₀ , con lo que no se modifican los valores del LS en la zona donde la grieta ya existe.

Dependiendo de los valores de ₀ y de  ' , para contemplar aquella zona en que la rotación ha introducido el LS  ' en la región en que la grieta ya existía, como se puede ver en la siguiente figura, es decir, si  '  ₀ , se modifica de nuevo el valor de  ' antes de aceptarlo como correcto. 

El LS  actualizado es finalmente: 

  ' v_ (40)

En el resto del dominio, si ₀ 0 y  '  0 , entonces

 '  ₀ , con lo que no se modifican los valores en la zona donde la grieta ya existe.

Posteriormente se reinicia el LS  en el caso de que se haya alejado demasiado de una función distancia. Por último, se ortogonaliza el LS  y se reinicializa a partir de su nivel cero.

No se han observado variaciones importantes al usar otros tamaños de malla o usando mallas no estructuradas en los procesos relacionados con la reinicialización y actualización de los LS. A continuación, se presentan los resultados obtenidos al aplicar las diversas opciones barajadas en la actualización mediante razonamientos geométricos.

Usando para actualizar el LS , el valor suavizado conforme varían los campos LS –Ec. (39)– se obtienen los LS mostrados en las figuras 10 a 15. Primero actualizamos los LS, figura 10, posteriormente se ortogonaliza el LS  respecto al LS , figura 11, a continuación se reinicializa el LS , figura 12, el resultado final se puede ver en la figuras 13. En las figuras 14 a 15 se presentan loa resultados finales para los dos siguientes pasos de tiempo. 

Figura 11. LS  ortogonalizado en paso 1:  en rojo y  en azul.  

azul.  

Figura 13. LS en paso 1:  en rojo,  en azul, r del sistema LS en verde y θ del sistema LS en negro.

Figura 14. LS en paso 2:  en rojo,  en azul, r del sistema LS en verde y θ del sistema LS en negro. 

Figura 15. LS en paso 3:  en rojo,  en azul, r del sistema LS en verde y θ del sistema LS en negro.  

8.          CONCLUSIONES  

Se ha desarrollado un método que permite la construcción de funciones distancia para el modelado de grietas en XFEM, partiendo del FMM. Este método es fácilmente programable, no incluye procesos iterativos – que hace que sea rápido en comparación con otras técnicas similares–, es extrapolable directamente a tres dimensiones y permite la implementación de narrow band.

Para aumentar la efectividad del método se podría mejorar el proceso de selección de qué elemento se ha de calcular en cada instante, ya que es el proceso más crítico desde el punto de vista de coste computacional y, además, puede influir en la solución.

9. REFERENCIAS

[1] J. A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science, Cambridge University Press.

Cambridge. R.U. 1999. 

[2] J. A. Sethian, “Evolution, Implementation, and application of level set and fast marching methods for advancing fronts”, Journal of Computational Physics, vol. 169, pp. 503-555, 2001.

[3] J. A. Sethian, “A Fast Marching Level Set Method for Monotonically Advancing Fronts”, Proc. Nat. Acad, vol. 93, nº. 4, pp.1591-1595, 1996.

[4] J. A. Sethian, “Fast marching methods”, Society of Industrial and Applied Math.. Review, vol 41, nº 2, pp.

199-235,1999.

[5] T. J. Barth, J.A. Sethian, ”Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains”, Journal of Computational Physics. 1998, vol.

145, pp. 1-40, 1998.

[6] N. Moës, J. Dolbow, T. Belytschko, “A Finite Element Method for Crack Growth Without Remeshing”,

International       Journal      for Numerical      Methods in

Engineering, vol. 46, nº 1, pp. 131-150, 1999. 

[7] M. Stolarska, D. Chopp, N. Moës, T. Belytschko, “Modelling crack growth by level sets in the extended finite element method”, International Journal for Numerical Methods in Engineering, vol. 51, pp. 943-960, 2001. 

[8] N. Moës, A. Gravouil, T. Belytschko, “Non-planar 3D crack growth by the extended finite element and level sets-part I mechanical model”, International Journal for Numerical Methods in Engineering, vol. 53, pp. 25492568, 2002.

[9] A. Gravouil, N. Moës, T. Belytschko, “Non-planar 3D crack growth by the extended finite element and level sets-part II: Level set update”, International Journal for Numerical Methods in Engineering, vol. 53, pp 25692586, 2002.

[10] M. Duflot, “A study of the representation of cracks with level sets”, International Journal for Numerical Methods in Engineering, vol. 70, nº11, pp. 1261-1302,  2006.

[11] N. Sukumar, D. L. Chopp, B. Moran, “Extended finite element method and fast marching method for three-dimensional fatigue crack propagation”, Engineering Fracture Mechanics, vol.70, pp. 29-48, 2003.

[12] D. L. Chopp, N. Sukumar, ”Fatigue crack propagation of multiple coplanar cracks with the coupled extended finite element / fast marching method”, International Journal of Engineering Science, vol. 41, pp.

845-869, 2003.

[13] N. Sukumar, D. L Chopp, E. Béchet, N. Moës, “Three- dimensional non-planar crack growth by a coupled extended finite element and fast marching method“, International Journal for Numerical Methods in Engineering, vol. 76, pp. 727-748, 2008.

[14] R. Kimmel, J. A. Sethian, “Computing geodesic paths on manifolds”, Proceedings of the national academy of Sciences , vol 95, nº15, pp.8431-8435, 1998. 

[15] J. A. Sethian, A. Vladimirsky, “Fast methods for the eikonal and related Hamilton-Jacobi equations on unstructured meshes”, Proceedings of the National Academy of Sciences, vol. 97, nº 11, pp. 5699-5703, 2000.

[16] D. L. Chopp, ”Some improvements of the fast marching method”, Society of Industrial and Applied Math, vol. 23, nº1,  pp 230-244, 2001. 

[17] R. N. Elias, M. A. D. Martins, A. L. G. A. Coutinho, “Simple finite element-based computation of distance functions in unstructured grids”, International Journal for Numerical Methods in Engineering, vol.72, pp. 10951110, 2007.

[18] D. Adalsteinsson, J. A. Sethian, ”The fast construction velocities in level set methods”, The Journal of Computational Physics, vol. 148, pp. 2-22, 1999.