Simulación estructural de espumas de aluminio a partir de imágenes 2D mediante la combinación de técnicas de homogeneización y machine learning

Structural simulation of 2D aluminium foam images by the use of homogenization and machine learning techniques

Borja Ferrandiz^1a, Manuel Tur^1b, Enrique Nadal^1c

Centro de Investigación en Ingeniería Mecánica (CIIM), Departamento de Ingeniería Mecánica y de Materiales ( DIMM), Universitat Politècnica de València, Camino de Vera, s/n 46022 Valencia, España. Email: ^aborferca@upvnet.upv.es, ^bmanuel.tur@mcm.upv.es, ^cennaso@upvnet.upv.es

RESUMEN

En la industria actual, el uso de materiales resistentes, rígidos, de bajo peso y con buenas propiedades tanto acústicas como térmicas es de gran interés. Entre estos materiales encontramos las espumas de aluminio. Para su uso, es necesario conocer su comportamiento estructural. Para la obtención de la geometría de una espuma de aluminio se pueden plantear diversas técnicas, todas ellas basadas en que la información inicial proviene de una imagen obtenida mediante una Tomografía Axial Computarizada (TAC). Una posible metodología, conocida comúnmente como segmentación, consiste en generar un CAD a partir de la imagen y de ahí el modelo de Elementos Finitos (EF). Otra opción es usar técnicas como el CellFEM o el cgFEM, donde cierta cantidad de píxeles, que definen las propiedades del material, son embebidos en cada elemento. De entre los diversos métodos que existen para evaluar la matriz de propiedades del material, en este trabajo se propone el uso de técnicas de homogeneización aceleradas mediante técnicas de machine learning. Dicha técnica se ha aplicado a problemas reales obteniendo un elevado speed up sin sacrificar la precisión.

PALABRAS CLAVE: homogeneización; espuma de aluminio; red neuronal; machine learning.

ABSTRACT

The use of resistant, rigid, low-weight materials with good both acoustic and thermal properties is very interesting in today’s industry. Among these materials, one can find aluminium foams, whose mechanical behaviour is necessary for their application. In order to obtain the geometry of an aluminium foam, several techniques can be applied, and all of them are based in the fact that information is initially obtained by a Computed Axial Tomography (CAT). One of these techniques, known as segmentation, involves a CAD being generated from an image in order to build the Finite Element (FE) model. Another option is to use techniques such as CutFEM or cgFEM, in which a certain amount of pixels, which define the properties of the material, are embedded in each element. Among the existing methods for evaluating the material properties matrix, this study proposes the use of homogenization techniques, sped up by the use of machine learning techniques. This method has been applied to real problems obtaining a high speed up, conserving precision.

KEYWORDS: homogenization; aluminium foam; neural network; machine learning.

1.    INTRODUCCIÓN

Las espumas de aluminio constituyen un nuevo material con excelentes propiedades sobre los materiales convencionales, tales como gran rigidez, baja densidad, capacidad de absorción de impactos, capacidad de absorción de ondas electromagnéticas, buena estabilidad térmica, así como gran capacidad de amortiguación y aislamiento acústico [1]. En los últimos años se han llevado a cabo diversos estudios para evaluar la idoneidad del uso de las espumas de aluminio en la industria del transporte, bien en la ámbito automotriz, o en el sector aeronáutico. Al mismo tiempo, la industrialización de algunos procesos de fabricación ha conseguido incrementar la aplicabilidad de las espumas de aluminio en la ingeniería [2]. Para el análisis mecánico de espumas de aluminio utilizando herramientas de análisis como el método de los elementos finitos, se puede partir de imágenes obtenidas mediante Tomografía Axial Computarizada (TAC). Esto permite caracterizar el material sin necesidad de ensayos experimentales. Un método tradicional para generar el modelo CAD a partir de las imágenes proporcionadas por el TAC pasa por un proceso de segmentación que define de manera explícita una frontera entre el aire y el aluminio, en este caso. La definición de dicha frontera de manera explícita implica la asunción que a partir de cierto valor de color se considera aluminio, siendo el resto aire. Este proceso es laborioso ya que en espumas de aluminio existen una infinidad de contornos, casi del tamaño del píxel, y hace muy difícil esta labor requiriendo la asistencia de personal cualificado.

Para evitar esta problemática, en una primera aproximación surgieron métodos de generación de modelos de elementos finitos que asociaban un píxel por elemento cuyas propiedades elásticas se obtenían según el color del píxel. Esta aproximación basada en la fuerza bruta, aunque sencilla y precisa, resulta ineficiente ya que obliga a usar una cantidad de elementos excesiva. Por ejemplo en una imagen 3D de 512 píxeles por dirección deberíamos resolver un problema de más de 130 millones de elementos. Para evitar esta problemática aparecen técnicas más novedosas como el CellFEM [3, 4, 5] o el cgFEM

[6, 7, 8]. Estas técnicas permiten eliminar el tedioso proceso de segmentación. Se basan en incluir cierto número de píxeles dentro de cada elemento y asociar propiedades de material a cada uno de los píxeles en función del color del mismo [5]. A partir de esta información en el cgFEM se utiliza una cuadratura de integración especial que asocia un punto de integración a cada píxel. Esta metodología, aunque precisa, resulta ineficiente ya que utiliza un número innecesario de puntos de integración. Para evitar este coste se debería realizar un proceso de homogeneización y utilizar una cuadratura de integración convencional. No obstante la realización de un proceso de homogeneización en cada uno de los elementos resultaría a todas luces ineficiente.

En aras de mejorar la eficiencia de la técnica cgFEM, en el presente trabajo se tratará de obtener un método de análisis mecánico de espumas de aluminio que proporcione resultados similares a la solución con fuerza bruta que considera un píxel por elemento, con un menor coste computacional. Para esto se utilizarán técnicas de homogeneización y machine learning aplicadas en el entorno cgFEM donde considera un conjunto de píxeles por elemento. La técnica propuesta consiste en realizar un proceso de aprendizaje previo mediante redes neuronales que permita relacionar de manera eficiente el ”color” de los píxeles (de la imagen del TAC) embebidos dentro de un elemento con la matriz de propiedades de material homogeneizada. Una vez realizado este proceso, durante la resolución del problema elástico, se leerá el color de los píxeles en cada uno de los elementos y se obtendrá a través de la red neuronal la matriz de propiedades elásticas homogeneizadas. Esto permite usar una cuadratura de integración estándar y no un punto de integración por píxel, sin tener que recurrir a realizar un proceso de homogeneización en cada uno de los elementos. El trabajo queda estructurado como sigue. Después de esta introducción, en la sección 2 se define el problema elástico lineal, en la sección 3 se describe la metodología usada: análisis de la imagen, entrenamiento y uso de la red neuronal y resolución del problema elástico. Finalmente en la sección 4 se describen los resultados obtenidos seguidos de las conclusiones en la sección 5.

2. PLANTEAMIENTO DEL PROBLEMA ELÁSTICO LINEAL

En esta sección se presenta brevemente el modelo del problema elástico lineal en 2 dimensiones. Denotamos las tensiones de Cauchy como σ, los desplazamientos como u, las deformaciones como ε, todos estos campos siendo definidos en el dominio Ω ⊂ R², con el contorno definido por ∂Ω. Las tracciones prescritas, denotadas por t son impuestas sobre la parte Γ_N del contorno, mientras los desplazamientos denotados por u¯ están prescritos sobre la parte complementaria Γ_D del contorno. Las cargas por unidad de volumen se denotan por b. El problema elástico toma la forma siguiente. Buscamos (σ,u) tal que se cumpla: admisibilidad estática:

L ^Tσ + b	=	0	en Ω	(1)
Gσ	=	t	en ΓN	(2)

donde L es el operador diferencial definido como:

L   (3)

y G es el operador proyección que proyecta el campo de tensiones en tracción sobre el contorno. El operador G es la forma matricial de la ley de Cauchy

n       o^T considerando el vector unitario normal n = n_x,n_y a Γ_N tal que:

" n G =   x 0

0 ny

ny# nx

(4)

admisibilidad cinemática:

u = u¯

en ΓD

(5)

relación constitutiva:

σ = Dε(u),  con ε(u) = Lu en Ω      (6)

donde la matriz D contiene los coeficientes de elasticidad de la ley constitutiva isotrópica lineal usual que relaciona tensión con deformación.

El problema de arriba toma la forma variacional primal:

Hallar u ∈ (V + {w}) a(u,v) = l(v)   donde:

Z

a        (7)

Z      Z

l(v) = b^TvdΩ +     t^TvdΓ

Ω      Γ_N

²

donde V = v | v ∈ ^h H¹(Ω) ⁱ , v|_ΓD = 0 y w es una cam-

po de desplazamientos particular que satisface las condiciones de contorno de Dirichlet.

2.1. Discretización en elementos finitos

Introduzcamos un esquema de discretización en elementos finitos clásico para el problema de elasticidad. El campo de desplazamientos aproximado u^h es buscado en el espacio de dimensión finita (V^h +{w}) ⊂ (V+{w}) tal que V^h se construye con las funciones de forma de elementos finitos con soporte local.

Usando el marco de referencia de Galerkin, la formulación variacional primal se transforma en:

Hallar u^h ∈ V^h V^h a(u^h,v) = l(v)   (8)

3.    METODOLOGÍA

La resolución del problema elástico (8) que se plantea en este trabajo tiene la peculiaridad que la distribución de propiedades de material varía en función del color de gris proporcionado por el TAC de la espuma de aluminio. Basándonos en la metodología cgFEM, cada uno de los elementos contendrá un conjunto de píxeles. En este trabajo vamos a tomar que cada elemento contiene 8 × 8 = 64 píxeles, ya que es el número de píxeles que toma la tecnología JPEG para realizar la compresión. En la figura 1 se muestra un ejemplo de 8x8 píxeles cuyas propiedades mecánicas se desean homogeneizar.

Figura 1: Conjunto de volúmenes de referencia (en azul) tomado de la imagen del TAC de una espuma de aluminio real (abajo). Problema de homogeneización de un volumen de referencia con sus 64 píxeles (arriba). En escala de grises se representa el módulo de Young de cada píxel de la imagen del TAC.

A esta unidad básica (volumen de referencia) se le aplica un proceso de homogeneización que permitirá extraer una matriz D cuyas propiedades sean equivalentes a las que se obtendrían tomando toda la información a nivel de píxel. Esto se hace resolviendo tres problemas ( como se explica en la sección 3.2) con la máxima precisión posible, que es utilizando un elemento finito por píxel, es decir, un modelo con 162 grados de libertad. Dicho proceso de homogeneización se debe realizar a cada elemento, dando lugar a un proceso de resolución del problema lento y tedioso.

Para ello, y teniendo la información proporcionada por la Transformada Coseno que usa el algoritmo JPEG, se entrena una red neuronal que sea capaz de relacionar la información que proporciona la compresión JPEG con la matriz D de propiedades de material homogeneizada, evitando así realizar el proceso de homogeneización en cada elemento una vez la red neuronal haya sido debidamente entrenada.

En esta sección se explican los tres ingredientes fundamentales: i) obtención de la compresión JPEG del conjunto de 64 píxeles, ii) obtención de la matriz de propiedades de material D homogeneizada a partir de los coeficientes JPEG y iii) entrenamiento de la red neuronal. Finalmente se describe el proceso que se sigue para resolver el problema global de elasticidad.

3.1. JPEG y transformada discreta de coseno (DCT)

JPEG es un estándar de compresión y codificación de imágenes fijas creado por el comité de expertos del que toma nombre, Joint Photographic Experts Group. Usa un algoritmo de compresión con pérdida (lossy compression) llamado transformada de coseno discreta ( DCT ). Esta operación matemática convierte el valor de cada píxel de la imagen al dominio de frecuencia (dominio de la transformada).

Un modelo perceptual (basado en el sistema visual humano) descarta las frecuencias más altas de variación de intensidad y tono de cada píxel, pues los coeficientes de las frecuencias más altas en el dominio de la transformada son normalmente valores más pequeños que afectan menos a la imagen. Este procesado se llama cuantificación y permite reducir el tamaño del archivo, lo cual afecta generalmente a la calidad de la imagen. La mayoría de aplicaciones de JPEG permiten al usuario controlar el factor de compresión.

La transformada de coseno discreta (DCT) descompone una imagen en una serie de sinusoides de diferentes frecuencias y amplitudes. En una imagen típica, la mayor parte de la información significante está contenida en unos cuantos coeficientes de su DCT, por lo que es usada a menudo en aplicaciones de compresión de imágenes así como en métodos espectrales para la resolución de ecuaciones en derivadas parciales[9]. La transformada de coseno discreta bidimensional de una imagen (matriz) A de N×M píxeles se define como:

Por tanto en el dominio de la frecuencia existen tantos coeficientes C_uv como píxeles tiene la imagen. Los factores de normalización α_u y α_v son meramente una convención y han sido elegidos de tal manera que el primer peso C₁₁ representa el valor medio de la imagen en su dominio (valor medio de la matriz A). Por tanto, cualquier matriz de N×M puede ser expresada a partir de una serie de funciones base sinusoidales multiplicadas por su correspondiente peso C_uv. En la figura 2 se muestran las funciones base para una imagen de 8 x 8 píxeles.

Figura 2: 64 funciones base para una imagen de 8x8 píxeles.

3.2. Homogeneización

La homogeneización es un proceso de inspiración física que consiste en substituir un material fuertemente heterogéneo por uno homogéneo equivalente. Este proceso es importante cuando se intenta, por ejemplo, estudiar deformaciones o temperaturas en cuerpos con heterogeneidades debidas a la presencia de impurezas distribuidas de una cierta manera en los mismos [10].

En cuanto a la homogeneización de propiedades elásticas de un volumen de referencia, puede ser llevada a cabo con diferentes métodos, como la imposición de condiciones de contorno de Dirichlet, Neumann y periódicas o la aplicación de las condiciones de restricción de Voigt y Reuss, todos ellos descritos en [11].

Figura 3: Ejemplo de volumen de referencia con 64 píxeles cuyo comportamiento se desea homogeneizar. En verde se representa el nuevo elemento finito obtenido.

Para la homogeneización de propiedades mecánicas de espumas de aluminio se resolverá un problema de elementos finitos aplicando las condiciones de contorno periódicas, pues se ha comprobado que en un material con inclusiones repartidas periódicamente en su dominio como el de la figura 4, el comportamiento del elemento (volumen de referencia) obtenido al resolver el problema de homogeneización es independiente de la posición de estas inclusiones dentro del volumen de referencia elegido. Esto no ocurre, por ejemplo, con la aplicación de condiciones de contorno de Dirichlet [11].

Figura 4: Matriz contínua (blanco) con inclusiones periódicamente distribuidas (amarillo). En negro se representa el volumen de referencia tomado para realizar el proceso de caracterización del material. Inclusión céntrica (izquierda) y excéntrica (derecha) dentro del volumen de referencia a homogeneizar.

3.2.1.    Condiciones de contorno periódicas

Para la homogeneización de propiedades elásticas del volumen de referencia, se resuelven tres problemas con deformación promedio unitaria (ε¯_x = 1, ε¯_y = 1 y γ¯_xy = 1). Los desplazamientos en los contornos enfrentados (inferiorsuperior, izquierdo-derecho) tendrán la misma variación a lo largo del contorno, ya que se aplican condiciones de contorno periódicas.

Tomando como esquema del volumen de referencia a optimizar la figura 5, se restringen los dos grados de libertad del nodo 1, u = 0,v = 0 (en los tres casos de deformación), y a partir de aquí se pueden prescribir los desplazamientos del nodo 9 como:

donde u denota el desplazamiento horizontal del nodo y v el desplazamiento en dirección vertical.

Figura 5: Esquema del conjunto 8x8 elementos finitos cuadriláteros lineales(uno por píxel) y nodos correspondientes.

Se puede de igual manera prescribir los desplazamientos de los nodos 73 y 81. En total, se restringen 8 desplazamientos correspondientes a los cuatro nodos esquina del conjunto de elementos.

En cuanto a los 7 nodos interiores de cada contorno, para garantizar la aplicación de las condiciones de contorno periódicas se deben formular 56 restricciones de la siguiente forma:

Se utilizará el método de los multiplicadores de Lagrange para encontrar la solución a las ecuaciones que minimice el funcional Π_p = u^TKu − u^TF, obtenido a partir de

(8). Se obtienen 56 ecuaciones. A modo de ejemplo, si se toma como nodo a optimizar el 10 se obtiene:

Estas ecuaciones se organizan de la siguiente manera:

Por tanto, para cada caso de deformación media, se construyen las matrices C y el vector g y se resuelve el sistema mostrado en la ecuación (14).

3.3. Redes neuronales

Una red neuronal es un procesador distribuido en paralelo que tiene la capacidad natural de almacenar conocimientos experimentales y hacerlos disponibles para su uso [12].

El aprendizaje es un proceso por el cual los parámetros libres de la red neuronal se adaptan mediante la estimulación del entorno en el que está embebida la red. El presente trabajo se centrará en el uso de redes neuronales entrenadas mediante aprendizaje supervisado. Éste asume un conjunto de datos de aprendizaje consistentes en N pares de datos de entrada-salida (objetivo):

T = (xi,di)iN=1             (15)

donde x_i es el vector de datos de entrada i-ésimo, d_i es el vector deseado de salida (objetivo) al vector de entrada x_i, y N es el tamaño de la muestra.

Por tanto, se desea, dada una muestra T, calcular los parámetros libres de la red neuronal para reducir una función error que mide la diferencia entre la salida y _i a un vector de entrada x_i y la salida deseada d_i para todos los i de una manera estadística. En el presente trabajo se ha utilizado el error cuadrático medio (ECM) [13]:

En la figura 6 se muestra un ejemplo de una red neuronal prealimentada con p elementos de entrada, una capa de neuronas oculta consistente de 3 neuronas de tipo función logística, y una única salida.

Figura 6: Ejemplo de red neuronal prealimentada con p entradas, una capa oculta de N (en este caso tres) neuronas y S (una única) salidas.

Se observa que el número de conexiones (sinapsis) entre la capa de entrada y la capa oculta es de Np, siendo N el número de neuronas ocultas. Por tanto, en esta primera sinapsis los coeficientes libres a optimizar serán N · p coeficientes llamados pesos sinápticos w, así como N niveles de activación b (uno por cada neurona de la capa oculta).

Los coeficientes libres se organizan de forma matricial, y se obtiene (17), donde sigmoide es la función de transferencia sigmoidea:

aNx1 = sigmoide(W ₁ Nxp · xpx1 + b₁Nx1)        (17)

En la segunda sinapsis entre la capa oculta y la capa de salida, los coeficientes libres son N · s pesos sinápticos w, así como s constantes; sea s el número de salidas ( en la figura 6 sólo una).

La salida se obtiene de la siguiente manera:

ysx1 = W₂sxN · aNx1 + b₂sx1  (18)

Por tanto, en una red neuronal prealimentada con una única capa oculta, el número de coeficientes libres a ajustar es N · p + N + N · s + s = N · (p + s + 1) + s, es decir, el número de coeficientes a ajustar aumenta de forma lineal con el número de neuronas de la capa oculta en este ejemplo, donde N es el número de neuronas de la capa oculta, y p y s son el números de elementos de cada par de datos de entrada y salida respectivamente.

3.3.1.    Entradas. Matriz C

Una imagen de 64 píxeles como en el caso que nos ocupa, tiene 64 pesos de su transformada DCT (Discrete Cosine Transform). En este primer análisis utilizaremos solamente tres de estos coeficientes (C₁₁, C₁₂ y C₂₁) que representaran la media y la información contenida en la frecuencia más baja en sentido horizontal y vertical. Para reducir el número de entradas, se utilizarán como entradas C₁₂/C₁₁ así como C₂₁/C₁₁, y la matriz de constantes elásticas resultante se multiplicará por el peso C₁₁ para obtener la matriz buscada. Por tanto, la capa de entrada tendrá 2 neuronas.

3.3.2.    Salidas. Matriz D

En cuanto a la salida, se buscarán 6 coeficientes de la matriz D (D₁₁, D₁₂, D₁₃, D₂₂, D₂₃ y D₃₃), al ser ésta simétrica. Por tanto, la capa de salida tendrá 6 neuronas.

3.3.3.    Aprendizaje

Desde que el algoritmo de propagación hacia atrás (Backpropagation algorithm) fue obtenido por Henry J. Kelley en 1960, han aparecido numerosos métodos para acelerar la convergencia del algoritmo. Estos se pueden clasificar en técnicas ad hoc y métodos quasi-Newton.

Los métodos quasi-Newton se consideran más eficientes, si bien los requisitos computacionales aumentan con el cuadrado del tamaño de la red neuronal [12]. Sin embargo, en el presente trabajo, la red neuronal apenas tendrá unos cuantos coeficientes libres, por lo que este segundo criterio será el adecuado. El algoritmo de MarquardtLevenberg usado por el software está descrito en [14]. La red neuronal utilizada en Matlab tendrá por defecto un aprendizaje de tipo Marquardt-Levenberg con ratios de entrenamiento, validación y test serán del 70% - 15% - 15% de los pares de datos de entrada-salida [14]. En la figura 7 se muestra un esquema de la red neuronal generada, con una capa oculta compuesta de 4 nodos o neuronas cuya función de transferencia es en este caso la tangente hiperbólica.

Para la resolución de un problema de elementos finitos mediante la técnica propuesta, una vez entrenada la red neuronal se seguirá el siguiente esquema:

1. Se agrupan los píxeles en conjuntos de 8×8.

2. Se aplica la transformada DCT al conjunto de píxeles para obtener los 64 coeficientes C_uv.

Figura 7: Esquema de la red neuronal generada.

3. Se evalúa la red neuronal con los ratios {C₁₂/C₁₁ y C₂₁/C₁₁} para obtener la matriz D del nuevo elemento. Este proceso es mucho más rápido que resolver el problema de homogeneización.

4. Se obtiene con ella la matriz de rigidez K_e del elemento.

5. Ensamblado de las matrices de rigidez. Obtención de KGlobal

6. Obtención del vector de fuerzas nodales f adaptado a los nodos de los nuevos elementos.

7. Resolución del sistema de ecuaciones.

4. RESULTADOS

En este apartado se ha probado el uso de las redes neuronales para acelerar el proceso de homogeneización. En primer lugar, se ha realizado un análisis para escoger la configuración de red neuronal más apropiada para la homogeneización de propiedades elásticas de 8x8 píxeles de una imagen de espumas de aluminio (sección 4.1). A continuación, se ha resuelto un problema a partir de una imagen real (sección 4.2).

4.1. Caracterización mecánica de espumas de aluminio usando 3 pesos de la transformada DCT

En esta sección se ha tratado de escoger la configuración de la red neuronal más ventajosa para la caracterización mecánica de espumas metálicas. Para esto diferentes configuraciones han sido analizadas, variando el número de capas ocultas así como el número de neuronas en cada capa.

Figura 8: Coeficiente D₁₁. Real vs red neuronal con 1 neurona en la capa oculta.

4.1.1. Entrenamiento A: Muestras de todo el espacio de entrada. Influencia del número de neuronas en la capa oculta y número de muestras sobre el ECM

En la Tabla 1 se muestra el error cuadrático medio (ECM) obtenido en todo el espacio {C₁₂/C₁₁,C₂₁/C₁₁} para varias redes neuronales con una capa oculta y nodos desde 1 hasta 8, así como dos capas ocultas y nodos desde 1/1 hasta 8/8. Como entradas a la red neuronal se proporcionan 100 conjuntos entrada-salida. Las mismas salidas (deseadas) han sido comparadas con la salida de la red neuronal para hallar el ECM. Cada caso se ha repetido 3 veces y se ha hallado la media del ECM.

Tabla 1. Influencia del Número de Capas Ocultas y Neuronas.

1 Capa Oculta		2 Capas Ocultas
Nodos ECM		Nodos ECM
1	1.979 ·10^−3	1/1	2.099 ·10^−3
2	1.237 ·10^−3	2/2	1.136 ·10^−3
3	6.773 ·10^−4	3/3	3.115 ·10^−4
4	8.313 ·10^−6	4/4	1.500 ·10^−6
5	3.161 ·10^−6	5/5	1.218 ·10^−6
6	2.497 ·10^−6	6/6	1.425 ·10^−6
7	2.594 ·10^−6	7/7	2.029 ·10^−6
8	2.211 ·10^−6	8/8	4.312 ·10^−6

Como vemos en la Tabla 1, a partir de 4 neuronas, el rendimiento de la red neuronal es muy bueno, y el ECM es de un orden despreciable.

En segundo lugar, se observa las ventajas del uso de una red neuronal de 1 sola capa oculta ya que con menos parámetros libres y un menor tiempo de cálculo se ha obtenido una precisión semejante a modelos con más capas ocultas.

Como se ve en las figuras 8, 9, 10, 11 y 12, para las salidas D₁₁, D₁₂, D₂₂ y D₃₃, una capa oculta con una, dos o tres neuronas no modela bien las funciones deseadas de los coeficientes de salida. Sin embargo, a partir de cuatro neuronas en la capa oculta la red neuronal es capaz de modelar la curva de manera adecuada. Como se ve en las figuras 13, 14 y 15, para las salidas D₁₃ y D₂₃, el uso de redes neuronales con una capa oculta de una, dos o tres neuronas tampoco simula correctamente las funciones de los coeficientes de salida. Para cuatro neuronas, la red neuronal, como se ve en la figura 16, es capaz de proporcionar una función similar a la del coeficiente de salida, mientras que para un número mayor de neuronas, no se ha obtenido una mayor precisión ( figura 17). Esto puede ser debido en parte a un exceso de coeficientes libres, además de al hecho de que el valor de los coeficientes D₁₃ y D₂₃ oscila en un rango de valores pequeño y muy cercano a cero (en el caso de materiales isótropos así como ortótropos, D₁₃ y D₂₃ es cero por definición, y en este caso no estamos considerando variaciones cruzadas de las propiedades del material). Por tanto, la red neuronal trata de ajustar un valor prácticamente nulo, y se produce un error de tipo numérico y sin una gran aportación al error cuadrático medio total.

Figura 9: Coeficiente D₁₁. Real vs red neuronal con 2 neuronas en la capa oculta.

Figura 10: Coeficiente D₁₁. Real vs red neuronal con 3 neuronas en la capa oculta.

4.1.2. Entrenamiento B: Muestras aleatorias. Influencia del número de neuronas en la capa oculta y número de muestras sobre el ECM

En esta sección, se han generado conjuntos de 8x8 píxeles aleatorios, para después obtener los pesos de la transformada DCT (entradas) y generar conjuntos de entradasalida para el entrenamiento. El error cuadrático medio se ha evaluado en todo el espacio al igual que en el caso anterior. Cada caso se ha evaluado tres veces.

En la Tabla 2 se observa que para un número de neuronas en la capa oculta mayor, la mayor precisión se alcanza con un mayor número de muestras.

Figura 11: Coeficiente D₁₁. Real vs red neuronal con 4 neuronas en la capa oculta.

Figura 12: Coeficiente D₁₁. Real vs red neuronal con 5 neuronas en la capa oculta.

Como se observa en la Tabla 2, en el caso de 3 pesos de la DCT, el entrenamiento con muestras aleatorias ha resultado menos eficiente que el realizado con muestras de todo el espacio muestral. Las muestras generadas aleatorias tienen más contenido en altas frecuencias de variación del color que las imágenes reales. Estas frecuencias no se toman en cuenta en nuestro análisis. Sin embargo, como se ve en las figuras 18 y 19, las muestras aleatorias tienden a presentar unas variaciones de color del píxel en las bajas frecuencias más suaves (C₁₂/C₁₁, C₂₁/C₁₁), por lo que las muestras de entrenamiento no cubren todo el espacio muestral y no son apropiadas para entrenar la red neuronal.

Se ha observado que el ECM crece por tanto, en las regiones donde los pesos asociados a las bajas frecuencias (C₁₂/C₁₁, C₂₁/C₁₁) son más altos.

4.2.   Imagen real

En esta sección se ha planteado un problema de elementos finitos, modelando una espuma de aluminio a partir de una imagen 2D de una tomografía axial computarizada (TAC). Las soluciones de referencia y la obtenida con el uso de redes neuronales han sido comparadas en cuanto a coste computacional y precisión de los resultados.

Figura 13: Coeficiente D₁₃. Real vs red neuronal con 1 neurona en la capa oculta.

Figura 14: Coeficiente D₁₃. Real vs red neuronal con 2 neuronas en la capa oculta.

4.2.1.    Planteamiento del problema

En la figura 20 se muestra una imagen de 512x512 píxeles obtenida con tomografía axial computarizada para obtener una solución de referencia. El modelo ha sido construido asociando a cada píxel un elemento finito cuyo módulo de Young varía linealmente con el color del píxel. El color más claro ha sido modelizado con E = 72GPa, mientras que el color más oscuro, asociado a los huecos de la espuma, ha sido modelizado con E = 1000Pa. El coeficiente de Poisson ha sido supuesto de 0.3 en todo el dominio de la geometría.

Figura 15: Coeficiente D₁₃. Real vs red neuronal con 3 neuronas en la capa oculta.

Figura 16: Coeficiente D₁₃. Real vs red neuronal con 4 neuronas en la capa oculta.

Se ha tomado como dimensiones de la espuma 1 × 1 unidades de longitud ya que esto no influye sobre el análisis de datos. El problema a resolver consta en la aplicación de una fuerza distribuida de forma lineal sobre la cara superior de la espuma. Se han restringido los desplazamientos verticales en la cara inferior así como los horizontales en el nodo inferior izquierdo para evitar movimientos de sólido rígido, tal y como se muestra en la figura 21.

4.2.2.    Modelo 0. Solución de referencia

En la figuras 22, 23 y 24 se representa la distribución de deformaciones para la solución de referencia (asociando un elemento a cada píxel).

Esta solución, por tanto, ha sido la utilizada para validar el modelo con redes neuronales. En la Tabla 3 se ha calculado la energía de deformación de referencia.

4.2.3.    Modelo 1. Solución con 1 coeficiente JPEG

En general el tratamiento de imágenes píxel a píxel no es viable. Para el caso en 2D, se han tenido que resolver 262.144 elementos. Para analizar un TAC 3D de resolución 512 píxeles por dimensión se requerirían un total de 134.217.728 elementos lo que hace esto inviable para los PC’s actuales de uso doméstico.

Figura 17: Coeficiente D₁₃. Real vs red neuronal con 5 neuronas en la capa oculta.

Figura 18: Coeficiente D₁₁. Distribución de ECM con 4 neuronas en la capa oculta tomando 50 muestras aleatorias.

Una de las técnicas para aliviar esta dificultad es la de utilizar el promedio del valor del color en una malla cuyos elementos contengan gran cantidad de píxeles. A continuación se ha resuelto el problema agrupando los píxeles en elementos de 8x8. A cada nuevo elemento le ha sido asociado el módulo de elasticidad medio del conjunto de 8x8 píxeles, es decir, el primer peso de la compresión JPEG, como se muestra en la imagen 25. En las figuras 26, 27 y 28 se representa la solución obtenida.

Como se observa, la homogeneización de constantes elásticas y posterior discretización del problema en elementos cuadriláteros de primer orden de 8x8 píxeles no permite obtener valores máximos de las deformaciones. Esto es debido, entre otros motivos, a que la aproximación se realiza en una malla más basta y con una reducción considerable del coste computacional (con nuestro portátil Asus K53SJ, y con nuestro código, el speed up fue de 240 veces). Esto provoca que la energía de deformación se subestime en todos los casos analizados, lo cual no nos pone en el lado de la seguridad.

Figura 19: Coeficiente D₁₁. Distribución de ECM con 4 neuronas en la capa oculta tomando 400 muestras aleatorias.

Tabla 3. Solución de Referencia

Energía de deformación (MJ)
Solución de referencia	1.3673

mente la media del color en cada elemento formado por 8x8 píxeles,se han considerado también las variaciones horizontales y verticales del color, a través de los coeficientes C₁₂ y C₂₁. Con respecto a la anterior sección esto no lleva ningún coste computacional en la fase online ya que las matrices D están precalculadas gracias al entrenamiento en la fase offline de la red neuronal. Notar que

Figura 22: Solución de referencia. Deformación en dirección x (ε_x).

Figura 23: Solución de referencia. Deformación en dirección y (ε_y).

este entrenamiento es independiente del problema y se ha realizado una única vez. El objetivo final es obtener una precisión mayor que el método anterior manteniendo el coste computacional gracias al uso de una red neuronal. Se ha usado para ello una red neuronal con una capa oculta de 4 neuronas, que ha sido entrenada con N muestras, que o bien han sido tomadas de la imagen (entrenamiento A), o generadas como un conjunto de 8x8 píxeles aleatorios (entrenamiento B).

Entrenamiento A

En este primer modelo de entrenamiento, se han elegido

N conjuntos de 8x8 píxeles de la imagen. Los resultados en cuanto a energía de deformación del problema se muestran en la Tabla 5. Cada caso ha sido evaluado tres veces y se ha hallado la media de la energía de deformaFigura 24: Solución de referencia. Deformación angular

(γ_y).

Figura 25: Imagen resultante tras la compresión JPEG con el primer peso.

ción.

Tabla 5. Energ´ia de Deformacion´      (MJ). Influencia del Numero´    de Muestras.

Energía de deformación ( MJ )
N=10	1.2238 (-10.50%)
N=50	1.2110 (-11.43%)
N=100	1.2485 (-8.69%)
N=500	1.3217 (-3.33%)

Como vemos en la Tabla 5, tomando 3 coeficientes de la DCT, el número de muestras óptimo ha sido el mayor posible, ya que ello permite cubrir el mayor espacio muestral posible.

Al igual que en la sección anterior, todos los experimentos tienden a subestimar la energía de deformación. En las figuras 29, 30 y 31 se representa las distribuciones

Figura 28: Solución con 1 peso de la DCT. en cuanto a la energía de deformación del problema se muestran en la Tabla 6. Cada caso ha sido evaluado tres veces y se ha hallado la media de la energía de deforma-

de deformaciones, para 3 coeficientes de la DCT y 500 ción.

Figura 31: Solución con 3 pesos de la DCT y 500 muestras.

Tabla 6. Energía de Deformación (MJ). Influencia del Número de Muestras.

Energía de deformación (MJ)
N=10	1.1908 (-12.91%)
N=50	1.1913 (-12.87%)
N=100	1.1929 (-12.76%)
N=500	1.1966 (-12.48%)

Como era previsible, este método para hallar los coeficientes ha resultado menos eficiente ya que no controlamos directamente el espacio de muestra de nuestras magnitudes de interés (C₁₁, C₁₂ y C₂₁), por ello no se ha conseguido mejorar los resultados con respecto al método clásico (media del color).

4.2.5.    Modelo 3. Solución con 6 coeficientes JPEG

En esta sección se resolverá el problema tomando 6 pesos de la DCT.

Entrenamiento A

Los resultados en cuanto a la energía de deformación del problema se muestran en la Tabla 7.

Tabla 7. Energía de deformación (MJ). Influencia del número de muestras.

Energía de deformación (MJ)
N=10	1.2925 (-5.47%)
N=50	1.2341 (-9.74%)
N=100	1.2500 (-8.58%)
N=500	1.2711 (-7.04%)

En este caso, la inclusión de más parámetros no ha aportado una mejora en la calidad de la solución. Este hecho podría ser debido a la necesidad de un tamaño muestral mayor. No obstante, si el elemento usado es lineal, la inclusión de variaciones de material de orden superior podría no tener influencia alguna en los resultados, o a lo sumo contribuir en una fuente de error numérico.

Entrenamiento B

Los resultados en cuanto a la energía de deformación del problema se muestran en la Tabla 8.

Tabla 8. Energ´ia de Deformacion´      (MJ). Influencia del Numero´    de Muestras.

Energía de deformación ( MJ )
N=10	1.1955 (-12.56%)
N=50	1.1978 (-12.40%)
N=100	1.2052 (-11.86%)
N=500	1.1971 (-12.45%)

Al igual que en el caso con 3 coeficientes JPEG, se han obtenido valores parecidos a los obtenidos con el método clásico. Esto es debido a que el entrenamiento con muestras aleatorias es menos eficiente y solo consigue capturar con precisión los valores medios.

5.    CONCLUSIONES

Los resultados mostraron que el uso de redes neuronales puede agilizar el proceso de homogeneización de propiedades elásticas de la imagen en conjuntos de 8x8 píxeles. Además, se ha observado que las frecuencias bajas en las propiedades del material son aquellas que contienen una mayor información de interés para evaluar las propiedades elásticas del elemento.

En cuanto al uso de 6 coeficientes JPEG, no se ha observado una mejora respecto al método con 1 coeficiente JPEG, confirmando la hipótesis anterior.

En el caso 3D el uso de mallas a nivel de píxel es inviable y la mejora de las técnicas de homogenización existentes es una necesidad evidente. Como trabajo futuro se plantea la implementación del método para una imagen en 3 D.

En definitiva, en el presente trabajo se ha conseguido plantear un método de homogeneización de las propiedades elásticas a partir de imágenes que permitirá incrementar la eficiencia de la metodología cgFEM en el análisis de estructuras cuya geometría se obtiene a partir de imágenes.

AGRADECIMIENTOS

Los autores agradecen la ayuda recibida por el Ministerio de Economía y Competitividad a través del proyecto DPI2013-46317-R así como el apoyo recibido por la Generalitat Valenciana mediante el proyecto Prometeo/2016/007.

REFERENCIAS

[1] C. Meola, S. Boccardi, and G. M. Carlomagno, “Composite materials in the aeronautical industry,” Elsevier Ltd, 2015.

[2] J. Banhart, “Aluminum foams for transport industry,” Materials and Design, vol. 18, pp. 217–220, 1997.

[3] J. Parvizian, A. Düster, and E. Rank, “Finite cell method,” Computational Mechanics, vol. 41, no. 1, pp. 121–133 , sep 2007. [Online]. Available: http://link.springer.com/10.1007/ s00466-007-0173-y

[4] A. Düster, J. Parvizian, Z. Yang, and E. Rank, “The finite cell method for three-dimensional problems of solid mechanics,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 45-48, pp. 3768–3782, aug 2008. [Online]. Available: http://www.sciencedirect.com/science/article/pii/ S0045782508001163

[5] L. Giovannelli, J. J. Ródenas, J. M. Navarro, and M. Tur, “Direct medical image-based finite element modelling for patient-specific simulation of future implants,” Elsevier Ltd, 2017.

[6] E. Nadal, “Cartesian grid FEM (cgFEM): High performance hadaptive FE analysis with efficient error control. Application to structural shape optimization,” Ph.D. dissertation, Universitat Politècnica de València, 2014.

[7] M. Tur, J. Albelda, O. Marco, and J. J. Ródenas, “Stabilized method of imposing Dirichlet boundary conditions using a recovered stress field,” Computer Methods in Applied Mechanics and Engineering, vol. 296, pp. 352–375, aug 2015. [Online]. Available: http://www.sciencedirect.com/science/article/pii/ S0045782515002467

[8] M. Tur, J. Albelda, E. Nadal, and J. J. Ródenas, “Imposing Dirichlet boundary conditions in hierarquical cartesian meshes by means of stabilized Lagrange multipliers,” International Journal for Numerical Methods in Engineering, vol. Accepted, 2013.

[9] [Online].    Available:  https://es.mathworks.com/help/images/ discrete-cosine-transform.html

[10] M. Lobo and M. E. Pérez, “Problemas de homogeneización en la ingeniería. Una experimentación numérica,” Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, vol. 7, no. 4, pp. 455–472, 1991.

[11] J. Schröder, “Constraint/boundary conditions on the microscale,” A numerical two-scale homogenization scheme: the FE²method, Udine, 2014.

[12] M. Hagan and M. B. Menhal, “Training feedforward networks with the marquardt algorithm,” IEEE transactions on neural network, vol. 5, no. 6, pp. 989–993, 1994.

[13] I. W. Sandberg, “Feedworward neural networks: An introduction,” Nonlinear Dynamical Systems: Feedforward Neural Network Perspectives, Wiley and Sons, 2001.

[14] H. Damuth and M. Beale, “Faster training,” Neural Network Toolbox User Guide, The MathWorks, Inc, 2002.