Análisis numérico de uniones
abulonadas en materiales compuestos
Cristhian David Rodríguez Flores
Tesis de Maestría en Ciencias de la Ingeniería – Mención Aeroespacial
Universidad Nacional de Córdoba – Instituto Universitario Aeronáutico
Director: Dr. Ing. Walter B. Castelló
Co Director: MSc. Ing. Julio C. Massa
Noviembre de 2014
ANÁLISIS NUMÉRICO DE UNIONES ABULONADAS
EN MATERIALES COMPUESTOS
Resumen. La combinación de bajo peso y alta resistencia que brindan los materiales
compuestos, los posiciona convenientemente en aplicaciones estructurales exigentes, sin
embargo es conocido que la vinculación o unión entre elementos estructurales fabricados
en compuestos es un punto crítico. En este trabajo se desarrolla un modelo numérico
para estudiar una unión abulonada de solape simple de dos placas de material compuesto.
El análisis del problema se hace empleando un código de elementos finitos y el modelo
desarrollado contempla el múltiple contacto entre las placas de compuesto, el cuerpo del
bulón, la tuerca y las arandelas. Las placas se consideran fabricadas con varias capas de
fibra de carbono unidireccional en una matriz de epoxi, las arandelas son de acero, mientras
que el bulón y la tuerca son de titanio. Los resultados de este análisis computacional
concuerdan muy bien con los obtenidos experimentalmente por otros autores, no solo en
cuanto a nivel tensional sino también en las deformaciones y las distorsiones mostradas
por las placas. A partir de los resultados satisfactorios obtenidos con el modelo de un
solo bulón, se extendió el estudio a la simulación numérica del comportamiento de otras
configuraciones de abulonado cuyos resultados muestran buena concordancia con
datos obtenidos experimentalmente por otros autores en situaciones similares.
NUMERICAL ANALYSIS OF COMPOSITES BOLTED JOINTS
Abstract. It is well known that combination of low weight and high strength of composite
materials is a great advantage while designing demanding structural applications; however
it is also well known that connections between structural elements made of composite
are a critical issue. In this thesis a numerical model is developed to study a single lap
bolted union of two composite plates. The problem is addressed using a commercial
finite element code and developing a model that takes into account multiple contacts
between the plates of composite, the body of the bolt, the nut and the washers. The plates
contain several layers of unidirectional carbon fiber in an epoxy matrix considered as
sub-laminates; the washers are made of steel, while the bolt and nut are constructed of
titanium. The results of the computational analysis agree very well with those obtained
experimentally by other authors, not only in terms of stress level but also on the deformations and distortions shown by the plates. Based on the satisfactory results obtained
with the model of a single bolt, the work continued with numerical simulations of the
behavior of other configurations of bolted unions. The results of the simulation in those
cases also showed very good agreement with experimental data obtained by other authors
for similar situations.
i
ii
A Xime por estar siempre ahí, con su apoyo, amor y comprensión
en los días tristes y alegres, soleados y lluviosos
A Fabiola, mi madre por su amor, apoyo y confianza permanente,
rompiendo la barrera de la distancia
A Reinaldo, mi padre, y a Gaby, María y Raquel, mis hermanas,
por su amor y su fe
iii
iv
Agradecimientos
Cuando inicié esta etapa, lejos de mi patria y de mi familia, un país generoso, grande
y amable me acogió; con gente solidaria, confiable y cariñosa que supo tenderme una mano
cuando lo necesité: Argentina, gracias, te llevo atesorada en mi corazón.
Al regresar a los libros luego de casi una década, es inevitable encontrarse con que la
ciencia ha avanzado de modo que estás perdido en algún punto de ese avance. Es entonces
cuando se necesita la guía de quienes han seguido ese avance y han hecho de su vida una
religión de servicio a la formación de los jóvenes, capitalizando en su forma de ser, paciencia,
sabiduría y demás virtudes que los caracterizan como verdaderos maestros. Mi admiración
y gratitud para los que fueron los guías en esta hermosa etapa de nuevo aprendizaje: José
Inaudi, Gustavo Scarpin, Walkiria Schultz, Patricia Silvetti, Andrea Costa, Mario D’errico,
Gustavo Wolfman y Eduardo Zapico. Una especial mención merecen dos excelentes
maestros: Carlos Sacco y Sergio Elaskar, siempre con una palabra de aliento en el camino,
gracias por estar prestos a brindarme su ayuda y apoyo incondicional más allá de lo académico.
Gracias también a mis compañeros de clase: Luis Soria, Carlos Contreras, Andrés Cimino,
Dino Antonelli, Mario Rabale, José Alé, Marco Gómez y a todos los que compartieron este
interesante camino, que el éxito los cobije siempre. A mis amigos presentes: Juan Carlos,
Angelita, Ivanna, Pablo, Jeremías, Josué, Verónica, Enrique, Jorge, Luis L., Andrés, Luis
Germán, Gaby, Paola, Andrea C., Fernanda, Liz, Carlita y la pequeña Rocío, gracias por
su apoyo. A Doña. Marina, mi gratitud y mi cariño por su forma de ser, su preocupación,
solidaridad y por hacerme sentir parte de su familia, compartiendo hermosos momentos. A
mis amigos de siempre: Fredy, Paúl, Edward y Fátima gracias por su aliento permanente
desde el lugar donde se encuentren.
A Marco, Raúl, Ramiro, Raquel Alicia, Carlota, Silvia, Marco Jr., Tania, Rodrigo, Susy,
Sophia, y a todos quienes conforman mi familia, su apoyo siempre ha sido, es y será importante.
A mi padre Reinaldo, a Gaby, María y Raquel mis hermanas, gracias por su cariño,
confianza y preocupación que impulsan siempre mi caminar. A mi madre, todo mi amor,
por todos los sacrificios hechos, porque siempre ha confiado en mí, incluso cuando yo
mismo me he perdido la fe.
A Xime, todo el amor del mundo, por estar conmigo siempre, por darme su amor,
comprensión, ternura, apoyo incondicional y poner su hombro para sacar adelante este
proyecto y nuestro hogar.
A mis directores, Walter Castelló y Julio Massa, mi gratitud eterna y mi admiración por
llevar adelante el reto de dirigir este trabajo, superando las barreras físicas de la distancia,
por su disposición a contestar mis inquietudes y enfrentar los altibajos de este proyecto, por
buscar soluciones para salir de los problemas y por sacarme adelante.
Finalmente, mi gratitud a mi Ecuador y la Fuerza Aérea Ecuatoriana por creer en mí, a
la Fuerza Aérea Argentina, la UNC y el IUA por su apoyo y prestancia. Gracias finalmente
a mis compañeros de trabajo Galo Álvarez, José Salcedo, Raúl Banderas, Hernán Salazar y
a todos quienes confiaron en que era posible salir adelante con este proyecto y que hoy por
la fragilidad de la memoria no constan en estos párrafos.
Cristhian Rodríguez
v
vi
INDICE
Página
Capítulo 1: INTRODUCCIÓN
1.1 Motivación............................................................................................................................. 1
1.2 Justificación del proyecto...................................................................................................... 3
1.3 Objetivos................................................................................................................................ 4
1.4 Alcance del proyecto............................................................................................................. 4
1.5 Metodología de análisis......................................................................................................... 5
1.6 Contenido de la tesis.............................................................................................................. 6
Capítulo 2: MECÁNICA DE MATERIALES COMPUESTOS
2.1 Introducción............................................................................................................................ 7
2.2 Compuestos laminados........................................................................................................... 7
2.2.1 Secuencia del laminado............................................................................................... 8
2.3 Deformaciones y tensiones en compuestos........................................................................... 9
2.4 Matriz de rigidez de un laminado delgado........................................................................... 13
2.4.1 Matriz de rigidez del laminado usado en la tesis...................................................... 20
Capítulo 3: MODELACIÓN NUMÉRICA DE COMPUESTOS
3.1 Introducción.......................................................................................................................... 23
3.2 Estado actual del análisis numérico de materiales compuestos............................................23
3.3 Análisis por elementos finitos.............................................................................................. 27
3.4 Ecuación de gobierno........................................................................................................... 29
3.5 Selección del elemento finito y de las funciones de forma.................................................. 30
Capítulo 4 ANÁLISIS NUMÉRICO
DE UNIONES DE MATERIAL COMPUESTO ABULONADAS
4.1 Introducción.......................................................................................................................... 35
4.2 Descripción del problema..................................................................................................... 36
4.3 Modelo de elementos finitos................................................................................................. 37
4.3.1 Malla de elementos finitos......................................................................................... 37
4.3.2 Tipo de elementos utilizados..................................................................................... 39
4.3.3 Condiciones de borde y de carga............................................................................... 39
4.3.4 Modelo constitutivo del material................................................................................41
vii
4.4 Validación numérica y experimental del modelo propuesto................................................ 42
4.4.1 Deformaciones en la superficie de la placa............................................................... 43
4.4.2 Rigidez y deformaciones de la unión......................................................................... 48
4.4.3 Desplazamiento del plano medio de la placa............................................................. 48
4.4.4 Distribución de tensiones........................................................................................... 50
4.4.5 Deformación de la junta............................................................................................ 53
4.5 Resultados numéricos para uniones abulonadas dobles....................................................... 54
4.5.1 Modelo de dos bulones axiales.................................................................................. 54
4.5.1.1 Geometría de la probeta con dos bulones axiales......................................... 55
4.5.1.2 Deformaciones en la superficie de la placa y rigidez de la unión................. 56
4.5.1.3 Distribución de tensiones en la unión con dos bulones axiales..................... 56
4.5.1.4 Distorsión del plano medio de la placa para dos bulones axiales.................. 58
4.5.1.5 Deformación de la unión con dos bulones axiales.........................................59
4.5.2 Modelo de dos bulones transversales........................................................................ 59
4.5.2.1 Geometría de la probeta con dos bulones transversales................................ 60
4.5.2.2 Deformaciones en la superficie de la placa y rigidez de la unión................. 61
4.5.2.3 Distribución de tensiones en la unión con dos bulones transversales........... 61
4.5.2.4 Distorsión del plano medio de la placa para dos bulones transversales........ 63
4.5.2.5 Deformada final de la unión con dos bulones transversales......................... 64
Capítulo 5: CONCLUSIONES
5.1 Introducción......................................................................................................................... 65
5.2 Conclusiones........................................................................................................................ 66
5.3 Contribuciones del presente trabajo..................................................................................... 67
5.4 Líneas de trabajos futuros.................................................................................................... 68
ANEXO: MATRICES DE RIGIDEZ DE LAMINADOS Y SUBLAMINADOS
1
Introducción........................................................................................................................... 69
2
Matrices de transformación tensión-deformación.................................................................. 69
2.1 Transformación de tensiones..........................................................................................70
2.2 Transformación de deformaciones................................................................................. 72
3
Condición de tensión plana.....................................................................................................73
4
Obtención de la matriz de flexibilidad de sublaminados........................................................76
4.1 Elementos de la matriz S debido a tensiones en el plano.................................................. 77
4.2 Elementos de la matriz S debido a las tensiones fuera del plano...................................... 79
4.3 Elementos de S debido a tensiones cortantes fuera del plano........................................ 82
4.4 Conformación de la matriz de rigidez............................................................................ 83
BIBLIOGRAFÍA........................................................................................................................ 85 viii
CAPÍTULO 1
1
INTRODUCCIÓN
1.1 Motivación
En las últimas décadas la industria ha generado nuevos materiales, que debido a sus
características se adaptan mejor a las nuevas condiciones de trabajo a las que son
sometidos. Entre los materiales que han alcanzado un amplio desarrollo en los últimos
tiempos se encuentran los materiales compuestos, caracterizados como su nombre lo
indica por ser elaborados con dos o más componentes. Habitualmente uno de los
componentes cumple la función de reforzar al otro, para lograr una determinada
característica deseada como la resistencia mecánica. Entre estos materiales se puede
identificar, entre otros, a los compuestos reforzados por fibras.
Reedy [1] define como material compuesto a aquel construido con dos o más
componentes que juntos producen propiedades deseables que no pueden alcanzar
ninguno de ellos por sí solos. En el presente trabajo, se tratan materiales compuestos por
una matriz o base epóxi que recibe un refuerzo de fibras continuas, estas últimas
poseedoras de propiedades mecánicas que superan incluso a las de los metales más
resistentes. Los materiales de matriz epóxi con refuerzo de fibras, que en el presente trabajo
se denominan materiales compuestos o simplemente compuestos, han tenido una amplia
difusión en la industria automotriz de avanzada y particularmente en la aeronáutica.
La aplicación de los materiales compuestos en el campo aeronáutico ha ido creciendo
a pasos agigantados. Esto se debe esencialmente a la versatilidad de estos materiales,
principalmente para proveer características de resistencia óptimas para diversos tipos de
aplicaciones. Así mismo, el uso de estos materiales permite una reducción de peso
importante, que se logra reemplazando componentes estructurales metálicos de
aeronaves por otros de compuesto. Existen varios casos de aviones modernos que
explotan las características de los materiales compuestos, entre ellos un caso singular es
el Boeing 787 [2] donde un 50 % de su peso corresponde a materiales compuestos.
En la Figura 1-1 se muestran los materiales usados en el fuselaje del Boeing 787, su
ubicación y su peso dado como porcentaje del peso total. Se observa la participación
mayoritaria de los compuestos del orden del 50 %, seguido por aluminio ( 20 % ) ,
titanio (15 %) y acero (1 0 %). Para visualizar el progreso en el uso de materiales
compuestos basta mencionar que un modelo anterior, el Boeing 777 usa sólo un 12 %
de materiales compuestos y 50 % de aluminio.
1
Figura 1-1: Materiales usados en el fuselaje del avión Boeing 787, ( fuente www.boeing.com )
Para hacer un uso correcto de los materiales compuestos, es importante conocer
acabadamente las propiedades del material resultante. Esto es necesario a fin de
garantizar que la utilización del material cumplirá los requerimientos de diseño,
interacción y uso de los elementos finales. La determinación de propiedades puede
hacerse mediante mediciones experimentales, sin embargo en la mayoría de los casos
las mediciones resultan costosas y difíciles de realizar.
Más allá de la capacidad técnica requerida para implementar ensayos mecánicos para
determinar ciertas características y comportamientos, los elevados costos de estos procesos
pueden hacer inviable un proyecto. En contraposición las ecuaciones que gobiernan el
comportamiento de materiales compuestos pueden ayudar a estimar las características
del material. Sin embargo en la mayoría de los casos no son de fácil resolución a través
de métodos analíticos y eso lleva a aplicar métodos numéricos para encontrar resultados
aceptables y con bajo costo.
El método de los elementos finitos es una herramienta muy potente para el análisis
de estructuras. Mas aún, en la actualidad se cuenta con una gran cantidad de software
comercial con algoritmos computacionales robustos, capaces de simular de una manera
acertada el comportamiento de los materiales compuestos. Estos códigos han llegado a
un nivel de desarrollo tal que, además de presentar los resultados numéricos del
problema, brindan herramientas de reproducción gráfica en las que pueden ingresarse
las geometrías motivo de análisis y observarse su comportamiento a lo largo del estudio.
Dentro de este tipo de software se encuentra ABAQUS®, un programa basado en el
método de los elementos finitos que tiene la capacidad de resolver problemas que van
desde un simple análisis lineal hasta complejas simulaciones no lineales ( contacto entre
partes, interacción termo-mecánica, etc. ) [3].
2
En lo que respecta a la simulación computacional de elementos estructurales fabricados
en compuesto, a principios de la década pasada se inició el proyecto BOJCAS ( Bolted
Joints in Composite Aircraft Structures) [4]. En ese esfuerzo se reunieron: tres empresas
constructoras de aeronaves, cuatro laboratorios aeroespaciales, dos universidades y dos
compañías de investigación; de ocho países europeos.
El objetivo del proyecto BOJCAS fue analizar las juntas abulonadas de materiales
compuestos, muy utilizadas en la construcción de estructuras aeronáuticas modernas.
Entre otras configuraciones de unión, se analizaron placas de fibra de carbono unidas
mediante una solapa simple y con un solo bulón. Este esquema de unión fue sometido a
ensayos de tracción de manera experimental, y posteriormente estos ensayos fueron
replicados de manera numérica por McCarthy y otros [5].
En este trabajo de tesis se busca reproducir a nivel computacional, mediante el uso del
software comercial de análisis por elementos finitos ABAQUS®, una unión abulonada
de solapa simple de dos placas de fibra de carbono. Particularmente se pretende estudiar
el mismo esquema propuesto en el proyecto BOJCAS [4] y [5], debido a que los resultados
experimentales disponibles permiten validar la simulación computacional. Luego, una
vez replicados los resultados de BOJCAS, se busca extender el modelo numérico al estudio
de otras configuraciones de abulonado para estudiar el comportamiento de las mismas.
1.2
Justificación del proyecto
El desarrollo de simulaciones computacionales del comportamiento de uniones de partes
construidas con materiales compuestos, representa un importante aporte para ampliar el
campo de aplicación de estos materiales en la construcción de elementos estructurales en
todas las ramas de la ingeniería. Es de amplio conocimiento que los materiales compuestos
muestran un excelente comportamiento mecánico, con la excepción de las zonas donde
no existe continuidad de fibras como por ejemplo agujeros y bordes de unión.
Por ese motivo el desarrollo y validación de un modelo numérico, que sea capaz de
brindar resultados satisfactorios en la solución de ese tipo de problemas, resulta crucial
en el estudio de estructuras construidas con material compuesto.
Por otra parte la metodología numérica de análisis se presenta como una importante
herramienta de apoyo para complementar los estudios experimentales de uniones de piezas
construidas con compuestos. La simulación computacional abre un camino al desarrollo
de nuevos experimentos ( numéricos ) , disminuyendo los costos de análisis experimental
de piezas, esquemas de vinculación, tipo de unión, etc. Dicho esto, queda justificada la
conveniencia de explorar alternativas para simular por la vía de elementos finitos el
problema de las uniones abulonadas de piezas construidas con materiales compuestos.
3
1.3
Objetivos
El presente trabajo de tesis tiene como objetivo general, reproducir a nivel computacional los ensayos experimentales de uniones de placas laminadas en material compuesto
que fueron desarrollados dentro del proyecto BOJCAS. Para cumplir este objetivo se ha
empleado el software comercial de análisis por elementos finitos ABAQUS®.
Después de validar el modelo de unión simple, el siguiente objetivo es desarrollar y
estudiar otras configuraciones de abulonado a fin de analizar el comporta-miento de las
mismas ( deformaciones, tensiones y modos de falla).
Para alcanzar los objetivos generales enunciados en el párrafo anterior se plantean
los siguientes objetivos específicos:
•
Estudiar las características del programa de elementos finitos ABAQUS® con la
finalidad de familiarizarse con su uso.
•
Analizar el trabajo desarrollado por los participantes del proyecto BOJCAS [4], en
lo que refiere a la junta abulonada de solapa simple con un bulón.
•
Definir un modelo numérico apropiado para replicar el comportamiento del abulonado
de solapa simple de material compuesto.
•
Validar el modelo comparando los resultados numéricos del presente trabajo con
sus similares experimentales y numéricos del proyecto BOJCAS [4] y del trabajo de
McCarthy [5].
•
Sobre la base de los resultados obtenidos y usando las características del modelo
validado, desarrollar mediante el mismo software modelos para otras configuraciones
de abulonado con el fin de anticipar los resultados que se obtendrían al realizar un
ensayo de tracción sobre los mismos.
1.4 Alcance del proyecto
En este trabajo se ha considerado que los elementos constituyentes del material
compuesto tienen un comportamiento lineal, de este modo se puede inferir la zona de
falla del material pero no se ha incluido ningún modelo de daño. Por otro lado se ha
estudiado la interacción entre las placas de compuesto y también la interacción entre las
placas y el bulón. Para ello fue necesario considerar los contactos entre las diferentes
partes involucradas en la unión. El modelo empleado tiene en cuenta la fricción en las
zonas de contacto entre las partes, como así también el ajuste previo debido al torque
aplicado sobre el bulón.
4
1.5 Metodología de análisis
Con el fin de alcanzar los objetivos planteados para este trabajo, se sigue la siguiente
metodología:
1. Se analiza el proceso experimental llevado a cabo en el Proyecto BOJCAS [4] y los
trabajos realizados por McCarthy [5], para conocer en detalle: i) los criterios utilizados,
ii) los ensayos experimentales de tracción realizados y iii) las simulaciones numéricas
llevadas a cabo para estudiar el comportamiento de las uniones de placas de material
compuesto por medio del abulonado.
2. Se estudia la mecánica de los materiales compuestos, haciendo énfasis en las relaciones
tensión-deformación, las diferentes características de simetría de propiedades que se
puede lograr a partir de la combinación fibra-matriz epóxi y otros temas de interés
para el trabajo. Luego se identifica la característica que mejor se adapta al material,
experimentando con el fin de seleccionar un modelo matemático satisfactorio que lo
represente numéricamente.
3. Se utiliza el software ABAQUS® para el análisis numérico, validando un modelo que
resulta adecuado para simular los ensayos experimentales, con materiales, tipos de
elementos finitos y mallados que representan de forma adecuada las características
del experimento y de las componentes estructurales que motivan el análisis.
4. Se comparan los resultados numéricos, obtenidos en el presente trabajo para uniones
de solapa simple y con un solo bulón, con los datos experimentales reportados en el
Proyecto BOJCAS y los valores numéricos publicados en el trabajo de McCarthy
5. Se analizan otras configuraciones de abulonado con dos bulones para determinar el
comportamiento que se observaría al realizar un ensayo sobre las mismas.
6. Se utiliza el trabajo experimental desarrollado por Morales y Quadri [6] en la
Universidad Nacional de Córdoba para determinar coincidencias que permitan
reforzar la validez del modelo que busca predecir el comportamiento de otras
configuraciones diferentes a la analizada por McCarthy.
5
1.6 Contenido de la tesis
El contenido del presente trabajo se detalla a continuación. En el Capítulo 2, se
expone la mecánica de los materiales compuestos, haciendo hincapié en las principales
relaciones entre las tensiones, los desplazamientos y las deformaciones, y se formula la
matriz de rigidez de un material laminado.
En el Capítulo 3, se expone la mecánica de los materiales compuestos, abarcando las
características que mejor se ajustan a las necesidades del presente trabajo, así como la
teoría de los laminados. Eso facilita después conocer e identificar los tipos de comportamiento que se han presentado durante el análisis experimental, y permite plantear
posteriormente las ecuaciones matemáticas que los representan.
En el Capítulo 4 se desarrolla el modelo numérico computacional presentado en esta
tesis, el cual permite reproducir las condiciones bajo las cuales se llevaron a cabo las
experimentaciones del proyecto BOJCAS y los trabajos numéricos de McCarthy. En
este capítulo se describe el modelo desarrollado para replicar los resultados experimentales,
se valida el modelo, y se desarrollan modelos similares para otras configuraciones de
abulonado que permiten anticipar el comportamiento que se observaría al experimentar
sobre las mismas.
En el Capítulo 5 se exponen las conclusiones y recomendaciones derivadas del
presente trabajo. Se enuncian las contribuciones de esta tesis al ámbito académico y se
listan posibles líneas de trabajo futuras que se pueden derivar de la misma.
En el Anexo se detalla la metodología para obtener las matrices de rigidez de
laminados cuyas capas tienen diferentes orientaciones, así como también para el caso de
sub-laminados.
6
CAPÍTULO 2
2
MECÁNICA DE
MATERIALES COMPUESTOS
2.1 Introducción
Los compuestos laminados, de acuerdo a Paris y colab. [7], son considerados desde el
punto de vista macroscópico como la combinación de dos o más materiales con interfaces
de separación entre ellos. El material así obtenido muestra habitualmente propiedades
mecánicas que no pueden ser alcanzadas por sus constituyentes actuando en forma aislada.
Dado que los materiales compuestos presentan como ventaja elevadas relaciones
resistencia-peso y rigidez-peso, resultan adecuados para ser utilizados en estructuras en
las cuales el peso es una variable fundamental en el proceso de diseño. Sus aplicaciones,
principalmente en la industria aeronáutica, han sido notorias debido a esas características,
aunque también se observa un crecimiento de su campo de aplicación en otras ramas de
la industria como lo afirman Car y colab. [8]
En este capítulo se cubre la mecánica de los compuestos laminados. Se definen las
relaciones aplicables a un laminado sometido a cargas, en lo que refiere a los campos de
desplazamientos, deformaciones y tensiones; como así también las relaciones entre
ellos. Se han introducido además las relaciones asociadas al estado de tensión plana, que
es el más generalizado para el análisis de compuestos. Por último a partir de las
relaciones mencionadas, se obtienen las matrices de rigidez de las capas y la matriz de
rigidez global del laminado.
2.2 Compuestos laminados
Habitualmente los materiales compuestos de fibra continua están constituidos por
capas de material, pegadas entre si de modo que forman un laminado. Las capas que
conforman este tipo de compuestos, pueden estar constituidas de fibras largas, continuas
unidireccionales, tejidos multidireccionales ó de otro tipo; embebidas en una matriz
inicialmente líquida para facilitar la cohesión. Esta matriz luego de su curado a
temperatura ambiente, ó mediante tratamiento térmico, se endurece dando el soporte
necesario a las fibras del material laminado.
7
En la Figura 2-1 se muestran dos láminas de un compuesto laminado de fibra de
carbono en matriz epoxi donde se observa claramente la orientación de las fibras.
Figura 2-1: Compuesto laminado de fibra de carbono en matriz epoxi
Fuente http://www.fibretechsolutions.co.nz
En un compuesto las capas adyacentes, constituidas con el mismo material y con la
misma orientación, se denominan grupo de capas. Debido a que sus propiedades y sus
orientaciones son las mismas, pueden ser tratadas como una capa única.
En este trabajo se analiza un compuesto laminado, conformado por 40 capas de fibra de
carbono unidireccional con orientaciones diversas, preimpregnadas en una matriz epoxi.
2.2.1
Secuencia del laminado
Con el fin de analizar a los laminados, se utiliza un sistema de coordenadas ortogonal
x-y-z, en el que su eje z es perpendicular al plano del laminado. La orientación de las
capas unidireccionales es especificada por el ángulo θ con respecto al eje x contenido en
el plano local del laminado. El ángulo se considera positivo en el sentido antihorario.
El número de capas dentro de un laminado se especifica por medio de un subíndice,
que no figura cuando se trata de una capa única. En la Figura 2-2, a modo de ejemplo se
caracteriza la secuencia de laminado de un compuesto de 10 capas.
Figura 2-2: Descripción de la conformación de un laminado
8
Cuando existe un laminado simétrico con respecto al plano medio ( x - y ) , se indica
simetría mediante el subíndice “s”, en el exterior del corchete, describiendo únicamente
en el interior del mismo la configuración de la parte superior al plano medio.
Como se indicó previamente, en este trabajo se analiza un laminado construido con 40
capas de fibra de carbono unidireccional dispuestas en una configuración 45º,0º,−45º,90º,
secuencia que se repite hasta alcanzar la vigésima capa, luego de lo cual se invierte la
configuración, y se continúa hasta la cuadragésima capa. Utilizando el código descrito,
la secuencia del laminado estudiado en el presente trabajo es [45º/0º/−45º/90º] 5S . Esta
combinación de ángulos conforma una configuración denominada “cuasi-isótropa”,
debido a que la orientación y la secuencia de apilamiento es tal que el material resultante se
comporta como un material isótropo. Los trabajos de Beckwith [9] y Kóllar y Springer [10]
enuncian varios criterios para considerar a un material compuesto laminado como cuasiisótropo, los cuales se cumplen en el laminado analizado en el presente trabajo:
•
El ángulo de las fibras (θ en radianes ) u orientación de cada capa es θ = i π / n , en
donde i es un entero y n es el número total de orientaciones de la fibra.
•
Deben existir al menos tres direcciones diferentes de las fibras ( n ≥ 3 ) .
•
La combinación de capas debe ser balanceada, es decir, por cada ángulo positivo
+θ, sobre el plano medio del laminado, debe existir un ángulo −θ, bajo dicho plano.
•
A cada capa de la parte superior al plano medio del laminado, debe corresponderle
exactamente la misma capa en la parte inferior a dicho plano, las cuales deben tener
iguales los siguientes parámetros:

Material.

Fracción volumétrica de fibra y resina.

Espesor de la capa.

Tipo y geometría de la fibra.

Distancia al plano medio.
Bajo estas consideraciones, y dadas las condiciones del laminado que se analiza en el
presente estudio, se confirma la característica de cuasi isotropía del mismo.
2.3 Deformaciones y tensiones en compuestos
La versatilidad que se puede obtener al elaborar materiales compuestos de fibras en
matriz epóxi, proviene del hecho de que el material se comporta de manera diferente en
cada dirección según sea la orientación de las fibras.
9
Las deformaciones ε están definidas a partir del cambio en los desplazamientos u. De
esta manera teniendo en cuenta pequeñas deformaciones, las componentes del tensor de
deformaciones se pueden obtener a partir de la siguiente relación:
=
ε
εx 
ε y 
 
εz 
∂
γ  =
 yz 
γ xz 

γ xy 

 ∂ / ∂x
 0
 0

 0
 ∂ / ∂z
 ∂ / ∂y
0
∂ / ∂y
0
∂ / ∂z
0
∂ / ∂x
0
0
∂ / ∂z
∂ / ∂y
∂ / ∂x
0



 =
u



u
 
v

 w

→=
ε
∂u
(1)
Luego la relación constitutiva entre las tensiones y las deformaciones para un
material anisótropo y linealmente elástico, en un sistema coordenado x-y-z, puede
evaluarse a partir de la expresión:
σ x = C11 ε x + C12 ε y + C13 ε z + C14 γ yz + C15 γ xz + C16 γ xy
σ y = C21 ε x + C22 ε y + C23 ε z + C24 γ yz + C25 γ xz + C26 γ xy
σ z = C31 ε x + C32 ε y + C33 ε z + C34 γ yz + C35 γ xz + C36 γ xy
τ yz = C41 ε x + C42 ε y + C43 ε z + C44 γ yz + C45 γ xz + C46 γ xy
(2)
τ xz = C51 ε x + C52 ε y + C53 ε z + C54 γ yz + C55 γ xz + C56 γ xy
τ xy = C61 ε x + C62 ε y + C63 ε z + C64 γ yz + C65 γ xz + C66 γ xy
donde las Cij son constantes que definen al material.
La relación general (2) puede ser expresada matricialmente como:
σ x 
σ 
 y
σ z 
 =
τ yz 
τ xz 
 
τ xy 
 C11
C
 21
C31

C41
C
 51
C61
C12
C22
C32
C13
C23
C33
C14
C24
C34
C15
C25
C35
C42
C52
C43
C53
C44
C54
C45
C55
C62
C63
C64
C65
C16 
C26 
C36 

C46 
C56 

C66 
εx 
ε 
 y
 ε z 
 
γ yz 
γ xz 
 
γ xy 
→ σ =
Cε
(3)
donde la matriz C es la matriz de rigidez en el sistema coordenado x-y-z. Al invertir la
relación expresada previamente, obtenemos la siguiente:
εx 
ε 
 y
 ε z 
 =
γ yz 
γ xz 
 
γ xy 
 S11
S
 21
 S31

 S 41
 S51

 S61
S12
S13
S14
S15
S 22
S32
S 23
S33
S 24
S34
S 25
S35
S 42
S52
S62
S 43
S53
S63
S 44
S54
S64
S 45
S55
S65
10
S16 
S 26 
S36 

S 46 
S56 

S66 
σ x 
σ 
 y
σ z 
 
τ yz 
τ xz 
 
τ xy 
→ ε =
Sσ
(4)
En la ecuación (4), S es la matriz de flexibilidad del material, la cual resulta ser la
inversa de la matriz C. Para el caso de un material elástico se cumple que tanto la matriz
de flexibilidad como la de rigidez son simétricas con lo que únicamente 21 de los 36
elementos que las componen son independientes.
En el caso de láminas, se observa que las tensiones en dirección transversal al plano
de la lámina de material compuesto σ z , τ yz y τ xz se consideran nulas. Esto conlleva a un
estado de tensión plana, donde (4) se transforma en:
εx 
 
ε y 
γ 
 xy 
 S11 S12 S16  σ x 
εz 
 
 
=

γ yz 
 S12 S 22 S 26  σ y 


γ 
 S16 S 26 S36  τ xy 
 xz 
 S13
S
 14
 S15
S 23
S 24
S 25
S36 
S 46 
S56 
σ x 
 
σ y 
τ 
 xy 
(5)
A partir de las relaciones planteadas, y bajo consideraciones atribuibles a los arreglos
de las fibras, los compuestos pueden caracterizarse en general como anisótropos cuando
no hay planos de simetría respecto al alineamiento de las fibras. Los compuestos
monoclínicos son aquellos donde hay un plano de simetría respecto al alineamiento de
las fibras. El caso más habitual es el de los compuestos ortótropos donde hay tres planos
de simetría mutuamente perpendiculares respecto al alineamiento de las fibras. Los
transversalmente isótropos tienen las mismas características del ortótropo pero uno de
sus planos de simetría puede tratarse como isótropo. Finalmente para el caso de los
metales se pueden tratar como materiales isótropos dado que todos los planos son
planos de simetría.
Las diferentes características que presentan los compuestos hacen que se modifique
la matriz de flexibilidad para adaptarse a esas particularidades. Para el caso de un
material ortótropo con planos de simetría x-y, x-z y y-z, donde σ x que está en el plano
γ=
0 ),
de simetría x-y que no provoca deformaciones cortantes fuera del plano ( γ=
xz
yz
la tensión de corte τ xy también es nula debido a que σ x también está en el plano de
simetría x-z con lo que los elementos S 14 , S 15 y S 16 son nulos. Con argumentos
similares, se puede demostrar que para un material ortótropo los elementos S 24 , S 25 , S 26,
S 34 , S 35 , S 36, S 45 , S 46 , y S 56 también son nulos. De este modo la relación deformacióntensión (4) queda establecida por:




S = 




0 
S 22 S 23 0 0 0 
S33 0 0 0 

S 44
0 0 
sim
S55
0 

S66 
S11
S12
S13
11
0 0 (6)
Consecuentemente la matriz de rigidez se reduce a:
0 0  C11 C12 C13

C22 C23 0 0 

C33 0 0 C= 
 C44 0  sim
C55

 0 

0 
0 

0 
0 

C66 
(7)
En las matrices (6) y (7) que contienen 12 elementos diferentes de cero, sólo 9
elementos son independientes entre sí. Luego la matriz C puede representarse de la
siguiente forma:




C = 
 0
 0

 0
0
0
0
L
0
0
0
0
0
0
0
0
0
M
0
0
0









(8)
donde las submatrices L y M, pueden escribirse como:
L
=
  E3 2 
 E1 1 − ν 23 

  E2
1 

D 


sim



E
E2 ν 12 + 3 ν 13 ν 23 
E2


 E

E2 1 − 3 ν 132 
 E1 




 
E2
E3 ν 23 + ν 13 ν 12  
E1

 

 E


E3 1 − 2 ν 122 
E



1
E3 (ν 13 +ν 12 ν 23 )
0 
 G23 0


G13 0 
M = 
 sim
G12 

(9)
(10)
siendo el denominador D igual a:
D=
2
E1 E2 E3 −ν 23
E1 E32 −ν 122 E22 E3 −ν 13 ( 2ν 12ν 23 −ν 13 ) E2 E32
E1 E2 E3
(11)
Para el caso del laminado HTA/6376 que se analiza en esta tesis, cuyas propiedades
mecánicas se listan en la Tabla 2-1 la matriz de flexibilidad para la capa orientada a 0º
está dada por la ecuación (12):
12
Tabla 2-1: Propiedades de rigidez unidireccional [GPa] para una lámina de carbono/epóxi HTA/6376
E 11
E 22
E 33
G 12
G 13
G 23
ν 12
ν 13
ν 23
140
10
10
5,2
5,2
3,9
0,3
0,3
0,5
0
0
0
0
0
0
0,3
 1
−
 140 140

1


10



S = 



sim





0,3
140
0,5
−
10
1
10
−
1
3,9
0
1
5, 2

0 

0 


0 
2

−9 m
x
1
0

N
0 


0 

1 
5, 2 
(12)
Invirtiendo la matriz (12) se obtiene la matriz de rigidez para una capa del material
compuesto estudiado, la cual es:
 143, 7



C = 





6,16
6,16
0
0
13, 60
6,93
0
0
13, 6
0
0
3,9
0
sim
5, 2
0 

0 
0 
m2
 x1
0 9
N
0 

0 

5, 2 
(13)
2.4 Matriz de rigidez de un laminado delgado
En este trabajo se aplican los conceptos y desarrollos de Kóllar y Springer [10]. El
análisis que se presenta está basado en la teoría de placas laminadas asumiendo que los
laminados son planos. Por conveniencia, se ubica al laminado en un sistema de
coordenadas referencial x-y-z, y se usa un plano de referencia en la mitad del laminado
como se muestra en la Figura 2-3.
Se supone que el laminado experimenta pequeñas deformaciones y que las deformaciones varían linealmente a través del laminado en la dirección z. De este modo las
deformaciones cortantes son despreciables, dado que la tensión fuera del plano ( normal
σ z y las cortantes τ xz , τyz ) son pequeñas en comparación con las tensiones en el plano:
σ x , σ y y τxy. Estas aproximaciones implican que las relaciones tensión-deformación
pueden aplicarse bajo consideraciones de tensión plana.
13
Figura 2-3: Plano de referencia y ejes coordenados
Las deformaciones membranales se dan en la ecuación (14) donde u y v representan las
componentes del desplazamiento en x y y respectivamente y el superíndice o variables
asociadas al plano medio como se indica en la Figura 2-4.
=
ε xo
∂u o
∂x
=
ε yo
∂v o
∂y
=
γ xyo
1  ∂u o ∂v o 
+


∂x 
2  ∂y
(14)
Adoptando la hipótesis de Kirchhoff, que indica que las normales a la superficie de
referencia permanecen normales y rectas, y suponiendo pequeñas deflexiones, el ángulo
de rotación de la normal al plano de referencia es:
=
β xz
ángulo de rotación
=
∂ wo
∂x
donde w o es el desplazamiento fuera del plano de referencia.
Figura 2-4: Deformación de una placa en el plano x-z
14
(15)
El desplazamiento total en la dirección x e y está dado por:
∂ wo
u =u − z β x z =u − z
∂x
∂wo
v=
vo − z β y z = vo − z
∂y
o
o
(16)
Sustituyendo (16) en las definiciones de deformación (14) se obtiene:
∂u
=
∂x
∂v
=
∂y
=
εx
=
εy
∂ uo
∂ 2 wo
−z
∂x
∂ x2
∂ vo
∂ 2 wo
−z
∂y
∂ y2
1∂u ∂v
+

=
2 ∂y ∂x
γ x y=
(17)
∂ 2 wo
1  ∂ u o ∂ vo 
z
+
−


∂x 
∂x∂ y
2 ∂y
Las derivadas segundas en la ecuación (17) miden las curvaturas del plano de referencia,
∂ 2 wo
∂ 2 wo
κ
=
y
∂ x2
∂ y2
de modo que (17) se puede escribir como:
κx =
εx 
 
εy 
=
γ 
 xy 
κxy =
∂ 2 wo
∂x∂ y
(18)
 ε xo 
κx 
 o
 
ε y  + z κ y 
γ o 
 
κ xy 
 xy 
(19)
Las fuerzas y momentos actuando sobre un elemento, se indican en la Figura 2-5 y se
obtienen integrando las tensiones membranales en el espesor según (20):
hs
hs
=
Nx
σ dz
M
z σ dz
V ∫
∫=
∫=
=
Ny
σ y dz
My
∫=
=
N xy
τ dz
M
∫=
∫
−h i
x
x
hs
−h i
hs
− hi
xy
xy
x
−h i
hs
x
z σ y dz
Vy
∫=
− hi
hs
−h i
hs
− hi
∫
hs
− hi
τ xz dz
τ yz dz
z τ xy dz
Figura 2-5: Fuerzas y momentos actuando en el plano de referencia
15
(20)
En la ecuación (20) los esfuerzos N i representan las fuerzas en el plano distribuidas
sobre el contorno de la placa, M i son los momentos distribuidos en el contorno, y V i son
las fuerzas cortantes transversales distribuidas en el contorno. Los límites de integración
h s y h i , son las distancias desde el plano de referencia a la superficie superior e inferior
del laminado, respectivamente. Recordando la relación tensión-deformación para tensión
plana (5), y denominando Q a la matriz de rigidez, se tiene:
σ x 
 εx 
 
 
σ y  = Q  ε y 
τ 
γ 
 xy 
 xy 
(21)
En la ecuación (21), Q es la matriz de rigidez de la capa en el sistema coordenado
x-y, y la misma está relacionada con la matriz de rigidez Q en un sistema de
coordenadas x 1 -x 2 -x 3 ( donde el eje x 1 se encuentra alineado con la dirección
longitudinal de la fibra ). Mediante la transformación (22), cuya deducción está en el
Anexo de esta tesis, se puede establecer que:
 Q11 Q12 Q16 


Q22 Q26 
=
Q =
 sim
Q66 

Tσ −1
 Q11 Q12

Q22

 sim

Q16 

Q26  Tε
Q66 
(22)
Las matrices T σ y T ε son las matrices de transformación, relacionadas con las
funciones trigonométricas del ángulo de orientación de las fibras, θ, dadas por:

cos 2 θ

Tσ 
=
sen 2 θ
 − cos θ senθ


cos 2 θ

=
Tε 
sen 2 θ
 −2 cos θ sen θ

sen 2 θ
2 cos θ sen θ 

−2 cos θ sen θ 
cos 2 θ − sen 2 θ 
cos θ
2
cos θ senθ
sen 2 θ
(23)
cos θ sen θ


− cos θ sen θ 
cos 2 θ − sen 2 θ 
cos θ
2
2 cos θ sen θ
(24)
Reemplazando (21) en las ecuaciones de fuerzas y momentos (19) se obtiene:
 ε xo 
 
Q d z  ε yo  +
ε o 
 xy 
 Nx 


=
 Ny 
N 
 xy 
∫
 Mx 


=
My 
M 
 xy 
ε

∫− hi z Q d z  ε
ε

hs
− hi
 κx 
 
∫− hi z Q d z  κ y 
κ 
 xy 
hs
(25)
hs
o
x
o
y
o
xy


+


16
 κx 
 
2
∫− hi z Q d z  κ y 
κ 
 xy 
hs
Partiendo de (25) se pueden evaluar las matrices de rigidez del laminado, que se
definen entonces como:
=
A
hs
hs
−h i
−h i
Q dz
B ∫
∫=
∫
=
z Q dz
D
hs
−h i
z 2 Q dz
(26)
siendo los elementos de las submatrices (26) :
=
Ai j
hs
Q dz
B ∫
∫=
ij
−h i
ij
hs
=
z Qi j d z
Di j
−h i
hs
z Q dz
i
∫=
2
ij
−h i
1, 2, 6
(27)
Debido a que Q es constante en cada capa, las integrales de las ecuaciones que
anteceden pueden reemplazarse por sumatorias:
Ai j=
1 m
1 m
2
2
Q
z
−
z
D
=
Qi j ) ( zk3 − zk3−1 ) (28)
(
)
(
(
)
∑
∑
ij k
k
ij k
k
k −1
ij
k
2 k 1=
3k 1
1 =
∑ (Q ) ( z
m
k
− zk −1 )
Bi j=
en donde m es el número total de capas o grupos de capas en el laminado, z k y z k-1 son las
distancias desde el plano de referencia a las dos superficies de la k-ésima capa y ( Qi j ) k
son los elementos de la matriz de rigidez de la k-ésima capa.
Con las definiciones (28), las expresiones para las fuerzas y momentos en el plano se
convierten en:
 Nx 
 A11
N 
 A
 y 
 12
 N xy 
 A16
=



 B11
 Mx
 B
M 
 12
 y
 B16
 M xy 
A12
A22
A16
A26
B11
B12
B12
B22
A26
A66
B16
B26
B12
B22
B26
B16
B26
B66
D11
D12
D16
D12
D22
D26
o
B16   ε x
 o
B26   ε y
o
B66   γ xy
 
D16   κ x
D26   κ y

D66  κ xy










(29)
En la ecuación (29) los vectores a izquierda y derecha representan los esfuerzos y
deformaciones generalizadas. Invirtiendo la ecuación (29), se obtienen las deformaciones
y curvaturas en términos de las fuerzas y momentos en el plano:
 ε xo
 o
 εy
o
 γ xy


 κx
κ
 y
 κ xy

 α11

 α

 12

 α16


=

 β11



 β12

 β16


α12
α 22
α 26
α16
α 26
α 66
β11
β12
β16
β12
β 22
β 26
β16
β 26
β 66
β12
β 22
β 26
β16
β 26
β 66
δ11
δ12
δ16
δ12
δ 22
δ 26
δ16
δ 26
δ 66
de modo que resulta la relación:
17



















Nx 
N y 
N xy 

Mx 
My 

M xy 
(30)
 α11 α12 α16
α
 12 α 22 α 26
 α16 α 26 α 66

 β11 β12 β16
β
 12 β 22 β 26
 β16 β 26 β 66
β11 β12
β12 β 22
β16 β 26
β16 
β 26 
β 66 
δ11 δ12 δ16
δ12 δ 22 δ 26
δ16 δ 26 δ 66




 = 








A11
A12
A16
B11
B12
A12
A22
A26
B12
B22
A16
A26
A66
B16
B26
B11
B12
B16
D11
D12
B12
B22
B26
D12
D22
B26
B66
D16
D26
B16
B16 
B26 
B66 

D16 
D26 
D66 
−1
(31)
Las submatrices A, B y D representan la rigidez del laminado, y permiten obtener la
respuesta del mismo a fuerzas y momentos en el plano medio o de referencia:
• A ij es la rigidez en el plano que relaciona las fuerzas N con las deformaciones en
el plano ε.
• D ij es la rigidez a la flexión que relaciona los momentos M con las curvaturas κ.
• B ij es la rigidez de acople en el plano y fuera del mismo, que relaciona las fuerzas N
con las curvaturas κ y los momentos M con las deformaciones en el plano ε.
La existencia de la matriz B implica el acoplamiento entre la flexión y la elongación,
debido a que tanto fuerzas y curvaturas, momentos y deformaciones existen de manera
simultánea. En un laminado simétrico como el caso del analizado en este trabajo, la
capa ubicada en la posición +z es idéntica a la ubicada en −z, por lo que su matriz de
rigidez es la misma, es decir:
Q( z ) = Q ( − z )
(32)
Sustituyendo (32) en (26) la matriz B se anula, por lo que para un laminado simétrico
no hay acoplamiento flexión-elongación y la ecuación (29) se reduce a:
 Nx

 Ny
N
 xy


=


 A11


 si m

 Mx 


My  =


 M xy 
 D11


 si m

A12
A22
D12
D22
A16 

A26 
A66 
 ε xo
 o
εy
γ o
 xy





(33)
D16 

D26 
D66 
 κx 
 
κy 
 
κ xy 
(34)
Para este caso, la matriz de flexibilidad se puede expresar como:
 A11


 sim
A12
A22
A16 
A26 
A66 
−1
 a11

= 
 sim
18
a12
a22
a16 
a26 
a66 
(35)
 D11


 sim

D12
D22
D16 

D26 
D66 
−1
 d11

= 
 sim

d16 

d 26 
d 66 
d12
d 22
(36)
Con lo que la relación entre las deformaciones y curvaturas con las fuerzas y
momentos se puede reescribir como:
 ε xo 
 o
εy  =
γ o 
 xy 
 a11


s i m

 κx 
 
κy  =
κ 
 xy 
 d11


 si m

a12
a22
d12
d 22
a16 

a26 
a66 
 Nx 


 Ny 
N 
 xy 
(37)
d16 

d 26 
d 66 
 Mx 


My 
M 
 xy 
(38)
Para el caso de capas ortótropas los elementos Q16 y Q 26 en la matriz Q son nulos,
con lo cual los valores de los elementos A16 , A26 , D16 , D26 de las matrices A, B y D
también se anulan, y esto resulta en:
 Nx 
 
 Ny  =
N 
 xy 
 A11

 A12
0

 Mx 


My  =


 M xy 
 D11

 D12
 0

A12
A22
0
D12
D22
0
0 

0 
A66 
 ε xo 
 o
εy 
γ o 
 xy 
(39)
0 

0 
D66 
 κx 
 
 κy 
κ 
 xy 
(40)
y por lo tanto la relación entre las deformaciones y curvaturas con las fuerzas y
momentos termina siendo:
 ε xo   a11
 o 
 ε y  =  a12
γ o   0
 xy  
 κ x   d11
  
 κ y  =  d12
κ   0
 xy  
a12
a22
0
d12
d 22
0
19
0 

0 
a66 
 Nx

 Ny
N
 xy





(41)
0 

0 
d 66 
 Mx 


My 
M 
 xy 
(42)
2.4.1
Matriz de rigidez del laminado usado en la tesis
Para este análisis, bajo las consideraciones de tensión plana, la matriz de flexibilidad S
de una capa orientada a 0º fabricada con el material HTA/6376 es:
0,3
 1
−
 140 140

1
S = 

10

 sim


0 

m2
0  0− 9
x1

N

1 
5, 2 
(43)
luego la matriz de rigidez Q para esta capa con orientación 0º, resulta de invertir la
matriz de flexibilidad (44), siendo:
140,91

Q = Q0 = 

 sim
3, 02
10, 07
0 

N
0  x 1 0 9 2
m

5, 20 
(44)
Aplicando sobre la matriz (44), la transformación de coordenadas prevista en (22),
las matrices de rigidez para las capas orientadas a 45º, −45º y 90º resultan:
 44, 45

Q 45 = 

 sim
44, 45
 44, 45



 sim
34, 05
 10, 07

= 

 sim
3, 02
Q −45
=
Q90
34, 05
44, 45
140,91
32, 71

N
32, 71 x 1 0 9 2
m

36, 23
(45)
−32, 71

N
−32, 71 x 1 0 9 2
m

36, 23 
(46)
0 

N
0  x 1 09 2
m

5, 20 
(47)
Basado en (28) y considerando que la matriz B es nula, las matrices A y D para un
laminado HTA/6376 con una configuración [45º/0º/−45º/90] 5s, resultan:
20
 311,84

A = 

 sim
96,39
311,84
 743, 65

D = 

 sim


9 N
0  10
x
m

107, 73 
(48)
30,18 

30,18  x 10 9 Nm

256,37 
(49)
0
230,82
634, 41
La matriz de fuerzas y momentos en el plano para el compuesto que se estudia en el
presente trabajo, según lo expuesto en (29), es:
 N x  311.84

 
 Ny  

 
 N xy  

 

=
 Mx  

 
My  

 
 M xy  
96,39
0
0
311,84
0
0
107, 73
0
743, 65
sim
21
 ε xo



 ε yo
0
0 


 γ xyo
0
0 
9 
 x 10 
230,82 30,18 
 κx


634, 41 30,18 
 κy


 κ xy
256,37 
0
0












(50)
22
CAPÍTULO 3
3
MODELACIÓN NUMÉRICA DE COMPUESTOS
3.1
Introducción
El método de elementos finitos es una técnica poderosa para la solución de
ecuaciones diferenciales e integrales que gobiernan una gran variedad de problemas de
ingeniería y ciencias aplicadas. Actualmente existe una gran cantidad de software
comercial con capacidad de resolución de diferentes problemas de ingeniería a través
del método de elementos finitos.
En este capítulo se describe el estado del arte del análisis de materiales compuestos
aplicando métodos numéricos. Se exponen las teorías que han propuesto algunos
autores para analizar este tipo de materiales, en particular la de Reedy [1]. También se
describe el procedimiento general para tratar numéricamente las ecuaciones de
gobierno, la ecuación constitutiva adoptada, y el tipo de elemento finito utilizado. Se
define la aproximación empleada, sus funciones de forma y el desarrollo de las matrices
típicas que permiten el análisis por elementos finitos.
3.2
Estado actual del análisis numérico de materiales compuestos
Zienkiewicz y colab. [11] afirman que las técnicas numéricas, fundamentalmente el
método de elementos finitos, facilitan el análisis de estructuras elaboradas con materiales
compuestos. Car y colab. [7] clasifican a las teorías más utilizadas en el ámbito del
método de elementos finitos para el análisis de materiales compuestos, de la siguiente
manera:
1. Teoría de capa única equivalente.
2. Teoría de elemento sólido bi y tridimensional.
3. Teoría de aproximación bidimensional por capas.
Entre los primeros trabajos que se ocuparon de desarrollar un análisis tridimensional de
compuestos laminados por el método de los elementos finitos está el de Barker y colab.
[12]. Estos investigadores llegaron a la conclusión de que es necesario contar con
elementos individuales en cada capa del compuesto cuando se tienen diferentes
orientaciones de las fibras.
En la Figura 3-1 se bosquejan los dos tipos de análisis por elementos finitos de placas
laminadas: i ) con elementos de placa ó ii ) modelando con elementos tridimensionales.
23
Figura 3-1: Placa de material compuesto
( a ) Análisis como placa. ( b ) Análisis con elementos tridimensionales
En el estado actual del desarrollo de los materiales compuestos, es práctica habitual
utilizar una gran cantidad de capas, pero de un espesor relativamente reducido. Bajo las
consideraciones expresadas por Barker y colab. [12], una formulación de elementos
finitos que cumpla con esas consideraciones requeriría un número considerable de
elementos tridimensionales, lo cual conllevaría un costo computacional elevado en el
análisis de estructuras de grandes dimensiones.
Con el paso de los años se ha desarrollado el concepto de sub-laminado para subsanar
este problema. De esa manera se reduce el costo computacional y se pueden usar elementos
tridimensionales en el caso de grandes estructuras de compuesto. Según exponen Kóllar
y Springer [10], un sub-laminado permite modelar varias capas de material compuesto
como un elemento finito único a través del espesor. En un sub-laminado las propiedades
del material se obtienen mediante integración de las propiedades de las láminas a través
del espesor del sub-laminado. En la Figura 3-2 se esquematiza la discretización en sublaminados y la malla de elementos finitos.
Figura 3-2: Laminado, sub-laminado y malla de elementos finitos
La matriz de rigidez C del sub-laminado define la relación entre la tensión promedio
σ y la deformación promedio ε :
24
σ x 
σ 
 y
σ z 
=
τ 
 yz 
τ xz 
 
τ xy 
εx 
ε 
 y
εz 
C 
γ yz 
γ xz 
 
γ xy 
=
→ σ
Cε
(1)
donde los términos que contienen una barra son las tensiones y deformaciones
promedio. Denotando los elementos de la matriz de flexibilidad S se tiene:
 ε x   S1 1
ε  
 y   S2 1

 εz 
  S3 1
 =
γ y z   S 4 1
γ x z   S5 1
  

γ x y 
  S6 1
S1 2
S1 3
S1 4
S1 5
S2 2
S2 3
S2 4
S2 5
S3 2
S3 3
S3 4
S3 5
S4 2
S4 3
S4 4
S4 5
S5 2
S5 3
S5 4
S5 5
S6 2
S6 3
S6 4
S6 5
S1 6 
S 2 6 
S3 6 

S4 6 
S5 6 

S6 6 
σ x 
σ 
 y

σ z 

 
τ y z 
τ x z 
 

τ x y 

(2)
Las tensiones promedio en (1) y (2) están definidas por las siguientes ecuaciones:
=
σx
=
σy
=
τ xy
1
1
dz
Nx =
=
σ x σz σz
∫
hs hs
hs
1
1
dz
Ny =
=
σ y τ xz τ xz
∫
h
hs s
hs
1
1
dz
N x y=
=
τ x y τ yz τ yz
∫
h
s
hs
hs
(3)
donde h s es el espesor del sub-laminado, ( N x , N y , N xy ) son los esfuerzos membranales
por unidad de longitud y (σ x , σ y , τ xy ) son las tensiones membranales.
Las deformaciones promedio son:
1
=
ε z εx εx
dz
hs ∫hs
1
dz
=
γ yz =
γ y z εy εy
hs ∫hs
1
dz
γ xz =
γ x z γ xy = γ xy
hs ∫hs
=
εz
(4)
En las ecuaciones (3) y (4) se observa que las tensiones σ z , τ xz , τ yz y las deformaciones
ε x , ε y , γ xy no varían a través del espesor.
Kóllar y Springer [10] indican que los elementos de la matriz S se obtienen
asumiendo que el sub-laminado se encuentra sometido a las siguientes condiciones:
•
Tensiones en el plano ( σ x , σ y , τ xz ), que permiten obtener los elementos de la
primera, segunda y sexta columna de la matriz.
25
•
Tensión normal al plano del laminado (σ z ) que permite obtener los elementos
de la tercera columna de la matriz.
•
Tensiones cortantes fuera del plano ( τ yz , τ xz ), con las que se determinan los
elementos de la cuarta y quinta columna de la matriz.
A partir de la obtención de la matriz de flexibilidad S 6 x6 , es posible obtener la matriz
de rigidez C del sub-laminado. Esta última matriz permite considerar al laminado como
un sólido tridimensional. El desarrollo de la formulación que permite obtener los
elementos de la matriz se detalla en el Anexo de esta tesis, y su resultado final es:



S= 




S11
S 21
S31
0
0
S61
S12
S 22
S32
0
0
S62
S13
S 23
S33
0
0
S63
0
0
0
S 44
S54
0
0
0
0
S 45
S55
0
S16
S 26
S36
0
0
S66








(5)
donde:
 S11
 S 21
S
 61
S12
S 22
S62
S16 
S 26  = hs A
S66 
=
[ S31 S32 S36 ]
∑ (  S
M
k =1
 S13 
 S11
 S 23  = −  S 21
S 
S
 63 
 61
(5-a)
S12
S 22
S62
31
S32
)
S36  k ( zk − zk −1 ) Q k A −1
 C13 

S16 
M
z
z
−
1




1
k
k
−
S 26  ∑
 C23  C


h
1
k
=
s
 C63  ( 33 )k 
S66 
k


 C 

1 M  zk − zk −1 
1 M   13  zk − zk −1 
=
[ S33 ]
∑  ( C )  − [ S31 S32 S36 ] h ∑
 C23  C

hs k 1 =
(
=
k 1
s
33
33 ) k 

k 

C
  63 k

 S 44 S 45 
=
 S54 S55 
S
1 M 
∑  ( zk − zk −1 )  44
hs k =1 
 S54
S 45  

S55  k 
(5-b)
(5-c)
(5-d)
(5-e)
En la ecuación (5) los términos S ij no listados son nulos. Como ya se indicó h s es el
espesor del sub-laminado, y A 3x3 es la matriz de rigidez que relaciona las fuerzas con
las deformaciones generalizadas, definida en las ecuaciones (25) y (26) del Capítulo 2 .
z k es la distancia de cada capa al plano de referencia. S y C son las respectivas
matrices de flexibilidad y rigidez de la capa única de fibra unidireccional. Estas matrices
deben ser calculadas para cada capa del laminado, según la orientación que tienen las
fibras dentro de la placa. M es el número total de capas que conforman la placa.
26
Para el caso de este trabajo, la placa de material compuesto de fibra de carbono en matriz
epóxi HTA/6376 está laminada en la disposición [45º/ 0º/ –45º/ 90] 5s. A este material se
le ha aplicado el concepto expuesto en los párrafos precedentes, dividiéndolo en 10 sublaminados de modo que cada uno de ellos tiene una disposición [45º/0º/–45º/90]. Esta
subdivisión permite tener diez elementos (en lugar de 40) en el espesor de la placa en el
análisis por elementos finitos. Aplicando el concepto de sub-laminado y las ecuaciones (5),
se obtiene la matriz de flexibilidad S y la matriz de rigidez C del compuesto que se utiliza
en esta tesis.
S
0
0
0 
18, 44 −5, 70 −6,13
 −5, 70 18, 44 −6,13
0
0
0 


0
0
0 
 −6,13 −6,13 79, 43
×10−12
 0
0
0
224,36
0
0 


0
0
0
224,36
0 
 0
0
0
0
0
48, 27 
 0
(6)
C
 63,11
 21, 68

 6,54
 0

 0
 0
(7)
21, 68
63,11
6,54
0
0
0
6,54
6,54
13, 60
0
0
0
0
0
0
4, 46
0
0
0
0
0
0
4, 46
0
0 
0 

0 
×109
0 

0 
20, 72 
Las propiedades ortótropas equivalentes del sublaminado se pueden obtener a partir
de las relaciones de la ley de Hooke aplicadas a la matriz de flexibilidad. Estas
propiedades se usan como una opción para modelar las placas de compuesto en el
análisis como elemento sólido continuo. La Tabla 3-1 presenta las constantes
ingenieriles (elementos de la matriz C ) para el laminado HTA/6376.
Tabla 3-1: Propiedades de rigidez ortótropa equivalente para el material HTA/6376 en GPa obtenidas a partir
de la matriz C
E xx
E yy
E 33
G xy
G xz
G yz
ν xy
ν xz
ν yz
54,23
54,23
12,59
20,72
4,46
4,46
0,309
0,332
0,332
3.3
Análisis por elementos finitos
El análisis por elementos finitos como lo describe Kóllar y Springer [10] comprende
los siguientes pasos:
a. Generación de una malla que abarque la estructura, con lo que se obtiene un conjunto
de subdominios simples que no se intersecan entre sí, y cuyo comportamiento se
especifica mediante parámetros asociados a puntos (nodos ).
27
e
b. Determinación de la matriz de rigidez K de cada elemento,
K e δe = fe
(8)
siendo f e el vector de fuerzas actuando en los nodos del elemento y δ
e
los
desplazamientos incógnitas nodales del elemento.
La matriz de rigidez del elemento puede obtenerse a partir de principios de formulación variacional, tales como el principio de la mínima energía potencial, principio de la
mínima energía complementaria o el principio de la energía estacionaria de Reissner.
Además pueden aplicarse métodos de residuos ponderados como el de Galerkin. La
conformación de la matriz de rigidez global K de la estructura se obtiene por medio
e
de ensamblaje de las matrices de rigidez de los elementos K .
c. Las cargas aplicadas en la estructura son reemplazadas por un sistema de fuerzas
equivalentes f que actúan sobre los nodos y se obtiene.
Kδ= f
(9)
donde δ es el vector de desplazamientos nodales.
d. Los desplazamientos nodales son calculados a partir de la solución del sistema (9), lo
que es equivalente a:
δ = K −1 f
(10)
e. El vector δ es subdividido en subvectores δe , cada uno de los cuales representa los
desplazamientos de los puntos nodales de un elemento particular.
f. Los desplazamientos en un punto dentro de un elemento cualquiera se calculan a
partir de la expresión:
u = N δe
(11)
siendo u el vector de los desplazamientos y N la matriz de funciones de forma.
g. Las deformaciones en un punto dentro de un elemento se calculan por :
ε=Βδ
(12)
donde B es la matriz deformación-desplazamiento (a veces denominada matriz
cinemática dado que relaciona deformaciones con desplazamientos ), y se define más
adelante.
h. Las tensiones en un punto dentro de los elementos se calculan por medio de la
relación:
(13)
σ =Εε
siendo E ≡ C la matriz de rigidez característica del material.
28
3.4
Ecuación de gobierno
El principio de los trabajos virtuales provee una alternativa a las ecuaciones de equilibrio
para la formulación de problemas por elementos finitos, entre los autores que aplican
esta metodología se puede nombrar a Barbero [13], Liu y Quek [14] y Oñate [15]. El uso
del principio de trabajos virtuales es la metodología empleada por el software que se ha
aplicado en el presente trabajo, tal como lo indica el Manual de Teoría de ABAQUS® [3].
El principio de los trabajos virtuales permite obtener la forma débil de las ecuaciones
de movimiento del sólido, y en condiciones estacionarias (despreciando los efectos de
aceleraciones) permite recuperar las ecuaciones de equilibrio. En la formulación típica
de elementos finitos el principio de trabajos virtuales resulta en:
∫
V
σ i j δ ε i j 0
d V − ∫ ti δ ui d S − ∫ fi δ ui d V =
S
V
(14)
Siendo t i las fuerzas de tracción por unidad de área actuando sobre la superficie, f i
son las fuerzas de cuerpo por unidad de volumen, δε ij son las deformaciones virtuales y
δ u i los desplazamientos virtuales. El primer término de la ecuación (14), representa el
trabajo virtual interno, desarrollado por las tensiones y los otros términos el trabajo
virtual externo, desarrollado por las fuerzas externas.
Expandiendo la ecuación en (14) al estado tridimensional de deformaciones, se obtiene
que el trabajo virtual interno está dado por:
δ W=
int
∫ (σ
V
xx
δ ε x x + σ y y δ ε y y + σ z z δ ε z z + τ y z δ γ y z + τ x z δ γ x z + τ x y δ γ x y ) d V
(15)
que en notación compacta se escribe como
δ Wi n t =
∫
V
σ T δ ε d V
(16)
mientras que el trabajo virtual externo está dado por:
=
δ We x t
∫
S
t T δ u dS + ∫ f T δu dV
V
(17)
Teniendo presente las expresiones del Capítulo 2 para las deformaciones, ecuación
(1), y para las tensiones, ecuación (3), se tiene:
ε = ∂u
(18)
σ = Cε
(19)
Partiendo de (18), y dado que los desplazamientos virtuales son desplazamientos
compatibles con los vínculos y con la continuidad estructural, las deformaciones virtuales
están dadas por:
δε = ∂ δu
29
(20)
por tanto, llevando (19) y (20) a (16), el equilibrio entre el trabajo virtual externo (17) y
el trabajo virtual interno (16) puede escribirse como:
∫
V
∫
εT C
=
∂ δu dV
S
t T δu dS + ∫ f T δu dV
(21)
V
Las integrales sobre el volumen y la superficie del cuerpo, pueden ser representadas
como una sumatoria extendida a los m elementos en que se divide al cuerpo, de modo que:
m
m
 =
 t T δu dS + f T δu dV 
ε T C ∂ δu dV 
∑
∑
∫
∫V


∫S

V


e 1=
e 1 
=
(22)
Los desplazamientos dentro del elemento finito se aproximan como una combinación
lineal de las funciones de forma N y los desplazamientos nodales d:
u = Nd
(23)
Reemplazando (23) en (20) y ese resultado a su vez en la componente de trabajo
virtual interno para un elemento extraída de (22) tenemos:
=
ε
δ Winte
=
∂ Nd
=
B
∂N
→=
ε
Bd
(24)
d B CB
δd dV d  ∫ B CB dV  δd = d K δd
∫=


T
T
T
T
Ve
T
e
Ve
(25)
e
donde K es la matriz de rigidez del elemento.
dV
K e = ∫ B TCB
(26)
Ve
Trabajando sobre el trabajo virtual externo del elemento, tenemos:
(
)
δ Weex t =∫ t T δ u d S + ∫ f T δ u d V = ∫ t T N d S + ∫ f T N d V δ d =P e δ d
Se
Ve
Se
Ve
(27)
e
donde P representa el vector de fuerzas del elemento.
=
Pe
∫
Se
t T N dS + ∫ f T N dV
Ve
(28)
Ensamblando las matrices de rigidez de los elementos K e y los vectores de fuerzas
nodales P e de los elementos, se conforma el sistema global para el cuerpo:
Kd = P
3.5
(29)
Selección del elemento finito y de las funciones de forma
El análisis que se realiza en este trabajo busca obtener resultados que se aproximen
aceptablemente a los reportados por el proyecto BOJCAS [4] y a los obtenidos con los
modelos numéricos desarrollados por McCarthy [5], Ekh [16] y otros autores dentro del
mismo.
30
El software ABAQUS® permite modelar placas de material compuesto como láminas
(shell ) o como elementos sólidos tridimensionales. Song [17] indica que una placa de
compuesto generalmente se analiza mediante el uso de elementos de lámina (bidimensionales) pero que en ciertos casos se requiere el uso de elementos tridimensionales, eso
ocurre cuando:
• los efectos cortantes transversales son importantes.
• las tensiones normales no pueden ser ignoradas.
• se necesita mejorar las definiciones de contacto.
Sobre esa base, se ha optado por modelar el conjunto analizado en el presente trabajo
con elementos tridimensionales ya que, como se muestra más adelante, se requiere
establecer definiciones de contacto entre sus componentes, además de obtener desplazamientos en el sentido transversal en puntos específicos del espesor de las placas.
El elemento seleccionado para el modelado es un hexaedro cuadrático, ya que se adapta
mejor a la geometría que el hexaedro lineal y además su uso está recomendado en el
Manual de Usuario de ABAQUS® [3]. Se adoptó ese elemento (hexaedro cuadrático)
principalmente porque brinda mejores resultados a un menor costo computacional, y
además se adapta de mejor manera a geometrías complejas.
El hexaedro contiene 20 nodos, con tres grados de libertad de desplazamiento en
cada uno, dando un total de 60 grados de libertad. Se trata de un elemento cuadrático,
que además cuenta con 8 puntos de Gauss para la integración en su interior (integración
reducida). Se adoptó esa opción debido a que capta de mejor manera los desplazamientos
transversales, la flexión y la torsión existentes en el modelo motivo de análisis, reduciendo
el costo computacional del análisis. El elemento es llamado C3D20R en el programa
ABAQUS® y se muestra en la Figura 3-3:
Figura 3-3: Elemento hexaédrico cuadrático, C3D20R
31
Las coordenadas de los nodos en el sistema local del elemento ( ξ - η - ζ ) cuyo
origen está en el centro del mismo, están dadas en la Tabla 3-2.
Tabla 3-2: Coordenadas de los nodos en el elemento hexaédrico serendípito cuadrático
Nodo
1
2
3
4
5
Coordenada
( 1, –1,–1)
( 1, 1, –1)
(–1, 1,–1)
(–1,–1,–1)
( 1, –1, 1)
Nodo Coordenada
6
( 1, 1, 1 )
7
(–1, 1, 1)
8
(–1,–1, 1)
9
( 1, 0,–1)
10
( 0, 1,–1)
Nodo Coordenada
11
(–1,0,–1)
12
(0,–1,–1)
13
( 1, 0, 1)
14
( 0, 1, 1)
15
(–1, 0, 1)
Nodo Coordenada
16
( 0,–1, 1)
17
( 1,–1, 0)
18
( 1, 1, 0)
19
(–1, 1, 0)
20
(–1,–1, 0)
Las funciones de forma de este elemento hexaédrico serendípito (elemento cuadrático
donde se condensan los nodos de las caras y el nodo centroidal ), son las estándar que se
listan en (30), empleadas en todos los códigos de elementos finitos.
N1 =
(1 + ξ ) (1 − η ) (1 − ξ ) ( ξ − η − ζ − 2 ) / 8 N 2 =+
(1 ξ )(1 + η )(1 − ζ ) ( ξ + η − ζ − 2 ) / 8
N 3 =−
(1 ξ )(1 + η ) (1 − ζ ) (−ξ + µ − ζ − 2) / 8
N 4 =−
(1 ξ )(1 − η )(1 − ξ ) (−ξ − η − ζ − 2) / 8
N 5 =+
(1 ξ )(1 − η ) (1 + ζ ) ( ξ − η + ζ − 2 ) / 8 N 6 =+
(1 ξ )(1 + η )(1 + ζ ) ( ξ + η + ζ − 2 ) / 8
N 7 =−
(1 ξ )(1 + η )(1 + ζ ) (−ξ + η + ζ − 2) / 8
N8 =−
(1 ξ )(1 − η )(1 + ζ ) (−ξ − η + ζ − 2) / 8
N 9 =+
(1 ξ ) (1 − η 2 ) (1 − ζ ) / 4
N10 =−
(1 ξ 2 ) (1 + η )(1 − ζ ) / 4
N11 =−
(1 ξ ) (1 − η 2 ) (1 − ζ ) / 4
N12 =
(1 − ς 2 ) (1 − η )(1 − ζ ) / 4
N13 =
(1 + ξ ) (1 − η 2 ) (1 + ζ ) / 4
N14 =
(1 − ξ 2 ) (1 + η )(1 + ζ ) / 4
N15 =
(1 − ξ ) (1 − η 2 ) (1 + ζ ) / 4
N16 =
(1 − ξ 2 ) (1 − η )(1 + ζ ) / 4
N17 =+
(1 ξ )(1 − η ) (1 − ζ 2 ) / 4
N18 =+
(1 ξ )(1 + η ) (1 − ζ 2 ) / 4
N19 =−
(1 ξ )(1 + η ) (1 − ζ 2 ) / 4
N 20 =−
(1 ξ )(1 − η ) (1 − ζ 2 ) / 4
(30)
Estas funciones de forma también son usadas para la interpolación de las coordenadas
de los nodos del elemento, dado que es un elemento isoparamétrico, lo cual para el presente
caso queda de la siguiente manera:
=
x
20
N i (ξ ,η , ζ ) xi
y
∑=
20
=
N i (ξ ,η , ζ ) yi
z
∑
20
∑ N (ξ ,η , ζ ) z
i
i
(31)
=i 1 =i 1 =i 1
La matriz de interpolación N para el elemento hexaédrico de 20 nodos con 3 grados de
libertad de desplazamiento en cada uno, es de dimensión (3 x 60) y tiene la siguiente forma:
N
0 
 N1 0 0 N 2 0 0 N 3 0 0 N 4 0 0 …… N 20 0
 0 N 0 0 N 0 0 N 0 0 N 0 …… 0 N
0 
1
2
3
4
20

 0 0 N1 0 0 N 2 0 0 N 3 0 0 N 4 …… 0
0 N 20 
32
(32)
El vector de desplazamientos nodales d de la ecuación (23), el cual conjuntamente
con la matriz N permite determinar los desplazamientos en los elementos finitos, tiene
dimensión [60 x 1] y está dado por:
 u1 
 v1 
 w1 
 u2 
 v2 
 
d =  w2 
  
  
 u20 
 v20 

 w20 

(33)
La matriz B definida en (24) desarrollada a partir de las matrices ∂ y N , tiene la
siguiente forma:
B
∂

 ∂x 0 0 


0 ∂ 0


∂y


0 0 ∂

∂z 

∂ ∂
0

∂z ∂y 

∂
∂
0


∂x 
 ∂z
∂ ∂ 0
 ∂y ∂x



0 
 N1 0 0 N 2 0 0 N 3 0 0 …… N 20 0
 0 N 0 0 N 0 0 N 0 …… 0 N
0  (34)
1
2
3
20


0 N 20 
 0 0 N1 0 0 N 2 0 0 N 3 …… 0
Efectuando el producto indicado en (34) se obtiene:
∂N 3
∂N 2
 ∂N1
0
0
0
0
0
0
 ∂x
∂x
∂x

∂N 3
∂N 2
 0 ∂N1 0
0
0
0
0

∂y
∂y
∂y

∂N 3
∂N1
∂N 2
 0
0
0
0
0
0
∂z
∂z
∂z
B=
∂N 3 ∂N 3
∂N1 ∂N1
∂N 2 ∂N 2

0
0
 0
∂z ∂y
∂z ∂y
∂z ∂y

 ∂N1 0 ∂N1 ∂N 2 0 ∂N 2 ∂N 3 0 ∂N 3
 ∂z
∂x ∂z
∂x ∂z
∂x
 ∂N ∂N
∂N 3 ∂N 3
∂N 2 ∂N 2
1
 1
0
0
0
∂y ∂x
∂y ∂x
 ∂y ∂x
33
...
∂N 20
∂x
...
0
∂N 20
∂y
...
0
0
...
0
∂N 20
∂z
0
∂N 20
0
∂z
∂N 20 ∂N 20
...
∂y
∂x
...

0 

0 

∂N 20 

∂z 
(35)
∂N 20 
∂y 

∂N 20 
∂x 

0 

Dado que las funciones de forma mostradas en (30) están definidas en términos de
las coordenadas naturales ‘ξ ’, ‘η’ y ‘ζ ’, sus derivadas respecto a ‘x’, ‘y’ y ‘z’ se obtienen
aplicando la regla de la cadena para derivadas parciales, lo cual se expresa matricialmente
como:
=
G








∂x
∂ξ
∂x
∂η
∂x
∂ζ
∂y
∂ξ
∂y
∂η
∂y
∂ζ
∂z
∂ξ
∂z
∂η
∂z
∂ζ
∂N i 
 ∂ξ 


 ∂N i 
=
→


 ∂η 
 ∂N i 
 ∂ζ 










 ∂N i 
 ∂x 
∂ 
 N 
G  i
 ∂y 
 ∂N i 

 ∂z 

(36)
donde G es la matriz jacobiana de la transformación, la cual generalmente se denota con
la letra J. La matriz jacobiana también puede relacionarse con las funciones de
interpolación a través de reemplazar las coordenadas x, y y z por las relaciones
planteadas en (30), quedando de la forma:
20
20
 20 ∂N i
∂N i
∂N i
x
y
zi
∑ i
∑
∑
i
∂ξ
∂ξ
∂ξ
=
i 1=
i 1
 i 1=
20
20
 20 ∂N
∂N i
∂N i
i
G =  ∑ xi
yi
zi
∑
∑
∂η
∂η
∂η
 i 1=
=
i 1
i =1
 20
20
20
∂N i
∂N i
∂N i

x
y
zi
∑
∑
∑
i
i
 i 1=
∂ζ
∂ζ
∂ζ
=
i 1
i =1











(37)
Invirtiendo la matriz G dada en (37) y reemplazando en (36) se obtiene:
 ∂N i
 ∂x

 ∂N i

 ∂y
 ∂N
 i
 ∂z
 ∂N i

 ∂ξ





− 1  ∂N i
 = G 

 ∂η

 ∂N

 i

 ∂ζ












(38)
Esta expresión permite desarrollar la matriz B mostrada en (34) en función de las
coordenadas naturales, a partir de lo cual se puede determinar la matriz de rigidez K e
del elemento tridimensional que figura en las ecuaciones (25) y (26) y que se repite a
continuación como ecuación (39):
Ke =
1 1 1
K =
e
∫∫ ∫
∫
Ve
B TCB
dV
B TCB
det G d ξ dη d ζ
−1−1−1
34
(39)
(40)
CAPÍTULO 4
4
4.1
ANÁLISIS NUMÉRICO DE UNIONES DE
MATERIAL COMPUESTO ABULONADAS
Introducción
La vinculación de piezas de material compuesto es un ítem que requiere especial
atención durante el diseño de estructuras, especialmente en la industria aeroespacial.
Esas vinculaciones suelen tener alta influencia sobre la capacidad de soportar cargas y
transmitirlas. En el caso específico de las uniones abulonadas, debido a la presencia del
bulón y del agujero en el material, las tensiones y deformaciones pueden variar tridimensionalmente. Los factores que intervienen en ese tipo de vinculaciones son: la flexión de
las partes unidas y la inclinación que experimenta el bulón causada por la aplicación de la
carga excéntrica. Esos dos aspectos afectan directamente al compuesto, en el que
pueden aparecer modos de falla asociados a los mismos.
El análisis de uniones de material compuesto se ha encarado por métodos analíticos y
por medio del análisis numérico, esto último ha sido realizado en el pasado mayoritariamente con modelos bidimensionales. El desarrollo computacional de los últimos años
ha hecho que sea más habitual llevar a cabo análisis tridimensionales, los cuales
inicialmente se realizaban apilando elementos sólidos ortótropos que representaban las
capas del compuesto. Además, por supuesto, se empleaban procedimientos experimentales
que permitían obtener propiedades equivalentes de las placas de material. En la actualidad,
a través de módulos específicos para materiales compuestos, estas propiedades equivalentes son incorporadas en los códigos comerciales de resolución por elementos finitos.
En el año 2001 se desarrolló el proyecto “BOJCAS - Bolted joints in composite aircraft
structures” [4], en el que se analizó el comportamiento de uniones abulonadas de materiales
compuestos. En ese proyecto se realizaron procedimientos tanto experimentales como
numéricos, aplicando cargas de tracción sobre conjuntos de placas abulonadas con solape
simple. Posteriormente y aun como parte de ese proyecto, McCarthy, et al. [5]
desarrollaron modelos numéricos para contrastarlos con los resultados experimentales
del proyecto BOJCAS. En el trabajo de McCarthy se utilizó como base las mismas
formas, disposiciones y componentes de los materiales ensayados experimentalmente en
el proyecto BOJCAS.
El modelo que se desarrolla en este capítulo está basado en los materiales y características
geométricas utilizados en los trabajos antes mencionados. El modelo se basa en el empleo
del software comercial de análisis por elementos finitos ABAQUS®, y los resultados del
mismo se contrastan con los obtenidos en el proyecto BOJCAS.
35
Para la comparacion se han empleado los resultados experimentales y los numéricos.
El objeto de la comparación es validar el modelo numérico desarrollado en este trabajo,
para luego en base al mismo implementar modelos de mayor cantidad de bulones y con
diferente disposición de los mismos. El objetivo final es realizar una predicción del
comportamiento y eficiencia de las configuraciones habituales de abulonado.
4.2
Descripción del problema
En el año 2011, Morales y Quadri [6] realizaron un trabajo experimental en la
Universidad Nacional de Córdoba para caracterizar uniones mecánicas de materiales
compuestos. En ese trabajo se estudiaron los efectos del cambio en la cantidad y
disposición de bulones en la unión de materiales compuestos. Morales y Quadri
emplearon en su estudio compuestos de fibra de carbono y fibra de vidrio en matriz
epoxi. Ese trabajo despertó interés por estudiar en esta tesis la influencia del número y
disposición de bulones sobre las características de la unión. Sin embargo por su carácter
experimental, ese trabajo no presenta todos los datos necesarios para desarrollar y
validar un modelo numérico.
El trabajo de McCarthy en cambio, si bien no tiene el mismo objetivo que el de
Morales y Quadri, provee suficiente información para elaborar el modelo numérico
adecuado para el presente trabajo. La unión analizada en el trabajo de McCarthy, es de
solape simple con un único bulón centrado en uno de los extremos. La geometría es la
mostrada en la Figura 4-1, y está basada en la norma ASTM D5961 [18] que regula las
relaciones del diámetro del bulón (d ), con el ancho (a), la distancia al extremo (e) y
espesor (t); de la manera indicada en la figura:
a
=
d
48
= 6
8
e
=
d
24
= 3
8
d
=
t
8
= 1, 6
5, 2
Figura 4-1: Geometría de la probeta con un solo bulón, unidades en milímetros
36
Estas relaciones entre las dimensiones, según lo indica McCarthy tienen como
objetivo inducir a una falla por aplastamiento en la junta.
El material usado por McCarthy es un compuesto de fibra de carbono en resina epoxi
elaborada por Hexcel Composites, de amplia aplicación en la industria aeroespacial, cuyo
nombre comercial es HTA/6376. La fibra es del tipo unidireccional, con un espesor
nominal de 0,13 mm y cuyas propiedades de rigidez se describen en la Tabla 4-1.
Tabla 4-1: Propiedades de rigidez unidireccional para fibra de carbono-epoxi HTA / 6376 en GPa
E 11
E 22
E 33
G 12
G 13
G 23
ν 12
ν 13
ν 23
140
10
10
5,2
5,2
3,9
0,3
0,3
0,3
Entre las secuencias de capas de compuesto utilizadas por McCarthy se encuentra la
disposición cuasi isótropa [45,0,–45,90] 5s , la cual se ha empleado en este trabajo.
Los resultados experimentales para el caso descripto en los párrafos anteriores se
encuentran bien documentados y sirven luego para contrastar los resultados numéricos
obtenidos en este trabajo de tesis.
4.3
Modelo de elementos finitos
El modelo de elementos finitos empleado en este trabajo comprende cinco elementos:
dos placas de material compuesto laminado, dos arandelas y una combinación bulón–
tuerca; similar al modelado realizado en el trabajo de McCarthy.
La placa fue modelada según las dimensiones mencionadas en la Sección 4.2, sin
omitir ninguna parte de la misma dentro del modelo. La arandela fue modelada de
acuerdo a las dimensiones utilizadas por el autor del trabajo de referencia. El bulón y la
tuerca fueron modelados en un solo cuerpo, del mismo modo que se lo hace en [5]. Para
el caso de dimensiones no especificadas dentro del trabajo [5], correspondientes al
bulón y arandela, se han tomado las mismas de un catálogo normalizado basado en las
normas ISO-4017 y norma ISO-7091 respectivamente.
4.3.1
Malla de elementos finitos
Cada componente del modelo tiene una malla característica. En la arandela, se
emplearon elementos hexaédricos, considerando un solo elemento en el espesor y treinta
y dos elementos igualmente espaciados en el perímetro como se muestra en la Figura 4-2.
El conjunto bulón-tuerca, es modelado como un solo cuerpo como se puede observar
en la Figura 4-3. Se emplearon elementos de forma hexaédrica y tetraédrica. Las
secciones que representan la cabeza del bulón y la tuerca, tienen treinta y dos elementos
en su perímetro circular, tres elementos en la dirección radial y un solo elemento en el
espesor de la cabeza del bulón y la tuerca.
37
Figura 4-2: Mallado de la arandela
La sección que representa el cuerpo del bulón, fue mallada con treinta y dos elementos
igualmente espaciados en el perímetro circular, dos elementos radialmente colocados y
diez elementos a lo largo del vástago. Se ha aplicado una malla más densa al vástago del
bulón dado que es la parte que entra en contacto con las placas de laminado durante el
ensamblaje del modelo.
Figura 4-3: Mallado y corte transversal del conjunto bulón-tuerca
En el caso de la placa que representa al material compuesto se utilizó un mallado de
mayor densidad en algunos lugares específicos, como se muestra en la Figura 4-4. Esto
se debe a la necesidad de contar con valores de deformación en puntos definidos en el
trabajo de McCarthy. Este mallado más denso se encuentra en el sector aledaño al
orificio de abulonado, que es donde se presenta el contacto entre los cuerpos
involucrados en el análisis.
38
Figura 4-4: Mallado de la placa de material compuesto
4.3.2
Tipo de elementos utilizados
Basado en experiencias anteriores se adoptó para las placas un elemento cuadrático
hexaédrico de veinte nodos e integración reducida denominado C3D20R porque permite
obtener buenos resultados.
Por otra parte, para las arandelas y las secciones que usan elementos hexaédricos en el
bulón y la tuerca, se emplearon elementos del tipo C3D8R. Este último elemento es un
hexaedro tridimensional de ocho nodos con interpolación lineal e integración reducida.
Mientras que en el núcleo del vástago del bulón se usaron elementos prismáticos
triangulares denominados C3D6, de seis nodos con aproximación lineal e integración
reducida, dado que permite captar mejor la geometría del núcleo del bulón.
4.3.3
Condiciones de borde y de carga
En el ensayo de tracción llevado a cabo como parte del proyecto BOJCAS, se empleó
una máquina cuyas mordazas actúan sobre un área de 3600 mm2 (75 x 48mm) en los
extremos de la probeta. Esta característica es importante dado que a través de estas
superficies se produce la transmisión de cargas desde la máquina de tracción hacia la
probeta ensayada. El modelo incorpora las zonas de enganche como se muestra en la
Figura 4-5Figura 4-5. En esas áreas se han impuesto las siguientes condiciones de
borde:
i) Restricción de los desplazamientos en los tres grados de libertad (u=v=w=0) para
la placa que se acopla al cabezal fijo de la máquina de ensayos (izquierda); y
ii) Restricción de los desplazamientos en las direcciones ‘y’ y ‘z’ en la placa acoplada
al cabezal móvil (que tiene libertad de desplazarse según el eje longitudinal de la
máquina, eje ‘x’) como se indica a la derecha de la Figura 4-5.
39
Figura 4-5: Condiciones de borde impuestas al modelo
A fin de simular adecuadamente la unión, el bulón tiene una precarga axial (carga de
7,2 MPa) tal como se ha indicado en el trabajo de McCarthy. Para ello se ha aplicado
una opción prevista en el software que está basado en la publicación realizada por
Milligan [19]. Esta carga fue colocada en el centro del bulón, en dirección axial y como
se muestra en la Figura 4-6.
Figura 4-6: Detalle de la zona de solapamiento y precarga del bulón
La carga de tracción en la dirección ‘x’ del conjunto estudiado, se aplica sobre la
placa que está vinculada al cabezal móvil y varía de 0 a 5 kN. Se adoptó ese valor de
carga con el objetivo de mantener a la probeta en el rango elástico durante todo el
ensayo. La carga es aplicada en la placa derecha, que tiene libertad de movimiento en la
dirección x, como se muestra en la Figura 4-7.
40
Figura 4-7: Ubicación de la carga de tracción coloreada en rojo
La probeta a ensayar tiene varias partes que entran en contacto: i) las placas de
material compuesto entre sí, ii ) el bulón con las placas laminadas en el orificio de
abulonado, iii) las arandelas con las placas de laminado, y iv) las arandelas con el
conjunto bulón-tuerca. Se definieron las superficies de contacto identificando en cada
pieza del conjunto las zonas donde efectivamente se producirá la interferencia.
La condición de contacto para las superficies involucradas contempla un coeficiente
de rozamiento µ = 0,5. Este valor se estimó partiendo de las características de los materiales
en contacto, según los valores típicos descriptos en el trabajo de Dominique y colab. [20].
Para el caso del bulón se consideró la sugerencia de definición de contacto provista por
Milligan [21]. Es necesario considerar que una buena definición de los contactos entre
las partes involucradas permite obtener resultados más próximos a la realidad, así como
la correcta transmisión de las fuerzas entre las partes en contacto.
4.3.4
Modelo constitutivo del material
Como se mencionó previamente, el material analizado en el trabajo de McCarthy y
en el proyecto BOJCAS está fabricado a partir de fibra de carbono en matriz epoxi
(HTA/6376). Las propiedades de la fibra unidireccional se presentan en la Tabla 4-1. En
el modelo numérico presentado en esta tesis, de manera similar a lo propuesto en [5], se
utilizaron dos métodos para representar el material: i ) un material compuesto laminado
donde la disposición y orientación de las placas es parte de los datos del problema, y ii )
un material sólido con características de ortotropía donde los datos del problema son las
propiedades aparentes del material.
41
Al definir el material laminado se contemplan las capas de material dispuestas de
acuerdo a lo realizado en el experimento. Se han considerado cuarenta capas con las
propiedades dadas en la Tabla 4-1 y una configuración cuasi isótropa [45,0,–45,90] 5s
utilizando las opciones previstas por ABAQUS® para ese efecto. En el caso de
considerar las propiedades mecánicas equivalentes, al igual que en [4] y [5], se han
aplicado las propiedades que se muestran en la Tabla 4-2.
Tabla 4-2: Propiedades de rigidez ortótropa equivalente para
la configuración [45,0,-45,90] 5s de fibra HTA/6376 en [GPa]
E xx
E yy
E 33
G xy
G xz
G yz
ν xy
ν xz
ν yz
54,25
54,25
12,59
20,72
4,55
4,55
0,309
0,332
0,332
Los otros componentes de la junta fueron modelados con materiales isótropos, el
conjunto bulón-tuerca se considera fabricado en titanio y las arandelas construidas con
acero. Las propiedades utilizadas para estos elementos se muestran en la Tabla 4-3.
Tabla 4-3: Propiedades de los materiales isótropos utilizados en el modelo
Material
E [GPa]
ν
Acero
Titanio
210
110
0,30
0,29
Resulta importante destacar que las propiedades de rigidez equivalente que se
presentan en la Tabla 4-2, son muy similares a las que se pueden obtener al despejarlas
de los elementos de la matriz de flexibilidad del sub-laminado [45º/0º/–45º/90º] que se
presentan en la Tabla 3-1 del Capítulo 3. Las pequeñas diferencias que se aprecian pueden
ser atribuidas a que las propiedades equivalentes deducidas por McCarthy en su trabajo
fueron obtenidas por medios analíticos y afinadas mediante procesos experimentales
donde pueden influir factores externos.
4.4
Validación numérica y experimental del modelo propuesto
Para realizar la validación del modelo propuesto en este trabajo, se compararon los
resultados obtenidos (este trabajo) con los experimentales [4] y numéricos [5]. Los
resultados mostrados en [5] han sido cotejados contra los dados en BOJCAS, siguiendo
en particular dos grupos de variables:
1. Las deformaciones en puntos específicos de la superficie de la placa de compuesto.
2. La rigidez de la unión.
En la experimentación se instrumentaron ocho galgas extensiométricas, ubicadas en
el sector del abulonado y dispuestas como se muestra en la Figura 4-8. La deformación
experimental medida en las 8 galgas extensiométricas se ha graficado en la Figura 4-9.
42
Figura 4-8: Disposición de las galgas extensiométricas ( dimensiones en milímetros )
Esas mismas ubicaciones de las galgas han sido consideradas en el modelo numérico
desarrollado en el presente trabajo, asegurando la existencia de un nodo en esos puntos.
El objetivo de eso es validar los modelos desarrollados en este trabajo, del mismo modo
que se hace en [5].
4.4.1
Deformaciones en la superficie de la placa
El modelo numérico desarrollado fue solicitado hasta un nivel de carga de 5 kN, y se
calcularon los desplazamientos de los nodos ubicados en los lugares especificados para
las galgas extensiométricas y en las direcciones correspondientes (transversal para la
galga Nº 7 y axial para las demás). A partir de los resultados se obtienen las gráficas
carga-deformación que son contrastadas con las gráficas experimentales (Figuras 4-10
hasta 4-13) y con las obtenidas por McCarthy (Figuras 4-14 hasta 4-18).
Teniendo en cuenta lo observado en el trabajo realizado en la Universidad Nacional
de Córdoba por Morales y Quadri [3], y ratificado en el modelo numérico ejecutado en
ABAQUS®, es la placa sobre la que se aplica la carga la que presenta el mayor nivel de
deformación. Por este motivo los cálculos se realizan sobre la placa en la que se aplica
la carga y es donde se analizan los resultados.
En las placas modeladas en este trabajo como un material sólido con propiedades
ortótropas equivalentes, se obtienen los resultados que se muestran en la Tabla 4-. En
esa tabla también se han listado los resultados de los modelos experimentales [4] y los
numéricos obtenidos en [5].
43
Tabla 4-4: Micro deformaciones obtenidas para un material solido ortótropo con propiedades equivalentes
Galga
Modelo experimental [4]
Modelo numérico [5]
Este trabajo
1
2
3
4
5
6
7
8
–1.8
760
–349
–488
–400
–218
–367
–353
149
633
–244
–438
–346
–182
–414
–346
84 ,7
717 ,1
–315 ,5
–491 ,8
–363 ,0
–227 ,7
–410 ,4
–363 ,0
Figura 4-9: Deformación experimental medida en las galgas extensiométricas
Como se muestra en los gráficos superpuestos de las Figuras 4-10 hasta 4-13, y
también se evidencia en los datos de la Tabla 4-4, el modelo de elementos finitos
desarrollado en este trabajo capta adecuadamente los desplazamientos de las galgas
extensiométricas. La mayor diferencia respecto a los resultados experimentales
BOJCAS, se observa en la deformación medida por la galga Nº 1. Sin embargo también
puede verse, que de acuerdo con los resultados experimentales esa galga prácticamente
no presenta desplazamiento durante la aplicación de la carga. Esto último puede deberse
a un problema durante el proceso de ensayo.
Las diferencias también se observan en la comparación de los resultados de este
trabajo con los numéricos del trabajo de McCarthy. Sin embargo comparativamente los
resultados numéricos se muestran más cercanos entre sí, inclusive en la galga N° 1. En
las demás galgas la máxima diferencia relativa en valor absoluto resulta del orden del
12 %, y este valor está asociado a la galga Nº 7.
44
5
Galga Nº 1
Carga (KN)
Carga (KN)
5
4
3
2
Galga Nº 2
4
3
2
1
1
Experimental
Experimental
Numérico actual
0
-200
0
200
400
Numérico actual
0
600
0
200
Microdeformación
400
600
800
Microdeformación
Figura 4-10: Deformaciones en las galgas Nº 1 y Nº 2
5
Galga Nº 3
Carga (KN)
Carga (KN)
5
4
3
2
Galga Nº 4
4
3
2
Experimental
1
0
Experimental
1
Numérico actual
-800
-600
-400
-200
0
0
Numérico actual
-800
-600
-400
-200
0
Microdeformación
Microdeformación
Figura 4-11: Deformaciones en las galgas Nº 3 y Nº 4
5
Galga Nº 5
4
Carga (KN)
Carga (KN)
5
3
2
Galga Nº 6
4
3
2
1
1
Experimental
Experimental
Numérico actual
0
-800
-600
-400
Numérico actual
-200
0
0
-800
-600
Microdeformación
-400
-200
Microdeformación
Figura 4-12: Deformaciones en las galgas Nº 5 y Nº 6
45
0
5
Galga Nº 7
Carga (KN)
Carga (KN)
5
4
3
Galga Nº 8
4
3
2
2
Experimental
1
Experimental
1
Numérico actual
0
-800
-600
-400
Numérico actual
-200
0
0
-800
-600
Microdeformación
-400
-200
0
Microdeformación
Figura 4-13: Deformaciones en las galgas Nº 7 y Nº 8
Figura 4-14: Deformaciones en las galgas del modelo numérico de McCarthy
5
Galga Nº 1
4
Carga (KN)
Carga (KN)
5
3
2
Galga Nº 2
4
3
2
1
1
Numérico McCarthy
Numérico McCarthy
Numérico actual
0
-200
0
200
400
Numérico actual
0
600
Microdeformación
0
200
400
600
Microdeformación
Figura 4-15: Deformaciones obtenidas numéricamente en las galgas Nº 1 y Nº 2
46
800
5
5
Carga (KN)
Carga (KN)
Galga Nº 3
4
3
2
Galga Nº 4
4
3
2
1
1
Numérico McCarthy
Numérico McCarthy
Numérico actual
0
-800
-600
-400
Numérico actual
-200
0
0
-800
-600
Microdeformación
-400
-200
0
Microdeformación
Figura 4-16: Deformaciones obtenidas numéricamente en las galgas Nº 3 y Nº 4
5
5
4
4
3
3
2
2
1
1
Numérico McCarthy
Numérico actual
0
Galga Nº 6
Carga (KN)
Carga (KN)
Galga Nº 5
-800
-600
-400
-200
0
0
Numérico McCarthy
Numérico actual
-800
-600
-400
-200
0
Microdeformación
Microdeformación
Figura 4-17: Deformaciones obtenidas numéricamente en las galgas Nº 5 y Nº 6
5
Galga Nº 7
Carga (KN)
Carga (KN)
5
4
3
2
Galga Nº 8
4
3
2
1
1
Numérico McCarthy
Numérico McCarthy
Numérico actual
0
-800
-600
-400
Numérico actual
-200
0
0
Microdeformación
-800
-600
-400
-200
Microdeformación
Figura 4-18: Deformaciones obtenidas numéricamente en las galgas Nº 7 y Nº 8
47
0
4.4.2
Rigidez y deformaciones de la unión
En lo que refiere a la rigidez de la unión, los valores se muestran en la Tabla 4-5. En
esa tabla se observa que el modelo numérico desarrollado en este trabajo produce un valor
muy cercano al experimental, mostrando una diferencia de tan solo 5,8 %.
Tabla 4-5: Valores de rigidez de la junta abulonada en kN / mm
Experimental [4]
Numérico [5]
Diferencia %
Este trabajo
Diferencia %
28
31,5
12,5 %
29,64
5,8 %
En el modelo de McCarthy [5] y en el modelo experimental BOJCAS [4], se
determina el valor de la deformación axial, a través del promedio de los valores de
deformación obtenidos en las galgas Nº 1 y Nº 2, mientras que la deformación flexional
se obtiene restando al valor de deformación de la galga Nº 2, la deformación de la galga
Nº 1 y se divide ese resultado por dos. Matemáticamente se tiene:
ε axial
=
ε1 + ε 2
ε flexional
=
2
ε 2 − ε1
(1)
2
Aplicando la misma metodología a los valores obtenidos con el modelo numérico de
este trabajo, se obtienen los resultados que se muestran a continuación en la Tabla 4-6.
Tabla 4-6: Microdeformación de la probeta
Experimental [4]
Numérico [5]
Diferencia
Este trabajo
Diferencia
Axial
379
391
3,1 %
401
5,8 %
Flexional
381
242
36,5 %
316
17,1 %
Como se puede apreciar en la Tabla 4-6, a pesar de la diferencia que se observó en la
galga Nº 1, los valores de deformación obtenidos con el modelo de este trabajo guardan
una buena correlación con los experimentales. La diferencia en la deformación axial es
del 6 %, mientras que para la deformación flexional, se ubica en el 17 %.
4.4.3
Desplazamiento del plano medio de la placa
Para la comparación del modelo elaborado con capas apiladas, en el trabajo de
McCarthy [5], se emplean los trabajos numéricos llevados a cabo por Anderson [22] y
Ekh [16]. En esos trabajos los autores han empleado mallas altamente refinadas, con
elementos de cuarto orden, alcanzando hasta 1,2 millones de grados de libertad.
La comparación se realiza inicialmente entre los desplazamientos en la dirección ‘z’,
que sufre el plano central de la placa y son medidos en los lados de la misma. En los
modelos numéricos de cuarto orden se realiza esta medición cuando la junta sometida a
tracción, se desplaza 0,5 mm en la dirección ‘x’ (axial). El sistema coordenado usado
por el modelo de Anderson se presenta en la Figura 4-19.
48
Figura 4-19: Líneas en los bordes y arista en el orificio donde se han evaluado tensiones
La Figura 4-20 muestra las curvas obtenidas por Ekh [16], que reflejan el comportamiento del plano medio de la probeta abulonada. En esa misma figura se comparan esos
resultados con los obtenidos en esta tesis:
Desplazamiento Z [mm]
0.6
Ver esquema aclaratorio en la Figura 4-19
0.5
Lado 1 - Ekh
0.4
Lado 1 actual
0.3
Lado 2 - Ekh
0.2
Lado 2 actual
0.1
0
-0.1
-0.2
-0.3
0
10
20
30
40
50
60
70
80
Distancia sobre el plano medio [mm]
Figura 4-20: Modelo de Ekh vs modelo actual,
desplazamiento de los extremos del plano medio en la dirección z
La tendencia obtenida con el modelo numérico de este trabajo es la misma que
presenta el modelo de Ekh, presentándose un valor de desplazamiento diferente para
cada lado de la probeta, lo que confirma la presencia de torsión observada en [16] y
[22]. Esto también es mencionado en [5]. En la se presentan los resultados reportados
por [5], y se los compara con los calculados en este trabajo. Esta comparativa muestra
que los resultados guardan la misma tendencia, pero los resultados dados en [5] están
por debajo de los obtenidos en este trabajo.
49
Desplazamiento Z (mm)
0.6
Ver esquema aclaratorio en la Figura 4-19
0.5
Lado 1 - McCarthy
0.4
Lado 1 actual
0.3
Lado 2 - McCarthy
0.2
Lado 2 actual
0.1
0
-0.1
-0.2
-0.3
0
10
20
30
40
50
60
70
80
Distancia sobre el plano medio (mm)
Figura 4-21: Modelo de McCarthy vs. modelo actual,
desplazamiento de los extremos del plano medio en la dirección z
Los desplazamientos de cada uno de los lados obtenidos por McCarthy, si bien
mantienen la tendencia están más apartados que los obtenidos a través del modelo de
Ekh [16] (ver Figura 4-20) y también del modelo actual.
4.4.4
Distribución de tensiones
El modelo numérico propuesto en este trabajo permite establecer la distribución de
tensiones en cada una de las capas a través del espesor, y a partir de las funciones de
aproximación para el modelo de propiedades equivalentes. Los resultados del trabajo
experimental [6], y de otros autores [23], [24], [25] y [27] que han estudiado uniones
abulonadas, muestran que se produce deformaciones por aplastamiento en la zona de
contacto del bulón con las placas. Con el fin de captar de manera adecuada los valores
de tensión en el orificio, la zona de abulonado cuenta con una malla más fina que el
resto de las partes involucradas en el modelo. La malla en la zona del bulón se muestra
en la Figura 4-22.
Figura 4-22: Detalle del mallado en la zona próxima al orificio de abulonado
50
Para referenciar la ubicación de las tensiones, McCarthy [5] y Anderson [22] en sus
respectivos trabajos han realizado una partición imaginaria de la probeta, paralela al eje
‘x’ pasando por el centro del orificio de abulonado. Esto define un plano de partición x-z
como se muestra en la Figura 4-23 (que reproduce la parte derecha de la Figura 4-19),
ubicando en este plano al orificio de abulonado a lo largo de la altura del semicilindro
formado por el corte. La altura corresponde al espesor de la probeta, tomando como 0
mm el plano de las caras de contacto entre placas y 5,2 mm la cara libre:
Figura 4-23: Referencia para la ubicación de las tensiones
En esta tesis se calcularon las tensiones de dos maneras distintas usando:
i ) el modelo de propiedades ortótropas homogéneas y
ii ) el modelo de laminado por capas (sub-laminados),
los resultados numéricos obtenidos con esos dos modelos se contrastaron con las curvas
obtenidas numéricamente por McCarthy [5] y Anderson [22]. Como en esos dos trabajos
las tensiones están reportadas en forma de gráficos, los valores necesarios para hacer
comparaciones se extrajeron usando el software “Engauge Digitizer ”.
En la Figura 4-24 se comparan las tensiones calculadas del modelo de propiedades
ortótropas homogéneas obtenidas con el programa ABAQUS®, con los resultados del modelo
de cuarto orden de Anderson [22] y los resultados numéricos provistos por McCarthy [5].
5
4
Espesor (mm)
Espesor (mm)
5
Numérico actual
Anderson
3
4
3
2
2
1
1
0
0
-500
-1000 -1500 -2000 -2500 -3000
Numérico actual
McCarthy
0
0
-500
-1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Tensión radial [MPa]
Figura 4-24: Distribución de tensiones en el espesor
51
Las curvas de la Figura 4-24 muestran la misma tendencia y los valores de las tensiones
son similares. Sin embargo los resultados de Anderson, presentan una tendencia hacia
valores del orden 3000 MPa en el plano de contacto entre las placas de compuesto.
En las Figuras 4-25 y 4-26 se comparan los resultados de las tensiones calculadas
numéricamente con el modelo de sublaminados usando el programa ABAQUS®, con los
resultados del modelo numérico de cuarto orden de Anderson [22]. Se presenta la
distribución de tensión en las capas con la misma orientación de fibras.
5
4
Espesor (mm)
Espesor (mm)
5
Anderson 0º
Numérico actual
3
4
3
2
2
1
1
0
0
-500
0
-1000 -1500 -2000 -2500 -3000
Anderson 45º
Numérico actual
0
Tensión radial [MPa]
-500 -1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Figura 4-25: Distribución de tensiones en las capas a 0º y 45º
5
4
Espesor (mm)
Espesor (mm)
5
Anderson -45º
Numérico actual
3
4
3
2
2
1
1
0
0
-500 -1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Anderson 90º
Numérico actual
0
0
-500 -1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Figura 4-26: Distribución de tensiones en las capas a –45º y 90º
En las Figura 4-27 y 4-28 se comparan las tensiones calculadas numéricamente en este
trabajo con el modelo de sublaminados usando el programa ABAQUS®, con las tensiones
calculadas numéricamente por McCarthy [5]. En ambas figuras, todas las orientaciones de
las capas presentan la misma tendencia y valores similares para la tensión radial.
52
5
Espesor (mm)
Espesor (mm)
5
McCarthy 0º
4
Numérico actual
3
2
1
1
0
-500
0
-1000 -1500 -2000 -2500 -3000
Numérico actual
3
2
0
McCarthy 45º
4
0
-500
-1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Tensión radial [MPa]
Figura 4-27: Distribución de tensiones en las capas a 0º y 45º
5
Espesor (mm)
Espesor (mm)
5
McCarthy -45º
4
Numérico actual
3
Numérico actual
3
2
2
1
1
0
0
-500
-1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
McCarthy 90º
4
0
0
-500 -1000 -1500 -2000 -2500 -3000
Tensión radial [MPa]
Figura 4-28: Distribución de tensiones en las capas a –45º y 90º
Como conclusión de las comparaciones realizadas queda demostrado que el modelo
numérico propuesto en este trabajo resulta adecuado para calcular las tensiones en una
unión abulonada de materiales compuestos.
4.4.5
Deformación de la junta
En la Figura 4-29-a se muestra la deformación observada en la unión durante los
ensayos experimentales para un nivel de carga de 5 kN. Además, en la Figura 4-29-b
se presenta la deformada del modelo numérico del trabajo de McCarthy [5]. En la
Figura 4-30 por otra parte se muestra la geometría deformada que se obtuvo con el
modelo numérico desarrollado en este trabajo.
Hay que remarcar que con el modelo de elementos finitos presentado en este trabajo,
se determinaron deformaciones similares a las medidas en el proyecto BOJCAS. Además,
se observa similitud en la inclinación que toma el bulón respecto a su eje longitudinal, y
la flexión causada por la carga excéntrica aplicada a la junta. El trabajo de McCarthy [5]
53
no provee una evaluación cuantitativa de la deformación de la unión, sino que únicamente
la menciona de modo cualitativo (ver Figura 29-b).
Figura 4-29: Deformación de la junta al aplicar 5 kN de carga
(a) Modelo experimental. (b) Modelo numérico de McCarthy [5]
Figura 4-30: Deformación de la junta al aplicar 5 kN en este trabajo
4.5
Resultados numéricos para uniones abulonadas dobles
A partir del modelo desarrollado y validado en la Sección 4.4, se pretende elaborar
dos modelos numéricos similares para estudiar uniones mecánicas con dos diferentes
configuraciones de abulonado: i) dos bulones dispuestos en dirección longitudinal
(axial) y ii) dos bulones dispuestos en dirección transversal. Sabiendo que el modelo de
un solo bulón ha dado resultados acordes a los obtenidos experimentalmente, se espera
que los resultados obtenidos con estos modelos tengan un grado similar de certeza.
4.5.1
Modelo de dos bulones axiales
Al igual que en los estudios del proyecto BOJCAS, la unión es sometida a una carga
de 5 kN. Se han verificado: la rigidez, la deformación axial y la deformación flexional
54
de la unión abulonada. Se estudia además la distribución de tensiones en los agujeros
de la placa, para un desplazamiento prefijado de 0,5 mm del extremo derecho móvil
respecto al extremo izquierdo fijo (ver Figura 4-32), y la distorsión por torsión que sufre
el plano medio de la placa. El dimensionamiento de esta nueva configuración, está
basado en la norma ASTM D5961 [18], y en el trabajo desarrollado por McCarthy [28]
para el caso de una unión de solape simple con múltiples bulones.
4.5.1.1. Geometría de la probeta con dos bulones axiales
La probeta con dos bulones axiales tiene las mismas dimensiones que la probeta de
un solo bulón, se han adicionando 36 mm a su longitud original por la existencia del
segundo bulón. La configuración de las capas del laminado, como la orientación de
fibras es la misma que para el modelo de un solo bulón. La geometría de esta probeta se
muestra en la Figura 4-31. Además la Figura 4-32 presenta una vista en perspectiva del
modelo de elementos finitos empleado en el estudio de esta configuración.
Figura 4-31: Geometría de la probeta con dos bulones axiales (unidades en milímetros)
Figura 4-32: Modelo de elementos finitos de la unión con dos bulones axiales
55
4.5.1.2. Deformaciones en la superficie de la placa y rigidez de la unión
Se han seleccionado dos puntos, similares a la posición de las galgas 1 y 2, que
servirán para determinar la deformación axial y flexional de la junta. A partir de los
resultados correspondientes a las deformaciones axiales en cada lado de la probeta, y
aplicando las ecuaciones (1) de la sección 4.4.2, se obtienen las deformaciones mostradas
en la Tabla 4-7. El máximo desplazamiento sufrido por el conjunto al aplicar la carga de
5 kN es de 0,1253 mm. La rigidez de unión se presenta también en la Tabla 4-7, donde
además se han reproducido las deformaciones y la rigidez asociadas a bulón centrado
(ver Tabla 4-5 y Tabla 4-6).
Tabla 4-7: Deformación y rigidez de la unión con dos bulones axiales (provistos por este trabajo)
Valores comparados
1 Bulón centrado
2 Bulones axiales
Diferencia %
Deformación axial [µε]
400,9
383,94
4,24
Deformación flexional [µε]
316,2
66,63
78,92
Rigidez [kN/mm]
29,64
39,92
34,68
Cuando se comparan los valores de la Tabla 4-7, se observa que si bien la unión con
dos bulones axiales mantiene casi inalterada su deformación axial respecto a un solo
bulón centrado, la deformación flexional con dos bulones es mucho menor (se reduce
aprox. 80 %). Por otra parte comparando la rigidez de la unión con dos bulones axiales
con el caso de un bulón centrado, se observa que la rigidez aumenta en el orden de 35 %.
Se puede decir entonces que el uso de dos bulones axiales mejora el desempeño mecánico
de la unión.
4.5.1.3. Distribución de tensiones en la unión con dos bulones axiales
Se han extraído los valores de tensiones en los orificios de abulonado, tanto para el
modelo representado por capas así como del modelo representado por propiedades
equivalentes. La zona de abulonado presenta una malla más fina que el resto del modelo
de la placa, como se muestra en la Figura 4-33.
Figura 4-33: Detalle de la malla en proximidades de los orificios de la unión con dos bulones axiales
56
Para referenciar la ubicación de las tensiones, se utiliza la misma representación
considerada por McCarthy [5] y Anderson [22], y mostrada en la Figura 4-23. Además para
identificar los orificios de abulonado se ha establecido el esquema de numeración que se
muestra en la Figura 4-34, notar que el bulón 1 es el que está próximo al extremo.
Figura 4-34: Esquema de numeración de los orificios de la unión con dos bulones axiales
En la Figura 4-35 se graficó la distribución de tensiones en los orificios de abulonado
a través del espesor cuando el conjunto es sometido a tracción y se desplaza 0.5 mm,
Espesor (mm)
para el caso del modelo numérico con propiedades homogéneas equivalentes.
5
4
3
Orificio Nº 1
2
Orificio Nº 2
1
0
0
-500
-1000
-1500
-2000
Tensión [MPa]
Figura 4-35: Distribución de tensiones del modelo numérico para la unión con dos bulones axiales
De acuerdo con los valores obtenidos del análisis, el orificio Nº 2 alcanza una tensión
de compresión máxima de 932,95 MPa, un 40 % menos que el caso de un solo bulón
centrado. Por otro lado, el orificio Nº 1 se encuentra en compresión con una tensión
máxima de 529,77 MPa, lo cual representa un 79 % menos que el caso de un solo bulón
centrado. Ambas tensiones son medidas en la superficie que está en contacto con la otra
placa.
57
La distribución de tensiones que se presentan al modelar las probetas por capas
muestra una distribución de tensiones a través del espesor que ha sido representada en
los gráficos de la Figura 4-36.
5
5
o
Orificio N 1
Espesor (mm)
Espesor (mm)
o
4
45º
3
0º
-45º
2
Orificio N 2
4
1
1
0
0
-500
-1000
-1500
0º
-45º
2
90º
0
45º
3
-2000
90º
0
-500
-1000
-1500
-2000
Tensión [MPa]
Tensión [MPa]
Figura 4-36: Distribución de tensiones por capas a través del espesor en el orificio Nº 1 y Nº 2 para la unión
con dos bulones axiales
Las distribuciones por capas muestran resultados similares a los obtenidos cuando se
considera material homogéneo ortótropo, en particular cuando se observa el comportamiento de las placas con orientación longitudinal (0º). Resulta significativo en ambos
casos la descarga de tensiones, que se produce al emplear dos bulones axiales en lugar de
uno solo centrado. Se observa que en el caso del orificio N° 1 la reducción de tensiones en
las fibras de las capas con orientación longitudinal es del orden de 60 %. Esta reducción
en el nivel de tensiones implica una menor interacción entre los bulones y las placas, lo
cual lleva a pensar que habrá una disminución del daño por aplastamiento en las placas
vinculadas con el bulón.
4.5.1.4. Distorsión del plano medio de la placa para dos bulones axiales
Los trabajos desarrollados por McCarthy y otros autores dentro del proyecto BOJCAS,
han utilizado el desplazamiento del plano medio para determinar la presencia de torsión
en las probetas. En el caso de una unión con dos bulones axiales, donde se emplea
nuevamente el sistema de orientación propuesto por Ekh [16] y que se muestra en la
Figura 4-19, se obtienen los desplazamientos que se indican en la Figura 4-37.
Los resultados mostrados en el gráfico de la Figura 4-37 confirman la presencia de
torsión, al igual que sucedió en las probetas con un solo bulón, pero con valores algo
superiores cuando se consideran dos bulones axiales.
58
Desplazamiento Z (mm)
1.2
1
Modelo de dos bulones axiales
0.8
0.6
Lado 1
0.4
Lado 2
0.2
0
-0.2
-0.4
0
20
40
60
80
100
120
Distancia sobre el plano medio (mm)
Figura 4-37: Desplazamiento del plano medio en la dirección z para el modelo de dos bulones axiales
4.5.1.5. Deformación de la unión con dos bulones axiales
Si se observa la Figura 4-38, se ve que la deformada obtenida del análisis numérico
es coincidente con las deformaciones axiales y flexionales mostradas previamente. Se
observa también que la distorsión de los bulones es mucho menor a la que se presenta
con un solo bulón centrado, esto implica menor interferencia con las placas y por lo
tanto también menores tensiones en el espesor de la misma.
Figura 4-38: Deformada final de la unión con dos bulones longitudinales
4.5.2
Modelo de dos bulones transversales
El dimensionamiento de esta nueva configuración, se basa en las normas ASTM D5961
[18] y ASTM D3039 [26]. La distancia entre los bulones se establece de manera que la
probeta mantenga en lo posible proporciones similares a las que se han usado en el caso
de bulones axiales. El ancho de las probetas empleadas previamente no permite incluir
otro bulón, dado que la norma ASTM D5961 exige que la distancia a los extremos debe
ser de 24 mm. Por este motivo se mantiene esta distancia entre los bulones, lo cual
modifica el ancho de la probeta empleada con un solo bulón. Además esta configuración
cumple la norma y también brinda una probeta con el menor ancho posible para esta
59
configuración de dos bulones transversales. En este caso nuevamente la unión con dos
bulones transversales es sometida a una carga de 5 kN, y se verifican los mismos
parámetros y en las mismas condiciones que en la configuración anterior.
4.5.2.1. Geometría de la probeta con dos bulones transversales
La probeta con dos bulones transversales tiene el mismo largo y espesor que la
probeta de un solo bulón centrado, y se ha modificado únicamente su ancho. La
disposición de capas es la misma que para los modelos mostrados previamente. La
geometría de esta probeta se observa en la Figura 4-39. En la Figura 4-40 se muestra
una vista en perspectiva del modelo de elementos finitos empleado en el estudio de la
configuración de dos bulones transversales.
Figura 4-39: Geometría de la probeta con dos bulones transversales (unidades en milímetros)
Figura 4-40: Modelo de elementos finitos de la unión con dos bulones transversales
60
4.5.2.2. Deformaciones en la superficie de la placa y rigidez de la unión
Al igual que en el caso de dos bulones axiales, se han tomado dos puntos equivalentes
a las posiciones de las galgas 1 y 2. Estos puntos sirven para determinar la deformación
axial y flexional de la junta. Se ha considerado que esos puntos están ubicados en ambas
caras de la probeta, a 3 mm del solapamiento en dirección longitudinal y sobre la línea
central de la probeta en dirección transversal. En estos puntos se han evaluado las deformaciones axiales en cada lado de la probeta, y aplicando las ecuaciones detalladas en (1)
en la sección 4.2.2, se obtuvieron los resultados presentados en la Tabla 4-8. También
se han incorporado los resultados de las deformaciones y la rigidez obtenidas para el
caso de un bulón centrado (ver Tabla 4-5 y Tabla 4-6).
Tabla 4-8: Deformación y rigidez de la unión con dos bulones transversales
Valores comparados
Deformación axial [µε]
Deformación flexional [µε]
Rigidez [kN/mm]
1 Bulón Centrado
2 Bulones Transv.
Dif. %
400,9
316,2
29,64
453,49
340,26
32,49
13,1
7,6
9,6
Al contrastar los resultados de la Tabla 4-8, se observa una mayor deformación axial.
Además de manera inversa al caso de dos bulones axiales, las deformaciones por flexión
son levemente mayores. Esto último no es conveniente dado que a mayor deformación
de flexión se presupone una mayor interacción bulón-probeta y por consiguiente mayor
daño. En cuanto a la rigidez de la unión, dado que el máximo desplazamiento sufrido
por el conjunto es 0.1539 mm al aplicar la carga de 5 kN, el valor de rigidez es mayor
que en el caso de un bulón centrado. Sin embargo el aumento en la rigidez de la unión
es del orden de un tercio del que se consigue con dos bulones axiales.
4.5.2.3. Distribución de tensiones en la unión con dos bulones transversales
En la zona de abulonado en la Figura 4-41 se ha empleado nuevamente una malla
más fina para lograr captar adecuadamente las tensiones y deformaciones, provocadas
por la interacción entre el bulón y la probeta.
Figura 4-41: Detalle de la malla en proximidades de los orificios en la unión con dos bulones transversales
61
Nuevamente se utiliza la misma representación considerada por McCarthy [5] y
Anderson [22], que se ha mostrado en la Figura 4-23. En este caso la identificación de
los orificios se estableció con el esquema de numeración indicado en la Figura 4-42.
Figura 4-42: Esquema de numeración de los orificios de la unión con dos bulones transversales
Si se consideran propiedades ortótropas homogéneas, la distribución de tensiones
resultantes en los orificios a través del espesor (cuando el conjunto sometido a tracción
se desplaza 0.5 mm) es la que se observa en la Figura 4-43.
Espesor (mm)
5
4
3
Orificio Nº 1
2
Orificio Nº 2
1
0
0
-500
-1000
-1500
-2000
Tensión [MPa]
Figura 4-43: Distribución de tensiones del modelo numérico con dos bulones transversales
Los resultados graficados en la Figura 4-43 quedan sobrepuestos, lo que muestra que
las tensiones son iguales en ambos orificios. Además se ve que el pico de tensiones
897,5 MPa en la superficie de contacto entre probetas, es levemente menor (del orden
del 10 %) al caso de un solo bulón centrado.
Por otra parte la distribución de tensiones a través del espesor en la Figura 4-44,
corresponden al modelo laminado por capas. En este modelo se observa que las capas se
descargan con respecto al caso de un bulón centrado (ver Figuras 4-25 hasta 4-28),
disminuyendo las tensiones para el caso de las láminas orientadas axialmente (0°) en el
62
orden del 35 %. Sin embargo ambos orificios tienen la misma magnitud de tensiones a
diferencia del caso de dos bulones longitudinales donde uno de los orificios tiene una
mayor descarga de tensiones.
5
Orificio 1
Espesor (mm)
Espesor (mm)
5
4
45º
3
0º
2
0
90º
0
-500
4
45º
3
0º
2
-45º
1
Orificio 2
-45º
90º
1
0
-1000 -1500 -2000 -2500 -3000
0
-500
-1000 -1500 -2000 -2500 -3000
Tensión [MPa]
Tensión [MPa]
Figura 4-44: Distribución de tensiones por capas a través del espesor en el orificio Nº 1 y Nº 2
para la unión con dos bulones transversales
4.5.2.4. Distorsión del plano medio de la placa para dos bulones transversales
Se emplea nuevamente el sistema de orientación que usó Ekh [16] para representar
esta distorsión (ver Figura 4-19). Los resultados del desplazamiento del plano medio se
muestran en la figura siguiente.
Desplazamiento Z (mm)
0.8
0.7
Modelo de dos bulones transversales
0.6
0.5
0.4
Lado 1
0.3
Lado 2
0.2
0.1
0
-0.1
-0.2
-0.3
0
10
20
30
40
50
60
70
80
Distancia sobre el plano medio (mm)
Figura 4-45: Desplazamiento del plano medio en la dirección z para el modelo de dos bulones transversales
63
Los resultados presentados en el Figura 4-45 muestran la distorsión típica de torsión
en la placa, que también había sido detectada en los ensayos de tracción. Al igual que
sucede en las probetas con un solo bulón y dos bulones axiales. Sin embargo se puede
notar que en este caso, el desplazamiento relativo entre los lados es significativamente
mayor que lo observado en los otros casos estudiados en este trabajo. Si bien esta
distorsión no afecta los parámetros mecánicos de la unión, si puede suceder que la
misma afecte el normal comportamiento de la unión entre partes estructurales.
4.5.2.5. Deformada final de la unión con dos bulones transversales
De acuerdo con los resultados del análisis numérico, la deformada final para una
carga de 5 kN es la que se observa en la Figura 4-46. Se puede ver que la distorsión por
flexión es similar al caso de un solo bulón centrado, lo cual permite inferir que la
interacción entre los bulones y la placa será similar al caso mencionado.
Figura 4-46: Deformada final de la unión con dos bulones transversales
64
CAPÍTULO
5
5
CONCLUSIONES
5.1 Introducción
En este trabajo se han simulado computacionalmente diferentes configuraciones para
uniones abulonadas de elementos fabricados con materiales compuestos. La validación
del modelo numérico propuesto se realizó a partir de los ensayos experimentales del
Proyecto BOJCAS [4], que fue llevado a cabo en la Unión Europea. Además se han
replicado los trabajos numéricos realizados por McCarthy y colab. [5], como parte del
desarrollo de un modelo numérico eficiente y robusto. De esta manera se ha obtenido y
validado un modelo adecuado para el análisis de uniones de materiales compuestos de
solape simple vinculados por medio de bulones.
Las simulaciones requirieron estudiar los desarrollos realizados dentro de los trabajos
mencionados e identificar:
• Los materiales utilizados.
• Las configuraciones del compuesto analizado.
• Las configuraciones de abulonado seleccionadas.
• El método de ensayo mecánico aplicado.
• Los resultados obtenidos experimentalmente.
• Los resultados obtenidos numéricamente.
• Las conclusiones a que se arribó.
En este trabajo se han expuesto los principales tópicos referidos al comportamiento
mecánico de materiales compuestos que se relacionan directamente con el análisis realizado.
Se han identificado las ecuaciones que modelan la relación tensión-deformación, así
como la mejor metodología general de análisis para la probeta abulonada a través del
método de elementos finitos.
Se enfrentó la complejidad de modelar el problema tridimensional, con el fin de
poder captar fielmente las interacciones debido al contacto existente entre los
componentes involucrados en el modelo. Para ello se realizó un análisis de las
características de las placas bajo el concepto de sublaminado, obteniéndose resultados
satisfactorios con el modelo propuesto en el presente trabajo que fueron contrastados
con resultados de trabajos experimentales y numéricos de otros autores.
65
La buena concordancia observada al comparar los resultados de este trabajo con los
resultados experimentales y numéricos reportados por BOJCAS [4], McCarthy y colab.
[5], Anderson [22] y Ekh [16], permiten concluir que el modelo elaborado representa de
una manera realista el comportamiento mecánico de las uniones abulonadas de materiales
compuestos. Esa concordancia da un respaldo adecuado a la modelación numérica de
otras configuraciones de abulonado y solapamiento, presentadas en esta tesis, basada en el
modelo desarrollado en este trabajo para un abulonado simple.
5.2 Conclusiones
El análisis de las uniones abulonadas de placas de materiales compuestos realizado
en este trabajo arrojó resultados convincentes, lo cual sugiere que es posible la
simulación computacional de estas piezas que tienen un comportamiento complejo. En
líneas generales las conclusiones de este trabajo son:
• La simulación computacional por elementos finitos del abulonado de placas de
materiales compuestos, proporcionó resultados concordantes con los que se
obtuvieron en los ensayos experimentales del proyecto BOJCAS. En particular en
lo referente al comportamiento mecánico del conjunto, es decir se obtuvieron
desplazamientos y deformaciones similares a los medidos experimentalmente y
plena concordancia en el valor de la rigidez de la unión.
• La distribución de tensiones en el espesor, sobre la línea de bulones, obtenida
numéricamente en este trabajo presenta similitudes importantes con los
resultados numéricos reportados en la literatura para casos similares. Igual
situación se repite con los desplazamientos del plano medio de las placas de
compuesto.
• La selección de la aproximación por elementos finitos, empleando en el mallado
de la geometría elementos tridimensionales, resulta en un todo de acuerdo con la
naturaleza del problema. Los elementos que mejor se adaptan y captan de una
forma eficiente los cambios bruscos en la distribución espacial de las tensiones,
han sido los hexaedros cuadráticos de 20 nodos.
• La extensión del modelo, desarrollado para el análisis la unión simple de un bulón,
al estudio del problema de dos bulones con diferentes configuraciones, mostró un
comportamiento cuantitativamente equivalente a los resultados experimentales
obtenidos por Morales y Quadri [6]. Además se reprodujeron adecuadamente
los estados tensionales que producen el aplastamiento de las placas de material
compuesto, lo que posteriormente provoca el daño observado por dichos autores.
66
• El estudio de configuraciones más generales de abulonado, muestra que las
configuraciones donde los elementos de unión están dispuestos en tándem (en la
dirección de carga) presentan mejores resultados. Particularmente, la ventaja de
esas configuraciones está asociada a la mejora de los parámetros principales de
la unión: rigidez, deformaciones globales y tensiones de aplastamiento en los
orificios.
• Los resultados obtenidos permiten concluir que el método de elementos finitos es
una potente herramienta para simular condiciones reales de trabajo en uniones
abulonadas de materiales compuestos. Se ha comprobado además que los
resultados obtenidos son bastante fiables, esto puede aprovecharse a los fines de
ahorrar tiempo y recursos en procesos industriales durante la fase de diseño de
uniones abulonadas.
5.3 Contribuciones del presente trabajo
A continuación se enumeran las principales contribuciones de este trabajo:
Se ha desarrollado una metodología para replicar situaciones reales de vinculación de
piezas de materiales compuestos aplicando software comercial de elementos finitos. Lo
importante de esta metodología, es que se han establecido los parámetros adecuados
para que los resultados de la simulación sean acordes a la realidad del problema.
La complejidad del problema estudiado y los resultados precisos obtenidos, muestran
que es posible emplear esta metodología en forma alternativa y/o complementaria a los
ensayos experimentales a realizarse sobre uniones abulonadas de partes estructurales de
compuesto. En función de esto se puede afirmar que los detalles de modelación que se
presentan en esta tesis son un aporte importante para la industria.
La extensión del modelo propuesto a configuraciones más generales, a partir de la
validación de la metodología de cálculo sobre un esquema simple de unión, ha permitido
establecer que configuración de unión es la más adecuada.
Por último y a nivel personal, este trabajo brinda una excelente oportunidad para
establecer nuevas líneas de investigación en el campo de los compuestos dentro del
Área de Ciencias Aeronáuticas del Departamento de Seguridad y Defensa (Universidad
de las Fuerzas Armadas) en Ecuador y en el Centro de Investigación y Desarrollo de la
Fuerza Aérea Ecuatoriana.
67
5.4
Líneas de trabajos futuros
A partir del análisis numérico de una unión mediante bulones para placas de material
compuestos presentado en estas tesis, es posible delinear algunas líneas futuras de trabajo
como por ejemplo:
•
Reemplazar la representación realizada a través de sólidos tridimensionales,
desarrollada en el presente trabajo, por un modelo de lámina equivalente (con
elementos de lámina planos), a fin de comparar los resultados obtenidos y
verificar las diferencias que provee cada modelo.
•
Extender el análisis numérico, desarrollado en el Proyecto BOJCAS y en esta
tesis, a la unión de materiales compuestos con otras configuraciones de
apilamiento, otros tipos de solapamiento y diferentes tipos de elementos de
vinculación (remaches, pegado, etc.)
•
Realizar trabajos simultáneos de experimentación y modelación numérica
para el análisis del comportamiento mecánico de materiales compuestos
empleando para ello distintos tipos de ensayos mecánicos, que simulen
condiciones reales de operación de este tipo de materiales y comparar los
resultados obtenidos con las simulaciones computacionales. Esto permitirá
sin duda hacer aportes significativos a la industria para la solución de
problemas reales en el diseño de estructuras de materiales compuestos.
•
Analizar la influencia del torque aplicado por el bulón sobre las placas de
material compuesto.
•
Investigar cómo influye la interferencia (control de huelgos) en los orificios
de abulonado, en las uniones abulonadas de materiales compuestos.
68
ANEXO
MATRICES DE RIGIDEZ DE LAMINADOS
Y SUBLAMINADOS
1 Introducción
En este anexo se presentan las consideraciones pertinentes a la mecánica de materiales
compuestos, que son un compendio de los desarrollos realizados por varios autores, en
particular Reddy [1] y Kóllar [10], que han estudiado el comportamiento y la formulación
matemática de los materiales compuestos en profundidad.
2. Matrices de transformación tensión-deformación
Consideramos dos sistemas coordenados cartesianos: el p-q-r al que en adelante
llamaremos “original” y el p'-q'-r' al que trataremos como “primo”. Esos dos sistemas
tienen como punto de coincidencia invariable su origen pero el primo está rotado con
respecto al original. La relación entre la orientación del sistema p'-q'-r', con la del
sistema referencial p-q-r, puede expresarse en base a tres ángulos: θ p , θ q y θ r que son
las rotaciones del sistema coordenado primo como se muestra en la Figura 1.
Figura 1: Rotaciones consecutivas para alcanzar la posición final de p'-q'-r'
69
Los ángulos θ p , θ q y θ r son respectivamente las rotaciones alrededor de los ejes p,
q y r, y para definirlos se debe recurrir a los cosenos directores (r 11 , r 21 , r 31 ), (r 12 , r 22 ,
r 32 ) y
(r 13 , r 23 , r 33 ), que dan la orientación del sistema coordenado primo respecto al original y
que se definen del modo presentado en la Figura 2.
=
r11 cos
=
Ω pp
r12 cos
=
Ω pq
r13 cos Ω pr
r21 cos
Ωqp
r22 cos
Ωqq
r23 cos Ωq r
=
=
=
(1)
=
r31 cos
=
Ωrp
r32 cos
=
Ωrq
r33 cos Ωrr
Figura 2: Definición de los cosenos directores
La relación de los ángulos θ p , θ q y θ r , con los cosenos directores está dada por las
expresiones:
1/2
θq =
f  −r31 , ( r112 + r212 ) 


 r
r 
θr =
f  21 , 11 
 cos θ q cos θ q 
 r
r 
θp =
f  32 , 33 
 cos θ q cos θ q 
(2)
donde f es una función de dos argumentos, que está definida de la siguiente manera:
 tan −1 ( x / y ) cuando y > 0

f ( x, y ) = 
 tan −1 ( x / y ) + π cuando y < 0
Cuando θ q = 90o → θ q = f (r 12 , r 22 ) y θ r = 0.
70
(3)
2.1 Transformación de tensiones
Las tensiones en el sistema coordenado primo, son calculadas a partir de las
tensiones en el sistema coordenado original, empleando para ello la transformación:
σ ′p   Tσ 11
 ′ 
σ q   Tσ 21
 σ r′   Tσ 31
τ ′  =  T
 q r   σ 41
τ ′pr   Tσ 51
  T
τ ′pq   σ 61
Tσ 12
Tσ 22
Tσ 32
Tσ 42
Tσ 52
Tσ 62
Tσ 13 Tσ 14
Tσ 23 Tσ 24
Tσ 33 Tσ 34
Tσ 43 Tσ 44
Tσ 53 Tσ 54
Tσ 63 Tσ 64
Tσ 15
Tσ 25
Tσ 35
Tσ 45
Tσ 55
Tσ 65
Tσ 16
Tσ 26
Tσ 36
Tσ 46
Tσ 56
Tσ 66









σ p 
 
σ q 
 σ r 
τ  → σ′ = Tσ σ
 qr 
τ pr 
 
τ pq 
(4)
En la ecuación (4), la matriz de transformación Tσ puede ser expresada como:
Tσ = Tσp Tσq Tσr
(5)
siendo:
1
0

0
Tσp = 
0
0

0
0
cos 2 θ p
0
sen 2 θ p
2sen θ p cos θ p
0
0
2
sen θ p
2
cos θ p
−2sen θ p cos θ p
0
−sen θ p cos θ p
0
0
senθ p cos θ p
0
0
cos θ p − sen θ p
0
0
0
cos θ p
sin θ p
 cos 2 θ q

0

 sen 2 θ q
q
Tσ = 
0

 − sen θ q cos θ q

0

 cos 2 θ r

2
 sen θ r

0
Tσr = 
0


0

 − senθ r cos θ r
0
2
2
0
1
0
0
sen 2 θ q
0
cos 2 θ q
0
0
0
0
cos θ q
2sen θ q cos θ q
0
−2sen θ q cos θ q
0
0
0
sen θ q cos θ q
0
0
senθ q
cos 2 θ q − sen 2 θ q
0
sen 2 θ r
cos 2 θ r
0
0
0
0
0
0
0
0
0
1
0
0
0
cos θ r
sen θ r
0
−sen θ r
cos θ r
sen θ r cos θ r
0
0
0



0 
 (6)
0 
− sen θ p 

cos θ p 
0
0
0 

0 
0 
(7)
− sen θ q 
0 

cos θ q 
2 sen θ r cos θ r 

−2sen θ r cos θ r 

0
 (8)
0


0

cos 2 θ r − sen 2 θ r 
Dado que en el presente trabajo estamos interesados en una rotación del sistema
primo alrededor del eje r, donde se mantienen coincidentes r y r', se establece que los
71
ángulos θ p y θ q son nulos. De esta forma Tσp y Tσq se convierten en matrices identidad,
y la matriz Tσr resulta igual a Tσ .
En el caso de sistemas cartesianos bidimensionales p-q y p'-q', la relación entre las
tensiones en el sistema primo y el sistema original esta dado por:
σ ' p 


=
σ 'q 
τ ' pq 


 cos 2 θ r
 sen 2 θ
r

sen
θ
cos
θr
−
r

sen 2 θ r
cos 2 θ r
senθ r cos θ r
2 senθ r cos θ r 
− 2senθ r cos θ r 

cos 2 θ r − sen 2θ r 
σ p 
 
σ q 
 
τ pq 
(9)
2.2 Transformación de deformaciones
Similarmente a lo mostrado previamente, la deformación en el sistema coordenado
primo con respecto al sistema coordenado original está dado por:
 ε ' p   Tε 11 Tε 12
ε '  T
 q   ε 21 Tε 22
 ε 'r   Tε 31 Tε 32
γ '  =  T
Tε 42
 qr   ε 41
γ ' pr   Tε 51 Tε 52
γ ' pq   Tε 61 Tε 62

 
Tε 13
Tε 23
Tε 33
Tε 43
Tε 53
Tε 63
Tε 14 Tε 15 Tε 16
Tε 24 Tε 25 Tε 26
Tε 34 Tε 35 Tε 36
Tε 44 Tε 45 Tε 46
Tε 54 Tε 55 Tε 56
Tε 64 Tε 65 Tε 66
 εp 
 ε 
 q 
  ε r  → ε′ = T ε
ε
  γ qr 


 γ
  pr 
 
 γ pq 
(10)
En la ecuación (10) el vector ε representa las deformaciones ingenieriles. La matriz
Tε que se aplica a este tipo de transformaciones, que es distinta de Tσ , puede ser
expresada como:
Tε = Tεp Tεq Tεr
(11)
siendo estas matrices:




p
Tε = 




1
0
0
0
0
0
0
cos 2 θ p
sen 2 θ p
− 2sen θ p cos θ p
0
0

cos 2 θ q

0

sen 2 θ q

Tεq = 
0

 −2sen θ q cos θ q

0

0
1
0
0
0
0
0
0
0
0
cos θ p
sen θ p
0 
0 

0 
(12)
0 

− sen θ p 
cos θ p 
sen θ q cos θ q
0
− sen θ q cos θ q
0
2
cos θ q − sen 2 θ q
0
0 
0 

0 
(13)
− senθ q 

0 
cos θ q 
0
sen 2 θ p
cos 2 θ p
2senθ p cos θ p
0
0
0
sen θ p cos θ p
− sen θ p cos θ p
cos 2 θ p − sen 2 θ p
0
0
sen 2 θ q
0
cos 2 θ q
0
2sen θ q cos θ q
0
0
0
0
cos θ q
0
sen θ q
72

cos 2 θ r

sen 2 θ r

0
Tεr = 
0


0

 − 2sen θ r cos θ r
sen 2 θ r
cos 2 θ r
0
0
0
2sen θ r cos θ r
0
0
1
0
0
0
0
0
0
cos θ r
sen θ r
0
0
0
0
−sen θ r
cos θ r
0
sen θ r cos θ r 
− sen θ r cos θ r 

0
 (14)
0


0

2
2
cos θ r − sen θ r 
Al igual que en la matriz de transformación de tensiones, el giro alrededor del eje r
hace que las matrices Tεp y Tεq sean iguales a la matriz identidad, y Tεr sea igual a Tε .
Este hecho al ser aplicado a un caso bidimensional da lugar a:
 ε 'p 


 ε 'q 


γ ' pq 

cos 2 θ r

sen 2 θ r

 −2sen θ r cos θ r

sen 2 θ r
cos 2 θ r
2sen θ r cos θ r
sen θ r cos θ r
− sen θ r cos θ r
cos 2 θ r − sen 2 θ r





εp 
 
 εq 
 
γ pq 
(15)
3 Condición de tensión plana
Una de las principales suposiciones para el análisis mecánico de placas laminadas de
materiales compuestos, es la condición de tensión plana bajo la cual se asume que las
tensiones normales y las tensiones de corte fuera del plano son nulas. Considerando que
la dirección z es normal a la placa de compuesto, se tiene:
σz = 0
τ xz = 0
τ yz = 0
(16)
La condición dada en (16) aproxima adecuadamente las tensiones de una placa
delgada de materiales compuestos laminados cuando las fibras son paralelas al plano x-y
y la placa se carga a lo largo de sus extremos. De este modo las fuerzas resultan
paralelas a dicho plano y se distribuyen uniformemente en el espesor. Si bien esta
suposición no provee las tensiones exactas en la placa, sirve para obtener valores con
una precisión aceptable. En virtud de esto, la relación deformación-tensión, mostrada en
la ecuación (5) del Capítulo 2 referida al sistema coordenado x-y-z es:
 ε x   S11
  
ε y  = 
γ   sim
 xy  
S12
S 22
S16 
S 26 
S36 
σ x 
 
σ y 
τ 
 xy 
Invirtiendo la ecuación (17) se obtiene la matriz de rigidez Q que está dada por:
73
(17)
σ x   Q11
  
σ y  = 
τ   sim
 xy  
Q12
Q22
Q16 
Q26 
Q66 
εx 
 
ε y 
γ 
 xy 
(18)
Esta relación se adapta a la condición que presenta el material, y de acuerdo a la
orientación de sus fibras puede ser tratado como monoclínico, ortótropo, transversamente isótropo o completamente isótropo. En esta tesis se trabaja con un material
ortótropo, lo cual conlleva a que los elementos S 16 y S 26 de la matriz de flexibilidad S
sean nulos. Entonces la matriz de flexibilidad tiene la siguiente forma:
−ν xy / Ex
 1/ Ex

S = 

 sim
1/ E y


0 

1 / Gxy 
0
(19)
Luego, expresando los elementos de la matriz de rigidez Q en términos de las
constantes ingenieriles se tiene:
−ν xy Ex / D
 Ex / D

Q= 

 sim
Ey / D
0 

0 

Gxy 
(20)
Ey 2
1−
ν xy =
1 −ν xyν yx
donde D =
Ex
Con la ecuación (20) se puede obtener la matriz de rigidez, para una capa de material
laminado de fibras unidireccionales que están alineadas con el eje x. A partir de esta, y
mediante el uso de las matrices de transformación, se puede determinar la matriz de
rigidez para otras capas cuya orientación de fibras no esté alineada con el eje x.
Las relaciones tensión-deformación en los sistemas coordenados original y primo
definidos anteriormente, están dadas por:
σ = Qε
σ ′ = Q′ ε ′
(21)
Multiplicando a la primera de las relaciones descritas en (21) por la matriz de
transformación de tensiones Tσ resulta:
Tσ σ = Tσ Q ε

σ′
74
(22)
Para que en la ecuación (22) el vector de tensiones en el sistema coordenado primo
σ′ esté relacionado con el vector de deformaciones ε′ en el mismo sistema, es necesario
introducir la relación descrita en (10). Entonces se multiplica por ( Tε −1 Tε ), que resulta
en la matriz identidad y permite llegar al objetivo buscado:
σ′ = Tσ Q Tε−1 Tε ε


 
(23)
ε′
Q′
A partir de esta relación se deduce que la matriz de rigidez en el sistema coordenado
primo está dada por:
Q′ = Tσ Q Tε−1
(24)
Si tenemos un compuesto laminado de varias capas y un sistema de coordenadas
global x-y-z, cada capa tiene un sistema cartesiano propio x'-y'-z' donde sus fibras están
alineadas en la dirección longitudinal del eje x'. De modo que para referirlas al sistema
global del compuesto deberemos usar la ecuación (24) invertida de modo que:
Q = Tσ−1 Q′ Tε
(25)
En el caso de un material compuesto laminado conformado por varias capas del
mismo material, la matriz Q′ es igual en todos las capas y corresponde a la inclinación
de 0º de las fibras respecto al eje x'. Mientras que la matriz Q tiene en cuenta la
orientación de las fibras mediante la asociación de su ángulo de inclinación, algo que
está incluido en las matrices de transformación.
Para ejemplificar, se considera un compuesto de fibra de carbono unidireccional
HTA/6376 formado por una capa orientada a 0º y una a 45º con las propiedades
mecánicas descritas en la Tabla 2-1 del Capítulo 2. La matriz de flexibilidad para la
capa orientada a 0º, es decir en el sistema asociado a la dirección longitudinal de la fibra
resulta:
 1
0,3
−
 140 140

1

S = 
10

 sim


0 

m2

0  x1
0− 9
N

1 
5, 2 
(26)
Invirtiendo esta relación se tiene:
140,91

Q′ = 

 sim
3, 02
10, 07
75
0 

N
0  x 1 0 9 2
m

5, 20 
(27)
Utilizando la ecuación (25) la fibra alineada a 0º llevada al sistema global tiene una
matriz de transformación:
−1
 cos 2 0
sen 2 0
2sen 0 cos 0 


2
cos 2 0
0 cos 0 
− 2sen=
 sen 0

2
2 
 − sen 0 cos 0 sen 0 cos 0 cos 0 − sen 0 
=
Tσ−10o
=
Tε 0o
=
Q
0
1

0
0


cos 2 0
sen 2 0
sen 0 cos 0 


− sen 0 cos 0 
sen 2 0
cos 2 0 =


2
2 
 − 2sen 0 cos 0 2sen 0 cos 0 cos 0 − sen 0 
 1
 0

 0
−1
0 

0=
Q′

1 
0
1
0
 1
 0

 0
0
1
0
0
1
0
1

0
0

0

0
1 
0
1
0
−1
(28)
0

0
1 
(29)
0 
0  =
⇒ Q 0 Q′
1 
(30)
Para la fibra alineada a 45º se tiene:
−1
Tσ-145o
 cos 2 45
sen 2 45
2sen 45 cos 45 


2
=
−2sen 45 cos 45 
cos 2 45
 sen 45


2
2
 − sen 45 cos 45 sen 45 cos 45 cos 45 − sen 45
Tε 45o

cos 2 45
sen 2 45
sen 45 cos 45   0,5 0,5 0,5 

 

=
−sen 45 cos 45  =
sen 2 45
cos 2 45

 0,5 0,5 −0,5 

 
2
2
0 
 −2sen 45 cos 45 2sen 45 cos 45 cos 45 − sen 45  −1 1
Q 45o
0 
 0,5 0,5 1  140,91 3, 02
 0,5 0,5 0,5 

 



N
x1
=
0 9 2  0,5 0,5 − 0,5
10, 07
0   0,5 0,5 −1 
m 
 − 0,5 0,5 0   sim
5, 20 
0 

 
 −1 1
 0,5 0,5 1 


=
 0,5 0,5 −1
 − 0, 5 0,5 0 


−1
(31)
(32)
−1
Q 45o
 44, 45

= 
 sim

34, 05
44, 45
32, 71 

N
32, 71  x 1 09 2
m
36, 23 
(33)
(34)
A partir de las ecuaciones mostradas previamente es posible determinar las matrices
de rigidez A, B y D de un compuesto laminado, como se muestra en el Capítulo 2 del
presente trabajo en las ecuaciones (25) y (26).
76
4 Obtención de la matriz de flexibilidad de sublaminados
Los sublaminados son divisiones que se realizan al laminado, para facilitar el análisis
por elementos finitos usando mallas de elementos tridimensionales. Para esto, se emplea
la relación entre las deformaciones y las tensiones según la ecuación (5) del Capítulo 2
que debe ser cumplida por el sublaminado. Trasladando esta expresión a un sistema
coordenado x-y-z, utilizando S para referirse a la matriz de flexibilidad del sublaminado,
y denotando deformaciones y tensiones promedios mediante ε , γ , σ y τ tenemos:
 εx
ε
 y
 ε z

 γ yz
 γ xz

 γ xy
  S11
 
  S 21
  S31
= 
  S 41
  S51
 
  S61
S12 S13
S 22 S 23
S32 S33
S 42 S 43
S52 S53
S62 S63
S14
S 24
S34
S 44
S54
S64
S15
S 25
S35
S 45
S55
S65
S16
S 26
S36
S 46
S56
S66









σx
σ
 y
 σ z

 τ yz
 τ xz

 τ xy









(35)
En el presente desarrollo hay que tener en cuenta que las relaciones (3) y (4) del Capítulo 3,
mostradas a continuación, que se refieren a las tensiones y deformaciones promedio:
1
1
=
σ x dz
Nx =
σz σz
∫
hs hs
hs
1
1
=
σy =
σ y dz
Ny =
τ xz τ xz
∫
h
hs s
hs
1
1
=
τ xy =
τ xy dz
N xy =
τ yz τ yz
∫
hs hs
hs
(36)
1
dz
=
ε z εx εx
hs ∫hs
1
dz
=
γ yz =
γ yz εy εy
hs ∫hs
1
dz
γ xz =
γ xz γ xy = γ xy
hs ∫hs
(37)
=
σx
=
εz
4.1 Elementos de la matriz S debido a tensiones en el plano
Para obtener los elementos de la matriz S es necesario imponer la condición de que
las tensiones fuera del plano σ z , τ xz ,τ yz sean nulas como se muestra en la Figura 3. De
este modo la ecuación (35) se puede escribir como se muestra en (38), (39) y (40).
77
Figura 3: Tensiones por capas y tensiones promedio en la condición dada
 ε x   S11 S12
  
 ε y  =  S 21 S 22
γ   S
 xy   61 S62
εz
=
[
S31
S32
γ yz   S 41 S 42
 = 
γ xz   S51 S52
S16 

S 26 
S66 

σ x 
 
σ y 
τ 
 xy 
(38)
]
σ x 
 
σ y 
τ 
 xy 
(39)
σ 
S 46   x 
 σ 
S56   y 
τ xy 
(40)
S36
Basados en las ecuaciones (21) y (36) de este anexo, y con las condiciones previas,
las curvaturas son nulas por lo que se tiene:
 Nx 


=
 Ny 
N 
 xy 
 ε xo 
 o 
A
=
 εy 
 o 
 γ xy 
 εx 


A  εy 


 γ xy 
(41)
donde A es la matriz de rigidez del sublaminado. Tomando las ecuaciones (36) y (37)
resulta:
 Nx 


=
 Ny 
N 
 xy 
 σx 


h=
s  σy 
τ 
 xy 
 εx 
 εx 




A
=
 ε y  A  εy 




 γ xy 
 γ xy 
de donde:
78
(42)
 εx 
 
=
 εy 
 
 γ xy 
σ x 
 
−1
h=
σ y 
s A
 
τ xy 
 S11

 S 21
S
 61
S12
S 22
S62
S16   σ x 
 
S 26   σ y 
 
S66  τ xy 
(43)
Bajo la condición de tensión plana, la deformación transversal puede escribirse según
lo descrito en la ecuación (5) del Capítulo 2:
ε z = [ S13 S23 S36 ]
σ x 
 
σ y 
τ 
 xy 
(44)
donde σ x , σ y ,τ xy son las tensiones en las capas como lo muestra la Figura 3.
Recordando la relación tensión-deformación para tensión plana presentada en la
ecuación (21) de Capítulo 2, se tiene:
σ x 
εx 
 
 
σ y  = Q  ε y 
τ xy 
γ xy 
 
 
(45)
Utilizando las ecuaciones previas (44) y (45), junto a (37), la ecuación (44) se
convierte en:
ε z = [ S13 S23 S36 ]
σ x 
 
Q hs A σ y 
τ 
 xy 
−1
(46)
Utilizando nuevamente las definiciones dadas en (37) y reemplazando la integral con
la sumatoria extendida al número de capas, M, resulta:
=
εz
M
∑ [S
k =1
13
S 23
σ x 
 
S36 ] ( zk − zk −1 ) Q k A −1 σ y 
τ 
 xy 
ε z = [ S31 S32
σ x 
 
S36 ] σ y 
τ 
 xy 
(47)
(48)
En la condición analizada en el presente caso, la ecuación (40) puede escribirse como:
 0   S 41 S 42
 0  = S
   51 S52
Esta condición se satisface sólo si:
79
σ 
S 46   x 
σ
S56   y 
τ xy 
(49)
 S41 S 42
 S51 S52
S 46   0
=
S56   0
0
0
0
0 
(50)
4.2 Elementos de la matriz S debido a las tensiones fuera del plano
Aquí se considera que en el laminado se tienen únicamente tensiones transversales,
lo cual se consigue considerando:
a. Un sublaminado restringido en sus extremos de modo que las deformaciones
ε x , ε y , y γ xy sean nulas, al que se le aplica una tensión uniforme σ z = σ z que
genera las tensiones σ x , σ y , y τ xy en el plano.
b. Un sublaminado al que se aplica una tensión plana − σ x , − σ y , y − τ xy .
c. Se superponen las condiciones a) y b) con lo que se logra un laminado en el
que σ x , σ y , τ xz , τ yz y τ xy son nulas.
(a)
(b)
(c)
Figura 4: (a) Laminado restringido en sus extremos, (b) Laminado sujeto a una tensión
y (c) Sublaminado no restringido sujeto a un tensión σ z
σx
Las relaciones tensión-deformación para la condición ‘a.’ se convierten en:
εx 
 
ε y  =
γ 
 xy 
 S13 
S  σ
 23  z
 S63 
(51)
ε z = S33 σ z
(52)
γ yz 
 =
γ xz 
 S 43 
S  σ z
 53 
(53)
Aplicando las ecuaciones anteriores en la ecuación (3) del Capítulo 2, se obtiene:
σz
=
C33 ε z
σ x   C13 
   
σ y  = C23  ε z
τ  C 
 xy   63 
de donde:
80
(54)
εz
1
σz
C33
=
σ x   C13 
    1
σz
σ y  = C23 
τ  C  C33
 xy   63 
(55)
Combinando estas ecuaciones en las correspondientes de (36) y (37), y reemplazando
la integral por una sumatoria:
εz =
1
hs
1
=
hs
σ x 
 
1
σ y  =
hs
τ 
xy
 
M
∑
k =1
∫
hs
1
σ z ds
C33
M
zk − zk −1
∑ (C )
k =1
(56)
σz
33 k
 C13 
  zk − zk −1
C23  ( C ) σ z
33 k
C 
 63 
(57)
Trabajando sobre la condición ‘b.’ usando las ecuaciones (43), (47) y (57) se llega a:
σ x 
 
1
εz =
− [ S31 S32 S36 ] σ y  =
− [ S31 S32 S36 ]
hs
τ 
xy
 
M
∑
k =1
 εx 
 S11 S12 S16  σ x 
 S11 S12 S16 
 

 

1 M
−  S 21 S 22 S 26  σ y  =
−  S 21 S 22 S 26  ∑
εy  =
h
γ 
 S S S  τ 
 S S S  s k =1
 61 62 66   xy 
 61 62 66 
 xy 
 C13 
  zk − zk −1
C23  ( C ) σ z
33 k
C 
 63 
(58)
 C13 
  zk − zk −1
C23  ( C ) σ z
33 k
C 
 63 
(59)
Con la condición ‘c.’ de superposición, operando sobre las ecuaciones (56), (58) y
(59) se obtiene:


 C13 
 1 M zk − z K −1
1 M   zk − zk −1 
εz  ∑
=
− [ S31 S32 S36 ] ∑ C23 
σ z
h
C
h
C
(
)
(
)
1
1
k
k
=
=
33
33
s
s

 
k
k 
C63 






 J13 
 
 J 23 
 
 J 63 
81
(60)
 εx 
 S11 S12 S16 
 C13 
 

 1 M   zk − zk −1
σz
 ε y  = −  S 21 S 22 S 26  ∑ C23 
 

 hs k =1   ( C33 )k
 S61 S62 S66 
C63 
γ xy  
(61)
J 33
Debido a que τ xz y τ yz son nulas, y a la condición de ortotropía mostrada en la matriz
de rigidez en la ecuación (6) del Capítulo 2, las deformaciones γ xz y γ yz resultan nulas,
y con la definición de (37), se tiene:
0 
 =
0 
 S 43 
  σz
 S53 
(62)
esta ecuación se satisface únicamente si:
 S 43   0 
 =  
 S53   0 
(63)
4.3 Elementos de S debido a tensiones cortantes fuera del plano
Para este desarrollo se aplican las tensiones cortantes τ xz y τ yz en el sublaminado. Las
tensiones σ x , σ y , σ z y τ xy son nulas, de modo que las relaciones tensión-deformación
son:
γ yz   S 44 S 45  γ yz 
(64)
 = 
  
γ xz   S54 S55  γ xz 
εx 
 
ε y 
 =
εz 
γ 
 xy 
 S14

 S 24
S
 34
 S64
S15 

S 25 
S35 

S65 
τ yz 
 
τ xz 
(65)
S 45  τ yz 
 
S55  τ xz 
(66)
Para una sola capa tenemos:
γ yz   S 44
 = 
γ xz   S54
Combinando las ecuaciones (37) y (66) y reemplazando las integrales por sumatorias
el resultado es:
γ yz 
=
 
γ xz 
M 
 S 44
1
∑  ( zk − zk −1 ) 
hs k =1 
 S54

82
S 45    τ yz 
 

S55  k   τ xz 

(67)
lo cual permite determinar los elementos correspondientes de S por inspección simple
de (64) y (67):
 S 44 S 45 
=

 S54 S55 
M 
S
1
∑  ( zk − zk −1 )  44
hs k =1 
 S54
S 45    τ yz 


S55  k   τ xz 

(68)
Para sublaminados sometidos a tensiones que provocan únicamente deformaciones
γ yz y γ xz con todas las demás deformaciones nulas, la relación (65) se torna:
0 
0 
0  =
 
0 
 S14
 S 24
S
 34
 S64
S15
S 25
S35
S65





 S14
 S 24
S
 34
 S64
S15
S 25
S35
S65

0

0
 = 0



0
τ yz 
 
τ xz 
(69)
0
0
0

0
(70)
de modo que:
4.4 Conformación de la matriz de rigidez
Combinando los resultados obtenidos en los tres pasos previos, se obtiene la relación
deformación-tensión siguiente:
εx 
 S11
 

ε y 
 S 21
ε 
S
 z
 31
=
 
 S 41
γ yz 

γ 
 S51
xz
 


 S61
γ xy 

S12
S13
S14
S15
S 22
S 23
S 24
S 25
S32
S33
S34
S35
S 42
S 43
S 44
S 45
S52
S53
S54
S55
S62
S63
S64
S65
S16 

S 26 
S36 

S 46 

S56 

S66 
σ x 
 
σ y 
σ 
 z
 
τ yz 
τ 
 xz 

τ xy 

(71)
La matriz de rigidez se obtiene invirtiendo la matriz de flexibilidad S. Para el caso de
la condición de tensión plana, la relación (71) se transforma en:
 ε x   S11 S12
  
 ε y  =  S21 S 22
  
γ xy   S61 S62
83
S16  σ x 
 
S 26  σ y 
S66  τ xy 
(72)
84
BIBLIOGRAFÍA
[1] Reedy, J.N., Mechanics of Laminated Composite Plates and Shells, Second
Edition, CRC Press, Boca Raton, USA, 2004.
[2] The Boeing Company, Composites in the airframe and primary structure, Boeing
787, From the Ground Up, Aeromagazine, 04, 2006.
[3] Getting started with Abaqus: Interactive Edition, Dassault Systèmes, Providence,
USA., 2010.
[4] McCarthy, M.A., BOJCAS: Bolted joints in composite aircraft structures. Air and
Space Europe, 3, 3/4, pp.139-142, 2001.
[5] McCarthy, M.A., McCarthy, C.T., Lawlor, V.P. and Stanley, W.F., Three-dimensional
finite element analysis of single-bolt, single-lap composite bolted joints,
Composite Structures, 71, pp.140-158, 2005.
[6] Morales, M. y Quadri, T., Caracterización de uniones mecánicas en materiales
compuestos, Proyecto de Grado de la FCEFyN, U.N.C., 2011.
[7] Paris, F., Cañas, J. y Marín, J.C., Introducción al Análisis y Diseño con Materiales
Compuestos, Primera Edición, Universidad de Sevilla, Sevilla, España, 2006.
[8] Car, E., Oller, S. y Oñate, E., Tratamiento Numérico de Materiales Compuestos,
Primera Edición, Centro Internacional de Métodos Numéricos en Ingeniería,
Barcelona, España, 2000.
[9] Beckwith, S.W., Quasi-Isotropic Composite Materials/Structures Design - Part I,
SAMPE Journal, 48, pp 48, 2012.
[10] Kóllar, L.P. and Springer, G.S., Mechanics of Composite Structures, First Edition,
Cambridge University Press, Cambridge, United Kingdom, 2003.
[11] Zienkiewicz, O., Zhu, J. and Taylor, R., The Finite Element Method, Vol. 1 y 2,
Sexta Edición, Elsevier, 2006.
[12] Barker, R.M., Lin, F.T. and Dana, J.R., Three Dimensional Finite-Element Analysis
of Laminated Composites, Computers & Structures, 2, pp. 1013-1029, 1972.
[13] Barbero, E., Finite Element Analysis of Composite Materials, First Edition, CRC
Press, Boca Raton, USA, 2007.
[14] Liu G. and Quek S., The Finite Element Method, A Practical Course, First Edition,
Butterworth-Heinemann, 2003.
[15] Oñate, E., Structural Analysis with the Finite Element Method, Vol. 1, First
Edition, Springer, Barcelona, España, 2009.
85
[16] Ekh, J., Three-dimensional Stress Analysis. BOJCAS Interim Report, 2001.
[17] Song, G., Composites modeling: analysis procedures and techniques, Conferencia
virtual, http://www.simulia.com/download/Webinars/Central_Composite_Series/, 2012.
[18] ASTM Standard. D5961/D5961M-13 Standard Test Method for Bearing Response
of Polymer Matrix Composite Laminates, 2013.
[19] Milligan, D., Abaqus Composite Tutorial: Creating a Bolt Preload in 3 Easy Steps,
Converging on Composites - The Autodesk Composites Blog: http://info.firehole.com,
2012.
[20] Dominique, F., André, P. and André, Z., Mechanical Behaviour of Materials:
Volume II: Fracture Mechanics and Damage. Springer, París, Francia, 2012.
[21] Milligan, D., Abaqus Composite Tutorial – Bolted Joint Contact Definition, Converging
on Composites - The Autodesk Composites Blog: http://info.firehole.com/, 2012.
[22] Anderson, B., A Splitting Method for Fast Solution of 3D Contact Problems in
Bolted Joints. BOJCAS interim report, 2001.
[23] Ireman, T., Three-dimensional stress analysis of bolted composite single lap joints.
Composite Structures, 43, pp.195-216, 1998.
[24] Ramkumar, R.L. and Tossavainen, E., Bolted Joints in Composite Structures:
Design, Analysis and Verification. Air Force Wright Aeronautical Laboratories,
Hawthorne, USA. 1985.
[25] Rowlands, R.E., Rahman, M.U., Wilkinson, T.L. and Chiang, Y.I., Single and multiple
bolted joints in orthotropic materials, Composites, 13, pp. 273-279, 1982.
[26] ASTM Standard. D3039/D3039M-08 Standard Test Method for Tensile Properties
of Polymer Matrix Composite Materials, 2008.
[27] Ekh, J., Multi-Fastener Single-lap Joints in Composite Structures, Tesis Doctoral
del Royal Institute of Technology, Estocolmo, Suecia, 2006.
[28] McCarthy, C.T. and McCarthy, M.A., Finite element analysis of the effects of
clearance on single-shear, composite bolted joints, journal of plastics, rubber and
composites, The Institute of Materials, 32, Nº 2, London, United Kingdom, 2003.
86

examen 2.2

modelo_examen_trigonometria

Resueltos - FerMates

Integrales inmediatas y casi inmediatas

Tarea 8 – CII (2015)

Dr. José Emiliano Martínez Ordaz jose.martinez

Anexo 3-C Campo Nombre de campo Descripción Tipo Decimales

Solución Tarea #2 Problema #1

Mathcad - Trabajo Practico 1

Ejercicios resueltos del primer examen parcial del curso 2014–2015