1

DESARROLLO DE UNA APLICACIÓN PARA EL ANÁLISIS DE ESTRUCTURAS PIEZOELÉCTRICAS USANDO EL MÉTODO DE ELEMENTOS FINITOS

ALEJANDRO LÓPEZ RIVERA

UNIVERSIDAD AUTÓNOMA DE OCCIDENTE FACULTAD DE INGENIERÍA

PROGRAMA DE MAESTRÍA EN INGENIERÍA SANTIAGO DE CALI

2017

DESARROLLO DE UNA APLICACIÓN PARA EL ANÁLISIS DE ESTRUCTURAS PIEZOELÉCTRICAS USANDO EL MÉTODO DE ELEMENTOS FINITOS

ALEJANDRO LÓPEZ RIVERA

Proyecto de gado para optar al título de Magíster en Ingeniería

Director EDINGUER ENRIQUE FRANCO G.

Doctor en Ingeniería Mecánica

UNIVERSIDAD AUTÓNOMA DE OCCIDENTE FACULTAD DE INGENIERÍA

PROGRAMA DE MAESTRÍA EN INGENIERÍA SANTIAGO DE CALI

2017

Nota de aceptación:

Ing. PhD. MAURICIO BARRERA Jurado

Ing. PhD. HECTOR JARAMILLO Jurado

Santiago de Cali, 30 de marzo del 2017.

Aprobado por el Comité de Grado en cumplimiento de los requisitos exigidos por la Universidad Autónoma de Occidente para optar por el título de magíster en ingeniería con énfasis en Mecatronica.

A mis padres Leider López y María Elena Rivera, a los que amo profundamente.

AGRADECIMIENTOS

Agradecimientos a mi familia, a mi padre Leider López, mi madre María Elena Rivera, hermano Héctor Daniel López y mi hijo Diego Alejandro López, quien con su apoyo y motivación son motores de mi vida, a mi director de proyecto de grado, Doctor Edinguer Enrique Franco, por su guía y enseñanza a lo largo de mi proceso educativo, a Karol Vanessa Lenis por su amor y compañía.

CONTENIDO

pág. ABSTRACT 11

RESUMEN 12

INTRODUCCIÓN 13

2. OBJETIVOS 16

2.1. OBJETIVO GENERAL 16

3. JUSTIFICACIÓN 17

4. ANTECEDENTES 19

5. MARCO TEÓRICO 21

5.1 FORMULACIÓN DEL ELEMENTO FINITO PIEZOELÉCTRICO 21

5.2 IMPLEMENTACIÓN NUMÉRICA 24

5.3 ANÁLISIS ARMÓNICO 29

5.4 CALCULO DE LA IMPEDANCIA ELÉCTRICA 30

6. METODOLOGÍA EXPERIMENTAL 32

6. RESULTADOS 35

6.1 VERIFICACIÓN DEL SOLUCIONADOR DEL PROBLEMA MECÁNICO 35

6.2 RESPUESTA EN FRECUENCIA DE LOS ELEMENTOS PIEZOELÉCTRICOS 39

6.2.1 APC-181. 39

6.2.2 APC 613. 43

6.2.3 APC-1029. 46

6.3 ANÁLISIS DE UN ACTUADOR TIPO LANGEVIN 50

6.4 ANÁLISIS DE LA DEPENDENCIA DE LA MALLA 53

6.5 COMPARACIÓN CON EL SOFTWARE COMERCIAL COMSOL 57

7 CONCLUSIONES 59

7.1 TRABAJOS FUTUROS 61

BIBLIOGRAFÌA 62

LISTA DE FIGURAS

pág.

Figura 1. Hexaedro de 8 nodos con sus coordenadas naturales ( ηyμξ, ). 24 Figura 2. Montaje experimental usado en la medición de la 31 impedancia eléctrica de elementos piezoeléctricos. Figura 3. Desplazamiento (arriba) y esfuerzo normal (abajo) 36 de la viga sometida a carga axial. Figura 4. Desplazamiento (arriba) y esfuerzo normal (abajo) 37 de la viga sometida a carga transversal. Figura 5. Cerámicas APC de prueba. 39 Figura 6. Malla usada para el análisis del elemento 39 APC-181 (3250 nodos y 2520 hexaedros). Figura 7. Impedancia eléctrica (ohm) en función de la 40 frecuencia para el elemento APC-181. Figura 8. Principales modos de vibración observada en el elemento APC-181. 42 Figura 9. Impedancia eléctrica en función de la frecuencia calculada 42 por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-181. Figura 10. Impedancia eléctrica (ohm) en función de la 43 frecuencia para el elemento APC-613. Figura 11. Principales modos de vibración observada en el elemento APC-181. 45 Figura 12. Impedancia eléctrica en función de la frecuencia 45 calculada por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-613. Figura 13. Impedancia eléctrica (ohm) en función de la 46 frecuencia para el elemento APC-1029. Figura 14. Principales modos de vibración observada en 47 el elemento APC-1029.

Figura 15. Impedancia eléctrica en función de la frecuencia calculada 48 por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-1029. Figura 16. Aspecto del actuador tipo Langevin analizado. 51 Figura 17. Impedancia eléctrica del actuador tipo Langevin 51 calculada por elementos finitos. Figura 18. Modos de vibración observados a las frecuencias de 52 antiresonancia (izquierda) y resonancia (derecha) del actuador tipo langevin. Figura 19. Viga en cantilever de material piezoeléctrico. 54 Figura 20. Impedancia eléctrica en función de la frecuencia de la 54 viga en cantiliver para los tres diferentes tamaños de malla. Figura 21. Modo de vibración fundamental de la viga empotrada 56 obtenido con las diferentes mallas. Figura 22. Comparación de la impedancia eléctrica en función 58 de la frecuencia para un disco piezoeléctrico obtenida por diferentes métodos.

LISTA DE TABLAS

pág.

Tabla 1. Coeficientes para el método de cuadratura de Gauss 27 con 3=n . Tabla 2. Propiedades físicas de la cerámica PZT usadas en 33 este trabajo. Tomadas de literatura de Boucherd. Tabla 3. Comparación de los resultados analíticos y por elementos 37 finitos del problema lineal elástico analizado. Tabla 4. Dimensiones de los elementos piezoeléctricos analizados. 39 Tabla 5. Propiedades del actuador analizado. 51 Tabla 6. Frecuencia del modo fundamental para tres malla 55 con diferente número de elementos.

ABSTRACT The work describes the development of a code in Matlab, which, it use the Finite Element Method (MEF), it is intended to provide support for the development of piezoelectric structures. It explained how by means of the (MEF) the piezoelectric continuum is divided, in subdomains, and as the solution of each of them, allow to obtain the behavior of the piezoelectric structures. Later In order to corroborate the operation of the code, tests are performed comparing the results of this, with a case study, the case of the cantilevered beam, then compared with the experimental results of two discs and ring of piezoelectric material, with a sonotrode, it is an object that brings together passive elements with an active or piezoelectric element, the work probe the relation of the mesh with the results and analyzed and the comparison with a commercial software.

RESUMEN El trabajo describe el desarrollo de un código en Matlab, que usando el método de Elementos Finitos (MEF), pretende brindar un apoyo para el desarrollo de estructuras piezoeléctricas. Se explicará como por medio del (MEF) se divide el continuo piezoeléctrico, en subdominios, y como la solución de cada uno de estos, permiten obtener el comportamiento de las estructuras piezoeléctricas a diseñar. Para corroborar el funcionamiento del código, en primer lugar, se compara los resultados del código con los resultados teóricos, del estudio de una viga en voladizo sometida a una carga. En segundo lugar, se comparan los resultados del código, con resultados experimentales de dos discos y un anillo piezoeléctricos. En tercer lugar, se analiza la respuesta del comportamiento de un sonotrodo. Finalmente se revisa la dependencia de la malla con la respuesta obtenida y se realiza la comparación con un software comercial. El desarrollo de este código facilitara el diseño de estructuras piezoeléctricas. Palabras clave: Ingeniería Mecánica, Método de Elementos Finitos, Diseño de elementos piezoeléctricos, anillo piezoeléctrico, Análisis Estructural. Sonotrodo.

INTRODUCCIÓN La piezoelectricidad es la característica que poseen algunos materiales de generar una tensión eléctrica en respuesta a una tensión mecánica aplicada. Este comportamiento electromecánico es reversible, pues al aplicar una tensión eléctrica se genera deformación mecánica. La existencia de la piezoelectricidad fue sugerida a mediados del siglo XIX, a partir de resultados teóricos y el estudio de la piroelectricidad1. Sin embargo, fue hasta 1880 cuando se demostró experimentalmente el efecto piezoeléctrico en compuestos naturales como la turmalina y el cuarzo2. La primera aplicación de los materiales piezoeléctricos fue el sonar, desarrollado en Francia en los años de la primera guerra mundial para contrarrestar la amenaza de los submarinos alemanes3. El material usado para los primeros sensores fue el cuarzo. A partir del sonar fueron encontradas muchas aplicaciones industriales, la mayoría relacionadas con la inspección no destructiva de materiales y la generación de vibraciones mecánicas de alta intensidad y alta frecuencia. En la década de 1940 se realizaron las primeras ecografías, mostrando su potencial para aplicaciones médicas4. En la década de 1950 se descubrió en Japón el material cerámico PZT (zirconato titanato de plomo), que es el material con mayor efecto piezoeléctrico conocido. Actualmente los materiales piezoeléctrico tienen una gran variedad de aplicaciones, como transductores para la generación y recepción de ondas acústicas, micromecanismos (MEMS), sistemas de disparo de bolsas de aire (airbags), sistemas de ignición en motores de combustión, actuadores para el posicionamiento fino en óptica, memorias no volátiles de computador, acelerómetros, entre muchas otras. También se destacan las aplicaciones industriales de potencia, como limpieza, soldadura y atomización ultrasónicas5. Muchos dispositivos y equipos tienen componentes piezoeléctricos incorporados, por tanto, el diseño y mantenimiento de estos equipos requieren dominio del tema. 1 ERHART Mgr Jiri. Piezoelectricity and ferroelectricity Phenomena and properties. Department of Physics FP TUL. Presentación de clase, año no disponible. 1-65 p 2 CURRIE J, Currie P. Développement, par pression, de electricité polaire dans les cristaux hémièdres à faces inclinées. Comptes rendus. Académie des Science. 1880. 294-295 p. 3 KINO Gordon S. Acoustic Waves: Devices, Imaging, and Analog Signal Processing. Vol 2, Prentice Hall, 1987. 531 p 4 CHEEKE J. David, Zagzebski J. Fundamentals and Applications of Ultrasonic Waves. Physics Department Concordia University. Montreal: CRC Press LLC 2004. 72-719 p. 5 FRANCO Edinguer E, Cuervo Andrés M, Barrera Helver M. Desarrollo de un programa por elementos finitos para analizar la respuesta vibratoria de estructuras piezoeléctricas. Facultad de Ingeniería, Universidad Autónoma de Occidente, Cali, 2016. 1-14 p.

La respuesta electromecánica de un dispositivo piezoeléctrico depende de la excitación, las características de los materiales y la geometría del dispositivo. La principal dificultad en el análisis dinámico es la geometría, porque en la mayoría de los casos es complicada y los métodos analíticos no son útiles para representarlas. Por tanto, se hacen muchas simplificaciones, obteniendo resultados más cualitativos que cuantitativos6. El método de los elementos finitos (MEF) se ha empelado ampliamente para el análisis de la respuesta en frecuencia de estructuras piezoeléctricas debido a que permite representar las diferentes geometrías de manera eficiente7. El MEF divide el dominio en subdomminios de geometría simple. La solución de la ecuación diferencial constitutiva para esta geometría simple es conocida. El problema global se soluciona sumando la contribución de cada subdominio, haciendo que el campo de respuesta sea suave y cumpla las condiciones de frontera. Cada uno de estos subdomminios y su respectiva formulación matemática, se denomina elemento finito. El MEF permite la simulación numérica de la respuesta dinámica del dispositivo piezoeléctrico. Como se tienen en cuenta la geometría del dispositivo, los resultados son más precisos y proporcionan más información que los métodos analíticos. La simulación numérica permite predecir el comportamiento dinámico y evaluar el efecto de cambios en la geometría o en el material. Esto reduce el costo de desarrollo, pues el trabajo experimental requerido es menor. Además le proporciona al diseñador más libertad, pues se pueden evaluar geometrías más elaboradas, lo que beneficia la innovación.

6 PIEFORT Vincent, Preumont André. Finite Element Modelling of Piezoelectric Active Structures. Active Structures Laboratory, ULB-CP 165/42, Brussels, 2001. 1-17 p 7 Ibid, p 1-17.

En este trabajo se propone el desarrollo de un código por elementos finitos para el análisis dinámico de estructuras piezoeléctricas. El código consiste en una serie de rutinas en Matlab que implementan la solución del problema dinámico armónico con elementos tipo hexaedro con interpolación lineal. Las mallas de cálculo serán elaboradas en el software libre Salome y exportadas en un formato que pueda ser leído por Matlab. El postprocesamiento de los resultados será realizado en el software Paraview

2. OBJETIVOS 2.1. OBJETIVO GENERAL Desarrollar un código, en el programa Matlab, el cual, aplicando el método de elementos finitos (FEM), permita analizar el comportamiento electromecánico de estructuras piezoeléctricas. 2.2 OBJETIVOS ESPECÍFICOS Desarrollar la formulación de un elemento finito piezoeléctrico tipo hexaedro. Programar en Matlab un código de análisis que permita obtener el desplazamiento temporal, cuando se aplica un campo eléctrico. Probar el código, analizando el comportamiento vibratorio eléctrico de elementos piezoeléctricos comunes (discos y anillos). Simular un comportamiento vibratorio de un sistema que incluya un material piezoeléctrico activo y un material pasivo acoplado.

3. JUSTIFICACIÓN En el proceso de optimización de diseño de los elementos piezoeléctricos, por lo general, estos se someten a ciclos de diseño, que componen pruebas de ensayo y error, los cuales son procesos de diseño costosos y demorados. Durante la última década se ha avanzado en la investigación del modelado por el método de elementos finitos de elementos piezoeléctricos, lo que ha logrado mejorar la calidad de los dispositivos y reducir el tiempo de desarrollo. Sin embargo, el software disponible para este fin es comercial y deben adquirirse licencias de software muy costosas. El análisis piezoeléctrico es un tema específico que hace parte de varios programas comerciales denominados multifísicos. Por tanto, se requiere adquirir un paquete grande y costoso de software, que queda muchas veces subutilizado. Siendo consecuente con esta idea, si se requiere un software para análisis piezoeléctrico, es mejor desarrollarlo, la propuesta de este trabajo, es el desarrollo de un código, usando las herramientas de software libre de enmallado y post procesamiento disponibles. Por otro lado, dispositivos piezoeléctricos tales como sensores para ensayos no destructivos, actuadores de potencia para aplicaciones industriales, acelerómetros, entre muchos otros, son productos de alto valor agregado, donde el conocimiento y experiencia representan la mayor parte de su precio. El desarrollo local de este tipo de productos puede ayudar a mejorar la competitividad de la industria local y ahorrar dinero en la construcción de los mismos. Las licencias comerciales de programas de elementos finitos, como Ansys o Comsol, que realizan análisis piezoeléctricos son caros, y muchas empresas pequeñas o ingenieros independientes no tienen recursos financieros para adquirirlos. Aunque Matlab también es un software licenciado, el programa desarrollado puede traducirse fácilmente a otro entorno de programación, como por ejemplo Python Numeric, C/C++ o Fortran. Por otro lado, como el código en Matlab muestra todos los detalles de la implementación del cálculo, se tiene un conocimiento profundo del problema y la posibilidad de cambiar el código, haciéndolo muy flexible. Este enfoque busca generar un impacto en el gasto que se genera en el momento de adquisición de una licencia comercial, esto se logra implementando el trabajo con herramientas de software libre, convirtiéndose en una herramienta de gran utilidad, para ingenieros en el campo industrial o académico, además de aportar a

la apropiación del conocimiento personal y el de las personas que aportan en este campo, permitiendo así un trabajo colaborativo en pro de nuevos desarrollos en este campo. Es así, como el código desarrollado será usado para analizar los dispositivos piezoeléctricos requeridos en un proyecto de investigación sobre atomización ultrasónica que está desarrollando el grupo de investigación PAI+ del Departamento de Energética y Mecánica de la UAO. Por otro lado, dispositivos piezoeléctricos tales como sensores para ensayos no destructivos, actuadores de potencia para aplicaciones industriales, acelerómetros, entre muchos otros, son productos de alto valor agregado, donde el conocimiento y experiencia representan la mayor parte de su precio. El desarrollo local de este tipo de productos puede ayudar a mejorar la competitividad de la industria local y ahorrar dinero en la construcción de los mismos. Teniendo un modelo electromecánico, basado en el método de elementos finitos, este se podría ser usado para obtener un posible comportamiento de la salida del elemento piezoeléctrico. Permitiendo mejorar el rendimiento del transductor. El diseño es estudiado con respecto a la magnitud de salida y una respuesta de frecuencia plana. Los resultados del modelado serían utilizados para facilitar la fabricación de un prototipo. El diseño y análisis de dispositivos piezoeléctricos requiere el modelado numérico de su respuesta en frecuencia. Los métodos analíticos son útiles con geometrías simples, donde es viable asumir que las deformaciones mecánicas varían en una sola dirección. Esta simplificación del problema elimina la mayoría de modos de vibración, proporcionando resultados que pueden ser incompletos o inexactos. El análisis por medio del MEF proporciona soluciones más detalladas. Esta información adicional le permite al diseñador tomar decisiones más acertadas y reducir el trabajo experimental requerido en el desarrollo de un dispositivo, lo que significa reducción de costos y mejora de la calidad del producto.

4. ANTECEDENTES El rango de audición del humano se encuentra en el ancho de banda comprendido, entre los 20 Hz y los 20 kHz. Con aplicaciones como los sonares, se comenzó a trabajar en el rango de los 20 Khz a 100 Khz, y con el uso y descubrimiento de las características de los elementos piezoeléctricos, las aplicaciones han permitido trabajar en rangos muy altos, como los Mhz. El uso de elementos que tienen características piezoeléctricas es muy normal en el desarrollo de componentes tecnológicos, un ejemplo de ello, son los acelerómetros, que en la industria son ampliamente usados para realizar pruebas no destructivas, en el análisis de vibraciones mecánicas, en la industria aérea, se encuentran en sistemas de control de vuelo y en misiles, en la industria automotriz, en los sistemas de accionamiento principal de bolsas de aire, y en productos de consumo masivo, como; celulares y tablets. Como actuadores se encuentran en posicionadores de precisión en telescopios, microscopios, cabezales de lectura de discos duros, micro mecanismos, entre otros. En el campo de la microelectrónica son la base de las memorias no volátiles FeRAM, osciladores y filtros, y en el campo del ultrasonido de potencia, los materiales piezoeléctricos, se usan para crear vibraciones de alta potencia, permitiendo aplicaciones en soldadura y corte de materiales, limpieza de piezas, atomización, eliminación de espuma en procesos de fermentación, entre otros. En el campo de la medicina, el uso del ultrasonido, se observa en la adquisición de imágenes, con el ecógrafo en 3D. Los ejemplos mencionados, conllevan a notar la relevancia de los materiales piezoeléctricos en los procesos de la industria moderna y con el paso del tiempo nuevas aplicaciones van haciendo su aparición. Para el desarrollo de elementos piezoeléctricos en el pasado, el desarrollo de transductores piezoeléctricos se basaba en investigaciones experimentales, las cuales involucraban varias etapas de diseño y ciclos de fabricación, repitiéndose, hasta que se cumplieran las especificaciones óptimas8. Posteriormente se comienza a trabajar en procesos de optimización, buscando el modelado de los elementos piezoeléctricos, con esto, se pretende cambiar el procedimiento basado en ensayo y error, para obtener un ciclo de desarrollo más corto, buscando el menor consumo de tiempo. 8 SHIARI Behrouz, Wild Peter M. Coupled finite element analysis of a piezo-ceramic force transducer. Deparment of Mechanical & Aerospace Engineering, Carlenton University, Ottawa, Canada, 2006. 167-181 p.

Por desgracia, la mayor dificultad técnica en la realización de una exitosa simulación de transductores es la falta de comprensión de las interacciones entre la parte estructural y campos eléctricos. Por ello, aunque el método de elementos finitos es una técnica que se usa en macro estructuras, se revisa y aplica la misma técnica en micro estructuras, para el diseño de sistemas que usan piezas piezoeléctricas9. El análisis por elementos finitos de estructuras piezoeléctricas es un tema relativamente nuevo que ha sido abordado por diferentes autores. Avdiaj et al. Dedujeron la formulación del elemento tipo hexaedro lineal y los usaron para analizar los modos de vibración de geometrías simples (paralelepípedos)10. Piefort hizo un estudio detallado del tema, analizando la deducción de las ecuaciones constitutivas a partir de la termodinámica, la obtención de una ecuación de movimiento usando el principio de Hamilton y particularizando esta ecuación de movimiento para varios elementos finitos bidimensionales y tridimensionales, y reportando varias aplicaciones11. Hasta donde se sabe, no hay un software libre desarrollado específicamente para analizar estructuras piezoeléctricas por elementos finitos. Code-Aster es un software libre para análisis estructural para ingenierías mecánica y civil. Este software no tiene elemento piezoeléctrico, sin embargo, haciendo ciertas manipulaciones se pueden solucionar problemas piezoeléctricos con elementos térmicos, pero es un proceso complicado y del cual hay muy poca documentación.

9 Ibid., p 167-181. 10 AVDIAJ Sefer, Setina Janez, Syla Naim, Modeling of the piezoelectric effect using the finite – element method, University of Prishtina, Faculty of Natural Science and Mathematics, Prishtina, Kosovo, 2009. 1-9 p. 11 PIEFORT, Op. cit. p 1-17.

5. MARCO TEÓRICO 5.1 FORMULACIÓN DEL ELEMENTO FINITO PIEZOELÉCTRICO Las ecuaciones constitutivas para un material piezoeléctrico, con comportamiento lineal, son las siguientes:

EeSc=T TE (1) Eξ+Se=D S (2)

Donde TTTTTTT=T 121323332211 es el vector de esfuerzo, TTTTTSS=S 121323332211 el vector deformación 321 EEE=E el vector de campo eléctrico, 321 DDD=D el vector de desplazamiento eléctrico, Ec es la matriz de rigidez del material medida en un campo eléctrico constante, ξ es la matriz dieléctrica medida a deformación constante, e la matriz de efecto piezoeléctrico12. Las ecuaciones dinámicas de un continuo piezoeléctrico pueden obtenerse usando el principio de Hamilton, definiendo un funcional, que consiste en la diferencia entre las energías deformación y electrostática, y aplicando la ecuación de Euler-Lagrange para establecer el principio variacional que gobierna el material piezoeléctrico. La ecuación resultante es la que se discretiza, usando el método de elementos finitos y luego es simplificada verificando cualquier variación arbitraria de las variables compatible con las condiciones de frontera esenciales13. En la discretización usando MEF, el campo vectorial de desplazamiento {u } y el campo escalar de potencial eléctrico φ se expresan como una combinación lineal de los valores nodales y las funciones de forma o interpolación,

12 FRANCO, Op. cit. p. 1-14. 13 PIEFORT, Op. cit. p 1-17.

iu uN=u (3)

iφ φN=Φ . (4) Aquí, uN y φN son las matrices de funciones de forma y iu y iφ son los vectores de desplazamiento mecánico y potencial eléctrico nodales, respectivamente. La deformación mecánica S y el campo eléctrico E se relacionan con los desplazamientos y el potencial eléctrico, como se observa en las siguientes ecuaciones:

iuiu uB=uNL=S (5)

iφiφ φB=φN=E ▽ (6) Donde L es el operador de derivación que relaciona el desplazamiento y la deformación y es obtenida de la teoría de la elasticidad. La ecuación de movimiento para un elemento finito tiene la forma14.

iiuφiuuii f=φK+uKuC+uM }]{[}]{[...

(7)

iiφφiφu g=φK+uK (8) La ecuación (7) corresponde a la ecuación de movimiento de un sistema elástico, discreto de varios grados de libertad. La naturaleza del modelo determina la forma de las matrices de masa ( M ), rigidez ( uuK ) y amortiguamiento ( C ). El acoplamiento electro-mecánica se realiza por medio del cuarto término del lado izquierdo, que contiene la matriz ( uφK ) que depende de las constantes elásticas y piezoeléctricas y el vector de potencial eléctrico ( φ ). La ecuación (8) modela el problema eléctrico, que también está acoplado con la parte mecánica por medio del primer término del lado izquierdo. El segundo término modela la carga eléctrica generada por el potencial eléctrico aplicado. Es decir, esta ecuación es la suma de las cargas eléctricas generadas por el desplazamiento mecánico y el potencial 14 Ibid., p 1-17.

eléctrico. ig Es un término fuente que se puede usar para fijar un valor de carga eléctrica. Las matrices de las ecuaciones (7) y (8) se obtienen a partir de la solución del problema variacional y sus valores son los siguientes15.

uu

Tuφφu

φV

Tφφφ

φT

V

Tuuφ

uE

V

Tuuu

uT

uV

Kβ+Mα=CK=K

dVBεB=K

dVBeB=K

dVBcB=K

dVNNρ=M

(9)

Donde ρ es la densidad, α y β son los coeficientes del modelo de amortiguamiento de Rayleigh, los términos if y ig representan la fuerza mecánica y la carga eléctrica aplicadas, cuyo valor es

cTuS

sT

uV

bT

ui PN+dSPN+dVPN=f (10)

QφNσdSN=g TuS

Tφi (11)

Donde bP es una fuerza distribuida en el volumen (fuerza de cuerpo), sP es la fuerza distribuida (presión) aplicada en la frontera y cP son las fuerzas nodales externas. En la ecuación (11), σ y Q son las cargas eléctricas distribuidas y puntuales aplicadas, respectivamente.

15 FRANCO, Op. cit. p. 1-14.

5.2 IMPLEMENTACIÓN NUMÉRICA La implementación numérica del problema comienza con la particularización de las matrices definidas en el conjunto de ecuaciones (9) para un elemento finito determinado. En este trabajo se trata de un hexaedro con interpolación lineal, como se muestra en la Figura 1. Cada uno de los nodos de la figura, tiene 4 grados de libertad, los cuales son; tres de desplazamiento y el de potencial eléctrico, teniendo en conjunto, un total de 32 grados de libertad16. Figura 1. Hexaedro de 8 nodos con sus coordenadas naturales ( ηyμξ, ).

Fuente: FRANCO Edinguer E, Cuervo Andrés M, Barrera Helver M. Desarrollo de un programa por elementos finitos para analizar la respuesta vibratoria de estructuras piezoeléctricas. Facultad de Ingeniería, Universidad Autónoma de Occidente, Cali, 2016. 1-14 p. Los grados de libertad mecánicos y eléctricos son organizados en dos vectores separados, como se muestra a continuación:

Twvuwvuwvuwvu=u 888333222111 .... (12)

Tφφφφφφφφ=φ 87654321 (13) La interpolación del desplazamiento y el potencial eléctrico en una posición arbitraria dentro del elemento se hace por medio de la matriz de interpolación, definida por las funciones de forma en coordenadas naturales. Las funciones de forma para un hexaedro con interpolación lineal son: 16 Ibid., p 1-14.

iiii μμ+ηη+ξξ+=N 11181 (14)

Donde iii yμη,ξ son las coordenadas naturales del nodo i . Se tienen ocho funciones de forma, una por cada nodo, las cuales son polinomios de primer grado que deben cumplir las siguientes condiciones:

nodos otros0 nodo elen 1 i

=N i (15)

1=N i Para todo iii yμη,ξ . (16) Las matrices de funciones de forma quedan definidas por las Ecuaciones (3) y (4) y su forma es:

821

821

821

00...0000

00...0000

00...0000

NNN

NNN

NNN=N u (17)

87654321 NNNNNNNN=N φ (18) La matriz 243xuN interpola los desplazamientos, que es un campo vectorial tridimensional, mientras 81xφN interpola el potencial eléctrico, que es un campo escalar.

La matriz de derivación L se obtiene de la teoría de la elasticidad17 y está dada por:

δxδ

δzδ

δyδ

δzδ

δxδ

δyδ

δzδ

δyδ

δxδ=L

0

0

0

00

00

00

(19)

Al aplicar el operador diferencial L a uN se obtiene la matriz 246xuB . Por otro

lado, el operador gradiente está dado por, T

δzδ

δyδ

δxδ=

▽ y al aplicar este operador

diferencial a la matriz φN se obtiene la matriz 83XφB . Queda claro que las matrices 246xuB y 83XφB contienen derivadas de las funciones de forma con respecto a las coordenadas globales. Estas derivadas se ven afectadas por la transformación de coordenadas y usando la regla de la cadena se tiene:

δzδNδyδNδxδN

J=

δzδNδyδNδxδN

zyx

zyx

zyx

=

δμδNδηδNδξδN

i

i

i

i

i

i

i

i

i

(20)

17 AULD B. A. Acoustic Fields and Waves in Solids, Volumes I and II. Vol 8. Florida: Krieger Publishing Company, 1975. 405 p.

Donde la matriz J se conoce como el Jacobiano y define la transformación entre coordenadas naturales (locales) y las globales (cartesianas). Las derivadas de las funciones de forma con respecto a las coordenadas cartesianas se obtienen invirtiendo el sistema (pre multiplicando por 1J ), cuyos los términos están dados por:

ii

ii

ii

ii

ii

ii

ii

ii

ii

zδμδN

=δμδzz

δηδN

=δηδzz

δξδN

=δξδz

yδμδN

=δμδyy

δηδN

=δηδyy

δξδN

=δξδy

xδμδN

=δμδxx

δηδN

=δηδxx

δξδN

=δξδx

(21)

El determinante del Jacobiano es usado para transformar las integrales desde el sistema de coordenadas globales al sistema de coordenadas locales, de la siguiente manera:

J=dxdydz=dV , (22) Logrando cambiar los límites de integración y facilitando la aplicación del método de cuadratura de Gauss18. Las integrales mostradas en (9) son calculadas como se muestra a continuación:

n

i=

n

=j

n

=kkjikji

V

ωωωζη,ξf=dζdηdξηζξ,f=dVzy,x,f=I1 1 1

1

1

1

1

1

1

(23)

18 BATHE K.-J. Finite Element Procedures, New Jersey: Prentice Hall,1996. 1037 p.

Los iω corresponde a los valores de ponderación para las posiciones ix y de forma análoga para las otras direcciones, La integral calcula en una nube de nnn puntos dentro del elemento hexaédrico. Por ejemplo, para la matriz de rigidez el integrando en coordenadas locales es Jdetηζξ,Bcηζξ,B=ηζξ,f uETu . Los otros integrandos se manejan de forma análoga. Como el método de cuadratura de Gauss calcula el valor exacto de la integral de un polinomio de grado 12 n , donde n es número de punto de la cuadratura, se usó un valor de 3=n para precisión completa. Los valores de ix y iω se muestran en la Tabla 1. Tabla 1. Coeficientes para el método de cuadratura de Gauss con 3=n .

x i ωi

5/3 5/9 0 8/9

5/3 5/9 Las cerámicas piezoeléctricas PZT (circonato-titanato de plomo) corresponden a un material cristalino con simetría mmC6 . Suponiendo que la dirección de la polarización es el eje z , las matrices de rigidez, piezoelectricidad y dieléctrica, están dadas por19:

E

E

E

EEE

EEE

EEEE

c

c

c

ccc

ccc

ccc=c

44

44

66

331313

131112

131211

00000

00000

00000

000

000

000

19 AULD. Op. cit.

T

e

e

e

e

e=e

15

15

33

13

13

00

00

000

00

00

00

S

S

S

ε

ε

ε=ε

33

11

11

00

00

00 (24)

5.3 ANÁLISIS ARMÓNICO Se pueden realizar varios tipos de análisis dinámicos (transitorio, armónico y modal) con diferentes tipos de excitación. La mayoría de las veces, al trabajar con materiales piezoeléctricos se desea obtener la respuesta en frecuencia para localizar las resonancias, que son determinadas frecuencias a las cuales el piezoeléctrico vibra con mayores amplitudes. El análisis consiste en determinar los desplazamientos para un voltaje aplicado en los terminales del piezoeléctrico. Como el voltaje de entrada es conocido y se supone que no hay fuerzas externas aplicadas, la Ecuación de movimiento (1) queda:

iuφiuuii φK=uK+CM uu }]{[}{][...

, (25)

Donde se puede ver que le término del acoplamiento electro-mecánico pasa a convertirse en un término fuente que modela la excitación. La Ecuación (2) no es necesaria. El análisis armónico se realiza mediante la aplicación de un potencial eléctrico (voltaje) armónico, para obtener la respuesta (desplazamientos) en el estado

estable del sistema, suponiendo que la fase de la señal de excitación es cero, el

potencial eléctrico, está dado por jωωeφ=φ 0 y la respuesta es im

jωωte uj+u=u=u 0* , donde el superíndice (∗) significa que se trata de una

cantidad compleja. De esta manera, el desplazamiento complejo para una frecuencia angular kω se obtiene solucionando el siguiente sistema lineal:

0

0

²

² φK=

u

u

MωKCω

CωMωK uφ

kuuk

kkuu (26)

La señal de excitación debe aplicarse de una manera especial. Los electrodos en ambas caras de las cerámicas son superficies equipotenciales, es decir, se encuentran al mismo voltaje. Si los voltajes son diferentes en ambas caras, se genera un campo eléctrico E constante, de igual forma que en un condensador de placas paralelas, donde el voltaje varía linealmente a través del espesor de la cerámica. Por tanto, el voltaje de los nodos que está fuera de los electrodos se puede calcular interpolando los voltajes de los electrodos20. 5.4 CALCULO DE LA IMPEDANCIA ELÉCTRICA La curva de impedancia eléctrica en función de la frecuencia ωZ es la manera más usual de representar la respuesta en frecuencia de dispositivos piezoeléctricos. En esta curva, los puntos de valor mínimo corresponden a las resonancias del sistema, relacionada cada una con un modo específico de vibración. La impedancia eléctrica se obtiene a partir de la carga eléctrica en el electrodo de la cerámica que no se encuentra conectado a tierra21, por tanto:

e

ee

jωωV

=IV

=ωZ (27)

20 FRANCO, Op. cit. p. 1-14. 21 ANDRADE, Op. cit. p. 1674-1683.

Donde eV corresponde a la amplitud del voltaje de excitación y eQ es la carga eléctrica del electrodo, calculada como:

S

Te dSnD=ωQ (28)

Donde n corresponde a un vector unitario normal al electrodo. El desplazamiento eléctrico D se obtiene al evaluar la Ecuación (2) usando los desplazamientos calculados en el problema armónico:

iuφSru uBej+φBεuBe=D . (29)

En la Ecuación (29) se realiza el cálculo en coordenadas locales (a nivel del elemento) Por otro lado, eQ es una cantidad compleja, igualmente Z , como se espera. La integral de la ecuación ωQe fue aproximada calculando el promedio en el electrodo del desplazamiento eléctrico en la dirección z y multiplicando por el área.

6. METODOLOGÍA EXPERIMENTAL La Figura 2 muestra el montaje experimental usado para medir la impedancia eléctrica de los elementos piezoeléctricos. La excitación es una señal sinusoidal a bajo voltaje (menor a 10 Vpp) y frecuencia variable, proveniente del generador de señales (Keysight 33210A). Un osciloscopio (Keysight DSO-X 2022A con ancho de banda de 200 MHz) se usa para adquirir las señales de excitación 0a y un voltaje proporcional a la corriente que circula por la cerámica 1a . Un circuito electrónico especialmente diseñado permite la adquisición simultánea de las dos señales, protegiendo eléctricamente los instrumentos. El sistema es controlado por un computador a través de una red LAN, de tal manera que se pueden almacenar ambas señales para cada frecuencia en un rango deseado y de manera automática. Las señales almacenadas son procesadas posteriormente para obtener la impedancia eléctrica. Figura 2. Montaje experimental usado en la medición de la impedancia eléctrica de elementos piezoeléctricos.

El procesamiento de las señales medidas consiste en calcular la Amplitud de la señal de salida 1a , obteniendo el valor de voltaje real, según los parámetros de adquisición proporcionados por el osciloscopio. Esta señal de voltaje se convierte en corriente, al dividir por la resistencia de la carga. Después, usando la transformada de Fourier, se calcula la diferencia entre las fases de las señales 1yaa0 , permitiendo el cálculo de la parte real ri e imaginaria ii de la corriente. La impedancia eléctrica se calcula usando la siguiente ecuación:

ir

ee

ji+iV

=IV

=Z

(30)

Como la fase de la señal de excitación no varía y se usa como referencia, el valor de eV es siempre una constante real. A pesar de que las cerámicas fueron compradas a un fabricante reconocido, las especificaciones técnicas que ellos proporcionan no son completas. Para no realizar una simplificación del modelo de material, se optó por usar las propiedades físicas reportadas en la literatura de Boucher22. La Tabla 2 muestra los valores correspondientes al material PZT-5H, que se trata de una cerámica PZT dura, como las usadas para aplicaciones de potencia, que es precisamente el tipo de cerámicas existentes.

22 BOUCHER Didier, Lagier Michel, Maerfeld C. Computation of the Vibrational Modes for Piezoelectric Array Transducers using a Mixed Finite Element-Perturbation Method. IEEE Transactions Sonics and Ultrasonics. 1981. 318-329 p.

Tabla 2. Propiedades físicas de la cerámica PZT usadas en este trabajo. Tomadas de literatura de Boucherd23

Propiedad Valor Unidad

c11E

13,7 × 1010

N/m² c33E

12,4

c44E

3,14

c66E

3,4

c12E

6,97

c13E

7,16

e33 23,24 C/m²

e13 -6,62

e15 17,03

ε11S /ε0 1704 -

ε33S /ε0 1434

ρ 7600 Kg/m³

Fuente: BOUCHER Didier, Lagier Michel, Maerfeld C. Computation of the Vibrational Modes for Piezoelectric Array Transducers using a Mixed Finite Element-Perturbation Method. IEEE Transactions Sonics and Ultrasonics. 1981. 318-329 p.

23 Ibid., p 318-329.

6. RESULTADOS Inicialmente se realizó la evaluación del solucionador mecánico por medio de un problema de carga estática. Fue analizado el problema clásico de una viga en voladizo sometida a dos condiciones diferentes de carga, axial y transversa. Después fue analizada la respuesta en frecuencia de dos anillos y un disco piezoeléctrico (códigos del fabricante: APC-181, APC-613 y APC-1029), mostrando los resultados obtenidos para la impedancia eléctrica y los modos de vibración. La impedancia eléctrica fue comparada con los resultados experimentales. Posteriormente se analiza un sonotrodo, que consiste en un actuador de potencia que contiene, además de la cerámica piezoeléctrica, materiales pasivos en su construcción. y se muestran los resultados para este elemento, que combina elementos activos y pasivos. Posteriormente se analiza cómo varía la solución según cambia el tamaño de la malla, para esto se analiza una viga piezoeléctrica en voladizo. Finalizando se hace una comparación con los resultados obtenidos con un software comercial (Comsol) y un modelo analítico unidimensional, mostrando aceptable concordancia. 6.1 VERIFICACIÓN DEL SOLUCIONADOR DEL PROBLEMA MECÁNICO En la práctica, la mayor parte del código de elementos finitos se dedica a la solución del problema elástico, pues debido al análisis armónico realizado, la parte eléctrica es pequeña. En la Ecuación de movimiento (7) y (8) se puede ver que la parte eléctrica se reduce al término fuente que modela la excitación y consiste en una matriz de valor constante. Por tanto, se puede decir que la base del código es la parte mecánica, siendo necesario corroborar que proporcione resultados coherentes. La validación de la parte mecánica del código se realizó abordando un problema con solución analítica conocida. Se trata del análisis de una viga esbelta de sección constante a dos condiciones diferentes de carga: axial y transversal. La viga consiste en un elemento de sección rectangular de 0,15x0,06 m² y una longitud de 1,8 m, fabricado de acero estructural A36, empotrado en un extremo y cargado en el otro. Fue aplicada una carga de 5,2 kN en las dirección axial a tensión y transversalmente. En ambos casos la carga fue distribuida en el área transversal.

Para el caso de carga axial, el desplazamiento en extremo se calcula como:

AEPL=δz (31)

y el esfuerzo normal en un plano transversal en la mitad de la viga como:

AP=σ zz (32)

Donde P es la carga aplicada, L y A son la longitud y el área transversal de la viga y E es el módulo de elasticidad del material. Para el caso de carga transversal, el desplazamiento en el extremo de la viga se calcula como:

EIPL=δ

3

x 3 (33)

Y el esfuerzo normal máximo es:

IMy=σ zz

(34)

Donde I es el momento de inercia de la sección transversal, M es el momento flector máximo que se presenta en el empotramiento y (y) es la distancia desde el eje neutro hasta la superficie de la viga. La Figura 3 muestra el desplazamiento y el esfuerzo normal para el caso de carga axial. Se puede ver que la viga se deforma de la manera esperada y se presenta un esfuerzo en la dirección Szzz uniforme, con excepción en las esquinas de los extremos. Este esfuerzo mayor en las esquinas es resultado de la concentración de esfuerzos, ocasionada por la menor área de material en estas zonas. Este resultado es el esperado.

La Figura 4 muestra el desplazamiento y el esfuerzo normal para el caso de carga transversal. Se puede ver a deflexión esperada de la viga en la Tabla 3. La distribución del esfuerzo también es coherente, pues se pueden ver los valores máximos en el empotramiento, que es el sitio con mayor momento flector, además de presentarse a tensión en la parte superior y a compresión en la inferior. La Tabla 3 muestra los resultados numéricos y analíticos. Se puede ver que el error porcentual es menor al 2,16% en todos los casos. En el caso axial los resultados deben ser exactos, como efectivamente se obtiene en el caso del esfuerzo. En el caso del desplazamiento se muestra el máximo desplazamiento, que ocurre en las esquinas del extremo cargado. Al inspeccionar el centro del área cargada, se pudo ver que el valor es más exacto. Por tanto, la diferencia se presenta por el carácter tridimensional del análisis. Figura 3. Desplazamiento (arriba) y esfuerzo normal (abajo) de la viga sometida a carga axial.

Figura 4. Desplazamiento (arriba) y esfuerzo normal (abajo) de la viga sometida a carga transversal.

Tabla 3. Comparación de los resultados analíticos y por elementos finitos del problema lineal elástico analizado.

Solución analítica

Elementos finitos

Error porcentual

Axial zδ 5,07x10-6 m 5,15x10-6 m 1,58%

Axial σ zz 5,778x105 Pa 5,778x105 Pa 0,00%

Transversal xδ 0,0029 m 0,0029 m 0,00%

Transversal zzσ

41,6x106 Pa 42,5x106 Pa 2,16%

Los resultados mostrados corroboran que la parte mecánica del código desarrollado proporciona resultados coherentes. 6.2 RESPUESTA EN FRECUENCIA DE LOS ELEMENTOS PIEZOELÉCTRICOS Se usaron tres elementos piezoeléctricos para permitir la comparación entre los resultados experimentales y los obtenidos por medio del código de Matlab. Fueron seleccionadas tres geometrías diferentes. El procedimiento experimental, se menciona en el Capítulo 3. Donde se describe como se obtienen las curvas de impedancia eléctrica. Estos resultados se comparan con la respuesta del código por elementos finitos. Los elementos piezoeléctricos son mostrados en la Figura 5, y corresponden a las referencias APC-181, APC-613 y APN-1029 del fabricante American Piezo. Las dimensiones se muestran en la Tabla 4. Figura 5. Cerámicas APC de prueba.

Tabla 4. Dimensiones de los elementos piezoeléctricos analizados

Ítem Radio externo (mm)

Radio interno (mm)

Espesor (mm)

181 30,00 5,25 5,00 613 25,40 12,60 6,40

1029 10,00 No tiene 2,00 6.2.1 APC-181. El procedimiento para obtener una respuesta del código desarrollado, comienza con la creación de la malla, que, para todos los casos, fue generada por medio del programa Salome. En la Figura 6 se puede observar la malla para el ítem 181. A continuación, la malla se exporta en formato de texto,

donde se encuentran las matrices de las coordenadas de los nodos y de conectividad. Haciendo modificaciones simples en el archivo de texto, se logra que sea leída en Matlab. El archivo de malla también proporciona los valores del número de nodos y elementos. Las condiciones de frontera son aplicadas en grupos de nodos que son creados en Salome y también exportados a un archivo de texto. Las matrices de propiedades de los materiales se incluyen también en un archivo de texto. El código genera las matrices de las Ecuaciones (9) a partir de la malla y las matrices de material, usando cuadratura de Gauss de 3 puntos para la integración. Las matrices de cada elemento luego son ensambladas y aplicadas las restricciones, en caso de que las hubiera. Finalmente, se hace el análisis armónico en la banda de frecuencia de interés, guardando en cada paso frecuencia el vector de solución, que contiene los desplazamientos complejos. Figura 6. Malla usada para el análisis del elemento APC-181 (3250 nodos y 2520 hexaedros).

La Figura 7 muestra la curva de impedancia eléctrica calculada en el rango entre 10 kHz y 2 MHz, con pasos de 1 kHz. Se muestran la magnitud y fase de la impedancia eléctrica. Se evidencia la presencia de múltiples resonancias, donde el valor de a impedancia se reduce considerablemente, por tanto, se transfiere más energía eléctrica. Las primeras resonancias corresponden a los modos radial y de pared (thickness wall) y cerca a los 1,4 MHz se presenta el pico de menor impedancia, que corresponde a un modo espesor. Estos resultados concuerdan con los reportados en la literatura para un elemento piezoeléctrico con características similares24.

24 Ibid., p 318-329.

Figura 7. Impedancia eléctrica (ohm) en función de la frecuencia para el elemento APC-181.

En la Figura 7, los puntos de menor impedancia en la respuesta del código, representan los modos de vibración. Cada modo de vibración corresponde a una manera específica de vibración de la estructura. Los modos pueden ser animados haciéndolos variar sinusoidalemnte, usando las amplitudes y fases correspondientes a cada grado de libertad. Para este propósito se creó un script en Matlab que lee el vector desplazamiento a una determinada frecuencia y genera una serie de vectores de desplazamiento que cubren un ciclo de vibración. Estos vectores son guardados en formato VTK para ser leídos en Paraview y poder visualizar y animar el modo de vibración. En la Figura 8 se muestra los primeros 10 modos, es decir, la vibración ocurrida en los 10 puntos de menor impedancia eléctrica. El modo 1 es el radial fundamental, los modos 2 y 3 son los modos de pared, pues se nota que la superficie de la pared se mueve alrededor de un nodo de desplazamiento casi nulo. Los modos 5, 6 y 7 son armónicos de los modos de los modos anteriores. El modo 10 es el fundamental, pues presenta la menos impedancia eléctrica, y es espesor. Es decir que los mayores desplazamientos se presentan en la dirección del espesor. Si bien los gráficos proporcionan una idea del modo de vibrar, lo más recomendable es implementar la animación, que muestra con más detalle y sin ambigüedades de interpretación el movimiento.

Figura 8. Principales modos de vibración observada en el elemento APC-181.

Modo 1 Modo 2

Modo 3 Modo 4

Modo 5 Modo 6

Modo 7 Modo 8

Modo 9 Modo 10 La Figura 9 muestra la impedancia eléctrica simulada y experimental del anillo piezoeléctrico analizados (APC-181). Teniendo en cuenta las limitaciones del sistema de medición, se puede concluir que las curvas de la Figura 9, muestran buena concordancia, pues la cantidad y posición de los picos es muy similar. En ambos casos se pueden ver los modos radiales y de pared esperados. En el caso del elemento APC-181 se puede observar un pico cerca de 1,4 MHz, que como se mostró anteriormente se trata de un modo espesor. Este modo se observó inicialmente por elementos finitos y era inesperado. Para corroborar su existencia se realizó la caracterización experimental hasta 2 MHz, siendo encontrado a una frecuencia cercana. Este modo de vibración es muy interesante, porque podría

permitir la fabricación de transductores de alta frecuencia sin tener los espesores tan pequeños requeridos con geometría tipo disco. Figura 9. Impedancia eléctrica en función de la frecuencia calculada por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-181.

6.2.2 APC 613. El elemento APC-613 fue analizado en un rango de frecuencia entre 10 kHz y 1 MHz El resultado obtenido después del procesamiento del código, arroja como resultado la Figura 10, en este caso, se observa que el punto de mínima impedancia se encuentra en la frecuencia de 1000 Khz, este el modo de resonancia. Los otros puntos de mínima impedancia antes de su frecuencia de resonancia, son los modos de vibración, se evalúa de nuevo el código, haciendo el análisis para los 7 modos de mínima impedancia, obteniendo la Figura 11, donde se observa como los modos 2 y 3, son los modos de pared, mientras que los modos 6 y 7 son los armónicos de estos primeros, el modo 8 es el modo de resonancia que se usa para el trabajo de esta geometría.

La comparación de la respuesta en frecuencia, obtenida con el código de elementos finitos y el montaje experimental de una geometría, se observa en la Figura 12, en donde se observa que los puntos de mínima impedancia en la prueba experimental, se encuentra en una frecuencia de Khz=f r 1097 , mientras que con el código de elementos finitos, nos indica que la frecuencia de resonancia tiene un valor de

Khz=f rfe 983 teniendo una diferencia de 114=ff=f rferdesv kHz, esto indica que el error relativo del procedimiento tiene un valor de 10.39=erel , es importante tener en cuenta que las propiedades físicas no son aproximadas, ocasionando la desviación entre los valores reales y los valores del código. Los modos de vibración restantes se presentan una desviación importante. Figura 10. Impedancia eléctrica (ohm) en función de la frecuencia para el elemento APC-613.

Figura 11. Principales modos de vibración observada en el elemento APC-181.

Modo 1 Modo 2

Modo 3 Modo 4

Modo 5 Modo 6

Modo 7 Modo 8 La Figura 12 muestra la impedancia eléctrica simulada y experimental del anillo piezoeléctrico (APC-613). Como en el caso del (APC 181), se observa la buena concordancia de las curvas por la cantidad y posición de los picos, siendo similares entre sí. En ambos casos se pueden ver los modos radiales y de pared esperados. También en el caso del elemento APC-613 se observa claramente cómo el problema de la sensibilidad en el sistema de medición afecta la forma de la curva de impedancia, ocasionando en muchos casos que los picos más pequeños no puedan ser visualizados.

Figura 12. Impedancia eléctrica en función de la frecuencia calculada por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-613.

6.2.3 APC-1029. La curva de impedancia eléctrica del elemento APC-1029 (disco) se muestra en la Figura 13. De ella se obtienen siete puntos de mínima impedancia, correspondientes a siete modos de vibración. El modo fundamental se encuentra a

1026=f kHz, que corresponde al modo espesor. Este resultado es el esperado, tal como se reporta en la literatura25. También se pueden observar los modos de pared y radiales. La Figura 14 muestra gráficamente los siete modos de vibración más importantes obtenidos. Se pueden ver claramente los modos radiales (1 y 2), los modos de pared (3 y 4) y el modo espesor fundamental (7). Los modos de vibración 5 y 6 son modos acoplados, es decir, son la sobreposición de dos o más modos. En este caso se debe tratar de la sobreposición de un modo espesor con modos de pared, pues se puede ver que incluyen características de ambos.

25 Ibid., p 318-329.

Figura 13. Impedancia eléctrica (ohm) en función de la frecuencia para el elemento APC-1029.

La Figura 15 compara los resultados entre la respuesta obtenida por el código de elementos finitos, y el montaje experimental, se observa una desviación de 9,6=erelpara las frecuencias de resonancia del modo fundamental. Este error es bajo, teniendo en cuenta la incertidumbre en las propiedades físicas del material. En general, se pueden ver similitudes claras en las curvas, como el número de picos de mínima amplitud y su tamaño relativo, siendo observados claramente los modos radiales y el de espesor. Sin embargo, la frecuencia exacta a la cual se presentan los picos es diferente. Estas diferencias deben ser causadas por el uso de unas propiedades de material genéricas y limitaciones del sistema de medición. A pesar de las diferencias, se puede concluir que el análisis por elementos finitos reproduce en buena medida el comportamiento en frecuencia observado experimentalmente.

Figura 14. Principales modos de vibración observada en el elemento APC-1029. Modo 1

Modo 2

Modo 3

Modo 4

Modo 5

Modo 6

Modo 7 El análisis anterior es importante en el diseño de muchos dispositivos piezoeléctricos. Por ejemplo, en el caso de ensayos no destructivos por ultrasonido, donde se usan transductores que generan ondas longitudinales (ondas de presión), un disco cerámico vibrando en el modo espesores el elemento activo que genera y recibe las ondas acústicas. Para esta aplicación el disco cerámico es excitado con un pulso de alto voltaje y un ancho variable, de tal manera que la frecuencia fundamental del espectro de excitación sea la misma frecuencia del modo espesores de la cerámica. Si las frecuencias no coinciden, la respuesta del transductor puede ser muy pobre.

Figura 15. Impedancia eléctrica en función de la frecuencia calculada por elementos finitos (arriba) y medida experimentalmente (abajo) para el elemento APC-1029.

En el análisis fueron usadas propiedades genéricas del material piezoeléctrico porque el fabricante (American Piezo) solamente reporta algunas de las propiedades elásticas. El cálculo por elementos finitos requiere la matriz de rigidez completa, que para una material con simetría hexagonal son 5 constantes. El fabricante solamente reporta 2 constantes, los módulos de Young en la dirección del espesor ( 33Y ) y en la dirección transversal ( 11Y ). Por tanto, se usaron las propiedades genéricas para un material piezoeléctrico usado en aplicaciones de potencia (PZT-5H) halladas en la literatura de ANDRADE. Por lo general, la obtención de las propiedades elásticas exactas de un material piezoeléctrico implica un proceso de caracterización elaborado que requiere equipos de laboratorio y el

ajuste de los datos experimentales a modelos teóricos usando algoritmos de optimización. El sistema de medición implementado en el laboratorio está en proceso de desarrollo y aún tiene limitaciones que impiden la obtención exacta de la curva de impedancia eléctrica. Actualmente, el sistema puede levantar la curva, identificando con precisión la localización (frecuencia) y tamaño relativo de los picos, pero los valores reales de la amplitud aún no pueden ser obtenidos, en especial los de mínima impedancia. Se trata de un problema de software que afecta la sensibilidad del sistema, pues los valores de impedancia muy baja no son escasamente detectados, y será solucionado en el futuro. 6.3 ANÁLISIS DE UN ACTUADOR TIPO LANGEVIN Un actuador tipo Langevin es un dispositivo piezoeléctrico donde los elementos activos, generalmente anillos de cerámica tipo PZT, son presionados entre dos cilindros de material pasivo (no piezoeléctrico). Dependiendo de la longitud de los elementos pasivo, el actuador puede ser sintonizado en una frecuencia determinada de trabajo. Estos dispositivos tienen muchas aplicaciones industriales y son los más usados en ultrasonido de potencia. La Figura 16, muestra el actuador analizado. Consiste en un amplificador mecánico de tipo escalonado con dos segmentos de igual longitud, pero diferente área. Al tener la misma longitud, los dos segmentos presentan el modo flextensional a la misma frecuencia, pero la diferencia de área causa que el desplazamiento sea mayor en el segmento de menor área, que en este caso cumple el propósito de herramienta. Ambos segmentos son resonadores de media longitud de onda ( 2/λ ). El segmento de mayor área contiene las cerámicas piezoeléctricas, tipo PZT-5H, en la mitad de su longitud. Actuadores de este tipo son usados en aplicaciones como atomización ultrasónica, aceleración de reacciones químicas, entre otras.

Figura 16. Aspecto del actuador tipo Langevin analizado.

El principio de funcionamiento del resonador mostrado es el siguiente: al aplicar un voltaje a la cerámica, esta se deforma generando una onda de presión o tensión mecánica que se propaga en el material pasivo. La onda se refleja al final del material pasivo y regresa a la cerámica. Existe una frecuencia determinada a la cual la onda regresa en fase con el desplazamiento de la cerámica y se produce interferencia constructiva. Esta frecuencia corresponde al modo de vibración flextensional, donde el sonotrodo vibra aumentando y reduciendo su longitud, presentándose un alto desplazamiento. La frecuencia depende de la velocidad de la onda longitudinal en el material, por tanto, la longitud del material pasivo depende del tipo de material utilizado. La tabla 5 muestra las propiedades más importantes del actuador. Tabla 5. Propiedades del actuador analizado.

Material piezoeléctrico PZT-5H (geometría del ítem 613)

Material pasivo Aluminio

Dimensiones Longitud 122,4 mm (61,2 + 61,2 mm) Radio mayor 12,7 mm Radio menor 3,175 mm

Frecuencia de diseño 40 kHz

El análisis por elementos finitos fue realizado usando el código desarrollado. Sin embargo, algunos cambios adicionales debieron ser realizados. El material pasivo fue modelado como un material piezoeléctrico con constantes piezoeléctricas y dieléctricas nulas, y las constantes elásticas correspondientes a un material isotrópico. Esto se hizo generando grupos de elementos que permiten aplicar materiales diferentes en la malla. Finalmente se realizó un análisis armónico en el rango entre 10 y 80 kHz. Figura 17. Impedancia eléctrica del actuador tipo Langevin calculada por elementos finitos.

La Figura 17 muestra la curva de impedancia eléctrica en función de la frecuencia en el rango entre 10 y 80 kHz. Se puede claramente una resonancia alrededor de los 40 kHz, exactamente la frecuencia de diseño. Sin embargo, hay un detalle que no se ha aclarado completamente, el cual tiene que ver con la forma de la curva. Generalmente se tiene primero el pico de resonancia (mínima impedancia eléctrica)

y después el de antiresonancia (máxima impedancia eléctrica). El resultado muestra el comportamiento contrario. Figura 18. Modos de vibración observados a las frecuencias de antiresonancia (izquierda) y resonancia (derecha) del actuador tipo langevin.

Para evaluar la manera a la que vibra el sonotrodo a la frecuencia de diseño se hizo el gráfico de los modos de vibración. La Figura 18 muestra los modos de vibrar en las frecuencias de antiresonancia (39 kHz) y resonancia (40,5 kHz). Se puede ver que en la antiresonancia se tiene un modo que no es flextansional, con desplazamiento lateral. El modo observado en la frecuencia de resonancia es claramente flextensional, pues se puede ver cómo la longitud del sonotrodo se alarga y se acorta, con desplazamiento laterales pequeños, que no se notan a la escala del gráfico. También se puede ver que las amplitudes de los desplazamientos son mayores en el caso de la resonancia, a pesar de que fue usado un factor de escala para los desplazamientos 10 veces menor. 6.4 ANÁLISIS DE LA DEPENDENCIA DE LA MALLA En toda simulación por elementos finitos debe hacerse un análisis que permita asegurar que el resultado no depende del tamaño de los elementos finitos. Este

análisis se denomina dependencia de la malla y consiste en solucionar el problema para mallas cada vez más refinada, hasta obtener un resultado constante. Cuando se trabaja con elementos finitos con interpolación de primera orden, el análisis de dependencia de la malla es imprescindible para asegurar la calidad de los resultados. Para analizar cómo el tamaño de la malla afecta el resultado, fue calculada la respuesta en una viga en cantiliver de material piezoeléctrico. La Figura 19 muestra un esquema del problema. Los electrodos donde se aplica el voltaje corresponden a las superficies superior e inferior, aquellas de mayor área. Cuando el espesor de la viga es pequeño comparado con el ancho, esta configuración presenta un modo de vibración característico, con desplazamientos principalmente en el plano de los electrodos. Figura 19. Viga en cantilever de material piezoeléctrico.

La viga fue modelada aplicando una restricción de desplazamiento cero en la cara izquierda. Como se trata de un problema dinámico, además de modificar la matriz de rigidez, también se deben modificar la matriz de masa y el vector de cargas para modelar el empotramiento de las vigas. Fue realizado un análisis armónico en un rango de frecuencia entre 10 kHz y 1 MHz, con pasos de 1 kHz. Este procedimiento fue realizado para tres mallas con diferente número de elementos. La Figura 20 muestra las curvas de impedancia eléctrica para las tres mallas usadas. Se puede ver claramente una resonancia (modo fundamental), donde la impedancia eléctrica es mínima. Para la malla 1 (la menos refinada) se obtuvo una frecuencia de resonancia mayor que para las otras dos, y a medida que se refina la malla, la frecuencia disminuye, mostrando un claro comportamiento asintótico. Queda claro que la malla menos refinada introduce un error elevado en el cálculo

de la frecuencia. La Tabla 6 muestra el tamaño de las mallas y la frecuencia del modo fundamental. Figura 20. Impedancia eléctrica en función de la frecuencia de la viga en cantiliver para los tres diferentes tamaños de malla.

Tabla 6. Frecuencia del modo fundamental para tres mallas con diferente número de elementos.

Número de elementos

Número de grados de

libertad

Frecuencia del modo

fundamental

Malla 1 52 208 486 kHz

Malla 2 400 1600 452 kHz

Malla 3 3200 12800 442 kHz Las diferencias también se manifiestan en la forma como la estructura vibra a la frecuencia de resonancia. La Figura 21 esquematiza el modo de vibración para las tres mallas usadas. Si bien la vibración es similar, pues las mayores amplitudes de

desplazamiento se presentan en el plano de los electrodos, hay diferencia en los desplazamientos laterales. Sin embargo, en la práctica estas pequeñas diferencias pueden no ser relevantes. Las amplitudes mostradas en las barras de color (colorbar) son diferentes porque no se muestran en el mismo instante de tiempo. El máximo desplazamiento es casi el mismo en los tres casos. Figura 21. Modo de vibración fundamental de la viga empotrada obtenido con las diferentes mallas.

El modo de vibración mostrado, con desplazamientos de mayor amplitud en el plano del electrodo, puede tener aplicaciones en la medición de viscosidad de líquidos. Esta estructura piezoeléctrica inmersa puede generar cizallamiento en un líquido, ocurriendo una pérdida de energía que varía con la viscosidad. Analizando la curva de resonancia de los casos inmerso en aire e inmerso en el líquido, puede determinarse la viscosidad. La determinación de la viscosidad, y posiblemente también la elasticidad (viscoelasticidad), puede realizarse por el ajuste de un modelo

o por calibración del dispositivo. Esta es una aplicación interesante que vale la pena evaluar. 6.5 COMPARACIÓN CON EL SOFTWARE COMERCIAL COMSOL Fue realizada una comparación entre los resultados obtenidos con el código desarrollado en este trabajo, el software comercial Comsol y un modelo analítico. En este análisis se usó una disco piezoeléctrico de 50,6 mm de diámetro y 2,8 mm de espesor (American Piezo, referencia APC-612), vibrando en el vacío, es decir, sin ningún tipo de restricción de los desplazamientos. Fueron usadas las mismas propiedades físicas reportadas anteriormente. Para este caso, donde el diámetro es mucho mayor que el espesor, existe una solución analítica que modela la vibración en la dirección del espesor (unidimensional). Este modelo calcula la impedancia eléctrica de la siguiente manera26:

2

2tan

11

33

233

33

33

kl

kl

εe

+c

cjωωAl=Z

SE

E

S33

e (34)

Donde l y A son el espesor y el área del disco, ω es la frecuencia angular, c33

E,

e33 y ε33S

son las constantes elástica, piezoeléctrica y dieléctrica, respectivamente, y:

S

E

c

εe

+c

ρω=k

33

233

33

(35)

26 NEGREIRA Carlos, Adamowski Julio C, Buiochi Flavio, Andrade Marco A. Analysis of 1-3 Piezocomposite and Homogeneous Piezoelectric Rings for Power Ultrasonic Transducers, University of Sao Paulo – USP, Escola Politecnica, Sao Paulo, 2009. 1-7 p.

Donde ρc es la densidad de la cerámica. La Figura 22 muestra la impedancia eléctrica obtenida de las tres maneras: solución analítica, software Comsol y código desarrollado en este trabajo. La solución analítica muestra claramente el pico del modo espesor, que también es obtenido por elementos finitos. En el modelo de Comsol no fue incluido amortiguamiento, por tal razón se presenta una curva tan poluida, llena de múltiples resonancias. Este resultado es correcto numéricamente, sin embargo, en la realidad estos modos son atenuados por el amortiguamiento del material. El código desarrollado tiene implementado amortiguamiento de Rayleigh, por tanto, se muestra una curva más suave, donde se puede ver que el modo espesor está muy cercano al proporcionado por el modelo analítico. En el caso de Comsol, el modo espesor no se visualiza correctamente. Las primeras resonancias, correspondientes a los modos radiales, se pueden observar en los resultados por elementos finitos, mostrando una buena concordancia entre los resultados de este trabajo y Comsol. Figura 22. Comparación de la impedancia eléctrica en función de la frecuencia para un disco piezoeléctrico obtenida por diferentes métodos.

7 CONCLUSIONES Se implementó un código en Matlab para el análisis por elementos finitos de la respuesta en frecuencia de estructuras piezoeléctricas. Fueron usados elementos tipo hexaedro con 8 nodos e interpolación lineal, para un total de 32 grados de libertad (24 desplazamientos y 8 voltajes). El código fue diseñado para realizar el análisis armónico de una estructura piezoeléctrica, incluyendo materiales pasivos (no piezoeléctricos). Sin embargo, permite otro tipo de análisis, tales como estático, transiente y modal, con aplicaciones en resistencia de materiales o acústica. Las mallas fueron elaboradas en Salome, que es un software de código abierto (licencia GPL) desarrollado para el pre y post procesamiento en cálculos de ingeniería. Si bien Salome trabaja muy bien con mallas de tetraedros, el enmallado usando hexaedros es poco automatizado, siendo difícil generar mallas de geometría complicada. El post procesamiento fue realizado en Paraview, que también se trata de un software libre (licencia BSD) para la visualización de datos de ingeniería. Paraview está diseñado para trabajar con ficheros de datos muy grandes y posee una gran cantidad de tipos de gráficos y opciones para la manipulación y edición de gráficos científicos. El código desarrollado en este trabajo exporta los resultados al formato VTK (Visualization Toolkit) para ser leído por Paraview. Matlab es el único software no libre que fue usado en este trabajo. La principal ventaja de Matlab es la facilidad del modelado matemático, lo que permite el rápido desarrollo de los métodos numéricos involucrados. Sin embargo, la licencia académica disponible solo puede ser usada para fines académicos y de investigación. Para uso comercial se requiere una licencia costosa. Por esta razón, un trabajo futuro consiste en implementar el código en Python, o una mezcla de Python y C para eliminar la dependencia de software propietario, reducir los tiempos de procesamiento y aumentar el tamaño de los modelos. Aunque el código permite resolver miles de elementos, modelos más elaborados requieren un mejor desempeño, que no puede lograrse con matlab. La parte más importante del código es el solucionador del problema elástico (parte mecánica). Por esta razón fue analizado un problema de simple de resistencia de materiales, consistente en una viga prismática empotrada en un extremo y cargada en el extremo otro. Fueron analizados dos casos: carga axial y carga transversal. Los resultados mostraron buena concordancia con la solución analítica conocida, con desviaciones en los desplazamientos con un error porcentual de 1,58%.

También los valores de esfuerzo calculados mostraron concordancia con la teoría. Estos resultados mostraron que la parte mecánica del código funciona correctamente. La parte del acoplamiento electro-mecánico es relativamente simple. Como la excitación es un voltaje sinusoidal y no hay cargas mecánicas aplicadas, el término de acoplamiento electro-mecánico de la ecuación de movimiento se convierte en un término fuente constante. El análisis armónico se realizó entonces suponiendo desplazamientos también sinusoidales, expresados como números complejos. Los desplazamientos fueron almacenados para cada frecuencia y los modos de vibración animados en Paraview a partir de los desplazamientos en formato de magnitud y fase. Un disco y dos anillos de material piezoeléctrico fueron analizados para determinar su respuesta en frecuencia. Fue realizado el análisis armónico en una banda de frecuencia lo suficientemente ancha para capturar el modo fundamental y otros modos adicionales. Las curvas de impedancia eléctrica mostraron una serie de resonancias, cada una relacionada con un modo de diferente de vibración, tal como era esperado. La animación en Paraview de la respuesta armónica de los desplazamientos permitió identificar los modos de vibración correspondientes a cada resonancia. Fue posible identificar los modos radiales, de pared y espesor de los elementos piezoeléctricos, similares a mostrados en la literatura. Las mediciones experimentales realizadas en el laboratorio mostraron buena concordancia con las curvas obtenidas por elementos finitos. Las diferencias son explicadas por las propiedades físicas de los materiales, las cuales no son conocidas con precisión. Esto se debe a que el fabricante (APC) proporciona sólo algunos valores, los necesarios para el cálculo analítico de los dispositivos más comunes. Para la simulación por elementos finitos se requieren todas las constantes elásticos, piezoeléctricas y dieléctricas. Estos valores fueron obtenidos de la literatura, donde se reportan para materiales genéricos. Por otro lado, el sistema de caracterización debe ser mejorado, pues aún presenta problemas en el manejo de la escala vertical. En definitiva, las curvas simuladas y medidas tienen razonable concordancia con relación a la posición (frecuencia) de las resonancias, lo que en la práctica es la información más importante. Se sabe que la interpolación de primera orden induce errores grandes en determinados análisis, dependiendo principalmente de la geometría. Por esta razón se realizó un ensayo de dependencia de la malla. Fue analizada la respuesta vibratoria de una viga en cantilever de PZT con diferente discretización. Resultó claro que mallas poco refinadas llevan a errores grandes en la frecuencia

fundamental, además de ligeros cambios en el modo de vibrar del elemento. Es decir, se requieren mallas grandes para obtener resultados precisos. Esto muestra la necesidad de implementar elementos con interpolación cuadrática. Fue comparada la respuesta en frecuencia de un disco piezoeléctrico de PZT calculado por medio de un modelo analítico clásico, el software comercial de MEF Comsol y el código desarrollado en este trabajo. El resultado mostró buena concordancia en la determinación de la frecuencia del modo espesor fundamental. También se mostró concordancia entre los resultados obtenidos con Comsol y el código desarrollado para los modos radiales. Sin embargo, para una mejor concordancia se debe implementar el amortiguamiento en el modelo de Comsol. En conclusión, el código desarrollado en este trabajo proporciona valores coherentes de las frecuencias de resonancia y la forma de los modos de vibración. Por tanto, puede ser usado como una herramienta de análisis en el estudio o diseño de dispositivos piezoeléctricos. 7.1 TRABAJOS FUTUROS Para el diseño y análisis de dispositivos piezoeléctricos, las siguientes mejoras deben ser implementadas en el código (en orden de importancia): Incluir el elemento tipo tetraedro para una mejor representación de geometrías complejas. Aumentar el orden de interpolación para obtener mejor precisión. Mejorar el desempeño programando el núcleo de cálculo en C o C++ y técnicas de procesamiento en paralelo, como el uso de la librería Openmp, y de esta manera aumentar el número de elementos del modelo a cientos de miles. Crear una interfaz gráfica que le facilite el trabajo al usuario.

BIBLIOGRAFÌA ANDRADE Marco, Pérez Nicolas, Buiochi Flavio, Adamowski Julio. Matrix method for acoustic levitation simulation. IEEE Transactions Ultrasonics Ferroelectrics and Frequency Control. 2011. 1674-1683 p. AVDIAJ Sefer, Setina Janez, Syla Naim, Modeling of the piezoelectric effect using the finite – element method, University of Prishtina, Faculty of Natural Science and Mathematics, Prishtina, Kosovo, 2009. 1-9 p. AULD B. A. Acoustic Fields and Waves in Solids, Volumes I and II. Vol 8. Florida: Krieger Publishing Company, 1975. 405 p. BATHE K.-J. Finite Element Procedures. New Jersey: Prentice Hall,1996. 1037 p. BOUCHER Didier, Lagier Michel, Maerfeld C. Computation of the Vibrational Modes for Piezoelectric Array Transducers using a Mixed Finite Element-Perturbation Method. IEEE Transactions Sonics and Ultrasonics. 1981. 318-329 p. CURRIE J, Currie P. Développement, par pression, de electricité polaire dans les cristaux hémièdres à faces inclinées. Comptes rendus. Académie des Science. 1880. 294-295 p. CHEEKE J. David, Zagzebski J. Fundamentals and Applications of Ultrasonic Waves. Physics Department Concordia University. Montreal: CRC Press LLC 2004. 72-719 p. FRANCO Edinguer E, Cuervo Andrés M, Barrera Helver M. Desarrollo de un programa por elementos finitos para analizar la respuesta vibratoria de estructuras piezoeléctricas. Facultad de Ingeniería, Universidad Autónoma de Occidente, Cali, 2016. 1-14 p. KINO Gordon S. Acoustic Waves: Devices, Imaging, and Analog Signal Processing. Mèxico: Prentice Hall, 1987. 531 p

NEGREIRA Carlos, Adamowski Julio C, Buiochi Flavio, Andrade Marco A. Analysis of 1-3 Piezocomposite and Homogeneous Piezoelectric Rings for Power Ultrasonic Transducers, University of Sao Paulo – USP, Escola Politecnica, Sao Paulo, 2009. 1-7 p. PIEFORT Vincent, Preumont André. Finite Element Modelling of Piezoelectric Active Structures. Active Structures Laboratory, ULB-CP 165/42, Brussels,2001. 1-17 p. SHIARI Behrouz, Wild Peter M. Coupled finite element analysis of a piezo-ceramic force transducer. Deparment of Mechanical & Aerospace Engineering, Carlenton University, Ottawa, Canada, 2006. 167-181 p. ERHART Mgr Jiri. Piezoelectricity and ferroelectricity Phenomena and properties. Department of Physics FP TUL. Presentación de clase, año no disponible, 1-65 p.