Shape Optimization of 2D parametric models based on geometric weighting and evolutionary rules

==Resumen==

En este trabajo se presenta un método para la optimización de forma de modelos bidimensionales sometidos a casos de carga simple o múltiple. La optimización se realiza aplicando iterativamente reglas evolutivas, basadas en el nivel de esfuerzos. Dichas reglas determinan zonas de la frontera en las que el material es subutilizado o sobreutilizado con el objetivo de que el modelo evolucione a un modelo más ligero y con un nivel de esfuerzos alto y homogéneo. La evolución se lleva a cabo modificando la frontera del modelo lentamente. Dado que la frontera se define con curvas paramétricas B-spline los cambios generados resultan en una frontera suave. La principal propuesta del método es que la modificación de la frontera está dada por un conjunto de desplazamientos que toman en cuenta, en su magnitud y dirección, la información geométrica de la vecindad, a diferencia de otros métodos que simplemente generan desplazamientos perpendiculares a la frontera. Al final del trabajo se presenta el diseño de un marco de bicicleta para mostrar el buen desempeño del método.

==Abstract==

This paper presents a method for shape optimization of two-dimensional models subjected to simple or multiple load cases. The optimization is performed iteratively using evolutionary rules, based on the stress level. These rules determine regions on the boundary where the material is underused or overused; the objective is that the model evolves to a minimum weight model with a high and homogeneous stress level. The evolution is performed by modifying the boundary of the model slowly. Since the boundary is defined by parametric B-spline curves, generated changes result in a smooth boundary. The main proposal of the method is that boundary modifications are given by a set of displacements that, its magnitude and direction, take into account the geometrical information of the neighborhood, unlike other methods that simply generate displacement perpendicular to the boundary. Finally, the paper presents the design of a frame bike in order to show the good performance of the method.

==Palabras clave==

Diseño mecánico ; Diseño óptimo ; Optimización estructural ; Optimización de forma ; Métodos evolutivos

==Keywords==

Mechanical design ; Optimal design ; Structural optimization ; Shape optimization ; Evolutionary methods

==1. Introducción==

Tradicionalmente, la calidad de la solución a un problema de diseño depende de la experiencia diseñador. El tiempo de diseño y sus costos son una función de la experiencia del diseñador y su equipo de trabajo, así como de su conocimiento de metodologías efectivas de diseño. En este sentido puede pensarse que, entre mayor sea la experiencia mejor será el diseño; sin embargo, es imposible que un diseñador pueda proponer y analizar todas las posibles soluciones a un problema de diseño pues, en general, estos son problemas con múltiples soluciones. Por otra parte, el número de variables de diseño que pueden ser consideradas por una persona y su equipo es limitado. Lo anterior trae como consecuencia que es casi imposible obtener de manera tradicional la solución óptima a un problema de diseño, entendida como la mejor solución posible al problema.

Los métodos de optimización se caracterizan por su capacidad de abordar problemas con múltiples: variables, escenarios y soluciones. Dadas estas características, estos métodos son ideales para la solución de problemas de diseño. El primero en aplicar optimización en el diseño de estructuras fue Maxwell [[#bib0005|[1]]]  a finales del siglo xix , seguido por Michell [[#bib0010|[2]]]  con su Optimal Layout Theory. Con la segunda mitad del siglo xx  comenzó el desarrollo de métodos de búsqueda basados en Inteligencia Artificial (AI por sus siglas en inglés), métodos poderosos de análisis como el Método de los Elementos Finitos (FEM) y el Método de los Elementos Frontera (BEM), y de equipos de cómputo cada vez más veloces y con mayor capacidad. Esta conjunción de desarrollos ha permitido en los últimos años la aparición de poderosos métodos de optimización aplicados en o creados para diferentes áreas de la ingeniería.      

Una de estas áreas es la ingeniería de estructuras para la cual se han desarrollado varios métodos que se engloban con el nombre de métodos de optimización estructural. Los métodos más sobresalientes son aquellos que simulan la manera en la que los seres se adaptan a su medio en la naturaleza, esto es, métodos que lenta e iterativamente cambian las características de una estructura de tal forma que ésta mejora su desempeño en cada iteración. En este sentido, se puede decir que dichos métodos emplean reglas evolutivas para tomar decisiones durante el proceso de diseño lo cual, en cierta medida, sustituye la experiencia de los diseñadores

Cada uno de los métodos emplea reglas o estrategias y criterios similares para provocar que las estructuras evolucionen hacia un óptimo. Lo que verdaderamente distingue a los métodos es la manera en que modifican a las estructuras. Por ejemplo, en el método Evolutionary Structural Optimization o ESO (Xie [[#bib0015|[3]]] , Steven [[#bib0020|[4]]] , Querin [[#bib0025|[5]]] ) la evolución se logra quitando y añadiendo segmentos de material en zonas sub y sobre utilizadas, respectivamente. Métodos que emplean motores de búsqueda inteligentes como Algoritmos Genéticos (GA por sus siglas en inglés) y que tiene muchos exponentes (Velázquez [[#bib0030|[6]]]  y Annicchiarico [[#bib0035|[7]]]  solo por mencionar un par) combinan las características de soluciones diversas a lo largo de muchas iteraciones hasta lograr la convergencia en el desempeño de la estructura. En estos no hay una adición o remoción de material como tal. No así en los desarrollados por Mattheck [[#bib0040|[8]]]  and [[#bib0045|[9]]] , Sethian [[#bib0050|[10]]] , Ramm [[#bib0055|[11]]] , Maute [[#bib0060|[12]]] , Cerver a[[#bib0065|[13]]] , entre otros, en los cuales se modifica la forma de una estructura con base en cambios en su frontera, definidos por la aplicación de reglas evolutivas en los puntos que la constituyen.      

En estos últimos casos, en los que propiamente se habla de optimización de forma, la frontera de una estructura cambia con base en su respuesta ante una situación dada de cargas. Dicha respuesta puede estar dada en términos de esfuerzo, deformación, energía de distorsión, características geométricas como curvatura, etc., calculados sobre los puntos de la frontera. En general, la optimización de la forma se genera aplicando, sobre dichos puntos, desplazamientos en dirección perpendicular a la frontera y cuya magnitud es calculada de manera diferente por cada método.

Aunque la manera descrita de llevar a cabo la optimización ha demostrado ser efectiva pues genera buenos resultados, puede ser mejorada. El desplazar perpendicularmente punto a punto la frontera, en general, pierde la información geométrica importante de los puntos vecinos. Es decir, la respuesta en un punto de una estructura depende en gran medida de las características de su vecindad, por lo que desplazar la frontera solo en dirección perpendicular no toma en cuenta dicha dependencia.

En este trabajo se presenta un método de optimización de forma que iterativamente modifica la forma de una estructura, modelada con curvas paramétricas, tomando en cuenta la respuesta de la estructura en los puntos de su frontera y su relación geométrica con su vecindad. Como se describe más adelante, la información de las relaciones geométricas se logra a través de las funciones de interpolación que definen las curvas paramétricas con las que se modela el contorno de la estructura. Dichas funciones se utilizan para calcular valores de ponderación que miden la influencia de la vecindad sobre un punto dado, para calcular la dirección que debe tener el desplazamiento aplicado durante cada iteración del proceso de optimización.

==2. Modelado paramétrico==

Una característica que es requerida en la forma de una estructura es que su contorno no presente cambios bruscos de dirección que puedan resultar, por ejemplo, en concentradores de esfuerzos. Por lo tanto, durante el proceso de optimización es necesario evitar que el contorno de la estructura forme zonas de curvatura elevada. El uso de curvas paramétricas adecuadas puede ayudar para tal situación, en particular aquellas diseñadas especialmente para comportarse «suavemente» ante los cambios.

En el diseño el uso de las curvas paramétricas, tales como las B-splines, es muy socorrido debido a que presentan propiedades deseables en dicha actividad (Shene [[#bib0070|[14]]] ):
* Intuitivas- los cambios realizados en ellas tienen una interpretación geométrica que puede intuirse
* Flexibles- la forma de crearlas y editarlas es muy sencilla y se tiene gran control sobre dichos cambios
* Manejo unificado- la creación de diferentes curvas no requiere técnicas distintas de manipulación
* Invariantes- la curva no cambia su geometría bajo transformaciones tales como traslación y rotación
* Eficiencia y estabilidad numérica- la ejecución de cálculos sobre la curva no provocan distorsión de su forma

Sean los puntos ''P''<sub>''1''</sub>'',P''<sub>''2''</sub>'',…,P''<sub>''m''</sub>  de control de un polígono y ''0=u''<sub>''1''</sub>''≤u''<sub>''2''</sub>''≤…≤ u''<sub>''m''</sub>''=1''  valores del parámetro ''u''  que define ''m-1''  segmentos. La curva B-spline ''C(u)'' , formada por ''m-1''  segmentos polinomiales de grado ''p''  cada uno, contenida por el polígono de control se define como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>C(u)=\sum_{i=1}^mN_{i,p}{\left(u\right)}^\ast P_i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

donde ''N''<sub>''i,p''</sub>''(u)''  son funciones de interpolación de grado ''p''  definidas por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{ll}
N_{i,0}(u)= & \begin{array}{cc}
1 & si\mbox{   }u_i\leq u\leq u_{i+1}\\
0 & \mbox{para}\mbox{   }\mbox{otro}\mbox{   }\mbox{caso}
\end{array}\\
N_{i,p}(u)= & \frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}(u)+\\
\mbox{   } & \frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u)
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

De las ecuaciones anteriores y de la [[#fig0005|figura 1]]  se destaca que el cambio en la forma de una curva B-spline es resultado del cambio de posición de los puntos de control y que dicho cambio genera cambios solo en una región bien definida de la curva.

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr1.jpg|center|367px|Curva B-spline y su polígono de control. a) curva original, b) el punto de ...]]


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

Figura 1.

Curva B-spline y su polígono de control. a) curva original, b) el punto de control 6 es desplazado lo cual genera que la forma de una región bien definida cambie; fuera de dicha región la curva permanece sin cambios.

</span>
|}

En cuanto a la optimización de forma, si el análisis en un punto indica que la forma de una estructura debe cambiar, entonces dicho cambio afecta una región bien definida de la curva B-spline que modela la frontera y, por lo tanto, dicho cambio no debe estar basado únicamente en la información del punto de análisis.

==3. Método==

===3.1. Regiones modificables y reglas evolutivas===

El método propuesto en este trabajo emplea reglas evolutivas para determinar aquellas regiones, con base en el análisis de los nodos, que deben ser modificadas. En particular, se presenta el método aplicado a problemas de estructuras de mínimo peso con nivel de esfuerzo homogéneo; sin embargo, el método no está restringido solo a este tipo de problemas. Las reglas evolutivas para determinar las regiones modificables son dos: la primera identifica aquellas regiones en las que la estructura está subutilizada, ecuación (3); la segunda identifica las regiones que son sobreutilizadas, ecuación (4).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\sigma }_i\leq TR^\ast {\sigma }_r</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\sigma }_i\geq TA^\ast {\sigma }_r\mbox{ }o\mbox{ }{\sigma }_i\geq {\sigma }_y</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

donde ''σ''<sub>''i''</sub>  es el esfuerzo de von Mises calculado en un punto ''i''  de la frontera, ''σ''<sub>''r''</sub>  es un esfuerzo de referencia (el máximo o el medio de von Mises en la estructura), ''σ''<sub>''y''</sub>  es el esfuerzo de fluencia del material del modelo, ''TR''  y ''TA''  son tasas de rechazo y adición de material respectivamente. ''TR''  y ''TA''  marcan los límites en el valor de esfuerzo inferior y superior respectivamente; esto es, regiones subutilizadas son aquellas en que el material es excesivo para la situación de esfuerzos generada y por ello debe retirarse material de ellas, mientras que regiones sobreutilizadas representan una amenaza para la integridad del modelo y deben ser reforzadas agregando material.      

La optimización de una estructura se lleva a cabo en un proceso iterativo, evaluando en cada iteración lo siguiente: si alguna de las reglas o ambas se cumplen, la estructura es modificada antes de iniciar una nueva iteración; si las reglas no son satisfechas entonces se dice que se ha llegado a un estado permanente, es decir se ha generado una estructura que no será modificada para los valores dados de ''TR''  y ''TA.''  En este punto ''TR''  y ''TA''  cambian de manera que la evolución continúe hacia una estructura con menor peso y un nivel de esfuerzos más alto y homogéneo. El cambio en los valores de ''TR''  y ''TA''  está dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>TR_{j+1}=TR_j+TER</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>TA_{j+1}=TA_j+TEA</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

donde ''j''  representa el ''j'' -ésimo estado permanente, ''TER''  y ''TEA''  son las tasas evolutivas de rechazo y adición respectivamente.      

En el caso de estructuras sometidas a N casos de carga, las ecuaciones (3) y (4) se modifican de la siguiente forma:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\sigma }_{i,k}\leq TR^\ast {\sigma }_{r,k};\mbox{   }\forall \mbox{   }k=</math><math>1,2,...,\mbox{   }\mbox{N}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\exists \mbox{   }k:\mbox{   }{\sigma }_{i,k}\geq TA^\ast {\sigma }_{r,k};\mbox{   }k=</math><math>1,2,...,\mbox{   }\mbox{N}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)
|}

donde ''σ''<sub>''i,k''</sub>  y ''σ''<sub>''r,k''</sub>  son los esfuerzos de von Mises en el punto ''i''  de la frontera para el ''k'' -ésimo caso de carga.      

La ecuaciones (7) y (8) representan reglas con lógica AND y OR respectivamente (Li [[#bib0075|[15]]] ); esto es, se remueve material siempre y cuando se cumpla la ecuación (7) para todos los casos de carga mientras que se añade material cuando al menos para un caso de carga se cumple (8).      

===3.2. Relación entre nodos y puntos de control===

La optimización de forma se lleva a cabo alterando la geometría de la frontera del modelo o estructura a optimizar, lo cual se logra desplazando los nodos que componen a la curva B-spline que la modela. El desplazamiento de los nodos es resultado de cambiar la posición de los puntos de control de la curva, de acuerdo con la ecuación (1). Como fue mencionado anteriormente, en este trabajo se propone un método para cambiar la forma de la frontera de un modelo empleando la información no solo de los nodos en particular sino también de su vecindad.

La posición de un nodo ''C(u''<sub>''i''</sub>'')''  de una curva B-spline depende de la posición de ''p+1''  puntos de control, es decir: ''P''<sub>''i-1''</sub>'', P''<sub>''i''</sub>'', P''<sub>''i+1''</sub>'',…,P''<sub>''i+p-1''</sub>  (Shene [[#bib0070|[14]]] ). En la [[#fig0010|figura 2]]  se muestra, para un polinomio de grado ''p=3'' , cuáles son los puntos de control que definen la posición de un punto ''C(u''<sub>''i''</sub>'')'' . Para la curva mostrada en la figura, la ecuación (1) adquiere la forma de la ecuación (9).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{cccccccccccccccccccccccccc}
C(u_1) & = & N_1(u_1)P_1 & + & N_2(u_1)P_2 & + & N_3(u_1)P_3 & + & 0 & + & 0 & + & 0 & + & 0 & + & \cdots  & + & 0 & + & 0 & + & 0 & + & 0 & ;\\
C(u_2) & = & N_1(u_2)P_1 & + & N_2(u_2)P_2 & + & N_3(u_2)P_3 & + & N_4(u_2)P_4 & + & 0 & + & 0 & + & 0 & + & \cdots  & + & 0 & + & 0 & + & 0 & + & 0 & ;\\
C(u_3) & = & 0 & + & N_2(u_3)P_2 & + & N_3(u_3)P_3 & + & N_4(u_3)P_4 & + & N_5(u_3)P_5 & + & 0 & + & 0 & + & \cdots  & + & 0 & + & 0 & + & 0 & + & 0 & ;\\
C(u_4) & = & 0 & + & 0 & + & N_3(u_4)P_3 & + & N_4(u_4)P_4 & + & N_5(u_4)P_5 & + & N_6(u_4)P_6 & + & 0 & + & \cdots  & + & 0 & + & 0 & + & 0 & + & 0 & ;\\
C(u_5) & = & 0 & + & 0 & + & 0 & + & N_4(u_5)P_4 & + & N_5(u_5)P_5 & + & N_6(u_5)P_6 & + & N_7(u_5)P_7 & + & \cdots  & + & 0 & + & 0 & + & 0 & + & 0 & ;\\
\vdots  & = & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & + & \vdots  & ;\\
C(u_{m-2}) & = & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & \cdots  & + & N_{m-3}(u_{m-2})P_{m-3} & + & N_{m-2}(u_{m-2})P_{m-2} & + & N_{m-1}(u_{m-2})P_{m-1} & + & N_m(u_{m-2})P_m & ;\\
C(u_{m-1}) & = & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & \cdots  & + & 0 & + & N_{m-2}(u_{m-1})P_{m-2} & + & N_{m-1}(u_{m-1})P_{m-1} & + & N_m(u_{m-1})P_m & ;\\
C(u_m) & = & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & 0 & + & \cdots  & + & 0 & + & 0 & + & N_{m-1}(u_m)P_{m-1} & + & N_m(u_m)P_m & ;
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)
|}

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr2.jpg|center|374px|La posición de un nodo sobre la curva paramétrica está definida por un conjunto ...]]


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

Figura 2.

La posición de un nodo sobre la curva paramétrica está definida por un conjunto de puntos de control.

</span>
|}

En forma matricial, la ecuación (9) es:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\boldsymbol{\mbox{C}}=\underset{\sim}{\boldsymbol{\mbox{N}}}\bullet \boldsymbol{\mbox{P}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

Por otra parte, las funciones de interpolación tienen la propiedad de unicidad, ecuación (11). la cual indica que los puntos de control ''P''<sub>''i-1''</sub>'',P''<sub>''i''</sub>'',P''<sub>''i+1''</sub>'',…,P''<sub>''i+p-1''</sub>  contribuyen en la posición del nodo ''C(u''<sub>''i''</sub>'')''  con porcentajes ''N''<sub>''i-1''</sub>''(u''<sub>''i''</sub>''),N''<sub>''i''</sub>''(u''<sub>''i''</sub>''),N''<sub>''i+1''</sub>''(u''<sub>''i''</sub>''),…,N''<sub>''i+p-1''</sub>''(u''<sub>''i''</sub>'')'' .

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>N_{i-1}(u_i)+N_i(u_i)+N_{i+1}(u_i)+...+N_{i+p-1}(u_i)=</math><math>1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

En este mismo sentido, un punto de control ''P''<sub>''i''</sub>  influye en la posición de los nodos ''C(u''<sub>''i-p+1''</sub>''),C(u''<sub>''i-p+2''</sub>''),…,C(u''<sub>''i''</sub>''),C(u''<sub>''i+1''</sub>'')''  siendo ''N''<sub>''i''</sub>''(u''<sub>''i-p+1''</sub>''),N''<sub>''i''</sub>''(u''<sub>''i-p+2''</sub>''),…,N''<sub>''i''</sub>''(u''<sub>''i''</sub>''),N''<sub>''i''</sub>''(u''<sub>''i+1''</sub>'')''  los valores de dicha influencia, respectivamente, tal como se muestra en la  [[#fig0015|figura 3]] . Es importante señalar que en este caso la propiedad de unicidad no se satisface, es decir:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>N_i(u_{i-p+1})+N_i(u_{i-p+2})+...+N_i(u_i)+N_i(u_{i+1})\not =</math><math>1</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr3.jpg|center|287px|La posición de un punto de control está definida por un conjunto de nodos sobre ...]]


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

Figura 3.

La posición de un punto de control está definida por un conjunto de nodos sobre la curva paramétrica.

</span>
|}

De las relaciones entre nodos, puntos de control y las ecuaciones (11) y (12) se concluyen dos cosas: en primer lugar, la posición de un nodo está dada por un conjunto de puntos de control los cuales contribuyen de manera porcentual a dicha posición; en segundo lugar, la posición de un punto de control está dada por un conjunto de nodos sobre una curva B-spline cuyas funciones de interpolación no cumplen con la propiedad de unicidad y, por lo tanto, no pueden considerarse como valores porcentuales o fraccionarios de influencia. Esto último es de gran importancia para las ideas que se presentan a continuación.

Si el nodo ''C(u''<sub>''i''</sub>'')  ''  debe cambiar su posición por un desplazamiento <math display="inline">{\overset{\rightarrow}{D}}_i</math> , entonces ''P''<sub>''i-1''</sub>'',P''<sub>''i''</sub>'',P''<sub>''i+1''</sub>'',…,P''<sub>''i+p-1''</sub>  cambiarán la suya por los desplazamientos <math display="inline">{\overset{\rightarrow}{d}}_{i-1},{\overset{\rightarrow}{d}}_i,{\overset{\rightarrow}{d}}_{i+1},\ldots ,{\overset{\rightarrow}{d}}_{i+p-1}</math> , respectivamente. Ahora bien, si dos nodos vecinos ''C(u''<sub>''i''</sub>'')''  y ''C(u''<sub>''i+1''</sub>'')  ''  deben cambiar su posición, respectivamente, por desplazamientos <math display="inline">{\overset{\rightarrow}{D}}_i</math>  y <math display="inline">{\overset{\rightarrow}{D}}_{i+1}</math>  independientes entre sí, entonces para los puntos de control relacionados con estos nodos deben calcularse los conjuntos de desplazamientos, también independientes entre sí, <math display="inline">{\overset{\rightarrow}{d}}_{i-1,i},{\overset{\rightarrow}{d}}_{i,i},{\overset{\rightarrow}{d}}_{i+1,i},\ldots ,{\overset{\rightarrow}{d}}_{i+p-1,i}</math>  y <math display="inline">{\overset{\rightarrow}{d}}_{i,i+1},{\overset{\rightarrow}{d}}_{i+1,i+1},\ldots ,{\overset{\rightarrow}{d}}_{i+p-1,i+1,}{\overset{\rightarrow}{d}}_{i+p,i+1}</math> , donde el primer subíndice corresponde al punto de control a desplazar y el segundo al nodo que provoca el desplazamiento.      

Se observa que, para que ambos desplazamientos nodales sean posibles, varios puntos de control deberían sufrir dos desplazamientos independientes lo cual, claramente, no produciría las posiciones deseadas de los nodos. La explicación de esto es simple: una de las características principales de las curvas B-spline, como se explicó anteriormente, es que al desplazar un punto de control se cambia la geometría de una región y no solo de un nodo. Como consecuencia, satisfacer desplazamientos nodales independientes y simultáneos no es posible.

En el caso general, es decir que todos los nodos de la curva deban cambiar de posición, los desplazamientos simultáneos que debería sufrir un punto de control ''P''<sub>''i''</sub>  son: <math display="inline">{\overset{\rightarrow}{d}}_{i,i-p+1,}{\overset{\rightarrow}{d}}_{i,i-p+2},\ldots ,{\overset{\rightarrow}{d}}_{i,i},{\overset{\rightarrow}{d}}_{i,i+1}</math> . Cada uno de estos desplazamientos, todos ellos independientes entre sí, indica que el punto de control debe cambiar su posición para que un nodo cambie parcialmente la suya.      

En este punto es importante señalar lo siguiente. A pesar de que la ecuación (10) aparentemente indica una manera de calcular los puntos de control para un conjunto de nodos dado, simplemente premultiplicando por la inversa de la matriz, en realidad no es lo mejor en el caso descrito. Lo anterior se debe a que, de acuerdo con la ecuación (2), las funciones de interpolación ''N''<sub>''i,p''</sub>''(u)''  dependen de relaciones entre valores discretos del parámetro (''u''<sub>''1''</sub>'', u''<sub>''2''</sub>'', …, u''<sub>''m)''</sub> , los cuales a su vez dependen de la geometría de la curva B-spline. Dada una curva, se pueden calcular los parámetros con base en los métodos: espaciado uniforme, longitud de cuerda, centrífugo, entre otros, obteniendo resultados muy diferentes con cada uno de ellos. En este trabajo se empleó el método Centrífugo debido a que es el que genera menores “rizos” en las curvas B-spline (Shene[[#bib0070|[14]]] ).      

Si la posición de cada nodo es modificada directamente, en lugar de hacerlo a través de los puntos de control, las relaciones entre los valores ''u''<sub>''1''</sub>'', u''<sub>''2''</sub>'', …, u''<sub>''m''</sub>  cambia y, por lo tanto, las funciones de interpolación también deberían cambiar, esto es, ser recalculadas. Como consecuencia, la matriz de funciones de interpolación cambiaría. Es por esto que en los sistemas CAD la forma de una curva B-spline, en la mayoría de los casos, es modificada desplazando los puntos de control.      

De lo anterior, resulta evidente que es imposible cumplir exactamente con el desplazamiento requerido sobre cada punto de control que, de manera independiente, cada nodo requiere. Dicho de otra manera, el optimizar la forma de un modelo requiere que su frontera sea modificada con base en el valor de una o varias características de la misma; por ejemplo, en el caso de este trabajo, el nivel de esfuerzos. Dichas características no solo dependen de la geometría del modelo, sino de muchas otras variables. Esta complejidad inherente provoca que los desplazamientos nodales necesarios en la optimización sean independientes entre sí y, por lo tanto, imposibles de satisfacer simultáneamente si se desea conservar las bondades de las curvas B-spline.

Para remediar lo anterior, en este trabajo se propone utilizar la información de la influencia del punto de control sobre la posición de los nodos, es decir, los valores de las funciones de interpolación, para determinar la prioridad de los desplazamientos. Así, las funciones de interpolación se convierten en valores de ponderación en el cálculo de los desplazamientos que un punto de control debería realizar en cada iteración del proceso de optimización. En la siguiente sección se describe detalladamente este punto.

===3.3. Cálculo de los desplazamientos ponderados===

El cálculo de los desplazamientos de los puntos de control se realiza en tres pasos:
* Determinar los desplazamientos nodales <math display="inline">{\overset{\rightarrow}{D}}_i</math>
* Calcular la dirección de los desplazamientos de los puntos de control
* Calcular la magnitud de los desplazamientos de los puntos de control.

''Paso 1.''  Los desplazamientos nodales se calculan para aquellos nodos que satisfacen las ecuaciones (3) y (4) en el caso de un solo caso de carga, o las ecuaciones (7) y (8) para múltiples casos de carga. Su magnitud está dada por el valor del esfuerzo calculado en el nodo y el valor del esfuerzo de referencia. En el caso de carga simple, la ecuación (13) corresponde a los nodos que cumplen con la ecuación (3), mientras que la (14) corresponde a los que satisfacen la (4).

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\|{\overset{\rightarrow}{D}}_i\right\|=Abs\left(\frac{{\sigma }_i-TR^\ast {\sigma }_r}{TR\ast {\sigma }_r}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\|{\overset{\rightarrow}{D}}_i\right\|=Abs\left(\frac{{\sigma }_i-TA^\ast {\sigma }_r}{TA\ast {\sigma }_r}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

Para múltiples casos de carga, a las ecuaciones (7) y (8) les corresponden las ecuaciones (15) y (16) respectivamente:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\|{\overset{\rightarrow}{D}}_i\right\|=Min\left(Abs\left(\frac{{\sigma }_{i,k}-TR^\ast {\sigma }_{r,k}}{TR^\ast {\sigma }_{r,k}}\right)\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\|{\overset{\rightarrow}{D}}_i\right\|=Max\left(Abs\left(\frac{{\sigma }_{i,k}-TA^\ast {\sigma }_{r,k}}{TA^\ast {\sigma }_{r,k}}\right)\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

donde ''k''  indica cada uno de los casos de carga.      

La dirección de los desplazamientos nodales está dada por el vector unitario perpendicular y externo de la frontera en el nodo (<math display="inline">{\overset{\mbox{ˆ}}{n}}_i</math> ). La dirección es:
* <math display="inline">-{\overset{\mbox{ˆ}}{n}}_i</math> , si las ecuaciones (3) o (7) se satisfacen                          
* <math display="inline">+{\overset{\mbox{ˆ}}{n}}_i</math> , si las ecuaciones (4) u (8) se cumplen                  

''Paso 2.''  La dirección de los desplazamientos de los puntos de control se basa en los desplazamientos nodales calculados y las funciones de interpolación. La posición de un punto de control ''P''<sub>''i''</sub>  cambia si los nodos ''C(u''<sub>''i-p+1''</sub>''),C(u''<sub>''i-p+2''</sub>''),…,C(u''<sub>''i''</sub>''),C(u''<sub>''i+1''</sub>'')  ''  son modificados por desplazamientos <math display="inline">{\overset{\rightarrow}{D}}_{i-p+1,}{\overset{\rightarrow}{D}}_{i-p+2},\ldots ,{\overset{\rightarrow}{D}}_i,\mbox{   }{\overset{\rightarrow}{D}}_{i+1}</math> . La relación que existe entre nodos y el punto de control está dada por los valores de las funciones de interpolación ''N''<sub>''i''</sub>''(u''<sub>''i-p+1''</sub>''),N''<sub>''i''</sub>''(u''<sub>''i-p+2''</sub>''),…,N''<sub>''i''</sub>''(u''<sub>''i''</sub>''),N''<sub>''i''</sub>''(u''<sub>''i+1''</sub>'')  '' , las cuales son empleadas como valores de ponderación para calcular la dirección del vector <math display="inline">{\overset{\rightarrow}{d}}_i</math> , [[#fig0020|figura 4]] , de acuerdo con la siguiente ecuación:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overset{\mbox{ˆ}}{d}}_i=\frac{\sum_{k=-p+1}^1\left(N_i{\left(u_{i+k}\right)}^\ast {\overset{\rightarrow}{D}}_{i+k}\right)}{\left\|\sum_{k=-p+1}^1\left(N_i{\left(u_{i+k}\right)}^\ast {\overset{\rightarrow}{D}}_{i+k}\right)\right\|}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)
|}

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr4.jpg|center|371px|La dirección del vector de desplazamiento necesario para un punto de control se ...]]


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

Figura 4.

La dirección del vector de desplazamiento necesario para un punto de control se calcula con base en los valores de las funciones de interpolación y los desplazamientos nodales de la frontera.

</span>
|}

''Paso 3.  ''  El cálculo de la magnitud del vector <math display="inline">{\overset{\rightarrow}{d}}_i</math>  depende de la magnitud del desplazamiento (<math display="inline">{\overset{\rightarrow}{D}}_q</math> ) que corresponde al nodo que tiene mayor influencia en la posición del punto de control ''P''<sub>''i''</sub> , es decir, ''q''  corresponde al subíndice del parámetro ''u''  de:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>N_i(u_q)=Max\left(N_i(u_{i-p+1}),N_i(u_{i-p+2}),...,N_i(u_i),N_i(u_{i+1})\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

De tal forma que la magnitud de <math display="inline">{\overset{\rightarrow}{d}}_i</math>  está dada por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left\|{\overset{\rightarrow}{d}}_i\right\|=\frac{\left\|{\overset{\rightarrow}{D}}_q\right\|}{Max\left(\left\|{\overset{\rightarrow}{D}}_1\right\|+\left\|{\overset{\rightarrow}{D}}_2\right\|+...+\left\|{\overset{\rightarrow}{D}}_m\right\|\right)}{\,}^\ast esize^\ast f</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

donde ''esize''  es el tamaño promedio de los elementos finitos en la malla del modelo y ''f''  es un factor de corrección para evitar que la frontera del modelo se distorsione debido a los desplazamientos acumulados (Cervera  [[#bib0065|[13]]] , Yang [[#bib0080|[16]]] ). Por lo tanto, combinando los resultados de (17) y (19), el vector de desplazamiento de cada punto de control, en una iteración dada del proceso de optimización, se calcula como:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\overset{\rightarrow}{d}}_i=\vert \vert {\overset{\rightarrow}{d}}_i\vert \vert ^\ast {\overset{\mbox{ˆ}}{d}}_i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

===3.4. Criterio de paro===

El criterio de paro empleado busca la convergencia en el nivel de esfuerzos del modelo (Cervera [[#bib0065|[13]]] ). El nivel de esfuerzos ''NS''  en una malla con ''ne''  elementos finitos está dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>NS=\sum_{i=1}^{ne}{\left({\sigma }_i-{\sigma }_r\right)}^2</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

Con base en este valor, el criterio de paro está dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Abs\left(\frac{NS_{t+1}-NS_t}{NS_t}\right)\leq {10}^{-5}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

donde ''t''  y ''t+1''  corresponden al número de las iteraciones consideradas en el cálculo.      

==4. Caso de estudio==

En este apartado se aplica el algoritmo de optimización propuesto en la optimización de un marco de bicicleta, modelado con un material con módulo de Young 200 GPa y coeficiente de Poisson de 0,3, sujeto a tres casos de carga definidos de acuerdo con el tipo de esfuerzo que realiza el ciclista (ver [[#tbl0005|tabla 1]] ). El peso considerado del ciclista fue 75 kg (∼ 750 N).

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

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

Tabla 1.

Casos de carga considerados en el diseño óptimo del marco de bicicleta (Tom Compton [[#bib0085|[17]]] )                  

|-

! Caso
! Descripción
! Cond. Avance
! Velocidad m/s
! Dirección de la carga
! Fuerza Asiento (N)
! Fuerza Manubrio (N)
! Fuerza pedal (N)
|-

| LC1
| Avance en una superficie plana de pendiente 0̊ a velocidad moderada
| Ligero
| 10
| x
| 0
| 0
| 0
|-

| 
| 
| 
| 
| y
| −800
| 0
| −104
|-

| LC2
| Avance con pendiente ligera (1/10) y velocidad moderada
| Moderado
| 8
| x
| 80
| 0
| 52
|-

| 
| 
| 
| 
| y
| −796
| 0
| −512
|-

| LC3
| Avance con pendiente grande (1/25) y baja velocidad. El ciclista no está sentado en el asiento
| Fuerte
| 4
| x
| 0
| 400
| 414
|-

| 
| 
| 
| 
| y
| 0
| 0
| −1.655
|}

El diseño inicial del marco se muestra en la [[#fig0025|figura 5]] , el cual está modelado con 10 curvas: L1, L3, L4, L5, L7, L9 y L10 son líneas fijas que pueden modificarse, mientras que L2, L6 y L8 son curvas B-spline de grado 3.

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr5.jpg|center|376px|Diseño inicial.]]


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

Figura 5.

Diseño inicial.

</span>
|}

Los valores de los parámetros de optimización empleados son TR = 0,2, TA = 0,95, TRE = 0,01 y TAE = 0,001. Los valores del análisis por elementos finitos calculados para el diseño inicial son los mostrados en la [[#tbl0010|tabla 2]] .

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

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

Tabla 2.

Esfuerzos calculados por medio de análisis por elementos finitos en el diseño incial

|-

! Caso de carga
! Esfuerzo máximo de von Mises (kPa)
! Esfuerzo promedio de von Mises (kPa)
|-

| LC1
| 623,7
| 37,6
|-

| LC2
| 658,7
| 53,7
|-

| LC3
| 600,5
| 86,9
|}

El proceso de optimización se llevó a cabo en 732 iteraciones durante las cuales se eliminó material lentamente del marco mientras la distribución de esfuerzos se homogeneizaba para emplear más eficientemente el material, tal como se muestra en la [[#fig0030|figura 6]] . El diseño final del marco de bicicleta generado por el proceso de optimización es el mostrado en la [[#fig0035|figura 7]] , el cual presenta los valores de esfuerzo de von Mises mostrados en la [[#tbl0015|tabla 3]] . En esta misma tabla se muestra una comparación entre los valores de los diseños inicial y final.

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr6.jpg|center|367px|Distribuciones de esfuerzos de von Mises de cada caso de carga para diferentes ...]]


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

Figura 6.

Distribuciones de esfuerzos de von Mises de cada caso de carga para diferentes iteraciones.

</span>
|}

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr7.jpg|center|301px|Marco de bicicleta óptimo generado por el método propuesto.]]


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

Figura 7.

Marco de bicicleta óptimo generado por el método propuesto.

</span>
|}

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

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

Tabla 3.

Esfuerzos calculados por medio de análisis por elementos finitos en el diseño final (óptimo)

|-

! Caso de carga
! Esfuerzo máximo de von Mises (kPa)
! Incremento respecto del inicial (%)
! Esfuerzo promedio de von Mises (kPa)
! Incremento respecto del inicial (%)
|-

| LC1
| 627,7
| 0,6
| 188,9
| 402,4
|-

| LC2
| 1183,4
| 79,7
| 274,2
| 410,6
|-

| LC3
| 1.602,9
| 166,9
| 442,6
| 409,3
|}

El marco obtenido muestra una mejor utilización del material lo cual se refleja en los valores y porcentajes comparativos mostrados en la [[#tbl0015|tabla 3]] , sobre todo en cuanto a los datos del esfuerzo promedio. Por otra parte, este diseño es similar al generado por Mevdev [[#bib0090|[18]]]  ([[#fig0040|fig. 8]] ) para un marco de bicicleta empleado para una competencia de velocidad.

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

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


[[Image:draft_Content_868989460-1-s2.0-S0213131512000235-gr8.jpg|center|301px|Bicicleta comercial cuyo diseño fue generado por Mevdev.]]


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

Figura 8.

Bicicleta comercial cuyo diseño fue generado por Mevdev.

</span>
|}

==5. Conclusiones==

Existe un número considerable de métodos de optimización de forma en la literatura. Sin embargo, la manera de modificar la geometría de los modelos o estructuras optimizadas de estos métodos no considera la información completa de la región que se está optimizando, sobre todo aquella relacionada con la geometría de la vecindad y su relación con los puntos que se modifican. El método presentado en este trabajo solventa esta deficiencia empleando las características de las curvas paramétricas, en particular las B-spline, las cuales están dadas por las funciones de interpolación que les dan forma. El método así concebido demostró ser capaz de optimizar una estructura sometida a múltiples casos de carga, al grado que el diseño obtenido es similar a uno existente comercialmente.

==Agradecimientos==

A la Universidad Nacional Autónoma de México por el apoyo proporcionado para la realización de este trabajo a través del proyecto IN108909 «Herramientas Computacionales para el Diseño Óptimo en Proyectos de Ingeniería e Innovación Tecnológica» del Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (PAPIIT) de la Dirección General de Asuntos de Personal Académico (DGAPA).

==Bibliografía==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[1]]] J.C. Maxwell; On Reciprocal Figures, Frames and Diagrams of Force; Trans Roy Soc Edinburgo, 26 (1872)</li>
<li><span id='bib0010'></span>
[[#bib0010|[2]]] A.G.M. Michell; The limits of economy of material in frame structures; Philosophical Magazine, 8 (1904), pp. 589–597</li>
<li><span id='bib0015'></span>
[[#bib0015|[3]]] Y.M. Xie, G.P. Steven; Evolutionary Structural Optimization (ESO);  (First Edition)Springer, London (1997)</li>
<li><span id='bib0020'></span>
[[#bib0020|[4]]] G.P. Steven, Y.M. Xie; Evolutionary Structural Optimization with FEA; Computational Mechanics, 1 (1993), pp. 27–34</li>
<li><span id='bib0025'></span>
[[#bib0025|[5]]] O. Querin; Evolutionary Structural Optimisation: Stress Based Formulation and Implementation. Ph. D. Dissertation; Department of Aeronautical Engineering, University of Sydney, Australia (April 1997)</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] F. Velázquez, S. Santillán; Diseño Óptimo Evolutivo; Ingeniería, Investigación y Tecnología, VII 3 (2006), pp. 139–150</li>
<li><span id='bib0035'></span>
[[#bib0035|[7]]] W. Annicchiarico, M. Cerrolaza; Structural Shape Optimization 3D Finite-Element Models Based on Genetic Algorithms and Geometric Modeling; Finite Elements in Analysis and Design, 37 (2001), pp. 403–415</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] C. Mattheck, S. Burkhardt; A new method of structural shape optimization based on biological growth; Intn. J. Fatigue, 12 (1990), pp. 185–190</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] C. Mattheck; Design in Nature- Learning from Trees;  (First Edition)Springer Verlag, Germany (1998)</li>
<li><span id='bib0050'></span>
[[#bib0050|[10]]] J.A. Sethian, A. Wiegmann; Structural Boundary Design via Level Set and Immersed Interface Methods; J. Comp. Phys., 163 (Sep 2000), pp. 489–528</li>
<li><span id='bib0055'></span>
[[#bib0055|[11]]] E. Ramm, K. Maute, S. Schwarz; Adaptative Topology and Shape Optimization; Computational Mechanics – New Trends and Applications, CIMNE, Barcelona, Spain (1998)</li>
<li><span id='bib0060'></span>
[[#bib0060|[12]]] K. Maute, E. Ramm; Adaptative Topology Optimization; Structural Optimizaction, 10 (1995), pp. 100–112</li>
<li><span id='bib0065'></span>
[[#bib0065|[13]]] E. Cervera, J. Trevelyan; Evolutionary structural optimization based on boundary representation of NURBS: Part II: 3D algorithms; Computers & Structures, 83 (23–24) (2005), pp. 1917–1929</li>
<li><span id='bib0070'></span>
[[#bib0070|[14]]] C. Shene, Introduction to computing with geometry. Course Notes. [WWW document; consultado 15 Ene 2010]. Disponible en: http://www.cs.mtu.edu/∼shene/COURSES/cs3621/NOTES/ , 1998, July 1                                    </li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] Q. Li, G.P. Steven, Y.M. Xie; On equivalence between stress criterion and stiffness criterion in evolutionary structural optimization; Structural Multidisciplinary Optimisation, 18 (1999), pp. 67–73</li>
<li><span id='bib0080'></span>
[[#bib0080|[16]]] H. Yang, W. Wang, J. Sun; Control Point adjustment for B-Spline Curve Approximation; Cumputer – Aided Design, 36 (2004), pp. 639–652</li>
<li><span id='bib0085'></span>
[[#bib0085|[17]]] Tom Compton, “Analytic Cycling – Forces on Rider” [documento WWW; consultado 10 Jul 2010], Disponible en: http://www.analyticcycling.com/ForcesPower_Page.html , 2001.                                    </li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] J.J. Mevdev; Design and Optimization of a Carbon Monocoque Bicycle Frame; 4th year Thesis, Department of Aeronautical Engineering, University of Sydney, Australia (1995)</li>
</ol>