Validation of a numeric approach for modelling delamination damage caused by low energy impact

1 Introduction

Este paper recoge el desarrollo de un modelos numerico para predecir el comportamiento de materiales compuesto frente a impactos de baja energía. El modelo trata de predecir el comportamiento elástico así como el daño interlaminar producido a consecuencia del impacto.

Este tipo de modelos puede tener aplicación para optimizar de manera virtual el diseño de un producto fabricado en material compuesto sin tener que recurrir a prototipos funcionales que incrementarían tanto el tiempo como el coste de desarrollo del producto.

2 Descripción de los test de validación

2.1 Material

Para validar el modelo se han considerado 3 tipos de laminado fabricados a partir de infusión con 6, 10 y 14 láminas de epoxy (SiComin’s SR 8100/SD 8822) respectivamente reforzados con un tejido de fibra de carbono (Hexcel’s Ref. 43199). Las placas de laminados tienen un diámetro de 60 mm y un espesor de 1,4 mm 2,3 mmy 3,2 mm respectivamente.

2.2 Propiedades de la lamina

La tabla 1 recoje las propiedades mecánicas de la lamina siendo las direcciones 1 y 2 las principales en el plano del laminadoy la 3 la normal al propio plano [1].

Tabla 1. Propiedades de la lámina
E1 [GPa]	E2 [GPa]	E3 [GPa]	G12 [GPa]	G23 [GPa]	G31 [GPa]
57	57	2	3.1	2.9	2.9
ν ¸'© Ù?¸'© Ù?	ν ¸'© Ù?¸'© Ù?	ν ¸'© Ù?¸'© Ù?	ρ [kg/m3]	X1 & X2 [MPa]	X13& X23 [MPa]
0.08	0.08	0.08	1430	784	81

Los módulos elásticos en el plano E₁ y E₃ (ASTM D3039)_, el módulo de cortadura fuera de plano (G₂₃ y G_31, ASTM D 3518), la densidad (ρ¸'© Ù? ASTM 792¸'© Ù? y la resistencia en el en plano (X₁ y X₂)y a cortadura X₁₃ y X₂₃ fueron ensayados in-situ.

La fracción de volumen de fibra del material caracterizado y el porcentaje de porosidad fueron caracterizados in-situ y correponden a 49% (ASTM D3171) y 3.2% (ASTM D2734) respectivamente.

2.3 Propiedades interlaminares

Para simular el daño por delaminación se considera la energía de fractura interlaminar en modo I y II. La energía de fractura en modo I fue caracterizada a través de un ensayo DCB (Double Cantilever Beam) conforme a la norma ASTM D5528-13 [2]. La energía de fractura en modo II fue caracterizada a través de un ensayo ENF (End Notched Fracture) de 4 puntos [3].

La tabla 2 muestra los valores obtenidos y empleados en este proceso de validación.

Tabla 2. Propiedades de la lámina [2-3]
GIC [J/m2]	GIIC [J/m2]
500.0	3300.0

2.4 Caracterización de impacto del material

Los ensayos de impacto con los que se ha validado el ensayo han sido caracterizados a través de una máquina de impacto de dardo Fractovius plus con un impactador de cabeza hemiesférica de 20 mm de diámetro y una masa de 2 kg. Está equipada con un modulo que evita las multiples colisiones posteriores debidos al rebote entre el impactador y la probeta.

Gracias al transductor de fuerza equipado sobre el impactador se registra el historial de la fuerza de contacto.El especimen se apoya sobre un anillo anular con un diámetro interno y externo de 40 mm y 60 mm respectivamente.

Para los tres tipos de láminado se han considerado 2 niveles de energía distintos. Un primer nivel de energía que denominaremos subcrítica con una forma sinusoidal caracteristica fuerza-tiempo (sobre el que no se ha podido identificar daño alguno tras el impacto) y un segundo nivel de energía que denominaremos supercrítica de delaminación con una forma de curva fuerza-tiempo que presenta una primera fase linear y una fase posterior en el que aparece cierta oscilación fruto del daño interno provocado por el laminado.

La figura 1 muestra el patron de ambos tipos de curva fuerza –tiempo.

La tabla 3 muestra el nivel de energía de impacto empleado en cada caso.

Tabla 3. Propiedades de la lámina
	Nº de láminas	6	10	14
Sub-crítico	Energía [J]	1,0	1,0	2,0
	V impactador [m/s]	1,00	1,00	1,41
Super-crítico	Energía [J]	4,0	12,0	8,0
	V impactador [m/s]	2,00	3,46	4,47

3 Modélo numérico

Para caracterizar el dañó por delaminación se ha empleado LS-Dyna V971 un código basado en el método de elementos finitos inicialmente desarrollado para trabajar con métodos de integración temporal explicito. Para post-procesar los resultados se ha utilizado el módulo de post-procesado LS_PREPOST V4.3.

(a)

(b)

(b)

Figura 1. Curva caracteristica fuerza de contacto frente a tiempo: laminado con a) 6 láminas, b) 10 láminas y c) 12 láminas.

3.1 Modelo geometrico

La figura 2 muestra el modelo geométrico empleado. El modelo considera el soporte anular sobre el que descansa la placa, la cabeza del impactador y el laminado que se caracteriza lámina a lámina.

(a)

(b)

Figura 2. Modelo geometrico y mallado. a) aspecto general y b) detalle del laminado.

1.1 Caracterización del daño por delaminación

Para caracterizar el daño por delaminación se han considerado dos tipos de aproximación. Por un lado, se ha considerado el daño a través de la aproximación CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK con el parámetro OPTION 9-Discrete Crack Model cuyo criterio de delminación se basa en la energía a fractura de modo I y modo II. La aproximación considera las leyes Power Law o Benzeggagh y Kenane para el modo mixto de fractura.

Por otro, se considera el daño a través de un elemento sólido cohesivo de espesor cero intercalado entre elementos de láminas adyacentes a través de la ley de material MAT_COHESIVE_MIXED_MODEL. La aproximación que considera este comando comporta una formulación de ley cohesiva admitiendo también las leyes de Power Law o Benzegah Kenane para el modo mixto de fractura. En este último caso, también es preciso aplicar el comando ERODING_SINGLE_SURFACE para que, en caso de que el elemento cohesivo se abra totalmente, se genere de manera automática una ley de contacto entre láminas adyacentes que evite que se penetren entre sí.

3.2 Mallado

Todos los modelos se basan en un mallado con elementos de tipo sólido con integración reducida (ELFORM=1: constant stress solid element).

Para generar la malla se ha utilizado el mallador de Ansys Workbench^® y su extensión a LS-Dyna^®. En el caso del impactador se ha aplicado un tamaño medio de malla de 2,5 mm y para el caso del soporte anular se ha utilizado un elemento de malla de 4 mm.

El laminado se ha mallado lámina a lámina poniendo un único elemento por espesor. Cada una de las láminas se ha dividido en dos zonas, mallando aquella zona susceptible de ser delaminada con un tamaño de elemento medio de 0.5 mm y aquella otra zona fuera de dicha zona a través de un tamaño de elemento medio de 1.5 mm (ver figura 2).

El mallado del laminado se ha realizado fuera de Ansys^® a través de un procedimiento desarrollado internamente en Matlab^® para este cometido. El algoritmo desarrollado parte de un mallado 2D de la cara superior y la extruye capa a capa. En el caso de caracterizar la cohesión entre láminas a través de elementos de contacto, el mallado extruye láminas sin conectar los nodos de una lámina con la de la lámina inmediatamente superior (ver figura 3). El caso de emplear elementos cohesivos, los mallados de dos láminas adyacentes son conectados por el elemento cohesivo (ver figura 3).

Figura 3. Esquema de generación del mallado en el caso de emplear elementos cohesivos y elementos de contacto.

3.3 Materiales

Para el caso del impactador y del soporte anular se ha utilizado una ley de material isotrópica elástica y lineal definiéndola a través de MAT001-ELASTIC. La tabla 4 presenta las propiedades introducidas.

Tabla 4 Propiedades del acero
E [GPa]	ν	ρ [kg/m3]
210	0.3	7850

Puesto que únicamente se considera la cabeza del indentador, es preciso incorporar masa adicional para que la cabeza adquiera la inercia total del impactador (2 kg). Este aumento de masa se ha incorporado a través del comando MASS_PART en lugar de incrementar artificialmente la densidad de la cabeza. De esta manera se consigue que parámetros como el contacto, etc. no se vean afectados por el incremento artificial de la densidad de la cabeza.

En el caso del laminado ,se emplea la formulación MAT002-ORTHOTROPIC ELASTIC con las propiedades que se recogen en la tabla 2. Las direcciones principales del material se introducen a través de la opción AOPT=2.

En lo que respecta a la delaminación, tanto en la aproximación basada en el contacto como en la de la ley cohesiva se utilizan las propiedades consideradas en la tabla 2. En ambos casos se aplica la ley de Power Law con el exponente 1.21.

3.4 Contactos

Aparte del ya mencionado contacto entre láminas se han definido dos contactos usando AUTOMATIC_SURFACE_TO_SURFACE. Se ha definid un primer contacto entre el soporte anular (contacto maestro) y la lámina inferior (contacto esclavo) del laminado y otro entre el laminado (contacto maestro) e impactador (contacto esclavo).

En ambos casos se ha empleado un coeficiente de rozamiento tanto estático como dinámico de 0.2.

3.5 Condiciones de contorno

A todos los cuerpos se les aplica aceleración debido a la gravedad a través del comando BODY_Z.

La caracterización del impactador desde su altura total requeriría un tiempo de computación muy alto. En su lugar se considera el instante t=0 de la simulación aquel en el que el impactador entra en contacto con el laminado. En este caso, la energía acumulada por el impactador se aplica en forma de velocidad, de acuerdo a los valores recogidos en la tabla 3.

Tanto la aceleración por gravedad como la velocidad inicial se aplican durante la fase de relajación dinámica y se mantienen durante la simulación transitoria posterior (véase apartado 3.5).

3.6 Relajación dinámica

Antes de simular el transitorio se inicializa la aceleración de gravedad de manera que el laminado se encuentre en el instante t=0 s en equilibrio de fuerzas considerando todas las fuerzas actuantes. Con este objetivo se aplica la relajación dinámica. En este caso se ha realizado una relajación dinámica empleando un aproximación explícita a través del comando DYNAMIC_RELAXATION con IDRFLG=1.

3.7 Amortiguamiento

El modelo considera un amortiguamiento global de un 2% respecto del amortiguamiento crítico a través del comando GLOBAL_DAMPING.

3.8 Control de hourglass

Al emplear elementos de integración reducida, se ha aplicado una ley que minimiza el modo de deformación espúrea de elementos de integración reducida conocido como hourglass o efecto de reloj de arena. En este caso, al tratarse de elementos sólidos, se ha optado por el comando HOURGLASS con la opción IHQ =4 y un valor de QH = 0.1.

En todas las simulaciones se ha optado por un tiempo de simulación de 0.004 s. El tiempo de integración temporal inicial es de 1.38.10^-7 s y, posteriormente, se deja al propio código modificarlo de manera automática según criterios internos relacionados con el tamaño mínimo de elemento. En este cálculo el criterio basado en el tamaño crítico se ha minorado 0,7 gracias a la variable TSSFA del comando TIMESTEP.

4 Resultados

4.1 Validación del modelo para el caso subcrítico

En la siguiente figura se comparan los resultados obtenidos a partir de los modelos numérico con los ensayos presentados en el apartado 2 para el caso subcrítico. Para obtener el historial de fuerza frente al tiempo entre el impactador y el laminado se emplea el comando RCFORCE con una resolución de 1.10-6 s.

Ajustado el coeficiente de penalti del contacto, los resultados son practicamente identicos para el caso de elementos cohesivos y elementos de contacto entre láminas en todos los casos subcríticos. En ningún momento se aprecia un daño significativo por delaminación en ninguno de los laminados modelados.

1.1 Inestabilidad en delaminación por contacto

A diferencia del subcrítico, con el supercrítico, en los modelos considerados se produce una inestabilidad sobre los elementos de contacto. Como consecuencia de esta inestabilidad, se detectan una serie de delaminaciones puntuales distribuidos sin un patrón claro sobre la zona de contacto entre el impactador y el laminado (véase figura 5).

Este comportamiento aparentemente es muy similar al descrito por [4].

Tras revisar la literatura, para tratar de mitigar dicho efecto se han realizado los siguientes ajustes:

• Cambiar la formulación de contacto entre el impactador y la primera capa a la opción SOFT =1 y SOFT=2

• Reducir el valor de penalti normal y tangencial de contacto bajo la opción de formulación SOFT=0, SOFT=1 y SOFT=2.

• Incorporar un amortiguamiento local en el contacto contacto a través del parámetro VDC.

Con ninguna de las acciones anteriores se ha llegado a reducir de manera significativa este comportamiento inestable del contacto.

El modelo basado en elemento cohesivo proporciona un comportamiento más coherente con los resultados experimentales, sin la presencia de inestabilidad alguna, tal y como muestra la figura 6. Es por ello que se decidió continuar el análisis únicamente con el modelo basado en elementos cohesivos.

1.1 Caso supercrítico

La figura 7 muestra la relación fuerza de contacto frente al tiempo para el caso del laminado de 10 láminas bajo un impacto de 12 J para el caso experimental y su predicción numérica.

(a)

(b)

(c)

Figura 4. Curva caracteristica fuerza de contacto frente a tiempo: laminado con a) 6 láminas, b) 10 láminas y c) 12 láminas.

Figura 5. Rotura por delaminación entre la lamina superior y la lámina inmediatamente inferior (rojo indica la apertura total de la ley cohesiva).

Figura 6. Indice de delamiación entre la lamina superior y la lámina inmediatamente inferior (rojo indica la apertura total del elemento cohesivo).

Figura 7. Curva caracteristica fuerza de contacto frente a tiempo: experimental frente a predicción numérica para el laminado de 10 láminas e impacto de 12 J.

El resultado muestra que el modelo desarrollado predice de manera precisa el efecto de la delaminación en el comportamiento macroscópico del laminado, en especial, en la fase inicial de la delaminación. Predice también la oscilación caracteristica posterior al comienzo de la delaminación. No obstante el comportamiento es algo distintos en la segunda parte de la curva. Esto podría ser debido a que únicamente se está considerando el fallo por delaminación mientras que en el experimental se puede producir daño intralaminar.

La figura 8 muestra la zona delaminada ( zona sin elementos) a lo largo de todas las capas del laminado. La tabla 5 muestra la comparación entre el área dañada experimental y su predicción.

Figura 8. Rotura por delaminación entre la lamina superior y la lámina inmediatamente inferior (rojo indica la apertura total de la ley cohesiva).

Tabla 1. Área delaminada media entre láminas
	Numérico [mm2]	Experimental [mm2]
10 láminas-12 J	191	235

Conclusiones

A continuación se enumeran las principales conclusiones obtenidas de este estudio:

• Se ha desarrollado un modelo numérico basado en elementos finitos que permite caracterizar el daño por delaminación en laminados sujetos a impacto de baja energía.

• El modelo ha sido validado para el caso de energías subcríticas para laminados de 6, 10 y 14 láminas.

• El modelo también ha sido validado para el caso de energías supercríticas para el laminado de 10 láminas. Para este caso, el modelo predice de manera precisa la influencia de la delaminación sobre el comportamiento macroscópico del laminado. Predice también de manera razonable el área de delaminación y las oscilaciones posteriores al comienzo de daño por delaminación.

A continuación se citan los siguientes pasos de que se pretenden llevar a cabo en el futuro relacionados con esta actividad:

• Extender la validación del modelo numérico a los laminados de 6 y 14 láminas para el caso supercrítico.

• Incorporar el daño intralaminar para tratar de mejorar el comportamiento tras el inicio de la delaminación.

• Validar la resistencia residual post impacto de laminados a través de la validación del comportamiento a tracción de probetas dañada.

Agradecimientos

Este trabajo de investigación ha sido desarrollado gracias al proyecto ACTIMAT encuadrado en el programa ELKARTEK del Gobierno Vasco.

Referencias

[1] H. Zabala, L. Aretxabaleta, G. Castillo, J. Urien, J. Aurrekoetxea, Impact velocity effect on the delamination of woven carbon–epoxy plates subjected to low-velocity equienergetic impact loads. Composites Science and Technology, 94, pág. 48–53 (2014).

[2] H. Zabala, L. Aretxabaleta, G. Castillo,J. Aurrekoetxea, Loading rate dependency on mode I interlaminar fracture toughness of unidirectional and woven carbon fibre epoxy composites, Composite Structures 121, pág. 75–82 (2015).

[3] H. Zabala, L. Aretxabaleta, G. Castillo,J. Aurrekoetxea. Dynamic 4 ENF test for a strain rate dependent mode II interlaminar fracture toughness characterization of unidirectional carbon fibre epoxy composites, Polymer Testing, 55, pág. 212–218 (2016).

[4] Lstc Composite models. URL: http://www.dynasupport.com/howtos/material/composite-models. (2017-05-01).