Novel Discrete Multiscale Model (DM) for non-linear analysis in composite materials

1 Introducción

El aumento de complejidad de la microestructura interna de los materiales compuestos hace que en los últimos años los métodos numéricos basados en elementos finitos representen una solución al problema de como analizar a nivel microscópico el comportamiento de estos materiales. De hecho, la sustitución de análisis experimentales por simulaciones virtuales permite obtener resultados comparables de forma más económica. Para el análisis del comportamiento de estos materiales se han desarrollado varias técnicas como la Simulación Numérica Directa (DNS) que utiliza el Método de Elementos Finitos a una escala para modelar la microestructura. Este tipo de simulación representa el nivel más alto de precisión, pero debido a la discretización explícita de la microestructura, la simulación de estructuras complejas de grandes dimensiones tiene costes computacionales prohibitivos. En el caso particular de los materiales compuestos, se desarrollaron formulaciones simplificadas con las que reducir el coste del análisis, entre ellas está la homogenización fenomenológica, o teoría de mezcla clásica [1], donde la hipótesis fundamental es que para todos los componentes del compuesto la deformación sea la misma. Otro tipo de homogeneización fenomenológica es la teoría de mezclas Serial/Paralelo (SP) [2,3], que considera la distribución topológica de cada uno de los componentes que componen el compuesto. Oller [4] y Otero [5] extendieron esta teoría a los compuestos con matriz reforzada y nanotubos. Otras técnicas permiten de analizar a nivel microscópico el comportamiento de los materiales que componen una estructura son los Métodos Multiescala [6,7]. Estos métodos se basan en la solución de varios análisis anidados, cada uno en una escala diferente, para obtener la respuesta global del sistema que se quiere analizar. Estos métodos han llegado a acoplar escalas atómicas con escalas estructurales [8].

En el caso de estructura de materiales compuestos, el comportamiento de la microestructura se define como la solución de un Elemento de Volumen Representativo (RVE) o microestructura del material compuesto. En la mayoría de los casos, esto consiste en obtener las tensiones en la macro-escala a partir de un modelo numérico micro del material compuesto, donde las condiciones de borde se definen desde la macro. Sin embargo, a pesar del aumento de las capacidades computacionales actuales, realizar un análisis multiescala de una estructura real, teniendo en cuenta las no linealidades materiales, sigue siendo una tarea muy cara en términos de coste computacional [9]. Varios modelos de reducción se han desarrollado para la optimización de ese tipo de análisis como Response Surface Method (RSM) [10], Proper Ortogonal Decomposition (POD) [11] y otros [12] [13].

Analizando el nivel macro, podemos observar que en los casos más comunes la fractura empieza en una zona muy pequeña de la estructura. Esto significa que la mayor parte del dominio sigue en régimen lineal-elástico y no requiere el cálculo de un RVE. En esta línea, el modelo propuesto por Otero et al. [14] se basa en la definición de una función de activación no lineal (NLAF). Los autores definen una función de umbral, basada en la energía de deformación requerida por el RVE para alcanzar su tensión de fallo. El procedimiento requiere un primer análisis de los diferentes RVE de la estructura, bajo las diferentes cargas que se les aplican. El coste computacional de esta operación se reduce gracias a un algoritmo Smart First Step (SFS). Posteriormente, el método sólo analiza los RVE que alcanzan la non-linealidad. La principal ventaja de este método es que permite realizar el análisis no lineal de la estructura, con un coste computacional razonable. En [14] se demuestra que puede reducir el coste computacional del análisis en más de un 98% en comparación con un procedimiento Multiescala Completo. El problema de este método surge cuando queremos cambiar la estructura macro, o las cargas aplicadas, utilizando el mismo RVE. En ese caso el método no puede utilizar ningún dato de los calculados anteriormente y se debe determinar nuevamente la función de activación no lineal, con el consiguiente coste computacional.

Otra metodología para la optimización de este tipo de análisis es el procedimiento de simulación a ruptura. En este caso, la macroestructura se simula linealmente utilizando un tensor de rigidez obtenido a partir de un análisis inicial del RVE. Una vez que ha terminado la macro-simulación, se pueden localizar los elementos con campos de esfuerzos de deformación más grandes y se puede analizar el RVE en esas condiciones. Este análisis proporcionará la carga en la que se produce el fallo. La principal ventaja de este enfoque es que en la mayoría de los campos de ingeniería el interés principal es obtener la carga de fallo, y no el rendimiento de la estructura en campo no-lineal.

El propósito de este trabajo es proporcionar un procedimiento capaz de predecir la respuesta mecánica no lineal de los materiales compuestos con el fin de reducir el coste computacional necesario para el análisis. A las formulaciones desarrolladas les llamaremos Modelos Multiescala Discretos (Discrete Multiscale, DM). Éstos utilizan un Data Base (DB) calculado a priori que define el comportamiento del RVE. A través de la definición de un parámetro de daño equivalente ( ${\textstyle {d}_{eq}}$), función de la tensión global resultante desde la microescala, se realizan una serie de ensayos virtuales del RVE, controlados por el campo de deformaciones, y los resultados se almacenan en el DB. Los resultados a guardar son el estado de tensiones y deformación homogenizados alcanzados en ciertos niveles de ${\textstyle {d}_{eq}}$. Finalmente, la solución de la estructura macroescala a través del Método Multiescala de Primer Orden (FE2) será reemplazada por la interpolación de los datos almacenados en el DB. El método propuesto, el Discrete Multiscale Threshold Surface (DMTS) obtiene el estado tenso-deformacional para el que el material daña. En este momento el modelo DM es sustituido por un RVE.

En este trabajo se ha prestado especial atención a la creación y validación de la base de datos. Posteriormente se describe el modelo DMTS y su aplicación al estudio de una estructura completa comparando los tiempos de cálculo y los resultados con los que se obtendrían con un método FE2.

2 Formulación Multiescala de primer orden

Esta sección resume los principales conceptos y ecuaciones de la formulación de Multiescala de primer orden [14]. Ésta se basa en la obtención de la respuesta de la estructura macroscópica a partir del análisis de un elemento de volumen representativo (RVE) que define el material a nivel de microescala.

En un procedimiento de múltiple escala, la deformación obtenida al analizar la estructura macroscópica se utiliza para definir las condiciones de contorno en el problema en la microescala. La solución del problema en la microescala, en tales condiciones, actúa como una ley constitutiva equivalente para la macroescala, ya que proporciona parámetros materiales como la rigidez y las tensiones. Esta ley constitutiva equivalente se utiliza en todos los puntos de Gauss del modelo macroscópico para obtener la respuesta global de la estructura. En caso de un comportamiento no lineal del modelo microscópico, el RVE debe resolverse para diferentes condiciones de contorno hasta que ambas escalas alcancen el equilibrio, asegurando la consistencia entre las soluciones de micro y macroescala.

La técnica de homogeneización de primer orden utilizada asume una separación de escala de la micro a la macro-escala, donde la longitud característica de la microescala será mucho menor que la longitud de los elementos macroscópicos. El procedimiento supone pequeños desplazamientos y un comportamiento estructural cuasi estático en ambas escalas, macro y micro.

El algoritmo propuesto puede definirse en tres fases principales. El primer paso es la transición desde macro a micro, donde la deformación macroscópica en cualquier punto de la malla macroscópica se transfiere a la microescala. Posteriormente, el problema de condición de frontera (BVP) se resuelve en la microescala para obtener su respuesta mecánica. Finalmente, los parámetros de macroescala se obtienen a través de la transición desde micro a macro (Figura 1). Se utilizará el subíndice, para distinguir las dos escalas, μ para el micro dominio (RVE) y m para el macro dominio.

En el nivel macro, el punto de partida para un método de homogeneización computacional basado en la cinemática es la suposición de que el tensor de deformación mecánica, ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{m}}$, en cada punto del dominio en la macroescala, ${\textstyle {\Omega }_{m}}$ (donde la posición se define a través del vector ${\textstyle {\mathit {\boldsymbol {x}}}_{m}}$) Se puede obtener por un cierto instante t como el promedio de volumen del campo de deformación mecánica microscópica, ${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{\mu }}$, definido en cada punto del dominio de la microescala, ${\textstyle {\Omega }_{\mu }}$ (donde la posición se define a través del vector ${\textstyle {\mathit {\boldsymbol {x}}}_{\mu }}$) y en el mismo instante t como:

${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{m}\left({\mathit {\boldsymbol {x}}}_{m},t\right)=}$${\frac {1}{{V}_{\mu }}}\int _{{\Omega }_{\mu }}^{}{\mathit {\boldsymbol {\epsilon }}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu },t\right)dV$

(1)

El campo microscópico de deformación puede expresarse como el gradiente simétrico del campo de desplazamiento microscópico, ${\textstyle {\mathit {\boldsymbol {u}}}_{\mu }=}$$({{u}^{x}}_{\mu },{{u}^{y}}_{\mu })$:

${\textstyle {\mathit {\boldsymbol {\epsilon }}}_{m}\left({\mathit {\boldsymbol {x}}}_{m},t\right)=}$${\frac {1}{{V}_{\mu }}}\int _{{\Omega }_{\mu }}^{}{{\nabla }^{s}{\mathit {\boldsymbol {u}}}}_{\mu }dV$

(2)

Sin pérdida de generalidad, podemos descomponer el desplazamiento en la microescala como:

${\textstyle {\mathit {\boldsymbol {u}}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu },t\right)=}$${\mathit {\boldsymbol {\epsilon }}}_{m}\left({\mathit {\boldsymbol {x}}}_{m},t\right)\cdot \left({\mathit {\boldsymbol {x}}}_{\mu }-\right.$$\left.{\mathit {\boldsymbol {x}}}_{g}\right)+{\overline {\mathit {\boldsymbol {u}}}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu },t\right)$

(3)

Donde ${\textstyle {\overline {\mathit {\boldsymbol {u}}}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu },t\right)}$ son los desplazamientos y sus fluctuaciones con respecto a los campos medios y ${\textstyle {\mathit {\boldsymbol {x}}}_{g}}$ (el centroide de la microescala) en cada instante ${\textstyle t}$.

Finalmente, también tenemos que definir las condiciones de contorno correspondientes en la microescala. Para los campos de fluctuación de desplazamiento, aplicamos condiciones de contorno periódicas, ya que generalmente proporcionan una respuesta intermedia y más exacta en comparación con otro tipo de condiciones de contorno, tal como se describe en [15] [16] [17] [18]:

${\textstyle {\overline {\mathit {\boldsymbol {u}}}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu }^{+},t\right)=}$${\overline {\mathit {\boldsymbol {u}}}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu }^{-},t\right)\,\forall \,pair\left\{{\mathit {\boldsymbol {x}}}_{\mu }^{+},{\mathit {\boldsymbol {x}}}_{\mu }^{-}\right\}\,\in \,{\partial \Omega }_{\mu }$

(4)

Siendo ${\textstyle {x}_{\mu }^{+}}$ y ${\textstyle {x}_{\mu }^{-}}$ dos puntos opuestos en el límite RVE ${\textstyle {\partial \Omega }_{\mu }}$.

Para el problema mecánico, el principio Hill-Mandel de macro-homogeneidad establece que la potencia de tensión macroscópica, producto del tensor de tensión macro, ${\textstyle {\mathit {\boldsymbol {\sigma }}}_{m}}$, y la tasa de deformación macro, ${\textstyle {\overline {\mathit {\boldsymbol {\epsilon }}}}_{m}}$, debe ser igual a la potencia de tensión microscópica sobre el promedio RVE en el dominio de volumen. Por lo tanto:

${\textstyle {\mathit {\boldsymbol {\sigma }}}_{m}:{\overline {\mathit {\boldsymbol {\epsilon }}}}_{m}=}$${\frac {1}{{V}_{\mu }}}\int _{{\Omega }_{\mu }}^{}{\mathit {\boldsymbol {\sigma }}}_{\mu }:{\overline {\mathit {\boldsymbol {\epsilon }}}}_{\mu }dV$

(5)

El tensor de esfuerzo macroscópico homogeneizado puede obtenerse como el campo de esfuerzo microscópico del RVE promediado en el volumen como:

${\textstyle {\mathit {\boldsymbol {\sigma }}}_{m}={\frac {1}{{V}_{\mu }}}\int _{{\Omega }_{\mu }}^{}{\mathit {\boldsymbol {\sigma }}}_{\mu }\left({\mathit {\boldsymbol {x}}}_{\mu },t\right)dV}$

(6)

La relación macroscópica constitutiva definida por las propiedades homogeneizadas del RVE puede obtenerse después de la solución de la microescala.

Suponiendo que el equilibrio de la microescala se expresa como:

${\textstyle \int _{{\Omega }_{\mu }}^{}{\mathit {\boldsymbol {\sigma }}}_{\mu }:{{\nabla }^{s}{\overline {\mathit {\boldsymbol {u}}}}}_{\mu }dV=}$$0$

(7)

Como se describe en [19] y [20], el tensor constitutivo homogeneizado ${\textstyle {\mathit {\boldsymbol {C}}}^{H}}$ puede definirse como:

${\textstyle {\mathit {\boldsymbol {C}}}^{H}={\frac {1}{{V}_{\mu }}}\int _{{\Omega }_{\mu }}^{}{\mathit {\boldsymbol {C}}}_{\mu }dV}$

(8)

Donde ${\textstyle {\mathit {\boldsymbol {C}}}_{\mu }}$ es el tensor constitutivo del RVE.

La evaluación del tensor constitutivo homogeneizado se realiza con un método de perturbación. Para cada columna j del tensor constitutivo se aplica al RVE una pequeña perturbación de deformación (δ ${\textstyle {\overline {\epsilon }}_{j}}$) para obtener, junto con la Ecuación 7, un tensor de tensión perturbado ( ${\textstyle {\delta }^{j}{\overline {\sigma }}}$). Las j columnas del tensor constitutivo homogeneizado pueden obtenerse como:

${\textstyle {\mathit {\boldsymbol {C}}}_{j}^{H}\equiv {\frac {{\delta }^{j}{\overline {\mathit {\boldsymbol {\sigma }}}}}{\delta {\overline {\mathit {\boldsymbol {\epsilon }}}}_{j}}}}$

(9)

Figura 1. Esquema de formulación multiescala de primer orden.

3 Implementacion de un Modelo Discreto Multiescala

Como se ha descrito antes, proponemos un método que representa la solución del problema de la macroescala y sustituye el FE² por un Modelo Discreto Multiescala (DM). Éste consiste en analizar el comportamiento del RVE en un número discreto de estados tenso-deformacionales para, posteriormente, utilizar los mismos para predecir el comportamiento del material. Un caso particular de un modelo discreto es el Discrete Multiscale Threshold Surface (DMTS) en el que la información discreta se utiliza para evaluar si el RVE ha alcanzado su comportamiento no lineal. Si lo ha alcanzado, el programa resolverá el micro-modelo RVE cada paso de tiempo. En la Figura 2 se muestra el procedimiento del método DMTS.

Figura 2. Algoritmo DMTS.

Tal como se ha comentado, un modelo discreto multiescala se basa en la realización de un estudio preliminar del RVE para diferentes casos de carga, como tracción, compresión, cortante, etc... Estos análisis, que denominaremos historial de deformaciones, pueden definirse en los espacios multidimensionales ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\gamma }_{xy})}$ para 2D o ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{zz},{\gamma }_{xy},{\gamma }_{yz},{\gamma }_{xz})}$ para problemas mecánicos en 3D. En este capítulo, por simplicidad, sólo describiremos el caso 2D, pero es perfectamente extensible a 3D.

La Base de Datos de tensiones homogeneizadas, cuya estructura será descrita en el siguiente capítulo, almacenará la información de la microescala para hacer accesibles estos datos durante el análisis de macroescala. Por razones obvias de las dimensiones de la base de datos y para obtener el comportamiento del RVE en un determinado estado de esfuerzo-deformación se introducirá un parámetro, denominado daño equivalente, que proporciona la relación del esfuerzo homogeneizado y el correspondiente elástico. Este parámetro se definirá como:

${\textstyle {d}_{eq}={\frac {\parallel {\sigma }_{el}-{\sigma }_{real}\parallel }{\parallel {\sigma }_{el}\parallel }}\quad ,\quad {d}_{eq}\in [0,1]}$

(10)

Donde ${\textstyle {\mathit {\boldsymbol {\sigma }}}_{real}}$ es la tensión homogeneizada, ${\textstyle {\mathit {\boldsymbol {\sigma }}}_{el}=}$${\mathit {\boldsymbol {C}}}_{el}^{H}:{\mathit {\boldsymbol {\epsilon }}}$ es la tensión elástica y ${\textstyle {\mathit {\boldsymbol {C}}}_{el}^{H}}$ es el Tensor Elástico Constitutivo Homogeneizado.

La discretización de la respuesta RVE con un parámetro único es necesaria para la interpolación de las informaciones alamcenadas en el DB en caso de direcciones de deformación diferentes de las calculadas anteriormente.

3.1 Discretización espacio de deformaciones

Para obtener una definición completa del comportamiento de la microescala, se impondrá como historial de deformaciones el valor del punto equispacio proyectado a la esfera de radio unitario centrada en el origen de los ejes ${\textstyle ({\epsilon }_{xx},{\epsilon }_{yy},{\gamma }_{xy})}$. Podemos definir de forma única estos puntos en el sistema de coordenadas esféricas 3D como la combinación de tres parámetros ${\textstyle \theta ,\varphi ,r}$, como podemos ver en la Figura 3. Donde los ángulos ${\textstyle \theta }$ y ${\textstyle \varphi }$ representan la dirección de las cargas de deformación que aplicamos a la microescala y ${\textstyle r}$ es la intensidad de deformación unitaria.

De coordenadas esféricas a cartesianas podemos reconstruir el vector de deformación como:

${\textstyle \left\{{\begin{matrix}{\epsilon }_{xx}=r\cdot cos(\theta )\\{\epsilon }_{yy}=r\cdot sin(\theta )\cdot cos(\varphi )\\{\gamma }_{xy}=r\cdot sin(\theta )\cdot sin(\varphi )\end{matrix}}\right.}$

(11)

${\textstyle r={\sqrt {{\epsilon }_{xx}^{2}+{\epsilon }_{yy}^{2}+{\gamma }_{xy}^{2}}}}$

(12)

De la Ecuación 11 podemos observar que la deformación es periódica con funciones ${\textstyle sen}$ y ${\textstyle cos}$ de los ángulos ${\textstyle \theta \in [0,\pi ]}$ y ${\textstyle \varphi \in [0,2\pi ]}$. La subdivisión de este intervalo determina el número de análisis a realizar y la precisión de la discretización.

Figura 3. Esfera de radio unitario en el espacio de deformaciones.

(b)

Figura 4. Evolución del número de análisis necesarias.

Introduciendo el parámetro ${\textstyle m\in [1,inf[}$ como subdivisión del intervalo ${\textstyle [0,2\pi ]}$ podemos definir de forma única la dirección de deformación con dos parámetros enteros ${\textstyle \in [-}$$m,m]$, que llamaremos tag. Con este método cada punto será unitario espaciado uniformemente del otro y los ángulos ${\textstyle \theta }$ y ${\textstyle \varphi }$ variarán entre ${\textstyle [{\frac {-\pi }{m}},{\frac {\pi }{m}}]}$. De esta manera, podemos observar que la cantidad total de análisis es ${\textstyle (2\cdot m+}$$1{)}^{2}$. Considerando la superposición de direcciones de deformación debido al uso de ${\textstyle [-}$$\pi ,\pi ]$ como intervalo para ambos ángulos, la cantidad única de historias de deformaciones se reducirá a ( ${\textstyle 2\cdot m}$)·( ${\textstyle m-}$$1$)+2. A pesar de esta reducción, la definición de un valor alto de m implica un aumento exponencial del coste de análisis, como se puede ver en la Figura 4 anterior.

3.2 Análisis de la microestructura y construcción de la Base de Datos

Las diferentes deformaciones definidas anteriormente serán aplicadas al RVE haciendo un Análisis Clásico de Multiescala de Primer Orden, como se describe en el Capítulo 1. Para simplificar utilizamos como geometría en la macroescala un triángulo con un punto de integración donde las deformaciones se imponen como desplazamiento a los nodos. De esta forma, todos los grados de libertad del triángulo serán fijos y para cada paso de tiempo se resolverá el problema de la microescala, evaluando el esfuerzo homogeneizado y el tensor constitutivo. Con estas informaciones podemos determinar el valor correspondiente de daño equivalente en cada paso de tiempo y guardar en el database el correspondiente estado tenso-deformacional. En este trabajo el ${\textstyle {d}_{lim}=}$${d}_{eq}=0.01$

En el caso de un análisis 2D, la siguiente figura puede proporcionar una mejor comprensión de cómo funciona el algoritmo propuesto con una representación gráfica.

Figura 5. Representación gráfica de la Base de Datos en el espacio de deformaciones.

La construcción de la Data Base (DB) se realiza aplicando sobre el modelo del RVE los distintos campos de deformaciones definidos por los tags. Éstos corresponden a los puntos A, B, C y D de la Figura 5. La superficie definida por estos puntos para ${\textstyle {d}_{lim}}$ corresponde al Threshold en el que empieza el daño. La superficie definida por los puntos A’, B’, C’, D’ correspondería a un daño equivalente mayor.

Cuando se realiza una simulación, la deformación que se aplique sobre el elemento no tiene por qué coincidir con ninguno de los tags analizados. Los puntos E y E' representan la intersección de la dirección de la tensión y la deformación con las superficies correspondientes distintos daños equivalentes.

El valor de tensión y deformación para el punto E de la Figura 5 representa el límite elástico y se determina interpolando los datos almacenados en el DB correspondientes a los tags cercanos. De esta forma podemos reconstruir la respuesta de deformación y tensión del RVE para cualquier dirección posible en el espacio de deformaciones.

4 Definición de la Discrete Multiscale Threshold Surfaces para distintos materiales

4.1 Material de daño a Tracción y Compresión

Para validar la construcción de la base de datos con un material simple analizamos un RVE con una ley constitutiva de daño en tensión y compresión, de acuerdo con la ley constitutiva implementada en [21].

Tabla 1. Damage Traction/Compression Properties
Young Modulus	900	MPa
Poisson Ratio	0.15
Stress Traction Limit	0.25	MPa
Traction Fracture Energy	0.016	J/mm2
Stress Compression Limit	10.5	MPa
Compression Fracture Energy	40.0	J/mm2

A continuación, se presenta la sección de la superficie límite elástica de deformación para ${\textstyle {{d}_{lim}=d}_{eq}=}$$0.01$.

Las superficies obtenidas con el modelo reducido propuesto en Figura 6 representan la deformación limite elástica del modelo constitutivo de daño. Posteriormente, será posible reproducir una superficie de fluencia de daño a través de una serie de puntos discretos obtenidos desde los análisis del RVE.

(a)	(b)
(c)	(d)

Figura 6. Superficie limite elástica en el espacio de deformaciones.

4.2 Material compuesto

El material compuesto estudiado se compone de una matriz de resina epoxi (F155) y fibra de carbono. La resina epoxy se simulará mediante la ley constitutiva utilizada en el apartado 4.1 y [21], mientras que la fibra de carbono se simulará como material elástico. La geometría del RVE se compone de 5 inclusiones simétricas como en la Figura 7 y la malla utilizada para el análisis está compuesta por 272 elementos lineales en pequeñas deformaciones con 305 nodos.

Figura 7. Geometría y malla RVE.

Tabla 2. Propiedades Mecánicas Resina Epoxy
Young Modulus	3.24e3	MPa
Poisson Ratio	0.32
Stress Traction Limit	80	MPa
Traction Fracture Energy	0.73	J/mm2
Stress Compression Limit	240	MPa
Compression Fracture Energy	2.19	J/mm2

Tabla 3. Propiedades Mecánicas fibras de carbono
Young Modulus	235e3	MPa
Poisson Ratio	0.21

El análisis del RVE descrito anteriormente proporciona una serie de puntos discretos con los que definir la superficie de daño o, lo que es lo mismo, una Discrete Multiscale Threshold Surface. A continuación, se presenta la sección de la superficie límite elástica para ${\textstyle {d}_{eq}=}$$0.01$ en el espacio de las tensiones normales.

Figura 8. Superficie límite elástica en el espacio de tensiones para el material compuesto.

5 Validatión

5.1 Test a cortante

Se ha realizado una prueba de cortante por una geometría cuadrada, de 1 m de lado, compuesta por elementos lineales en pequeñas deformaciones utilizando el material compuesto definido en la sección 4.2. En caso de DMTS, sólo los elementos de la macroescala que exceden el límite elástico generan el RVE mientras los que todavía siguen en el régimen lineal no necesitan la generación RVE.

En la siguiente figura se muestra el comportamiento del compuesto si este se analiza mediante un procedimiento FE2 o si se analiza con la formulación DMTS propuesta. Como se puede ver, los resultados son prácticamente idénticos.

Figura 9. Curvas Fuerza-Desplazamiento para los diferentes métodos analizados.

La Figura 9 y la Tabla 4 a continuación proporcionan respectivamente la respuesta de Fuerza-Desplazamiento de la macroescala bajo la prueba de cortante y los resultados en términos de tiempo y memoria para todos los métodos: FE2, MCM completo y adaptivo.

Tabla 4. Comparación de tiempo y memoria
Type	Time [s]	Memory [kB]	Number of Active RVEs	Time Speedup
Full Multiscale	844.72	169788	484/484	-
DMTS	646.39	92140	388/484	1.307

En las siguientes figuras podemos ver la evolución de ${\textstyle {u}_{y}}$, ${\textstyle {\sigma }_{min}}$, ${\textstyle {\epsilon }_{max}}$, ${\textstyle {d}_{eq}}$ para cada uno de los métodos analizados a diferentes niveles de deformación en la dirección y. Podemos apreciar cómo los resultados obtenidos con DMTS reflejan los valores de referencia calculados con FE2.

5.2 Refuerzo en material compuesto

En este capítulo analizamos una viga realizada con el material compuesto descrito en la sección 4.2 para validar el método DMTS propuesto.

La viga está discretizada con una malla de elementos lineales en pequeñas deformaciones (1075 elementos y 678 nodos) con una geometría como en Figura 11.

Figura 11. Geometría y malla.

La viga está sometida a un desplazamiento vertical impuesto ${\textstyle {u}_{y}}$ y fija en toda su base inferior. En este ejemplo, igual que en el anterior, se analizarán los resultados obtenidos a través de la simulación DMTS, con los resultados de un modelo multiescala de primer orden (FE2). Para reducir el coste computacional de ésta última simulación, únicamente se analizarán con RVEs los 84 elementos marcados en rojo en la Figura 11.

Como se puede observar en la Figura 10, las curvas Fuerza-Desplazamiento obtenidas mediante la formulación propuesta (DMTS) coinciden con el método FE2. Se puede apreciar también que la fuerza máxima alcanzada sea la misma para los dos casos.

La comparación de los resultados de tensiones, deformaciones y daño obtenidos con FE2 y DMTS muestra también valores muy parecidos. En las Figuras 13 y 14 se proporcionan el estado final del refuerzo y la distribución de daño en la estructura en el RVE. Mientras en la Figura 15 se observan los elementos que se han activados, generando un RVE, después de alcanzar el valor límite de ${\textstyle {d}_{eq}}$.

Figura 12. Comparación de curvas Fuerza-Desplazamiento para el refuerzo.

Figura 13. Distribución de ${\textstyle {d}_{eq}}$ en el refuerzo para el caso FE2.

Figura 14. Distribución de ${\textstyle {d}_{eq}}$ en el refuerzo para el caso DMTS.

	0.5 mm	2.25 mm	4.5 mm	${u}_{y}$
FM2
DMTS

	0.5 mm	2.25 mm	4.5 mm	${\sigma }_{min}$
FM2
DMTS

	0.5 mm	2.25 mm	4.5 mm	${\epsilon }_{max}$
FM2
DMTS

	0.5 mm	2.25 mm	4.5 mm	${d}_{eq}$
FM2
DMTS

Figura 10. Resultados de desplazamiento vertical ( ${\textstyle {u}_{y}}$), tensión ( ${\textstyle {\sigma }_{min}}$), deformación ( ${\textstyle {\epsilon }_{max}}$) y daño equivalente ( ${\textstyle {d}_{eq}}$).

Figura 15. Elementos activos, donde se requiere RVEs, para el caso DMTS.

Las figuras siguientes proponen una comparación entre FE2 y DMTS, para el mismo valor de ${\textstyle {u}_{y}=}$$2.27$ respectivamente en términos de daño equivalente, tensión máxima y mínima, deformación máxima y mínima. Podemos observar como los resultados en el caso DMTS reproducen de manera óptima el modelo FE2 de referencia.

(a)

(b)

(c)

(d)

(e)

Figura 16. Comparación de Resultados de daño equivalente ( ${\textstyle {d}_{eq}}$), tensión ( ${\textstyle {\sigma }_{max},\,{\sigma }_{min}}$) y deformación ( ${\textstyle {\epsilon }_{max},\,{\epsilon }_{min}}$) para el método FE2 (izquierda) y el método DMTS (derecha).

Finalmente, analizamos los tiempos de cálculo y la memoria necesarios para los análisis del refuerzo en los dos casos.

Gracias al método DMTS el tiempo de cálculo se ha reducido de 2 veces para el modelo adaptivo y la memoria utilizada de un 30% respecto al caso FE2, que hemos considerado simplificado analizando 84 de los 1075 elementos del refuerzo.

Table 5. Comparación de tiempo y memoria
Type	Time [s]	Memory [kB]	Number of Active RVEs	Time Speedup
Full Multiscale	40267.5	152212	84/1075	-
DMTS	19183.3	108696	70/1075	2.1

6 Conclusions

En este documento se ha analizado un nuevo método de cálculo mediante modelos multiescala para el estudio de estructuras en materiales compuestos, que optimiza el coste computacional de la simulación. El Discrete Multiscale Threshold Surface (DMTS) se basa en la definición de una base de datos que contiene el comportamiento del material, la cual se obtiene a partir del análisis de un RVE en un número discreto de estados tenso-deformacionales.

Gracias a una previa validación del método con un material simple y la caracterización de una base de datos especifica de una microestructura con fibras largas, se ha procedido al estudio de un refuerzo en fibras de carbono.

El resultado de este estudio muestra como con el DMTS se obtienen resultados comparables a los obtenidos con Métodos Multiescala de primer orden (calculado pare 84 RVEs), obteniendo ventajas relevantes en términos de coste computacional. Obviamente, un incremento de RVEs para el caso FE2 produciría un speed-up aun mayor, entendiéndose constante el tiempo de cálculo obtenidos con el método DMTS.

Los autores desean destacar que, teniendo en cuenta el elevado número de análisis necesario para la discretización del espacio de deformación, si el número total de puntos de integración de la macroescala, es decir, el número de RVEs, es menor que el número de historias de deformación necesarias para la construcción del database del material, el uso de un modelo FE2 tendrá un coste computacional más bajo. Sin embargo, esto es cierto sólo si analizamos la macroescala sólo una vez (situación muy limitada). Es suficiente que cambiemos las condiciones de borde o cualquier parámetro del análisis que necesitaremos recalcular todos los RVEs, mientras que utilizando el método DMTS no es necesario.

Agradecimientos

El presente trabajo ha sido apoyado por el European Research Council a través de la Advanced Grant: ERC-2012-AdG 320815 COMP-DES-MAT ‘‘Advanced tools for computational design of engineering materials’’, por el Ministerio de Economía y Competividad a través del proyecto: MAT2014-60647-R ‘‘Multi-scale and multi-objective optimization of composite laminate structures (OMMC)’’, por el 7º Programa Marco de la Unión Europea bajo IRSES Marie Curie Action: PIRSES-GA-2013-612607 TCAiNMaND, por la colaboración entre EU-H2020 (Acuerdo nº 690638) y la República Popular China (Acuerdo nº [2016] 92) "ECOCOMPASS", y por el Centro Internacional de Métodos Numéricos en Ingeniería (CIMNE). Todo este apoyo es reconocido con gratitud.

Referencias

[1] S. Oller et al., An anisotropic elastoplastic constitutive model for large strain analysis of fiber reinforced composite materials, Computer Methods in Applied Mechanics and Engineering, 185, pág. 245-277 (2000).

[2] F. Rastellini et al., Composite materials non-linear modelling for long fibre-reinforced laminartes: Continuum basis, computational aspects and validation, Computers & Structures, 86, pág. 879-896 (2008).

[3] X. Martínez et al., A numerical procedure simulating RC structures reinforced with FRP using serial/parallel mixing theory, Computers & Structures, 86, pág. 1604-1618 (2008).

[4] S. Oller, Numerical simulation of mechanical behavior of composite materials, 1st Edition, CIMNE-Springer (2014).

[5] F. Otero et al., Study and prediction of the mechanical performance of nanotube-reinforced composite, Composite Structures, 94, pág. 2920-29300 (2012).

[6] E. Sanchez-Palencia, Non-homogeneous media and vibration theory, Lecture Notes in Physics, 127, (1980).

[7] P. Suquet, Local and global aspects in the mathematical theory of plasticity, A. Sawczuk, G. Bianchi (Eds.), (1985).

[8] F. Montero et al., Multiscale thermo-mechanical analysis of multi-layered coatings in solar thermal applications, Finite Elements in Analysis and Design, 27, pág. 31-43 (2017).

[9] F. Otero, X. Martinez, S. Oller and O. Salomón, Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations, Composite Structures, 122, pág. 405-416 (2015)

[10] Q.-C. H. J. Yvonnet, The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains, J. Comput. Phys., 223, pág. 341-368 (2007).

[11] J.A. Hernández et al., High-performance model reduction techniques in computational multiscale homogenization, Comput. Methods Appl. Mech. Eng., 176, pág. 149-189 (2014).

[12] P. W. I. Temizer, An adaptive method for homogenization in orthotropic nonlinear elasticity, Comput. Methods Appl. Mech. Eng., 35–36, pág. 3409-3423 (2007).

[13] F. Otero et al., An efficient multi-scale method for non-linear analysis of composite structures, Comp. Struct., 131, pág. 707-719 (2015).

[14] F. Otero et al., Multiscale computational homogenization: review and proposal of a new enhanced-first-order method, Archives of Computational Methods in Engineering, pág. 1-27 (2016).

[15] T. Kanit et al., Determination of the size of the representative volume element for random composites: statistical and numerical approach., Int. J. Solids Struct., 13-14, pág. 40 (2003).

[16] M. O.-S. a. I. J. M. Jiang, Scale-dependent bounds on effective elastoplastic response of random composites, J. Mech. Phys. Sol., 49, pág. 655-673 (2001).

[17] M. O.-S. a. I. J. M. Jiang, Apparent elastic and elastoplastic behavior of periodic composites, Int. J. Solids Struct., 39, pág. 199-212 (2002).

[18] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Prob. Eng. Mech., 21, pág. 112-132 (2006).

[19] S. Ghosh, K. Lee and S. Moorthy, Two scale analysis of heterogeneous elastic-plastic materials with asymptotic homogenization and Voronoi cell finite element model, Computer Methods in Applied Mechanics and Engineering, 132, pág. 63–116 (1996).

[20] F. Otero, X. Martinez, S. Oller and O. Salomón, An efficient multi-scale method for non-linear analysis of composite structures, Composite Structures, 131, pág. 707-719 (2015).

[21] M. P. L. R. R. e. a. Petracca, Regularization of first order computational homogenization for multiscale analysis of masonry structures, Comput. Mech., 57, pág. 257 (2016).

Referencias

[1] S. Oller et al., An anisotropic elastoplastic constitutive model for large strain analysis of fiber reinforced composite materials, Computer Methods in Applied Mechanics and Engineering, 185, pág. 245-277 (2000).

[2] F. Rastellini et al., Composite materials non-linear modelling for long fibre-reinforced laminartes: Continuum basis, computational aspects and validation, Computers & Structures, 86, pág. 879-896 (2008).

[3] X. Martínez et al., A numerical procedure simulating RC structures reinforced with FRP using serial/parallel mixing theory, Computers & Structures, 86, pág. 1604-1618 (2008).

[4] S. Oller, Numerical simulation of mechanical behavior of composite materials, 1st Edition, CIMNE-Springer (2014).

[5] F. Otero et al., Study and prediction of the mechanical performance of nanotube-reinforced composite, Composite Structures, 94, pág. 2920-29300 (2012).

[6] E. Sanchez-Palencia, Non-homogeneous media and vibration theory, Lecture Notes in Physics, 127, (1980).

[7] P. Suquet, Local and global aspects in the mathematical theory of plasticity, A. Sawczuk, G. Bianchi (Eds.), (1985).

[8] F. Montero et al., Multiscale thermo-mechanical analysis of multi-layered coatings in solar thermal applications, Finite Elements in Analysis and Design, 27, pág. 31-43 (2017).

[9] F. Otero, X. Martinez, S. Oller and O. Salomón, Numerical homogenization for composite materials analysis. Comparison with other micro mechanical formulations, Composite Structures, 122, pág. 405-416 (2015)

[10] Q.-C. H. J. Yvonnet, The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains, J. Comput. Phys., 223, pág. 341-368 (2007).

[11] J.A. Hernández et al., High-performance model reduction techniques in computational multiscale homogenization, Comput. Methods Appl. Mech. Eng., 176, pág. 149-189 (2014).

[12] P. W. I. Temizer, An adaptive method for homogenization in orthotropic nonlinear elasticity, Comput. Methods Appl. Mech. Eng., 35–36, pág. 3409-3423 (2007).

[13] F. Otero et al., An efficient multi-scale method for non-linear analysis of composite structures, Comp. Struct., 131, pág. 707-719 (2015).

[14] F. Otero et al., Multiscale computational homogenization: review and proposal of a new enhanced-first-order method, Archives of Computational Methods in Engineering, pág. 1-27 (2016).

[15] T. Kanit et al., Determination of the size of the representative volume element for random composites: statistical and numerical approach., Int. J. Solids Struct., 13-14, pág. 40 (2003).

[16] M. O.-S. a. I. J. M. Jiang, Scale-dependent bounds on effective elastoplastic response of random composites, J. Mech. Phys. Sol., 49, pág. 655-673 (2001).

[17] M. O.-S. a. I. J. M. Jiang, Apparent elastic and elastoplastic behavior of periodic composites, Int. J. Solids Struct., 39, pág. 199-212 (2002).

[18] M. Ostoja-Starzewski, Material spatial randomness: From statistical to representative volume element, Prob. Eng. Mech., 21, pág. 112-132 (2006).

[19] S. Ghosh, K. Lee and S. Moorthy, Two scale analysis of heterogeneous elastic-plastic materials with asymptotic homogenization and Voronoi cell finite element model, Computer Methods in Applied Mechanics and Engineering, 132, pág. 63–116 (1996).

[20] F. Otero, X. Martinez, S. Oller and O. Salomón, An efficient multi-scale method for non-linear analysis of composite structures, Composite Structures, 131, pág. 707-719 (2015).

[21] M. P. L. R. R. e. a. Petracca, Regularization of first order computational homogenization for multiscale analysis of masonry structures, Comput. Mech., 57, pág. 257 (2016).