Modelización y Análisis de Vigas Piezoeléctricas para la Captación de Energía

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería Civil Modelización y Análisis de Vigas Piezoeléctricas para la Captación de Energía Autor: María Dolores Lorenzo Espejo Tutor: Luis Rodríguez de Tembleque Solano Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 iii Trabajo Fin de Grado Ingeniería Civil Modelización y Análisis de Vigas Piezoeléctricas para la Captación de Energía Autora: María Dolores Lorenzo Espejo Tutor: Luis Rodríguez de Tembleque Solano Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023 v Proyecto Fin de Carrera: Modelización y Análisis de Vigas Piezoeléctricas para la Captación de Energía Autor: María Dolores Lorenzo Espejo Tutor: Luis Rodríguez de Tembleque Solano El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2023 El Secretario del Tribunal vii A mis padres y a mis hermanos. ix Agradecimientos Con este trabajo finaliza mi etapa de estudiante de grado. Estoy profundamente agradecida a estos años por haberme permitido encontrar mi pasión en la ingeniería civil y por el crecimiento académico y personal que me han aportado. Llegar hasta aquí ha sido posible gracias a la ayuda de mi tutor Luis, a quien agradezco su colaboración y su implicación desde principio a fin de la elaboración de este trabajo. También quiero agradecer a mis amigos y a mi familia la confianza que han depositado en mí durante todo este proceso. Gracias por apoyarme y por guiarme a lo largo de este camino. María Dolores Lorenzo Espejo Sevilla, 2023 ÍNDICE DE FIGURAS Figura 1 (a): Efecto piezoeléctrico directo 2 Figura 1 (b): Efecto piezoeléctrico inverso 2 Figura 2: Técnica de aumento de rango de excitación del recolector piezoeléctrico 3 Figura 3: Empleo de sensores piezoeléctricos en aerogeneradores 3 Figura 4: Morfologías de los recolectores unimorfo y bimorfos. 7 Figura 5: Dispositivos piezoeléctricos de modos 33 y 31. 8 Figura 6: Geometría del recolector piezoeléctrico a estudiar. 9 Figura 7: Alzado de la viga piezoeléctrica. 10 Figura 8: Geometría de los elementos SOLID186 y SOLID226. 14 Figura 9: Viga piezoeléctrica 3D mallada. 16 Figura 10: Aplicación de la carga en el extremo libre. 16 Figura 11: Geometría de los elementos PLANE183 y PLANE223. 17 Figura 12: Modelo de la viga piezoeléctrica en 2D 19 Figura 13: Secuencia de carga escalón aplicada 20 Figura 14: Influencia del parámetro β en el voltaje producido según α 25 Figura 15 (a): Geometría mallada en modelo 2D 26 Figura 15 (b): Geometría mallada en modelo 3D 26 Figura 16: Relación fuerza-flecha (𝑡/𝐿=0,1, 𝑡/𝑊=4 y β=3) 27 Figura 17: Relación fuerza-voltaje (𝑡/𝐿=0,1, 𝑡/𝑊=4 y β=3) 27 Figura 18: Influencia de β en el voltaje según el modelo empleado (𝑡/𝐿=0,1, 𝑡/𝑊=4) 28 Figura 19 (a): Geometría para t/W=0,1, 𝑡/𝐿=0,1 y β=3 30 Figura 19 (b): Geometría para t/W=4, 𝑡/𝐿=0,1 y β=3 30 Figura 20: Influencia de 𝑡/𝑊en la relación fuerza-flecha 30 Figura 21: Influencia de 𝑡/𝑊en la relación fuerza-voltaje (𝑡/𝐿=0,1 y β=3 ) 31 Figura 22 (a): Deformada para t/W= 0,1, 𝑡/𝐿=0,1 y β=3 31 Figura 22 (b): Deformada para t/W= 4, 𝑡/𝐿=0,1 y β=3 31 Figura 23 (a): Voltaje para t/W= 0,1, 𝑡/𝐿=0,1 y β=3 32 Figura 23 (b): Voltaje para 𝑡/𝑊= 4, 𝑡/𝐿=0,1 y β=3 32 Figura 24 (a): Influencia de efectos de 2º orden en la flecha para t/W=0,1 34 Figura 24 (b): Influencia de efectos de 2º orden en el voltaje para t/W=0,1 34 Figura 24 (c): Influencia de efectos de 2º orden en la flecha para t/W=1 34 Figura 24 (d): Influencia de efectos de 2º orden en el voltaje para t/W=1 34 Figura 24 (e): Influencia de efectos de 2º orden en la flecha para t/W=4 34 Figura 24 (f): Influencia de efectos de 2º orden en el voltaje para t/W=4 34 Figura 25: Influencia de 𝑡/𝐿 en la relación fuerza-flecha (𝑡/𝑊=4 y β=3) 36 Figura 26: Influencia de 𝑡/𝐿 en la relación fuerza-voltaje (𝑡/𝑊=4 y β=3) 36 xvii Figura 27 (a): Influencia de efectos de 2º orden en la flecha para 𝑡/𝐿=0,1 38 Figura 27 (b): Influencia de efectos de 2º orden en la flecha para 𝑡/𝐿=0,05 38 Figura 27 (c): Influencia de efectos de 2º orden en la flecha para 𝑡/𝐿=0,03 38 Figura 27 (d): Influencia de efectos de 2º orden en la flecha para 𝑡/𝐿=0,02 38 Figura 27 (e): Influencia de efectos de 2º orden en la flecha para 𝑡/𝐿=0,01 38 Figura 28 (a): Influencia de efectos de 2º orden en el voltaje para 𝑡/𝐿=0,1 39 Figura 28 (b): Influencia de efectos de 2º orden en el voltaje para 𝑡/𝐿=0,05 39 Figura 28 (c): Influencia de efectos de 2º orden en el voltaje para 𝑡/𝐿=0,03 39 Figura 28 (d): Influencia de efectos de 2º orden en el voltaje para 𝑡/𝐿=0,02 39 Figura 28 (e): Influencia de efectos de 2º orden en el voltaje para 𝑡/𝐿=0,01 39 Figura 29 (a): Deformada para 𝑡/𝐿=0,01, 𝑡/𝑊=4 y β=3 en primer orden 40 Figura 29 (b): Deformada para 𝑡/𝐿=0,01, 𝑡/𝑊=4 y β=3 en segundo orden 40 Figura 30 (a): Deformada para 𝑡/𝐿=0,1, 𝑡/𝑊=4 y β=3 en primer orden 40 Figura 30 (b): Deformada para 𝑡/𝐿=0,1, 𝑡/𝑊=4 y β=3 en segundo orden 40 Figura 31 (a): Voltaje para 𝑡/𝐿 =0,01, 𝑡/𝑊 =4 y β=3 en primer orden 41 Figura 31 (b): Voltaje para 𝑡/𝐿 =0,01, 𝑡/𝑊 =4 y β=3 en segundo orden 41 Figura 32 (a): Voltaje para 𝑡/𝐿 =0,1, 𝑡/𝑊 =4 y β=3 en primer orden 41 Figura 32 (b): Voltaje para 𝑡/𝐿 =0,1, 𝑡/𝑊 =4 y β=3 en segundo orden 41 Figura 33: Influencia de los efectos de 2º orden en la flecha según β (𝑡/𝐿 =0,02 y 𝑡𝑊=0,1) 43 Figura 34: Influencia de los efectos de 2º orden en el voltaje según β (𝑡/𝐿 =0,02 y 𝑡𝑊=0,1) 43 Figura 35: Influencia de los efectos de 2º orden en la flecha según β (𝑡/𝐿 =0,1 y 𝑡𝑊=4) 44 Figura 36: Influencia de los efectos de 2º orden en el voltaje según β (𝑡/𝐿 =0,1 y 𝑡𝑊=4) 45 Figura 37: Influencia de β en el voltaje según 𝑡/𝑊 (𝑡/𝐿 =0,1) 46 Figura 38: Influencia de β en el voltaje según 𝑡𝐿 (𝑡/𝑊=4) 47 Figura 39 (a): Influencia de β en el voltaje según 𝑡/𝐿 para t/W=0,1 en 1er orden 49 Figura 39 (b): Influencia de β en el voltaje según t/L para t/W=0,1 en 2º orden 49 Figura 39 (c): Influencia de β en el voltaje según t/L para t/W=0,8 en 1er orden 49 Figura 39 (d): Influencia de β en el voltaje según t/L para t/W=0,8 en 2º orden 49 Figura 39 (e): Influencia de β en el voltaje según t/L para t/W=4 en 1er orden 49 Figura 39 (f): Influencia de β en el voltaje según t/L para t/W=4 en 2º orden 49 Figura 40: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3, 𝑡/𝑊=4, 𝑡/𝐿 =0,1 y 𝛽𝑅=0,00001) 51 Figura 41: Influencia de los efectos de 2º orden en el voltaje en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,1 y 𝛽𝑅=0,00001) 51 Figura 42: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,03 y 𝛽𝑅=0,00001) 52 Figura 43: Influencia de los efectos de 2º orden en el voltaje en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,03 y 𝛽𝑅=0,00001) 52 Figura 44: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,1 y y 𝛽𝑅=0,01) 54 Figura 45: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,1 y y 𝛽𝑅=0,01) 54 Figura 46: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,03 y 𝛽𝑅=0,01) 55 Figura 47: Influencia de los efectos de 2º orden en el voltaje en análisis transitorio (𝛽=0,3, 𝑡/𝑊 =4, 𝑡/𝐿 =0,03 y 𝛽𝑅=0,01) 55 xix NOTACIÓN 𝐝 Tensor de coeficientes piezoeléctricos d31 coeficiente piezoeléctrico para el modo de carga 31 𝐞 Tensor de coeficientes piezoeléctricos t Espesor total tb Espesor de la placa base tp Espesor de la placa piezoeléctrica C Capacitancia C Matriz de amortiguamiento 𝐂𝐄 Tensor de rigidez 𝐃 Desplazamiento de carga dieléctrica E Campo eléctrico Eb Módulo elástico del material base Ep Módulo elástico del material piezoeléctrico F Fuerza aplicada en el extremo libre K Matriz de rigidez L Longitud de la viga M Matriz de masa Q Carga 𝐒𝐄 Tensor de compliancia 𝑆11 E Componente del tensor de compliancia en el eje vertical U Energía total del sistema V Voltaje W Ancho de la viga α Relación entre el módulo elástico del material base y el módulo elástico del material piezoeléctrico αR Parámetro de amortiguamiento de Rayleigh β Relación entre el espesor del material base y el espesor del material piezoeléctrico βR Parámetro de amortiguamiento de Rayleigh ε Deformación mecánica ξ Factor de amortiguamiento σ Tensión mecánica ϵ33 T Componente del tensor de permitividad en el eje horizontal ϵT Permitividad 1 Introducción 2 1 INTRODUCCIÓN n este capítulo se expondrán las razones que motivan el estudio realizado en este trabajo y se llevará a cabo una revisión de la literatura publicada en el campo de los dispositivos piezoeléctricos. Posteriormente, se expondrán los objetivos perseguidos y la organización que sigue el contenido presentado. 1.1. Motivación El desarrollo tecnológico de los últimos años ha dado lugar a la proliferación de numerosos dispositivos inalámbricos y autónomos en ámbitos tan variados como la medicina, la industria o la ingeniería civil. La presencia de estos dispositivos plantea nuevos retos, como el uso de energías alternativas que permitan la independencia de baterías convencionales y evitar las desventajas que estas conllevan: duración limitada, necesidad de reemplazo o aporte de peso entre otras. Una posible fuente de energía para estos dispositivos podría basarse en la conversión de energía mecánica a energía eléctrica mediante el uso de materiales y sistemas piezoeléctricos. Los materiales piezoeléctricos permiten la generación de carga eléctrica en sus superficies al experimentar una deformación mecánica, es el llamado efecto piezoeléctrico directo. Del mismo modo, estos materiales se deforman al hacer pasar por ellos un voltaje eléctrico, el efecto piezoeléctrico inverso (ver Figura 1.) (a) (b) Figura 1: (a) Efecto piezoeléctrico directo, (b) efecto piezoeléctrico inverso (extraído de [1]) Este efecto hace a los recolectores 1 energéticos piezoeléctricos útiles ante cualquier tipo de excitación externa susceptible de crear deformaciones mecánicas en el material: vibraciones acústicas, incidencia del viento o dilataciones térmicas. Este sistema de recolección de energía tiene numerosas ventajas sobre otro tipo de fuentes de energía, como ser ligeros o no generar dispersión de calor [2]. Además, en los últimos años se ha investigado acerca de diversos aspectos relacionados con su rendimiento, obteniendo resultados que avalan su utilidad. Por ejemplo, se han estudiado técnicas para aumentar el rango de frecuencias de excitación del recolector. Algunas de estas son el empleo de múltiples vigas en voladizo (ver Figura 2), la adición o eliminación de una masa puntual en el extremo libre de la viga o el uso de imanes para provocar la repulsión magnética en la viga [3]. También se ha verificado el buen comportamiento de los recolectores piezoeléctricos bajo fatiga, demostrando una tasa limitada de degradación del rendimiento para ciclos de carga equivalentes hasta a 3 años de carga de tráfico [4]. Incluso se ha demostrado su viabilidad económica para la alimentación de dispositivos de demandas de energía bajas o medias, como en puertas automáticas [5]. Este conjunto de ventajas permite que los recolectores piezoeléctricos se empleen en multitud de ámbitos. Ejemplo de ello es su uso en aplicaciones en altas temperaturas, como en el ámbito aeronáutico, automóvil o en industrias energéticas [6]. En este contexto, los dispositivos piezoeléctricos son una de las alternativas 1 A lo largo de este trabajo se denominará a los generadores piezoeléctricos “recolectores piezoeléctricos”, por su traducción del término en inglés, harvester. E 3 3 preferidas de entre las técnicas resistentes a las altas temperaturas, por su bajo coste y su pequeño tamaño. También han resultado de especial utilidad en el contexto de la monitorización de la salud estructural en ingeniería civil. La creciente atención al mantenimiento de las infraestructuras civiles ha ocasionado la aparición de estándares, procedentes de administraciones nacionales u organizaciones, que recomiendan o imponen el uso de técnicas de monitorización de estructuras [7]. El uso de dispositivos como sensores piezoeléctricos resulta de utilidad en este contexto, dado que permiten medir parámetros como deformaciones, corrosión o fisuras con sensores baratos y de bajo mantenimiento [8]. Finalmente, poseen un gran potencial para la monitorización de aerogeneradores. La sensorización de los aerogeneradores requiere, tradicionalmente, un gran porcentaje de los costes asociados a la generación de energía eólica, dado que las localizaciones de estos aerogeneradores pueden ser remotas y el acceso a sus sensores es complicado. Por ello, esta industria se beneficia de los sensores piezoeléctricos por su larga vida útil y su escaso mantenimiento [9]. De la investigación llevada a cabo hasta ahora se puede extraer que los recolectores piezoeléctricos resultan de utilidad en contextos muy variados. Sin embargo, esto supone que las condiciones de excitación para cada uno de estos recolectores podrían ser exclusivas y específicas a cada aplicación. Del mismo modo, las propiedades y dimensiones de los recolectores podrían variar adaptándose a las condiciones y requisitos de cada entorno. Por lo tanto, resulta interesante evaluar el rendimiento de los recolectores en distintos escenarios y con distintas morfologías. Con este fin, se ha estudiado el efecto de parámetros como la longitud total de la viga piezoeléctrica o de la placa de material piezoeléctrico [10]. En el contexto de excitación aerodinámica, se ha evaluado la influencia de la superficie de exposición al viento y del peso de una masa puntual añadida al extremo de la viga [11]. También se ha comprobado la respuesta de un recolector piezoeléctrico bajo cargas puntuales estáticas y dinámicas [12]. Los estudios llevados a cabo en este trabajo pretenden continuar en esta línea de comprobaciones y averiguar la influencia de diversos parámetros geométricos, métodos de análisis y tipos de excitación en la eficiencia de un recolector piezoeléctrico. Asimismo, estos estudios son una continuación de trabajos previamente elaborados en el Dpto. Mecánica de Medios Continuos y Teoría de Estructuras, como el referenciado en [13]. Figura 2: Técnica de aumento de rango de excitación del recolector piezoeléctrico (extraído de [3]) Figura 3: Empleo de sensores piezoeléctricos en aerogeneradores (extraído de [9]) Introducción 4 1.2. Objetivos Este trabajo pretende modelizar una viga piezoeléctrica y diversas variaciones de ella para identificar parámetros relacionados con la cantidad de energía recolectada en distintos planteamientos. Así, se ensayarán diferentes morfologías de viga, tipos de solicitación y regímenes de funcionamientos con el fin de conocer la susceptibilidad del voltaje producido a la variación de alguno de estos parámetros. Para ello, en este trabajo se ha comprobado la relevancia en la recolección de energía piezoeléctrica de los siguientes aspectos: • Monitorización según modelos de elementos finitos (MEF) 2D o 3D. • Variaciones en las siguientes relaciones de aspecto: relación espesor-ancho, relación espesor-longitud y relación entre el espesor del material piezoeléctrico y el espesor del material sustrato. • Consideración de efectos de segundo orden. • Análisis estático o transitorio. La comparación de los resultados obtenidos permitirá reconocer las discordancias entre distintos casos de análisis y los contextos en los que son de utilidad los distintos modelos. Más importantemente, los resultados nos permitirán nos permitirán extraer conclusiones con respecto a los valores óptimos de las variables de diseño. 1.3. Organización del documento En este apartado se clarifica la organización del contenido presentado en los siguientes capítulos. • Capítulo 2: Modelización de vigas recolectoras. En este capítulo se presenta el esquema general de viga piezoeléctrica que se va a ensayar a lo largo de todo el trabajo. Se expone también el marco teórico que rige el comportamiento piezoeléctrico, que nos permitirá llevar a cabo comprobaciones analíticas y comprender los resultados arrojados por los modelos de elementos finitos. Se explicará el proceso de elaboración de los modelos de elementos finitos. • Capítulo 3: Análisis de los resultados. Se validarán los modelos de elementos finitos elaborados, comprobando que son coherentes y equivalentes. Posteriormente, se someterá al modelo a la variación de sus dimensiones para comprobar la relevancia que las relaciones de aspecto tienen sobre la eficacia del recolector. También se compararán los resultados obtenidos bajo el análisis de primer orden y el de segundo orden, extrayendo conclusiones acerca de su impacto en los resultados. Por último, se obtendrá la respuesta del recolector piezoeléctrico bajo un análisis transitorio. • Capítulo 4: Resumen y conclusiones. Se sintetizará el contenido desarrollado a lo largo de todo el documento, así como los resultados de mayor interés. Se expondrán las conclusiones obtenidas que pretenden satisfacer los objetivos inicialmente planteados. Finalmente, se indicarán algunos de los posibles trabajos futuros que podrían contribuir a la investigación del comportamiento de recolectores piezoeléctricos. 5 5 13 U=2(1+𝛼𝛽)𝑠11 𝐸𝑡𝑝𝐿3 𝛼2𝐵𝑊 𝐹2−3𝑑31𝛽(1+𝛽)𝑠11 𝐸 2𝑡𝑝3𝐿2 𝛼𝐵 (𝑉 𝑡𝑝)𝐹 +𝜖33 𝑇𝑊𝐿𝑡𝑝 2[1+(3𝛽2(1+𝛽)2𝑡𝑝4𝑠11 𝐸 2 (1+𝛼𝛽)𝐵 −1)𝑑31 2 𝑠11 𝐸𝜖33 𝑇](𝑉 𝑡𝑝)2 (23) La carga eléctrica se puede obtener derivando la energía total, en (23), con respecto a la tensión: 𝑄=𝑑𝑈 𝑑𝑉=−3𝑑31𝛽(1+𝛽)𝑠11 𝐸2𝑡𝑝2𝐿2 𝛼𝐵 𝐹+𝜖33 𝑇𝑊𝐿𝑡𝑝 2[1+(3𝛽2(1+𝛽)2𝑡𝑝4𝑠11 𝐸2 (1+𝛼𝛽)𝐵 −1) 𝑑31 2 𝑠11 𝐸𝜖33 𝑇]V (24) Si consideramos únicamente la carga producida por la actuación de la fuerza: 𝑄=−3𝑑31𝛽(1+𝛽)𝑠11 𝐸 2𝑡𝑝2𝐿2 𝛼𝐵 𝐹 (25) La capacitancia obtenida es: 𝐶=𝜖33 𝑇𝑊𝐿 𝑡𝑝[1+(3𝛽2(1+𝛽)2𝑡𝑝4𝑠11 𝐸 2 (1+𝛼𝛽)𝐵 −1)𝑑31 2 𝑠11 𝐸𝜖33 𝑇] (26) Por último, conociendo que el voltaje se define como la carga entre la capacitancia: 𝑉=−3 𝑑31𝛽(1+𝛽)𝑠11 𝐸𝑡𝑝2𝐿 𝜖33 𝑇𝑊𝑡𝑝𝐵[1+(3𝛽2(1+𝛽)2𝑡𝑝4𝑠11 𝐸 (1+𝛼𝛽)𝐵 −1)𝑑31 2 𝑠11 𝐸𝜖33 𝑇]𝐹 (27) 2.2.3. Descripción del modelo teórico Las expresiones obtenidas durante el desarrollo teórico llevado a cabo en los apartados anteriores permiten confeccionar un modelo de cálculo analítico. Este se empleará para la validación de los modelos, comprobando la coherencia entre sus resultados y aquellos obtenidos mediante el método de elementos finitos. Para simulaciones de mayor complejidad se usarán únicamente los de elementos finitos. Para la elaboración de este modelo, se han integrado algunas de estas ecuaciones en un programa del software MATLAB, cuyo código se puede consultar en el Anexo 1. Así, se han usado las ecuaciones (16), (25) y (26), para obtener la flecha máxima, la carga y la capacitancia respectivamente. El voltaje se ha obtenido relacionando los resultados de carga y capacitancia. Como se puede observar en las expresiones (16), (25) y (26), los parámetros implicados son los parámetros geométricos presentados en el apartado 2.1 y las componentes de los tensores materiales que resultan de interés para nuestro caso de carga. Para una mayor aclaración, se explican estos últimos: • 𝑠11 𝐸 es la componente del tensor de compliancia en el eje vertical • 𝜖33 𝑇 es la componente del tensor de permitividad en el eje horizontal • 𝑑31 es el coeficiente piezoeléctrico para el modo de carga 31 • 𝐸𝑝, incluido en 𝛽, se obtiene según 1 𝑠11 𝐸 A diferencia de los modelos de elementos finitos, este modelo no permite visualizar las distribuciones de estas propiedades a lo largo de la viga ni los puntos en los que tienen lugar los valores máximos cuantificados. No Modelización de vigas recolectoras de energía 14 obstante, es evidente por el esquema de carga y condiciones de contorno presentado en el apartado 2.1, que la flecha máxima tendrá lugar en el extremo libre de la viga. Asimismo, el voltaje máximo se producirá en el extremo empotrado, donde las tensiones serán mayores como consecuencia de los grados de libertad restringidos. 2.3. Modelo de elementos finitos 3D Para realizar algunas de las simulaciones se ha construido un modelo de elementos finitos 3D en el software ANSYS. En este apartado se indicarán las fases necesarias para la creación del modelo, así como los parámetros a introducir. El código de este modelo se puede consultar en el Anexo 2. En primer lugar, será necesario seleccionar los elementos a utilizar para la discretización de la geometría. Para el caso 3D se han elegido los elementos tipo SOLID186 y SOLID226. El elemento SOLID186 es un elemento cuadrático 3D de 20 nodos. Cada nodo posee 3 grados de libertad: translación en x, translación en y y translación en z. Este elemento ha sido seleccionado para el material sustrato. Para el material piezoeléctrico se han elegido elementos SOLID226 porque permiten modelizar el comportamiento en caso de campos acoplados. Al igual que el anterior, este elemento cuadrático 3D posee 20 nodos. No obstante, en este caso cada nodo tendrá 4 grados de libertad: las tres translaciones según los ejes y el voltaje producido. Figura 8: Geometría de los elementos SOLID186 y SOLID226 (extraída de [17]) A continuación, se introducen las propiedades de ambos materiales. Mientras que las del acero inoxidable se introducen mediante valores escalares, habrá que introducir las del material piezoeléctrico de manera matricial. El programa se basa en una formulación de las ecuaciones constitutivas piezoeléctricas equivalente a las introducidas en (1), pero que emplea variables alternativas: {𝝈𝑫}=[𝑪𝑬−𝒆𝑻 𝒆 𝝐𝑺]{𝜺𝑬} (28) Las variables empleadas en esta formulación se definen en la Tabla 3. Las relaciones entre las variables empleadas en las dos formulaciones de las ecuaciones constitutivas son las siguientes: 𝑺𝑬=[𝑪𝑬]−𝟏 (29) 𝒆=𝒅𝒄𝑬 (30) 15 𝝐𝑺=𝝐𝑻−𝒅𝒆𝑻 (31) La matriz de rigidez posee una dimensión de 6x6, atendiendo a sus 6 grados de libertad: 3 en las direcciones principales y 3 en las tangenciales: 𝑪𝑬= [ 𝑐11 𝑐12 𝑐13 0 0 0 𝑐12 𝑐11 𝑐13 0 0 0 𝑐13 𝑐13 𝑐33 0 0 0 0 0 0 𝑐44 0 0 0 0 0 0 𝑐44 0 0 0 0 0 0 𝑐66 ] El tensor de coeficientes piezoeléctricos tiene una dimensión de 3x6: 𝒆=[0 0 0 0 𝑒15 0 0 0 0 𝑒15 0 0 𝑒31 𝑒31 𝑒33 0 0 0] Por último, la matriz de permitividad es de 3x3: 𝝐𝑺=[𝜖11 0 0 0 𝜖11 0 0 0 𝜖33] El superíndice “S” acompañando al tensor de permitividad indica que la propiedad es medida a deformación constante. Además de estos valores de permitividad absoluta, habrá que especificar el valor de permitividad en el vacío, para poder computar la permitividad relativa. Posteriormente, se define la geometría, creando volúmenes a partir de puntos. Resulta evidente que el modelo 2D distará del 3D en este apartado, dado que en el primero tendremos que modelizar áreas. La geometría resultante es mallada y se puede proceder a la aplicación de condiciones de contorno y fuerzas actuantes en base a los nodos de la malla. La geometría mallada se puede consultar en la Figura 9. Para simular el empotramiento se restringen los 3 grados de libertad de los nodos en una de las caras extremas del material sustrato. La solicitación se modeliza como una serie de cargas puntuales actuando sobre los nodos del extremo no empotrado. El valor de fuerza impuesto se reparte, por tanto, a lo largo del ancho de la viga. La disposición de la carga se representa en la Figura 10. Finalmente, se puede resolver el modelo. En los casos en los que se haya llevado a cabo un análisis de segundo orden, se activa en este punto el comando de análisis no lineal, “NLGEOM”. Este comando habilita la consideración de grandes desplazamientos y tiene en cuenta el efecto de rigidización por tensión en la matriz de rigidez. Propiedad Símbolo Unidad Tensión mecánica 𝝈 Pa Desplazamiento de carga dieléctrica 𝑫 C/𝑚2 Tensor de rigidez 𝑪𝑬 𝑃𝑎 Tensor de coeficientes piezoeléctricos 𝒆 C/𝑚2 Deformación mecánica 𝜺 N/𝑚2 Campo eléctrico E N/C Tabla 3: Variables implicadas en las ecuaciones constitutivas Modelización de vigas recolectoras de energía 16 Figura 9: Viga piezoeléctrica 3D mallada Figura 10: Aplicación de la carga en el extremo libre 2.4. Modelo de elementos finitos 2D A lo largo del desarrollo de las simulaciones se presentan en diversos casos morfologías para las que el ancho es significativamente menor que el espesor total. Consecuentemente, resulta de interés plantear un modelo 2D para estudiar el problema como tensión plana, para lograr un menor esfuerzo computacional. El código de este modelo se puede consultar en el Anexo 2. 17 Esto requerirá, en primer lugar, cambiar los tipos de elementos seleccionados. En este caso se ha recurrido a los elementos PLANE183 y PLANE223. El elemento PLANE183 es un elemento 2D cuadrático de 8 nodos con dos grados de libertad: traslación en x y traslación en y. Este elemento se ha usado para modelizar la placa sustrato. Al igual que en el caso 3D, para la modelización del material piezoeléctrico se requiere un elemento que permita el acoplamiento de campos. Por ello, se ha empleado el elemento PLANE223, elemento 2D cuadrático con 8 nodos con 3 grados de libertad: dos traslaciones según los ejes x e y y el voltaje. Figura 11: Geometría de los elementos PLANE183 y PLANE223 (extraída de [17]) Los valores de las propiedades materiales no cambian para el modelo 2D, dado que se emplearán lo mismos materiales. Sin embargo, sí varía la estructura matricial de estos. En primer lugar, porque al tener menos grados de libertad, todos los tensores implicados tendrán dimensiones más reducidas. En segundo lugar, porque como se puede observar en la Figura 11, los elementos están definidos por los ejes x e y, de modo que no es posible crear áreas 2D en el plano XZ. Dado que, atendiendo a los ejes definidos en la Figura 6, nos interesa proyectar el volumen sobre nuestro plano XZ, será necesario rotar los tensores materiales de modo que el eje y pase a representar el eje vertical. El eje x continuará siendo el eje longitudinal. Solo a efectos de explicar la conversión de las matrices materiales desde el caso 3D al caso 2D, se denotarán las variables 2D con el indicador “ * ”. • Tensor de rigidez: El tensor de rigidez pasará a ser de dimensión 4x4, debido a las tensiones en las direcciones normales x, y, z y la componente tangencial xy. El tensor tendrá la forma: 𝑪∗=[𝑐11∗ 𝑐12∗ 𝑐13∗ 0 𝑐21∗ 𝑐22∗ 𝑐23∗ 0 𝑐31∗ 𝑐32∗ 𝑐33∗ 0 0 0 0 𝑐44∗] Y sus componentes se corresponderán con las componentes del tensor 3D según: Modelización de vigas recolectoras de energía 18 Tensor de rigidez Valor en 2D Valor en 3D 𝑐11∗ 𝑐11 𝑐12∗ 𝑐13 𝑐13∗ 𝑐12 𝑐21∗ 𝑐13 𝑐22∗ 𝑐33 𝑐23∗ 𝑐13 𝑐31∗ 𝑐12 𝑐32∗ 𝑐13 𝑐33∗ 𝑐11 𝑐44∗ 𝑐66 Tabla 4: Equivalencia entre las componentes del tensor de rigidez 3D y 2D • Tensor de coeficientes piezoeléctricos: Pasa a ser de dimensiones 2x4. Tendrá la siguiente estructura: 𝒆∗=[0 0 0 0 𝑒21∗ 𝑒22∗ 0 0] La equivalencia entre los coeficientes 2D y 3D será según: Tensor piezoeléctrico Valor en 2D Valor en 3D 𝑒21∗ 𝑒13 𝑒22∗ 𝑒33 Tabla 5: Equivalencia entre las componentes del tensor piezoeléctrico 3D y 2D • Matriz de permitividad: La matriz de permitividad queda de dimensiones 2x2, por las dos posibles direcciones de aplicación de campos eléctricos y las dos posibles direcciones de polarización: 𝝐∗=[𝜖11∗ 0 0 𝜖22∗] La correspondencia entre los coeficientes en 2D y en 3D se muestra en la siguiente tabla: 19 Tensor de permitividad Valor en 2D Valor en 3D 𝜖11∗ 𝜖11 𝜖22∗ 𝜖33 Tabla 6: Equivalencia entre las componentes del tensor de permitividad 3D y 2D Una vez definidas las propiedades materiales, se puede confeccionar la geometría de la viga. Como se ha indicado previamente, en este caso se representan áreas en el plano XY, considerando el ancho despreciable frente al resto de dimensiones. El resto del proceso será paralelo al proceso para el modelo 3D. Se ejecutará el mallado con los elementos 2D y se impondrán las condiciones de contorno y las solicitaciones. En este caso, el empotramiento de una de las caras extremas solo requerirá restringir 2 grados de libertad. El valor de la solicitación será, al igual que en el caso anterior, el correspondiente a la fuerza nominal distribuida a lo largo del ancho de la viga. Figura 12: Modelo de la viga piezoeléctrica en 2D 2.5. Modelización de la respuesta dinámica Los modelos mostrados hasta este punto del documento son modelos estáticos, en los que las solicitaciones y las respuestas se consideran independientes del tiempo. No obstante, posee interés considerar la respuesta del recolector en función de una carga transitoria o dinámica, de modo que se puedan simular cargas más realistas, como aquellas que actúan de manera intermitente. Para ello, se ha adaptado el modelo de elementos finitos 2D a un análisis transitorio en el que la solicitación es una carga tipo escalón. Se ha seleccionado una carga escalón para lograr excitar todas las frecuencias naturales. Las tareas del preprocesado de este modelo son idénticas a aquellas del preprocesado estático, ya que se trata de la misma morfología y los mismos materiales. Por ello, se describirá únicamente el proceso de elaboración de la sección de solución. En primer lugar, se definirá el análisis en tipo transitorio y se activará el cálculo en segundo orden en los casos en los que sea necesario. Como se ha indicado previamente, la carga se ha aplicado en forma de escalón, de modo que la secuencia de carga es: Modelización de vigas recolectoras de energía 20 Figura 13: Secuencia de carga escalón aplicada Como se observa, la simulación se divide en tres pasos. El primero de ellos acaba en el instante 0 y almacena las condiciones iniciales. Estas condiciones iniciales corresponden a la ausencia de carga aplicada y consecuentemente, a flecha y voltaje nulos. En el segundo paso, se aplica una carga de 700 N en el extremo libre. Este paso comienza donde finaliza el anterior y acaba en 0.001 s. El tercer instante acaba en 0.1 s y en él se mantiene la carga. Sirve para finalizar la simulación. En estos pasos de carga intervienen varios parámetros que cabe destacar: • Interviene un diferencial de tiempo que fija el número de subpasos que se toman en cada paso de carga. Este parámetro resulta determinante para la convergencia de la solución, dado que para cada uno de estos subpasos se integrará un valor de la respuesta. Para aquellos casos en los que el comportamiento sea no lineal, la curva poseerá puntos singulares. Si el subpaso es pequeño, la curva integrada podrá adaptarse mejor al entorno de estos puntos y proveer una solución adecuada. De lo contrario, la respuesta será no convergente. Se debe tener en cuenta que mientras menor sea el subpaso, más precisa será la respuesta, pero mayor será el esfuerzo computacional. El objetivo será encontrar una respuesta suficientemente precisa sin malgastar esfuerzo computacional. • También intervienen los parámetros α y β de amortiguamiento de Rayleigh. Este tipo de amortiguamiento supone que este es proporcional a la masa y a la rigidez del sistema. 𝐂=αR𝐌+βR𝐊 (32) En este caso se ha determinado αR=0 y se ha estudiado la respuesta para dos valores de βR, de modo que el amortiguamiento será independiente de la masa y proporcional en mayor o menor medida a la rigidez. En los resultados se comprobará la influencia del parámetro βR en la respuesta. Finalmente, se solucionan los tres pasos y se representarán las respuestas con respecto al tiempo. 21 3 ANÁLISIS DE LOS RESULTADOS n el presente capítulo se describirán los estudios ensayados con los modelos introducidos en el Capítulo 2. Para cada uno de estos casos de estudio se darán a conocer los valores de los parámetros de entrada, los resultados obtenidos y se extraerán conclusiones al respecto del comportamiento y la eficacia del recolector. 3.1. Propiedades de los materiales A diferencia de las dimensiones del recolector, los materiales empleados son los mismos en todas las simulaciones. Consecuentemente, no cambiarán los valores de los parámetros materiales, aunque sí su estructura para cada modelo, según los expuesto en el Capítulo 2. Por ello, se exponen aquí los valores de estos parámetros, extraídos del artículo referenciado en [18]. • Propiedades materiales del modelo MEF 3D Propiedades PZT-4 Símbolo Valor Unidad 𝑐11 139 GPa 𝑐12 77,8 GPa 𝑐13 74,3 GPa 𝑐33 115 GPa 𝑐44 25,6 GPa 𝑐66 30,6 GPa 𝑒31 -5,2 C/𝑚2 𝑒33 15,1 C/𝑚2 𝑒15 12,7 C/𝑚2 𝜖11/𝜖0 763 - 𝜖33/𝜖0 658 - 𝜌𝑝 7495 𝑘𝑔/𝑚3 ν𝑝 0,3 - Tabla 7: Propiedades materiales del material PZT-4 en el modelo MEF 3D E 3 Análisis de los resultados 22 Propiedades de acero inoxidable Símbolo Valor Unidad ρb 7750 kg/m3 Eb 195 GPa νb 0,3 − Tabla 8: Propiedades materiales del acero inoxidable en el modelo MEF 3D Así, las matrices quedarían: 𝑪𝑬= [ 𝑐11 𝑐12 𝑐13 0 0 0 𝑐12 𝑐11 𝑐13 0 0 0 𝑐13 𝑐13 𝑐33 0 0 0 0 0 0 𝑐44 0 0 0 0 0 0 𝑐44 0 0 0 0 0 0 𝑐66 ] = [ 139 77,8 74,3 0 0 0 77,8 139 74,3 0 0 0 74,3 74,3 115 0 0 0 0 0 0 25,6 0 0 0 0 0 0 25,6 0 0 0 0 0 0 30,6 ] 𝒆=[0 0 0 0 𝑒15 0 0 0 0 𝑒15 0 0 𝑒31 𝑒31 𝑒33 0 0 0]=[0 0 0 0 12,7 0 0 0 0 12,7 0 0 −5,2 −5,2 15,1 0 0 0] 𝝐𝑺=[𝜖11 0 0 0 𝜖11 0 0 0 𝜖33]=[763 0 0 0763 0 0 0 658] • Propiedades materiales del modelo MEF 2D Propiedades de acero inoxidable Símbolo Valor Unidad ρb 7750 kg/m3 Eb 195 GPa νb 0,3 − Tabla 9: Propiedades materiales del acero inoxidable en el modelo MEF 2D 29 se exponen en las tablas 14 y 15. Asimismo, para una mejor comprensión de la variación de la geometría de la viga al variar su relación espesor-ancho, se representan aquí los modelos resultantes de imponer 𝑡 𝑊=0,1 y 𝑡 𝑊=4 en las figuras 19 (a) y (b) respectivamente. Cabe destacar que las simulaciones de este apartado se han realizado con el modelo 3D, pues para los tres escenarios con menor relación espesor-ancho, el ancho no es despreciable frente al resto de dimensiones y no es válida la hipótesis de viga delgada. Aplicando una sucesión de carga de 0 a 350 N con un paso de 50 N, se obtiene la relación de flecha frente a carga presentada en la Figura 20 y la relación de voltaje frente a carga representada en la Figura 21. Se ha normalizado la carga con respecto al valor máximo aplicado, 𝐹0=350 𝑁 y el voltaje, con respecto al voltaje máximo obtenido en el caso de 𝑡 𝑊=0,1, 𝑉0=715,57 V. En las gráficas mostradas se puede observar que mientras mayor sea la relación 𝑡 𝑊, mayores son la flecha y el voltaje generado a igualdad de carga impuesta. La relación que guardan la fuerza y la flecha no es más que la rigidez de la viga, de modo que la pendiente de la curva es la rigidez. La rigidez a flexión de una viga se define según la siguiente expresión: 𝐾𝑓𝑙𝑒𝑥,𝑧=12⋅𝐸⋅𝐼𝑧 𝐿3, donde 𝐸 es el módulo de Young, 𝐼𝑧 es el módulo de inercia de la sección resistente a la flexión y L es la longitud de la viga. Este módulo de inercia corresponde al de un rectángulo, que atendiendo a nuestras medidas será: 𝐼𝑧=𝑊⋅𝑡3 12. A mayor ratio entre espesor y ancho, menor el ancho y consecuentemente, menor el módulo de inercia de la sección. Así, la viga posee menor rigidez a flexión y flecta más para un mismo valor de fuerza. Por lo tanto, tiene lugar una mayor deformación en la viga y se genera mayor voltaje. Se puede concluir que la viga es más eficiente mientras mayor sea la relación espesor-ancho. Para recalcar el impacto que el ancho de la viga tiene en la flecha y en el voltaje, se muestran en las figuras 22 (a) y (b) las deformadas de la viga para el caso 𝑡 𝑊= 0,1 y 𝑡 𝑊= 4 a igualdad de carga y con un factor de escala 5. Asimismo, se muestran en las figuras 23 (a) y (b) la distribución del voltaje generado en las superficies del material piezoeléctrico para los casos 𝑡 𝑊= 0,1 y 𝑡 𝑊= 4 respectivamente. Valores de relación 𝑡 𝑊 0,1 0,5 1 2 4 t (m) 0,008 0,008 0,008 0,008 0,008 𝑊 (m) 0,08 0,016 0,008 0,004 0,002 Tabla 14: Valores del espesor y el ancho según los ratios considerados Dimensiones Símbolo Valor Unidad 𝐿 0,08 m 𝑡𝑝 0,002 m 𝑡𝑏 0,006 m Tabla 15: Valores de los parámetros geométricos para la comprobación de la influencia de 𝑡 𝑊 3 Análisis de los resultados 30 Figura 19: (a) Geometría para 𝑡 𝑊=0,1,𝑡𝐿=0,1 y β=3 (b) Geometría para 𝑡 𝑊=4,,𝑡𝐿=0,1 y β=3 Figura 20: Influencia de 𝑡 𝑊 en la relación fuerza-flecha (𝑡𝐿=0,1 y β=3 ) 31 Figura 21: Influencia de 𝑡 𝑊 en la relación fuerza-voltaje (𝑡𝐿=0,1 y β=3 ) (a) (b) Figura 22: (a) Deformada para 𝑡 𝑊= 0,1, 𝑡𝐿=0,1 y β=3 (b) Deformada para 𝑡 𝑊= 4, 𝑡𝐿=0,1 y β=3 3 Análisis de los resultados 32 (a) (b) Figura 23: (a) Voltaje para 𝑡 𝑊= 0,1, 𝑡𝐿=0,1 y β=3 (b) Voltaje para 𝑡 𝑊= 4, 𝑡𝐿=0,1 y β=3 3.2.2.1. Influencia de efectos de segundo orden según 𝒕 𝑾 Llegados a este punto del documento, todos los resultados se han obtenido según un análisis de primer orden. Cabe revisar el fundamento del análisis de primer orden para identificar las implicaciones que esto haya podido tener en nuestros resultados. La teoría de primer orden supone que es correcto calcular las deformaciones con respecto a la configuración indeformada, por su similitud con la situación deformada. Esta hipótesis es válida en el caso de que los desplazamientos sean pequeños. En el caso de que esto no se cumpla, la solución en primer orden arrojará unos resultados incorrectos. Existe una segunda fuente de no linealidad en el cálculo de los desplazamientos de la viga, que también se desprecia como consecuencia de la hipótesis de pequeños desplazamientos. Se trata de la rigidización por tensión. Podemos asumir un sólido que tras la aplicación de una carga se distorsione significativamente, de modo que en ciertos puntos de su dominio hayan tenido lugar desplazamientos importantes. Inmediatamente deja de ser aplicable la hipótesis de pequeños desplazamientos. Estos desplazamientos locales provocarán tensiones locales y consecuentemente, esfuerzos locales. Estos esfuerzos forman parte de las fuerzas internas que la matriz de rigidez ha de relacionar con los desplazamientos. De este modo, aparte de la componente de rigidez relacionada con el material y la geometría de la viga, aparece una componente de rigidez relacionada con los esfuerzos internos que surgen como consecuencia de una distorsión geométrica importante. Esta componente es la matriz elemental de rigidez geométrica, que proporciona un aumento de la capacidad portante de la estructura como consecuencia de su estado tensional. Este efecto resulta de mayor relevancia en estructuras con poca rigidez a flexión y mucha rigidez axial, como vigas delgadas o placas [17]. Los métodos de elementos finitos tienen en cuenta la rigidización por tensión mediante la obtención de la matriz elemental de rigidez geométrica en cada punto. Para ello se llevan a cabo en cada punto sucesivas iteraciones dentro de un mismo paso de carga, con el fin de extraer el valor de las fuerzas internas en función de los desplazamientos locales y así actualizar el valor de la matriz de rigidez geométrica [18]. Como se ha razonado previamente, el requisito indispensable para poder adoptar la hipótesis de un comportamiento lineal es que los desplazamientos sean pequeños con respecto a la geometría del sólido. En este apartado se pretende comprobar si esta hipótesis arroja resultados verosímiles para nuestro recolector piezoeléctrico y en qué casos lo hace, de modo que se ejecutarán análisis de segundo orden y se compararán con los resultados según el análisis de primer orden. El análisis de segundo orden tendrá en cuenta los dos aspectos de no linealidad que se han introducido: los términos de segundo orden en el tensor de deformaciones y la rigidización del comportamiento por efecto de la tensión. El método de cálculo no lineal que sigue el software es el de Newton-Raphson. De este modo, se han repetido las simulaciones llevadas a cabo en este apartado, pero considerando los efectos de segundo orden. En la Figura 24 se representan los resultados obtenidos en comparación con sus homólogos 33 en primer orden. Como se puede observar, la diferencia entre los resultados según el cálculo en primer orden y el cálculo en segundo orden es muy pequeña. Esta diferencia aumenta mientras mayor sea la relación 𝑡 𝑊 y es más evidente en el caso del voltaje que en la flecha, donde el voltaje en segundo orden es mayor que en primer orden. No obstante, esta discordancia es de un 3% en caso más desfavorable, el voltaje para 𝑡 𝑊=4. Consecuentemente, la diferencia se puede considerar despreciable en todos los casos. Estos resultados ponen en evidencia que, si bien anteriormente se ha demostrado que la relación espesor-ancho influye en gran medida en la flecha obtenida, en ningún caso de los estudiados esta flecha es lo suficientemente grande para invalidar la hipótesis de pequeños desplazamientos. Es decir, los términos no lineales implicados en los tensores de deformaciones serán despreciables. 3 Análisis de los resultados 34 (a) (b) (c) (d) (e) (f) Figura 24: (a) Influencia de 2º orden en fuerza-flecha para 𝒕𝒘=0,1 (b) Influencia de 2º orden en fuerza-voltaje para 𝒕𝒘=0,1 (c) Influencia de 2º orden en fuerza-flecha para 𝒕𝒘=1 (d) Influencia de 2º orden en fuerza-voltaje para 𝒕𝒘=1 (e) Influencia de 2º orden en fuerza-flecha para 𝒕𝒘=4 (f) Influencia de 2º orden en fuerza-voltaje para 𝒕𝒘=4 35 3.2.3. Influencia de la relación de aspecto 𝒕𝑳 y los efectos de segundo orden Una vez expuesta la relevancia de la relación de aspecto del espesor de viga frente al ancho, cabe plantearse la influencia de otras dimensiones sobre el voltaje generado, como la longitud de la viga. En primera instancia, este análisis se efectuará en primer orden. Al igual que en el apartado previo, se estudiarán los resultados de flecha y voltaje para distintas relaciones del espesor total frente a la longitud, estas son: 0,1; 0,05; 0,03; 0,02 y 0,01. A igualdad de canto, estas relaciones resultan en las medidas expuestas en la tabla 16. Estas simulaciones se han llevado a cabo para un ratio de 𝑡 𝑊= 4 y β=3, de modo que el resto de dimensiones quedan como se expone en l Tabla 17. De este modo, ensayando el modelo para unas cargas de 0 a 350 N, con un paso de 50 N, obtenemos la evolución carga-flecha presentada en la Figura 25 y una evolución carga-voltaje presentada en la Figura 26. En estas gráficas se ha normalizado el valor de carga con respecto a la carga máxima aplicada (𝐹0=350 𝑁) y el voltaje, con respecto al voltaje máximo obtenido en el caso de 𝑡𝐿=0,1 (𝑉0=15390 V). Si evaluamos las gráficas, se observa que a menor ratio de 𝑡𝐿, más flecha y más voltaje se obtiene. Sin embargo, la diferencia entre los resultados para distintas relaciones es considerablemente más drástica en el caso de la flecha que en el del voltaje. Esta gran influencia se puede explicar de nuevo considerando la aparición de L en la rigidez de la viga. Así, consultando la expresión en (32) se observa que la longitud aparece en el denominador elevada a la tercera potencia. Del mismo modo, se puede consultar la expresión (16) de la que se obtiene la flecha máxima. La longitud aparece en el numerador, de nuevo elevada a la tercera potencia. Es decir, un incremento pequeño de la longitud aumenta significativamente la flecha y con ello, la deformación y voltaje generados. La relación directa entre la longitud y el voltaje se puede consultar en (27). En dicha expresión se observa que la longitud aparece en el numerador elevada a 1, de modo que el voltaje aumenta linealmente con la longitud. Esto explica que la diferencia entre los voltajes obtenidos para distintos valores de 𝑡𝐿 sea menos drástica que en el caso de la flecha. En conclusión, se demuestra que a igualdad de canto y de ancho, un recolector piezoeléctrico es más eficiente mientras mayor sea su longitud. Valores de relación 𝑡𝐿 0,1 0,05 0,03 0,02 0,01 t (m) 0,008 0,008 0,008 0,008 0,008 𝐿 (m) 0,08 0,16 0,267 0,4 0,8 Tabla 16: Valores del espesor y la longitud para distintas relaciones de 𝑡𝐿 Dimensiones Símbolo Valor Unidad 𝑡𝑝 0,002 m 𝑡𝑏 0,006 m 𝑡 0,008 m 𝑊 0,002 m Tabla 17: Valores de las variables geométricas para la comprobación de la influencia de 𝑡𝐿 3 Análisis de los resultados 36 Figura 25: Influencia de 𝑡𝐿 en la relación fuerza-flecha ( 𝑡 𝑊=4 y β=3) Figura 26: Influencia de 𝑡𝐿 en la relación fuerza-voltaje (𝑡 𝑊=4 y β=3) 37 3.2.3.1. Influencia de efectos de segundo orden según 𝒕𝑳 Al igual que para la relación 𝑡 𝑊, cabe cuestionar la validez de estos resultados en los que se ha obviado el comportamiento no lineal del recolector. Para ello, se ha repetido el cálculo de la viga para los valores de 𝑡𝐿 impuestos en las simulaciones previas y para las mismas cargas, pero bajo el análisis de segundo orden. Los resultados obtenidos se han comparado con los resultados homólogos en primer orden y se exponen a continuación. En primer lugar, se muestra la comparativa de la flecha con respecto a la carga aplicada y en segundo lugar, el voltaje generado con respecto a la carga aplicada. Los valores de las gráficas han sido normalizados con respecto a la fuerza máxima aplicada (F0=350 N) y al voltaje máximo que se obtiene para 𝑡𝐿=0,1 según el análisis de primer orden (V0 = 15390 V). La primera conclusión que se puede extraer es que a menor relación 𝑡𝐿, más distan los resultados obtenidos en el análisis no lineal y el análisis lineal. Como ya se ha razonado previamente, la longitud tiene un gran peso en la rigidez y la flecha de la viga, dado que interviene en sus expresiones elevada a la tercera potencia. Esto provoca que para longitudes grandes, se obtengan flechas muy grandes, de modo que la hipótesis de pequeños desplazamientos no se cumple. Consecuentemente, habría que tener en cuenta los términos de deformaciones de segundo orden y la rigidización tensional, que puede ser de gran importancia en vigas delgadas. Como el análisis lineal desprecia estos dos aspectos, los resultados del análisis en primer orden distan en exceso de los resultados en segundo orden, que representan con más verosimilitud el problema. La discordancia entre los resultados de la flecha en análisis de primer orden y segundo orden tiene también se pone de manifiesto en el voltaje obtenido, aunque se puede observar que la variación es ligeramente más lenta. El voltaje generado está directamente relacionado con las deformaciones mediante el coeficiente piezoeléctrico e31. Si aplicamos la teoría de pequeños desplazamientos se desprecian los términos de segundo orden de la deformación, conduciendo a un valor erróneo de la deformación y, consecuentemente, del voltaje. Por el contrario, podemos observar como para las relaciones 𝑡𝐿 más pequeñas las soluciones lineales y no lineales son muy similares. Para elementos menos esbeltos la flecha es menor, la hipótesis de pequeños desplazamientos es aceptable y se pueden despreciar los aspectos no lineales sin que ello suponga un error importante en los resultados obtenidos de flecha y voltaje. 3 Análisis de los resultados 38 (a) (b) (c) (d) (e) Figura 27: (a) Influencia de efectos de 2º orden en la flecha para 𝑡𝐿=0,1 (b) Influencia de efectos de 2º orden en la flecha para 𝑡𝐿=0,05 (c) Influencia de efectos de 2º orden en la flecha para 𝑡𝐿=0,03 (d) Influencia de efectos de 2º orden en la flecha para 𝑡𝐿=0,02 (e) Influencia de efectos de 2ºorden en la flecha para 𝑡𝐿=0,01 45 Figura 36: Influencia de los efectos de 2º orden en el voltaje según β (𝑡𝐿=0,1 y 𝑡 𝑊=4) 3.2.5. Influencia de la combinación de relaciones de aspecto A lo largo del Capítulo 3 se ha evaluado el impacto que ciertas relaciones de aspecto tienen en el rendimiento del recolector piezoeléctrico, estudiando la relación entre los resultados obtenidos y las propiedades mecánicas y piezoeléctricas del recolector. En este último apartado del análisis estático se pretende comprobar la influencia global de estas relaciones de aspecto mediante la combinación de todas ellas. 3.2.5.1. Influencia de β y 𝒕 𝑾 en análisis de segundo orden Para comprobar el efecto de β sobre el voltaje generado según la relación de espesor-ancho, se ha analizado el voltaje producido para las mismas relaciones 𝑡 𝑊 que en el apartado 3.2.2, una carga de 350 N y diversos valores de β. Los resultados se exponen en la Figura 37. El voltaje se ha normalizado con respecto al voltaje máximo producido en la curva de la relación, V0=3211,5 V. De este análisis se extrae que, para las dimensiones estudiadas, el β óptimo es independiente de la relación espesor-ancho, puesto que el máximo de las curvas permanece en 0,3 para todos los casos. No obstante, en el apartado 3.2.5.3 se comprobará que esta cualidad no es extensiva a todos los casos. Además, se comprueba que la curva es más tendida mientras menor es el ratio. Es decir, para mayores ratios, el entorno del máximo en el que los valores de voltaje producido son cercanos al óptimo es menor. A menor ancho, más susceptible es la generación de energía a un cambio en la relación entre espesores. Asimismo, se observa que lejos del β óptimo, la diferencia en voltaje producido para distintas relaciones 𝑡 𝑊 es menos significativa. 3 Análisis de los resultados 46 Figura 37: Influencia de β en el voltaje según 𝑡 𝑊 (𝑡𝐿=0,1) 3.2.5.2. Influencia de β y 𝒕𝑳 en análisis de segundo orden En este apartado se pretende comprobar el impacto que la razón de dimensiones 𝑡𝐿 y el análisis no lineal puedan tener sobre el óptimo de la relación entre los espesores de los materiales. Con este fin se ha repetido el análisis de voltaje obtenido en función de 𝑡𝐿 variando los valores de 𝛽 y para una carga de 350 N. Los resultados se han representado en la Figura 38. El voltaje se ha normalizado con respecto al voltaje máximo obtenido para 𝑡𝐿=0,1, V0=85473 V. Se observa que, para una misma longitud, el cambio de β tiene un efecto limitado sobre el voltaje óptimo: una desviación pequeña del β óptimo provoca una disminución pequeña del voltaje obtenido. Esto se cumple para todas las relaciones 𝑡𝐿 contempladas, al contrario de lo observado en la correlación de β con la relación 𝑡 𝑊. Es decir, para todos los valores de 𝑡𝐿, el entorno en el que los valores de voltaje son cercanos al óptimo posee aproximadamente la misma amplitud, desde β=0,2 a β=0,4. Este resultando es distinto del obtenido para la respuesta a la variación de 𝑡 𝑊, en la que el entorno del óptimo cambiaba significativamente para distintas relaciones de 𝑡 𝑊, dando lugar a curvas más o menos tendidas. Se puede concluir que la susceptibilidad del voltaje ante al cambio de β es independiente de 𝑡𝐿. Se observa que el β óptimo se mantiene en 0,3 para la relación de 𝑡𝐿 más baja, cambiando hasta 0,4 para longitudes hasta 10 veces mayores. Así, se concluye que, al contrario de lo que sucedía para la relación 𝑡 𝑊 , el β óptimo sí varía según el ratio 𝑡𝐿. También se puede observar que las curvas con valores de 𝑡𝐿 pequeños se agrupan, de modo que son casi coincidentes en el entorno del β óptimo. La distancia entre ellas se hace mayor en valores más lejanos del óptimo. Por lo tanto, se concluye que para una razón de los espesores de los materiales óptima o cercana a la óptima, se obtienen rendimientos muy similares para 𝑡𝐿 desde 0,01 a 0,05. En este aspecto tiene mucho peso la consideración de los efectos de segundo orden. Si observamos los resultados de la Tabla 18, realizados para un β=3, observamos que los voltajes obtenidos para 𝑡𝐿 igual a 0,05, 0,03 y 0,01 distaban mucho entre sí en el 47 análisis de primer orden. Sin embargo, los voltajes obtenidos en segundo orden mostraban valores mucho más similares. Esto último explica la superposición de las curvas que se observa en la gráfica en torno al β óptimo y su escasa distancia en valores fuera del óptimo. El cálculo lineal habría aportado curvas mucho más distantes entre sí. Figura 38: Influencia de β en el voltaje según 𝑡𝐿 (𝑡 𝑊=4) 3.2.5.3. Influencia de β, 𝒕 𝑾 y 𝒕𝑳 en análisis de 2º orden En este apartado se lleva a cabo el análisis de la respuesta del recolector piezoeléctrico midiendo la influencia de todas las relaciones de aspecto contempladas en este trabajo. Para ello se ha evaluado la respuesta del recolector para las razones de espesor-ancho, 0,1, 0,8 y 2, y tres valores de 𝑡𝐿 : 0,1; 0,05 y 0,02. Asimismo, para cada una de estas combinaciones se ha analizado la respuesta para una carga de 350 N y los siguientes valores de β: 0,1; 0,2; 0,3; 0,4; 0,5 y 1. Se han encontrado los resultados expuestos en la Figura 39. El voltaje se ha normalizado con respecto al voltaje máximo obtenido para 𝑡𝐿=0,02 y 𝑡 𝑊=2 según el análisis de primer orden, V0= 188590 V. Es evidente por la comparación de los resultados de primer orden y segundo orden que para la viga más esbelta es esencial considerar el análisis en segundo orden independientemente de la relación espesor-ancho. La diferencia entre las respuestas de voltaje obtenidas según los cálculos lineal y no lineal es grande para todos los casos de 𝑡 𝑊 e incrementa con 𝑡 𝑊. Como se ha explicado previamente, esto se debe a que un ancho menor provee un módulo de inercia de la sección resistente a flexión menor, provocando una mayor flecha. Consecuentemente, los desplazamientos se alejarán cada vez más de cumplir la hipótesis de pequeños desplazamientos. En cuanto a las vigas con menor 𝑡𝐿, se aprecia que para las dos menores relaciones espesor-ancho, la respuesta en primer orden y en segundo orden es muy similar. Las deformaciones no lineales tienen poca relevancia en estos casos. 3 Análisis de los resultados 48 Por otra parte, se puede observar que para la relación 𝑡 𝑊 más baja, la diferencia de voltaje producido para distintas longitudes es menor que para la 𝑡 𝑊 más alta, considerando la respuesta en segundo orden. Es decir, el voltaje obtenido en la viga con menor ancho es más susceptible a un cambio de longitud. Esto se debe a que el módulo de inercia de la sección resistente a flexión es menor porque el ancho es menor, de modo que la disminución de la rigidez provocada por un aumento de longitud es aún más notable. Así, la deformación incrementa sustancialmente y el voltaje producido también. Por último, se puede observar que el β óptimo cambia según la situación, salvo para 𝑡𝐿=0,1, donde este es 0,3 independientemente del espesor-ancho. Para el resto de los ratios de 𝑡𝐿 contemplados, el óptimo cambia en función de 𝑡 𝑊. Asimismo, para una misma relación 𝑡 𝑊, el óptimo cambia según 𝑡𝐿. En este aspecto posee mucha relevancia la consideración de los efectos de segundo orden, pues aunque en el análisis de primer orden los óptimos también cambian, no se sitúan en los mismos valores que para el caso de segundo orden. Esto arroja una conclusión relevante para el diseño del recolector: la combinación de los espesores de los materiales óptima depende tanto de los materiales seleccionados, como de la relaciones que guarda el espesor de la viga con el ancho y la longitud. 49 (a) (b) (c) (d) (e) (f) Figura 39: (a) Influencia de β en el voltaje según 𝑡𝐿 para 𝑡 𝑊=0,1 en 1er orden (b) Influencia de β en el voltaje según 𝑡𝐿 para 𝑡 𝑊=0,1 en 2º orden (c) Influencia de en el voltaje según 𝑡𝐿 para 𝑡 𝑊=0,8 en 1er orden (d) Influencia de en el voltaje según 𝑡𝐿 para 𝑡 𝑊=0,8 en 2º orden (e) Influencia de en el voltaje según 𝑡𝐿 para 𝑡 𝑊=2 en 1er orden (f) Influencia de en el voltaje según 𝑡𝐿 para 𝑡 𝑊=2 en 1er orden 3 Análisis de los resultados 50 3.3. Respuesta dinámica de vigas unimorfas En este apartado se exponen los resultados obtenidos de los análisis transitorios llevados a cabo, cuyo proceso de modelización se ha descrito en el apartado 2.5. 3.3.1. Influencia de la relación de aspecto t/L, factor de amortiguamiento ξ y los efectos de segundo orden En el apartado 3.2.3 se ha demostrado la influencia que la relación de aspecto 𝑡𝐿 posee sobre el comportamiento y el rendimiento de la viga piezoeléctrica ante una carga estática. En concreto, se ha demostrado que a menor relación 𝑡𝐿, mayor es el rendimiento del recolector. Asimismo, se ha comprobado que a mayor esbeltez de la viga, mayor influencia poseen los efectos de segundo orden, provocando que los resultados calculados en primer orden no describan correctamente el comportamiento del recolector. Del mismo modo, resulta interesante comprobar si estas conclusiones se mantienen para el análisis transitorio. Para comprobar la influencia de estos aspectos se ha obtenido la respuesta de la viga piezoeléctrica con 𝑡 𝑊=4, β=0,3 y 𝑡𝐿=0,1 y 0,03. El valor de la carga escalón aplicada ha sido de 700 N. Se ha empleado un factor de amortiguación de Rayleigh β𝑅=0,00001. Esto permite calcular el factor de amortiguamiento según la expresión: ξ =𝑐 𝑐𝑐=βR k 2√𝑘𝑚 . En esta expresión la rigidez y la masa variarán según las dimensiones y el régimen de funcionamiento que estemos evaluando, de modo que los factores de amortiguamiento para cada caso se exponen en la Tabla 21. Así, se representan en las figuras 40 y 41 la respuesta para 𝑡𝐿=0,1. El desplazamiento se ha normalizado según el espesor del material piezoeléctrico, tp=2 mm, y el voltaje se ha normalizado según el máximo voltaje obtenido con el cálculo de 1er orden, V0=171892 V. En primer lugar, cabe destacar que las respuestas dinámicas obtenidas se estabilizan en los valores de las respuestas estáticas correspondientes, una vez finalizada la fase transitoria. De este modo, se obtiene que la flecha estática remanente es mayor en el análisis de primer orden que en el análisis de segundo orden. Paralelamente, el voltaje final es mayor en primer orden que en segundo orden. En los resultados se observa un efecto relevante que no se contemplaba en el análisis estático: si bien el valor final de voltaje es el mismo, existe un periodo de tiempo en el que se alcanza un voltaje ligeramente superior. Es decir, la producción de voltaje del recolector piezoeléctrico es mayor de lo que arrojaban los análisis estáticos. A pesar de que la flecha en el análisis de primer orden es mucho mayor que aquella en el análisis de segundo orden, los valores finales de voltaje no difieren en la misma medida, teniendo valores mucho más similares. Esto se puede razonar observando la influencia de los efectos de segundo orden en β para las relaciones de aspecto 𝑡 𝑊=4 y 𝑡𝐿=0,1, expuesta en las figuras 35 y 36, la variación de β influye más en los efectos de segundo orden de la flecha que del voltaje. Por otra parte, se puede apreciar que para el valor de amortiguamiento proporcionado se producen varias fluctuaciones en la respuesta tanto en primer orden como en segundo orden, aunque el amortiguamiento será mayor en la respuesta de segundo orden, dado que este es proporcional a la rigidez y la respuesta en segundo orden es más rígida. Se puede apreciar que mientras que la respuesta en primer orden describe en su fase transitoria una sinusoide perfecta, el cálculo no lineal provoca que la respuesta de segundo orden sea más irregular. Además, la rigidización tensional que tiene lugar en la respuesta de segundo orden provoca un aumento en su frecuencia natural y por lo tanto, que las oscilaciones en la respuesta tengan un menor periodo. Se han ejecutado las mismas comprobaciones para 𝑡𝐿=0,03 y su resultado se expone en las figuras 42 y 43. La flecha está normalizada según el espesor del material piezoeléctrico, tp=2 mm, y el voltaje se ha normalizado según el máximo voltaje obtenido con el cálculo en primer orden, 1030000 V. Al igual que para el análisis estático, se verifica que mientras más esbelta sea la viga, mayor voltaje produce y mayor es la influencia de los efectos de segundo orden. Como consecuencia de su esbeltez, la viga fluctúa durante más tiempo y se aprecia que en el periodo de tiempo estudiado la respuesta aún no se ha estabilizado 51 en su valor estático. Asimismo, se vuelven a apreciar las irregularidades en las curvas de segundo orden. En este caso es muy evidente que el periodo de la vibración en segundo orden es menor que en primer orden. 𝛽𝑅 K (N/kg) ξ 𝑡𝐿=0,1 1er orden 0,00001 20245,26 0,06 2º orden 0,00001 37975,37 0,08 𝑡𝐿=0,03 1er orden 0,00001 552,5734 0,005 2º orden 0,00001 8625,47 0,02 Tabla 21: Factor de amortiguamiento (β=0,3, 𝑡 𝑊=4 𝑦 𝛽𝑅=0,00001) Figura 40: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,1 y 𝛽𝑅=0,00001) Figura 41: Influencia de los efectos de 2º orden en el voltaje en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,1 y 𝛽𝑅=0,00001) 3 Análisis de los resultados 52 Figura 42: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,03 y 𝛽𝑅= 0,00001) Figura 43: Influencia de los efectos de 2º orden en el voltaje en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,03 y 𝛽𝑅= 0,00001) A continuación, cabe cuestionarse la relevancia del factor de amortiguamiento en las conclusiones extraídas. Si bien el valor del parámetro de Rayleigh 𝛽𝑅 igual a 0,00001 aporta factores de amortiguamiento razonables, se comprobará la variación de la respuesta para 𝛽𝑅=0,01, dando lugar a unos factores de amortiguamiento excesivos, expuestos en la Tabla 22. Se han ejecutado los análisis con las mismas relaciones de aspecto que en el caso anterior (𝑡 𝑊=4, β=0,3 y 𝑡𝐿=0,1 y 0,03) y la misma carga, 700 N. Así, se representan en las figuras 53 44 y 45 la respuesta para 𝑡𝐿=0,1. El desplazamiento se ha normalizado según el espesor del material piezoeléctrico, tp=2 mm, y el voltaje se ha normalizado según el máximo voltaje obtenido con el cálculo de 2º orden, V0=160517 V. A diferencia del caso anterior, se comprueba que la respuesta está sobreamortiguada, es decir, no tiene lugar ni una oscilación completa antes de que la respuesta se estabilice en su valor estático. En el caso de la flecha, la curva se estabiliza de manera inmediata, mientras que para el voltaje, oscila ligeramente en su fase transitoria, de modo que el voltaje arrojado es mayor que el encontrado para la comprobación estática. Asimismo, se observa una mayor diferencia entre las curvas de primer y segundo orden en la flecha que en el voltaje, que de nuevo se explica consultando las gráficas 35 y 36. Cabe destacar que tanto el valor de la flecha máxima como el voltaje son menores en este caso sobreamortiguado que en el caso anterior. Se han ejecutado las mismas comprobaciones para 𝑡𝐿=0,03 y su resultado se expone en las figuras 42 y 43. La flecha está normalizada según el espesor del material piezoeléctrico, tp=2 mm, y el voltaje se ha normalizado según el máximo voltaje obtenido con el cálculo en primer orden, V0= 509975 V. Análogamente a la repuesta para 𝑡𝐿=0,1, se observa que existe una fluctuación temporal en el voltaje generado, previamente a la estabilización del valor del voltaje. Sin embargo, se puede comprobar que para esta viga más esbelta, el vértice de la gráfica es aún más pronunciado, existiendo una mayor diferencia entre el voltaje máximo alcanzando y el voltaje remanente. Este fenómeno tiene dos posibles explicaciones. La primera de ellas es que se trate de un error numérico producido en la integración de la respuesta. Si bien esta explicación parece descartable por el pequeño paso que se ha empleado en el cálculo (0.00000005 s), el valor del vértice es mucho mayor que el valor remanente de voltaje, lo que lo hace susceptible de ser incoherente. La segunda opción es que en efecto, la respuesta mostrada sea real. Esto aportaría una conclusión muy relevante para el diseño del recolector: para conocer el valor máximo de energía eléctrica generada ante una carga intermitente o tipo escalón sería necesario considerar la respuesta transitoria. En efecto, al evaluar la respuesta transitoria, se conoce que se genera un voltaje mayor que aquel obtenido para la respuesta estática, de modo que esta última subestima el rendimiento del recolector. Asimismo, al tratarse de un voltaje casi 3 veces mayor que el que aporta el análisis estático, diseñar la instalación eléctrica para el valor de diseño estático puede provocar un fallo de la instalación eléctrica por efecto Joule, dado que el diámetro de cable necesario será mayor del existente. 𝛽𝑅 K (N/kg) ξ 𝑡𝐿=0,1 1er orden 0,01 20245,26 55,23 2º orden 0,01 37975,37 75,63 𝑡𝐿=0,03 1er orden 0,01 552,5734 5 2º orden 0,01 8625,47 19,74 Tabla 22: Factor de amortiguamiento (β=0,3, 𝑡 𝑊=4 𝑦 𝛽𝑅=0,01) 3 Análisis de los resultados 54 Figura 44: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,1 y y 𝛽𝑅=0,01) Figura 45: Influencia de los efectos de 2º orden en la flecha en análisis transitorio (𝛽=0,3,𝑡 𝑊=4, 𝑡𝐿=0,1 y y 𝛽𝑅=0,01) 61 ANEXO 1: CÓDIGOS DE MATLAB • Modelo analítico clc clear all format long %introduzco propiedades materiales del material 1, base E1=195e9; %Young %introduzco propiedades materiales del material 2, piezoeléctrico E2=81e9; perm33=5.82e-9; d31=-124e-12; s11=1/E2; %introduzco parámetros geométricos L=0.1; t1=0.6e-2; %grosor base t2=0.2e-2; %grosor piezoeléctrico W=0.36e-2; %parámetros relativos alpha=E1/E2; beta=t1/t2; %valor de la fuerza F=-350; %cálculo del desplazamiento máximo en el caso mecánico (expresionesq %extraídas de “Modelling and analysis of a piezoelectric unimorph %cantilever for energy harvesting application") wmax=-(4*F*L^3*(1+alpha*beta))/(W*E2*t2^3*(1+4*alpha*beta+6*alpha*beta^2+... 4*alpha*beta^3+alpha^2*beta^4)) %cálculo del voltaje producido B=(t2^4/((alpha^2)*(E2^2)))*(1+4*alpha*beta+6*alpha*beta^2+... 4*alpha*beta^3+alpha^2*beta^4); %B ESTÁ BIEN aux1=(3*beta^2*(1+beta)^2*t2^4*s11^2)/((1+alpha*beta)*B) aux=1+(aux1-1)*(d31^2/(s11*perm33)); Q=(-3*(d31*beta*(1+beta)*s11^2*t2^2*L^2)*F)/(alpha*B) C=(perm33*W*L/t2)*aux V=Q/C Anexo 1: Códigos de MATLAB 62 62 • Influencia de β en el voltaje en función de α clc close all clear all format long %introduzco propiedades materiales del material 1, base E1=195e9; %Young %introduzco propiedades materiales del material 2, piezoeléctrico E2=81e9; perm33=5.82e-9; d31=-124e-12;% s11=1/E2; %introduzco parámetros geométricos t2=0.2e-2; %grosor piezoeléctrico %parámetros relativos alpha1=E1/E2; rel_t_W=4; %ratio t/W rel_t_L=0.1; %valor de la fuerza F=-350; %cálculo del desplazamiento máximo en el caso mecánico (expresionesq %extraídas de “Modelling and analysis of a piezoelectric unimorph %cantilever for energy harvesting application") %cálculo del voltaje producido voltaje2=[]; carga2=[]; beta=0.1:0.1:4 alpha=[0.25 0.5 1 alpha1] for j=1:length(alpha) voltaje1=[]; carga1=[]; for i=1:length(beta) t1=beta(i)*t2; t=t1+t2; W=t/rel_t_W; L=t/rel_t_L; B=(t2^4/((alpha(j)^2)*(E2^2)))*(1+4*alpha(j)*beta(i)+6*alpha(j)*(beta(i))^2+... 4*alpha(j)*(beta(i))^3+(alpha(j))^2*(beta(i))^4); aux1=(3*(beta(i))^2*(1+beta(i))^2*t2^4*s11^2)/((1+alpha(j)*beta(i))*B); aux=1+(aux1-1)*(d31^2/(s11*perm33)); Q=(-3*(d31*beta(i)*(1+beta(i))*s11^2*t2^2*L^2)*F)/(alpha(j)*B); 63 C=(perm33*W*L/t2)*aux; V=Q/C; voltaje1=[voltaje1 V]; carga1=[carga1 Q]; end voltaje2=[voltaje2; voltaje1]; % carga2=[carga2; carga1]; end Vmax=max(abs(voltaje2), [], 'all'); grid ON plot(beta,-(voltaje2(1,1:end))/Vmax,'-',beta,-(voltaje2(2,1:end))/Vmax,'-x',beta,- (voltaje2(3,1:end))/Vmax,'-o',beta,-(voltaje2(4,1:end))/Vmax,'-+') grid on axis square legend('\alpha=0.25','\alpha=0.5','\alpha=1','\alpha=2.4','fontsize',15) xlabel('\beta (t_b/t_p)','fontsize',15) ylabel('V/V_o','fontsize',15) title('Influencia de la relación entre espesores de materiales','fontsize',15) Anexo 2: Códigos de ANSYS 64 64 ANEXO 2: CÓDIGOS DE ANSYS • Modelo 3D estático en análisis de primer orden finish /clear !------------------- ! Preprocesador !------------------- /PREP7 ET,1,SOLID186 ET,2,SOLID226 KEYOPT,2,1,1001 !DEFINO PROPIEDADES MATERIALES !------------------------------ !* Prop. mat.1 !------------- MP,EX,1,195E9 MP,PRXY,1,0.3 MP,DENS,1,7750 !------------ !*Prop. mat.2 !------------- MP,DENS,2,7495 EMUNIT,EPZRO,8.854E-12 !DEFINO PERMITIVIDAD EN EL VACÍO *SET,PER11,763 *SET,PER22,763 *SET,PER33,658 MP,PERX,2,PER11 MP,PERY,2,PER22 MP,PERZ,2,PER33 !DEFINO MATRIZ DE RIGIDEZ 65 *SET,C11, 139E9 *SET,C12,77.8E9 *SET,C13,74.3E9 *SET,C22, 139E9 *SET,C23,74.3E9 *SET,C33, 115E9 *SET,C44,25.6E9 *SET,C55,25.6E9 *SET,C66,30.6E9 !* TB,ANEL,2,1,21,0 TBTEMP,0 TBDATA,,C11,C12,C13,0,0,0 TBDATA,,C22,C23,0,0,0,C33 TBDATA,,0,0,0,C44,0,0 TBDATA,,C55,0,C66,,, TB,PIEZ,2,,,0 !MATRIZ DE COEFICIENTES PIEZOELECTRICOS *SET,E11,0 *SET,E12,0 *SET,E13,-5.2 *SET,E21,0 *SET,E22,0 *SET,E23,-5.2 *SET,E31,0 *SET,E32,0 *SET,E33,15.1 *SET,E41,0 *SET,E42,12.7 *SET,E43,0 *SET,E51,12.7 *SET,E52,0 *SET,E53,0 *SET,E61,0 *SET,E62,0 *SET,E63,0 TBDATA,1,E11,E12,E13,E21,E22,E23 TBDATA,7,E31,E32,E33,E41,E42,E43 Anexo 2: Códigos de ANSYS 66 66 TBDATA,13,E51,E52,E53,E61,E62,E63 !------------ ! Geometria !------------ Rat_t_L = 0.1 Rat_t_W = 4 beta_espesor=0.1 *SET,TP,2E-3 *SET,TB,TP*beta_espesor *SET,T,TP+TB *SET,W, T/Rat_t_W *SET,L,T/Rat_t_L K,1, K,2,L,0,0, K,3,L,W,0 K,4,0,W, K,5,0,0,TP K,6,L,0,TP K,7,L,W,TP K,8,0,W,TP K,9,0,0,-TB K,10,L,0,-TB K,11,L,W,-TB K,12,0,W,-TB ! volumen 1 (base) V,1,2,3,4,9,10,11,12 ! volumen 2 (piezoelectrico) V,1,2,3,4,5,6,7,8 !MALLADO !-------- 67 NDIVX=120*0.1/Rat_t_L NDIVY=8 vsel,S,,,1 VATT,1,,1 allsel vsel,S,,,2 VATT,2,,2 allsel LESIZE,6,,,NDIVX LESIZE,8,,,NDIVY LESIZE,5,,,4 LESIZE,13,,,8 MSHKEY,1 !MALLADO MAPPED MSHAPE,0,3D !0 ES MALLADO HAEXAÉDRICO O CUADRILÁTERO Y 3D VMESH,ALL ! APLICO EMPOTRAMIENTO !-------------------- NSEL,S,LOC,X,0 D,ALL,UX,0 NSEL,R,LOC,Y,W/2 D,ALL,UY,0 NSEL,S,LOC,X,0 NSEL,R,LOC,Z,-TB D,ALL,UZ,0 ASEL,S,LOC,Z,0 DA,ALL,VOLT,0 !IMPONGO VOLTAJE NULO EN LA BASE DEL PIEZOELÉCTRICO FINISH !------------------- ! Resolucion !------------------- /SOLU Anexo 2: Códigos de ANSYS 68 68 F=350 NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 F,ALL,FZ,-F/(NDIVY*2) !elementos cuadraticos "NDIVY*2" allsel NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 NSEL,R,LOC,Y,0 F,ALL,FZ,-(F/(NDIVY*2))/2 allsel NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 NSEL,R,LOC,Y,W F,ALL,FZ,-(F/(NDIVY*2))/2 allsel ALLSEL !IMPORTANTE SOLVE FINISH !------------------- ! Post-Procesado !------------------- /POST1 PLNSOL, U,Z, 0,1.0 PRNSOL,VOLT PRNSOL,U,Z • Modelo 2D estático en análisis de primer orden 69 finish /clear !------------------- ! Preprocesador !------------------- /PREP7 ET,1,PLANE183 KEYOPT,1,3,0 ET,2,PLANE223 KEYOPT,2,1,1001 KEYOPT,2,3,0 !-------------------- !* Prop. mat.1 !-------------------- MP,EX,1,195E9 MP,PRXY,1,0.3 MP,DENS,1,7750 !------------------- !*Prop. mat.2 !------------------- MP,DENS,2,7495 EMUNIT,EPZRO,8.854E-12 !DEFINO PERMITIVIDAD EN EL VACÍO *SET,PER11,763 *SET,PER22,658 *SET,PER33,763 MP,PERX,2,PER11 MP,PERY,2,PER22 MP,PERZ,2,PER33 !DEFINO MATRIZ DE RIGIDEZ *SET,C11,139E9 *SET,C12,74.3E9 *SET,C13,77.8E9 *SET,C14,0 *SET,C21,74.3E9 Anexo 2: Códigos de ANSYS 70 70 *SET,C22,115E9 *SET,C23,74.3E9 *SET,C24,0 *SET,C31,77.8E9 *SET,C32,74.3E9 *SET,C33,139E9 *SET,C34,0 *SET,C41,0 *SET,C42,0 *SET,C43,0 *SET,C44,30.6E9 TB,ANEL,2,1,21,0 TBTEMP,0 TBDATA,,C11,C12,C13,0,0,0 TBDATA,,C22,C23,0,0,0,C33 TBDATA,,0,0,0,C44,0,0 TB,PIEZ,2,,,0 *SET,E11,0 *SET,E12,-5.2 *SET,E21,0 *SET,E22,15.1 *SET,E31,0 *SET,E32,0 *SET,E41,0 *SET,E42,0 TBDATA,1,E11,E12,0,E21,E22,0 TBDATA,7,E31,E32,0,E41,E42,0 BETAD, 0.01 !------------- ! Geometria !------------- Rat_t_L = 0.1 Rat_t_W=4 beta_espesor=0.1 *SET,TP,2E-3 77 !------------------- ! Resolucion !------------------- /SOLU F=350 NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 F,ALL,FZ,-F/(NDIVY*2) !elementos cuadraticos "NDIVY*2" allsel NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 NSEL,R,LOC,Y,0 F,ALL,FZ,-(F/(NDIVY*2))/2 allsel NSEL,S,LOC,X,L NSEL,R,LOC,Z,0 NSEL,R,LOC,Y,W F,ALL,FZ,-(F/(NDIVY*2))/2 allsel ALLSEL NLGEOM,1 SOLVE FINISH !------------------- ! Post-Procesado !------------------- /POST1 PLNSOL, U,Z, 0,1.0 Anexo 2: Códigos de ANSYS 78 78 PRNSOL,VOLT PRNSOL,U,Z • Modelo 2D estático en análisis de segundo orden finish /clear !------------------- ! Preprocesador !------------------- /PREP7 ET,1,PLANE183 KEYOPT,1,3,0 ET,2,PLANE223 KEYOPT,2,1,1001 KEYOPT,2,3,0 !------------------------------ !DEFINO PROPIEDADES MATERIALES !------------------------------ !* Prop. mat.1 !------------- MP,EX,1,195E9 MP,PRXY,1,0.3 MP,DENS,1,7750 !*Prop. mat.2 !------------- MP,DENS,2,7495 EMUNIT,EPZRO,8.854E-12 !DEFINO PERMITIVIDAD EN EL VACÍO *SET,PER11,763 *SET,PER22,658 *SET,PER33,763 MP,PERX,2,PER11 MP,PERY,2,PER22 79 MP,PERZ,2,PER33 !DEFINO MATRIZ DE RIGIDEZ *SET,C11,139E9 *SET,C12,74.3E9 *SET,C13,77.8E9 *SET,C14,0 *SET,C21,74.3E9 *SET,C22,115E9 *SET,C23,74.3E9 *SET,C24,0 *SET,C31,77.8E9 *SET,C32,74.3E9 *SET,C33,139E9 *SET,C34,0 *SET,C41,0 *SET,C42,0 *SET,C43,0 *SET,C44,30.6E9 TB,ANEL,2,1,21,0 TBTEMP,0 TBDATA,,C11,C12,C13,0,0,0 TBDATA,,C22,C23,0,0,0,C33 TBDATA,,0,0,0,C44,0,0 TB,PIEZ,2,,,0 *SET,E11,0 *SET,E12,-5.2 *SET,E21,0 *SET,E22,15.1 *SET,E31,0 *SET,E32,0 *SET,E41,0 *SET,E42,0 TBDATA,1,E11,E12,0,E21,E22,0 TBDATA,7,E31,E32,0,E41,E42,0 ! Geometria !------------ Anexo 2: Códigos de ANSYS 80 80 Rat_t_L = 0.01 Rat_t_W=4 beta_espesor=0.1 *SET,TP,2E-3 *SET,TB,TP*beta_espesor *SET,T,TP+TB *SET,W, T/Rat_t_W *SET,L,T/Rat_t_L K,1, K,2,L,0 K,3,L,tp K,4,0,tp K,5,0,-tb K,6,L,-tb !area 1 base A,1,2,6,5 !area 2 (piezoelectrico) A,1,2,3,4 !MALLADO !-------- NDIVX=120*0.1/Rat_t_L NDIVY=4 asel,S,,,1 AATT,1,,1 allsel asel,S,,,2 AATT,2,,2 allsel LPLOT LESIZE,1,,,NDIVX LESIZE,2,,,NDIVY*2 LESIZE,5,,,NDIVY 81 MSHKEY,1 !MALLADO MAPPED MSHAPE,0,2D !0 ES MALLADO HAEXAÉDRICO O CUADRILÁTERO Y 2D AMESH,ALL !---------------------- !APLICO EMPOTRAMIENTO !------------------------ NSEL,S,LOC,X,0 D,ALL,UX,0 allsel NSEL,S,LOC,X,0 NSEL,R,LOC,Y,-TB D,ALL,UY,0 allsel NSEL,S,LOC,Y,0 D,ALL,VOLT,0 !IMPONGO VOLTAJE NULO EN LA BASE DEL PIEZOELÉCTRICO allsel !------------------------ ! Resolución: !------------------------ /SOLU F = 350 NSEL,S,LOC,X,L,L+L/NDIVX NSEL,R,LOC,Y,0 F,ALL,FY,-F/W allsel NLGEOM,1 solve Anexo 2: Códigos de ANSYS 82 82 FINISH !------------------------ ! Postprocesado: !------------------------ /POST1 PLNSOL, U,Y, 0,1.0 PRNSOL,U,Y PRNSOL,VOLT • Modelo 2D dinámico en análisis de segundo orden finish /clear !------------------------- ! Preprocesador !------------------- /PREP7 ET,1,PLANE183 KEYOPT,1,3,0 ET,2,PLANE223 KEYOPT,2,1,1001 KEYOPT,2,3,0 !----------------------------- !DEFINO PROPIEDADES MATERIALES !------------------------------ !* Prop. mat.1 !------------- MP,EX,1,195E9 MP,PRXY,1,0.3 MP,DENS,1,7750 !------------ !*Prop. mat.2 !------------- MP,DENS,2,7495 83 EMUNIT,EPZRO,8.854E-12 !DEFINO PERMITIVIDAD EN EL VACÍO *SET,PER11,763 *SET,PER22,658 *SET,PER33,763 MP,PERX,2,PER11 MP,PERY,2,PER22 MP,PERZ,2,PER33 !DEFINO MATRIZ DE RIGIDEZ *SET,C11,139E9 *SET,C12,74.3E9 *SET,C13,77.8E9 *SET,C14,0 *SET,C21,74.3E9 *SET,C22,115E9 *SET,C23,74.3E9 *SET,C24,0 *SET,C31,77.8E9 *SET,C32,74.3E9 *SET,C33,139E9 *SET,C34,0 *SET,C41,0 *SET,C42,0 *SET,C43,0 *SET,C44,30.6E9 TB,ANEL,2,1,21,0 TBTEMP,0 TBDATA,,C11,C12,C13,0,0,0 TBDATA,,C22,C23,0,0,0,C33 TBDATA,,0,0,0,C44,0,0 TB,PIEZ,2,,,0 *SET,E11,0 *SET,E12,-5.2 *SET,E21,0 *SET,E22,15.1 *SET,E31,0 *SET,E32,0 Anexo 2: Códigos de ANSYS 84 84 *SET,E41,0 *SET,E42,0 TBDATA,1,E11,E12,0,E21,E22,0 TBDATA,7,E31,E32,0,E41,E42,0 ! Geometria !------------ Rat_t_L = 0.1 Rat_t_W=4 beta_espesor=3 *SET,TP,2E-3 *SET,TB,TP*beta_espesor *SET,T,TP+TB *SET,W, T/Rat_t_W ! 0.36E-2 *SET,L,T/Rat_t_L K,1, K,2,L,0 K,3,L,tp K,4,0,tp K,5,0,-tb K,6,L,-tb !area 1 base A,1,2,6,5 !area 2 (piezoelectrico) A,1,2,3,4 !MALLADO !-------- *SET,NDIVX,120*0.1/Rat_t_L *SET,NDIVY,4 asel,S,,,1 AATT,1,,1 allsel asel,S,,,2 85 AATT,2,,2 allsel LPLOT LESIZE,1,,,NDIVX LESIZE,2,,,NDIVY*2 LESIZE,5,,,NDIVY MSHKEY,1 !MALLADO MAPPED MSHAPE,0,2D !0 ES MALLADO HAEXAÉDRICO O CUADRILÁTERO Y 2D AMESH,ALL !APLICO EMPOTRAMIENTO NSEL,S,LOC,X,0 D,ALL,UX,0 allsel NSEL,S,LOC,X,0 NSEL,R,LOC,Y,-TB D,ALL,UY,0 allsel NSEL,S,LOC,Y,0 D,ALL,VOLT,0 !IMPONGO VOLTAJE NULO EN LA BASE DEL PIEZOELÉCTRICO allsel !------------------------ ! Solución: !------------------------ /SOLU ANTYPE,4 NLGEOM,1 TRNOPT,FULL Anexo 2: Códigos de ANSYS 86 86 LUMPM,0 OUTRES,ERASE OUTRES,BASIC,1 deltaTime=0.0001 Betad_fact = 0.01 ! EN EL INSTANTE T=0s !--------------------- TIME,0.0000000001 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,1, ! EN EL INSTANTE T=0.001s !------------------------ *SET,F,700 NSEL,S,LOC,X,L,L+L/NDIVX NSEL,R,LOC,Y,0 F,ALL,FY,-F/W allsel TIME,0.001 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,2, ! EN EL INSTANTE T=0.5s !------------------------ TIME,0.1 AUTOTS,-1 93 NSEL,R,LOC,Y,0 F,ALL,FY,-F/W allsel TIME,0.001 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,2, ! EN EL INSTANTE T=0.005s !------------------------ TIME,0.005 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,3, LSSOLVE,1,3,1 FINISH !------------------------ ! Postprocesado: Anexo 2: Códigos de ANSYS 94 94 !------------------------ /POST26 NSOL,2,2,U,Y, UY_2 STORE,MERGE PRVAR,2 PLVAR,2 • Modelo 2D transitorio en análisis de segundo orden finish /clear !------------------------- ! Preprocesador !------------------- /PREP7 ET,1,PLANE183 KEYOPT,1,3,0 ET,2,PLANE223 KEYOPT,2,1,1001 KEYOPT,2,3,0 !DEFINO PROPIEDADES MATERIALES !------------------------------ !* Prop. mat.1 !------------- MP,EX,1,195E9 MP,PRXY,1,0.3 MP,DENS,1,7750 !*Prop. mat.2 95 !------------- MP,DENS,2,7495 EMUNIT,EPZRO,8.854E-12 !DEFINO PERMITIVIDAD EN EL VACÍO *SET,PER11,763 *SET,PER22,658 *SET,PER33,763 MP,PERX,2,PER11 MP,PERY,2,PER22 MP,PERZ,2,PER33 !DEFINO MATRIZ DE RIGIDEZ *SET,C11,139E9 *SET,C12,74.3E9 *SET,C13,77.8E9 *SET,C14,0 *SET,C21,74.3E9 *SET,C22,115E9 *SET,C23,74.3E9 *SET,C24,0 *SET,C31,77.8E9 *SET,C32,74.3E9 *SET,C33,139E9 *SET,C34,0 *SET,C41,0 *SET,C42,0 *SET,C43,0 *SET,C44,30.6E9 TB,ANEL,2,1,21,0 Anexo 2: Códigos de ANSYS 96 96 TBTEMP,0 TBDATA,,C11,C12,C13,0,0,0 TBDATA,,C22,C23,0,0,0,C33 TBDATA,,0,0,0,C44,0,0 TB,PIEZ,2,,,0 *SET,E11,0 *SET,E12,-5.2 *SET,E21,0 *SET,E22,15.1 *SET,E31,0 *SET,E32,-5.2 *SET,E41,0 *SET,E42,0 TBDATA,1,E11,E12,0,E21,E22,0 TBDATA,7,E31,E32,0,E41,E42,0 ! Geometria !------------ Rat_t_L = 0.03 Rat_t_W=4 beta_espesor=0.3 *SET,TP,2E-3 *SET,TB,TP*beta_espesor *SET,T,TP+TB *SET,W, T/Rat_t_W ! 0.36E-2 *SET,L,T/Rat_t_L 97 K,1, K,2,L,0 K,3,L,tp K,4,0,tp K,5,0,-tb K,6,L,-tb !area 1 base A,1,2,6,5 !area 2 (piezoelectrico) A,1,2,3,4 !MALLADO !-------- *SET,NDIVX,50*0.1/Rat_t_L *SET,NDIVY,4 asel,S,,,1 AATT,1,,1 allsel asel,S,,,2 AATT,2,,2 !AATT,1,,1 allsel LPLOT LESIZE,1,,,NDIVX LESIZE,2,,,NDIVY*2 LESIZE,5,,,NDIVY MSHKEY,1 !MALLADO MAPPED Anexo 2: Códigos de ANSYS 98 98 MSHAPE,0,2D !0 ES MALLADO HAEXAÉDRICO O CUADRILÁTERO Y 2D AMESH,ALL !APLICO EMPOTRAMIENTO NSEL,S,LOC,X,0 D,ALL,UX,0 allsel NSEL,S,LOC,X,0 NSEL,R,LOC,Y,-TB D,ALL,UY,0 allsel NSEL,S,LOC,Y,0 D,ALL,VOLT,0 !IMPONGO VOLTAJE NULO EN LA BASE DEL PIEZOELÉCTRICO allsel !------------------------ ! Solución: !------------------------ /SOLU ANTYPE,4 NLGEOM,1 !NLGEOM,0 !PSTRES,1 NSUBST,5,0,0 TRNOPT,FULL 99 LUMPM,0 OUTRES,ERASE OUTRES,BASIC,1 deltaTime=0.00000005 Betad_fact = 0.01 ! EN EL INSTANTE T=0s !--------------------- TIME,0.0000000001 AUTOTS,-1 DELTIM,0.001, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,1, ! EN EL INSTANTE T=0.001s !------------------------ *SET,F,350*2 NSEL,S,LOC,X,L,L+L/NDIVX NSEL,R,LOC,Y,0 F,ALL,FY,-F/W allsel TIME,0.001 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 Anexo 2: Códigos de ANSYS 100 100 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,2, ! EN EL INSTANTE T=0.015s !------------------------ TIME,0.015 AUTOTS,-1 DELTIM,deltaTime, , , ALPHAD,0 BETAD,Betad_fact KBC,0 TSRES,ERASE LSWRITE,3, LSSOLVE,1,3,1 FINISH !------------------------ ! Postprocesado: !------------------------ /POST26 NSOL,2,2,U,Y, UY_2 STORE,MERGE PRVAR,2 PLVAR,2