Simulation of damage phenomena in metamaterials

<!---<sup>1</sup>Escuela Técnica Superior de Ingeniería Aeronáutica y del Espacio, Universidad Politécnica de Madrid, <sup>2,3</sup>Grupo Avanzado De Modelado Y Simulación No-Lineal De Sólidos, Escuela Técnica Superior de Ingeniería Aeronáutica y del Espacio, Universidad Politécnica de Madrid.

Pza. del Cardenal Cisneros, 3, 28040 Madrid-España
<!--<sup>*</sup>rodrigo.montero.mayoral@alumnos.upm.es-->
==Resumen==

Los metamateriales se generan a partir de una red interrelacionada de celdas (caracterizadas por su geometría y conectividad) y pueden presentar un comportamiento a nivel macroscópico que difiere del que caracteriza a sus constituyentes más elementales. La modelización de la influencia que el daño y fractura de los elementos constituyentes tienen en las propiedades macroscópicas del metamaterial es de relevancia para su estudio mecánico. Alterando el comportamiento resistente a nivel elemento (variando su módulo elástico longitudinal), es posible aproximar el efecto del daño en la estructura completa. Es decir, se considera válida la sustitución del proceso de dañado y deformación plástica por una serie de estados intermedios contenidos en el régimen elástico lineal. En estos estados, la geometría del elemento permanece inalterada, y cada uno de ellos queda caracterizado por un valor del módulo elástico. En un proceso iterativo, los elementos que soportan tensiones superiores a las impuestas como fronteras, progresan a través de los diferentes estados, disminuyendo su módulo elástico, hasta que se consideran rotos y son eliminados de la estructura, pudiéndose alcanzar el colapso de la propia estructura. Al calcular, se desprecian los efectos másicos y se parte de un estado libre de tensiones residuales en cada iteración. Se presenta un algoritmo sencillo para simular el efecto del daño en estructuras de metamaterial, aplicable a cualquier programa de elementos finitos.

'''Palabras clave''': Metamateriales, daño, fractura 

==Abstract==

Metamaterials are generated from an interrelated set of cells and can present a macroscopic behaviour that differs from the one that characterizes its basic constituents. Modelling the influence that the damage and fracture of the elemental constituents have in the macroscopic properties of the metamaterial is relevant for its mechanical analysis. Altering the resistant behaviour of the elements (changing their longitudinal elastic modulus), it is possible to approximate the effect of the damage in the complete structure.  It is considered as valid the substitution of the damage and plastic deformation phenomena by intermediate states contained in the linear elastic deformation regime. Each state is characterized by the elastic module of the element, whose geometry remains unmodified. In an iterative process, when the induced stress in the elements is greater than the one stablished as the limit, they progress through the different states, diminishing their elastic modulus until they are considered as fractured and are eliminated from the structure. Mass effects are ignored, and a stress-free structure is implemented to calculate each iteration. A simple algorithm is presented to simulate the effect of damage in metamaterial structures, applicable to any finite element software.

'''Keywords''': Metamaterials, damage, fracture 

==1. Introducción==

Los metamateriales son materiales formados a partir de una red interrelacionada de celdas, caracterizadas por su geometría y conectividad <span id='citeF-1'></span>[[#cite-1|[1]]] <span id='citeF-2'></span>[[#cite-2|[2]]]. Estas se ensamblan (tanto en el plano como en el espacio) para obtener las caras y/o vértices del sólido a generar. Las propiedades que presentan a nivel macroscópico pueden diferir de aquellas que caracterizan a sus constituyentes más elementales. El propio nombre sugiere este comportamiento, pues el prefijo ''meta-'', derivado de la preposición <math display="inline">\mu \varepsilon \tau \acute{\alpha }</math>, equivale a ''más allá'' <span id='citeF-3'></span>[[#cite-3|[3]]], lo que indica que a partir de determinadas configuraciones de celda, es posible obtener propiedades anómalas <span id='citeF-4'></span>[[#cite-4|[4]]]-<span id='citeF-5'></span>[[#cite-5|[5]]], como el comportamiento mecánico auxético <span id='citeF-6'></span>[[#cite-6|[6]]] o índices de refracción negativos <span id='citeF-7'></span>[[#cite-7|[7]]]. Este tipo de materiales no son extraños a la naturaleza, que presenta una amplia variedad, donde la distribución de celdas dentro del sólido y la propia estructura de cada celda varían desde disposiciones altamente regulares (''honeycomb'' <span id='citeF-8'></span>[[#cite-8|[8]]]) hasta completamente aleatorias (esponjas y espumas). Esto lleva a diferenciar entre metamateriales homogéneos (caracterizados por una distribución periódica de celdas) y no homogéneos (caracterizados por una distribución de celdas no periódica) <span id='citeF-9'></span>[[#cite-9|[9]]]. En ciertos campos, como el electromagnetismo, es posible emplear provechosamente las propiedades de los metamateriales no homogéneos (como, por ejemplo, las espumas <span id='citeF-10'></span>[[#cite-10|[10]]]), mientras que en otros, como el del diseño estructural, los homogéneos presentan cualidades que resultan más apropiadas <span id='citeF-11'></span>[[#cite-11|[11]]] <span id='citeF-12'></span>[[#cite-12|[12]]]. Históricamente, la posibilidad de fabricar estos materiales ha ido ligada al desarrollo de técnicas de impresión 3D <span id='citeF-13'></span>[[#cite-13|[13]]], que permite obtenerlos a partir de la idealización de la celda unidad y del uso de leyes topológicas que gobiernan las posiciones que toma. Esta tecnología está en continuo desarrollo, por lo que los métodos utilizados no están estandarizados <span id='citeF-14'></span>[[#cite-14|[14]]]. 

El creciente interés por los metamateriales artificiales se debe no solo a la dependencia de las propiedades macroscópicas de la celda de partida <span id='citeF-15'></span>[[#cite-15|[15]]]-<span id='citeF-16'></span>[[#cite-16|[16]]], sino también a la capacidad de disminuir la densidad que presentan las estructuras de metamaterial frente a sólidos homogéneos fabricados a partir del mismo material base <span id='citeF-17'></span>[[#cite-17|[17]]]. En el ámbito de la resistencia de materiales se utilizan de forma habitual metamateriales homogéneos, lo que facilita la modelización de la respuesta a estados de cargas con capacidad de dañar y romper un número significativo de elementos. 

Para conocer la respuesta mecánica y afrontar el problema de diseño, se precisa de un modelo que incluya los efectos derivados del fallo de los elementos constituyentes del metamaterial en su respuesta mecánica, ya que por construcción y debido a su dimensión respecto al metamaterial en conjunto, son más susceptibles de fallar. La evolución y acumulación de estos fallos lleva al colapso de la estructura. Por establecer una analogía con los materiales homogéneos, este proceso de evolución del efecto de daño desde el nivel fundamental hasta el macroscópico es similar al proceso de nucleación, crecimiento y propagación de grietas en roturas frágiles o al de formación, crecimiento y coalescencia de microvacíos en roturas dúctiles.

==2. Método==

===2.1 Simplificación del efecto del daño===

Tomando inspiración de diversos metamateriales presentes en la naturaleza, donde las aristas de cada celda son esbeltas, se toma como elemento de partida para generar las estructuras el elemento 1D barra-viga de dos nodos (''beam''), de sección cuadrada y constante.

Para simular el efecto del daño que genera una carga en la estructura de metamaterial, se acepta como válida la sustitución del proceso de dañado y deformación plástica a nivel elemento por una serie de estados intermedios en los que se implementan las hipótesis del régimen elástico lineal. Cada estado intermedio queda caracterizado por la geometría y el módulo de elasticidad del elemento <math display="inline">E</math>, al que se le asignan valores decrecientes, desde el valor característico del material base, para el primer estado (cuando el elemento se encuentra dentro del régimen de deformación elástico), hasta que se anula (cuando se considera que el elemento ha fallado). La geometría se mantiene constante e igual a la del estado nominal, por lo que se está ignorando el proceso real de deformación propio de un sólido homogéneo metálico con cierta capacidad dúctil.

La idea del efecto de daño que se ha incluido es considerar únicamente el mayor valor de la tensión principal en el elemento de metamaterial, una vez que dicha tensión principal supera un valor umbral, se supone que parte del elemento ha sido dañado, por tanto, se multiplica el valor del módulo de elasticidad del material en el elemento por un factor que intenta representar el comportamiento introducido por esa sección dañada, añadiendo flexibilidad a su comportamiento,  dando lugar al siguiente paso de la simulación. A partir de un valor, se supondrá que el elemento está totalmente dañado y se elimina de la simulación.

No se dispone de referencias en las que apoyarse para determinar con precisión los valores que debería tomar <math display="inline">E</math> en los diferentes estados, por ello se utiliza una distribución sencilla. 

===2.2 Evolución de los elementos===

Una vez descrito el proceso basado en estados intermedios que sustituye a la deformación plástica y al crecimiento de grietas o microrroturas, es necesario establecer el número de estados con el que trabajar y el criterio para pasar de un estado al siguiente. 

En el estudio se ha considerado que tras diez estados de modificación del módulo de Young, pasos que representan una flexibilización gradual del elemento debido al daño sufrido, este falla. Es claro que es una simplificación y que no se incluye efectos de grandes desplazamientos, plasticidad o fallo por pandeo elástico, modos que son muy importantes en la realidad, pero que no han sido incluidos en nuestro análisis para simplificar el cálculo. Esta simplificación del efecto del daño es excesiva para materiales dúctiles y en casos de pandeo,  pero puede emplearse como primera aproximación en metamateriales. Se elige 10 como el número de estados intermedios y un rango de valores de <math display="inline">E</math> decrecientes, que comprende desde el módulo de elasticidad característico del material base hasta el 80% de ese valor, aparte del valor nulo. Este número de estados es por cuestiones de limitar los pasos de cálculo y el valor límite es por cercanía al valor de fallo del material, aluminio.

Para establecer el criterio de paso de un estado al siguiente, se analizan las tensiones que soportan los elementos. Se supone el caso más sencillo, en el que las tensiones axiles son las más relevantes, siendo las cortantes despreciables frente a estas. Al trabajar con estructuras en el plano, el tensor de tensiones queda:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }  =  \begin{pmatrix}\sigma _x     & 0     & 0 \\   \\ 0     & 0     & 0 \\   \\ 0     & 0     & 0 \\\end{pmatrix}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

La tensión principal (<math display="inline">\sigma _1</math>), puede aproximarse como la suma de la tensión normal debida al momento flector y la tensión axil. A partir de los criterios de Rankine y de Tresca, y siendo únicamente la tensión principal <math display="inline">\sigma _1</math> distinta de cero, se establece que el máximo de esta tensión sobre el elemento es el valor con el que se debe comparar para determinar si el elemento a estudiar, perteneciente a un estado intermedio (asociado a un <math display="inline">E_{Nk}</math>), debería pasar al siguiente estado (asociado a <math display="inline">E_{Nk+1}</math>), acercándolo así al fallo. Para el paso desde el estado inicial al primer estado de daño, la tensión con la que comparar <math display="inline">\sigma _1</math> es la tensión de fluencia del material base. Para las tensiones que establecen la frontera entre el resto de estados intermedios, se dispone un rango que comprende desde la tensión de fluencia, hasta un 130% de ese valor. A partir de estos razonamientos se impone (en MPa):

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> E_N  = 10^3 \cdot  \begin{pmatrix} 70  &  68  &  67  &  66  &  65  &  64 &  63 &  62  & 61 & 60  & 57 \end{pmatrix} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \sigma _N  = \begin{pmatrix} 90  &  95  &  96  &  97  &  98  &  99  &  100  &  105  &  110  &  114  &  115 \end{pmatrix}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}


La precisión de los valores no es crítica para el proceso, ya que el objetivo es determinar si el método permite representar en primera aproximación el comportamiento de las estructuras de estudio. El objetivo es simular el comportamiento general del proceso de dañado en metamateriales, desde un punto de vista cualitativo, con tal de disponer de una herramienta que permita su simulación y ajustar los valores para casos concretos.

En la estructura de partida, todos los elementos están caracterizados por el mismo <math display="inline">E_{N1}</math> y la misma tensión <math display="inline">\sigma _{N1}</math>.  

Trabajar con más de un estado intermedio hace mandatorio el uso de un proceso iterativo, para que los elementos tengan la posibilidad de atravesar todos ellos. En cuanto a la transición entre diferentes estados, se aplican ciertas restricciones. Solo se puede pasar de un estado de daño, al inmediatamente posterior, independientemente de la tensión que esté soportando el elemento. No se puede ir hacia estados anteriores, ni se puede evitar pasar por alguno de los estados. Por lo tanto, para que un elemento falle, tiene que haber pasado antes por todos los estados.

===2.3 Generación de las estructuras===

La generación de las geometrías se lleva a cabo a partir del uso de simetrías y traslaciones de las posiciones de los nodos de la celda base. El material de partida utilizado es isótropo, con un módulo de elasticidad <math display="inline">E= 70  </math> GPa y un coeficiente de Poisson <math display="inline">\nu = 0.3 </math>, representativo del <math display="inline">Al1080A</math>. La sección de los elementos es cuadrada, de lado <math display="inline">l= 0.5 </math> mm constante. 

Intentando imitar los ensayos experimentales, el parte superior e inferior del dominio modelado, se introduce una superficie que actúa como interfaz donde se aplican los desplazamientos, equivalente a las mordazas o material donde se pega la probeta para ensayos. Por tanto, el estado de carga no se impone sobre los nodos del núcleo, sino sobre las interfaces (a través de la aplicación de desplazamientos nodales que provocan tracción). Como método de unión entre el dominio de estudio y las interfaces se utilizan <nowiki>Multi Point Constrains (MPCs)</nowiki>, estableciendo una unión entre mallas no conformes.

==3. Algoritmo==

Se comienza definiendo la geometría de la celda de partida, que se utiliza para generar la estructura. Una vez construida la estructura, se impone el estado de carga. Para cada elemento se obtiene el máximo de <math display="inline">\sigma _1</math> y el valor de <math display="inline">E</math> (perteneciente a <math display="inline">E_N</math>) que lo particulariza. Seguidamente, se compara <math display="inline">\sigma _1</math> máximo con el valor que se ha impuesto en el vector <math display="inline">\sigma _N</math> como frontera del estado intermedio que caracteriza al elemento. En caso de ser <math display="inline">\sigma _1</math> máximo mayor o igual al valor impuesto en <math display="inline">\sigma _N</math>, el elemento en cuestión progresa al siguiente estado (deja de estar identificado por <math display="inline">E_{Nk}</math>, y pasa a estarlo por <math display="inline">E_{Nk+1}</math>), o es eliminado de la estructura (en el caso de que se encontrara en el último estado). En el caso opuesto, se mantiene el estado en el que se encontraba el elemento (sigue estando asociado a <math display="inline">E_{Nk}</math>). Una vez analizados todos los elementos, se actualiza la estructura con las modificaciones pertinentes. 

En todas las iteraciones se parte de la geometría original, pero se modifican los valores de <math display="inline">E</math>  en aquellos elementos que se detectan como dañados y se eliminan de la estructura los que se definen como fallados. Aplicándose a esa estructura resultante el desplazamiento para volver a iterar.

Al modificar las propiedades de unos elementos y eliminar otros de la estructura, la redistribución de cargas provoca la activación de nuevos elementos, lo permite simular la evolución del efecto del daño. Este proceso se esquematiza mediante el diagrama mostrado en la [[#img-1|Figura 1]].

<div id='img-1'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-diagrama_7.png|320px|Diagrama del algoritmo.]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 1'''. Diagrama del algoritmo
|}

===3.1 Criterios de parada===

En caso de que la redistribución de cargas no genere modificaciones de <math display="inline">E</math> en ningún elemento, se considera que la estructura ha alcanzado el equilibrio, ya que las tensiones soportadas son insuficientes para hacer evolucionar a estructura hacia un nuevo estado de daño. Como comportamiento contrapuesto aparece el fallo global de la estructura, lo que se deduce de la aparición de mecanismos (pues ya no se tiene una sola estructura, sino varias) que impiden que se pueda seguir iterando. Dada esta situación, se considera que la estructura ha colapsado, al ser el estado de carga impuesto excesivo. Estos dos criterios se obtienen como respuesta de la estructura al caso de carga. Cabe mencionar la posible aparición de elementos subapoyados formando mecanismos en el interior de la estructura, debido a que los elementos adyacentes han fallado. En esta situación, la estructura no ha colapsado, por ello, se debe continuar iterando. Para distinguir entre el colapso total y la aparición de mecanismos en el interior de la estructura, se eliminan los elementos asociados a los nodos que forman los mecanismos (pues no están trabajando) y se comprueba si los elementos restantes forman una estructura  con todos los grados de libertad restringidos. Este proceso es a su vez iterativo, pues al quitar elementos, puede que otros que se encontraban apoyados, dejen de estarlo. Como último criterio, se limita el número máximo de iteraciones.

==4. Resultados y discusión==

===4.1 Análisis con estructuras jerárquicas===

Se plantea una estructura jerárquica de triángulos equiláteros, basada en la geometría mostrada en la [[#img-2|Figura 2]] y empleando dos niveles de jerarquización. Tanto el primer como segundo nivel constan de setenta y cuatro celdas de nivel anterior. La geometría y el caso de carga están basadas en <span id='citeF-18'></span>[[#cite-18|[18]]]. Para cada lado de cada triángulo equilátero, se utiliza un solo elemento, de longitud <math display="inline">1 </math> mm. 

<div id='img-2'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-t_jerarquia_1.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 2'''. Primer nivel de jerarquización
|}

<br>
<br>
Como estado de carga, se restringe el desplazamiento vertical en la interfaz inferior, se impone un desplazamiento según esa misma dirección en la interfaz superior y se empotra el nodo de la esquina inferior izquierda.

El modelo consta de <math display="inline">4.20224 \cdot 10^{5}</math> elementos en el núcleo, tal y como se muestra en la  [[#img-3|Figura 3]], y para un desplazamiento impuesto de <math display="inline">31 \cdot 10^{-2}</math> mm requiere <math display="inline">361</math> iteraciones hasta alcanzar el equilibrio.

<div id='img-3'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-t_4_4_iter_0.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 3'''. Geometría de la estructura inicial
|}


En la [[#img-4|Figura 4]] se aprecia el detalle de la zona de la que se eliminan elementos de la estructura de manera deliberada, entorno al plano de simetría horizontal.

<div id='img-4'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-t_4_4_iter_0_detalle_4.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 4'''. Detalle de la geomatría de la estructura
|}

<br>

En la [[#img-5|Figura 5]] se observa que las tensiones máximas se obtienen en el extremo derecho de la zona de la que se han eliminado elementos deliberadamente, coincidiendo con el resultado obtenido en <span id='citeF-18'></span>[[#cite-18|[18]]]. 
<br>
<div id='img-5'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[File:Montero-Mayoral_et_al_2023a_9544_img5_montero.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 5'''. Detalle de la distribución de tensiones en la primera iteración
|}


Hasta alcanzar el equilibrio ([[#img-6|Figura 6]]), se produce el fallo de numerosos elementos, principalmente en la parte superior derecha, y en menor medida en la inferior derecha. 

<div id='img-6'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-t_4_4_iter_386_stress.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 6'''. Distribución de tensiones una vez alcanzado el equilibrio
|}


En cuanto a los valores de <math display="inline">E</math> ([[#img-7|Figura 7]]), solo se han visto modificados en las inmediaciones de las posiciones en las que se ha producido el fallo de elementos. En el resto de la estructura, el valor de <math display="inline">E</math> se mantiene inalterado.

<div id='img-7'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-t_4_4_iter_386_E.png|282px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 7'''. Distribución de <math>E</math> una vez alcanzado el equilibrio
|}


En la región en la que aparece la concentración de tensiones se da la rotura del primer elemento (lo que se podría equiparar con la nucleación de una grieta). Posteriormente, los elementos van fallando a través de dos frentes, hacia el centro de la estructura, uno por encima y otro por debajo del plano de simetría horizontal. Este avance es semejante al que presentaría un sólido homogéneo, debido al crecimiento de grietas.

===4.2 Análisis con estructuras cuadrangulares===

La estructura se genera a partir de un cuadrado rotado 45º. Como caso de carga, en la interfaz inferior, los nodos quedan empotrados, mientras que en la superior se restringe el desplazamiento vertical, a la vez que se impone un desplazamiento horizontal, simulando así un caso de carga a cortadura. Para esta misma estructura, en <span id='citeF-19'></span>[[#cite-19|[19]]] se calcula el fallo del primer elemento mediante procedimientos numéricos y experimentales, dándose este en ambos casos en la esquina inferior izquierda. El modelo consta de 1.323 elementos en el núcleo, y para un desplazamiento impuesto de <math display="inline">11 \cdot 10^{-2}</math> mm requiere <math display="inline">72</math> iteraciones hasta alcanzar el equilibrio. 

El primer fallo aparece en la esquina inferior izquierda ([[#img-8|Figura 8]]), coincidiendo con el resultado obtenido en <span id='citeF-19'></span>[[#cite-19|[19]]].

<div id='img-8'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-r_20_E.png|186px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 8'''. Aparición del fallo de los primeros elementos
|}


<div id='img-9'></div>
{| class="wikitable" style="margin: 1em auto 0.1em auto;border-collapse: collapse;width:auto;" 
|-style="background:white;"
|style="text-align: center;padding:10px;"| [[Image:Draft_Montero_348638231-r_72_E.png|186px|]]
|-
| style="background:#efefef;text-align:left;padding:10px;font-size: 85%;"| '''Figura 9'''. Distribución de <math>E</math> una vez alcanzado el equilibrio
|}


Para la carga propuesta, hasta alcanzar el equilibrio, los elementos van fallando a través de un único frente, hacia el interior. Por otra parte, el efecto del daño no evoluciona significativamente en el resto de esquinas respecto de la iteración en la que se produjo el fallo de los primeros elementos ([[#img-9|Figura 9]]).

==5. Conclusiones==

El método propuesto, aunque está basado en hipótesis que ignoran parte de la física del problema (omite el proceso de deformación plástica y no tiene en cuenta las tensiones residuales) genera resultados que, en primera aproximación, permiten estimar a partir de un sencillo algoritmo la aparición y propagación del efecto de daño en estructuras de metamaterial, a través del modelado de las condiciones necesarias para que se produzca la fractura a nivel elemento. Estableciendo una serie de alternativas, relativas a la situación en la que se puede encontrar cada elemento, en función de la tensión a la que está sometido y las condiciones de daño que presenta, se recogen los diferentes escenarios que tienen lugar de manera recurrente, por lo que implementando un proceso iterativo es posible obtener la evolución de la estructura. Como parámetros de control, se pueden elegir el número de estados intermedios hasta la rotura y el valor del módulo elástico y la tensión máxima que caracterizan cada estado. La selección de dichos valores permite acomodar el algoritmo a requerimientos específicos, lo cual le hace bastante versátil en cuanto a la posibilidad de representar diferentes materiales y tipos de rotura.

Los resultados obtenidos son consistentes con lo calculado mediante otros procediminetos, ya sea para la distribución de tensiones a obtener (caso de las estructuras jerárquicas triangulares, al comparar con <span id='citeF-'></span>[[#cite-18|[18]]]), o para la zona en la que se da el fallo del primer elemento (caso de las estructuras cuadrangulares, al comparar con <span id='citeF-19'></span>[[#cite-19|[19]]]). En primera aproximación, se considera que el modelo genera resultados representativos del comportamiento observado en la bibliografía para este tipo de estructuras. 

En conclusión, se presenta un algoritmo que permite simular el efecto de daño de forma sencilla en estructuras de metamaterial, de forma rigurosa se trata de estructura celulares jerárquicas, pero el objetivo es aplicarlo al estudio de metamateriales, que se puede implementar en cualquier programa de elementos finitos (sin necesidad de modificar su código) y cuyos resultados coinciden con los encontrados en la bibliografía.

Si bien en mecánica del daño continuo se puede aproximar el daño mediante una transformación del área efectiva o de la rigidez efectiva,  se debe tener en cuenta la clara definición del criterio de inicio y evolución del daño, puntos que superan el objetivo de este estudio y cosntituyen un objetivo posterior. La idea es incluir una rigidez variable o quizás área efectiva variable como un resultado de la microestructura de sólidos celulares impresos en 3D, pues las porosidades son comunes en estos procesos de fabricación y permitir así una modelización sencilla del proceso de rotura de estos materiales, que permita implementar un modelo de daño mejor definido basado en datos experimentales y estudios de las propiedades de los materiales.

==Agradecimientos==

R.V. Montero-Mayoral, M.A. Sanz, V. F. González Albuixech y V. Yanes reconocen la financiación de European Union's Horizon 2020 Marie Sklodowska-Curie Actions-Innovative European Training networks del acuerdo de subvención  Nº 956401
<!--R.V. Montero-Mayoral, M.A. Sanz, V. F. Gonz\'alez Albuixech and V. Yanes aknoweledge the funding from the European Union's Horizon 2020 marie Sklodowska-Curie Actions-Innovative European Training networks under grant agreement Nº 956401-->

==Referencias==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">
<div id="cite-1"></div>
[[#citeF-1|[1]]] Wartak M.S., Tsakmakidis K.L.,  Hess O. Introduction to metamaterials. Physics in Canada, 67(1):30&#8211;34, 2011.

<div id="cite-2"></div>
[[#citeF-2|[2]]] Torres-Silva H. Una revisión sobre metamateriales. Publicaciones DYNA SL. DYNA New Technologies, 4:8&#8211;9, 2017.

<div id="cite-3"></div>
[[#citeF-3|[3]]] Kshetrimayum R.S. A brief intro to metamaterials. IEEE Potentials, 23(5):44&#8211;46, 2004.

<div id="cite-4"></div>
[[#citeF-4|[4]]] Lapine M.,  Tretyakov S. Contemporary notes on metamaterials.  IET Microwaves, Antennas & Propagation, 1(1):3&#8211;11, 2007.

<div id="cite-5"></div>
[[#citeF-5|[5]]] Bilotti F., Sevgi L. Metamaterials: Definitions, properties, applications, and FDTD-based modeling and simulation. International Journal of RF and Microwave Computer-Aided Engineering, 22(4):422&#8211;438, 2012.

<div id="cite-6"></div>
[[#citeF-6|[6]]] Castellanos L.M., Lopez, F.,  Reyes-Vera, E. Metamateriales: principales características y aplicaciones. Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, 40(156):395&#8211;401, 2016.

<div id="cite-7"></div>
[[#citeF-7|[7]]] Ren X.,  Das R.,  Tran P.,  Ngo T.D.,  Xie Y.M. Auxetic metamaterials and structures: a review. Smart Materials and Structures, 27(2):023001, 2018.

<div id="cite-8"></div>
[[#citeF-8|[8]]] Smith D.R,  Pendry J.B., Wiltshire M.C.K. Metamaterials and negative refractive index.  American Association for the Advancement of Science, 305(5685):788&#8211;792, 2004.

<div id="cite-9"></div>
[[#citeF-9|[9]]] Tancogne-Dejean T., Karathanasopoulos N.,  Mohr D.M. Stiffness and strength of hexachiral honeycomb-like metamaterials. Journal of Applied Mechanics, 86(11):111010, 2019.

<div id="cite-10"></div>
[[#citeF-10|[10]]] Cui T.J.,  Smith D.R., Liu R. Metamaterials: Theory, Design, and Applications. Springer, pp. 367, 2010.

<div id="cite-11"></div>
[[#citeF-11|[11]]] Longhui H.,  Lianwen D., Yuhan L., Heng L., Jun H., Shengxiang H., Hong C. Wide-angle microwave absorption performance of polyurethane foams combined with cross-shaped metamaterial absorber. Results in Physics, 11:769&#8211;776, 2018.

<div id="cite-12"></div>
[[#citeF-12|[12]]] Kadic M., Bückmann T.,  Schittny R., Wegener M. Metamaterials beyond electromagnetism.  Reports on Progress in Physics, 76(12):126501, 2013.

<div id="cite-13"></div>
[[#citeF-13|[13]]] Sanchez Torres L.X., Herrera C., Restrepo M. Estudio de la aplicación de los metamateriales en el diseño de antenas textiles. Revista Politécnica, 15:108&#8211;127, 2019.

<div id="cite-14"></div>
[[#citeF-14|[14]]] Yoon G.,  Kim I.,  Rho J. Challenges in fabrication towards realization of practical metamaterials. Microelectronic Engineering, 163:7&#8211;20, 2016.

<div id="cite-15"></div>
[[#citeF-15|[15]]] Gwanho Y., Inki K.,  Junsuk R. Challenges in fabrication towards realization of practical metamaterials. Microelectronic Engineering, 163:7&#8211;20, 2016.

<div id="cite-16"></div>
[[#citeF-16|[16]]] Kelly A. Why engineer porous materials? Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 364(1838):5&#8211;14, 2006.

<div id="cite-17"></div>
[[#citeF-17|[17]]] Liu Y.,  Zhang X. Metamaterials: a new frontier of science and technology.  Chemical Society Reviews, 40(5):2494&#8211;2507, 2011.

<div id="cite-18"></div>
[[#citeF-18|[18]]] Walser R.M. Electromagnetic metamaterials.  Proceedings: SPIE 4467 Complex Mediums II: Beyond Linear Isotropic Dielectrics, San Diego, 4467:1&#8211;15, 2001.

<div id="cite-19"></div>
[[#citeF-19|[19]]] Gibson L.J.,  Ashby M.F. Cellular solids: Structure and properties. Cambridge University Press, pp. 532, 1999.

<div id="cite-20"></div>
[[#citeF-20|[20]]] Shekhawat A. Toughness and damage tolerance of fractal hierarchical metamaterials. arXiv preprint arXiv:1611.01719, 2016.

<div id="cite-21"></div>
[[#citeF-21|[21]]] Turco E.,  Dell’Isola F.,  Rizzi N.L.,  Grygoruk R.,  Müller W.H, Liebold C. Fiber rupture in sheared planar pantographic sheets: numerical and experimental evidence. Mechanics Research Communications,  76:86&#8211;90, 2016.