SMEC 2009 2do SIMPOSIO REGIONAL EN MECÁNICA DE MATERIALES Y ESTRUCTURAS CONTINUAS: EXPERIMENTACIÓN, MODELADO NUMÉRICO Y TEÓRICO

21, 22 y 23 de Octubre de 2009 Cartagena de Indias - Colombia

SMEC 2009 2do SIMPOSIO REGIONAL EN MECANICA DE MATERIALES Y ESTRUCTURAS CONTINUAS: EXPERIMENTACIÓN, MODELADO NUMÉRICO Y TEÓRICO MEMORIAS DEL EVENTO Volumen 1, Octubre de 2009 Institución editora Facultad de Ingeniería – Universidad Tecnológica de Bolívar Los conceptos y opiniones de los artículos contenidos en estas memorias son responsabilidad de su autor; En ningún momento comprometen las orientaciones y políticas de la Facultad de Ingeniería de la Universidad Tecnológica de Bolívar.

Contacto: Prof. Jairo F. Useche, Ph.D. Departamento de Ingeniería Mecánica y Mecatrónica Universidad Tecnológica de Bolivar Parque Industrial Velez-Pombo, km.1 Tel/fax: +575 6535337/6619240 [email protected] Cartagena, Colombia

SMEC 2009 2do SIMPOSIO REGIONAL EN MECANICA DE MATERIALES Y ESTRUCTURAS CONTINUAS: EXPERIMENTACIÓN, MODELADO NUMÉRICO Y TEÓRICO MEMORIAS DEL EVENTO

Dirección Universitaria Patricia del Pilar Martínez Barrios Rectora Paola Amar Sepulveda Vicerrectora Académica Viviana Londoño Moreno Vicerrectora Administrativa Javier Campillo Jimenez Director Oficina de Investigaciones y Transferencia Tecnológica

Facultad de Ingeniería Jose Luis Villa Decano Facultad de Ingeniería Raul Padrón Carvajal Secretario Academico Justo Ramos Madrid Director Departamento de Ingeniería Mecánica y Mecatrónica

SMEC 2009 2do SIMPOSIO REGIONAL EN MECANICA DE MATERIALES Y ESTRUCTURAS CONTINUAS: EXPERIMENTACIÓN, MODELADO NUMÉRICO Y TEÓRICO MEMORIAS DEL EVENTO

Cuerpo Editorial Jairo Useche Vivero

Editor General Comité Editorial Universidad Tecnológica de Bolívar Universidad Tecnológica de Bolívar Universidad Tecnológica de Bolívar Universidad Tecnológica de Bolívar Universidad Tecnológica de Bolívar Universidad Tecnológica de Bolívar

Edgardo Arrieta Ortiz Jaiver Campillo Sharicar Mendez Villamizar Erick Guerrero Jose Luis Villa Ramirez Justo Ramos Madrid

Comité Científico Prof. Jairo Useche , Ph.D. Prof. Jairo Tovar , Ph.D. Prof. Edgardo Arrieta, M.Sc. Prof. E. L. Albuquerque, Ph.D. Prof. Renato Pavanello, Ph.D. Prof. Paulo Sollero, Ph.D. Prof. Alejandro Marañon, Ph.D. Prof. José Rafael Toro, M.Sc.

Universidad Tecnológica de Bolívar, Col. Universidad Tecnológica de Bolívar, Col. Universidad Tecnológica de Bolívar, Col. Universidade Estadual de Campinas, Brasil. Universidade Estadual de Campinas, Brasil. Universidade Estadual de Campinas, Brasil. Universidad de los Andes, Col. Universidad de los Andes, Col.

Diagramación y Diseño Sharicar Mendez Villamizar Erick Guerrero

Impresión JaveGraft Ltda.

Contenido pg. 1

Análisis por Elementos Finitos de Sólidos Incompresibles: Implementación del Modelo de Mooney-Rivlin – Useche, J.

11

2

Determinación de las propiedades mecánicas de materiales compuestos reforzados con fibras continuas utilizando el método de los elementos de contorno - Álvarez, H.

21

3

Implementación de esquemas de optimización en Ansys para el diseño de estructuras de nodos articulados - Rodríguez, W., Pallares, M.

31

4

Análisis numérico de problemas de contacto tipo Hertziano y no Hertiziano en uniones mecánicas - Jiménez, D.

51

Predicción de la transformación de fases durante un proceso de temple utilizando un modelado numérico - Caraballo, D.

72

5

Medición de deformaciones estáticas y dinámicas en laminados de fibra de vidrio utilizando strain gages y adquisición de datos - Mulford, E.

83

6

Optimización de un Transom de lanchas de planeo fabricadas en fibra de vidrio- García, J.

102

7

Factores de intensificación de esfuerzos utilizando elementos finitos y elementos de contorno - López, J., Diaz, J.

115

8 9

Large deflection of composite laminate thin plates by the boundary element method - Albuquerque, E.L. Analysis of crack propagation in hyperelastic materials - Sollero, P.

125 133

10

Local and global instabilities in nano-size rectangular prismatic gold specimens - Pacheco, A.

142

11

Diseño óptimo de estructuras: caso de vigas totalmente esforzadas Rodríguez W., Pallares, M.

163

12

Modelación y chequeo estructural de puente Kantinwrwa sobre el río Ariguani - Chang, G.

177

13

Caracterización de Compuestos Laminados a partir de la simulación numérica 3D a nivel micromecánica y mesomecánico.

189

Prefacio EL SIMPOSIO REGIONAL SOBRE MECANICA DE MATERIALES Y ESTRUCTURAS CONTINUAS-SMEC es un espacio destinado al encuentro de investigadores nacionales dedicados al desarrollo de investigación básica y aplicada en las áreas de experimentación y modelado matemático y numérico de problemas en mecánica de materiales y de estructuras, que busca la apropiación y difusión de conocimiento, aportando así al desarrollo científico regional y nacional. El presente documento contiene las memorias del SMEC 2009. Son trece trabajos presentados en el evento, que tratan tematicas que van desde el análisis y optimización estructural, modelado de materiales hiperelasticos, pasando por problemas en mecánica de fractura, micromecánica de materiales compuestos, hasta el modelado numérico a escala atomica y molecular de materiales. SMEC 2009 fué organizado por el Grupo de Investigación en Materiales y Estructuras – GIMAT de la Universidad Tecnológica de Bolívar, el cual tuvo lugar en Cartagena de Indias – Colombia los dias 21, 22 y 23 de octubre de 2010 en el campus de la universidad. Esperamos que este documento permita alcanzar uno de los objetivos fundamentales de SMEC: difundir conocimiento avanzado en el área de materiales estructuras, que aporte al desarrollo cientifico y tecnológico de la región y del pais.

Jairo F. Useche, Ph.D. Editor

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

Análisis por Elementos Finitos de Sólidos Incompresibles: Implementación del Modelo de Mooney-Rivlin

Jairo F. Useche Departamento de Ingeniería Mecánica Universidad Tecnológica de Bolívar, Colombia [email protected]

Resumen Los materiales incompresibles, tales como los elastómeros, son ampliamente utilizados en la industria debido a sus excelentes propiedades para absorber de energía, alta flexibilidad, resiliencia elevada y vida útil prologada entre otras. El análisis de componentes mecánicos fabricados con elastómeros requiere de modelos constitutivos especializados y de la utilización de modelos por elementos finitos que tengan en cuenta las características altamente no-lineales de este tipo de materiales. Este trabajo presenta el desarrollo de una formulación utilizando el método de los elementos finitos – MEF para el análisis de materiales hiper-elásticos isotrópicos basados en el modelo de Mooney-Rivlin, utilizando una formulación lagrangiana actualizada (Updated Lagrangian). Los resultados numéricos obtenidos muestran una alta correlación con aquellos encontrados en la literatura, lo cual demuestra una alta confiabilidad de la formulación desarrollada. Palabras claves: Elementos finitos, elastómeros, hiperelasticidad, MooneyRivlin, deformaciones finitas.

1. INTRODUCCION Los elastómeros son ampliamente utilizados en la industria debido a sus excelentes propiedades tales como su alta capacidad para absorber energía, alta flexibilidad y variable rigidez, resiliencia elevada y vida útil prologada, entre otras propiedades. Los elastómeros, o cauchos, son materiales poliméricos elásticos cuyas dimensiones pueden cambiar en gran medida cuando se someten a esfuerzos Segundo Simposio Regional en Mecánica de Materiales y Estructuras Continuas – SMEC 2009

11

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

(Smith y Hashemi, 2004). Mediante un proceso de vulcanización las moléculas del polímero se unen mediante enlaces entrecruzados, formando moléculas más largas que restringen el movimiento molecular, lo que le confiere a este material baja resistencia y alta flexibilidad. El análisis de componentes mecánicos fabricados con elastómeros requiere de modelos constitutivos especializados y de la utilización de modelos por elementos finitos que tengan en cuenta las características altamente no-lineales de este tipo de análisis. A diferencia del análisis de componentes fabricados con materiales metálicos, los elastómeros presentan propiedades especiales que generan una respuesta mecánica compleja (MSC Software Corporation, 2000): -

Un elastómero puede desarrollar grandes deformaciones bajo carga, que pueden a ser del orden de 500% en deformación unitaria. Las ecuaciones esfuerzo-deformación (ecuación constitutiva) que rigen a estos materiales es no-lineal. Son materiales casi-incompresibles, es decir, su volumen con cambia apreciablemente con el esfuerzo. Presentan amortiguamiento viscoso (respuesta visco-elásticas) lo cual genera una respuesta que depende del tiempo y de la temperatura.

Estas características hacen del análisis estructural por elementos finitos de elastómeros una labor bastante compleja y área de investigación activa hoy en día. Muchos han sido los trabajos de investigación desarrollados en las últimas décadas en este campo (Mackerle, 2004). Este trabajo presenta el desarrollo de una formulación por elementos finitos para el análisis en hiper-elasticidad basada en el modelo constitutivo de Mooney-Rivlin. Este trabajo hace parte del proyecto FEASY (A Finite Element-Based Computational Program for Static and Dynamic Analysis of Solids and Structures) desarrollado en la actualidad por el Grupo de Investigación en Materiales y Estructuras de la Universidad Tecnológica de Bolívar en Cartagena - Colombia.

2. FORMULACION POR ELEMENTOS FINITOS 2.1 Principio del trabajo virtual La formulación por elementos finitos parte de la formulación débil de las ecuaciones de equilibrio estático para un sólido deformable en el dominio Ω delimitado por una superficie Γ:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

12

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

∫ (∇σ + f )δ vd Ω = 0

(1)

Ω

En esta ecuación, σ representa el tensor de esfuerzos de Cauchy, f es el vector de fuerzas de cuerpo por unidad de volumen y δv es un campo virtual de velocidades (ver Habbit y Sorensen, 1996). Utilizando el teorema de la divergencia y manipulando matemáticamente se obtiene la formulación débil para la Ec. (1) en la descripción espacial (descripción Euleriana):

∫ σ : δ dd Ω − ∫ f ⋅δ vd Ω − ∫ t ⋅ δ vd Γ = 0 Ω

Ω

(2)

Γ

donde d representa la parte simétrica del tensor gradiente de velocidad y t = σ.n es el vector de tracción, siendo n un vector normal a la superficie Γ. 2.2 Tratamiento de la incompresibilidad Diferentes formulaciones por elementos finitos han sido propuestas para tratar el problema del modelado de materiales incompresibles. Entre ellos se pueden destacar los abordajes basados en multiplicadores de Lagrange o aquellos basados en métodos de penalización, en donde la presión hidrostática en el elemento es considerada una variable adicional en el problema, generando así formulaciones mixtas (desplazamientos y presiones como incógnitas nodales). Desafortunadamente ese tipo de abordaje resulta en formulaciones que en general no son simples y eficientes desde un punto de vista numérico (Bonet y Wood, 1997). Una forma de superar este inconveniente es utilizar una formulación integral, basada en el principio variacional de Hu-Washizu, la cual está dada por:

∫ σ ' : δ dd Ω + pv.div (δ v)− ∫ f ⋅ δ vd Ω − ∫ t ⋅ δ vd Γ = 0 Ω

Ω

(4)

Γ

donde se ha considerado la descomposición aditiva del tensor de Cauchy: σ = σ '+ pI , siendo σ ' la componente desviatorica del tensor de esfuerzos. En este abordaje se asume que p y J son constantes en un volumen v y se define

div (δ v) como:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

13

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

div (δ v ) =

1 ∇δ vd Ω v ∫Ω

(5)

2.3 Linealizacion de la formulación integral En este trabajo, el método de Newton-Raphson es utilizado para obtener una solución aproximada a la Ec. (2). La aplicación de este método requiere la linealizacion en desplazamiento para esta ecuación, la cual está dada por (Bonet y Wood, 1997):

∫ δ v : c : εd Ω + ∫ σ : (∇u)

T

Ω

Ω

∇v  d Ω + κ v (divδ v)(divδ u) − ∫ f ⋅ δ vd Γ  Ω

 ∂x  ∂u  ∂x  ∂u  −∫ p  ⋅ ×δ v − ⋅  ×δ v dξdη = 0  ∂ξ  ∂η  ∂η  ∂ξ   Γξ

(6)

En esta ecuación, u es el vector de desplazamiento, p es la presión, ∂x ∂ξ y ∂x ∂η , son vectores tangentes a Γ de tal manera que:

n = (∂x ∂ξ ×∂x ∂η ) ∂x ∂ξ ×∂x ∂η , ξ ,η son coordenadas normalizadas y

κ = v V , siendo V el volumen de un elemento diferencial medido en la configuración de referencia. 2.3 Discretizacion por elementos finitos La discretizacion por elementos finitos de la Ec. (6) se lleva a cabo de la manera habitual. En el presente trabajo se utiliza una formulación lagrangiano actualizado (Updated-Lagrangian formulation) en la solución. De esta manera, el sistema de ecuaciones discretizado puede ser escrito como: t

K T ∆ u = t F − t + ∆t R

(7)

donde t K T es la matriz de rigidez tangente, ∆u es el vector incremental de t +∆t

R es el vector de fuerzas externas aplicadas el instante de desplazamiento, tiempo ∆t + t; y es el vector de fuerzas internas en un instante de tiempo t. La matriz de rigidez tangente está dada por: t

K T = tt K L + tt K NL + tt K v

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

(8)

14

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

Las matrices tt K L , tt K NL , tt K v y el vector de fuerza interna definidos como (ver Bathes, 1996 y Bonet y Wood, 1997): t t

F se encuentran

K L = ∫ tt BTL t C tt B L d t v t

t t

t

(9)

v

K NL = ∫ tt BTNL t σ tt B NL d t v t

(8)

v

t t

1    1  K V = tκ t v  t ∫ tt gd t v t ∫ tt gT d t v  v t v  v t v  

t t

F = ∫ tt BTL t σˆ tt d t v t

(9)

(10)

v

3. MODELO DE MOONEY-RIVLIN 3.1. Funcional de energía En su forma más general, el modelo de Mooney-Rivlin es expresado matemáticamente como (Ogden, 1997):

Ψ (C) = ∑ µ rs ( I C − 3) ( II C* − 3) r

s

(11)

r ,s

donde:

I C = tr (C) = C : I; II C* = 12 ( I C 2 − II C ) ; II C = tr (CC) = C : C

(12)

IC e IIC representan el primer y segundo invariante del tensor derecho de deformaciones de Cauchy. En este trabajo se considera un modelo de dos parámetros. De esta manera la Ec. (1) se reduce a:

Ψ ( I C , II C ) = µ10 ( I C − 3) + µ01 ( II C* − 3)

(13)

Para pequeñas deformaciones se tiene que: E = 6 ( µ10 + µ01 ) y G = µ10 + µ 01 donde E es el modulo de Hooke y G el modulo de rigidez del material. La Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

15

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

incompresibilidad del material puede ser forzada en el contexto de una formulación por elementos finitos a través del concepto de materiales cuasiincompresibles. En este caso, el potencial elástico es considerado como la descomposición aditiva desacoplada de ψ(C) de un potencial distorsional y un potencial volumétrico (ver Holzapfel, 2000):

ˆ (C) + U ( J ) Ψ (C) = Ψ

( )

ˆ ˆ (C) = Ψ C donde Ψ

(14)

ˆ = J −2/3C y J = det (F ) siendo F el tensor con: C

gradiente de deformación. La expresión más simple para U(J) está dada por:

U ( J ) = 12 κ ( J −1)

2

(15)

Donde κ un parámetro de penalización que permite forzar la compresibilidad del material. Físicamente este parámetro corresponde con el modulo volumétrico (bulk modulus) del material.

3.2. Tensor de esfuerzos A partir de la Ec. (4), el tensor de esfuerzos de Cauchy puede escribirse como (ver Bonet y Wood, 1997):

ˆ b + 4 J −1Ψ ˆ b2 + 2 J Ψ ˆ I + pI σ = 2 J −1Ψ I II III

(6)

ˆ = ∂Ψ ˆ ∂I , Ψ ˆ = ∂Ψ ˆ ∂II y Ψ ˆ = ∂Ψ ˆ ∂III y p = κ ( J −1) donde: Ψ I C II C III C ˆ = III −1/3C siendo IIIC es la presión hidrostática en el material. Considerando C C el tercer invariante de C, la Ec. (3) puede escribirse como:

ˆ ( I , II , III ) = µ ( III −1/3 I − 3) + 1 µ  III −2/3 I 2 − III −2/3 II − 6 Ψ C C C 10 C C C C C C 2 01  

(7)

De esta manera se tiene que:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

16

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

ˆ ∂Ψ = µ10 III C−1/3 + µ01 III C−2/3 I C ∂I C ˆ ∂Ψ = − 12 µ01 III C−2/3 ∂II C ˆ ∂Ψ = − 13 µ10 III C−4/3 I C − 32 µ01 III C−5/3 II C* ∂III C (8a,b,c) Reemplazando en la Ec. (6) se tiene:

σ = 2µ10 J −5/3 (b − 13 IC I) + 2µ01 J −7/3 ( IC I − b) b − 23 IIC* I + pI

(9)

3.3. Tensor elástico El tensor elástico en la configuración espacial es obtenido a partir de tensor elástico en la configuración material, el cual está dado por (Sussman, T.D., 1987): 0

C = µ10 ( 0 C1 ) + µ01 ( 0 C2 ) + κ ( 0 C3 )

(10)

donde: 0

C1 = 4 J −2/3 ( 13 I C ϒ − 13 I ⊗ C−1 − 13 C−1 ⊗ I + 19 I C C−1 ⊗ C−1 )

0

C 2 = J −4/3 (4I ⊗ I − 2 ( ϒ + Λ ))

− 83 J −4/3  I C I ⊗ C−1 − C ⊗ C−1 + I C C−1 ⊗ I

(10a)

(10b)

− C−1 ⊗ C − 14 J −2 II C ϒ − 53 II C C−1 ⊗ C−1  0

C3 =

p

κ

J (C−1 ⊗ C−1 − 2 ϒ)

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

(10c)

17

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

0.15

0.1

0.05

0

-0.05 0

0.05

0.1

0.15

0.2

0.25

0.3

0.15

0.1

0.05

0

-0.05

-0.1

-0.15

-0.2 -0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Figura 1. Superior: Modelo por elementos finitos de membrana elástica con agujero. Inferior: Malla deformada por la acción del desplazamiento impuesto.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

18

Useche, J., Análisis por elementos finitos de sólidos incompresibles … 0.025

0.02

Fx [N]

0.015

0.01

κ = 50 κ = 100

0.005

Resultado MATLAB( κ = 1000) Resultado ANSYS (κ = 1000) 0

0

0.01

0.02

0.03

0.04 0.05 0.06 Desplazamiento [m]

0.07

0.08

0.09

0.1

Figura 2. Resultados para desplazamiento horizontal vs. Carga en los extremos de la membrana, obtenidos con diferentes valores del modulo κ y comparación con ANSYS®.

Utilizando una operación push-forward sobre el tensor 0 C , el tensor elástico en la configuración espacial, t C , puede ser obtenido (Bonnet y Wood, 1997).

4. IMPLEMENTACION COMPUTACIONAL La formulación por elementos finitos desarrollada fue implementada en un código computacional desarrollado en MATLAB®. Se implementaron elementos isoparamétricos de cuatro nodos con interpolación bi-lineal y cuadrática para elasticidad plana y bricks isoperimétricos de 8 nodos con interpolación lineal para discretización en el espacio. Se utilizo una formulación incremental utilizando el método de Newton-Raphson en la solución del sistema de ecuaciones. La siguiente sección presenta el análisis de problemas de grandes deformaciones en elasticidad plana y espacial.

5. EJEMPLOS NUMERICOS 5.1 Deformación de una membrana con agujero central En este ejemplo se presenta los resultados obtenidos en el análisis de una membrana elástica con un agujero central. El problema es analizado considerando un estado de esfuerzo plano considerando µ10 = 0.4556 N/mm2 y µ01 = 0.960 Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

19

Useche, J., Análisis por elementos finitos de sólidos incompresibles …

N/mm2. Las dimensiones de la membrana son: L/W = 3, D/W = 0.1, t/W = 1, siendo D el diámetro del agujero y t el espesor de la membrana. Son utilizados 500 elementos cuadráticos de cuatro nodos para la discretización del dominio (ver figura 1). En esta figura se muestra igualmente la deformación de la membrana cuando es aplicado en los extremos de la misma un desplazamiento horizontal uniforme igual a 0.1 m. La figura 2 muestra las curvas fuerza horizontal vs. deformación para diferentes valores del modulo κ y una comparación con los resultados obtenidos utilizando el paquete comercial para análisis por elementos finitos ANSYS 10.0®.

CONCLUSIONES Fue presentada una formulación por elementos finitos para el análisis de deformación en materiales hiperelásticos utilizando el modelo constitutivo de Mooney-Rivlin con dos parámetros. La formulación desarrollada está basada en el principio variacional de Hu-Washizu en conjunto con la técnica de dilatación media para tratar numéricamente la incompresilibilidad el material. Los resultados obtenidos a través de la herramienta computacional implementada en MATLAB® muestran una alta concordancia con aquellos obtenidos utilizando otros programas de elementos finitos, demostrando así la confiabilidad de la herramienta desarrollada. REFERENCIAS − Bonet, J. Wood, R.D., 1997, Nonlinear Continuum Mechanics for Finite Element Analysis, Cambridge Press, New York. − Bathes, K.J., 1996, Finite Elements Procedures, Prentice Hall, New Jersey. − Holzapfel, G.A., 2000, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, John Wiley & Sons, New York. − Mackerle, J., 2004, Rubber and Rubber-like Materials, Finite-Element Analyses and Simulations, an Addendum: A Bibliography (1997–2003, Modelling Simul. Mater. Sci. Eng. 12 1031-1053. − MSC.Software Corporation, 2000, Nonlinear Finite Element Analysis of Elastomers, Technical Paper: www.axelproducts.com/.../MARC_FEA_ELASTOMERS_2000.pdf. − Ogden, R.W., 1997, Nonlinear Elastic Deformations, Dover, Publications, Inc., New York. − Smith, W.F., Hashemi, J., 2004, Fundamentos de la Ciencia e Ingeniería de Materiales, 4ed., McGraw-Hill, México D.F. − Sussman, T.D., 1987, On Finite Element Analysis of Non-linear Incompressible Solids, Ph.D. Thesis, Massachusetts Institute of Technology. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

20

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

Determinación de las Propiedades Mecánicas de Materiales Compuestos Reforzados con Fibras Continúas Utilizando el Método de los Elementos de Contorno Hugo Alvarez Departamento de Ingeniería Mecánica Universidad Tecnológica de Bolívar, [email protected]

Resumen Este trabajo presenta una metodología basada en el método de elementos de contorno (BEM) para la determinación de las propiedades mecánicas de materiales compuestos reforzados con fibra, a partir de las propiedades mecánicas de la matriz y de las fibras. La literatura especializada ofrece modelos analíticos basados en el análisis micromecánico del compuesto que permiten predecir sus propiedades mecánicas promedio a partir de las propiedades de cada uno de los constituyentes del mismo. Sin embargo, al ser modelos basados en hipótesis simplificadoras, resultan ser imprecisos y por tanto de poca aplicación práctica en ingeniería. Buscando desarrollar modelos micro-mecánicos más robustos, se han propuesto modelos basados en elasticidad lineal y mecánica del continuo, cuya solución ha sido obtenida a través del método de los elementos finitos. En este trabajo, se desarrolla un modelo matemático basado en las ecuaciones de Navier en elasticidad plana, para describir el comportamiento micro-mecánico de un material compuesto de fibra continua y unidireccional. La solución numérica de este modelo es obtenida utilizando el método de los elementos de borde desarrollado a partir de la formulación por residuos ponderados de las ecuaciones de Navier. Es utilizado el método de sub-regiones para el tratamiento de la discontinuidad en las propiedades mecánicas del material sobre el dominio de solución del problema. Son utilizados elementos de borde discontinuos con interpolación cuadrática para la discretización del contorno del dominio. El modelo desarrollado es aplicado a la determinación de las propiedades mecánicas de varios compuestos de uso comercial. Estos resultados son contrastados con aquellos reportados en la literatura especializada Palabras Clave: Ecuaciones de Navier, materiales compuestos, método elementos de contorno, método subregiones, micromecánica. Segundo Simposio Regional en Mecánica de Materiales y Estructuras Continuas – SMEC 2009

21

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

1. INTRODUCCIÓN Los materiales compuestos han tenido gran acogida como materiales de ingeniería en los últimos años y para muchas aplicaciones incluso han llegado a reemplazar a materiales tradicionales como el acero, por esta razón se ha desarrollado mucho conocimiento en esta área. En la literatura existente podemos encontrar una gran variedad de modelos para la predicción de las propiedades de estos materiales en los diferentes niveles de análisis los cuales son micromecánico, macromecánico y estructural, en este trabajo solo se analizará el primer nivel. Dentro del nivel micromecánico a la vez podemos encontrar diferentes aproximaciones como las basadas en mecánica de materiales o modelos analíticos, modelos empíricos, modelos semi-empíricos y modelos numéricos, cada uno de estos modelos tiene sus ventajas y desventajas por ejemplo en los modelos analíticos encontramos diferentes reglas de mezclas para la predicción de las propiedades mecánicas del compuesto las cuales son ecuaciones sencillas pero están muy limitadas debido a los supuestos que se realizaron durante su derivación como una geometría de empaquetamiento cuadrada o rectangular de las fibras en la matriz lo cual es poco usual en los materiales compuestos de uso común ya que en estos es normal encontrar una geometría de empaquetamiento aleatoria. En cuanto a la obtención experimental de las propiedades mecánicas solo podemos decir que en la mayoría de los casos es muy costosa debido a los equipos necesarios para la obtención de estos datos, por esta razón se desarrollaron modelos semi-empíricos que requieren poca experimentación ya que las pruebas experimentales solo son realizadas para la medición de ciertas constantes con las que luego se alimenta el modelo semi-empírico, algunos de estos modelos resultan en ecuaciones extensas y poco practicas como los modelos de contigüidad presentados por Tsai [1] aunque posteriormente Halpin-Tsai [3] presentaron modelos más simples y prácticos para el cálculo de las mismas propiedades de todas forma estos modelos requieren experimentación por ser modelos semiempíricos, están limitados a compuestos de fibras continuas unidireccionales y no tienen en cuenta factores como la presencia de vacios o incrustaciones en el material. Por otro lado encontramos los modelos numéricos los cuales han obtenido gran popularidad debido a los pasos agigantados que ha dado la tecnología de la computación. Entre los más populares se encuentran los modelos basados en elementos finitos y los basados en elementos de contorno, estos modelos permiten eliminar muchos de los supuestos realizados durante la derivación de las aproximaciones mencionadas arriba y tienen en cuenta factores muy importantes como la concentración de esfuerzo observada en las zonas cercanas a la interface fibra matriz. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

22

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

En este trabajo se decidió el uso del método de elementos de contorno debido a su menor costo computacional y menor dificultad de mallado en comparación con el método de elementos finitos. Además se decidió utilizar elementos cuadráticos discontinuos debido a que al utilizar elementos de alto orden, se puede lograr una buena precisión con pocos elementos y al usar elementos discontinuos se facilita la implementación computacional del modelo. Así utilizando este método se solucionaron las ecuaciones de Navier para elasticidad plana y utilizando el método de subregiones se modelo una celda unitaria o elemento representativo de material compuesto sometido a cargas de tracción y tomando los desplazamientos obtenidos como variable primaria, mediante la ley de Hooke se calculo el modulo elástico efectivo del elemento representativo. La mayoría de los modelos basados en elasticidad utilizan celdas como elementos representativos del material para analizar materiales compuestos, esto se hace con el fin de llevar a cabo el proceso de “homogenización” del material, si observamos la sección transversal de una lamina compuesta reforzada con fibras continuas unidireccionales a una escala suficientemente pequeña podremos distinguir las dos fases que conforman el materia (Fibra y matriz), para lidiar con estas inhomogeneidades se debe tomar un elemento representativo a una escala δ como se muestra en la figura1 tal que cuando se mida alguna propiedad del material en cualquier zona geométrica del mismo en la misma escala δ, se debe obtener el mismo valor en la propiedad obtenida, así podemos decir que el material es macroscópicamente homogéneo. Son muchos los trabajos realizados en este área que utilizan aproximaciones basadas en celdas por ejemplo Yiwei Jiang et. al. [4] Presentan una aproximación basada en micromecánica para la derivación de las ecuaciones constitutivas para el análisis de compuestos de tejidos planos utilizando celdas de volumen representativas bajo el supuesto de deformación y esfuerzos uniformes, las ecuaciones constitutivas son promediadas a lo largo de la dirección del espesor. M. G. Knight et. al. [5] utilizan el método de elementos de contorno y una aproximación basada en celdas embebidas para analizar el efecto de las propiedades de los constituyentes del compuesto y la distribución espacial de las fibras en el comportamiento localizado de un compuesto reforzado con fibras continuas unidireccionales. Wulf et. al. [6] utilizaron este mismo método para estudiar la evolución del daño alrededor de una grieta en materiales compuestos. Es claro que este método es el más adecuado para lidiar con las inhomogeneidades presentes en ese tipo materiales por lo que este artículo está orientado a ilustrar una metodología de análisis para predecir las propiedades elásticas transversales de materiales compuestos utilizando celdas embebidas. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

23

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

Figura 1 Elemento representativo homogenizado

Figura 2 Celda representativa básica

2. GEOMETRÍA Y MATERIALES Se modelo una celda unitaria básica como la que se muestra en la figura 2 tomando una fibra de tamaño 20µm de diámetro y la matriz cuidando una fracción de área de 0.55 para la fibra. Para este caso la forma como se encuentran empaquetadas las fibras en la matriz corresponde a un arreglo cuadrado, que es el caso de análisis más sencillo lo cual posibilita la comparación de los resultados obtenidos con modelos analíticos o semi-empíricos disponibles en la literatura.

2.1 Análisis por el método de elementos de contorno Para este estudio fue desarrollado un programa en MatLab® el cual resuelve las ecuaciones de Navier presentadas en la Ec. (1), para problemas de elasticidad plana.
    0  Ω


(1)

La Ec. (1), satisface las siguientes condiciones naturales y esenciales mostradas en Ec. (2) y Ec. (3):

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

24

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

    Γ

(2)

t    n    Γ

(3)

Para obtener la formulación por elementos de contorno de Ec. (1) primero debemos llevarla a su forma integral lo cual se logra utilizando el método de residuos ponderados con funciones de ponderación en desplazamientos del tipo  con lo que obtenemos: Ω ,   Ω  0

(4)

Luego integrando por partes o aplicando el teorema de Gauss-Green a Ec. (4) y aplicando las condiciones de contorno descritas en Ec.(2) y Ec.(3) obtenemos Brebbia, [7].            Γ  Ω  Ω

(5)

Donde  y # representan las soluciones fundamentales en desplazamiento y tracciones de las ecuaciones de Navier. Ec. (5) puede ser expresada de manera discreta como sigue: '

'

&  

,

$ %  $ )   $ *+ (

(

 

+(

(6)

& ./ y 0./ son matrices de coeficientes de influencia que se obtienen de las Donde siguientes integrales: &   $    Φ2 Γ % Γ 3

1

G   $ Γ u Φ2 Γ 3

(7)

1

donde e representa el número de elementos sobre los cuales se realiza la integral. Si despreciamos las fuerzas de cuerpo el termino B desaparece de Ec. (6) y el sistema de ecuaciones final obtiene la siguiente forma:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

25

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

-  0#

(8)

Donde es un vector que contiene los desplazamientos de cada nodo,  es un vector que contiene las tracciones en los nodos y H, G representan las matrices de coeficientes de influencia resultantes de la formulación por BEM.

3. FORMULACIÓN POR SUB-REGIONES Debido a que Ec. (8), solo es válida para un solo dominio homogéneo, se debe desarrollar otra formulación válida para materiales compuestos de varias fases. Para lograr esto debemos (i) obtener un conjunto de ecuaciones valido para cada fase, (ii) obtener ecuaciones de acoplamiento que relaciones la interacción entre las diferentes fases del material, (iii) crear un sistema de ecuaciones global que incluya las variables y condiciones de contorno de todas las fases. Para el caso de un material compuesto conformado solo de fibras y resina, solo existen dos fases que intervienen en el comportamiento del material y si analizamos una celda representativa como la ilustrada en la figura 2, solo tendremos dos subregiones involucradas en el análisis. Así si dividimos un sistema en dos subregiones Ω1 y Ω2 como se muestra en la figura 3, el contorno del sistema se divide en dos conjuntos de contornos para cada subregión como se muestra en la figura 4, de esta forma podemos escribir un conjunto de ecuaciones para cada subregión. Así: Para Ω1 : % 67  ) 7

(9)

Para Ω2 : % 67  ) 7

(10)

Para una fácil implementación se deben arreglar los sistemas de ecuaciones expresados Ec. (9) y Ec. (10), de tal forma que en los vectores 67 y 7 se encuentren primero los valores correspondientes a los nodos no compartidos y luego los compartidos, de igual forma las matrices de coeficientes de influencia H y G deben ser arregladas de forma tal que coincidan con el orden dado a los vectores 67 y 7. Como se muestra en la figura 4, Γ8 y Γ9 representan el contorno no compartido y compartido de cada subregión respectivamente.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

26

Alvarez, H., Determinación de las propiedades mecánicas de materiales compuestos reforzados…

Figura 3. Sistema dividido en dos subregiones

Figura 4. Contornos asociados con (a) subregión 1 (b) subregión 2 Teniendo en cuenta lo anterior Ec. (9) y Ec.(10) pueden ser reescritas como sigue: %88 : 98 %

)88 %89 678 99 ; < 9 =  : 98 % ) 67

%88 : 98 %

)88 %89 678 ; < =  : %99 )98 679

)89 78 ;> ? )99 79 )89 78 ;> ? )99 79

(11)

(12)

Ahora debemos escribir las ecuaciones de acoplamiento. Para escribir estas ecuaciones nos basaremos en: (i) los desplazamientos en los nodos de la zona compartida deben ser los mismos para ambas subregiones (ii) las tracciones en la zona compartida deben cumplir la condición de equilibrio estático. Esto lo podemos escribir de forma matemática como sigue: @A

678 67 BA C < 9 =  9 ?  5.5Å , the summation in Eq. (1) is carried out for those values of j for which r (ij ) < 6Å to reduce the computational cost. We note that the TB potential and its partial derivatives with respect to r (ij ) are continuous at the cut off radius of 6Å within the accuracy of the machine. The potential energy V of a system of atoms equals the sum of the energy V(i) of all atoms in the system. That is, N

V = ∑V (i )

(2)

i =1

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

146

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

The interaction force vector f(ij) between atoms i and j equals the negative of the partial derivative of the potential energy with respect to r(ij), or

 ∂V (i ) ∂V ( j )  r (ij ) fα(ij ) = −  (ij ) + ( ji )  α(ij ) ∂r  r  ∂r

(3) (ij)

Here and below, the index α ranges from 1 to 3, and fα equals the component of f(ij) along the xα -coordinate axis of a rectangular Cartesian coordinate system.

2.2 Average Stresses For a system comprised of N atoms at 0 K, average values σ αβ of components of the Cauchy stress tensor are computed from the relation [35]

σ αβ =

1 2ΩT

N

N

∑∑ fα

( ij ) ( ij )

i =1 j =1 j ≠i

rβ

(4)

where ΩT equals the volume occupied by the system. Eq. (4) corresponds to the configurational part of the virial stress tensor. For a continuous body, the average value over volume Ω of the Cauchy stress tensor defined by

σ αβ =

1 σ αβ d Ω Ω Ω∫

(5)

can be written as

σ αβ =

1 1 σ αφ x β nˆφ dS = ∫ tα xβ dS ∫ Ω ∂Ω Ω ∂Ω

(6)

where nˆ is a unit outward normal to the boundary ∂Ω of Ω , t = nˆ ⋅ σ is the surface traction, and we have used the divergence theorem and the balance of linear momentum with null body forces. Thus the average Cauchy stress tensor multiplied by the volume of the region occupied by the body equals the first moment of tractions acting on the bounding surfaces of the body. For a discrete system, Eq. (5) can be written as,

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

147

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

σ αβ =

1 ΩT

Nb

∑ xβ

(i )

i =1

fα( i )

(7)

where Nb equals the number of atoms on the bounding surface of the region whose deformations are being studied. Assuming that the volume assigned to each atom is the same, Eq. (4) becomes

σ αβ =

1 N

N

1 ∑ (i ) i =1 2Ω

(i ) where ωαβ =

1 2Ω ( i )

N

∑ j =1 j ≠i

fα(ij ) rβ( ij ) =

1 N

N

∑ ωαβ

(i )

i =1

(8)

N

∑ fα

( ij ) ( ij )

j =1 j ≠i

rβ

is the dipole force tensor [36].

However, for a finite size specimen, Eq. (8) is approximately valid since the volume assigned to an atom on the bounding surface equals ½ of that assigned to an atom in the interior of the body, and the volume of an atom at a vertex of the region equals 1/8th of that for an interior atom. For a system having a large number of atoms, Eqs. (4) and (8) will give essentially the same values of average stresses. However, for a small system with a large number of atoms on the bounding surfaces compared to those in the interior, average stresses computed from these two equations will differ.

2.3 Strains We employ the MSPH method [23] to compute the spatial distribution of the deformation gradient F and spatial gradients, G, of F from positions of atoms in the current and the reference configurations. The Cauchy-Born rule [38-40] states that for a crystal with a simple Bravais lattice, relative position vectors r(ij) and R(ij) between atoms i and j in the current and the reference configurations are related by r(ij)= F(i)R(ij) where F(i) is the deformation gradient at the position of atom i in the reference configuration. We use components of the tensor G, the second-order spatial derivatives of displacements, to characterize local instabilities in an atomic system. A detailed description of the methodology is given in [16]. From F(i) at the point x(i), we evaluate there the Almansi-Hamel strain tensor ε(i) from

(

(i ) ε αβ = (1 / 2) δ αβ − ( F −1 ) ( F −1 ) ϕα ϕβ (i )

(i )

)

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

(9) 148

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

where δ αβ is the Kronecker delta. The volume averaged value, ε , of this tensor for the system is defined by

ε αβ =

1 ΩT

N

∫ ε αβ ( x )d Ω = ∑ i =1

ΩT

Ω(i ) (i ) ε αβ ΩT

(10)

where Ω (i ) and ΩT equal, respectively, the volume assigned to point i and the total volume of the system in the deformed configuration. We set Ω (i ) equal to the Voronoi volume associated with atom i. An approximation of the Voronoi volume is given by (e.g., see [33])

∑( r ) Ne

4π 3 Ω = ai , 3 (i )

j =1 j ≠i v Ne

ai = k

(ij ) −1

∑( r (ij ) )

−2

(11)

j =1 j ≠i

Here Ne equals the number of atoms in the neighborhood of atom i for which

r

( ij )

≤ ( 3 / 2)a0 , a0 is the lattice parameter, and we set the constant

k v = 0.55

found by computing the Voronoi volume of an atom at the centroid of the specimen, and equating it to the volume given by Eq. (11).

2.4 Molecular mechanics simulations We start numerical simulations by assigning the initial position vector XI(i) of each atom in the system in a perfect lattice configuration. Without applying any external force, each atom is allowed to move freely till the potential energy of the system has been minimized by using the conjugate gradient (CG) with warranted descent technique [42]. The minimization procedure is stopped when the magnitude of each component of the gradient of the internal energy at every atom in the system equals at most 1x10-8 eV/Ǻ. The position vector of an atom in this relaxed configuration is denoted by XR(i), and this configuration is taken as the reference configuration. Subsequently, after each increment in the prescribed displacements of atoms on the end faces of the specimen, the total potential energy is minimized. The change in the potential energy of the system from that in the reference configuration equals the strain energy required to deform the body or the system of atoms. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

149

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Table 5.1: Aspect ratios, number of atoms, and lengths of specimens tested in tension and compression.

L/H

No. of atoms

1 3 5 10 20

3430 9928 16787 32671 65883

Length (Å) ~37 ~110 ~188 ~367 ~742

The process is continued till atoms on the end faces have been given the desired axial displacement. The width H of specimens is ~37 Å and their lengths range from ~37 Å for L/H = 1 to ~742 Å for L/H = 20; the number of atoms in a specimen is listed in Table 5.1. In each case the specimen is oriented with the coordinate planes {1,0,0}, {0,1,0} and {0,0,1}.

3. INSTABILITIES IN TENSILE AND COMPRESSIVE DEFORMATIONS 3.1 Simple tension/compression tests For different L/H ratios, Fig. 1 depicts the variation with ε of the average value of σyy component of the Cauchy stress tensor found by using Eq. (4). It is observed that the variation with ε of the average value of σyy is the same for specimens with L/H ≥ 3. For the simple tension/compression simulations, the average values of σzz and σxx are negligible as compared to the average values of σyy. The applied boundary conditions allow atoms on the end faces to move freely in the X- and the Z-directions; consequently, edge effects are negligible and averaged values of σzz and σxx are very small. Average values of all shear stresses are negligible up to the discontinuity in the σyy vs. ε curves. In Table 2 we have listed, for different L/H ratios, values of the average σyy stress and the average axial strain at the yield point. Values of the axial yield stress and the average axial strain at yield for 10 ≥ L/H ≥ 3 for the simple tension case are ~ 5 GPa and ~ 8% respectively; the corresponding values for the tension test listed above are ~ 6.2 GPa and ~ 9.8%. In simple compression, a dependence of the average axial stress and the average axial strain at yield on the L/H ratio is also observed.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

150

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Simple tension / compression

6 5

Axial stress σyy [GPa]

4 3

L/H = 20 L/H = 10 L/H = 5 L/H = 3 L/H = 1

2 1 0 -1 -2 -3 -0.1

-0.05

0

[ Å/Å ] Axial strain ε [A/A]

0.05

0.1

Figure 1. For different values of L/H, evolution with the average axial strain ε of the average value of σyy component of the Cauchy stress tensor for the simple tension/compression tests.

For a given value of L/H, the yield stress in simple tension is higher than that in simple compression. Because of residual stresses in the reference configuration, the difference in the yield stress in simple tension and simple compression cannot be attributed to the Bauschinger effect. As evidenced from results exhibited in Fig. 2 a good agreement between three measures of average stresses is found but in this case there is more scatter in the average stresses derived from the first moment of forces on bounding surfaces than that from the other two methods. For L/H = 10, the evolution of the average normal components of the AlmansiHamel strain tensor with the average axial strain is plotted in Fig. 3. Curves for εxx vs. ε and εzz vs. ε overlap each other, and εyy is almost equal to ε till the specimen yields. Thus the effect of geometric nonlinearities can be ignored until the yield point. In Fig. 4 we have plotted the averagae axial stress vs. the average axial strain during unloading from two configurations – one just before the drop in the average axial stress vs. the average axial strain curve and the other just after this drop. When the specimen is unloaded from the configuration just before the average axial stress drops noticeably, the average axial stress vs. the average axial strain curve during unloading overlaps that during loading suggesting that the specimen deformed elastically. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

151

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Table 2: Values of the average axial stress and the average axial strain at the yield point for specimens with different L/H ratios deformed in simple tension and compression.

Simple tension L/H 1 3 5 10 20

Simple compression

σ yyyield

ε yyyield

σ yyyield

ε yyyield

(GPa) 5.134 5.050 4.996 4.990 4.618

(%) 7.928 8.119 8.117 8.086 7.584

(GPa) -2.498 -1.784 -1.810 -1.671 -1.387

(%) -8.345 -6.643 -6.897 -5.157 -3.874

Simple tension / compression

5

Axial stress σyy [GPa]

4

Eq.(4)(2.5) (2.7) Eq. Eq. Eq.(7)(3.10) Force/area FORCE/AREA

3 2 1 0 -1 -2 -0.06

-0.04

-0.02

0

0.02

0.04

Axial strain ε [[A/A] Å/Å ]

0.06

0.08

0.1

Figure 2. For L/H=10, comparison between different measures of the average Cauchy stress tensor in the simple tension/compression tests, (Eq. 2.7), (Eq. 3.10) and the mechanics of materials approach (force/area).

However, when the specimen is unloaded from the configuration just after the severe drop in the axial stress , there is a residual average axial strain at zero average axial stress. It confirms that the specimen deformed plastically during the instant the average stress dropped, and the permanent average axial strain equals 4.5%. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

152

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Simple tension / compression

Alm ansi - H am el strains [A/A] [ Å/Å ]

0.08 0.06 0.04

ε

yy

ε

xx

ε

zz

0.02 0 -0.02 -0.04 -0.06 -0.04

-0.02

0

0.02

0.04

Axial strain ε [A/A] [ Å/Å ]

0.06

0.08

Figure 3. For L/H=10, evolution with the average axial strain ε of the average values of normal components of the Almansi-Hamel strain tensor for the simple tension/compression tests. Curves for εxx vs. ε and εzz vs. ε overlap each other.

6 Loading Unloading-I Unloading-II

Normal stress σyy [GPa]

5 4 3 2 1 0 -1 -0.02

0

0.02

0.04

0.06

0.08

0.1

Axial strain ε [A/A]

Figure 4. For L/H=10, variation with the average axial strain ε of the average axial stress for loading and unloading paths in simple tension.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

153

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Minimum eigenvalue local Hessian

0.5

0

-0.5

-1

-1.5

-2 -0.1

L/H = 1 L/H = 3 L/H = 5 L/H = 10 L/H = 20

-0.05

0

Axial strain ε [A/A] [Å/ Å]

0.05

0.1

Figure 5. For simulations with different values of L/H, variation with the average axial strain ε of the minimum eigenvalue of the local Hessian H(i) in simple tensile/compressive deformations.

3.2 Analysis of local instabilities For different values of L/H, Fig. 5 exhibits the evolution of the minimum eigenvalue of the local Hessian H(i) among all atoms in the system for simple tensile/compressive deformations. For simple compression, the minimum eigenvalue continuously decreases and remains positive until the strain level where the sharp drop in the average axial stress - average axial strain curve occurs. The strain levels at which local instabilities, signified by the minimum eigenvalue of H(i) becoming negative, appear correspond to

ε yyyield (see Table 2). No local

instability occurred prior to this strain level. However, for simple tensile deformations, and for all values of L/H considered here, a group of atoms in each of the eight corners of the sample become unstable at an average axial strain of ~ 6% when curves in Fig. 5 exhibit the first discontinuity. The minimum eigenvalue remains negative and continuously decreases up to ~ 8% average axial strain when the minimum magnitude of the negative eignevalue drops noticeably; at this strain level the discontinuities in the average axial stress vs. the average axial strain curves also occur (cf. Fig. 1). Values of

ε yyyield listed in Table 2 correspond to the

initiation of the second sharp drop in the magnitude of the minimum eigenvalue of int s H(i). For both simple tensile and simple compressive deformations, ε yy ≤ Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

154

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens… int s ε yyyield for the five values of L/H considered here. The ε yy equals the average

axial strain when the minimum eighenvalue of H(i) becomes negative. These local instabilities do not cause, on the average, permanent deformations since overall deformations are reversible as should be clear from results exhibited in Fig. 4. For L/H = 10, Fig. 6 shows distributions of the minimum eigenvalue of H(i) , the CNP parameter and G for the simple compression at ε = ε yy

yield

greater than

ε

yield yy

and at ε slightly

. Negative eigenvalues of H(i) occur at points located in planes

{1,1,1} of high atomic density. The unstable atoms are not distributed symmetrically about the mid-section, Y = L/2, of the specimen (see Fig. 6b) possibly due to the asymmetry in the boundary conditions applied to the end faces and numerical truncation errors. Whereas atoms on one face are kept stationary, those on the other end face are displaced axially. Fig. 6d displays the distribution of the CNP at ε = ε yy

yield

. The CNP equals either zero or has very small values at

stable points but has large values at generally the same interior atomic positions where the local instability was predicted by the negative eigenvalues of H(i). At ε =

ε yyyield values of the CNP parameter are negligible everywhere except at some points located on the bounding surfaces. Fig. 6e,f exhibits the distribution of G at ε = ε yy

yield

. High values of G occur

at points close to the lateral surfaces and at points located near the end faces where displacements are prescribed. Values of G vanish at points in the interior of the specimen whose distance from the end faces exceeds ~37 Å, i.e., the width of the sample. For ε> ε yy

yield

high values of G

occur at the same atomic positions

(i)

where the minimum eigenvalue of H is negative (cf. Fig. 6f) implying that deformations are highly inhomogeneous in the neighborhoods of atoms that have become unstable. The sudden appearance of instabilities during an incremental axial strain of 0.01% suggests that either none of the three criteria used to characterize instability is robust enough to detect a gradual progression to an unstable state or the output time interval should have been considerably reduced. The latter option requires running the simulations several times with output at successively smaller intervals in order not to exceed the memory allocation in the computer. However, it has not been followed since for practical purposes instability strain within 0.01% is reasonably accurate.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

155

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

Fig. 7 depicts, for L/H=10, distributions of the minimum eigenvalue of H(i), the CNP and G for simple tensile deformations. At ε slightly less than

ε yyyield there

are very few atoms near the eight vertices of the specimen that have become unstable. However, at ε just greater than ε yy

yield

, several atoms have become

unstable. The distribution of unstable points has a serrated pattern along specimen’s Y-centroidal line formed by atoms located on planes of high atomic density. The instabilities under compression and simple compression do not propagate through all of the specimen length. A similar serrated pattern was reported by Liang and Zhou [42] who performed MD simulations of tensile deformations of Cu nanowires at 300 K with specimens having the same crystallographic orientations as in the present work. The distribution of unstable points predicted by the non-vanishing values of the CNP also coincides with that predicted by the minimum eigenvalue of the local Hessian (see Figs. 7b and 7d) becoming negative. The distribution of

G

depicted in Fig. 7e shows patterns symmetric with respect to the Y-centroidal line of the specimen. At ε= ε yy

yield

values of G

increase from almost zero at the

centroidal line of the specimen to ~0.07 at atoms located on the traction free lateral surfaces. The regions of high values of G located at points close to the end faces of the specimen observed in compression and simple compression are not present in specimens deformed in tension. At ε> ε yy

yield

high values of G

occur at

(i )

numerous points where the minimum eigenvalue of H is negative.

3.3 Global Instabilities Results exhibited in Fig. 8 reveal that the variation with the average axial strain of the strain energy density for 3 ≤ L/H ≤ 20 is essentially independent of the aspect ratio L/H in tension but not in compression. The avearage axial strain when the system becomes gloabally unstable in simple compression depends strongly upon the value of L/H.

CONCLUSIONS We have used molecular mechanics simulations with the tight-binding potential to study local and global instabilities in initially defect-free nano-specimens of gold deformed in simple tension/compression. The criteria used to delineate local Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

156

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

instabilities are: (i) a component of second order spatial partial derivatives of the displacement field having large values relative to its average value in the body, (ii) the minimum eigenvalue of the Hessian of the potential energy of an atom becoming non-positive, and (iii) a high value of the common neighborhood parameter. The system’s configuration is said to be globally unstable when its potential energy density changes significantly with an infinitesimal increase in the average axial strain. Conclusions from this work are: • • • • • •

•

The three criteria for the initiation of a local instability are met essentially simultaneously at the same atomic positions. The average values of the Cauchy stresses derived from different definitions of the Cauchy stress tensor agree well with each other. The response of a specimen in tension and compression is very different. This can be attributed to the presence of non-uniformly distributed stresses in the reference configuration. For specimens deformed in simple tension the average axial strain at the initiation of local instabilities is noticeably less than that when the specimen yields or that when it becomes globally unstable. Atoms on the end faces do not become unstable. Furthermore, atoms that become unstable are located away from the end faces where essential boundary conditions are imposed. The average axial stress vs. the average axial strain curve up to the yield point is nonlinear but deformations are elastic in the sense that the curve is reproduced during unloading. However, after the specimen has yielded (or has become globally unstable) there is a permanent axial strain induced in the specimen (~4.5% for L/H=10). The slope of the average axial stress vs. the average axial strain curve during unloading is the same as that of the curve during loading. The average axial stress at yield in simple compression is considerably less than that in simple tension.

ACKNOWLEDGEMENTS This work was partially supported by the ONR grant N00014-99-06-1-0567 to Virginia Polytechnic Institute and State University with Dr. Y.D.S. Rajapakse as the program manager, and Virginia Tech’s Institute of Critical Technologies and Sciences. AAP was also supported by the COLCIENCIAS – LASPAU scholarship, and the Universidad del Norte in Barranquilla, Colombia. Views expressed herein are those of authors, and neither of the funding agencies nor of their institutions. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

157

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

a

b

c

d

e

f

Figure 6. For the simple compression test, distribution of the minimum eigenvalue of the local Hessian, the CNP and

G

in the specimen with L/H=10; (a) minimum eigenvalue,

ε = -5.15%; (b) minimum eigenvalue, ε = -5.16%; (c) CNP on the mid-section X = 18 Å, ε = -5.15%; (d) CNP, ε = -5.16%; (e)

G

G

on the mid-section X = 18 Å, ε = -5.15%; and (f)

at unstable points, ε = -5.16%; In Fig. (d) atoms on the bounding surface have been

removed to clearly depict values of the CNP at unstable points in the interior of the specimen.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

158

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

c

d

a

b

e

f

Figure 7. For the simple tension test, distribution of the minimum eigenvalue of the local Hessian, the CNP and

G

in the specimen with L/H=10; (a) minimum eigenvalue, ε =

8.08%; (b) minimum eigenvalue, ε = 8.15%; (c) CNP on the mid-section X = 18 Å, ε = 8.08%; (d) CNP, ε = 8.15%; (e)

G

G

on the mid-section X = 18 Å, ε = 8.08%; and (f)

at unstable points, ε = 8.15%; in Fig. (d) atoms on the bounding surface have been

removed to clearly show values of the CNP at unstable points in the interior of the specimen.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

159

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

-4

12

x 10

Simple tension / compression

[ eV/Å33] ] Strain energy density [eV/A

10 L/H = 20 L/H = 10 L/H = 5 L/H = 3 L/H = 1

8 6 4 2 0 -2 -0.1

-0.05

0

[ Å/Å ] Axial strain ε [A/A]

0.05

0.1

Figure 8. Variation with the average axial strain ε of the strain energy density in simple tension/compression tests.

REFERENCES [1] J. Diao, K. Gall, M.L. Dunn, Atomistic simulation of the structure and elastic properties of gold nanowires, J. Mech. Phys. Solids 52 (2004) 1935-1962. [2] W. Zhang, T. Wang, X. Chen, Effect of surface stress on the asymmetric yield strength of nanowires, J. Appl. Phys. 103 (2008) 123527. [3] M.I. Baskes, Modified embedded-atom potentials for cubic materials and impurities, Phys. Rev. B 46 (1992) 2727-2742. [4] K. Gall, J. Diao, M.L. Dunn, The strength of gold nanowires, Nano Lett. 4 (2004) 2431-2436. [5] J. Diao, K. Gall, M.L. Dunn, J.A. Zimmerman, Atomistic simulations of the yielding of gold nanowires, Acta Mater. 54 (2006) 643-653. [6] Y. Liu, E. Van der Giessen, A. Needleman. An analysis of dislocation nucleation near a free surface, Int. J. Solids Struct. 44 (2007) 1719-1732. [7] R.E. Miller, D. Rodney, On the nonlocal nature of dislocation nucleation during nanoindentation, J. Mech. Phys. Solids 56 (2008) 1203–1223. [8] C. Truesdell, W. Noll W, The Non-Linear Field Theories of Mechanics, Second edition. Springer- Verlag (1992). [9] R. Hill, Acceleration waves in solids, J. Mech. Phys. Solids 10 (1962) 116. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

160

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

[10] K.J. Van Vliet, J. Li, T. Zhu, S. Yip, S. Suresh, Quantifying the early stages of Plasticity through nanoscale experiments and simulations, Phys. Rev. B 67 (2003) 104105. [11] P. Steinmann, A. Elizondo, R. Zunyk, Studies of validity of the CauchyBorn rule by direct comparison of continuum and atomistic modeling, Modell. Simul. Mater. Sci. Eng. 15 (2007) S271-S281. [12] J. Lu, L. Zhang, An atomistic instability condition and applications, J. Mech. Mater. Strut. 1 (2006) 663-648. [13] T. Kitamura, Y. Umeno, N. Tsuji, Analytical evaluation of unstable deformation criterion of atomic structure and its application to nanostructure, Comput. Mater. Sci. 29 (2004) 499-510. [14] T. Kitamura, Y. Umeno, R. Fushino, Instability criterion of inhomogeneous atomic system, Mater. Sci. Eng. A 379 (2004) 229-233. [15] S.V. Dimitriev, T. Kitamura, J. Li, Y. Umeno, K. Yashiro, N. Yoshikawa, Near-surface lattice instability in 2D fiber and half-space, Acta Mater. 53 (2005) 1215-1224. [16] A.A. Pacheco, R.C. Batra, Instabilities in Shear and Simple Shear Deformations of Gold Crystals, J. Mech. Phys. Solids 56 (2008) 3116-3143. [17] C. Kelchner, S. Plimpton, J. Hamilton, Dislocation nucleation and defect structure during surface indentation, Phys. Rev. B 58 (1998) 11085-11088. [18] J.A. Zimmerman, C.L. Kelchner, P.A. Klein, J.C. Hamilton, S.M. Foiles, Surface step effects on nanoindentation, Phys. Rev. Lett. 87 (2001) 165507. [19] Q.K. Li, M. Li, Atomic scale characterization of shear bands in an amorphous metal, Appl. Phys. Lett. 88 (2006) 241903. [20] H. Tsuzuki, P.S. Branicio, J.P. Rino, Structural characterization of deformed crystals by analysis of common atomic neighborhood, Comput. Phys. Commun. 177 (2007) 518–523. [21] R. Pasianot, Z.Y. Xie, D. Farkas, E.J. Savino, Representation of atomistic dislocation core structures, Scr. Metall. Mater. 28 (1993) 319-324. [22] C.S. Hartley, Y. Mishin, Representation of dislocation cores using Nye tensor distributions, Mater. Sci. Eng. A 400-401 (2005) 18-21. [23] G.M. Zhang, R.C. Batra, Modified smoothed particle hydrodynamics method and its application to transient problems, Comput. Mech. 34 (2004) 137146. [24] R.J. Hardy, Formulas for determining local properties in molecular dynamics simulations: shock waves, J. Chem. Phys. 76 (1981) 622-628. [25] M.W. Finnis, J.E. Sinclair, A Simple empirical N-body potential for transition-metals, Philos. Mag. A 50 (1984) 45-55. [26] M.I. Baskes, M.S. Daw, Embedded atom method: derivation and application to impurities, surfaces and other defects in metals. Phys. Rev. B 29 (1984) 6443-6453.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

161

Pacheco, A.A., Batra, R.C., Local and global instabilities in nano-size rectangular gold specimens…

[27] F. Voter, S.P. Chen, Characterization of Defects in Materials, In: Siegel R.W., Weertman J. R., Sinclair R. (Eds.), Mater. Res. Soc. Symp. Proc. 82. Materials Research Society, Pittsburgh (1987) 175. [28] K.W. Jacobsen, J.K. Norskov, M.J. Puska, Interatomic interactions in the effective-medium theory, Physical Review B 35 (1987) 7423-7442. [29] F. Ercolessi, M. Parrinello, E. Tosatti, Simulation of gold in the glue model, Philos. Mag. A 58 (1988) 213-226. [30] F. Cleri, V. Rosato, 1993. Tight-binding potentials for transition metals and alloys, Phys. Rev. B 48 (1993) 22-33. Ju, S.P., Lin J.S., Lee W.J., 2004. A molecular dynamics study of the tensile behavior of ultrathin gold nanowires. Nanotechnology 15, 1221-1225. [31] W.J. Lee, S.P. Ju, S.J. Sun, M.H. Weng, Dynamical behavior of 7-1 gold nanowires under different axial tensile strains, Nanotechnology 17 (2006) 32533258. [32] J.S. Lin, S.P. Ju, W.J. Lee, Mechanical behavior of gold nanowires with a multishell helical structure, Phys. Rev. B 72 (2005) 085448. [33] Q. Pu, Y. Leng, L. Tsetseris, H.S. Park, S.T. Pantelides, P.T. Cummings, Molecular dynamics simulations of stretched gold nanowires: The relative utility of different semiempirical potentials. J. Chem. Phys. 126 (2007) 144707. [34] G. Rubio-Bollinger, S.R. Bahn, N. Agrair, K.W. Jacobsen, S. Vieira, Mechanical properties and formation mechanisms of a wire of single gold atoms, Phys. Rev. Lett. 87 (2001) 026101 [35] J.A. Zimmerman, E.B. Webb, J.J. Hoyt, R.E. Jones, P.A. Klein, D.J. Bammann, Calculation of stress in atomistic simulation, Modell. Simul. Mater. Sci. Eng. 12 (2004) S319-S332. [36] G.P. Potirniche, M.F. Horstemeyer, G.J. Wagner, P.M. Gullett, A molecular dynamics study of void growth and coalescence in single crystal nickel, Int. J. Plast. 22 (2006) 257-278. [37] P.H. Mott, A.S. Argon, U.W. Suter, The atomic strain tensor, J. Comput. Phys. 101 (1992) 140-150. [38] M. Born, K. Huang, Dynamical Theory of Crystal Lattices, Oxford, Clarendon (1954). [39] J.L. Ericksen, On the Cauchy-Born rule, Math. Mech. Solids 13 (2008) 199-220. [40] I. Stakgold, The Cauchy relations in a molecular theory of elasticity, Quart. Appl. Math. 8 (1950) 169-186. [41] W.W. Hager, H. Zhang, A new conjugate gradient method with guaranteed descent and an efficient line search, SIAM Journal on Optimization 16 (2005) 170 – 192. [42] W. Liang, M. Zhou, Response of copper nanowires in dynamic tensile deformation, Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci. 218 (2004) 599606.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

162

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

Diseño óptimo de estructuras: caso de vigas totalmente esforzadas 1

Wilson Rodríguez Calderón Ing. Civil, Mag. Métodos Numéricos para Ingeniería - UPC Universidad de la Salle [email protected] 2

Myriam Rocío Pallares Muñoz Ing. Civil, Mag. Métodos Numéricos para Ingeniería - UPC Universidad Santo Tomás [email protected]

Resumen Para reducir al mínimo la cantidad de material y contar con una viga lo más liviana posible, se pueden variar las dimensiones de sus secciones transversales para que en cada una ocurra el esfuerzo de flexión máximo permisible. Una viga sometida a estas condiciones se denomina viga totalmente esforzada ó viga de resistencia constante. Ciertamente, casi nunca se alcanzan estas condiciones ideales, a causa de problemas prácticos de construcción y por la incertidumbre propia de las cargas de diseño. No obstante, el conocimiento de las propiedades de una viga totalmente esforzada tiene un papel importante en el diseño de estructuras para obtener pesos mínimos. En el presente trabajo se plantea la determinación de secciones de vigas totalmente esforzadas empleando técnicas de optimización, para lo cual es necesario identificar la función objetivo a minimizar, las variables de diseño y las restricciones del problema. Se reporta un caso de estudio de una estructuras diseñadas para mantener un esfuerzo máximo casi constante, estas son, las láminas de los muelles en los automóviles. La validación se realiza comparando los resultados de modelos de elementos finitos acompañados de los respectivos esquemas de optimización en Ansys y la solución teórica. Este trabajo pone de manifiesto que el objetivo de cada problema de ingeniería es obtener un diseño final, el mejor posible, el óptimo, y esto se realiza en dos etapas: el cálculo de la estructura y la modificación de las soluciones intermedias mediante procedimientos matemáticos, sin decisiones subjetivas.

1

Profesor investigador Programa de Ingeniería Civil Universidad de la Salle. Líder Centro de Investigación en Modelación Numérica y Desarrollo de Software CAE, “CIMON”. Bogotá D.C. Profesora investigadora Facultad de Ingeniería Civil Universidad Santo Tomás. Líder Grupo de Investigación en Simulación y Control Numérico, “SICON”. Bogotá D.C. 2

Segundo Simposio Regional en Mecánica de Materiales y Estructuras Continuas – SMEC 2009

163

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1. INTRODUCCION Cuando se quiere determinar las dimensiones de una viga prismática, normalmente se diseña para que resista el momento de flexión máximo Mmax, sin embargo, el momento generalmente varia a lo largo de la longitud de la viga y por tanto muchas secciones de la viga realmente estarán sobrediseñadas, sobre todo si se tiene cierta seguridad sobre las cargas que recibirá la estructura. Para optimizar el diseño y reducir el peso e inevitablemente los costos, es posible plantear un diseño de sección variable, tal que, a lo largo de la viga se presente en cada sección el esfuerzo máximo permisible o alguna condición de resistencia ajustada a normas o especificaciones propias del material y de la construcción misma de la viga. Este tipo de vigas son mas comunes en máquinas por la facilidad de mecanizado y de colocación, sin embargo, en estructuras de ingeniería civil también pueden emplearse con cierto cuidado, ya que, se sabe que las estructuras propias de las edificaciones o puentes presentan cierta incertidumbre en las cargas que recibirán sobre todo en lo que se refiere a cargas sísmicas y de viento. Algunos ejemplos de formas constructivas pueden ser las vigas acarteladas, vigas arriñonadas o de riñón de bóveda y vigas compuestas de acero reforzadas con cubreplacas soldadas. El análisis de esfuerzos en vigas no prismáticas, tienen cierta complejidad de cálculo, por tanto, es recomendable usar métodos numéricos como el método de los elementos finitos, dado que, en este es posible incluir una discretización suficiente que permite aproximar de manera adecuada el comportamiento estructural de las vigas de sección variable. Resultados experimentales y de teoría de elasticidad han mostrado que las hipótesis usadas para el cálculo de esfuerzos de flexión en perfiles no prismáticos, se cumplen siempre y cuando las variaciones de las fronteras superiores o inferiores no sean muy altas, en contraparte, los resultados de la formula de cortante son muy engañosos, ya que se sabe que la influencia de este varia de forma considerable con la altura en términos de la llamada distorsión de la sección o alabeo. Se advierte que la formula de flexión en vigas no prismáticas es solo un medio aproximado para el cálculo de esfuerzos de flexión, sin embargo, teóricamente puede determinarse la altura o la sección de una viga completamente esforzada, usando la clásica formula de flexión bajo un nivel permisible de esfuerzo σadm, así:

S=

M

σ adm

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

(1)

164

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

P L/2

L/2 h ho x

Figura 1. Viga empleada para la aproximación del problema de optimización de peso del muelle

A partir de esto es posible deducir el modulo de sección S como una función de x y de esta manera también obtener la altura variable de la sección si se deja el ancho constante. Dado que todas las secciones estarán sometidas de manera aproximada al esfuerzo admisible se dice que la viga está totalmente esforzada. Se debe poner también especial atención al cortante pero para esto se requieren técnicas numéricas avanzadas tanto para el análisis estructural, como para la optimización de la sección.

2. APROXIMACIÓN DEL PROBLEMA DE MUELLES Cuando se quiere estudiar el comportamiento de un resorte de hojas puede considerarse una viga simplemente apoyada totalmente cargada que soporta una fuerza concentrada en el centro de su luz. La viga de la Figura 1 tiene sección rectangular de ancho constante b y esfuerzo admisible σadm.

El momento interno de la viga de la Figura 1, expresado en función de la posición, 0≤x0. Empleando el segundo teorema del clásico método de área de momentos puede determinarse que la deflexión δ para el centro de la luz es:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

166

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

PL3 δ= 24EI 0

(8)

Si se compara esta formula con la de una viga de sección constante con momento de inercia I0, puede verse que esta es el doble de la de una viga de sección constante sometida a la misma carga P. Es decir la viga tiene la misma resistencia pero diferente rigidez a la de una viga de sección constante (la rigidez se reduce a la mitad). El volumen optimizado puede obtenerse por integración de la curva de altura h multiplicado por el valor de la base constante b.

Vol óptimo =

2bh0 L 3

(9)

3. MARCO TEÓRICO DE LA SOLUCIÓN NUMÉRICA Antes de empezar a describir los resultados es importante tener en cuenta algunas definiciones básicas que introducen el lenguaje propio del diseño de optimización. Definiciones Básicas Antes de describir los resultados del diseño de optimización, es necesario definir alguna de la terminología básica empleada en este tipo de análisis, como son, variable de diseño, variable de estado, función objetivo, diseño factible y diseño no factible, archivo de análisis, iteraciones, ciclos, conjunto de diseño, etc. El problema planteado sirve como un ejemplo representativo de optimización estructural. Variables de Diseño (DVs) son cantidades independientes que son iteradas en orden para alcanzar el diseño óptimo. Se especifican límites superiores e inferiores como "restricciones" en las variables de diseño. Estos límites definen el rango de variación de las Variables de Diseño (DV). En el problema de la viga, las alturas h2,h3,h4,…,h10, son candidatos obvios para definir las DVs. Las alturas h2,h3,h4,…,h10, no pueden ser inferiores a 1.2 cm que es primer valor asignado a h1, así, sus más bajos límites serían h2,h3,h4,…,h10>1.2 cm. También, tienen un límite superior de hmax = 15.0. Hasta 60 DVs pueden definirse en ANSYS para un problema de diseño de optimización.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

167

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

Variables de Estado (SVs) son cantidades que condicionan o restringen el diseño. Se conocen también como "variables dependientes", y son típicamente variables de respuesta que son funciones de las variables de diseño. Una variable de estado puede tener un límite máximo y un límite mínimo, o puede tener un “único límite”. El problema de la viga posee como variable de estado el esfuerzo máximo en la viga Smax. En total se hace un mapeo de nueve variables de estado SVs: Smax2,Smax3,Smax4,…,Smax10 acotadas superiormente e inferiormente con una variación respecto al esfuerzo admisible de ± 250 kgf/cm2. En el ANSYS se pueden definir hasta 100 SVs en un problema de diseño de optimización. La Función Objetivo es la variable dependiente que se intenta minimizar. Debe ser una función de las DVs, es decir, si cambian los valores de las DVs debe cambiar el valor de la función objetivo. En el problema tratado, el peso es la función objetivo (que se minimiza). En ANSYS solo puede ser definida una sola función objetivo. Las variables de diseño, las variables de estado, y la función objetivo son llamadas colectivamente variables de optimización. En una optimización realizada en ANSYS el usuario representa estas variables a través de parámetros. El usuario debe identificar en el modelo cuáles parámetros son representativos de las DVs, SVs, y de la función objetivo. Un conjunto de diseño (o diseño) es simplemente un conjunto único de valores de parámetros que representan una configuración particular de un modelo. Típicamente, un conjunto de diseño se caracteriza por los valores de las variables de optimización; sin embargo, todos los parámetros del modelo (incluso aquellos que no están identificadas como variables de optimización), se incluyen en el conjunto. Un diseño factible es uno que satisface bien todas las restricciones o condiciones especificadas tanto para las SVs así como para las DVs. Si una o ninguna de las restricciones no se satisface, el diseño es considerado como no factible. El mejor diseño es aquel que satisface todas las restricciones y produce el mínimo valor en la función objetivo. (Si todo el conjunto de diseño es no factible, el mejor conjunto de diseño es aquel que más se acerca a la factibilidad, independiente del valor que produzca en la función objetivo). El archivo de análisis es en ANSYS el archivo de entrada (el cual puede ser creado por varios caminos) que contiene una secuencia completa de análisis (preproceso, solución, postproceso). Este archivo debe contener un modelo definido paramétricamente usando parámetros que representen todas las entradas y las salidas que serán usadas como DVs, SVs, y la función objetivo.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

168

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

Figura 2. Flujo de datos de Optimización

Desde este archivo, se crea automáticamente un ciclo de archivos de optimización (Jobname.LOOP) que es usado por el optimizador para ejecutar o desarrollar los ciclos de análisis. Un ciclo es un paso a través del ciclo de análisis. (Esto es como un paso por el archivo de análisis). La salida para el último ciclo ejecutado se preserva en un archivo llamado Jobname.OPO. Una iteración de optimización (o simplemente iteración) es uno o más ciclos de análisis los cuáles resultan en un nuevo conjunto de diseño. Típicamente, una iteración equivale a un ciclo. Sin embargo, para el método de primer orden, una iteración representa más de un ciclo. La base de datos de optimización contiene el último entorno de optimización, e incluye las definiciones de las variables de optimización, los parámetros, todas las especificaciones de optimización y los conjuntos de diseño acumulados. Esta base de datos puede ser salvada (por Jobname.OPT) o resumida a cualquier tiempo en el optimizador. Algunos de los conceptos descritos anteriormente pueden ser entendidos mejor a través de una ilustración. La Figura 2 muestra el flujo de información durante un análisis de optimización. El archivo de análisis debe existir como una entidad separada, dado que la base de datos de optimización no hace parte de la base de datos del modelo de ANSYS. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

169

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

Figura 3. Elemento BEAM3 (Salida de resultados)

Elemento BEAM3 El elemento BEAM3 puede ser usado como un elemento 2D para la modelación de problemas de Esfuerzo-Deformación en Vigas, cuenta con dos nodos y 3 grados de libertad por nodo. La salida de resultados puede ser a través de fuerzas y momentos internos o a través de esfuerzos normales por flexión, como se ilustra en la Figura 3.

Método de solución Para la optimización se utiliza el método de primer orden. El número máximo de iteraciones es 32.

4. PLANTEAMIENTO DEL PROBLEMA, SOLUCIÓN Y DISCUSIÓN Optimización de sección de una viga completamente esforzada, simplemente apoyada y sometida a carga puntual en el centro de la luz (aproximación del problema del resorte de hojas metálicas)

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

170

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1

Y Z

X

OPTIMIZACION MASA DE ESTRUCTURA

Figura 4. Idealización por simetría de la mitad de la longitud original del muelle (carga y condiciones de contorno)

Considerando una viga de longitud L = 120 cm, ancho b=6 cm, peso especifico γ = 7.85*10-3 kgf/cm3, carga P= 6000 Kgf situada en el centro de la luz, σadm=5000 kgf/cm2, módulo elástico E=2*106 kgf/cm2 y altura para el extremo apoyado h1=1.2 cm, es posible optimizar la sección de la viga considerando la condición de viga completamente esforzada y para esto en el análisis y diseño se debe determinar la variación de la altura h a lo largo de la longitud y el valor de h y δ en el centro de la luz. Dado que la viga posee simetría no es necesario modelar toda la viga sino solo media viga empotrada en un extremo y con carga P/2=3000 Kgf, por tanto, es posible establecer el modelo simplificado de la Figura 4 obtenido en ANSYS. La discretización se realiza dividiendo la viga en 10 partes que permitan el cálculo aproximado de alturas de la viga para estos 10 tramos, sin embargo, el primer tramo es fijo y se le asigna un valor de espesor de 1.2 cm considerando que este cumple los requerimientos de cortante del extremo. En la figura 5 se muestra la discretización:

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

171

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1

Y Z

X1

2

3

4

5

6

7

8

9

10

OPTIMIZACION MASA DE ESTRUCTURA

Figura 5. Discretización del modelo

1

NODAL SOLUTION STEP=1 SUB =1 TIME=1 UY (AVG) RSYS=0 PowerGraphics EFACET=1 AVRES=Mat DMX =1.936 SMN =-1.936

Y Z

X

MX

MN

ZV =1 DIST=33 XF =30 YF =-1.5 Z-BUFFER -1.936 -1.721 -1.506 -1.291 -1.076 -.86048 -.64536 -.43024 -.21512 0

OPTIMIZACION MASA DE ESTRUCTURA

Figura 6. Isocontornos de desplazamiento en dirección y

La Figura 6 muestra los resultados de desplazamiento del muelle y se aprecia para el extremo un desplazamiento de 1.936 cm.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

172

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1

LINE STRESS STEP=1 SUB =1 TIME=1 VOLU VOLU MIN =43.2 ELEM=1 MAX =214.982 ELEM=10 ZV =1 DIST=33 XF =30 YF =-1.5 Z-BUFFER 43.2 62.287 81.374 100.461 119.548 138.635 157.722 176.809 195.895 214.982

Y Z

X

OPTIMIZACION MASA DE ESTRUCTURA

Figura 7. Volumen optimizado por tramos

La Figura 7 muestra la solución de volumen optimizado para la discretización de diez tramos, en ella se aprecia la transición suave esperada (comportamiento parabólico de la altura h). La curva de la figura 8 muestra una convergencia acelerada en las primeras 6 o 7 iteraciones y en las posteriores un refinamiento del valor del peso hacia el óptimo (11.56 kgf). También puede verse que la convergencia se da sin oscilaciones. La familia de curvas de la figura 9 muestra un comportamiento de la convergencia de las variables de diseño muy similar al de la función objetivo.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

173

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1 OPTIMIZATION

PLOT NO.

PESO PESO

1

48 44 40 36 32 28 24 20 16 12 8 1

7.4 4.2

13.8 10.6

20.2 17

26.6 23.4

33 29.8

OPTIMIZACION MASA DE ESTRUCTURA

Figura 8. Evolución de la función objetivo (peso del muelle)

Tabla 1. Comparación de resultados Propiedad

Sol teórica

Variación en la curva de convergencia Deflexión del extremo (cm) Volumen Optimo (cm2) H1 H2 H3 H4 H5 H6 H7 H8 H9 H10

---

M. Primer Orden (BEAM3) Suave

2.0

1.936

2880

2946.08

1.342 2.324 3.000 3.550 4.025 4.450 4.837 5.196 5.532 5.848

1.200 2.400 3.091 3.647 4.129 4.559 4.952 5.315 5.652 5.972

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

174

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

1 OPTIMIZATION H2 H3 H4 H5 H6 H7 H8 H9 H10 H2 H3 H4 H5 H6 H7 H8 H9 H10 ESPESORES H

PLOT NO.

2

22 20 18 16 14 12 10 8 6 4 2 1

7.4 4.2

13.8 10.6

20.2 17

26.6 23.4

33 29.8

NUMERO DE ITERACIONES

OPTIMIZACION MASA DE ESTRUCTURA

Figura 9. Convergencia de las variables de diseño (Alturas de los tramos H)

Comparación de resultados entre la solución numérica y la solución teórica. En la tabla 1 se resumen los resultados obtenidos a partir de las idealizaciones realizadas con elementos tipo viga (BEAM3). Es importante notar que existe un buen grado de aproximación a la solución teórica de las variables de diseño y de la función objetivo, así como también de la deflexión calculada.

CONCLUSIONES El diseño de Optimización representa una herramienta de competitividad extremadamente eficiente para las empresas dedicadas al diseño de maquinaria e infraestructura, dada la gran variedad de posibilidades de optimización flexibles y adaptables por medio de la parametrización de modelos. Sin embargo, hay que tener cuidado con fenómenos numéricos que afectan la solución de manera significativa al punto de no predecir el comportamiento real de las estructuras.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

175

Rodriguez, W., Pallares, M., Diseño óptimo de estructuras: caso de vigas totalmente esforzadas

Es importante que cuando se resuelva un problema por elementos finitos o cualquier otra técnica numérica, tener en cuenta que siempre estas conducen a soluciones aproximadas que deben ser evaluadas con criterio ingenieril.

REFERENCIAS 1. Ansys User Manual / Revision 5.0. Swanson Analysis Systems, Inc. Houston. Volume I (Procedures), Volume Ii (Comands), Volume Iii (Elements), Volume Iv (Theory). 2. COSMOS/M 2.6. Structural Research and Analysis Corporation. (SRAC), OPTSTAR (OPTIMIZATION). 3. Oñate Ibañez de Navarra, Eugenio, Cálculo de Estructuras por el Método de Elementos Finitos, CIMNE. Barcelona , (1995). 4. Zienkiewicz, O. C., Taylor, R. L. El Método de los Elementos Finitos, Edit. Mc. Graw - Hill. Barcelona, (1994).

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

176

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Modelación y Chequeo Estructural del Puente Kantinwrwa Sobre el Rio Ariguani Chang, Gustavo Facultad de Ingeniería - Departamento de Ingeniería Civil Unviersidad del Magdalena, Colombia [email protected]

Resumen En este documento se analizan las causas que generaron el colapso del tablero del puente sobre el río Ariguani, ubicado en las estribaciones de la Sierra Nevada entre los departamentos del Magdalena y Cesar. Inicialmente se realizó una visita técnica al sitio del siniestro con el fin de inspeccionar de forma directa el estado, forma y posición de la estructura caída, durante esta visita se tomaron la mayor cantidad de imágenes fotográficas para poder establecer las reales condiciones con las que se construyó el puente. Posteriormente se desarrolló un trabajo de oficina en la cual se elaboró un modelo matemático que nos permitió determinar la capacidad estructural de todos los elementos del sistema. Con toda esta información ordenada, procesada y analizada se establecieron las causas que ocasionaron la falla estructural, se chequearon las capacidades de los elementos y se dieron las recomendaciones a seguir en la reconstrucción. Palabras claves : Puente, chequeo estructural, modelación matemática.

1. INTRODUCCION El Gobierno Nacional a través de ACCION SOCIAL y de entidades como PRO SIERRA, se encuentran desarrollando proyectos de infraestructura física en zonas deprimidas con alto índice de pobreza, presencia de cultivos ilícitos y de grupos armados irregulares, con el fin de mejorar las condiciones socioeconómicas de la población colombiana. En este caso el gobierno nacional ha proyectado la construcción de nueve pueblos indígenas en las estribaciones de la Sierra Nevada de Santa Marta que beneficiaran a las cuatro comunidades Arhuaco, Kogui, Wiwa y Kankuamo.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

177

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Figura 1. Modelo computacional del puente.

Una de las obras de infraestructura corresponde al puente colgante denominado “KANTINWRWA”, ubicado sobre el río Ariguani el cual es el límite natural entre los departamentos de Magdalena y Cesar, el cual sufrió el pasado mes de abril-2009 la caída del sistema de piso (cables, pendolones, vigas metálicas, placas de concreto y lámina alfajor). El objetivo fundamental de este estudio es el de determinar las causas que llevaron a tal desenlace, por lo tanto el día 01 de mayo de 2009 se desarrolló una visita técnica al lugar del siniestro, en la cual se hizo una inspección visual y se tomaron distintas imágenes fotográficas a cada una de las partes importantes de la estructura colapsada con el fin de poder establecer las variables y las reales condiciones con las que se construyó el puente colgante. Posterior a la actividad de campo se consultaron diversas fuentes de información tales como: planos topográficos, memorias de cálculo y planos de construcción. Con toda esta información organizada, analiza se procedió a la elaboración de un modelo matemático tridimensional lo más real posible en donde se simulara su Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

178

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

geometría, masa y/o peso propio, rigidez, cargas externas (carga viva, fuerza de viento, fuerza sísmico, empujes de tierra entre otras.) y fue analizado a través del programa de análisis estructural SAP 2000 versión 14. El método empleado por el programa se basa en el método de los elementos finitos y más exactamente en el análisis matricial. Con los resultados obtenidos de esta modelación pudimos comprobar si los índices de sobre esfuerzo o flexibilidad de todos los elementos superaban el cien por ciento permitido. Por último se analizaron las fotografías tomadas, planos de construcción, los resultados de la modelación numérica y se procedió a dar las conclusiones y recomendaciones a seguir en el proceso de reconstrucción del puente.

2. METODOLOGIA 2.1 Geometria Del Modelo De los planos de topográficos y estructurales suministrados por el Arquitecto Daniel Varón se extrajeron las dimensiones que permitieron armar el modelo matemático 3D. En el apéndice A se relacionan las coordenadas empleadas a partir de las cuales se obtuvo el modelo 3D mostrado en la figura 1.

2.2 Evaluación De Cargas 2.2.1 Carga Muerta. Para las cargas muertas se tuvo en cuenta que el sistema de piso está conformado por una lámina de alfajor de 4 mm de espesor y una sobre placa de concreto reforzado de 40 mm de espesor y que fueron instaladas para darle mayor estabilidad al tablero del puente debido a su aumento del peso por carga muerta. Consultando los catálogos de los productos tenemos:

Elemento Lámina Alfajor Loseta concreto Total carga muerta

Peso (kgf/m2) 33.00 96.00 129.00

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

179

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Figura 2. Aplicación de cargas en el modelo.

El peso propio de las cerchas, pendolones, cables y de concreto como columnas, vigas y zapatas los incluye automáticamente el software. La carga muerta se distribuye hacia las cerchas de la siguiente forma: 113 kgf/m para la cercha central y de 56.5 kgf/m para las cerchas laterales.

2.2.2 Carga Viva. De acuerdo con el Código Colombiano de Diseño Sísmico de Puentes de 1995 y la AASHTO-2004 (American Association of State Highway and Transportation Officials), la carga viva para uso peatonal corresponde a 400 kgf/m2. Esta carga viva se distribuye hacia las cerchas de la siguiente forma: 350 kgf/m para la cercha central y de 175 kgf/m para las cerchas laterales. 2.2.3 Carga de Viento. Según el código de puentes para Vw = 160 KPH

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

180

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Figura 3. Cargas de barlovento (izquierda). Cargas de sotavento (derecha).

CARGA MÍNIMA DE VIENTO Barlovento= 450 kgf/m Sotavento= 225 kgf/m PARA 130KPH Barlovento: CSUP=CINF = Sotavento: CSUP=CINF =

450 x (130/160)² = 225 x (130/160)² =

297 kgf/m 149 kgf/m

Se tomaron estos valores pues la superficie de exposición para una presión de 370 Kgf/m2 no superaba la anterior carga por metro lineal. Las fuerzas de viento aplicadas al modelo son mostradas en la figura 3.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

181

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Espectro Elástico de Diseño

0.300 0.250 Sa (%g)

0.200 0.150 0.100 0.050 0.000 0.00

1.00

2.00

3.00

4.00

5.00

T (s)

Figura 4. Espectro elástico de diseño.

2.2.4 Movimientos Sísmicos de Diseños. La ubicación geográfica de el proyecto nos indica que este se encuentra ubicado en una zona de amenaza sísmica intermedia, el perfil estratigráfico del suelo clasifica se estima como un suelo tipo S-2: − − −

Art. A.2.3 Zonas de menaza sísmica Aa = 0.10 Art. A.2.4.1 Tipo de perfil de suelo S2, coeficiente de sitio S=1.2 Art. A2.5.1 Coeficiente de importancia, clasifica dentro de las estructuras de grupo de uso I, coeficiente de importancia I=1.0

Con base en la anterior información el espectro elástico de diseño tendría una forma como la que muestra la en figura 4:

2.2.5 Capacidad de Disipación de Energía. Se empleó un coeficiente de disipación de energía de R = 4.0 para este tipo de elementos metálicos.

3. COMBINACIONES DE CARGAS Las diversas cargas (viva, muerta, viento, fuerzas inerciales) deben combinarse con el objetivo de poder determinar las máximas solicitaciones sobre los elementos estructurales y estos puedan diseñarse de la manera más conveniente.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

182

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

La tabla 1 muestra la combinación de carga utilizada y los resultados obtenidos en esfuerzos para cada caso. Tabla 1. Combinación de carga utilizada y resultados obtenidos Combo #

4.

TIPO DE CARGA PP 1.4

CM 1.4

ENVCV

1 2

1.2

1.2

1.6

3 4

1.2 1.2

1.2 1.2

0.5 0.5

5 6 7 8

1.2 1.2 1.2 1.2

1.2 1.2 1.2 1.2

9 10 11 12

1.2 1.2 1.2 1.2

13 14

W

Sxx

Syy

0.5 0.5 0.5 0.5

1.0 -1.0 1.0 -1.0

0.3 0.3 -0.3 -0.3

1.2 1.2 1.2 1.2

0.5 0.5 0.5 0.5

0.3 -0.3 0.3 -0.3

1.0 1.0 -1.0 -1.0

1.2 1.2

1.2 1.2

0.5 0.5

15 16

0.9 0.9

0.9 0.9

1.3 -1.3

17 18

0.9 0.9

0.9 0.9

1.3 -1.3

Szz

Tmin

Tmáx

1.3 -1.3

1.0 -1.0 1.20 1.20 1.2 1.2

CHEQUEO DE ELEMENTOS

4.1 Chequeo de los Cables La máxima fuerza generada por carga muerta y viva que son los dos tipos de carga que mas generan fuerza axial en los cables fue de 35.75 toneladas (ver figura 5). Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

183

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Figura 5. Izquierda: resultados obtenidos en los cables. Derecha: Chequeo de las fuerzas en los pendolones

Los resultados del ensayo de tracción realizados por la empresa EMCOCABLES muestran que la probeta arrojo un resultado de 93.678 lbf que equivalen a 42.58 toneladas. Por lo tanto los cables están en capacidad de soportar una carga viva de 400 kgf/m2 y de la carga muerta sobre impuesta correspondiente de la losetas de concreto reforzado de 40 mm de espesor. Ya en el análisis de los pendolones Los niveles de sobreesfuerzo dan menor que la unidad (ver figura 5). Respecto a las columnas de concreto, de acuerdo con todas las combinaciones de carga se observa que estos elementos estructurales deben tener un área de acero mínima de 18.00 cm2. En los planos de construcción hechos por el diseñador original recomienda instalar 8 barras de 1.00 pulgada es decir un área de 40.48 cm2. Se observa que las columnas están sobradas. La figura 6 muestra los resultados obtenidos del modelo. Por otro lado, en el chequo de las platinas de anclaje Se asume que la fuerza axial que viaja por el cable será de 42.5 toneladas que corresponden a la máxima fuerza de tensión que resiste el cable. La forma en la que este elemento trabaja es a flexión y cortante (ver figura 6). Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

184

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

Figura 6. Izquierda: chequeo de las columnas de concreto. Derecha: Chequeo de las Platinas de Anclaje.

De los cálculos realizados la platina de soporte debería tener un espesor mínimo de 6.80 cm y la calidad del acero de dicha platina debe A36.

CONCLUSIONES Y RECOMENDACIONES − El puente colgante “KANTINWRWA” sobre el rio Ariguani, está la capacidad de en capacidad de soportar el peso propio de sus elementos constitutivos y el de las sobrecargas muertas compuestas por las losetas de concreto de 40 mm de espesor y el de una carga viva de servicio de 400 kgf/m2, que es la carga mínima recomendada por los códigos vigentes. ¡Es importante destacar la necesidad de limitar la carga viva a un valor de máximo de 480kgf/m2, con el fin de evitar llegar a la máxima carga de rotura de los cables.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

185

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

−

−

Que después de realizar el chequeo individual de los elementos más importantes de la estructura tales como cables, pendolones, cerchas y columnas estos pueden soportar las cargas y esfuerzos que les sean aplicados durante su vida útil. Es necesario que se cambie la platina que sirve para transmitir la fuerza de los cables a las dos barras ancladas en el muerto, actualmente los planos originales establecen que esta platina debe tener un espesor de 9.00 cm, la que está actualmente instalada es de ¾ pulgada es decir 1.90 cm y la que recomiendo debe tener como mínimo 6.80 cm.

−

Dado que el constructor no siguió las recomendaciones establecidas por el diseñador original en lo referente del anclaje de los cables al muerto, y decidió instalar 3 (tres) grapas (o perros) en cada extremo de cada cable. Es necesario se instalen como mínimo 4 (cuatro) grapas en cada extremo, esto basado en la “NORMA CHILENA DE EMERGENCIA OFICIAL NCH885. EOf72”.

−

Es también importante instalar contratuerca en cada grapa o perro y se debe apretar con el torque requerido por el fabricante de este tipo de elemento, controlar el nivel de grasa en las zonas de contacto de las grapas y el cable.

−

En el extremo del cable es imprescindible se instale algún dispositivo que permita ensanchar la punta del cable para que en un posible y futuro deslizamiento del cable por causa de la falla de las grapas solo descienda el puente y no colapse como así ocurrió. El ensanchamiento de la punta del cable no permitirá que este pase a través del ojal.

−

En cuanto a la prueba de carga que se debe realizar a la estructura reconstruida y de esta forma corroborar la capacidad de la estructura, se deben emplear sacos llenos de material pétreo de la zona que simulen la carga viva de 400 kgf/m2. El tablero del puente se debe dividir en 4 (cuatro) zonas; por lo tanto deberán demarcarse de la zona 1 a la zona 4, como sigue: 3 1

4 2

La forma en la que se debe colocar la carga es como se aprecia en las siguientes secuencias: •

Cargar solo la zona 1.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

186

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

• • • • • • • •

Cargar solo la zona 2. Cargar solo la zona 3. Cargar solo la zona 4. Cargar solo las zonas 1 y 2. Cargar solo las zonas 1 y 3. Cargar solo las zonas 1 y 4. Cargar solo las zonas 2 y 3. Cargar solo las zonas 2 y 4.

REFERENCIAS 1. Asociación Colombiana de Ingeniería Sísmica. Codigo Colombiano de Diseño Sismico de Puentes.1995. 2. American Association of State Highway and Transportation Officials. AASHTO.2004. 3. Wai-Fah Chen, Lian Duan. Bridge Engineering Handbook.2000. 4. Norma Chilena de Emergencia Oficial NCH885. EOf72

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

187

Chang, G., Modelación y chequeo estructural del puente Katinwrwa sobre el rio Ariguani

ANEXO A. Type POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT POINT

Name 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

X (m) -12 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 52

Y (m) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Z (m) 0.00 0.00 0.05 0.09 0.14 0.18 0.23 0.27 0.32 0.36 0.41 0.46 0.50 0.50 0.46 0.41 0.36 0.32 0.27 0.23 0.18 0.14 0.09 0.05 0.00 0.00

Nota: las coordenadas en el eje Z, representan la altura relativa que hay entre la cota del tablero y el inicio del puente que se toma como nivel 0.00.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

188

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Caracterización de Compuestos Laminados a Partir de la Simulación Numérica 3D a Nivel Micromecánico y Mesomecánico Edgardo W. Arrieta Ortíz Departamento de Ingeniería Mecánica Universidad Tecnológica de Bolívar [email protected]

Resumen La deformación bajo carga de laminados de fibra de vidrio para uso estructural en embarcaciones de la Armada Colombiana es evaluada mediante una serie de modelos computacionales por elementos finitos desde el nivel micromecánico incluyendo aproximaciones mesomecánicas. Bajo condiciones de laboratorio, se fabrica y prueba diferentes configuraciones de laminados para validar los resultados. Ante la ausencia de un modelo completo, que prediga el comportamiento mecánico de los laminados, han sido desarrolladas aproximaciones numéricas para simular estos materiales. También ha sido propuesta una relación constitutiva general algorítmica para caracterizar laminados. Sin embargo, para el diseño no hay a la mano un modelo de uso claro que prediga dicho comportamiento. Se prueba extensivamente modelos por elementos finitos simulando desde el nivel micromecánico el comportamiento del material y se compara los resultados contra pruebas de laboratorio. Palabras Clave: Mecánica de medios contínuos, Elasticidad, Métodos Numéricos, Elementos Finitos, Materiales Compuestos.

1. INTRODUCCION Las cada día mayores aplicaciones de los materiales compuestos, hacen no solo interesante sino muy productivo el estudio de los mismos. Construir herramientas de diseño que permitan a los fabricantes actuales agilizar sus procesos de manufactura usando materiales compuestos de alta eficiencia y bajo peso se convierte en un reto que exige de los más avanzadas técnicas de análisis teórico, experimental y computacional. Puesto que no existe un modelo completo de uso generalizado que permita predecir el comportamiento mecánico de los compuestos laminados de fibra de vidrio nos proponemos avanzar en la caracterización de Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

189

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

dichos laminados usando un modelo computacional basado en elementos finitos 3D aplicado a la micromecánica de una capa de material laminado usando Fibra de vidrio en una resina de poliéster. Modelos mixtos que requieren la calibración computacional de parámetros para una aproximación basada en micromecánica teórica han sido propuestos [1] para caracterizar compuestos basados en laminados de nanotubos de fibra de carbono. Principios físicos simples son aplicados a la determinación de las propiedades mecánicas de un laminado de fibra de vidrio experimentalmente y luego computacionalmente. La correlación entre los valores predichos por la teoría de mezclas [2] y los arrojados por la experimentación y por el modelo computacional, se muestra en los resultados.

2. MATERIALES, METODOS Y MEDICIONES Se utilzó fibra de vidrio provista por el astillero COTECMAR para realizar pruebas concernientes a sus propiedades fisicas y mecánicas. Usamos el tramado Woven Roving con densidad 800gr/m², al igual que el Matt 225 gr/m². Para las simulaciones numéricas se utilizó el software de código libre CODEASTER 9.4, ampliamente usado en universidades y centros de investigación en la Unión Europea, desarrollado por Electricitè de Fránce, y el generador de mallas Gmsh[3], también código libre. 2.1 Evaluación del diámetro de las fibras El análisis digital de las microfotografías tomadas con un microscopio óptico Leitz Wetzlar, permitió determinar que el diámetro promedio de las fibras usadas es de aproximadamente 20mm, En la Figura 1, el diámetro de la zona visible fue determinado y calibrado para el aumento de 100X, resultando de una medida de 1.2 mm. Las mediciones fueron tomadas en pixels sobre las imagenes digitalizadas, y corregidas por el ángulo de inclinación de la fibra con la horizontal: Si α es el ángulo de la fibra con el eje horizontal, d el diámetro de la fibra, y dh el ancho de la fibra medido en sentido horizontal, entonces: d = dhsen(α)

(1)

La evaluación digital de los diámetros de las fibras arrojó d = 23.61 mm, con desviación estándard de 3.18 mm es decir del 13.49%.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

190

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Figura 6. Medición del diámetro de las Fibras. Izquierda: Haz de fibras. Centro: Fibra individual. Derecha: fibras entrelazadas.

2.2 Evaluación del la densidad de las fibras Procedimos a evaluar el peso de un haz de fibras en balanzas de precisión y determinar el volumen desplazado usando un recipiente calibrado con agua, según muestran las gráficas. Los datos obtenidos mostraron que la densidad de las fibras es en promedio ρ = 2674Kg/m³, con desviación estándard del 12.38%.

2.3 Resistencia última y módulo de Young de las Fibras La resistencia última de las fibras, la matriz de poliéster usada, y de un laminado consistente en dos capas contiguas (una de Woven Roving 800 y la siguiente de Matt225), obtenidos en laboratorio, bajo temperatura y humedad fijas, fueron medidas en el laboratorio de ensayos de la UTB (Universidad Tecnológica de Bolivar), usando una máquina universal de 60 Ton (600KN), para un total de 15 probetas.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

191

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Habiendo determinado la densidad en peso por unidad de área del tramado WR800, y conociendo la densidad de haces de fibra por unidad de longitud podemos obtener el peso de cada haz de fibras por unidad de longitud. De igual modo, conociendo la densidad de las fibras y el peso por unidad de longitud de las mismas obtuvimos el area seccional promedio de un haz de dichas fibras. Si As es el área seccional de un haz de fibras, y se observa N haces de fibras horizontales y N verticales en un tramado de WR800, que debe pesar (800gr), cada uno de los cuales mide una longitud L de 1m, y notando que el peso de cada haz M, será el peso de la trama dividido el número de haces 2N, usando la densidad ρ: M = ρAsL

(2)

M = ρAs (1m) = 800g/(2N)

(3)

Notando que en el tramado, N es igual a 10 haces cada 5 cm es decir 200 haces por metro, entonces, hallamos As = 7.48x10-⁷ m². Aún mas midiendo el peso de un haz de fibras de longitud 1.45 m, que resultó ser de 3.23gr en cinco mediciones distintas, usando la Ec.(2), obtuvimos As con un valor de 8.33x10-7m². Usando el valor medido para el diámetro de las fibras (23.61 µm), esto nos conduce a aproximadamente 1900 fibras por haz. Verificamos, que usando un empaquetamiento hexagonal de fibras de sección circular, por geometría simple, el área ocupada es aproximadamente (1/η) veces el área de las fibras, con η igual a π/(4sen(60)). De este modo seis haces de fibras tendrán un diámetro visible aproximado de:

4 ⋅ (6A s ) ηπ

(4)

  π  = 0.906 η =   4 ⋅ sen(60) 

(5)

dH = con:

Usando el valor de As previamente obtenido, notamos que dH, debe ser de aproximadamente 2.6 mm, como en efecto se observa para un paquete de seis haces de fibras, salvo las diferencias entre un empaquetamiento hexagonal perfecto versus el arreglo que presentan las fibras al ser juntadas en el laboratorio. Procedimos entonces a evaluar el módulo de elasticidad y la resistencia a la tracción de los haces de fibra de vidrio.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

192

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Figura 7. Medición de Resistencia a la tensión y Módulo de Young. Izquierda: Montaje para grupos de seis. Centro: detalle de montaje. Derecha: Detalle montaje de fibras en maquina de ensayo.

Notando que los cabezotes de la máquina de tracción, marca SHIMADZU, con capacidad para 60Ton, prácticamente no se deforman durante los ensayos de los haces de fibras, podemos tomar la lectura de desplazamiento de los mismos (con precisión de 0.001mm) contra la curva de carga (con precisión de 0.01KN es decir de 10N o aprox 1.0kg-f), para obtener la constante de resorte de los haces de fibras considerados como barra o cable rígido a tensión. Siendo bien conocido que para una barra a tensión de área seccional A, módulo de Young E, y longitud L, la constante equivalente de resorte es:

k=

A⋅ E L

(6)

Podemos resolver entonces el Módulo de Young para la fibra de vidrio que se está usando. La tabla siguiente resume este proceso, habiendo descartado algunos ensayos por deslizamiento de las mordazas:

E=

L⋅k 6 ⋅ As

(7)

Arrojando entonces un valor del módulo de Young para la fibra de 37.5 GPa, con desviación de aproximadamente el 4%. Por otro lado, la evaluación de la carga máxima en cada ensayo, nos permite determinar el máximo esfuerzo de resistencia a la tensión de las fibras, resultando (Tabla 2.) Sut = 783 Mpa con deviación del 15% aprox. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

193

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Tabla 1. Cálculo Módulo de Young y Resistencia a la Tensión Fibra de Vidrio Probeta No. K (10⁵⁵ N/m) Fmax(KN) Sut(Mpa) E (Gpa) Datos 2 7,77 4,25 850,34 38,87 As= 8,33E-07 4 7,50 3,25 650,26 37,52 L= 0,25 6 7,22 4,25 850,34 36,11 7,496666667 3,916667 783,6467921 37,49833 Promedio 0,57735 115,5162603 1,375626 Desv.Est 0,275015151 3,67% 14,74% 14,74% 3,67%

2.4 Resistencia última y módulo de Young de los laminados Se preparó en laboratorio un laminado consistente de dos capas contiguas (una de Woven Roving 800 y la siguiente de Matt 225), bajo temperatura y humedad fijas, con relaciones de peso 45% resina, 55% fibra de vidrio. En el laboratorio de ensayos de materiales de la UTB se sometieron a tracción probetas rectangulares de 12 mm de ancho, 25 cm de largo y aprox. 3 mm de espesor, recortadas a 0º respecto a la dirección del WR. Despreciando la no linealidad de las curvas, y considerando el promedio del comportamiento lineal del laminado como producto de un material uniforme con módulo de elasticidad equivalente E, obtenemos que: Tabla 2. Modulo de Young y resistencia última de lamiando K (N/m) e (m) h (m) L (m) Fmáx (N) 615000 0,00115 0,01530 0,2130 2100 548000 0,00117 0,01570 0,2050 3125 581500 0,00116 0,01550 0,2090 2612,5 Promedio Desv 47376 0,00001 0,00028 0,0057 724,784451 8,15% 1,22% 1,82% 2,71% 27,74% Ensayo 1 2

de Fibra de Vidrio Sut (Mpa) 119,4 170,1 144,737833 35,9008642 24,80%

De modo que el laminado con 45% en peso de Matriz y 55% en peso de Fibra, con dos capas, Matt225 y WR800 posee un módulo de Young equivalente a 0° de 6.8 Gpa, con desviación del 14% aproximadamente.

2.5 Modelos usando Elementos Finitos Se desarrolló una serie de modelos 3D usando elementos finitos para determinar el módulo equivalente del laminado con fibras unidireccional, simulando fibras de 20 micras de diámetro inmersas en la matriz de poliéster con el módulo indicado por la literatura (3Gpa) [2], y dando a las fibras el módulo de 39 Gpa obtenido de la experimentación. Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

194

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Figura 8. Modelo Micromecánico 3D por Elementos Finitos. Izquierda: Sección transversal de la malla. Derecha: Vista en perspectiva de las Fibras

La malla obtenida usando Gmsh [3], disponiendo n fibras horizontalmente en dos hileras según el arreglo hexagonal. El bloque simulado mide 1mm de largo (relación de largo/diámetro = 50 para las fibras) y tiene dimensiones mínimas para alojar las n x m fibras en su interior. Se sometió a tensión el bloque aplicando 500Mpa de tensión en la dirección longitudinal.

3. RESULTADOS Es conocido que para un material compuesto con dos componentes dispuestos axialmente de modo uniforme [2] el módulo equivalente en términos de las proporciones volumétricas de fibra y matriz se puede expresar como:

E mf = α ∗ E f + (1 − α ) ∗ E m

(8)

donde,

E mf = Módulo Equivalente del laminado, E f , E m = Módulos de Young de fibra y matriz,

α, (1 − α ) = Fracciones volumétricas de fibra y matriz, v m , v f = Volumenes ocupados por fibra y matriz, ρ m , ρ f , ρ mf = Densidades de fibra, matriz y laminado,

β,(1 − β ) = Fracciones de masa de fibra y matriz,

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

195

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Tabla 4. Modulo Equivalente predicho por la Ec.(8), en Gpa ρfibra Efibra 2674 37,5 ρmatriz 1200 Ematriz 3 Capa Mfibra M resina β α Emf Matt225 54,2 115,35 0,31966971 0,1741 9,00793765 WR800 143 115,35 0,55351268 0,3575 15,3325707 2capas 197,2 230,7 0,46085534 0,2772 12,5650566 promedio 12,1702542 WR800* 71,5 115,35 0,38265989 0,2176 10,5082582 Matt225* 0,361333 115,35 0,00312271 0,0014 3,04843054

Lo anterior asumiendo que se genera un estado de esfuerzos uniaxial. El anterior modelo puede ser expresado en términos de las fracciones en peso, notando que:

ρ mf = αρ f + (1 − α )ρ m

(9)

y también

 ρmf β  ρ  f

  = α,  

 ρmf  ρm

(1 − β )

  = (1 − α ) 

(10)

Así entonces:

 ρ mf E mf = β   ρ  f

 ρ  ∗ E f + (1 − β ) mf  ρ   m 

 Em  

(11)

También resulta útil, escribir:

α=

1  1 − β   ρ f 1 +  ⋅   β   ρ  m 

   

(12)

La validez de las anteriores ecuaciones estará limitada por el hecho que en realidad, al poseer módulos distintos la deformación de las interfaces curvas entre fibra y matriz no producen una deformación geométricamente proporcional en el plano transversal, es decir las formas geométricas en el plano x-y de nuestro modelo son claramente distorsionadas, tal como era de esperarse.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

196

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

Tabla 5. Modulo Equivalente por Elementos Finitos versus el predicho por la Teoría Micromecánica. N2000

N4000

dL1

540,5405 775,51 957,8313 765,06 543,0108 774,19

dL2

17,12 23,83 33,07 25,13 18,41 24,52

def_Unit

SigmaEq Eequiv Alfa FEM 0,028549 500 17,51388 0,34 0,041168 500 12,14536 0,34 0,026402 500 18,93792 0,34

Eteor 14,73076 14,73076 14,73076

Al no tomar este tipo de suposiciones, esperamos que el modelo 3D arrojará valores más cercanos a los reales que el modelo micromecánico dado por la Ec. (11). Cada una de estas dos capas posee un modulo equivalente, y cada una tiene una proporción de fibra por unidad de peso total. Modelaremos cada una con el modelo micromecánico dado por la Ec.(8) y luego el total asumiremos que se comporta como un resorte equivalente a dos resortes en paralelo. El peso de material usado fue: 54.2gr de Matt, y 143gr de WR, (ambos en un área de 58x35 cm²), y se uso en total 230.7 gr de resina de poliéster. Con estos datos: Podemos observar que el modulo equivalente predicho por la teoría micromecánica es mucho mayor que el medido experimentalmente (6.78GPa). Esto puede ser debido a factores como la presencia de vacíos y la disparidad entre la suposición de deformación unidimensional y la real. Sin embargo, considerando que cada capa de laminado produce un material equivalente con constante equivalente propia, en el caso de tener capas del mismo espesor, asumiendo igual deformación en el sentido de la carga, el modulo equivalente para las dos capas juntas debe ser el promedio de los módulos de cada una de las capas, como se deduce usando la analogía que permite sumar resortes en paralelo como resistencias eléctricas en serie. Vemos que el modulo equivalente predicho por la teoría para las dos capas con base en su proporción en peso de fibra y resina (12,57GPa), difiere del módulo promedio de las dos capas (12,17 Gpa). Así pues según estas suposiciones mesomecánicas tal como en la Ec.(8):

Eeq = 0.5 ( E1 + E2 )

(13)

Adicionalmente, del material Matt225 no podemos a priori especificar qué proporción resiste carga, pero sí podemos calcularla de modo que el módulo Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

197

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

equivalente corresponda al obtenido en laboratorio. Así pues las dos ultimas líneas de la Tabla 4, corresponden a WR800 en el cual solo la mitad de masa de fibras (multiplicamos β or 0.5) actúa, y sólo 1/150 del material del Matt225 resistiría carga en la dirección axial. De los modelos generados por elementos finitos seleccionamos dos nodos en el interior ubicados entre 0.5 y 0.9 mm de la longitud total de 1mm del modelo. Se determina la deformación en z para dichos nodos y se evalúa la deformación unitaria. El material fue sujeto a 500 Mpa a tensión usando el software CodeAster[4]. Determinamos el valor del modulo de Young Equivalente basado en esta medida de deformación unitaria, variando un factor de escala en la configuración de la malla. Se determina la relación de área seccional de fibras dividida por el área seccional de la “probeta” en el modelo de elementos finitos (alfa en la Ec.(8)), para evaluar el Modulo equivalente usando la Ec.(8).

4. CONCLUSIONES De acuerdo con la Tabla 5, se observa que el modelo de elementos finitos arroja un modulo equivalente promedio de 16,2Gpa, para fibras con E=37,5Gpa, y matriz con E=3Gpa, con factor volumétrico del 34% para las fibras, lo cual representa el 10% de error respecto al predicho por la teoría correspondiente al modelo de mezcla unidimensional. Tal como lo muestra la Tabla 4, este valor corresponde muy bien a la condición de la capa de Woven Roving800 usada en la experimentación. Esto quiere decir que no es necesario asumir que solo el 50% del material del WR actúa en el ensayo a tensión, como sugerimos antes, y sugiere una interacción más pronunciada de las fibras transversales del tramado en la deformación del material. Simular esto último en modelos 3D, al igual que determinar las respuestas del material Matt225 con orientación aleatoria requiere de mayor capacidad computacional a la usada anteriormente. Aunque fue ensayada una colección extensa de modelos variando parámetros geométricos, para determinar rangos de convergencia, los modelos seleccionados indican claramente la habilidad de esta metodología para aportar a la evaluación y correlación de las propiedades de estos laminados para uso naval.

REFERENCIAS [1] Gang Huang et. Al, Material Characterization and Modeling of Single-Wall carbon Nanotube/Polyelectrolyte Multilayer Nanocomposites, Journal of Applied Mechanics, Sept. 2006, Vol. 73, pp. 737-744.

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

198

Arrieta, E., Caracterización de Compuestos Laminados de Fibra de Vidrio para Uso Naval…

[2] J.M. Berthelot, Composite Materials, Springer-Verlag, 1999. Págs. 10,17. [3] C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Meth. Eng., Accepted for publication, 2009. [4] Christophe Durand, Free software for computational Mechanics: EDF's choice, http://www.code-aster.org/V2/UPLOAD/DOC/Presentation/2007_nafems.pdf .

Segundo Simposio Regional en Mecánica No-Lineal de Estructuras Continuas – SMEC 2009

199