Optimización topológica bidireccional con esquema de penalización de material y restricciones de desplazamiento/Bidirectional topological optimization with material penalization scheme and displacement constraints

Abstract

This work presents the topology optimization of bi-dimensional structures in plane stress state, considering a bidirectional scheme of optimization, incorporating a penalization parameter that maximizes the stiffness, considering the kinematic constraints. An open source code was implemented in MATLAB, with capabilities to find the optimal structures for a certain volume fraction, modified to obtain the optimal structure considering a limited or admissible displacement. Two classical problems in the literature are analyzed by introducing the developed concepts. The results obtained indicate that, for a particular state of charge, it is possible to define different configurations for maximum rigidity, depending on the displacement considered. Each one of these configurations is structured on the basis of compression and tension elements that resemble a lattice, with almost uniform tensional distribution in each element. The relevance of this work is that it allows you to display the existence of maximum configurations, depending on the limit or admissible displacement and how this converges to a structure with uniform stress levels.

Headnote

RESUMEN

Se presenta la optimización topológica de estructuras bidimensionales en estado de tensión plana, considerando un esquema de optimización bidireccional que incorpora un parámetro de penalización para maximizar la rigidez de la estructura considerando restricciones del tipo cinemático. Un código abierto de optimización topológica escrito en MAIL AB® con capacidades de encontrar el óptimo para una fracción de volumen definida a priori, se modifica para obtener la estmctura óptima cuando se impone un desplazamiento límite o admisible. Dos problemas clásicos de la literatura se analizan introduciendo los conceptos desarrollados. Los resultados obtenidos permiten concluir que para un estado particular de carga es posible definir distintas configuraciones de máxima rigidez dependiendo del desplazamiento considerado. Cada una de estas configuraciones se estructura a base de elementos en compresión y tracción que se asemejan a un enrejado, con distribuciones tensionales casi uniformes en cada elemento. La relevancia del presente trabajo radica en que permite visualizar la existencia de configuraciones de máxima rigidez dependiendo del desplazamiento límite o admisible elegido y como este converge hacia una estructura con tensiones uniformes.

Palabras clave: Método de los elementos finitos, optimización topológica, rigidez máxima, estado tensional, desplazamiento.

ABSTRACT

This work presents the topology optimization of bi-dimensional structures in plane stress state, considering a bidirectional scheme of optimization, incorporating a penalization parameter that maximizes the stiffness, considering the kinematic constraints. An open source code was implemented in MATLAB, with capabilities to find the optimal structures for a certain volume fraction, modified to obtain the optimal structure considering a limited or admissible displacement. Two classical problems in the literature are analyzed by introducing the developed concepts. The results obtained indicate that, for a particular state of charge, it is possible to define different configurations for maximum rigidity, depending on the displacement considered. Each one of these configurations is structured on the basis of compression and tension elements that resemble a lattice, with almost uniform tensional distribution in each element. The relevance of this work is that it allows you to display the existence of maximum configurations, depending on the limit or admissible displacement and how this converges to a structure with uniform stress levels.

Keywords: Finite element method, topological optimization, maximum stiffness, stress state, displacement.

(ProQuest: ... denotes formulae omitted.)

INTRODUCCIÓN

La optimización estructural [1] nace en la segunda mitad del siglo XX inspirada en los trabajos sobre optimización de estructuras reticuladas realizados por Michell [2]. El objetivo principal de la optimización de estructuras es obtener un conjunto de valores para las variables de diseño, que minimicen (o maximicen, según sea el caso) una función objetivo, y satisfagan un conjunto de restricciones que dependen de estas mismas variables.

En función de las variables a optimizar, y de la formulación inicial del problema, la optimización estructural se clasifica usualmente en: optimización de dimensiones, optimización deformas y optimización topológica [3-4].

La optimización de dimensiones busca optimizar el tamaño de los elementos previamente definidos por el diseñador. Su aplicación más usual es en las estructuras reticuladas, donde se definen las longitudes y conectividades, para luego optimizar la sección de las barras, que corresponde a la variable de diseño. La optimización de dimensiones es la forma más simple de optimización estructural, y por ende suele ser la más común de encontrar en la literatura. En la optimización de formas el dominio inicial es una estructura ya definida, pero con ciertos parámetros modificables en su forma, que conesponden a las variables de diseño como las coordenadas de los nodos de un reticulado o el contorno de una estructura continua.

En la optimización topológica ( OT) se busca encontrar la mejor distribución de material dentro de un volumen finito (dominio de diseño) de forma tal que la estructura resultante maximice algún tipo de desempeño bajo ciertas restricciones [5]. La particularidad de este tipo de optimización es que los diseños no dependen del criterio del diseñador, como sucede en la optimización de formas y de dimensiones, en donde se restringen las posibles soluciones en función del diseño inicial.

Diferentes técnicas basadas en el método de los elementos finitos se plantean para abordar el problema de OT de estructuras, diferenciándose entre ellas en el tratamiento dado a la variable de optimización, pudiendo ser esta del tipo discreto o continuo.

Entre los esquemas de variable discreta destacan las técnicas de optimización estructural evolutiva (Evolutionary Structural Optimization, ESO) [6] y el esquema de optimización estructural evolutiva bi-direccional (Bi-directional Evolutionary Structural Optimization, BESO) [7].

En el método ESO, según criterios del tipo tensional y/o basados en la energía de deformación, se aplica un procedimiento que gradualmente va eliminando las partes ineficientes del dominio de diseño. En el método BESO, no solo se remueve el material ineficiente, sino que también se agrega en aquellas zonas en donde se necesite, es decir, posee un carácter bidireccional. Ambos esquemas ESO/ BESO se denominan del tipo hard kill puesto que remueven completamente los elementos ineficientes del dominio de diseño.

Dentro de los esquemas de variable continua destaca el modelo de material sólido isótropo con penalización (Solid Isotropic Material with Penalization, SIMP) [8]. La popularidad de este esquema radica en que utiliza una única variable de optimización (la densidad relativa) para controlar el problema.

El objetivo de la metodología SIMP es maximizar la rigidez de la estructura para una fracción de volumen/definida a priori (/= V*/V", razón entre el volumen final e inicial). Esto se hace indirectamente por medio de la minimización de la energía de deformación expresada en función de la variable de diseño penalizada.

Tanto en la aplicación del esquema ESO, BESO y SIMP la topología final presenta recurrentemente dos tipos de problemas [3]: i) existencia de zonas del dominio en donde la distribución de material se asemeja a un tablero de ajedrez, es decir, elementos contiguos que alternan entre sólido y vacío, y ii) dependencia de la malla de elementos finitos respecto de la topología final obtenida, es decir, se obtienen distintas soluciones dependiendo de la malla de elementos finitos utilizada.

Para eliminar ambos efectos los algoritmos de OT incorporan un "filtro de sensibilidades'! 10-12] que básicamente consiste en modificar la derivada de la función objetivo a partir de la media ponderada en un vecindario fijo al elemento, proyectando así este valor desde los elementos llenos a los vacíos generando de esta forma estructuras continuas y en donde la configuración final no depende del tamaño de la malla utilizada.

Huang y Xie [13] proponen un esquema híbrido BESO-SIMP (soft kill) del tipo bidireccional que utiliza un esquema de interpolación del material con penalización y que plantea que en lugar de remover totalmente los elementos ineficientes, estos sean reemplazados por elementos de baja densidad. Este esquema requiere de una menor cantidad de iteraciones y converge de forma estable hacia el óptimo [14].

Considerando las ventajas del esquema híbrido BESO-SIMP, el presente artículo aborda el problema de definir la distribución óptima de material, de forma tal que la estructura resultante satisfaga que el desplazamiento en un nodo de control uj, sea menor o igual al desplazamiento admisible uj para un volumen final V* objetivo.

Lo anterior tiene gran relevancia en problemas de la ingeniería civil, en donde el diseño estructural está condicionado tanto por las tensiones que se generan al interior del dominio, como por los desplazamientos que se desarrollan para un cierto nivel de cargas.

PLANTEAMIENTO DEL PROBLEMA

Para lograr el propósito de la presente investigación se utiliza un código abierto de optimización topológica escrito en MATLAB®con capacidades para encontrar el óptimo para una fracción de volumen definida a priori (Huang y Xie [13] http://www.isg.rmit.edu. au/), el cual se modifica para obtener la estructura óptima cuando se impone un desplazamiento límite o admisible.

El código resultante [17] se divide en seis rutinas principales:

1. Preproceso ( generación de malla, definición de las propiedades del material, cargas aplicadas y condiciones de contorno).

2. Ingreso de parámetros de optimización (fracción de volumen, tasa de evolución del volumen y desplazamiento límite).

3. Análisis estructural por elementos finitos.

4. Rutina de BESO-SIMP (cálculo de la derivada de la función objetivo, aplicación del esquema de filtrado de sensibilidades, actualización de densidades relativas y análisis de convergencia).

5. Nuevo ciclo de análisis estructural por elementos finitos hasta que se satisfaga la condición que el cambio en la función objetivo sea menor a un error [varepsilon] y que el desplazamiento del nodo de control uj sea menor o igual al desplazamiento límite ti/' definido para ese punto.

6. Posproceso (nivel de desplazamientos, distribución de tensiones principales y evaluación de tensión efectiva de Von Mises).

Para verificar la validez de la metodología propuesta e implementada [17], se desarrollan dos ejemplos de la literatura de OT [12,15], en donde se verifican y analizan los resultados obtenidos.

FORMULACIÓN MATEMÁTICA

La formulación general del problema de optimización topológica se expresa por [16]:

... (1)

... (2)

... (3)

... (4)

... (5)

En donde c1(x) corresponde a la función objetivo (energía de deformación), ki y ui son la respectiva matriz de rigidez y vector de desplazamiento asociados al i-ésimo elemento respecto del sistema de coordenadas globales.

Las ecuaciones (2), (3) y (4) expresan las restricciones del problema, en donde xi es la densidad relativa del elemento (variable de diseño) expresada como la razón entre la densidad actual respecto de la del material solido (ρ/ρo).

Cabe señalar que la variable de diseño tiene un límite inferior xmin para evitar la singularidad de la matriz de rigidez y un límite superior xi=1 que indica material lleno. Finalmente, uj y u* denotan el desplazamiento j de un nodo arbitrario y su respectivo valor de restricción.

La matriz de rigidez del i-ésimo elemento en coordenadas naturales se expresa por:

... (6)

En donde Bi es la matriz de deformación elemento finito rectangular Lagrangiano de cuatro nodos (denominado Q4), ver Figura 1:

D1 corresponde a la matriz constitutiva de un material elástico lineal e isótropo en estado de tensión plana:

... (7)

Siguiendo la metodología propuesta por SIMP el módulo de elasticidad se interpola, de acuerdo con:

... (8)

En donde E0 es el módulo de elasticidad del material totalmente sólido y el parámetro "p" actúa como un factor de penalización que castiga el módulo de elasticidad en todos aquellos elementos que presenten una densidad relativa intermedia (0<xi<1), generando de esta forma estructuras con una configuración vacía-llena.

Se demuestra [18] que el factor de penalización en elasticidad bidimensional es función del módulo de Poisson v0 y se estima de acuerdo con:

... (9)

(6) Reemplazando ecuación (8) en (7) se obtiene:

... (10)

... (11)

Finalmente, reemplazando en la ecuación (6), se obtiene la matriz de rigidez del i-ésimo elemento:

... (12)

Siendo k0 la matriz de rigidez inicial del elemento totalmente sólido. Debido a que la malla del dominio de diseño es estructurada k0 es igual para todos los elementos de la discretización.

La ecuación (2) establece que en todo momento debe existir equilibrio estático en la estructura y la ecuación (3), más que una restricción es un parámetro de evolución el que se debe satisfacer iteración a iteración. Por último, la ecuación (5) establece el dominio dentro del cual se mueve la variable de diseño.

El lagrangiano asociado al problema se reduce a:

... (13)

En donde las condiciones de optimalidad establecen que en el óptimo, se cumple que:

... (14)

Deducciones detalladas de la sensibilidad de la energía de deformación y el desplazamiento respecto de una variación en la variable de diseño pueden encontrarse en [12]; sin embargo, la determinación del valor de λ requiere especial atención. Si iij<uf, λ=0 y si uj<uj*λ tiende a infinito lo que quiere decir que la topología propuesta no cumple con la restricción de desplazamiento y por ende λ crece para redistribuir el material y tratar de generar una topología que sí cumpla con la restricción impuesta. Para este propósito λ se determina en función del desplazamiento de la siguiente iteración de acuerdo con el esquema presentado en la Figura 2.

A partir de la ecuación ( 14) se define el número de sensibilidad elemental:

... (15)

Valor que se utiliza como parámetro de comparación para la actualización de la densidad relativa xi en cada iteración [12].

FILTRO DE SENSIBILIDADES

Para evitar la formación de zonas con una distribución material similar a la de un tablero de ajedrez y la dependencia del diseño respecto de la discretización de elementos finitos, se utiliza un filtro de sensibilidades, el cual modifica el número de sensibilidad elemental iteración a iteración.

La modificación del número de sensibilidad elemental se inicia calculando el número de sensibilidad nodal, el cual se obtiene después de promediar ponderadamente las respectivas sensibilidades elementales de cada uno de los elementos de la vecindad que comparten un nodo común.

El algoritmo matemático que permite estimar el número de sensibilidad nodal se expresa por:

... (16)

Con:

... (17)

En donde M es la cantidad de elementos finitos conectados con el j-ésimo nodo y wi el factor de peso correspondiente al i-ésimo elemento. Siendo rij la distancia entre el centro del i-ésimo elemento al j-ésimo nodo.

En el caso particular de una malla estructurada de elemento finitos del tipo Q4, el número de sensibilidad nodal se obtiene después de aplicar la ecuación (16) a cada uno de los cuatro nodos del elemento, ver Figura 3.

La zona de influencia, denominada Ωi, queda definida por un círculo de radio ri = rmin centrado en el i-ésimo elemento. El valor del rmin se elige de forma tal que el tamaño del subdominio permita cubrir más de un elemento. Este esquema permite reevaluar el número de sensibilidad y obtener:

... (18)

Donde K es el número total de nodos contenidos en el subdominio y \v(r,i) es el factor de peso definido como:

w(rij = rB¿n-ry (19)

De este modo, aquellos nodos más cercanos al centro del elemento tendrán un mayor peso en el cálculo del número de sensibilidad ai que aquellos nodos más alejados del centro.

ACTUALIZACIÓN DE DENSIDADES

Para actualizar la densidad relativa xi de cada elemento, se ordena el número de sensibilidad ai de menor a mayor y se define una sensibilidad de comparación elemental ?com cual define un umbral a partir del cual se "eliminan" y/o se agregan o conservan elementos.

El esquema de actualización de la densidad relativa se obtiene por:

... (20)

En donde ?com se estima mediante un algoritmo de bisección, basándose en cumplir la restricción de volumen para esa iteración en particular. Esto porque la fracción de volumen también evoluciona a lo largo del proceso, comenzando usualmente con el valor máximo (todos los elementos con densidad relativa 1) y descendiendo en forma lineal a una la tasa de evolución, pero siempre mayor o igual a una fracción de volumen final definida a priori la que no puede ser inferior al 25% del volumen inicial [17].

Para valores menores a/0 = 0,25 la cantidad de material disponible en el dominio resulta ser muy baja y el algoritmo no logra la convergencia, oscilando la solución entre dos topologías. Adicionalmente resultan topologías con configuraciones a base de elementos (puntales y tensores) muy esbeltos, que como sistema estructural global genera estructuras con desplazamientos excesivos, sin sentido físico.

La fracción de volumen evoluciona de acuerdo con:

/ = máx {/(l - ER), f0} (21)

El ciclo de análisis estructural y la adición y sustracción de elementos continúa, hasta que se satisfaga la condición que el desplazamiento uj sea menor o igual al desplazamiento límite u *.

EJEMPLOS NUMÉRICOS

Con el fin de mostrar el potencial de la metodología de OT implementada, se presentan dos problemas de OT considerando su solución con y sin restricción local de desplazamiento.

Viga alta simplemente apoyada

Se busca definir la topología óptima de una viga alta simplemente apoyada (Figura 4), considerando que el desplazamiento horizontal de un punto ubicado sobre el apoyo móvil (Punto A) no supere 1 mm. La viga tiene una razón de aspecto H/L=1/2, un espesor de 1 ímn, una carga puntual es su centro de 100 N [12] y las propiedades mecánicas y parámetros de optimización se definen en Tabla 1. Además, se aprovecha la simetría del problema para aumentar la eficiencia computacional.

En la Figura 5 se presenta la evolución del desplazamiento del punto de control (nodo A) en función de la evolución de la fracción volumétrica, tanto para el caso sin restricción de desplazamiento como para la situación cuando se limita el valor del desplazamiento horizontal.

Se observa que para fracciones de volumen/> 0,56 ambos problemas (con y sin restricción) llegan a los mismos desplazamientos y por ende a las mismas topologías, lo que se debe a que la restricción de desplazamiento no se ha activado (restricción pasiva), sin embargo, para/< 0,56 las topologías empiezan a diferir, lo que se debe a que la restricción ahora es activa y afecta la topología final para cada fracción /establecida.

En la Figura 6 se presentan las topologías óptimas, para el caso sin y con restricción de desplazamiento, para una fracción volumétrica del 30%. En la Figura 7 se presentan las isolíneas de tensiones para ambos casos asociadas a la fracción de volumen antes indicada.

Se concluye que pese a la diferencia de desplazamientos horizontales del nudo A para una misma fracción de volumen (Figura 7), la magnitud de las tensiones principales alcanza valores similares al interior de la estructura optimizada.

Viga postensada simplemente apoyada

Considérese una viga de luz igual a 12 m y de altura de 1,2 m, la cual se optimiza topológicamente, considerando la acción de postensado longitudinal [15], de acuerdo con el modelo indicado en Figura 8 y considerando las propiedades indicadas en la Tabla 2.

Considerando las condiciones de simetría de la viga (cargas y condiciones de contorno) se analiza solamente la mitad de ella y se considera un borde de 0,15 m de ancho en todo el perímetro como no optimizable. El punto de control del desplazamiento vertical corresponde al ubicado sobre el borde superior e indicado como A en la Figura 8.

En la Figura 9 se presenta la evolución del desplazamiento del punto de control (punto A) en función de la evolución de la fracción volumétrica, sin considerar restricción de desplazamiento y considerando una restricción de desplazamiento vertical admisible de 40 nun. En la Figura 10 se presenta la topología óptima asociada a una fracción de volumen del 50% cuando a) no se limita el desplazamiento vertical y cuando b) se limite a un valor de 40 ímn.

Finalmente en Figura 11 se presentan para el mismo ejemplo las isolíneas de tensiones efectivas de Von Mises para ambos casos.

En resumen: en ambos casos (con/sin aplicar limitación en el desplazamiento vertical del punto A) la técnica de OT permite entender los mecanismos de distribución de tensiones al interior del elemento. CONCLUSIONES

Se ha presentado una metodología BESO-SIMP para la optimización topológica de estructuras con restricciones de desplazamiento. El esquema resulta altamente eficiente y permite estimar la estructura más rígida que satisface un determinado desplazamiento admisible impuesto para una respectiva fracción de volumen dada.

La relevancia del trabajo se traduce en que es posible definir configuraciones de máxima rigidez dependiendo del desplazamiento admisible que se considere, es decir, para una fracción de volumen/ se pueden obtener diferentes topologías dependiendo del desplazamiento admisible (o límite) que se establezca.

Se verifica además, que a medida que la fracción de volumen disminuye, la estructura evoluciona hacia configuraciones en forma de enrejado a base de elementos de "barra" sometidos a compresión y tracción con distribuciones tensionales cada vez más uniformes.

AGRADECIMIENTOS

Los autores agradecen el aporte financiero de CONICYT por medio del proyecto FONDECYT 1100374.

References

REFERENCIAS

[1] M. Victoria. "Optimización de forma y topología con malla fija y algoritmos genéticos". Departamento de Estructuras y Construcción. Universidad Politécnica de Cartagena. Cartagena, España. 2006.

[2] A.G.M. Michell. "The limit of economy of material in frame structures". Philosophical Magazine. Vol. 8, Issue 47, pp. 589-597.1904.

[3] M.P. Bendsøe and O. Sigmund. "Topology Optimization Theory, Method and Applications". Springer Verlag. Berlín, Alemania. 2003. ISBN: 3-540-42992-1.

[4] K. Ghabraie."Thesis for the degree of doctor of philosophy".School of civil, environmental and chemical engineering science, engineering and Health (SEH) portfolio.RMIT University. Melbourne, Australia. 2009.

[5] A. Tovar. "Optimización topológica con la técnica de los autómatas celulares híbridos". Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 21 N 4, pp. 365-383. 2005.

[6] Y.M. Xiey and G.P. Steven. "A simple evolutionary procedure for structural optimization". Computers and Structures. Vol. 49, Issue 5, pp. 885-896. 1993.

[7] X.Y. Yang, Y.M. Xie, G.P. Steven and O.M. Querin. "Bidirectional evolutionary method for stiffness optimization". AIAA Journal. Vol. 37, Issue 11, pp. 1483-1488. 1999.

[8] M. P.Bendsøe and N. Kikuchi. "Generating optimal topology in structural design using a homogenization method". Computer Methods in Applied Mechanics and Engineering. Vol. 71, Issue 2, pp. 197-224. 1988.

[9] A. Díaz and O. Sigmund. "Checkerboard patterns in layout optimization". Structural Optimization. Vol. 10, pp. 40-45. 1995.

[10] O. Sigmund. "On the design of compliant mechanisms using topology optimization". Mechanics of Structures and Machines. Vol. 25, Issue 3, pp. 493-524. 1997.

[11] O. Sigmund and J. Peterson. "Numerical instabilities in topological optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima". Structural and Multidisciplinary Optimization. Vol. 16, pp. 68-75. 1998.

[12] X. Huang and Y.M. Xie. "Evolutionary Topology Optimization of Continuum Structures". Methods and Applications. Wiley. Chichester. 2010. ISBN: 978-0-470-74653-0 (cloth).

[13] X. Huang and Y.M. Xie. "Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials". Computational Mechanics. Vol. 43, Issue 3, pp. 393-401. 2009.

[14] X. Huang and Y.M. Xie. "A revised BESO method for structures with design-dependent gravity loads". The fifth international structures engineering and construction conference, Challenges, Opportunities and Solutions in Structural Engineering and Construction. London, United Kingdom. September, 2009.

[15] P. Dombernowsky and A. Sondergaard. "Three-dimensional topology optimization inarchitectural and structural design of concrete structures". Proceedings of the International Association for Shell and Spatial Structures (IASS) Symposium. Valencia, España. Septiembre-Octubre, 2009.

[16] P. Tanskanen. "The evolutionary structural optimization method: theoretical aspects". Computer Methods in Applied Mechanics and Engineering. Vol. 191, pp. 5485-5498. 2002. "

[17] J. Sanhueza. "Optimización estructural evolutiva bi-direccional con penalización material y restricciones de desplazamiento". Memoria para optar al título de Ingeniero Civil. Departamento de Ingeniería Civil. Universidad de Concepción. Concepción, Chile. 2012.

[18] Z. Hashin and S. Shtrikman. "A variational approach to the theory of the elastic behavior of multiphase materials". Journal of Mechanics and Physics of Solids. Vol. 11, pp. 127-140. 1963.

[19] ACI318-08. "Requisitos de Reglamento para Concreto Estructural". American Concrete Institute. 2008.

[20] Q. Li, G.P. Steven and Y.M. Xie. "A simple checkerboard suppression algorithm for evolutionary structural optimization". Structural and Multidisciplinary Optimization. Vol. 22, pp. 230-239. 2001.

[21] R.B. Harber, C.S. Jog and M.P. Bendsøe. "A new approach to variable-topology shape design using a constraint on the perimeter". Structural Optimization. Vol. 11, pp. 1-11.1996.

[22] C.S. Jog. "Topology design of structures using a dual algorithm and a constraint on the perimeter". International Journal of Numerical Methods in Engineering. Vol. 54, pp. 1007-1019. 2002.