SIMULACION NUMERICA DE LA GENERACION DE ONDAS ACUSTICAS POR IRRADIACION LASER

J.I.Etcheverry^{1}, Mariana Mesaros^{2}

1 Departamento de Matemática, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires,

Pabellón I, Ciudad Universitaria, 1428 Buenos Aires, Argentina

2 Departamento de Física, Facultad de Ingeniería, Universidad de Buenos Aires, Paseo Colón 850, Buenos Aires, Argentina

En este trabajo se analiza la producción de ondas acústicas en materiales por medio de simulaciones de dinámica molecular. Se utilizan potenciales de interacción sencillos, para investigar cuáles de los fenómenos observados son reproducibles por medio de esta técnica.

In this work we analyze the production of acoustic waves by a simple molecular dynamics simulation, in order to investigate what of the observed phenomena are reproducible within this simple framework.

1 Introducción

Las alteraciones superficiales producidas por irradiación láser en un material pueden ser caracterizadas in situ a partir del análisis de las ondas acústicas que se producen en el mismo [1, 2, 3] . La ventaja de estas técnicas acústicas es que son no destructivas, y pueden ser llevadas a cabo en tiempo real.

La producción de ondas acústicas por acción del láser está vinculada a la dilatación térmica, y a grandes intensidades al impulso que se lleva el material evaporado. Algunos de estos fenómenos pueden ser descriptos cuantitativamente mediante las ecuaciones de termo-elasticidad, y describiendo el proceso de evaporación como un proceso de equilibrio [4, 5, 6] .

Sin embargo, la producción de ondas acústicas cuando el material sufre una transición de fase, por ejemplo, sólido-líquido, es mucho más compleja, debido a cambios de volumen específico, de las propiedades térmicas, elásticas, etc. Además, para pulsos láser en el rango de los picosegundos, es de esperar que haya importantes fenómenos de no equilibrio, y que no puedan ser descritos adecuadamente en un marco clásico.

Es en esta área donde consideramos que las simulaciones de dinámica molecular, aún cuando utilicen potenciales sencillos, son valiosas para obtener una comprensión más adecuada del fenómeno.

En el trabajo aquí presentado, se analiza la posibilidad de reproducir las ondas acústicas producidas en materiales mediante una simulación de dinámica molecular de una red inicial de partículas que interactúan con un potencial de Lénard-Jones. Se muestra que es factible obtener las ondas acústicas para muestras del orden de los 250000 átomos, con tiempos de ejecución en una PC del orden de las 12 horas.

2 Algoritmo

El problema básico considerado es el de resolver las ecuaciones de movimiento de un conjunto de partículas, sometidas a una interacción mutua, por ejemplo, mediante interacción de a pares. En este caso, las ecuaciones de Newton se escriben:

m_{i} \frac{d^{2}x_{i}}{dt^{2}} = ∑_{j≠i} F_{i,j}  (1)

donde m_{i} es la masa de la partícula i, F_{i,j} es la fuerza entre las partículas i y j, y x_{i} es la posición de la partícula i.

El cálculo de las aceleraciones de todas las partículas, utilizando el hecho de que F_{i,j} = −F_{j,i}, requiere del cálculo de N(N − 1)/2 interacciones de pares. Si a esto se le añade que este cálculo debe repetirse típicamente del orden de las 10000 o 100000 veces, el costo computacional es prohibitivamente grande ya para sistemas del orden de 1000 partículas.

En el caso de fuerzas de corto alcance, la situación mejora, puesto que para calcular la fuerza sobre una partícula, sólo es necesario calcular la interacción con aquellas partículas que estén suficientemente cerca. Esta situación se explota en los algoritmos de partículas en celdas, donde esencialmente hay una partición en celdas subyacente y donde las interacciones se calculan sólo con partículas de las celdas adyacentes a la de la partícula considerada. Esto agrega algún recargo a las operaciones que deben realizarse, pues debe llevarse cuenta de qué partículas están en cada celda a cada instante de tiempo, pero tiene la ventaja de que la cantidad de operaciones por partícula crece linealmente con N, con una constante que depende esencialmente de la densidad de partículas.

En las simulaciones que se muestran en este trabajo el potencial utilizado es un potencial de Lénard-Jones:

V (r) = \frac{ϵ}{4π} ((\frac{σ}{r} )^{12} − (\frac{σ}{r} )^{6})  (2)

con parámetros correspondientes al argón

ϵ = 119.8k(K) σ = 3.406A

(donde k es la constante de Boltzmann) ( [7] ).

En todos los casos que se muestran el problema es bidimensional, con una red inicial de 601 por 401 partículas. Las partículas se distribuyeron inicialmente en una red triangular de parámetro 1.12σ, que corresponde a la posición de equilibrio de dos partículas, y con una distribución de velocidades gaussiana. Para minimizar la influencia de pequeños transitorios iniciales sobre la posición de las partículas, se dejó evolucionar el sistema por un tiempo suficientemente largo como para que llegara a posiciones y velocidades de equilibrio.

Se implementaron condiciones de contorno libres, es decir, no hay restricciones ni sobre las partículas ni sobre la temperatura. De este modo, las simulaciones corresponden a la irradiación de un pequeño rectángulo de material en vacío.

La irradiación del láser se simuló suponiendo un pulso de perfil gaussiano, tanto espacial como temporalmente. El ancho espacial es de .04 μm, y el temporal de 5 ps. Para distribuir en profundidad la energía del pulso, se supuso que el coeficiente de absorción del material es proporcional a la densidad, de tal modo que el láser se atenuaba en 1/e en .017 μm.

3 Resultados numéricos

En el esquema de la Figura 1 se muestra la disposición de la muestra, del pulso de láser, y las posiciones iniciales de tres partículas P1, P2, P3 utilizadas como control. Estas están situadas en las capas atómicas 1, 101 y 201 respectivamente, contando desde la cara inferior de la muestra.

En la Figura 2 se muestra el movimiento de las partículas P1, P2, y P3, en función del tiempo.

Figura 1: Diagrama esquemático de la muestra, del pulso láser, y de los puntos de control P1, P2, P3 que se hallan en las capas atómicas 1, 101, 201, respecto de la cara inferior de la muestra

Figura 2: Desplazamiento de los puntos P1, P2, P3 en función del tiempo. La llegada de la onda acústica se señala con O, y su rebote en la cara inferior con R.

Notar los desplazamientos producidos por la llegada de la onda acústica, y el posterior rebote de la misma. Estos han sido señalados con O y R respectivamente.

En la Figura 3 se muestra el desplazamiento de la partícula P1 para varias intensidades de pulso de láser. Notar el crecimiento de la amplitud del desplazamiento producido por la onda acústica con la intensidad del pulso, y también el hecho de que el comienzo de la onda acústica llega primero a intensidades más elevadas. Esto último se debe a que el potencial es muy repulsivo para distancias interatómicas muy pequeñas, haciendo que el material sea prácticamente incompresible para compresiones grandes y por lo tanto la velocidad del sonido aumenta para ondas acústicas de gran intensidad.

En las Figuras 4, 5 se muestran los perfiles de energía cinética media y densidad obtenidos para varios tiempos consecutivos. Para obtener dibujos Figura 3: Desplazamiento del punto P1 para varias fluencias del láser. La amplitud del desplazamiento crece con la fluencia.

significativos fue necesario hacer un filtrado digital de los datos, sobre todo en el caso de la densidad, para eliminar las fluctuaciones en el número de partículas por celda. En los gráficos de densidad se ve claramente la onda acústica longitudinal, que se propaga inicialmente hacia abajo, y su posterior rebote en la cara inferior de la muestra. En los gráficos de energía cinética media se ve un incremento de energía cinética coincidente con el pasaje de la onda acústica, y un calentamiento remanente del material debido a la conversión de energía elástica en movimiento térmico aleatorio.

4 Conclusiones

El trabajo presentado demuestra la factibilidad de determinar las ondas acústicas mediante simulaciones de dinámica molecular, con una cantidad de partículas adecuada (aproximadamente 250000 en las simulaciones mostradas), y con tiempos de corrida razonables en una PC. En principio esta técnica permitirá analizar fenómenos de no equilibrio, como fusión y evaporación de material debidas a pulsos láser de duraciones en el orden de los picosegundos, y procesos de acumulación de defectos por irradiación múltiple de una misma muestra. Este análisis será reportado en futuros trabajos, debido en parte a la gran demanda de tiempo de cómputo que requiere un estudio paramétrico de las variables involucradas.

Referencias

[1] M. Mesaros, O. Martínez, G. M. Bilmes, J. O. Tocho, J. Appl. Phys. 81 (2), 1014-1016, Figura 4: Perfiles de energía cinética media obtenidos para tiempos t=20 ps, 40 ps, 60 ps, 80 ps, 100 ps.

Figura 5: Perfiles de densidad obtenidos para tiempos t=20 ps, 40 ps, 60 ps, 80 ps, 100 ps. Notar la pluma de material evaporado.

[2] M. Mesaros, O. Martínez, G. M. Bilmes, J. O. Tocho, Anales de la Asociación Química Argentina, vol 84, n2, 167-171, (1996).

[3] V. Gusev, A. A. Kolomensky, P. Hess, Appl. Phys A, 61, 285-298, (1995).

[4] C.B. Scruby, R.J. Dewhurst, D.A. Hutchins, S.B. Palmer. J. Appl. Phys. 51 (12), (1980).

[5] R.J. Dewhurst, D.A. Hutchins, S.B. Palmer, C.B. Scruby. J. Appl. Phys. 53 (6), (1982).

[6] F.V. Bunkin, Al.A. Kolomensky, V.G. Mikhalevich. Lasers in Acoustics, Harwood Acad. Publishers, (1991).

[7] P.W. Atkins, Físico-química. Tercera edición Addison Wesley Iberoamericana, (1991).