## Abstract

This work analyzes the influence of the discretization error contained in the Finite Element (FE) analyses of each design configuration proposed by the structural shape optimization algorithms over the behavior of the algorithm. The paper clearly shows that if FE analyses are not accurate enough, the final solution provided by the optimization algorithm will neither be optimal nor satisfy the constraints. The need for the use of adaptive FE analysis techniques in shape optimum design will be shown. The paper proposes the combination of two strategies to reduce the computational cost related to the use of mesh adaptivity in evolutionary optimization algorithms: (a) the use of the algorithm described by Bugeda et al. [1] which reduces the computational cost associated to the adaptive FE analysis of each geometrical configuration and, (b) the successive increase of the required accuracy of the FE analyses in order to obtain a considerable reduction of the computational cost in the early stages of the optimization process.

## Palabras clave

Optimización de forma estructural ; Estimación de error ; Refinamiento adaptativo ; Análisis de sensibilidades ; Algoritmos evolutivos

## Keywords

Structural shape optimization ; Error estimation ; Adaptive remeshing ; Sensitivity analysis ; Evolutionary algorithms

## 1. Introducción

Desde el punto de vista matemático, un problema de optimización se puede considerar como la minimización de una función f (v ) que depende de un conjunto de variables v y que está sujeta a una serie de restricciones. La forma general de un problema de este tipo es:

 ${\displaystyle {\begin{array}{lrcll}minimizar:&f({\boldsymbol {\mbox{v}}});{\boldsymbol {\mbox{v}}}&=&\lbrace v_{i}\rbrace &;i=1,\ldots ,n\\con:&{\boldsymbol {\mbox{g}}}({\boldsymbol {\mbox{v}}})&=&\lbrace g_{j}({\boldsymbol {\mbox{v}}})\rbrace &;j=1,\ldots ,m\\verificando:&g_{j}({\boldsymbol {\mbox{v}}})&\leq &0&;j=1,\ldots ,m\\&a_{i}&\leq &v_{i}\leq b_{i}&;i=1,\ldots ,n\\&&&&\end{array}}}$
( 1)

donde f   es la función objetivo (FO), ${\textstyle v_{i}}$ son las variables de diseño y gj son las restricciones en desigualdad, las cuales, en los problemas estructurales, normalmente se expresan en función de las tensiones y/o de los desplazamientos. Los valores ai y bi definen restricciones laterales. Cada individuo se caracteriza por un conjunto de valores de v que se corresponden con un diseño estructural específico. La definición de cada diseño en función de los valores de ${\textstyle v_{i}}$ se denomina la parametrización del problema de optimización. La resolución del problema de diseño óptimo consiste en hallar los valores de v que definen el mejor diseño.

Los algoritmos utilizados para la resolución de los problemas de optimización son, en general, iterativos. Con independencia del algoritmo utilizado, será necesario calcular los valores de f y g para cada uno de los distintos diseños propuestos durante el proceso iterativo de optimización. En este trabajo hemos considerado problemas de optimización de formas estructurales. Para este tipo de problemas los valores de f y g se obtienen normalmente mediante la utilización del método de los elementos finitos (MEF). Por tanto, es necesario construir una malla específica para cada uno de los distintos diseños a analizar y, a continuación, utilizar el MEF para obtener la respuesta estructural de cada diseño. Llegados a este punto se deben tener en cuenta dos aspectos principales, ligados a la evaluación de f y g mediante el MEF, que tienen una gran incidencia sobre el comportamiento global del proceso de optimización: el esfuerzo computacional requerido para la evaluación numérica de cada individuo (configuración geométrica) y la precisión de los resultados del MEF.

La importancia del esfuerzo computacional requerido para la evaluación de cada configuración geométrica es evidente. En los problemas de optimización como los que estamos considerando, la mayor parte del coste computacional está dedicado al análisis de los distintos individuos para la determinación de los valores de la FO y del grado de satisfacción de las restricciones.

Por otro lado, y en cuanto a la precisión de los resultados, se debe tener en cuenta que los métodos numéricos de análisis solo proporcionan valores aproximados para los datos que los algoritmos de optimización precisan (FO, restricciones,…). Si estos valores no son lo suficientemente precisos se puede introducir una cantidad de ruido excesivo en el proceso de optimización. Esto puede empobrecer la velocidad de convergencia del proceso de optimización, conducir la convergencia hacia una solución no óptima o incluso producir una solución final no factible que incumpla alguna de las restricciones. En el contexto del MEF, los métodos h -adaptables, p -adaptables y hp -adaptables se pueden utilizar para obtener soluciones con un nivel de precisión prescrito. No obstante, la utilización de estas técnicas implica un gran coste computacional que reduce la eficiencia computacional del proceso de optimización.

Este artículo mostrará que si no se garantiza una mínima calidad en los resultados del análisis de cada diseño, que son los que se utilizan para alimentar al algoritmo de optimización, no se puede asegurar un comportamiento correcto del proceso de optimización. Para ello, se mostrará la incidencia del valor máximo prescrito de la norma energética del error de discretización sobre los resultados finales obtenidos mediante un algoritmo evolutivo de optimización.

## 2. Motivación mediante el estudio de un caso sencillo: sección transversal de un tubo

 Figura 1. Sección transversal de un tubo. Datos para el análisis, modelo original, solución óptima analítica y variables de diseño.

La máxima tensión de von Mises admisible en el modelo se ha restringido a 2,0 × 106 .

Es bien conocido que la geometría analítica óptima del contorno exterior se corresponde con una forma circular. Las siguientes ecuaciones muestran los valores exactos de los desplazamientos y las tensiones para un punto dado en función de sus coordenadas cartesianas (x , y ), con k  = R0 /Ri , ${\textstyle r={\sqrt {x^{2}+y^{2}}}}$ y ϕ  = arctg(y /x ) situado en un cilindro de pared gruesa sujeto a una presión interior P . Por tanto, estas ecuaciones se pueden utilizar para evaluar la solución analítica de este problema de optimización:

 ${\displaystyle u_{r}={\frac {P(1+\nu )}{E(k^{2}-1)}}\left[(1-2\nu )r+{\frac {R_{0}^{2}}{r}}\right]}$
( 2)

- Tensiones en coordenadas cilíndricas y cartesianas:

 ${\displaystyle {\sigma }_{r}={\frac {P}{k^{2}-1}}\left(1-{\frac {R_{0}^{2}}{r^{2}}}\right)}$ ${\displaystyle {\sigma }_{t}={\frac {P}{k^{2}-1}}\left(1+{\frac {R_{0}^{2}}{r^{2}}}\right)}$ ${\displaystyle {\sigma }_{xx}={\sigma }_{r}{cos}^{2}(\phi )+{\sigma }_{t}{sin}^{2}(\phi )}$ ${\displaystyle {\sigma }_{yy}={\sigma }_{r}{sin}^{2}(\phi )+{\sigma }_{t}{cos}^{2}(\phi )}$ ${\displaystyle {\sigma }_{xy}=({\sigma }_{r}-{\sigma }_{t})sin(\phi )cos(\phi )}$ ${\displaystyle {\sigma }_{z}=2\nu {\frac {P}{k^{2}-1}}}$
( 3)

El radio externo de la solución analítica para este problema (1/4 de la sección transversal) evaluado a partir de (3) es R0  = 10,67033824461. Esto corresponde a un área Aopt  = 69,787307715081.

### 2.1. Algoritmo evolutivo

Para estos análisis numéricos se ha utilizado el algoritmo evolutivo de optimización denominado Differential Evolution (DE). DE es un algoritmo robusto que ha proporcionado buenos resultados en su aplicación a problemas de muy distinto tipo. DE fue desarrollado por Storn y Price [5] . Su idea fundamental es el operador diferencial, que obedece al mismo propósito que el parámetro de cruce en los algoritmos genéticos estándar, es decir, el intercambio de información entre individuos para producir las nuevas generaciones. De las dos versiones diferentes propuestas por Storn y Price para este algoritmo se ha escogido la DE1 para este trabajo.

Los valores iniciales de las variables de diseño y sus restricciones tanto laterales como geométricas para este problema se muestran en la tabla 1 .

Tabla 1. Sección transversal de un tubo. Valores de las variables de diseño.
Variable de diseño Valor inicial Rango Restricciones
V1 20 [5,2-50,0]
V2 19 [4,0-50,0]
V3 19 [4,0-50,0] V3  < V1  − 0,5
V4 20 [5,2-50,0] V4  > V2  + 0,5

Obsérvese que las restricciones entre los valores de las variables de diseño se han utilizado para minimizar la aparición de geometrías con formas extrañas.

### 2.2. Análisis mediante elementos finitos h -adaptables

 ${\displaystyle \Vert {\boldsymbol {\mbox{u}}}{\Vert }^{2}={\int }_{\Omega }{\boldsymbol {\sigma }}^{T}{\boldsymbol {\mbox{D}}}^{-1}{\boldsymbol {\sigma }}\quad d\Omega }$
( 4)

Por tanto, se utiliza la siguiente expresión para evaluar la estimación del error en la norma energética, η :

 ${\displaystyle {\eta }^{2}={\int }_{\Omega }{\left({\boldsymbol {\sigma }}^{\ast }-{\boldsymbol {\sigma }}^{h}\right)}^{T}{\boldsymbol {\mbox{D}}}^{-1}\left({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h}\right)d\Omega }$
( 5)

donde σh es el campo de tensiones obtenido directamente del análisis mediante el MEF, σ* es un campo de tensiones mejorado, D relaciona las deformaciones con las tensiones mediante σ  = Dɛ y Ω es el dominio de análisis. Para cuantificar la tolerancia en el nivel del error de discretización presente en un determinado análisis se ha usado el error relativo en norma energética γ definido como:

 ${\displaystyle \gamma ={\frac {\eta }{\Vert {\boldsymbol {\mbox{u}}}\Vert }}\approx {\frac {\eta }{\sqrt {\Vert {\boldsymbol {\mbox{u}}}^{h}{\Vert }^{2}+{\eta }^{2}}}}}$
( 6)

donde ||uh || es el valor de norma energética proporcionado por la solución de EF.

En los ejemplos numéricos se han utilizado elementos triangulares cuadráticos y la técnica del alisado global por mínimos cuadrados [6] para obtener el campo de tensiones mejorado requerido por el estimador de error de Zienkiewicz y Zhu. Las tensiones reconstruidas en nodos del contorno obtenidas mediante esta técnica han sido utilizadas para verificar las restricciones de tensión en el modelo. También se podrían utiliza otras técnicas como el Superconvergent Patch Recovery (SPR) [7] , [8]  and [9] o cualquiera de sus mejoras posteriores [10] , [11] , [12]  and [13] .

Para definir cada una de las configuraciones geométricas en función de las coordinadas de algunos puntos de definición (parametrización) se han utilizado B -splines cúbicos [14]

Para ver el efecto de la magnitud del error de discretización en norma energética contenido en el análisis de cada diseño sobre el comportamiento de los algoritmos de optimización se han estudiado 6 situaciones distintas correspondientes a diferentes niveles de error de discretización prescrito γ . Los valores del error prescrito para los primeros 5 análisis han sido 1, 2,5, 5, 10 and 20%. En el sexto caso se ha especificado una tolerancia γ  < 100 % lo cual, en la práctica, implica que no existe ningún control sobre la calidad de la solución y que el número de elementos generados para cada malla solo depende de criterios geométricos (ver los ejemplos en la figura 2 ).

### 2.3. Optimización utilizando el algoritmo evolutivo

La figura 3 muestra el efecto del valor máximo prescrito para la estimación del valor relativo de la norma energética del error sobre la evolución del área de la sección transversal del tubo para cada uno de los procesos de optimización puestos en marcha. Se puede observar que, al menos para este problema, todos los casos estudiados presentan aspectos globales bastante similares entre sí.

 Figura 3. Influencia de γ sobre la evolución de la función objetivo.

No obstante, y tal como se muestra en la figura 4 , el efecto de γ sobre la evolución del proceso de optimización es especialmente significativo si se comparan los resultados finales obtenidos para cada grado del error de discretización. La figura 4 muestra la evolución de la diferencia relativa entre el área proporcionada por el mejor individuo de cada generación y el área correspondiente a la solución óptima analítica. El gráfico muestra que la solución final asociada a cada proceso de optimización es significativamente distinta a la solución analítica para aquellos casos con valores altos del error prescrito γ , mientras que se aproxima hacia el valor exacto para los valores más bajos de γ .

 Figura 4. Influencia de γ sobre la evolución del error en la función objetivo con respecto a la solución analítica.

La figura 5 representa el efecto de γ sobre la forma de la solución final proporcionada por el algoritmo de optimización. Se puede observar claramente que para valores altos de γ el algoritmo converge hacia formas que son bastante distintas a la solución óptima analítica y solo se aproximan hacia esta última para valores decrecientes de γ . La diferencia en área entre la solución exacta y la obtenida para γ  = 1% es de solo un 0,38%. En cualquier caso, se debe tener en cuenta que la forma circular del contorno de la solución óptima analítica no se puede reproducir con total exactitud ya que los B -splines utilizados no son capaces de reproducir exactamente una forma circular.

 Figura 5. Influencia de γ sobre la solución óptima final obtenida (contorno en negro) y comparación con la solución óptima analítica (zona sombreada).

En el modelo geométrico se ha impuesto que, en sus extremos, la superficie externa del tubo sea perpendicular a las superficies de simetría. Cuando no se controla el error de discretización o para los valores más altos de γ , el tamaño de los elementos resulta demasiado grande para capturar dicha restrición geométrica, tal y como se muestra en el área sombreada de la figura 6 .

 Figura 6. Restricción geométrica y malla.

Tabla 2. Influencia de γ sobre la precisión de la función objetivo y el grado de satisfacción de las restricciones en tensión (tensiones de von Mises máximas de cada geometría final evaluadas con γ  ≤ 0,3%
γ A A  − Aopt σvm σvm  − σadm
Sin control 57,22 −18,00% 2222330 11,12%
20% 57,52 −17,58% 2229070 11,45%
10% 58,64 −15,97% 2265360 13,27%
5% 63,94 −8,39% 2074650 3,73%
2,5% 68,74 −1,50% 2025850 1,29%
1% 70,05 0,38% 2013710 0,69%

Este experimento numérico se ha reproducido utilizando programas de análisis MEF tanto comerciales como propios, modelos 2D y 3D de este experimento y algoritmos de optimización tanto genéticos como evolutivos, obteniéndose en todos los casos resultados muy similares.

## 3. Integración de los algoritmos evolutivos y los remallados h -adaptables mediante la proyección del estimador de error

Los resultados de la sección precedente han mostrado claramente que si no se asegura un nivel mínimo de precisión en los resultados utilizados para gobernar los procesos de optimización de formas estructurales, dichos procesos no convergerán hacia la verdadera solución óptima, produciendo soluciones que podrían violar notablemente las ecuaciones de restricción. No obstante, la utilización de técnicas h -adaptables para controlar la calidad de los análisis por elementos finitos implica un alto coste computacional adicional ya que cada configuración geométrica se debe evaluar más de una vez hasta alcanzar el nivel de precisión requerido. Debido al alto número de geometrías distintas que se deben analizar, el coste de la utilización de este tipo de técnicas puede ser crítico.

En el caso de la utilización de algoritmos evolutivos, y para reducir el coste computacional asociado con la generación de una malla adaptada para cada uno de los individuos (geometrías) a analizar, proponemos la utilización de la técnica descrita por Bugeda et al. [1] que a continuación se resume. El origen de esta técnica, y su aplicación a algoritmos de optimización deterministas, se puede ver en las referencias [2]  and [3] . En este trabajo hemos adaptado esta técnica para algoritmos evolutivos de optimización.

### 3.1. Algoritmo para la definición directa de mallas h -adaptadas para todos los individuos de una generación

Para cada generación de individuos, la estrategia propuesta se resume como sigue (ver su representación gráfica en la figura 7 ):

 Figura 7. Algoritmo para la definición directa de mallas h -adaptadas para todos los individuos de una generación.
• Selección de un individuo de referencia. Para cada generación se selecciona un individuo específico como referencia. Este individuo se puede fijar al inicio de todo el proceso de optimización como el diseño inicial. No obstante, los mejores resultados se han obtenido definiendo un individuo de referencia distinto para cada generación del proceso de optimización utilizando los valores medios de las variables de diseño de todos los individuos de la misma. Así, si se tiene una población con P individuos, los valores de vr que definen al de referencia se calculan como
 ${\displaystyle {\boldsymbol {\mbox{v}}}_{r}={\frac {1}{P}}\sum _{p=1}^{P}{\boldsymbol {\mbox{v}}}_{p}}$
( 7)

donde vp indica los valores de las variables de diseño v correspondientes al p -ésimo individuo utilizando la parametrización seleccionada

• Análisis de sensibilidad del individuo de referencia. Una vez que se ha obtenido una malla adaptada para el análisis del individuo de referencia, se realiza un análisis de sensibilidad completo de todas las magnitudes usadas en la estrategia de remallados adaptables. Esto incluye la evaluación de la sensibilidad del estimador de error, ver (5) , para cada elemento de la malla. La siguiente expresión fue obtenida en [2] para evaluar dicha magnitud:

 ${\displaystyle {\frac {\partial {\eta }_{e}^{2}}{\partial v_{i}}}={\int }_{{\Omega }_{\xi }}\left[{\left({\frac {\partial {\boldsymbol {\sigma }}^{\ast }}{\partial v_{i}}}-{\frac {\partial {\boldsymbol {\sigma }}^{h}}{\partial v_{i}}}\right)}^{T}{\boldsymbol {\mbox{D}}}^{-1}({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h})+({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h})^{T}{\frac {\partial {\boldsymbol {\mbox{D}}}^{-1}}{\partial v_{i}}}({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h})+({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h})^{T}{\boldsymbol {\mbox{D}}}^{-1}\left({\frac {\partial {\boldsymbol {\sigma }}^{\ast }}{\partial v_{i}}}-\right.\right.}$${\displaystyle \left.\left.{\frac {\partial {\boldsymbol {\sigma }}^{h}}{\partial v_{i}}}\right)+\right.}$${\displaystyle \left.({\boldsymbol {\sigma }}^{\ast }-{\boldsymbol {\sigma }}^{h})^{T}{\boldsymbol {\mbox{D}}}^{-1}({\boldsymbol {\sigma }}^{\ast }-\right.}$${\displaystyle \left.{\boldsymbol {\sigma }}^{h})tr\left({\boldsymbol {\mbox{J}}}^{-1}{\frac {\partial {\boldsymbol {\mbox{J}}}}{\partial v_{i}}}\right)\right]\vert {\boldsymbol {\mbox{J}}}\vert d{\Omega }_{\xi }}$
( 8)

donde Ωξ es el dominio del elemento en coordenadas locales y J es la matriz Jacobiana correspondiente a la transformación de coordenadas.

Fuenmayor et al. [4] extendieron el estimador de error de Zienkiewicz y Zhu [6] para el análisis de sensibilidad de formas para desarrollar un estimador del error de discretización en el análisis de sensibilidad, obteniendo también la expresión mostrada en (8) . Esto demuestra que la sensibilidad del error es equivalente al error en la sensibilidad.

La evaluación de ${\textstyle \partial {\eta }_{e}^{2}/\partial v_{i}}$ también requiere la evaluación del campo de tensiones mejorado σ* . En este trabajo se ha utilizado la técnica del alisado global por mínimos cuadrados [6] , aunque se podrían haber utilizado otros procedimientos como el SPR [7]  and [8] .

• Proyección de las magnitudes utilizadas en el proceso h-adaptable. Para cada uno de los individuos a analizar, los valores de todas las magnitudes involucradas en la estrategia de remallados adaptables, que se han obtenido en el paso anterior, se proyectan desde el individuo de referencia utilizando el correspondiente análisis de sensibilidad. Las magnitudes a proyectar son las coordenadas nodales (x e y ), el estimador de la norma energética del error η y la norma energética del error ||u ||. Para estas proyecciones se utilizan las siguientes expresiones:
 ${\displaystyle (x,y)_{p}=(x,y)_{r}+\sum _{i}^{n}(v_{p_{i}}-v_{r_{i}})\left({\frac {\partial x}{\partial v_{i}}},{\frac {\partial y}{\partial v_{i}}}\right)}$
( 9)
 ${\displaystyle {\eta }_{p}^{2}={\eta }_{r}^{2}+\sum _{i}^{n}(v_{p_{i}}-v_{r_{i}}){\frac {\partial {\eta }^{2}}{\partial v_{i}}}}$
( 10)
 ${\displaystyle \Vert {\boldsymbol {\mbox{u}}}{\Vert }_{p}^{2}=\Vert {\boldsymbol {\mbox{u}}}{\Vert }_{r}^{2}+}$${\displaystyle \sum _{i}^{n}(v_{p_{i}}-v_{r_{i}}){\frac {\partial \Vert {\boldsymbol {\mbox{u}}}{\Vert }^{2}}{v_{i}}}}$
( 11)

donde los subíndices r y p   están relacionados, respectivamente, con el individuo de referencia y el individuo hacia el cual se proyecta la información. ${\textstyle v_{i}}$ son las variables de diseño utilizadas para definir la geometría (ver [1] ).

Esta proyección proporciona, sin ningún cálculo previo adicional, una aproximación a los valores que se obtendrían para cada individuo específico si este se hubiese analizado utilizando la misma malla de elementos finitos que para el individuo de referencia. Se proporciona así, por tanto, la información necesaria para realizar un remallado adaptable sobre el nuevo diseño, incluso antes de efectuar ningún nuevo cálculo sobre el mismo.

La generación de cada nueva malla en el procedimiento de remallado requiere de la definición de un criterio de malla óptima. En este trabajo se considera que una malla es óptima cuando la densidad del error está uniformemente distribuida sobre todo el dominio (ver referencia [3] ).

La figura 8 muestra un ejemplo de análisis con una generación de 30 individuos. En la figura se muestra la malla h -adaptada del individuo de referencia junto con las obtenidas para cada uno de los individuos de la generación. Obsérvese que tan solo uno de los individuos ha requerido un segundo paso de refinamiento.

 Figura 8. Ejemplo de mallados de una generación de individuos obtenidos por proyección desde el individuo de referencia.

## 4. Ejemplos numéricos

Esta sección muestra cuatro ejemplos de optimización resueltos con la estrategia propuesta.

### 4.2. Optimización de la forma de un volante de inercia

Consideremos el problema de diseño óptimo mostrado en la figura 9 a) en el que se define un volante de inercia cuya zona interior es susceptible de recibir mejoras geométricas. Se ha aplicado una fuerza centrífuga correspondiente a una velocidad de rotación de 3 rad/s, y una fuerza tangencial distribuida sobre el contorno exterior de 1 N/mm2 . Las propiedades del material son E = 210 GPa, coeficiente de Poisson μ  = 0,29 y densidad ρ  = 7.860 kg/m3 .

 Figura 9. Optimización de un volante de inercia. Optimización topológica.

Este problema ya fue resuelto con anterioridad, ver [16] , utilizando una combinación de dos herramientas de optimización. En una primera etapa, se utilizó la malla representada en la figura 9 b) (6.000 cuadriláteros bilineales en tensión plana) para proceder a una optimización topológica, obteniendo la topología representada en la figura 9 c). A continuación, para obtener una forma más refinada, en una segunda etapa se utilizó una estrategia de optimización de forma basada en un algoritmo determinista con análisis de sensibilidad. La geometría inicial para esta segunda etapa se muestra en la figura 10 a). Este proceso redujo el peso total desde 1,53 kg hasta 1,45 kg.

 Figura 10. Optimización de un volante de inercia. Parametrización.

En este trabajo se presenta la resolución de este problema sustituyendo el algoritmo de optimización determinista por el algoritmo evolutivo con la integración de los remallados h -adaptables presentada.

La geometría de la solución obtenida mediante la optimización topológica presenta 4 ejes de simetría. No obstante, la carga tangencial aplicada no permite utilizar dichos ejes para simplificar el problema y, por ello, no es posible limitar la resolución del mismo a solo un octavo del total. De todos modos, teniendo en cuenta la presencia de los ejes de simetría se ha seguido el siguiente procedimiento:

• La geometría del problema se ha parametrizado utilizando un total de 60 variables de diseño que son algunas de las coordenadas de los puntos de definición cuyo movimiento está permitido.
• Se han aplicado restricciones lineales sobre dichas variables de diseño para conseguir que la solución del problema mantenga los ejes de simetría. De esta manera se han conseguido un total de solo 8 variables de diseño independientes (fig. 11 b).

 Figura 11. Optimización de un volante de inercia. Evolución de la función objetivo (peso) y de la diferencia con respecto a la solución final (error%).
• La función objetivo adoptada es el peso total del volante de inercia.
• El máximo valor permitido para la tensión de von Mises es de 100 N/mm2 .
• Para el algoritmo evolutivo de optimización se han utilizado 300 generaciones con 15 individuos cada una.
• El máximo error relativo admisible en la norma energética se ha fijado en un γ  = 5%.

La figura 11 muestra los valores del peso obtenidos a lo largo del proceso evolutivo de optimización, y la diferencia porcentual (error%) entre la función objetivo y el valor final obtenido. La figura 12 muestra la evolución de la geometría y las correspondientes mallas adaptadas. El peso total se reduce desde 1,53 kg hasta 1,4445 kg, lo cual es un resultado muy cercano al obtenido mediante el algoritmo determinista (1,45 kg). No obstante, si se tiene en cuenta solo la zona del dominio susceptible de ser modificada la reducción es de 0,25 a 0,17 kg.

 Figura 12. Optimización de un volante de inercia. Evolución de la geometría. Mallas adaptadas.

La figura 13 muestra una comparación entre la geometría original (línea discontinua) y la solución final obtenida (zona sombreada).

 Figura 13. Optimización de un volante de inercia. Comparación entre las geometrías original (trazo discontinuo) y optimizada (area sombreada).

### 4.3. Optimización de la forma de un gancho

 Figura 14. Optimización de un gancho. Geometría inicial y variables de diseño. Geometría optimizada.

Durante el proceso de optimización se han utilizado 400 generaciones con 20 individuos cada una. El máximo error relativo admisible en norma energética se ha fijado en un γ  = 2%.

La figura 14 permite comparar la solución final obtenida, con un área de 88,71 cm2 , con la geometría inicial que tenía un área de 181,22 cm2 . Puede observarse cómo el proceso de optimización desplaza el eje vertical del gancho para que su directriz coincida con la de la fuerza resultante de la carga aplicada, eliminando así la presencia de momento flector sobre esa región.

La figura 3 muestra que los resultados obtenidos durante las primeras etapas del proceso de optimización no están afectados de forma significativa por el nivel de precisión prescrito para el análisis. Esto ha permitido implementar una técnica muy simple, que también contribuye a reducir el coste computacional de la optimización, que ha sido utilizada en la resolución de este problema. Esta técnica consiste simplemente en definir el valor admisible para la estimación de la norma energética del error γ como una función del número de generación, haciendo decrecer su valor a medida que el proceso de optimización avanza. Esta técnica reduce considerablemente el coste computacional asociado a los análisis de las primeras generaciones de individuos sin afectar a la precisión de los resultados obtenidos en las últimas etapas del proceso. La evolución del error prescrito a lo largo de la optimización se puede ver en la figura 15 . Entre las generaciones 1 y 150 se ha definido un decrecimiento exponencial de γ desde un 20% hasta un 2%, habiéndose impuesto el límite más exigente de γ  = 2% solo a partir de la generación 150.

 Figura 15. Optimización de un gancho. Evolución del nivel de precisión prescrito (γ ) a medida que progresa el proceso de optimización.

La figura 16 muestra ejemplos de las distintas densidades de malla utilizadas durante el proceso de optimización. Obsérvese que dicho proceso es capaz de proporcionar aproximaciones razonables a la geometría final a partir de las primeras etapas del proceso utilizando mallas groseras (ver las geometrías para las generaciones 6 y 57) que requieren un bajo coste computacional. La figura 17 muestra la evolución del mejor individuo obtenido para una serie de generaciones a lo largo del proceso.

 Figura 16. Optimización de un gancho. Densidades de malla utilizadas durante el proceso de optimización. Las geometrías representadas corresponden a un individuo representativo de las generaciones 0, 3, 6, 57 y 375.

 Figura 17. Optimización de un gancho. Evolución del mejor individuo durante el proceso de optimización. Debajo de cada geometría se indica la generación en la que la misma se ha generado.

### 4.4. Optimización de la forma de la sección transversal de un túnel doble sumergido

 Figura 18. Optimización de un túnel sometido a presión. Modelo de análisis (en la versión electrónica de este artículo se puede consultar esta figura a color).

El máximo error relativo admisible en la norma energética para el análisis de la solución final se ha fijado en un γ  = 2%.

Al igual que en el ejemplo anterior, a lo largo del proceso se ha hecho variar en valor prescrito de γ con el número de generación. Se ha especificado una variación exponencial para el nivel de error admisible desde γ ≤ 30% para la primera generación hasta γ ≤ 2% para la generación 100, tal como se muestra en la figura 19 . A efectos de comparación, este problema se ha resuelto una segunda vez imponiendo un nivel de error constante γ ≤ 2% para todo el proceso.

 Figura 19. Optimización de un túnel sometido a presión. Evolución del nivel de precisión prescrito (γ ) a medida que progresa el proceso de optimización (en la versión electrónica de este artículo se puede consultar esta figura a color).

Hay que tener en cuenta que la variació de γ con el número de generación no ha de ser necesariamente una función exponencial. Los mayores beneficios del uso de esta ténica se obtendrían si dicha variación de γ está adecuadamente acoplada a la evolución del proceso de optimización. Cómo obtener este adecuado acoplamiento será materia de nuestras investigaciones en el futuro.

La figura 20 muestra la evolución de la función objetivo (área) durante ambos procesos de optimización. Esta figura muestra claramente que la utilización de valores altos de γ para las primeras generaciones y su progresiva reducción proporciona una evolución similar de los mejores individuos y reducciones del área similares después de 400 generaciones (de 4.105 m2 a 2.824 m2 para γ constante y a 2.772 m2 para γ variable). Las diferencias en los valores finales de la figura 20 son razonables teniendo en cuenta la naturaleza del algoritmo evolutivo y el hecho de que el proceso de optimización no se pueda considerar que haya convergido después de 400 iteraciones.

 Figura 20. Optimización de un túnel sometido a presión. Evolución del área (en la versión electrónica de este artículo se puede consultar esta figura a color).

La figura 21 muestra las diferencias entre la geometría original y las optimizadas en ambos casos. Se observa que la técnica de la reducción gradual de γ permite reducir sensiblemente el coste computacional sin perjudicar la evolución del proceso de optimización. Para este caso, con 400 iteraciones con error de discretización constante se ha necesitado un tiempo un 13% superior.

 Figura 21. Optimización de un túnel sometido a presión. Comparación entre las geometrías original (negro) y optimizada, con error (γ ) variable (línea a trazos), con error (γ ) constante (línea punteada) (en la versión electrónica de este artículo se puede consultar esta figura a color).

Dado que tanto las condiciones de contorno como la geometría son simétricas, la solución óptima también tendrá que ser simétrica. La figura 22 muestra en trazos continuos la geometría final de ambos túneles y en trazo discontinuo, la imagen especular del túnel de la derecha representada sobre el de la izquierda. Una solución inicial considerablemente asimétrica como la utilizada podría hacer que un algoritmo de optimización basado en gradiente se quedase estancado en un mínimo local proporcionando una solución no simétrica. Este tipo de comportamientos representa uno de los principales problemas de estos algoritmos. Frente a este problema, una de las principales ventajas de los algoritmos evolutivos es que son capaces de proporcionar el mínimo global del problema de optimización. Así, la figura 22 muestra que sin aportar ninguna información relativa a la simetría de la solución final, el algoritmo evolutivo ha sido capaz de evaluar una solución que tiende a reponer la simetría geométrica característica del mínimo global.

 Figura 22. Optimización de un túnel sometido a presión. Comparación entre los dos agujeros. Agujero izquierda linea negra. Agujero derecho línea continua. Imagen especular del agujero derecho sobre el izquierdo en la línea a trazos (en la versión electrónica de este artículo se puede consultar esta figura a color).

## 5. Conclusiones

Los resultados preliminares han demostrado que el nivel de error de discretización máximo impuesto a los análisis por el MEF tiene poca incidencia durante las primeras etapas del proceso de optimización. Por ello, se ha propuesto una técnica de reducción paulatina del mismo. A pesar de su simplicidad, esta técnica permite reducir de forma considerable el coste computacional asociado los análisis de las primeras generaciones de individuos sin afectar a la precisión de los resultados obtenidos al final de la optimización.

Los autores primero, tercero y quinto agradecen el soporte económico recibido para el desarrollo de este artículo a través de los proyectos de investigación DPI2007-66773-C02-01 y DPI2010-20542 del Ministerio de Educación y Ciencia (España) y los fondos recibidos de la Generalitat Valenciana y de la Universitat Politècnica de València.

El segundo autor agradece el soporte económico recibido para el desarrollo de este artículo a través del proyecto de investigación DPI2008-05250 del Ministerio de Educación y Ciencia (España).

Los autores también agradecen a MsC Santos López Real su ayuda con la utilización de los algoritmos evolutivos.

## References

1. [1] G. Bugeda, J.J. Ródenas, E. Oñate; An integration of a low cost adaptive remeshing strategy in the solution of structural shape optimization problems; Computers & Structures, 86 (2008), pp. 1563–1578
2. [2] G. Bugeda, J. Oliver; A general methodology for structural shape optimization problems using automatic adaptive remeshing; International Journal for Numerical Methods in Engineering, 18 (1993), pp. 3161–3185
3. [3] G. Bugeda, E. Oñate; A methodology for adaptive mesh refinement in optimum shape design problems; Computing Systems in Engineering, 5 (1994), pp. 91–102
4. [4] F.J. Fuenmayor, J.L. Oliver, J.J. Ródenas; Extension of the Zienkiewicz-Zhu error estimator to shape sensitivity analysis; International Journal for Numerical Methods in Engineering, 40 (1997), pp. 1413–1433
5. [5] R. Storn, K. Price, Differential Evolution - A simple and efficient adaptive scheme for global optimization over continuous space, International Computer Science Institute, Berkeley, CA, USA, Technical Report, March 1995 [consultado 10 Ene 2011], ftp://ftp.icsi.berkeley.edu/pub/techreports/1995/tr-95-012.pdf .
6. [6] O.C. Zienkiewicz, J.Z. Zhu; A simple error estimator and adaptive procedure for practical engineering analysis; International Journal for Numerical Methods in Engineering, 24 (1987), pp. 337–357
7. [7] O.C. Zienkiewicz, J.Z. Zhu; The superconvergent patch recovery and a posteriori error estimates. Part I: The recovery technique; International Journal for Numerical Methods in Engineering, 33 (1992), pp. 1331–1364
8. [8] O.C. Zienkiewicz, J.Z. Zhu; The superconvergent patch recovery and a posteriori error estimates. Part II: Error Estimates and Adaptivity; International Journal for Numerical Methods in Engineering, 33 (1992), pp. 1365–1382
9. [9] O.C. Zienkiewicz, J.Z. Zhu, J. Wu; Superconvergent Patch Recovery Techniques - some further Tests; Communications in Numerical methods in Engineering, 9 (1993), pp. 251–258
10. [10] N.E. Wiberg, F. Abdulwahab; Patch Recovery Based on Superconvergent Derivatives and Equilibrium; International Journal for Numerical Methods in Engineering, 36 (1993), pp. 2703–2724
11. [11] N.E. Wiberg, F. Abdulwahab, S. Ziukas; Enhanced superconvergent patch recovery incorporating equilibrium and boundary conditions; International Journal for Numerical Methods in Engineering, 37 (1994), pp. 3417–3440
12. [12] J.J. Ródenas, M. Tur, F.J. Fuenmayor, A. Vercher; Improvement of the Superconvergent Patch Recovery technique by the use of constraint equations: the SPR-C technique; International Journal for Numerical Methods in Engineering, 70 (2006), pp. 705–727
13. [13] P. Díez, J.J. Ródenas, O.Z. Zienkiewicz; Equilibrated Patch Recovery error estimates: simple and accurate upper bounds of the error; International Journal for Numerical Methods in Engineering, 69 (2007), pp. 2075–2098
14. [14] I.D. Faux, M.J. Pratt; Computational Geometry for Design and Manufacture; Ellis Horwood Limited (1987)
15. [15] J.J. Ródenas, F.J. Fuenmayor, J.E. Tarancón; A numerical methodology to asses the quality of the design velocity field computation methods in shape sensitivity analysis; International Journal for Numerical Methods in Engineering, 59 (2004), pp. 1725–1747
16. [16] J. Sienz, E. Hinton, G. Bugeda, S. Bulman; Some studies on Integrating Topology and shape optimization; M. Papadrakakis, B.H.V. Topping (Eds.), Innovative Computational Methods for Structural Mechanics, Saxe-coburg Publications (1999), pp. 223–256

### Document information

Published on 01/03/12
Accepted on 26/05/11
Submitted on 27/01/11

Volume 28, Issue 1, 2012
DOI: 10.1016/j.rimni.2011.11.005
Licence: Other

### Document Score

0

Times cited: 2
Views 32
Recommendations 0