A minimum weight formulation with stress constraints in topology optimization of structures

==Resumen==

El diseño óptimo de estructuras ha estado tradicionalmente orientado a la resolución de problemas de optimización de formas y dimensiones. Sin embargo, más recientemente ha surgido otra rama de investigación que propone modelos que proporcionan soluciones estructurales óptimas y que no requieren la definición previa de la tipología estructural: la optimización topológica de estructuras. Estas formulaciones proporcionan tanto la tipología estructural como la forma y dimensiones óptimas. Las formulaciones más habituales de estos planteamientos pretenden obtener una solución que maximice la rigidez de la estructura dadas unas limitaciones en la cantidad de material a utilizar. Estas formulaciones han sido ampliamente analizadas y utilizadas en la práctica pero presentan inconvenientes muy importantes tanto desde un punto de vista numérico como práctico. En este artículo se propone una formulación diferente a la de máxima rigidez para el problema de optimización topológica de estructuras que minimiza el peso e incorpora restricciones en tensión. Las ventajas de este tipo de planteamientos son muy importantes dado que se minimiza el coste de la solución, que es la situación más habitual en ingeniería, y además se garantiza la validez estructural de la misma. Además esta formulación permite evitar algunos de los problemas e inestabilidades numéricas que presentan las formulaciones de máxima rigidez. Finalmente, se resuelven algunos ejemplos prácticos para comprobar la validez de los resultados y las ventajas que ofrece la formulación propuesta.

'''Palabras clave''': Optimización topológica estructural, mínimo peso, restricciones en tensión, restricciones por bloques

==Abstract==

Optimum design of structures has been traditionally focused on the analysis of shape and dimensions optimization problems. However, more recently a new discipline has emerged: the topology optimization of the structures. This discipline states innovative models that allow to obtain optimal solutions without a previous definition of the type of structure being considered. These formulations obtain the optimal topology and the optimal shape and size of the resulting elements. The most usual formulations of the topology optimization problem try to obtain the structure of maximum stiffness. These approaches maximize the stiffness for a given amount of material to be used. These formulations have been widely analyzed and applied in engineering but they present considerable drawbacks from a numerical and from a practical point of view. In this paper the author propose a different formulation, as an alternative to maximum stiffness approaches, that minimizes the weight and includes stress constraints. The advantages of this kind of formulations are crucial since the cost of the structure is minimized, which is the most frequent objective in engineering, and they guarantee the structural feasibility since stresses are constrained. In addition, this approach allows to avoid some of the drawbacks and numerical instabilities related to maximum stiffness approaches. Finally, some practical examples have been solved in order to verify the validity of the results obtained and the advantages of the proposed formulation.

'''Keywords''': Structural topology optimization, minimum weight, stress constraints, block aggregated constraints


==1. Introducción==

El diseño óptimo de estructuras es una disciplina que nace como la prolongación natural de los métodos de análisis y cálculo estructural. Así, los planteamientos de diseño óptimo surgen de la necesidad de proponer técnicas objetivas que aborden la fase de concepción y definición de la solución estructural y que por tanto no se limiten a la realización de cálculos y comprobaciones de la capacidad resistente.

Los primeros modelos de diseño óptimo de estructuras aparecen a partir de los años 40 del siglo xx motivados por el gran desarrollo de la industria aeronáutica si bien ya Michell en 1904 [[#bib0005|[1]]]  propuso técnicas analíticas para obtener las soluciones estructurales óptimas para algunos problemas específicos. Sin embargo, esta disciplina se estudió de forma más rigurosa en los años 60 con los trabajos de Schmit [[#bib0010|[2]]], que sienta las bases de su análisis. A partir de estos primeros trabajos se han propuesto numerosos modelos para el diseño óptimo de estructuras hasta la actualidad. Así, aparecieron los modelos de optimización de dimensiones y, posteriormente, los modelos de optimización de formas, ampliamente desarrollados y estudiados en la actualidad.      

Más recientemente, se establece una nueva rama de estudio en esta disciplina: la optimización topológica de estructuras, a partir de los trabajos realizados por Bendsøe y Kikuchi en 1988 [[#bib0015|[3]]]. En esta disciplina se pretende obtener la distribución de material más adecuada que hay que disponer sobre un dominio predefinido para soportar las cargas que se imponen. Así, estos modelos pretenden obtener la solución estructural más adecuada sin necesidad de definir e imponer ni la tipología estructural ni la configuración de los elementos resistentes (geometría, dimensiones, sección,…). Se pretenden generalizar, por tanto, los planteamientos teóricos propuestos por Michell [[#bib0005|[1]]]  en 1904 y desarrollar modelos numéricos que permitan obtener soluciones óptimas a un mayor número de problemas en la práctica.      

La optimización topológica de estructuras pretende obtener la distribución de una cierta cantidad de material sobre un dominio predefinido que proporciona la estructura resistente más adecuada. De acuerdo con los planteamientos originales del problema de optimización topológica propuestos por Bendsøe [[#bib0015|[3-7]]], el objetivo principal de estos modelos consiste en obtener distribuciones de material adecuadas indicando las partes del dominio en las que debe o no debe existir material. En consecuencia, el planteamiento original del problema da lugar a modelos de diseño óptimo con variables de diseño binarias (material o vacío). En la práctica, los modelos más habituales que se proponen para abordar este problema plantean formulaciones en las que se asume que las variables de diseño son continuas para evitar el tratamiento de problemas de optimización con variables discretas. Asimismo, se asume que el estado material en cada elemento es uniforme en todos sus puntos. Además, con estas condiciones se evita que el problema esté mal puesto desde un punto de vista teórico, tal y como ocurre si se analiza el estado material en cada punto del dominio. Por lo tanto, en la práctica se establece una variable de diseño por cada elemento: la densidad relativa, que define el estado material del mismo. La densidad relativa de cada elemento puede tomar valores entre 0, que indica elemento vacío, y 1, que indica elemento con material sólido.      

Los modelos que utilizan variables de diseño continuas presentan dos inconvenientes fundamentales que es necesario tener en cuenta. El primero de ellos consiste en la necesidad de establecer un modelo de análisis que permita obtener el comportamiento estructural para valores intermedios de las densidades relativas de los elementos. Esta exigencia se ha resuelto en la práctica mediante la definición de modelos de microestructura resistente y la correspondiente homogeneización de las propiedades [[#bib0020|[4,5,7-13]]] para su empleo en modelos de cálculo de elementos finitos. El segundo inconveniente consiste en la necesidad de proporcionar distribuciones óptimas de material esencialmente binarias a partir de un modelo con variables de diseño continuas. Este punto se ha tratado habitualmente mediante el uso de funciones de penalización de las densidades relativas intermedias, o bien a través del modelo de microestructura propuesto o bien mediante funciones objetivo convenientemente definidas para esta finalidad.      

Las primeras formulaciones del problema de optimización topológica de estructuras plantean la obtención de la distribución de una cierta cantidad de material que maximiza la rigidez de la estructura resultante para el conjunto de cargas aplicadas. Estas primeras formulaciones del problema no corresponden con el tipo de planteamientos que tradicionalmente se analizan en diseño óptimo de estructuras ya que, desde un punto de vista ingenieril y económico, se pretende minimizar el coste o el peso de la estructura garantizando la capacidad estructural para soportar las cargas aplicadas. Además, las formulaciones de máxima rigidez presentan inestabilidades numéricas muy importantes como por ejemplo la disposición de material en damero o checkerboard [[#bib0025|[5,14]]] . Por otra parte, estas formulaciones de máxima rigidez ofrecen ventajas computacionales muy importantes que facilitan su uso en la práctica.      

Sin embargo, más recientemente se han propuesto formulaciones alternativas a esta de máxima rigidez que abordan los planteamientos más habituales en ingeniería (mínimo peso con restricciones en rigidez y mínimo peso con restricciones en tensión y/o en desplazamientos entre otros) [[#bib0075|[15-30]]]. En este artículo se propone y analiza una formulación que pretende obtener la solución estructural de mínimo peso teniendo en cuenta restricciones en tensión. Esta formulación presenta numerosas ventajas frente a las formulaciones de máxima rigidez dado que se evitan los inconvenientes más relevantes de estas formulaciones. Así, se evitan disposiciones de material en damero al imponer las restricciones en tensión y se aborda el problema tradicionalmente analizado en diseño óptimo de estructuras. Además, el uso de restricciones en tensión/desplazamientos garantiza la validez estructural de la solución óptima para las cargas y acciones externas que se aplican.      

El modelo de diseño que se propone en este trabajo utiliza como variables de diseño: las densidades relativas de cada elemento de la malla de elementos finitos con la que se resuelve el cálculo estructural. La función objetivo utilizada está basada en el peso de la estructura e incorpora una penalización sobre los valores intermedios de las densidades relativas para favorecer la obtención de soluciones óptimas con distribuciones de material esencialmente binarias. Las restricciones que se imponen son condiciones en tensión que comparan y limitan el valor de una cierta tensión de referencia, en la práctica la tensión de Von Mises, a la máxima tensión admisible de fallo del material. Para imponer estas restricciones en tensión se plantean tres formulaciones diferentes que presentan ventajas tanto desde un punto de vista computacional como desde un punto de vista teórico y numérico.

El artículo está estructurado en [[#sec0135|8]]  apartados. En este primer apartado de introducción se presenta el problema de optimización topológica, que se desarrolla de forma más específica en los siguientes apartados. Así, en el apartado [[#sec0010|2]]  se presenta el modelo de diseño óptimo haciendo especial hincapié en el planteamiento de la función objetivo y de las restricciones del problema. El apartado [[#sec0040|3]]  está destinado a abordar y analizar el modelo de cálculo de estructuras utilizado para obtener el comportamiento estructural de los diseños propuestos. Una vez planteado el problema de diseño óptimo y el modelo de cálculo, en el apartado [[#sec0065|4]]  se presenta el algoritmo de optimización utilizado para obtener la solución óptima. Los algoritmos de optimización propuestos requieren el cálculo de derivadas de primer y segundo orden de la función objetivo y las restricciones. Los algoritmos desarrollados para calcular estas derivadas se analizan en el apartado [[#sec0095|5]]. El apartado [[#sec0110|6]]  se dedica a establecer las técnicas de cálculo en paralelo que se pueden utilizar para reducir el tiempo de cálculo necesario para resolver el problema de optimización propuesto. Finalmente, en los apartados [[#sec0115|7]]  y [[#sec0135|8]]  se presentan algunos ejemplos de aplicación de los problemas de optimización planteados y se indican las conclusiones más importantes sobre este trabajo, respectivamente.      

<span id='sec0010'></span>

==2. Formulación del problema de optimización topológica de estructuras==

Tal y como se ha comentado en el apartado de introducción, el modelo que se propone en este trabajo para el problema de optimización topológica de estructuras se basa en un planteamiento que minimiza el peso y está sujeto a restricciones en tensión. Las variables de diseño de este problema indican la densidad relativa de cada uno de los elementos de la malla, en los cuales se asume que esta propiedad es uniforme.

Así, se puede plantear el problema de optimización topológica de estructuras de forma general como:

<span id='eq0005'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| <math>\hbox{hallar } \quad \boldsymbol{\rho }=\lbrace {\rho }_i\rbrace ,i=</math><math>1,\ldots ,n</math>
|-
|<math>\hbox{que minimice } \quad F(\boldsymbol{\rho })</math>
|-
|<math>\hbox{ verificando } \quad g_j\left(\boldsymbol{r}_j^o,\boldsymbol{\rho }\right)\leq 0,j=</math><math>1,\ldots ,m,</math>
|-
|<math>{\rho }_{m\acute{\imath}n}\leq {\rho }_i\leq {\rho }_{m\acute{a}x},i=</math><math>1,\ldots ,n.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

donde <math display="inline">\boldsymbol{\rho }=\left\{{\rho }_i\right\}</math>  es el vector de variables de diseño, ''F'' ('''''ρ''''') es la función objetivo a minimizar, ''g''<sub>''j''</sub>  son las restricciones en desigualdad del problema, ''n''  es el número de variables de diseño del problema y ''m''  es el número de restricciones de desigualdad. Los límites laterales de las variables de diseño se definen como ''ρ''<sub>''mín''</sub>  y ''ρ''<sub>''máx''</sub>  para el valor inferior y superior, respectivamente. Para el problema que se plantea las restricciones laterales toman el valor ''ρ''<sub>''máx''</sub>  = 1 (material sólido) para el límite superior y un valor inferior ''ρ''<sub>''mín''</sub>  ligeramente superior a 0 para evitar que la matriz de rigidez de la estructura sea singular. El valor habitualmente utilizado en la práctica es ''ρ''<sub>''mín''</sub>  = 0.001  [[#bib0015|[3,5]]].      

Los aspectos fundamentales del planteamiento del problema de optimización propuesto en [[#eq0005|(1)]]  residen en la definición de la función objetivo a minimizar y de las restricciones que se imponen. De estos aspectos depende el tipo de problema de optimización que se obtiene así como las técnicas que es necesario utilizar para resolverlo. Por lo tanto, procedemos a desarrollar y analizar la función objetivo propuesta en [[#eq0005|(1)]] y, posteriormente, al estudio de las restricciones en tensión.      

===2.1. Planteamiento de la función objetivo===

La función objetivo de la formulación que se propone en este trabajo se basa en la minimización del peso de la estructura, es decir:

<span id='eq0010'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F=\sum_{e=1}^{N_e}{\int }_{{\Omega }_e}{\rho }_e\quad {\gamma }_{mat}\quad d{\Omega }_e</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

siendo ''N''<sub>''e''</sub>  el número de elementos de la malla de elementos finitos, Ω<sub>''e''</sub>  el dominio del elemento ''e'' , ''ρ''<sub>''e''</sub>  la densidad relativa asociada al elemento ''e''  y ''γ''<sub>''mat''</sub>  el peso específico del material utilizado.      

Sin embargo, como ya se ha comentado en el apartado de introducción, es muy habitual incluir factores de penalización en la formulación para evitar disposiciones de material con un número elevado de zonas con densidades relativas intermedias. En este trabajo se incorpora una penalización de las densidades intermedias en la función objetivo dado que se considera que la obtención de distribuciones esencialmente binarias debe formar parte de la función objetivo del problema. Este planteamiento no es el habitual en otros trabajos sobre optimización topológica en los que la penalización de las densidades intermedias se incorpora en el modelo de microestructura que se utiliza. Estos planteamientos se explicarán en detalle más adelante en este artículo cuando se presente el modelo cálculo de estructuras por elementos finitos. Por lo tanto, para incorporar esta penalización de las densidades intermedias se plantea una función objetivo modificada, a partir de la propuesta en [[#eq0010|(2)]], como:

<span id='eq0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_p=\sum_{e=1}^{N_e}{\int }_{{\Omega }_e}({\rho }_e)^{1/p}\quad {\gamma }_{mat}\quad d{\Omega }_e</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

donde ''p''  es el factor de penalización de densidades intermedias. Este factor introduce una penalización en la función para las densidades intermedias cuando ''p''  > 1. Si, por el contrario, ''p''  < 1 se favorece la aparición de material con densidades intermedias. En el caso de que ''p''  = 1 se obtendrá la función objetivo inicial propuesta [[#eq0010|(2)]]  sin penalización alguna  [[#bib0120|[24,31-34]]].      

El efecto de este parámetro de penalización ''p''  sobre la función objetivo puede observarse fácilmente en la  [[#fig0005|figura 1]].

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr1.jpg|center|371px|Exponente de la función objetivo para distintos valores de penalización.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 1'''. Exponente de la función objetivo para distintos valores de penalización.

</span>
|}

Además, si la solución óptima presenta una distribución de material fundamentalmente binaria (0-1) la función de mínimo coste coincide con la función de mínimo peso dado que la penalización se establece para las zonas con densidades relativas intermedias. Por lo tanto, la formulación de mínimo coste representa, si la solución final es binaria, el peso de la estructura y además impone penalización para los elementos con densidad relativa intermedia.

Este planteamiento de la función objetivo es el que se utiliza habitualmente en la práctica para el problema que se plantea en este artículo. Sin embargo, en algunas aplicaciones es posible que no se alcancen soluciones esencialmente binarias debido a las restricciones en tensión que se imponen en el modelo. Este fenómeno es fácilmente explicable si se analiza la función objetivo conjuntamente con las restricciones en tensión. Así, las derivadas de la función objetivo ''F''<sub>''p''</sub>  propuesta en [[#eq0015|(3)]]  con respecto a cada una de las variables de diseño son siempre positivas, lo que indica que una reducción de la función objetivo requiere necesariamente que se reduzca el valor de la variable de diseño con respecto a la que se deriva. Si además se tiene en cuenta que las restricciones en tensión limitan inferiormente los valores de las densidades relativas para evitar que se incumplan las restricciones, se puede observar fácilmente que las soluciones binarias pueden no alcanzarse en la práctica. Esta afirmación se puede observar fácilmente en la [[#fig0010|figura 2]] , donde se analiza el efecto sobre una única variable de diseño.

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr2.jpg|center|365px|Ejemplo de solución óptima para distintos valores de penalización p cuando se ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 2'''. Ejemplo de solución óptima para distintos valores de penalización ''p''  cuando se consideran restricciones en tensión.                  

</span>
|}

En consecuencia, en algunas aplicaciones prácticas es posible que no se obtengan soluciones con la función objetivo propuesta en [[#eq0015|(3)]]. En estos casos, se puede forzar una solución binaria usando una función objetivo distinta destinada a este fin. Para ello se propone la función objetivo:

<span id='eq0020'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F_b(\boldsymbol{\rho })=\sum_{e=1}^{N_e}{\int }_{{\Omega }_e}{\Psi }_b({\rho }_e)\quad {\gamma }_{mat}d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Psi }_b({\rho }_e)={\rho }_e^{\frac{1}{(1+\beta )^2}}+</math><math>\beta {exp}^{-\beta {\rho }_e}sin(\pi {\rho }_e),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

siendo ''β''  un parámetro que toma valores en el intervalo [0,260, 7,404]. Valores de ''β  ''  que no pertenezcan al intervalo no conducen necesariamente a soluciones binarias dado que <math display="inline">{{\Psi }_b^{{'}}({\rho }_e)}_{{\rho }_e=1}>0\quad \forall e=</math><math>1,\ldots N_e</math>.      

En la [[#fig0015|figura 3]]  se muestran los valores del integrando Ψ<sub>''b''</sub> (''ρ''<sub>''e''</sub>) de la función objetivo propuesta en [[#eq0020|(4)]]  para distintos valores del parámetro ''β''  y de la densidad relativa ''ρ''<sub>''e''</sub>.

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr3.jpg|center|362px|Integrando de la función objetivo (Ψb(ρe)) (4) para distintos valores de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 3'''. Integrando de la función objetivo (Ψ<sub>''b''</sub> (''ρ''<sub>''e''</sub>)) [[#eq0020|(4)]]  para distintos valores de penalización.                  

</span>
|}

Esta modificación introduce cambios muy agresivos en los planteamientos de mínimo peso ya que las funciones resultantes presentan derivadas con respecto a las variables de diseño negativas para ''ρ''<sub>''e''</sub>  = 1. Por este motivo, las formulaciones resultantes no se comportan como formulaciones de mínimo peso sino que su objetivo es proporcionar soluciones binarias teniendo en cuenta el peso de la estructura. Es decir, las soluciones obtenidas no corresponderán a soluciones de mínimo peso (desde un punto de vista continuo de las variables de diseño) pero presentarán distribuciones de material binarias, lo que facilita en gran medida la utilización en la práctica de los diseños logrados. Para ello será necesario incrementar los valores de ''β''  de forma progresiva a lo largo del proceso de optimización hasta el límite superior del intervalo recomendado.      

Por el contrario, estos planteamientos también introducen mínimos locales y otros efectos no deseados. Por este motivo, la utilización de esta función objetivo se plantea como una segunda fase en el proceso de optimización que se utilizará solo en aquellos casos en los que las soluciones obtenidas con la función objetivo propuesta en [[#eq0015|(3)]]  no sean satisfactorias. Así, una vez alcanzada la solución óptima de mínimo peso utilizando la función objetivo propuesta en [[#eq0015|(3)]]  se puede plantear una segunda fase del proceso en la que se busquen las soluciones binarias utilizando la función objetivo modificada propuesta en [[#eq0020|(4)]]. Algunos ejemplos de aplicación práctica de este tipo de formulaciones pueden encontrarse en [[#bib0175|[35]]].      

Una vez planteadas y analizadas las diferentes formulaciones de la función objetivo que se han planteado para el problema de optimización topológica de estructuras propuesto se analizan las distintas alternativas para incorporar restricciones en tensión.

===2.2. Planteamiento de las restricciones en tensión===

La aplicación de restricciones en tensión es el planteamiento más habitual en otras disciplinas del diseño óptimo de estructuras como la optimización de formas y la optimización de dimensiones. Sin embargo, en optimización topológica de estructuras el objetivo más habitual es, como ya se ha comentado, la maximización de la rigidez fijado un volumen de material a emplear [[#bib0040|[3,5,8,10]]]. También hay que destacar algunas contribuciones que minimizan el peso y establecen como restricción un valor mínimo de rigidez para la estructura (por ejemplo, Liang ''et al''.  [[#bib0180|[36]]]). En otras aplicaciones la función objetivo consiste en minimizar las tensiones dada una cantidad de material a emplear (Allaire ''et al''.  [[#bib0075|[15]]]).      

En este trabajo se pretende abordar esta problemática desde el mismo punto de vista que en el resto de disciplinas de la optimización de estructuras minimizando el coste y estableciendo restricciones de tipo tensional. Con esta finalidad se han desarrollado tres formulaciones diferentes de las restricciones: una formulación local, una formulación global y una formulación agregada por bloques.

<span id='sec0025'></span>
===2.3. Restricciones en tensión de tipo local===

Generalmente, la forma más habitual de establecer restricciones en tensión consiste en definirlas en unos puntos predefinidos del dominio de la estructura desde un punto de vista local. En el caso de la optimización topológica de estructuras lo más habitual es introducir restricciones en tensión en el punto central de cada elemento finito [[#bib0080|[16,18,21,22,24,25,29-34,37]]]. Para ello, es habitual utilizar un criterio de rotura para plantear las restricciones en tensión. De este modo se pueden plantear las restricciones de tensión en cada punto como:

<span id='eq0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>g_{e,1}(\boldsymbol{\rho })=\overset{\mbox{ˆ}}{\sigma }({\boldsymbol{\sigma }}^h(\boldsymbol{r}_e^o,\boldsymbol{\rho }))-</math><math>{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}\leq 0,</math>
|-
|<math>g_{e,2}(\boldsymbol{\rho })={\overset{\mbox{ˆ}}{\sigma }}_{m\acute{\imath}n}-</math><math>\overset{\mbox{ˆ}}{\sigma }({\boldsymbol{\sigma }}^h(\boldsymbol{r}_e^o,\boldsymbol{\rho }))\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

donde <math display="inline">\overset{\mbox{ˆ}}{\sigma }({\boldsymbol{\sigma }}^h)</math>  es la tensión de comparación empleada, '''''σ'''''<sup>''h''</sup>  representa el tensor de tensiones calculado, <math display="inline">\boldsymbol{r}_e^o</math>  es el punto de aplicación de la restricción en tensión y '''''ρ'''''  es el vector de variables de diseño.      

Si el material que se utiliza presenta un comportamiento análogo a tracción que a compresión, como por ejemplo el acero o el aluminio, el criterio de rotura más habitual es el de Von Mises. En este caso, las dos restricciones en tensión propuestas en [[#eq0030|(6)]]  se reducen a una única restricción:

<span id='eq0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>g_e(\boldsymbol{\rho })={\overset{\mbox{ˆ}}{\sigma }}_{VM,e}({\boldsymbol{\sigma }}^h(\boldsymbol{r}_e^o,\boldsymbol{\rho }))-</math><math>{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

siendo <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{VM}({\boldsymbol{\sigma }}^h(\boldsymbol{r}_e^o,\boldsymbol{\rho }))</math>  la tensión de comparación de Von Mises obtenida a partir del tensor de tensiones calculado <math display="inline">({\boldsymbol{\sigma }}^h(\boldsymbol{r}_e^o,\boldsymbol{\rho }))</math>  y <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}</math>  el valor máximo de la tensión de comparación.      

En este trabajo se ha utilizado este criterio de comparación para imponer restricciones en tensión dado que algunos de los materiales utilizados habitualmente en ingeniería estructural se rigen por esta tensión de referencia. Además, este mismo método ha sido utilizado en otros trabajos sobre optimización topológica de estructuras por autores como Burger y Stainko [[#bib0080|[16]]], Duysinx y Bendsøe [[#bib0105|[21]]], Duysinx y Sigmund [[#bib0110|[22]]]  o Pereira et al. [[#bib0125|[25]]]. Otros criterios de rotura similares y aplicables a diferentes tipos de materiales pueden encontrarse por ejemplo en [[#bib0100|[20]]].      

Las formulaciones de las restricciones en tensión propuestas en las ecs.[[#eq0030|(6)]]  y [[#eq0035|(7)]]  están sujetas a fenómenos de singularidad debido al carácter discontinuo de la tensión cuando la densidad relativa del material ''ρ''<sub>''e''</sub>  tiende a cero. A medida que el proceso de optimización reduce el valor de la variable de diseño ''ρ''<sub>''e''</sub>  el valor que toma la restricción del elemento en el que se está reduciendo la densidad relativa va cambiando. Sin embargo, cuando la densidad relativa de este elemento se aproxima a cero es posible que la tensión a la que está sometido el punto central del mismo se incremente hasta llegar a incumplirse la restricción en tensión. Sin embargo, si se elimina por completo el material (''ρ''<sub>''e''</sub>  = 0) la restricción en tensión desaparece. En consecuencia, los algoritmos de optimización son en este caso incapaces de eliminar los elementos que presentan restricciones en tensión llevando al mínimo el valor de su densidad relativa aunque las soluciones óptimas así lo requieran. Algunos de estos fenómenos de singularidad aparecen reflejados en los trabajos de Cheng y Jiang [[#bib0090|[18]]], Cheng y Guo [[#bib0095|[19]]]  y Duysinx y Bendsøe [[#bib0105|[21]]], así como algunas de las posibles soluciones utilizadas para evitar este problema.      

En este trabajo se propone una relajación de la formulación de las restricciones basada en los trabajos de Cheng y Guo [[#bib0095|[19]]]  y de Duysinx y Bendsøe [[#bib0105|[21]]]  que consiste en la incorporación de un nuevo término que modifica el valor máximo de la tensión de rotura en función de la densidad relativa del elemento en el que se aplica la restricción. Así, la imposición de las restricciones en tensión incorporando este factor de relajación se puede escribir como:

<span id='eq0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>g_e=\left({\overset{\mbox{ˆ}}{\sigma }}_{VM,e}\left({\boldsymbol{\sigma }}_e^h(\boldsymbol{\rho })\right)-\right. </math><math>\left. {\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}{\varphi }_e\right)\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

donde ''φ''<sub>''e''</sub>  es el factor de relajación que se calcula como:

<span id='eq0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\varphi }_e=1-\varepsilon +\frac{\varepsilon }{{\rho }_e},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

siendo ''ɛ''  el parámetro de relajación de la formulación ([[#fig0020|fig. 4]]) y ''ρ''<sub>''e''</sub>  la densidad relativa del elemento ''e''  en el que se impone la restricción. El parámetro ''ɛ''  toma valores próximos a cero (''ɛ''  ∈ [0,001, 0,1]) para evitar una excesiva relajación de las restricciones.

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr4.jpg|center|360px|Factor de relajación para diferentes valores de ɛ.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 4'''. Factor de relajación para diferentes valores de ''ɛ'' .                  

</span>
|}

Además, este planteamiento de las restricciones en tensión se puede modificar para que tenga en cuenta el uso de tensiones reales o efectivas. Para ello se define una formulación alternativa a la propuesta en [[#eq0040|(8)]]  como:

<span id='eq0050'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>g_e=\left({\overset{\mbox{ˆ}}{\sigma }}_{VM,e}\left({\boldsymbol{\sigma }}_e^h(\boldsymbol{\rho })\right)-\right. </math><math>\left. {\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}{\varphi }_e\right){\rho }_e^q\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

donde el exponente ''q''  es un parámetro que toma los valores ''q''  = 0 cuando se impone la restricción en tensiones reales o ''q''  = 1 cuando se aplique la restricción en tensiones efectivas. Las ventajas e inconvenientes de estos planteamientos se han estudiado con detalle en  [[#bib0120|[24,31,32]]].      

La ventaja más importante de los planteamientos locales de las restricciones es su robustez. Algunos ejemplos de aplicación se han analizado con esta formulación obteniendo resultados satisfactorios. Además, con esta formulación se garantiza en mayor medida la validez de la solución estructural ya que se comprueba el estado tensional en el punto central de cada elemento de la estructura.

Sin embargo, la fiabilidad de la formulación propuesta y de los resultados obtenidos tiene como contrapartida el elevado coste computacional que se requiere para resolver los problemas de optimización resultantes. Los problemas de optimización topológica de estructuras utilizan en la práctica un número elevado de variables de diseño, normalmente tantas como elementos. Por otra parte, la formulación local del problema también obliga a imponer un número elevado de restricciones altamente no lineales en nuestro modelo. En consecuencia, los algoritmos de optimización desarrollados requieren un tiempo de cálculo muy elevado en la práctica. Por este motivo se plantean dos formulaciones adicionales para las restricciones en tensión con la finalidad de reducir las necesidades computaciones del problema de diseño óptimo: una formulación global y una formulación agregada por bloques de las restricciones.

<span id='sec0030'></span>
===2.4. Restricciones en tensión de tipo global===

Tal y como se ha comentado en el apartado anterior la formulación local de las restricciones en tensión requiere unas necesidades computacionales muy elevadas, motivo por el que se han propuesto nuevas formulaciones de las restricciones. La metodología más habitual para conseguir este tipo de planteamientos que requieren un menor coste computacional consiste en establecer una función que aproxime el comportamiento de todas las restricciones locales en una única restricción que recibe, en la terminología de optimización, el nombre de función global o restricción global. Esta función global presenta ventajas computacionales muy claras ya que se reduce considerablemente el número de restricciones del modelo pero, por el contrario, no garantiza de forma estricta el cumplimiento de las tensiones desde un punto de vista local, de modo que es posible que algunas tensiones locales excedan su valor máximo permitido. A pesar de esto, la formulación global de las restricciones actúa evitando que la violación de las restricciones locales sea excesivamente elevada de modo que las tensiones a nivel local satisfagan de forma aproximada el máximo valor permitido.

Como es obvio, la elección de esta función global es uno de los puntos claves para el buen comportamiento de la formulación. En este trabajo se propone la utilización de la función global de Kreisselmeier-Steinhauser como función de agregación de las restricciones locales. Esta función presenta un comportamiento logarítmico-exponencial y ha sido aplicada con éxito en otras disciplinas de la optimización de estructuras [[#bib0190|[38-40]]]. De acuerdo con este planteamiento se puede establecer una restricción global de tensiones como:

<span id='eq0055'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>G_{KS}=\frac{1}{\mu }ln\left[\sum_{j=1}^m{exp}^{\mu g_j}\right]-\frac{1}{\mu }\quad ln(m)\leq 0,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

donde ''m''  es el número de restricciones locales de desigualdad consideradas, ''g''<sub>''j''</sub>  es el valor que toma la restricción local del elemento ''j''  y ''μ''  es el factor o parámetro de agregación.      

Para el problema de optimización topológica que se plantea, la formulación propuesta en la ec.[[#eq0055|(11)]]  puede particularizarse de modo que se eviten, por ejemplo, fenómenos de escala con las unidades de medida de las restricciones así como incorporar factores de relajación para evitar los fenómenos de singularidad de las tensiones. Así, la formulación global propuesta en la ec.[[#eq0055|(11)]]  se reescribiría ahora como:

<span id='eq0060'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>G_{KS}=\frac{1}{\mu }\left[ln\left(\sum_{e=1}^{N_e}{exp}^{\mu ({\overset{\mbox{ˆ}}{\sigma }}_e^\ast -1)}\right)-\right. </math><math>\left. ln(N_e)\right]\leq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overset{\mbox{ˆ}}{\sigma }}_e^\ast =\frac{\overset{\mbox{ˆ}}{\sigma }({\boldsymbol{\sigma }}_e^h)}{{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}\quad {\varphi }_e},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

siendo <math display="inline">{\boldsymbol{\sigma }}_e^h</math>  la aproximación calculada del tensor de tensiones del material en el punto central del elemento ''e'' , ''σ''<sub>''máx''</sub>  la tensión máxima admisible y ''φ''<sub>''e''</sub>  el coeficiente de relajación de las tensiones propuesto en [[#eq0045|(9)]].      

El parámetro ''μ''  debe tomar el mayor valor posible dado que de esta forma la formulación global representará más fielmente el comportamiento de la formulación con restricciones locales. De hecho en el límite la restricción global propuesta en  [[#eq0055|(11)]]  queda como,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\underset{\mu \rightarrow \infty}{lim}G_{KS}=m\acute{a}x{\left\{g_j\right\}}_{j=1,\ldots ,m}\leq 0.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Sin embargo, valores muy elevados de este parámetro generan una función muy no-lineal y difícil de resolver con los métodos numéricos habituales. En la práctica se recomienda utilizar valores del parámetro de agregación en el intervalo ''μ''  ∈ [20, 40]. En [[#bib0155|[31]]]  se puede encontrar un estudio detallado de los valores más adecuados para el parámetro de agregación.      

La formulación global propuesta en [[#eq0060|(12)]]  consigue reducir enormemente las necesidades computacionales del problema de optimización topológica. Sin embargo, cuando el número de elementos y de restricciones locales es muy elevado la formulación global aglutina en una única función la información que proporcionan un número elevado de contribuciones locales. Por este motivo, cuando el número de elementos aumenta de forma importante, el cumplimiento aproximado de las restricciones desde un punto de vista local no se puede garantizar. Además, el agrupamiento de un número elevado de restricciones provoca que los algoritmos de optimización requieran de un número de iteraciones mucho mayor que para la formulación local dado que el efecto de agrupar restricciones hace que el algoritmo no disponga del mismo nivel de información que en la formulación local. En algunos casos, la agrupación de un número muy elevado de restricciones locales provoca que los óptimos obtenidos no sean adecuados. Por este motivo, se plantea una tercera formulación de las restricciones del problema: la formulación agregada por bloques.

===2.5. Reducción por bloques de restricciones de tipo local===

La formulación global puede presentar algunos inconvenientes cuando el número de restricciones aumenta de forma muy importante. Para evitar este efecto, se plantea una formulación por grupos o bloques de elementos. La idea principal de esta formulación consiste en agrupar los elementos en un número definido de bloques e imponer para cada bloque una restricción global similar a la propuesta en [[#eq0060|(12)]]. En la [[#fig0025|figura 5]]  se puede observar la composición de un bloque de elementos sobre el que se impondrá una función global de restricciones.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr5.jpg|center|214px|Definición de un bloque de elementos de la malla.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 5'''. Definición de un bloque de elementos de la malla.

</span>
|}

Con esta formulación el número de restricciones totales de cada problema es igual al número de bloques en que dividamos la malla de elementos finitos. Es obvio que este modelo incluye también a las formulaciones presentadas en los apartados [[#sec0025|2.3]]  y [[#sec0030|2.4]]  si se define un número de bloques adecuado. En el caso de la formulación con restricciones locales propuesta en el apartado [[#sec0025|2.3]]  bastará con tomar tantos bloques como elementos tenga la malla. Si, por el contrario, definimos un único bloque con todos los elementos de la malla estaremos utilizando la formulación global de las restricciones propuesta en el apartado [[#sec0030|2.4]]. Por lo tanto, un mayor número de bloques proporcionará soluciones mejores dado que el control sobre las tensiones es más estricto. Una vez fijado el número de bloques a utilizar en el problema, el número de elementos de la malla se distribuye de forma homogénea de modo que el número de elementos por bloque sea, aproximadamente, el mismo para todos.      

Esta formulación de las restricciones en tensión por bloques presenta las ventajas que ofrece la formulación local porque establece un control más efectivo sobre las restricciones en tensión. Por otra parte, debido al menor número de restricciones, también se consiguen reducir considerablemente las necesidades computacionales del problema de optimización, sobre todo en tiempo de ejecución.

Cada restricción agrupada por bloques se puede, por tanto, escribir como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>G_{KS,b}=\frac{1}{\mu }\left[ln\left(\sum_{e\in B_b}{exp}^{\mu ({\overset{\mbox{ˆ}}{\sigma }}_e^\ast -1)}\right)-\right. </math><math>\left. ln(N_e^b)\right]\leq 0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

siendo ''B''<sub>''b''</sub>  el conjunto de elementos de la malla que pertenecen al bloque ''b  '' , <math display="inline">N_e^b</math>  el número de elementos de la malla que pertenecen al bloque ''b''  y ''μ''  el parámetro de agregación propuesto en  [[#eq0055|(11)]]  y [[#eq0060|(12)]].      

Uno de los puntos claves de este planteamiento de las restricciones radica en el modo de definir la distribución y configuración de los bloques de elementos. Parece razonable suponer que la forma de generación de los bloques de elementos a partir de la malla original es un condicionante muy importante a tener en cuenta en esta formulación dado que afecta a las funciones globales y, en consecuencia, a las restricciones en tensión que definen. Esta distribución de los bloques, o lo que es lo mismo, la configuración de los elementos de cada bloque, puede realizarse de muy diversas formas: criterios de proximidad (por ejemplo, en la [[#fig0025|figura 5]]), geométricos (forma del bloque) o incluso de aleatoriedad (seleccionando los elementos de cada bloque de forma aleatoria). En este trabajo se ha utilizado una metodología para definir los bloques de elementos basada en la numeración de los elementos de la malla utilizada, de modo que los elementos de cada bloque se eligen de forma consecutiva siguiendo el orden numérico establecido para realizar el cálculo estructural. Esta técnica genera bloques con formas alargadas para las aplicaciones prácticas en optimización topológica. A pesar de esto, este método proporciona soluciones altamente satisfactorias ya que se ha comprobado en la práctica que la distribución de los bloques no supone un condicionante importante para el diseño final de la estructura. De hecho, en ningún ejemplo de los que se han estudiado se ha detectado una incidencia perceptible de la distribución de los bloques en las soluciones óptimas obtenidas.      

<span id='sec0040'></span>

==3. Cálculo de estructuras mediante el método de los elementos finitos==

Una vez planteado el problema de optimización es necesario desarrollar los puntos más importantes de la formulación. En primer lugar, será necesario realizar el cálculo de estructuras porque a partir de él se desarrollarán todas las fases posteriores del proceso de optimización.

El problema de optimización topológica de estructuras continuas que se propone requiere un método fiable de cálculo de estructuras para obtener el comportamiento resistente de las distintas soluciones planteadas durante todo el proceso de diseño. Esta técnica de cálculo debe, además, permitir la incorporación de las variables de diseño, las densidades relativas, en su formulación.

Una formulación del método de los elementos finitos será el modelo de cálculo del comportamiento resistente utilizado en este caso. El problema se estudiará bajo las hipótesis de elasticidad lineal, pequeños desplazamientos y pequeños gradientes de desplazamientos. Además, esta formulación nos permitirá obtener, además de la respuesta estructural, información acerca de otros aspectos vitales para el problema de optimización tales como el peso de la estructura, el valor de las restricciones e incluso el análisis de sensibilidad.

El modelo de cálculo que se propone parte de un planteamiento estándar del método [[#bib0120|[24,41,42,43]]]  sobre el que se aplican las modificaciones necesarias para incorporar el efecto de las variables de diseño del problema de optimización.      

A continuación, se propone un desarrollo general del análisis estructural sin considerar el efecto de la densidad relativa y el modelo propuesto para este problema de optimización con la finalidad de comparar ambas formulaciones y establecer las modificaciones que es necesario incorporar en una formulación convencional del problema.

<span id='sec0045'></span>
===3.1. Planteamiento del problema estructural===

Sea Ω<sup>''o''</sup>  un dominio del espacio material que está ocupado por un cuerpo que, bajo la acción de cargas exteriores aplicadas, se deformará ocupando un dominio Ω diferente. Entonces, cualquier punto ''P''<sup>''o''</sup>  del dominio Ω<sup>''o''</sup>  pasará a ocupar un punto ''P''  del dominio Ω una vez que se ha deformado el cuerpo. Si definimos como '''''r'''''<sup>''o''</sup>  y '''''r'''''  a los vectores de posición de ''P''<sup>''o''</sup>  y ''P''  podemos calcular los desplazamientos de cada punto según:

<span id='eq0080'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}(\boldsymbol{r}^o)=\boldsymbol{r}(\boldsymbol{r}^o)-</math><math>\boldsymbol{r}^o.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Una vez conocidos estos desplazamientos podemos obtener las deformaciones '''''ɛ''''' ('''''r'''''<sup>''o''</sup>) y las tensiones '''''σ''''' ('''''r'''''<sup>''o''</sup>) en cada punto ''P''<sup>''o''</sup>  como:

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\varepsilon }=\boldsymbol{Lu},\quad \boldsymbol{\sigma }=</math><math>\boldsymbol{D}\boldsymbol{\varepsilon }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

donde '''''L'''''  es la matriz de deformación-desplazamientos y '''''D'''''  es la matriz constitutiva del material. Los vectores '''''ɛ'''''  y '''''σ'''''  son los vectores de deformaciones y de tensiones respectivamente para cada punto (Hughes [[#bib0205|[41]]]).      

Para unas cargas externas aplicadas por unidad de volumen '''''b''''' ('''''r'''''<sup>''o''</sup>) sobre el cuerpo en el dominio Ω<sup>''o''</sup>  y unas cargas por unidad de superficie '''''t''''' ('''''r'''''<sup>''o''</sup>) en el contorno <math display="inline">{\Gamma }_{\sigma }^o</math>  el análisis estructural consiste en

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Hallar\quad \boldsymbol{u}\in H_u</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>tal\quad que\quad a(\boldsymbol{w},\boldsymbol{u})=</math><math>(\boldsymbol{w},\boldsymbol{b})_{{\Omega }^o}+(\boldsymbol{w},\boldsymbol{t})_{{\Gamma }_{\sigma }^o}\quad \forall \boldsymbol{w}\in H_w</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>siendo\quad a(\boldsymbol{w},\boldsymbol{u})=\iiint_{{\Omega }^o}(\boldsymbol{L}\boldsymbol{w})^T\boldsymbol{D}(\boldsymbol{Lu})d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\quad (\boldsymbol{w},\boldsymbol{b})_{{\Omega }^o}=</math><math>\iiint_{{\Omega }^o}\boldsymbol{w}^T\boldsymbol{b}d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18d)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\quad (\boldsymbol{w},\boldsymbol{t})_{{\Gamma }_{\sigma }^o}=</math><math>\iint_{{\Gamma }_{\sigma }^o}\boldsymbol{w}^T\boldsymbol{t}d\Gamma ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18e)
|}

donde las funciones de prueba '''''u'''''  y las funciones de test <math display="inline">\boldsymbol{w}</math>  satisfacen las condiciones contorno esenciales (desplazamientos prescritos) y las correspondientes condiciones de contorno homogéneas, respectivamente (Navarrina et al. [[#bib0120|[24]]]).      

<span id='sec0050'></span>

===3.2. Modelo numérico de elementos finitos===

La resolución analítica del problema estructural planteado anteriormente no es en general posible, lo que nos obliga a aproximar la solución exacta a través de una discretización finita del dominio. Para ello será necesario remplazar los espacios continuos de las funciones ''H''<sub>''u''</sub>  y <math display="inline">H_w</math>  por los correspondientes espacios discretizados <math display="inline">H_u^h</math>  y <math display="inline">H_w^h</math>  y por las respectivas funciones de prueba '''''u'''''<sup>''h''</sup>  y de test <math display="inline">\boldsymbol{w}^h</math>  de los mismos. Sean, por tanto, <math display="inline">\left\{{\Phi }_i(\boldsymbol{r}^o)\right\}</math>  y <math display="inline">\left\{w_j(\boldsymbol{r}^o)\right\}</math>  unas bases convenientemente elegidas de las funciones de prueba y de test en los correspondientes subespacios <math display="inline">H_u^h</math>  y <math display="inline">H_w^h</math> , que verifican las condiciones de contorno homogéneas del problema. La solución aproximada del problema estructural ([[#sec0045|3.1]]) se puede plantear como

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}^h(\boldsymbol{r}^o)=\boldsymbol{u}^p(\boldsymbol{r}^o)+</math><math>\sum_{i=1}^N{\boldsymbol{\Phi }}_i(\boldsymbol{r}^o){\boldsymbol{\alpha }}_i,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{w}^h(\boldsymbol{r}^o)=\sum_{j=1}^N\boldsymbol{W}_j(\boldsymbol{r}^o){\boldsymbol{\beta }}_j,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19b)
|}

siendo

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\Phi }}_i(\boldsymbol{r}^o)={\phi }_i(\boldsymbol{r}^o)\boldsymbol{I}_3,\quad \quad \boldsymbol{W}_j(\boldsymbol{r}^o)=</math><math>w_j(\boldsymbol{r}^o)\boldsymbol{I}_3,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

donde <math display="inline">{{\boldsymbol{\alpha }}_i}_{i=1,\ldots ,N}</math>  es el vector de desplazamientos nodales para cada nodo ''i''  y '''''u'''''<sup>''h''</sup> ('''''r'''''<sup>''o''</sup>) el vector de desplazamientos de un punto de coordenadas '''''r'''''<sup>''o''</sup> .      

La discretización espacial del dominio se realiza de modo que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\Omega }^o={\bigcup }_{e=1}^{N_e}{\Omega }_e,\quad {\Omega }_{e_1}\bigcap {\Omega }_{e_2}=</math><math>\phi \quad \forall \quad e_1\not =e_2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

donde Ω<sub>''e''</sub>  es el dominio ocupado por el elemento ''e''  y ''N''<sub>''e''</sub>  es el número de elementos finitos utilizados. Las bases de funciones de prueba y de test se plantean habitualmente, en problemas de cálculo de estructuras, mediante formulaciones isoparamétricas de tipo Galerkin, de modo que

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>w_j(\boldsymbol{r}^o)={\phi }_j(\boldsymbol{r}^o).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

Consecuentemente el modelo de elementos finitos con aproximación de tipo Galerkin consiste en:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Hallar\quad \boldsymbol{\alpha }=\left\{{\boldsymbol{\alpha }}_i\right\},\quad i=</math><math>1,\ldots ,N</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>tal~que\quad \sum_{i=1}^N\boldsymbol{K}_{ji}{\boldsymbol{\alpha }}_i=</math><math>\boldsymbol{f}_j,\quad j=1,\ldots ,N</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>siendo\quad \boldsymbol{K}_{ji}=a({\boldsymbol{\Phi }}_j,{\boldsymbol{\Phi }}_i),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\quad \boldsymbol{f}_j=({\boldsymbol{\Phi }}_j,\boldsymbol{b})_{{\Omega }^o}+</math><math>({\boldsymbol{\Phi }}_j,\boldsymbol{t})_{{\Gamma }_{\sigma }^o}-</math><math>a({\boldsymbol{\Phi }}_j,\boldsymbol{u}^p).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23d)
|}

donde los términos '''''K'''''<sub>''ji''</sub>  y ''f''<sub>''j''</sub>  pueden obtenerse ensamblando las contribuciones elementales de cada elemento como

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}_{ji}=\sum_{e=1}^{N_e}\boldsymbol{K}_{ji}^e,\quad \boldsymbol{f}_j=</math><math>\iint_{{\Gamma }_{\sigma }^o}{\boldsymbol{\Phi }}_j^T\boldsymbol{t}d\Gamma +</math><math>\sum_{e=1}^{N_e}\boldsymbol{f}_j^e,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

siendo

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}_{ji}^e=\iiint_{{\Omega }_e}(\boldsymbol{L}{\boldsymbol{\Phi }}_j)^T\boldsymbol{D}(\boldsymbol{L}{\boldsymbol{\Phi }}_i)d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{f}_j^e=\iiint_{{\Omega }_e}({\boldsymbol{\Phi }}_j^T\boldsymbol{b}-</math><math>(\boldsymbol{L}{\boldsymbol{\Phi }}_j)^T\boldsymbol{D}({\boldsymbol{Lu}}^p))d\Omega .</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25b)
|}

Una vez obtenida la solución al problema ([[#sec0050|3.2]]), se puede calcular para cualquier punto '''''r'''''<sup>''o''</sup>  ∈ Ω<sup>''o''</sup>  la aproximación de los desplazamientos, deformaciones y tensiones de acuerdo con la ec. [[#eq0085|(17)]]  y la Sección [[#sec0050|3.2]] como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{u}^h(\boldsymbol{r}^o)=\boldsymbol{u}^p(\boldsymbol{r}^o)+</math><math>\sum_{i=1}^N{\boldsymbol{\Phi }}_i(\boldsymbol{r}^o){\boldsymbol{\alpha }}_i,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\varepsilon }}^h(\boldsymbol{r}^o)=\boldsymbol{L}\boldsymbol{u}^h(\boldsymbol{r}^o),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\sigma }}^h(\boldsymbol{r}^o)=\boldsymbol{D}{\boldsymbol{\varepsilon }}^h(\boldsymbol{r}^o).</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26c)
|}

<span id='sec0055'></span>

===3.3. Análisis estructural de material con densidad relativa===

Para realizar el análisis estructural incorporando la densidad relativa de cada elemento seguiremos un esquema similar al planteado en el apartado anterior. Partimos de un dominio Ω<sup>''o''</sup>  que ahora estará ocupado por un material poroso y definimos la variable de diseño ''ρ'' ('''''r'''''<sup>''o''</sup>) como la densidad relativa del material en un punto ''P''<sup>''o''</sup>  de coordenadas '''''r'''''<sup>''o''</sup>.      

De acuerdo con [[#eq0080|(16)]] , el análisis estructural pretenderá obtener el campo de desplazamientos así como las deformaciones y tensiones asociadas para una configuración material dada.      

Si definimos como ''d'' Ω el volumen ocupado por una región diferencial de material entorno al punto ''P''<sup>''o''</sup> , el volumen de material sólido existente en esta región diferencial es ''ρ'' ('''''r'''''<sup>''o''</sup>) ''d'' Ω debido al efecto de la densidad relativa. De este modo, el análisis realizado en ([[#sec0050|3.2]]) puede reescribirse en la forma:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Dado\quad \rho ({\Omega }^o)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>hallar\quad \boldsymbol{u}\in H_u</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27b)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>tal\quad que\quad a(\boldsymbol{w},\boldsymbol{u})=</math><math>(\boldsymbol{w},\boldsymbol{b})_{{\Omega }^o}+(\boldsymbol{w},\boldsymbol{t})_{{\Gamma }_{\sigma }^o}\quad \forall \boldsymbol{w}\in H_w</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27c)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>siendo\quad a(\boldsymbol{w},\boldsymbol{u})=\iiint_{{\Omega }^o}(\boldsymbol{L}\boldsymbol{w})^T\boldsymbol{D}(\boldsymbol{Lu})\rho d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27d)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\quad (\boldsymbol{w},\boldsymbol{b})_{{\Omega }^o}=</math><math>\iiint_{{\Omega }^o}\boldsymbol{w}^T\boldsymbol{b}\rho d\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27e)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\quad (\boldsymbol{w},\boldsymbol{t})_{{\Gamma }_{\sigma }^o}=</math><math>\iint_{{\Gamma }_{\sigma }^o}\boldsymbol{w}^T\boldsymbol{t}d\Gamma ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27f)
|}

Como se puede observar, los cambios que se introducen en la formulación se reducen a tener en cuenta el efecto de la porosidad en la integración de los términos elementales correspondientes. Una vez calculadas estas contribuciones elementales, las deformaciones y tensiones se calculan del mismo modo que en ([[#sec0055|3.3]]). Sin embargo, es importante destacar que el análisis anterior carece de sentido físico cuando la densidad relativa es nula. Consecuentemente, los desplazamientos, deformaciones y tensiones tampoco existen.

===3.4. Modelo numérico de elementos finitos con densidad relativa===

El análisis estructural por elementos finitos para un material que incorpore el efecto de la densidad relativa se puede plantear a partir de ([[#sec0050|3.2]]). Según ([[#sec0055|3.3]]), una vez conocido el vector '''''ρ'''''  que contiene el valor de las variables de diseño, el problema a resolver consiste en

<span id='eq0220'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Hallar\quad \boldsymbol{\alpha }(\boldsymbol{\rho })</math>
|-
|<math>tal~que\quad \sum_{i=1}^N\boldsymbol{K}_{ji}(\boldsymbol{\rho }){\boldsymbol{\alpha }}_i(\boldsymbol{\rho })=</math><math>\boldsymbol{f}_j(\boldsymbol{\rho }),\quad j=1,\ldots ,N.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

Al igual que en la formulación sin densidad relativa los términos necesarios para resolver [[#eq0220|(28)]]  se pueden calcular ensamblando las contribuciones de cada elemento de forma independiente de modo que:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}_{ji}(\boldsymbol{\rho })=\sum_{e=1}^{N_e}\boldsymbol{K}_{ji}^e(\boldsymbol{\rho }),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{f}_j(\boldsymbol{\rho })=\iint_{{\Gamma }_{\sigma }^o}{\boldsymbol{\Phi }}_j^T\boldsymbol{t}d\Gamma +</math><math>\sum_{e=1}^{N_e}\boldsymbol{f}_j^e(\boldsymbol{\rho }),</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29b)
|}

siendo las contribuciones elementales

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}_{ji}^e({\rho }_e)=\iiint_{{\Omega }_e}(\boldsymbol{L}{\boldsymbol{\Phi }}_j)^T\boldsymbol{D}(\boldsymbol{L}{\boldsymbol{\Phi }}_i){\rho }_ed\Omega ,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30a)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{f}_j^e({\rho }_e)=\iiint_{{\Omega }_e}({\boldsymbol{\Phi }}_j^T\boldsymbol{b}-</math><math>(\boldsymbol{L}{\boldsymbol{\Phi }}_j)^T\boldsymbol{D}({\boldsymbol{Lu}}^p)){\rho }_ed\Omega .</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30b)
|}

Una vez que se ha obtenido la solución '''''α''''' ('''''ρ''''') del problema [[#eq0220|(28)]]  se pueden calcular para cualquier punto ''P''<sup>''o''</sup>  las aproximaciones de los desplazamientos, de las deformaciones y de las tensiones a partir de las expresiones propuestas en ([[#sec0050|3.2]]).      

<span id='sec0065'></span>

==4. Algoritmo de optimización==

La resolución de los problemas de optimización topológica de estructuras es, en general, muy complicada debido al elevado número de variables de diseño que presentan así como a la no-linealidad de la función objetivo y/o de las restricciones. Consecuentemente, los algoritmos y modelos que se utilizan para resolver este tipo de planteamientos deben proporcionar soluciones adecuadas con un coste computacional asumible. Ambos aspectos son fundamentales en la práctica.

Estos dos factores se acentúan en mayor medida en el problema que se presenta dado que, además de utilizar un número elevado de variables de diseño, incorpora una función objetivo no-lineal y un número elevado de restricciones en tensión altamente no-lineales.

En este trabajo se ha utilizado un algoritmo basado en programación lineal secuencial (SLP) dado que este tipo de métodos han demostrado ser eficaces en la resolución de problemas de diseño óptimo de estructuras con restricciones como los que se proponen en [[#bib0220|[44-47]]].      

===4.1. Algoritmo SLP-QLS (programación lineal secuencial con búsqueda unidireccional cuadrática)===

El algoritmo de programación lineal secuencial basa su funcionamiento en el aprovechamiento de las propiedades que ofrecen los espacios linealizados de las restricciones y de la función objetivo, sobre los que se aplican algoritmos como por ejemplo, el Simplex [[#bib0240|[48,49]]], para obtener la solución óptima. De forma general, la solución se obtiene de forma iterativa estudiando el comportamiento de la función objetivo y de las restricciones mediante una aproximación lineal. En cada iteración se pretende obtener una solución mejorada en dos fases. En primer lugar se obtiene mediante esta aproximación lineal la dirección de modificación del diseño que produce un mayor descenso en la función objetivo sin violar las restricciones para, a continuación en un segundo paso, determinar el factor de avance más adecuado para modificar el diseño en la dirección obtenida previamente. Por lo tanto, el vector de variables de diseño se actualiza en cada iteración como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\rho }}^{k+1}={\boldsymbol{\rho }}^k+</math><math>{\theta }^k\boldsymbol{s}^k,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

donde '''''ρ'''''<sup>''k''</sup>  es el vector de variables de diseño y ''θ''<sup>''k''</sup>  es el factor de avance en la dirección de modificación del diseño '''''s'''''<sup>''k''</sup>  en la iteración ''k''.      

La dirección de avance ('''''s'''''<sup>''k''</sup>) se calcula a través de diferentes métodos que utilizan la aproximación lineal del problema. En general, la obtención de esta dirección de modificación es la parte más importante y compleja del algoritmo y requiere de técnicas específicas y complejas debido al elevado número de variables de diseño y al elevado número de restricciones en tensión que pueden presentar los problemas de optimización topológica que se plantean. Para obtener esta dirección de avance se utiliza el método de programación lineal (Simplex) porque ha demostrado ser adecuado para este tipo de problemas. Sin embargo, este algoritmo debe ser complementado con otros más específicos. Así, por ejemplo, cuando el problema no presenta ninguna restricción activa se utilizará el método de máximo descenso y cuando existan restricciones fuertemente violadas se utilizará un algoritmo de entrada a la región factible. Por lo tanto, el algoritmo completo sigue un comportamiento diferente en función del diseño y del estado tensional de la estructura en cada iteración.      

===4.2. Algoritmo de máximo descenso===

Tal y como se ha comentado anteriormente, si el problema no presenta restricciones activas, el algoritmo Simplex no es aplicable y es necesario utilizar temporalmente otro algoritmo hasta que alguna restricción se active.

Debido a la simplicidad del problema cuando no hay restricciones en tensión activas y a que, en general, esta situación de optimización incondicionada es temporal, se ha optado por utilizar la dirección de máximo descenso de la función objetivo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_{i,o}^k=-{\frac{dF}{d{\rho }_i}}_{{\rho }_i^k}\Delta {\rho }_i^k,\quad i=</math><math>1,\ldots ,N_e,</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

donde <math display="inline">s_{i,o}^k</math>  es la componente ''i  ''  de la dirección de modificación del diseño sin normalizar (<math display="inline">\boldsymbol{s}_o^k</math>) en la iteración ''k''  y Δ'''''ρ'''''<sup>''k''</sup>  es el vector que indica la máxima modificación que puede experimentar el vector de variables de diseño '''''ρ'''''  en cada iteración teniendo en cuenta las restricciones laterales. La dirección de modificación del diseño '''''s'''''<sup>''k''</sup>  se obtiene como <math display="inline">\boldsymbol{s}^k=\boldsymbol{s}_o^k/\parallel \boldsymbol{s}_o^k\parallel </math>.      

Este algoritmo debe además complementarse con las correspondientes restricciones laterales para evitar que el algoritmo se detenga (debido a factores de avance nulos) antes de alcanzar el óptimo deseado.

===4.3. Algoritmo de programación lineal (Simplex)===

El algoritmo de programación lineal es el núcleo de todo el proceso de optimización dado que debe proporcionar una dirección de modificación del diseño que minimice la función objetivo pero que tenga en consideración todas y cada una de las restricciones en tensión y de las restricciones laterales.

El algoritmo basa su funcionamiento en la aplicación del algoritmo Simplex sobre el espacio linealizado. Para ello planteamos el desarrollo en serie de Taylor de primer orden de la función objetivo y las restricciones en torno a la solución actual ('''''ρ'''''<sup>''k''</sup>) como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F({\boldsymbol{\rho }}^{k+1})\approx F({\boldsymbol{\rho }}^k)+</math><math>{\frac{dF}{d\boldsymbol{\rho }}}_{{\boldsymbol{\rho }}^k}{\theta }^k\boldsymbol{s}^k</math>
|-
|<math>g_j({\boldsymbol{\rho }}^{k+1})\approx g_j({\boldsymbol{\rho }}^k)+</math><math>{\frac{dg_j}{d\boldsymbol{\rho }}}_{{\boldsymbol{\rho }}^k}{\theta }^k\boldsymbol{s}^k\quad j=</math><math>1,\ldots ,m.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

Además, también es necesario tener en cuenta las restricciones laterales de las variables de diseño, tanto las inferiores como las superiores. El algoritmo Simplex se aplica, por tanto, sobre este problema linealizado. El desarrollo completo de este algoritmo se recoge de forma detallada en [[#bib0155|[31,48]]].      

===4.4. Dirección de entrada a la región factible===

A pesar de la robustez del algoritmo de programación lineal es necesario recurrir en ocasiones a otros algoritmos. Así, por ejemplo cuando el diseño presenta restricciones fuertemente violadas es necesario disponer de un algoritmo que proporcione direcciones de avance que tiendan a devolver el diseño actual a la región factible de diseño.

Este algoritmo de dirección de entrada basa su funcionamiento en la utilización de los gradientes de las restricciones más violadas para obtener la dirección de avance más adecuada, que se obtiene como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>s_{i,o}^k=-\sum_{g_j\in g_v}g_j\frac{{\frac{dg_j}{d{\rho }_i}}_{{\rho }_i^k}}{{\left\|{\frac{dg_j}{d\boldsymbol{\rho }}}_{{\boldsymbol{\rho }}^k}\right\|}^2}\Delta {\rho }_i^k,\quad i=</math><math>1,\ldots ,N_e</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

donde <math display="inline">\boldsymbol{s}_{i,o}^k</math>  es la componente ''i''  del vector de dirección de modificación del diseño sin normalizar de la iteración ''k  '', <math display="inline">g_v</math>  es el conjunto de restricciones violadas del problema y Δ'''''ρ'''''<sup>''k''</sup>  es el vector que indica la máxima modificación que puede experimentar el vector de variables de diseño '''''ρ'''''  en cada iteración teniendo en cuenta las restricciones laterales. La dirección de modificación del diseño '''''s'''''<sup>''k''</sup>  se obtiene como <math display="inline">\boldsymbol{s}^k=\boldsymbol{s}_o^k/\parallel \boldsymbol{s}_o^k\parallel </math> .      

===4.5. Búsqueda del factor de avance===

Una vez calculada y normalizada la dirección de modificación del diseño ('''''s'''''<sup>''k''</sup>) es necesario calcular el factor de avance más adecuado (''θ''<sup>''k''</sup>). Para ello realizaremos una búsqueda unidireccional cuadrática en la dirección obtenida previamente. De este modo, se pretenden evitar los fenómenos oscilatorios que se producen habitualmente en torno a la solución cuando se utilizan aproximaciones lineales [[#bib0230|[46]]].      

Para llevar a cabo el cálculo del factor de avance necesitamos conocer, por tanto, las derivadas de segundo orden tanto de la función objetivo como de las restricciones. Dado que la dirección de modificación del diseño es conocida podremos evitar el cálculo de las derivadas completas de segundo orden realizando un cálculo direccional de estas derivadas en la dirección de modificación del diseño. A partir de estas derivadas planteamos el desarrollo en serie de Taylor de segundo orden de la función objetivo y de las restricciones como:

<span id='eq0265'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>F({\boldsymbol{\rho }}^{k+1})\approx F({\boldsymbol{\rho }}^k)+</math><math>{D_s\quad F}_{{\boldsymbol{\rho }}^k}{\theta }^k+{D_{ss}^2\quad F}_{{\boldsymbol{\rho }}^k}\frac{({\theta }^k)^2}{2}</math>
|-
|<math>g_j({\boldsymbol{\rho }}^{k+1})\approx g_j({\boldsymbol{\rho }}^k)+{D_s\quad g_j}_{{\boldsymbol{\rho }}^k}{\theta }^k+{D_{ss}^2\quad g_j}_{{\boldsymbol{\rho }}^k}\frac{({\theta }^k)^2}{2}j=1,\ldots ,m.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

Una vez planteado el desarrollo en serie de Taylor de segundo orden procedemos a calcular el factor de avance. Para ello, es necesario calcular previamente el valor mínimo y máximo que puede tomar el factor de avance que hace que no se incumplan las restricciones laterales de las variables de diseño. Una vez definidos estos límites, el algoritmo comprueba el valor de la función objetivo y de las restricciones en puntos equiespaciados suficientemente próximos dentro del rango admisible de valores del factor de avance y selecciona el más adecuado.

<span id='sec0095'></span>
==5. Análisis de sensibilidad de las restricciones==

Tal y como se ha planteado anteriormente los algoritmos de optimización propuestos requieren conocer las derivadas completas de primer orden y las derivadas de primer y segundo orden direccionales de las restricciones en tensión.

===5.1. Análisis de sensibilidad de primer orden de las restricciones===

El análisis de sensibilidad de las restricciones se desarrolla de forma analítica mediante el método de la variable adjunta. En este trabajo se desarrolla de forma específica este análisis para la formulación local de las restricciones en tensión. Un desarrollo más detallado de este análisis tanto para la formulación local como para las formulaciones global y agregada por bloques puede encontrarse en [[#bib0155|[31,50]]].      

La formulación local de las restricciones en tensión establece, de acuerdo con [[#eq0050|(10)]], una restricción en el punto central de cada elemento de la malla de elementos finitos utilizada. Esta formulación se puede plantear como:

<span id='eq0270'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>g_e(\boldsymbol{\rho })={{\tilde{g}}_e({\overset{\mbox{ˆ}}{\sigma }}_{VM,e},{\rho }_e)}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

donde <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{VM,e}</math>  es la tensión de comparación de Von Mises y '''''σ'''''<sub>''e''</sub>  es un vector que contiene las componentes del tensor de tensiones, ambos calculados en el punto central del elemento ''e''. El vector de desplazamientos nodales es '''''α''''' ('''''ρ'''''). Por simplicidad, se omite para el análisis de sensibilidad el uso del superíndice ''h''  introducido en el apartado ([[#sec0050|3.2]]) para indicar que las variables representan valores calculados mediante el método de elementos finitos.      

Las derivadas de primer orden de las restricciones en tensión con respecto a las variables de diseño se pueden plantear aplicando la regla de la cadena como:

<span id='eq0275'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{dg_e}{d{\rho }_i}={\frac{\partial {\tilde{g}}_e}{\partial {\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}{d{\boldsymbol{\sigma }}_e}}_{\begin{array}{l}
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\boldsymbol{\sigma }}_e}{d\boldsymbol{\alpha }}}_{\boldsymbol{\alpha }(\boldsymbol{\rho })}\quad \frac{d\boldsymbol{\alpha }}{d{\rho }_i}+</math><math>{{\delta }_{ie}\frac{\partial {\tilde{g}}_e}{\partial {\rho }_i}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

donde <math display="inline">\frac{dg_e}{d\boldsymbol{\rho }}=</math><math>\left\{\frac{dg_e}{d{\rho }_i}\right\}\quad i=1,\ldots ,N_e</math> , siendo ''N''<sub>''e''</sub>  el número de variables de diseño. La función ''δ''<sub>''ie''</sub>  es la delta de Kronecker.      

Una vez establecidas las derivadas de las restricciones de acuerdo con [[#eq0275|(37)]]  es necesario calcular las derivadas de los desplazamientos nodales '''''α''''' ('''''ρ''''') con respecto a la variable de diseño ''ρ''<sub>''i''</sub> . Estas derivadas también se pueden obtener de forma analítica derivando el sistema de ecuaciones lineales, definido en el cálculo de estructuras por el método de elementos finitos, propuesto en ([[#sec0055|3.3]]) como:

<span id='eq0280'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}\frac{d\boldsymbol{\alpha }}{d{\rho }_i}=</math><math>\frac{d\boldsymbol{f}}{d{\rho }_i}-\frac{d\boldsymbol{K}}{d{\rho }_i}\boldsymbol{\alpha }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

El cálculo de las derivadas de los desplazamientos requiere que se resuelva un sistema de ecuaciones para cada variable de diseño y caso de carga aplicado cuyas incógnitas son las derivadas buscadas. La matriz de coeficientes de este sistema de ecuaciones es la matriz de rigidez obtenida en el cálculo estructural. Las contribuciones que definen el vector de términos independientes del sistema de ecuaciones se pueden obtener de forma directa a partir del modelo de elementos finitos con densidad relativa. Si se tiene en cuenta que las cargas externas puntuales y las cargas externas por unidad de área no dependen para este problema de las variables de diseño, sus derivadas serán obviamente nulas. Por el contrario, las fuerzas por unidad de volumen, como por ejemplo el peso, presentan una dependencia lineal con respecto a las variables de diseño. Así, a nivel elemental las derivadas de las fuerzas por unidad de volumen dependen de la densidad relativa del elemento considerado [[#bib0155|[31,32,50]]]. Por tanto, las derivadas de las fuerzas nodales ('''''f'''''<sub>''k''</sub> , ''k''  = 1, …, ''N'') con respecto a la densidad relativa del elemento ''e''  pueden obtenerse ensamblando las contribuciones elementales de las fuerzas por unidad de volumen considerando que la densidad relativa del elemento ''e''  es ''ρ''<sub>''e''</sub>  = 1 y el resto son nulas (''ρ''<sub>''i''</sub>  = 0 ∀ ''i''  ≠ ''e'' , ''i''  = 1, …, ''N''<sub>''e''</sub>).      

Del mismo modo las derivadas de la matriz de rigidez '''''K'''''  con respecto a la variable de diseño ''ρ''<sub>''e''</sub>  se puede obtener ensamblando las matrices elementales '''''K'''''<sub>''i''</sub>  con ''ρ''<sub>''i''</sub>  = 1 si ''i''  = ''e''  y ''ρ''<sub>''i''</sub>  = 0 en cualquier otro caso [[#bib0250|[50]]]. En la práctica este cálculo se puede realizar a nivel elemental y ensamblar los resultados en el vector de términos independientes.      

Una vez planteada la obtención de los términos independientes de acuerdo con los planteamientos anteriores, las derivadas de los desplazamientos nodales se pueden calcular mediante dos métodos distintos: diferenciación directa analítica o diferenciación por el método de la variable adjunta [[#bib0155|[31,50]]].      

El método de diferenciación directa consiste en calcular las derivadas de los desplazamientos nodales resolviendo los sistemas de ecuaciones propuestos en [[#eq0280|(38)]]  y sustituir estas derivadas en [[#eq0275|(37)]], obteniendo así las derivadas de las restricciones. En la práctica este método requiere resolver tantos sistemas de ecuaciones adicionales como variables de diseño tenga el problema de optimización para cada caso de carga.      

En este trabajo el método de derivación utilizado es el método de la variable adjunta. Este planteamiento consiste en reemplazar las derivadas de los desplazamientos nodales obtenidas de acuerdo con [[#eq0280|(38)]]  en [[#eq0275|(37)]], con lo que se obtiene que:

<span id='eq0285'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{dg_e}{d{\rho }_i}={\boldsymbol{\lambda }}_e^T\left(\frac{d\boldsymbol{f}}{d{\rho }_i}-\right. </math><math>\left. \frac{d\boldsymbol{K}}{d{\rho }_i}\boldsymbol{\alpha }\right)+</math><math>{{\delta }_{ie}\frac{\partial {\tilde{g}}_e}{\partial {\rho }_i}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\boldsymbol{\lambda }}_e^T={\frac{\partial {\tilde{g}}_e}{\partial {\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}{d{\boldsymbol{\sigma }}_e}}_{\begin{array}{l}
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\boldsymbol{\sigma }}_e}{d\boldsymbol{\alpha }}}_{\boldsymbol{\alpha }(\boldsymbol{\rho })}\quad \boldsymbol{K}^{-1}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

la variable adjunta, que se obtiene como resultado de resolver el sistema de ecuaciones:

<span id='eq0295'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}^T{\boldsymbol{\lambda }}_e={\left({\frac{\partial {\tilde{g}}_e}{\partial {\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}{d{\boldsymbol{\sigma }}_e}}_{\begin{array}{l}
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\boldsymbol{\sigma }}_e}{d\boldsymbol{\alpha }}}_{\boldsymbol{\alpha }(\boldsymbol{\rho })}\right)}^T.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Todos los términos que se proponen en las expresiones anteriores se pueden obtener de forma directa mediante aplicación de reglas de derivación analíticas a partir de la formulación del problema estructural. Un desarrollo más extenso y detallado de las técnicas de análisis de sensibilidad que se proponen se puede encontrar en [[#bib0155|[31,50]]].      

La variable adjunta se puede obtener resolviendo el sistema de ecuaciones propuesto en [[#eq0295|(41)]]. De acuerdo con este planteamiento, con este método será necesario resolver tantos sistemas de ecuaciones como restricciones en tensión se incorporen en el algoritmo de optimización. En la práctica, este número es considerablemente más reducido que el número de variables de diseño ya que en el algoritmo de optimización solo se incorporarán las restricciones activas en cada iteración. Además, a pesar de que la matriz de coeficientes de [[#eq0295|(41)]]  aparece traspuesta, el cálculo de estos sistemas de ecuaciones se puede implementar de forma muy eficiente debido a las propiedades de simetría de la matriz de rigidez de la estructura. De este modo, en la práctica, la matriz de coeficientes de [[#eq0295|(41)]]  es la matriz de rigidez del análisis estructural al igual que ocurre en el método de diferenciación directa analítica.      

Una vez conocido el valor de la variable adjunta '''''λ'''''<sub>''e''</sub>  el cálculo de las derivadas de las restricciones se puede llevar a cabo a partir de [[#eq0285|(39)]]. Por lo tanto, el cálculo de estas derivadas también se realiza en dos pasos. En un primer paso se calcula el valor de la variable adjunta para cada restricción '''''λ'''''<sub>''e''</sub>  resolviendo el sistema de ecuaciones propuesto en [[#eq0295|(41)]]  y, a continuación, en un segundo paso se calculan las derivadas de las restricciones de acuerdo con [[#eq0285|(39)]]. Este método permite calcular las derivadas de las restricciones de forma selectiva, seleccionando por ejemplo aquellas restricciones que se encuentren activas en cada iteración, y de forma independiente para cada restricción. Esto también supondrá una ventaja muy importante si se pretende implementar este cálculo en paralelo en ordenadores de altas prestaciones.      

===5.2. Análisis de sensibilidad direccional de primer y segundo orden de las restricciones===

El análisis de sensibilidad de primer orden de las restricciones permite calcular la dirección de modificación del diseño más adecuada utilizando el algoritmo de programación lineal secuencial propuesto. Una vez calculada esta dirección se necesitan conocer las derivadas direccionales de primer y segundo orden de las restricciones con la finalidad de calcular el factor de avance más adecuado de acuerdo con [[#eq0265|(35)]].      

Para calcular las derivadas direccionales se utilizará un algoritmo de diferenciación directa analítica. En este caso no se utiliza el método de la variable adjunta ya que requeriría almacenar y calcular todas las variables adjuntas de todas las restricciones. El algoritmo de diferenciación directa permite calcular las derivadas de todas las restricciones, no solo de las restricciones activas. Las derivadas direccionales de primer orden de las restricciones en la dirección '''''s'''''  se pueden calcular a partir de [[#eq0270|(36)]]  como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>D_s\quad g_e=\frac{dg_e}{d\boldsymbol{\rho }}\boldsymbol{s},</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

y sustituyendo (<math display="inline">\frac{dg_e}{d\boldsymbol{\rho }}</math>) de acuerdo con [[#eq0275|(37)]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>D_s\quad g_e={\frac{\partial {\tilde{g}}_e}{\partial {\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\overset{\mbox{ˆ}}{\sigma }}_{VM,e}}{d{\boldsymbol{\sigma }}_e}}_{\begin{array}{l}
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}\quad {\frac{d{\boldsymbol{\sigma }}_e}{d\boldsymbol{\alpha }}}_{\boldsymbol{\alpha }(\boldsymbol{\rho })}\quad D_s\quad \boldsymbol{\alpha }+</math><math>{\frac{\partial {\tilde{g}}_e}{\partial {\rho }_e}}_{\begin{array}{l}
{\overset{\mbox{ˆ}}{\sigma }}_{\mbox{vm},e}({\boldsymbol{\sigma }}_e)\\
{\boldsymbol{\sigma }}_e(\boldsymbol{\alpha })\\
\boldsymbol{\alpha }(\boldsymbol{\rho })
\end{array}}s_e</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

donde ''D''<sub>''s''</sub>  '''''α'''''  se obtiene como:

<span id='eq0310'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{K}\quad (D_s\quad \boldsymbol{\alpha })=</math><math>(D_s\quad \boldsymbol{f})-(D_s\quad \boldsymbol{K})\quad \boldsymbol{\alpha }.</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

De un modo similar se pueden obtener las derivadas de segundo orden direccionales. Los algoritmos se omiten en este trabajo debido a su extensión. Un desarrollo completo de los métodos de análisis de sensibilidad propuestos se puede encontrar en [[#bib0155|[31,50]]]. Como se puede apreciar en [[#eq0310|(44)]]  se requiere resolver un sistema de ecuaciones adicional. De nuevo, la matriz de coeficientes de este sistema de ecuaciones es la matriz de rigidez del modelo de cálculo de estructuras. Por lo tanto, resolviendo solamente dos sistemas de ecuaciones adicionales (uno para las derivadas direccionales de primer orden y otro para las derivadas direccionales de segundo orden) se pueden calcular las derivadas direccionales de todas las restricciones del problema para cada caso de carga.      

<span id='sec0110'></span>
==6. Cálculo en paralelo en optimización topológica de estructuras==

Los modelos de diseño óptimo que se han propuesto involucran un número elevado de variables de diseño. En la práctica es fácil alcanzar varios miles de variables de diseño. Por otra parte, si se utiliza la formulación local el número de restricciones del modelo será igual o superior (si se aplican varios casos de carga) al número de variables de diseño. Por lo tanto, el problema de optimización resultante es muy complejo y en general el tiempo de cálculo necesario para resolverlo es muy elevado. Asimismo el análisis de sensibilidad de las restricciones también requiere un coste computacional elevado. Por este motivo se plantea incorporar las técnicas de cálculo en paralelo en el proceso de diseño óptimo.

Estas técnicas de cálculo en paralelo se aplican en este trabajo en los dos algoritmos que requieren la mayor parte del coste computacional: el análisis de sensibilidad de primer orden de las restricciones y el algoritmo de optimización basado en programación lineal secuencial.

Tal y como se ha comentado anteriormente en este trabajo el análisis de sensibilidad de primer orden se calcula mediante el método de la variable adjunta. Este método se desarrolla en dos fases: el cálculo de la variable adjunta y el cálculo de las derivadas a partir de ésta. Ambas fases se pueden desarrollar de forma independiente para cada restricción. Por lo tanto, la aplicación de técnicas de cálculo en paralelo en este caso es directa y los rendimientos obtenidos son, en la práctica, muy cercanos a los teóricos. En la [[#fig0030|figura 6]]  se muestran los resultados de aceleración de los cálculos del análisis de sensibilidad de primer orden obtenidos en un equipo de cálculo de altas prestaciones.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr6.jpg|center|392px|Aceleración obtenida por iteración en un problema con 7.200 variables de diseño ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 6'''. Aceleración obtenida por iteración en un problema con 7.200 variables de diseño y 7.200 restricciones locales en tensión, de las cuales 1.940 están activas. Equipo de cálculo formado por 4 procesadores Intel Xeon 7.120 M 3.0 GHz y 16 Gb de memoria.

</span>
|}

El algoritmo de programación lineal secuencial que se propone también se desarrolla en dos fases. En la primera fase se calcula la dirección más adecuada de modificación del diseño y en una segunda fase se calcula el factor de avance. El cálculo de la dirección de avance es la parte que requiere la mayor parte del tiempo CPU en la práctica. El tiempo de cálculo requerido para obtener el factor de avance es despreciable frente al cálculo de la dirección de modificación del diseño.

El cálculo de la dirección de modificación del diseño se lleva a cabo habitualmente mediante el algoritmo Simplex [[#bib0240|[48]]]  and [[#bib0245|[49]]]. Este algoritmo basa su funcionamiento en la aplicación iterativa de un algoritmo que modifica la matriz del problema linealizado. Esta transformación de la matriz del problema linealizado consiste en la aplicación de combinaciones lineales por filas de la misma [[#bib0155|[31]]]. Por lo tanto, estas combinaciones lineales por filas también se pueden llevar a cabo de forma independiente para cada fila, lo que permite desarrollar este cálculo en paralelo de forma sencilla. La aceleración obtenida en este caso como resultado de la utilización del cálculo en paralelo no es tan eficiente como para el cálculo del análisis de sensibilidad ya que es necesario desarrollar una parte del cálculo de forma secuencial [[#bib0155|[31]]]. En la [[#fig0030|figura 6]]  se muestran los resultados de aceleración de los cálculos obtenidos para este algoritmo en un equipo de cálculo de altas prestaciones.      

En conjunto la aceleración obtenida en cada iteración calculando en paralelo tanto el análisis de sensibilidad como el algoritmo Simplex se puede observar en la [[#fig0030|figura 6]] . Estos cálculos se han llevado a cabo utilizando directivas de paralelización OpenMP [[#bib0255|[51]]]  debido a su facilidad de implementación en un código secuencial ya existente y a su sencilla aplicación en equipos de cálculo de múltiples núcleos, que son los más habituales en la actualidad.      

<span id='sec0115'></span>
==7. Ejemplos de aplicación==

Como ejemplos de aplicación de la formulación de diseño óptimo así como de la metodología de optimización se analizan tres estructuras bidimensionales bajo la hipótesis de tensión plana. Se presentan tanto la evolución del proceso de optimización desde la solución inicial hasta la solución final como el estado tensional existente en la solución óptima final con el fin de verificar la validez de la formulación. También se tratan otros aspectos computacionales relacionados con el proceso de cálculo como por ejemplo, el tiempo de CPU total.

El modelo de elementos finitos propuesto para estos ejemplos utiliza elementos cuadráticos de 8 nodos con la finalidad de mejorar la aproximación tanto de los movimientos de la estructura como de las tensiones. Además, esta característica no supone un incremento de coste computacional relevante debido a que no se modifica ni el número de variables de diseño ni el número de restricciones del problema.

===7.1. Optimización topológica de una viga biapoyada de gran canto con carga inferior===

El primer ejemplo estudiado corresponde al problema de una viga biapoyada de gran canto sobre la que se aplica una carga vertical concentrada en la parte inferior de la misma. Un esquema de la geometría y dimensiones del dominio de estudio puede observarse en la [[#fig0035|figura 7]]. El peso propio de la estructura también se ha considerado en este ejemplo como una carga estructural. Por simetría solamente se analiza la mitad derecha del dominio, aunque se muestran las soluciones sobre el dominio completo.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr7.jpg|center|287px|Esquema y dimensiones de la viga biapoyada de gran canto con carga inferior. ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 7'''. Esquema y dimensiones de la viga biapoyada de gran canto con carga inferior. Unidades en metros.

</span>
|}

El dominio de la viga se ha discretizado usando una malla de 60 × 60 = 3.600 elementos cuadráticos de 8 nodos. El espesor de la viga analizada es de 1 cm.

La carga vertical se aplica sobre la cara inferior de los cuatro elementos del borde inferior más próximos al eje de simetría con la finalidad de evitar los efectos de concentración de tensiones. La carga vertical aplicada sobre la mitad derecha asciende a 50 ''k'' N, lo que supone un total de 100 ''k'' N si se considera el dominio completo. Esta carga se aplica de forma ligeramente distribuida ya que si se aplicasen cargas estrictamente puntuales los fenómenos de concentración de tensiones desvirtuarían los resultados del análisis estructural. Los fenómenos de concentración de tensiones conducirían a la no existencia de solución para este problema dado que no existiría ningún diseño que cumpliese las restricciones en tensión en las zonas afectadas. Además, este planteamiento no sería realista dado que en la práctica tanto las zonas de apoyo como las zonas de aplicación de cargas no son puntuales sino distribuidas.      

El material que se utiliza para este problema es acero con densidad ''γ''<sub>''mat''</sub>  = 7.850 kg/m<sup>3</sup>, módulo de Young ''E''  = 2, 1 × 10<sup>5</sup>  MPa, coeficiente de Poisson ''ν''  = 0, 3 y límite elástico <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}=</math><math>230\quad MPa</math>.      

En la [[#fig0040|figura 8]]  se puede observar la evolución de la solución para el problema propuesto desde la solución inicial. Para ello se muestran 6 soluciones intermedias correspondientes a números de iteración diferentes del proceso de optimización. Finalmente, en la [[#fig0045|figura 9]]  se puede observar la solución óptima alcanzada. La solución óptima obtenida corresponde con la solución esperada desde un punto de vista teórico. Para resolver este problema se ha utilizado la formulación local de las restricciones.

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr8.jpg|center|401px|Evolución de la solución hacia el óptimo para diferente número de iteraciones en ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 8'''. Evolución de la solución hacia el óptimo para diferente número de iteraciones en el problema de la viga biapoyada de gran canto.

</span>
|}

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr9.jpg|center|214px|Solución óptima obtenida para el problema de la viga biapoyada de gran canto.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 9'''. Solución óptima obtenida para el problema de la viga biapoyada de gran canto.

</span>
|}

El nivel de tensiones que presenta la solución final se muestra en la [[#fig0050|figura 10]]. Se muestran tensiones normalizadas por simplicidad. Como se puede observar el nivel de tensión de todos los elementos es correcto y no se superan los valores máximos permitidos para la tensión. Ello demuestra la validez de la formulación local de las restricciones para este problema.

<span id='fig0050'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr10.jpg|center|214px|Estado de tensiones normalizado para el problema de la viga biapoyada de gran ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 10'''. Estado de tensiones normalizado para el problema de la viga biapoyada de gran canto.

</span>
|}

En la tabla [[#tbl0005|1]]  se pueden observar y analizar los parámetros más importantes de este problema.

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1. Resumen de los parámetros más importantes del ejemplo de la viga biapoyada de gran canto

|-

! VIGA BIAPOYADA
! ''ɛ''
! ''p''
! ''q''
! Número de iteraciones
! Tiempo de CPU (s)
! Peso final (%)
|-

| Formulación local ([[#fig0045|fig. 9]])                                                    
| 0,01
| 4
| 0
| 564
| 202.360
| 6,18%
|}

===7.2. Optimización topológica de un gancho de grúa===

El segundo ejemplo que se propone corresponde al diseño y optimización de un gancho de grúa. Un esquema de la geometría y dimensiones del dominio de estudio puede observarse en la [[#fig0055|figura 11]]. El peso propio de la estructura también se ha considerado en este ejemplo como una carga estructural.

<span id='fig0055'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr11.jpg|center|233px|Esquema y dimensiones del gancho de grúa (unidades en centímetros).]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 11'''. Esquema y dimensiones del gancho de grúa (unidades en centímetros).

</span>
|}

El dominio del gancho se ha discretizado mediante 6.260 elementos cuadráticos de 8 nodos de acuerdo con la [[#fig0060|figura 12]]. El espesor del gancho es de 2 cm.

<span id='fig0060'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr12.jpg|center|232px|Malla de elementos finitos del gancho de grúa.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 12'''. Malla de elementos finitos del gancho de grúa.

</span>
|}

La carga vertical se aplica sobre el sector circular inferior abarcando un ángulo de 90<sup>''o''</sup> . Con el fin de simular el efecto las cargas reales no se ha utilizado una carga uniformemente distribuida. La carga aplicada sigue una ley senoidal que aplica un valor máximo en el eje vertical del arco de circunferencia sobre la que se aplican las cargas y un valor nulo en los extremos de aplicación de la carga, es decir a 45<sup>''o''</sup>  a cada lado del eje vertical del arco de circunferencia. La carga total aplicada es de 30 ''k'' N.      

El material utilizado en este problema es acero con densidad ''γ''<sub>''mat''</sub>  = 7.850 kg/m<sup>3</sup> , módulo de Young ''E''  = 2, 1 × 10<sup>5</sup>  MPa, coeficiente de Poisson ''ν''  = 0, 3 y límite elástico <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}=</math><math>230\quad MPa</math>.      

En la [[#fig0065|figura 13]]  se muestra la evolución de la solución a lo largo del proceso de optimización. Se muestran 6 soluciones intermedias que corresponden a diferente número de iteraciones del algoritmo de optimización desarrollado.

<span id='fig0065'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr13.jpg|center|395px|Evolución de la solución hacia el óptimo para diferente número de iteraciones en ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 13'''. Evolución de la solución hacia el óptimo para diferente número de iteraciones en el problema del gancho de grúa.

</span>
|}

La solución final obtenida se muestra en la [[#fig0070|figura 14]]. A juzgar por la distribución de material alcanzada, las soluciones que se utilizan en la práctica difieren sustancialmente de la solución óptima. Sin embargo, es necesario tener en cuenta que el modelo propuesto no incorpora restricciones por pandeo y este hecho puede afectar al diseño final propuesto. Por otra parte, también se muestra el nivel tensional de los elementos de la malla para la solución final en la [[#fig0075|figura 15]]. Como se puede observar el nivel de tensión de los elementos se satisface y no se superan los valores máximos permitidos para la tensión. En este caso también se ha utilizado la formulación local de las restricciones con la finalidad de garantizar la validez resistente de diseño final alcanzado.

<span id='fig0070'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr14.jpg|center|178px|Diseño óptimo obtenido para el problema del gancho de grúa.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 14'''. Diseño óptimo obtenido para el problema del gancho de grúa.

</span>
|}

<span id='fig0075'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr15.jpg|center|178px|Estado de tensiones normalizado para el problema del gancho de grúa.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 15'''. Estado de tensiones normalizado para el problema del gancho de grúa.

</span>
|}

En la tabla [[#tbl0010|2]]  se pueden observar y analizar los parámetros más importantes de este problema. Es importante en este punto destacar los elevados tiempos de cálculo que se requieren cuando se utiliza la formulación local de las restricciones. Este ejemplo se ha calculado en paralelo mediante directivas OpenMP en un servidor con dos procesadores Intel Xeon X5355 a 2,66 GHz y con 16 Gb de memoria. Para analizar este ejemplo se han utilizado los 8 núcleos de cálculo disponibles.

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 2. Resumen de los parámetros más importantes del ejemplo del gancho de grúa

|-

! GANCHO DE GRÚA
! ''ɛ''
! ''p''
! ''q''
! Número de iteraciones
! Tiempo de CPU (s)
! Peso final (%)
|-

| Formulación local ([[#fig0070|fig. 14]])                                                    
| 0,01
| 4
| 0
| 679
| 61.600
| 24,70%
|}

===7.3. Optimización topológica de la estructura de un puente colgante===

Este ejemplo corresponde al problema de un puente de gran luz con unas condiciones de apoyo que simulan las condiciones que presentan los puentes colgantes. Un esquema de la geometría y dimensiones del dominio de estudio puede observarse en la [[#fig0080|figura 16]]. Se han prescrito los movimientos en una zona localizada de la parte inferior del dominio que correspondería a un punto de apoyo (una pila). Asimismo también se fija la zona del tablero en el borde vertical derecho y se imponen los desplazamientos prescritos horizontales correspondientes en la zona del tablero de la izquierda para simular la continuidad del tablero en el centro de vano. El peso propio es muy importante en este tipo de estructuras y también se ha considerado en este ejemplo como una carga estructural. Solamente se analiza la mitad derecha del dominio por simetría, aunque se muestran las soluciones sobre el dominio completo.

<span id='fig0080'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr16.jpg|center|378px|Esquema y dimensiones de puente colgante. Unidades en metros.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 16'''. Esquema y dimensiones de puente colgante. Unidades en metros.

</span>
|}

El dominio de la viga se ha discretizado usando una malla de 6.030 elementos cuadráticos de 8 nodos. El espesor de la estructura es de 1 m.

La carga vertical se aplica de forma distribuida a lo largo del tablero del puente, que se establece a la cota 30 m y se representa fijando como ''ρ''  = 1 la densidad en una fila de elementos a esa altura. La carga vertical total aplicada sobre la mitad derecha asciende a 60.000 ''k'' N, lo que supone un total de 120.000 ''k'' N si se considera el dominio completo.      

El material que se utiliza para este problema presenta una densidad de ''γ''<sub>''mat''</sub>  = 2.400 kg/m<sup>3</sup>, módulo de Young ''E''  = 3 × 10<sup>5</sup>  MPa, coeficiente de Poisson ''ν''  = 0, 3 y límite elástico <math display="inline">{\overset{\mbox{ˆ}}{\sigma }}_{m\acute{a}x}=</math><math>24</math>  MPa.      

En la [[#fig0085|figura 17]]  se muestra la evolución de la solución a lo largo del proceso de optimización. Se muestran 6 soluciones intermedias que corresponden a diferente número de iteraciones del algoritmo de optimización desarrollado.

<span id='fig0085'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr17.jpg|center|376px|Evolución de la solución hacia el óptimo para diferente número de iteraciones en ...]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 17'''. Evolución de la solución hacia el óptimo para diferente número de iteraciones en el problema del puente colgante.

</span>
|}

La solución final obtenida se muestra en la [[#fig0090|figura 18]]. Como se puede observar el diseño final alcanzado no es el que se utiliza habitualmente en el diseño de puentes atirantados en Ingeniería. Este tipo de estructuras se definen en la práctica con un único mástil vertical (o pilona) que se establece como una prolongación de la pila que aparece desde el apoyo intermedio hasta el tablero del puente. Sin embargo, el diseño que se obtiene con la formulación propuesta define dos mástiles (o pilonas) colocadas de radialmente en forma de V desde la intersección del tablero con la pila situada sobre el apoyo intermedio. Este diseño es más eficiente, de acuerdo con la formulación propuesta, que el diseño con un único mástil vertical; si bien el modelo de diseño óptimo propuesto no incorpora restricciones por pandeo que podrían influir en la solución estructural alcanzada en este ejemplo. Por otra parte, también es importante destacar que para la longitud de vano que se analiza la solución de puente atirantado aparece como la solución estructural más adecuada frente a la disposición de una catenaria con péndolas verticales.

<span id='fig0090'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr18.jpg|center|376px|Diseño óptimo obtenido para el problema del puente colgante.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 18'''. Diseño óptimo obtenido para el problema del puente colgante.

</span>
|}

Por otra parte, también se muestra el nivel tensional de los elementos de la malla para la solución final en la [[#fig0095|figura 19]]. Como se puede observar, el nivel de tensión en los elementos se satisface si bien en algunos casos se superan ligeramente los valores máximos permitidos para la tensión dado que en este caso se utiliza una formulación por bloques. El número de bloques utilizado en este problema es de 603, lo que supone que cada bloque contiene 10 elementos de la malla.

<span id='fig0095'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_552266557-1-s2.0-S0213131511000368-gr19.jpg|center|376px|Estado de tensiones normalizado para el problema del puente colgante.]]


|-
| <span style="text-align: center; font-size: 75%;">

'''Figura 19'''. Estado de tensiones normalizado para el problema del puente colgante.

</span>
|}

En la [[#tbl0015|tabla 3]]  se pueden observar y analizar los parámetros más importantes de este problema. Es importante en este punto destacar la reducción en tiempo de cálculo obtenida con la formulación por bloques, incluso a pesar de que con esta formulación se necesitan en la práctica un mayor número de iteraciones del algoritmo de optimización. El valor del parámetro de agregación de las funciones globales utilizadas en la formulación por bloques también se puede observar en esta tabla.

<span id='tbl0015'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 3. Resumen de los parámetros más importantes del problema del puente colgante.

|-

! PUENTE COLGANTE
! ''ɛ''
! ''p''
! ''μ''
! Número de iteraciones
! Tiempo de CPU (s)
! Peso final (%)
|-

| Formulación por bloques ([[#fig0090|fig. 18]])                                                    
| 0,01
| 4
| 40
| 3.139
| 199.300
| 8,74%
|}

<span id='sec0135'></span>
==8. Conclusiones==

En este trabajo se propone una formulación de mínimo peso con restricciones en tensión para el problema de la optimización topológica de estructuras continuas. La formulación que se propone ha demostrado ser robusta y eficiente incluso cuando se resuelven modelos complejos con un elevado número de variables de diseño y de restricciones en tensión. La formulación propuesta ofrece ventajas importantes frente a las formulaciones habituales de máxima rigidez ya que se garantiza la validez de la solución y se evitan inestabilidades como disposiciones de material en damero. La metodología de diseño óptimo que se propone ha sido especialmente desarrollada e implementada para este problema teniendo en cuenta aspectos fundamentales como el tiempo de cálculo o la eficiencia de las operaciones. Así se han desarrollado algoritmos específicos para el cálculo de estructuras, el análisis de sensibilidad y la optimización de la solución. El uso de estas técnicas en los ejemplos de aplicación práctica ha demostrado la validez de la formulación y de la metodología de resolución propuesta.

==Agradecimientos==

Este trabajo ha sido parcialmente financiado por el Ministerio de Ciencia e Innovación (Proyectos DPI2009-14546-C02-01 y #DPI2010-16496), la ''Subdirección Xeral de I+D de la Xunta de Galicia''  (Proyectos PGDIT09MDS00718PR y PGDIT09REM005118PR) cofinanciado con fondos FEDER, la Universidade da Coruña y la Fundación de la Ingeniería Civil de Galicia.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] A.G.M. Michell; The limit of economy of material in frame structures; Philosophical Magazine, 8 (6) (1904), pp. 589–597</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] L.A. Schmit; Structural design by systematic synthesis; Proceedings of the Second ASCE Conference on Electronic Computation, Pittsburgh, Pensilvania, EE. UU. (1960), pp. 105–122</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] M.P. Bendsøe, N. Kikuchi; Generating optimal topologies in structural design using a homogenization method; Computer Methods in Applied Mechanics and Engineering, 71 (1988), pp. 197–224</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] M.P. Bendsøe; Optimal shape design as a material distribution problem; Structural Optimization, 1 (1989), pp. 193–202</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] M.P. Bendsøe; Optimization of Structural Topology, Shape, and Material; Springer-Verlag (1995)</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] M.P. Bendsøe; Variable-topology optimization: status and challenges; W. Wunderlich (Ed.), Proceedings of the European Conference on Computational Mechanics, ECCM’99, TUM, Munich (1999)</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] M.P. Bendsøe, O. Sigmund; Material interpolation schemes in topology optimization; Archive of Applied Mechanics, 69 (1999), pp. 635–654</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] P. Duysinx, Optimisation Topologique: du Milieu Continu a la Structure Elastique, PhD Thesis, Université de Liegé, Liegé, Bélgica, (1996).</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] H.A. Eschenauer, A. Schumacher; Simultaneous Shape and Topology Optimization of Structures; N. Olhoff, G.I.N. Rozvany (Eds.), Proc. on First World Congress of Structural and Multidisciplinary Optimization (WCSMO1), Pergamon (1995), pp. 177–184</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] H.A. Eschenauer, N. Olhoff; Topology optimization of continuum structures: A review; Applied Mechanics Review, 54 (2001), pp. 331–389</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] K.S. Park, S.K. Youn; Globally Convergent Topology Optimization unisng Level Set Method; Proceedings of the VI World Congress of Structural and Multidisciplinary Optimization (WCSMO6), Rio de Janeiro, Brasil (2005)</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] K. Suzuki, N. Kikuchi; A Homogeneization Method for Shape and Topology Optimization; Computer Methods in Applied Mechanics and Engineering, 93 (3) (1991), pp. 291–318</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] S.K. Youn, S.H. Park; A study on the shape extraction process in the structural topology optimization using homogenized material; Computers & Structures, 62 (3) (1997), pp. 527–538</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] A. Díaz, O. Sigmund; Checkerboard patterns in layout optimization; Structural Optimization, 10 (1995), pp. 40–45</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] G. Allaire, F. Jouve, H. Maillot; Topology optimization for minimum stress design with the homogenization method; Structural and Multidisciplinary Optimization, 28 (2004), pp. 87–98</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] M. Burger, R. Stainko; Phase-Field Relaxation of Topology Optimization with Local Stress Constraints; IUTAM Symposium on Topology Design Optimization of Structures, Machines and Materials, Rungsted, Dinamarca (2005)</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] I. Colominas, F. Navarrina, M. Casteleiro; Minimum weight with stress constraints topology optimization; Proceedings de la VIII International Conference on Computer Aided Optimum Design of Structures, Detroit, USA (2003)</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] G.D. Cheng, Z. Jiang; Study on topology optimization with stress constraints; Engineering Optimization, 20 (1992), pp. 129–148</li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] G.D. Cheng, X. Guo; ''ɛ'' -relaxed approach in structural topology optimization                                        ; Structural Optimization, 13 (1997), pp. 258–266</li>
<li><span id='bib0100'></span>
[[#bib0100|[20]]] P. Duysinx; Topology optimization with different stress limits in tension and compression; International report: Robotics and Automation, Institute of Mechanics, University of Liege, Lieja, Bélgica (1998)</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] P. Duysinx, M.P. Bendsøe; Topology optimization of continuum structures with local stress constraints; International Journal of Numerical Methods in Engineering, 43 (1998), pp. 1453–1478</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] P. Duysinx, O. Sigmund; New Developments in Handling Stress Constraints in Optimal Material Distributions; 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Design Optimization, Saint Louis, Missouri, EE. UU. (1998)</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] Q.Q. Liang, Y.M. Xie, G.P. Steven; Optimal selection of topologies for the minimum-weight design of continuum structures with stress constraints; Journal of Mechanical Engineering Science, Proceedings of The Institution of Mechanical Engineers, Part C, 213 (8) (1999), pp. 755–762</li>
<li><span id='bib0120'></span>
[[#bib0120|[24]]] F. Navarrina, I. Muíños, I. Colominas, M. Casteleiro; Topology Optimization of structures: A minimum weight approach with stress constraints; Advances in Engineering Software, 36 (2005), pp. 599–606</li>
<li><span id='bib0125'></span>
[[#bib0125|[25]]] J.T. Pereira, C.S. Barcellos, E.A. Fancello; Topology Optimization of Continuum structures with Material Failure Constraints; Structural and Multidisciplinary Optimization, 26 (1-2) (2004), pp. 50–66</li>
<li><span id='bib0130'></span>
[[#bib0130|[26]]] R. Stainko, M. Burguer; A one shot approach to topology optimization with local stress constraints; Proceedings of IUTAM-Symposium 2005 on Topology Design Optimization of Structures, Machines and Materials, publicado en Solid Mechanics and its Applications, 137 (2006), pp. 181–184</li>
<li><span id='bib0135'></span>
[[#bib0135|[27]]] M. Stolpe, T. Stidsen; A hierarchical method for discrete structural topology design problems with local stress and displacement constraints; International Journal for Numerical Methods in Engineering, 69 (5) (2007), pp. 1060–1084</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] Y.K. Sui, X.R. Peng, J.L. Feng, H.L. Ye; ICM Method of Topology Optimization for Continuum Structure with Globalization of Stress Constraints; Proceedings of the VI World Congress of Structural and Multidisciplinary Optimization (WCSMO6), Río de Janeiro, Brasil (2005)</li>
<li><span id='bib0145'></span>
[[#bib0145|[29]]] M. Werme, K. Svanberg; A hierarchical neighbourhood search method for stress-constrained topology optimization; Proceedings of the VI World Congress of Structural and Multidisciplinary Optimization (WCSMO6), Río de Janeiro, Brasil (2005)</li>
<li><span id='bib0150'></span>
[[#bib0150|[30]]] R.J. Yang, C.J. Chen; Stress-based topology optimization; Structural Optimization, 12 (1996), pp. 98–105</li>
<li><span id='bib0155'></span>
[[#bib0155|[31]]] J. París, Restricciones en tensión y minimización del peso: una metodología general para la optimización topológica de estructuras continuas, Tesis Doctoral, www.tdr.cesca.es , http://ruc.udc.es/dspace/handle/2183/1104 , Universidade da Coruña, España (2008).                                    </li>
<li><span id='bib0160'></span>
[[#bib0160|[32]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Topoloy optimization of continuum structures with local and global stress constraints; Structural and Multidisciplinary Optimization, 39 (4) (2009), pp. 419–437</li>
<li><span id='bib0165'></span>
[[#bib0165|[33]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Block aggregation of stress constraints in topology optimization of structures; Advances in Engineering Software, 41 (2010), pp. 433–441</li>
<li><span id='bib0170'></span>
[[#bib0170|[34]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Improvements in the treatment of stress constraints in structural topology optimization problems; Journal of Computational and Applied Mathematics, 234 (2010), pp. 2231–2238</li>
<li><span id='bib0175'></span>
[[#bib0175|[35]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Some Improvements in Minimum Weight Topology Optimization with Stress Constraints; Computer Aided Optimum Design of Structures XI, WIT Press, Southampton, Reino Unido (2009) pp. 71–82. ISBN: 978-1-84564-185-6</li>
<li><span id='bib0180'></span>
[[#bib0180|[36]]] Q.Q. Liang, G.P. Steven; A performance-based optimization method for topology design of continuum structures with mean compliance constraints; Computer Methods in applied Mechanics and Engineering, 191 (2002), pp. 1471–1489</li>
<li><span id='bib0185'></span>
[[#bib0185|[37]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Block Aggregation of Stress Constraints in Topology Optimization of Structures; Computer Aided Optimum Design of Structures X, WIT Press, Southampton, UK (2008) pp. 25-34. (ISBN: 978-1-84564-070-5)</li>
<li><span id='bib0190'></span>
[[#bib0190|[38]]] G. Kresisselmeier, R. Steinhauser; Systematic control design by optimizing a vector performance index; International Federation of Active Controls Symposium on Computer-Aided Design of Control systems, Zurich, Suiza (1979)</li>
<li><span id='bib0195'></span>
[[#bib0195|[39]]] J.R.R.A. Martins, N.M.K. Poon; On Structural Optimization Using Constraint Aggregation; Proceedings of the VI World Congress of Structural and Multidisciplinary Optimization (WCSMO6), Rio de Janeiro, Brasil (2005)</li>
<li><span id='bib0200'></span>
[[#bib0200|[40]]] H.P. Mlejnek, R. Schirrmacher; An engineers approach to optimal material distribution and shape finding; Computational Methods in Applied Mechanics and Engineering, 106 (1993), pp. 1–26</li>
<li><span id='bib0205'></span>
[[#bib0205|[41]]] T.J.R. Hughes; The finite element method: linear static and dynamic finite element analysis; Dover Publishers (2000)</li>
<li><span id='bib0210'></span>
[[#bib0210|[42]]] E. Oñate; Cálculo de Estructuras por el Método de Elementos Finitos: Análisis estático y lineal;  (Segunda edición)CIMNE, Barcelona (1995)</li>
<li><span id='bib0215'></span>
[[#bib0215|[43]]] O.C. Zienkiewickz, R.L. Taylor; El métodos de los elementos finitos, Vol. 1 y 2; CIMNE, Barcelona (2004)</li>
<li><span id='bib0220'></span>
[[#bib0220|[44]]] F. Navarrina, Una metodología general para optimización estructural en diseño asistido por ordenador, Tesis Doctoral, Universidad Politécnica de Cataluña, Barcelona, (1987).</li>
<li><span id='bib0225'></span>
[[#bib0225|[45]]] F. Navarrina, M. Casteleiro; A general methodologycal analysis for optimum design; International Journal of Numerical Methods in Engineering, 31 (1991), pp. 85–111</li>
<li><span id='bib0230'></span>
[[#bib0230|[46]]] F. Navarrina, R. Tarrech, I. Colominas, G. Mosqueira, J. Gómez-Calviño, M. Casteleiro; An efficient MP algorithm for structural shape optimization problems; S. Hernández, C.A. Brebbia (Eds.), Computer Aided Optimum Design of Structures VII, WIT Press, Southampton (2001), pp. 247–256</li>
<li><span id='bib0235'></span>
[[#bib0235|[47]]] T. Navarrina, R. Tarrech, I. Colominas, G. Mosqueira, J. Gómez-Calviño, M. Casteleiro; Efficient structural shape optimization: using directional high order sensitivity analysis to improve MP algorithms; Proceedings del 6th US National Congress on Computational Mechanics, Dearborn, USA (2001)</li>
<li><span id='bib0240'></span>
[[#bib0240|[48]]] G.B. Dantzig, M.N. Thapa; Linear Programming I:Introduction; Springer-Verlag, New York (1997)</li>
<li><span id='bib0245'></span>
[[#bib0245|[49]]] G.B. Dantzig, M.N. Thapa; Linear programming. II: Theory and extensions; Springer-Verlag, New York (2003)</li>
<li><span id='bib0250'></span>
[[#bib0250|[50]]] J. París, F. Navarrina, I. Colominas, M. Casteleiro; Stress constrains sensitivity Analysis in Structural Topology Optimization; Computer Methods in Applied Mechanics and Engineering, 199 (2010), pp. 2110–2122</li>
<li><span id='bib0255'></span>
[[#bib0255|[51]]] R. Chandra, R. Menon, L. Dagum, D. Kohr, D. Maydan, J. McDonald; Parallel Programming in OpenMP; Elsevier Science, Morgan Kauffmann Publishers (2000)</li>
</ol>