The influence of self‐weight of elastic 2D structures in topology optimization via numerical technique Smooth Evolutionary Structural Optimization (SESO)

Resumo

O presente artigo aborda a otimização topológica em problemas de elasticidade plana linear considerando a influência do peso próprio nos esforços em elementos estruturais. Utiliza‐se para este fim uma técnica numérica denominada Smooth ESO (SESO) que se baseia no procedimento de diminuição progressiva da contribuição de rigidez de elementos ineficientes com menores tensões até que ele não tenha mais influência. As aplicações do SESO são feitas com o método dos elementos finitos e considera‐se um elemento finito triangular e de alta ordem. Neste trabalho estende‐se a técnica SESO para a aplicação do peso próprio onde o programa, no cômputo de seu volume e peso específico, gera automaticamente uma força concentrada equivalente para cada nó do elemento. A avaliação é finalizada com a definição de um modelo de bielas e tirantes resultante das regiões de concentração de tensões. Nos exemplos de aplicação são apresentadas topologias ótimas de uma estrutura suspensa, de viga baixa e de viga parede considerando o peso próprio e obtendo‐se ótimas configurações e demonstrando que a consideração do peso próprio leva a maior robustez ao processo de otimização.

Abstract

This paper deals with topology optimization in plane elastic‐linear problems considering the influence of the self weight in efforts in structural elements. For this purpose it is used a numerical technique called SESO (Smooth ESO), which is based on the procedure for progressive decrease of the inefficient stiffness element contribution at lower stresses until he has no more influence. The SESO is applied with the finite element method and is utilized a triangular finite element and high order. This paper extends the technique SESO for application its self weight where the program, in computing the volume and specific weight, automatically generates a concentrated equivalent force to each node of the element. The evaluation is finalized with the definition of a model of strut‐and‐tie resulting in regions of stress concentration. Examples are presented with optimum topology structures obtaining optimal settings.

Palavras‐chave

Modelos de biela‐tirante ; Otimização evolucionária topológica ; Estruturas de concreto armado ; Peso próprio

Keywords

Strut and tie models ; Topology optimization ; Reinforced concrete structures ; Self weight

1. Introdução

Atualmente na engenharia estrutural os métodos numéricos de aplicação são ferramentas imprescindíveis e já nas 2 últimas décadas as técnicas de otimização vêm sendo empregadas no auxílio da modelagem estrutural. Dentre os métodos tradicionais de otimização, a otimização de forma e a paramétrica não são capazes de alterar a configuração da estrutura original e, desta maneira, não são possíveis de determinar os fluxos de tensões. A grande vantagem da otimização topológica (OT) é tornar possíveis estas alterações.

Os conceitos básicos da técnica OT foram estabelecidos há mais de um século, conforme descrevem Rozvany et al. [1] . Na análise topológica, 2 metodologias vêm se destacando, a abordagem micro e abordagem macro. A filosofia da técnica baseia‐se na remoção da parte do elemento que não for realmente necessário à estrutura. Esta remoção pode ser realizada de forma suave, atenuando os valores da matriz constitutiva do elemento.

Na abordagem microestrutural, o processo de otimização topológica baseado no estabelecimento de um relação entre a rigidez e a densidade associada ao domínio, que pode assumir qualquer valor entre 0 (vazio) e 1 (material sólido), sendo os valores intermediários correspondentes a um material poroso. A partir daí os algoritmos baseados em microestruturas propõem encontrar a melhor disposição do material, de modo a minimizar ou maximizar a função custo, Stump [2] . Como exemplo pode‐se citar: o método Simple Isotropic Material with Penalization (SIMP), utilizado por Bendsøe [3] , Rozvany et al. [4] e [1] .

Na abordagem macro, a topologia da estrutura é modificada mediante a inserção de furos no domínio. Como exemplo deste grupo de OT pode‐se citar o Evolutionary Structural Optimization (ESO), baseado no cálculo da função objetivo quando um elemento é removido da malha de elementos finitos, e Topological Sensitivity Analysis (TSA), que avalia a sensibilidade de uma função quando um furo é criado no corpo analisado, Labanowski et al. [5] .

Na análise estrutural de elementos em concreto armado a analogia apresentada por Ritter e Morsch no início do século XX , que associa uma viga a uma estrutura treliçada equivalente, ainda é utilizada para representar o comportamento estrutural deste elemento fletido. A denominada teoria de treliça clássica define barras representando os fluxos de tensão procurando criar o menor e o mais lógico caminho das cargas. Os elementos discretos (barras) representam os campos de tensão de tração (denominados tirantes) e de compressão (denominados bielas comprimidas) que surgem internamente ao elemento estrutural no efeito de flexão. Este analogia foi melhorada e é ainda utilizada pelas normas técnicas no dimensionamento de vigas em concreto armado à flexão e à força cortante e que estabelecem critérios diversos para a determinação dos limites de segurança nos seus procedimentos.

O trabalho de Schlaich et al. [6] descreve o modelo de bielas e tirantes de maneira mais generalizada, abrangendo o modelo de treliça de flexão e englobando regiões com descontinuidades e elementos estruturais especiais. Estas regiões, denominadas «Regiões de Descontinuidade D», onde as tensões de cisalhamento são significativas e a distribuição das deformações na seção transversal não é linear. Mesmo tratando de um modelo de idealização simples, é bastante utilizado na prática e eficiente para o conhecimento das regiões de concentração de tensões.

A utilização de ferramentas computacionais pode trazer maior confiabilidade e rapidez na definição de configurações geométricas e do fluxo de tensões dos elementos discretos no interior do elemento estrutural, de modo que a aplicação da OT no campo da análise estrutural em concreto armado nos modelos de bielas e tirantes vem sendo amplamente desenvolvido, destacando os trabalhos de Ali [7] , Liang et al. [8] , Liang et al. [9] , Liang et al. [10] , Liang e Steven [11] , Liang [12] e Brugge [13] .

Visando o desenvolvimento dos modelos de bielas e tirantes, faz‐se uma aplicação da técnica de otimização topológica Smooth ESO (SESO), seguindo os trabalhos de Simonetti et al. [14] , Simonetti et al. [15] na análise de elementos estruturais em concreto armado, realizando uma comparação das configurações ótimas considerando‐se a influência do peso próprio, o que muitas vezes pode levar a uma configuração final bem distinta quando da sua não consideração.

Trata‐se de uma variante do método ESO, cuja filosofia baseia‐se na observação de que se o elemento não for realmente necessário à estrutura, naturalmente sua contribuição de rigidez vai diminuindo progressivamente, até que ele não tenha mais influência. Isto é, sua heurística de remoção é feita de forma suave, atenuando os valores da matriz constitutiva do elemento, como se este estivesse em processo de danificação, sendo capaz de gerar os membros ideais de bielas e tirantes de forma mais natural. São apresentados 3 exemplos clássicos obtidos na literatura, que demonstram que o procedimento desenvolvido para o peso próprio é eficiente.

2. Otimização estrutural evolucionária

Xie e Steven [16] desenvolveram uma maneira bem simples de impor modificações na topologia da estrutura, feita mediante heurística de remoção gradual e sistemática de elementos finitos da malha, correspondentes a regiões que não contribuem efetivamente para o bom desempenho da estrutura.

Define‐se inicialmente uma malha de elementos finitos que circunscreva toda a estrutura, ou domínio estendido do projeto, de forma a incluir as condições de contorno em forças e deslocamentos, cavidades e demais condições iniciais. Em um processo iterativo, avaliam‐se os parâmetros de interesse de otimização, neste artigo em especial é a diminuição da massa mediante um critério de tensão máxima da estrutura. Assim, avaliam‐se as tensões de cada elemento, conforme o uso da inequação:

${\sigma }_{\mbox{e}}^{\mbox{vm}}<{\mbox{RR}}_{k}{\mbox{.}}\left({\sigma }_{\mbox{i,max}}^{\mbox{vm}}\right)$

( 1)

$RR_{k+1}=RR_{k}+RE$

( 2)

com ${\textstyle {\sigma }_{\mbox{e}}^{\mbox{vm}}{\mbox{e}}{\mbox{ }}{\sigma }_{i{\mbox{,max}}}^{\mbox{vm}}}$ sendo, respectivamente, as tensões principais de Von Mises do elemento «e» e máximas da estrutura na iteração efetiva «k», RR_k é raio de rejeição (RR) do k‐ésimo estado de equilíbrio.

Em cada iteração os elementos que satisfazem a inequação (1) são retirados da estrutura, figura 1 . O fator RR é aplicado para controlar o processo da remoção da estrutura (0,0 ≤ RR ≤ 1,0). O mesmo ciclo de remoção dos elementos usados pela desigualdade (1) é repetido até que não haja mais elementos que satisfaçam mais esta desigualdade (1). Quando esta situação ocorre, um estado de equilíbrio é alcançado. O processo evolucionário é definido adicionando a RR uma razão de evolução (RE). Assim, um novo ciclo de evolução inicia‐se, até que não existam mais elementos a serem eliminados com esta nova RR. A RR será atualizada conforme a equação (2), até àobtenção de uma configuração otimizada, alcançada pelo controle de um parâmetro de desempenho, denominado de índice de performance (IP). Este procedimento, também é conhecido como um método hard‐kill e pode ser interpretado como segue:

$D(j)={\begin{array}{c}D_{0}{\mbox{ }}{\mbox{if}}{\mbox{ }}{\mbox{ }}{\mbox{j}}\in \Gamma \\0{\mbox{ }}{\mbox{if}}{\mbox{ }}{\mbox{ }}{\mbox{j}}\in {\overline {\Gamma }}\end{array}}$

( 3)

Figura 1. Algoritmo evolucionário: baseado no elemento removido da malha.

onde D_i (j ) é a matriz constitutiva do elemento j  ∈ Ω na i‐ésima iteração, D₀ é a matriz constitutiva inicial do projeto, a qual é definida pra cada elemento finito, ${\textstyle \Omega ={\Gamma }_{i}+{\overline {\Gamma }}_{i}}$ é o domínio da estrutura, tal que ${\textstyle {\Gamma }_{i}=\lbrace \Omega /({\sigma }_{e}/{\sigma }_{VM}^{MAX}(\Omega ))\geq RR_{i}\rbrace }$ é o conjunto dos elementos que não serão removidos, e ${\textstyle {\overline {\Gamma }}_{i}=\Omega -{\Gamma }_{i}=\lbrace \Omega /({\sigma }_{e}/{\sigma }_{VM}^{MAX}(\Omega ))<RR_{i}\rbrace }$ é o conjunto dos elementos que serão removidos da estrutura (criação do vazio), todos na i ‐ésima iteração.

Na heurística de remoção via ESO, quando o elemento é removido do domínio de projeto durante o processo evolucionário, os elementos que permanecem na estrutura representam uma solução básica, conforme fora demonstrado matematicamente por Tanskanen [17] . Ele demonstra que a função discreta a ser minimizada com restrição por menor peso converge iterativamente, conforme a retirada de um elemento, que significa o encontro de um ponto de minimo.

3. Suavização da otimização estrutural evolucionária

Por outro lado, Tanskanem [17] destaca que a retirada de um elemento pode afetar a convergência do procedimento de otimização, pois o critério de retirada no ESO é indicado pelo atendimento da inequação (1), que muitas vezes pode ser radical, hard‐kill , uma vez que existem elementos que estão na vizinhança da esquerda dessa condição que numericamente são excluídos ou mantidos, mas que deveriam ainda compor ou ser retirados da estrutura, o que afeta o processo evolutivo.

Nesse sentido, é adequado aplicar sobre o método ESO uma condição de relaxação ou de suavização do ESO, apresentando em Simonetti et al. [15] um procedimento soft‐kill , de modo que o material, ou seja, os elementos que deveriam ser removidos pelo critério ESO – seguindo a inequação (1) – são organizados em n grupos e alocados em ordem crescente de tensões sendo ponderados segundo uma função 0 ≤ η (j ) ≤ 1.

Assim, p% destes n grupos são removidos, os grupos que contém os elementos com menores tensões (domínio Γ_LS ), e os grupos (1‐p%) são devolvidos a estrutura, domínio Γ_GS . Este processo de remoção e devolução de elementos à estrutura pode ser orientado por uma função, linear ou hiperbólica, que pondera a razão ${\textstyle {\sigma }_{e}^{VM}/{\sigma }_{\mbox{VM}}^{\mbox{MAX}}}$ dentro do domínio ${\textstyle {\overline {\Gamma }}}$ , isto é, permite que elementos que tenham suas tensões próximas e abaixo da tensão máxima, ${\textstyle {\sigma }_{\mbox{VM}}^{\mbox{MAX}}}$ sejam reintegrados a estrutura a cada iteração. A minimização da função objetivo é atingida mediante o encontro de um ponto ótimo estacionário, conforme [17] .

Uma forma alternativa, mais cara computacionalmente, de corrigir esse desvio seria a possibilidade de inserção do elemento na estrutura novamente, destacando que nesse sentido tem‐se uma variante do ESO, o Bidirectional ESO (BESO), Querin [18] , que possibilita a entrada e remoção de elementos no domínio de projeto. O SESO vem dessa filosofia, de forma consistente matematicamente e mais eficiente, que pondera a matriz constitutiva, fazendo com que o elemento que a priori seria eliminado no ESO é mantido e recebe uma relaxação, aplicando uma fração de sua rigidez inicial, permanecendo no projeto e naturalmente sua influência contribui para sua permanência ou fortalece a sua retirada por definitivo da malha inicial.

Os elementos que estão próximos ao limite pela esquerda desta tensão máxima são mantidos na estrutura, definindo assim um procedimento de retirada «não radical» e sim, de forma suave. O procedimento soft‐kill usado no SESO pode ser interpretado assim:

$D_{i}(j)=\left\{{\begin{array}{c}D_{0}{\mbox{ }}se{\mbox{ }}{\mbox{ }}j\in {\Gamma }_{i}\\D_{0}\cdot {\eta }_{j}({\overline {\Gamma }}){\mbox{ }}se{\mbox{ }}{\mbox{ }}j\in {\Gamma }_{GS}\\0{\mbox{ }}se{\mbox{ }}{\mbox{ }}j\in {\Gamma }_{LS}\end{array}}\right\}$

( 4)

onde ${\textstyle {\overline {\Gamma }}={\Gamma }_{LS}+{\Gamma }_{GS}}$ é o domínio dos elementos que atendem a inequação (1), ${\textstyle {\Gamma }_{LS_{i}}}$ é o domínio dos elementos que devem ser excluídos efetivamente, ${\textstyle {\Gamma }_{GS_{i}}}$ o domínio dos elementos que são devolvidos a estrutura, ${\textstyle 0\leq \eta ({\overline {\Gamma }})\leq 1}$ é uma função reguladora que pondera o valor da razão ${\textstyle {\sigma }_{e}^{VM}/{\sigma }_{VM}^{MAX}}$ dentro do domínio ${\textstyle {\overline {\Gamma }}}$ .

A proposta da suavização pode ser desempenhada por uma função linear do tipo ${\textstyle \eta ({\overline {\Gamma }})=\alpha j+\beta }$ ou uma função trigonométrica do tipo ${\textstyle \eta ({\overline {\Gamma }})=sen(\alpha j)}$ . Porque as 2 funções são contínuas e podem ser diferenciadas em todo o domínio ${\textstyle {\overline {\Gamma }}}$ e tem sua imagem variando de 0‐1, figura 2 .

Figura 2. Ilustração da suavização do volume dos elementos removidos na iteração i.

4. Índice de performance para a formulação Smooth Evolutionary Structural Optimization

O índice de performance (IP) é um parâmetro adimensional que mede a eficiência da performance da estrutura. O problema consiste na minimização da função objetivo em termos do peso, sujeito a uma restrição de tensão admissível (σ^projeto ), que é definida como:

${\mbox{minimize}}{\mbox{ }}W=\sum _{e=1}^{NE}w_{e}(t_{e})$ $sujeito{\mbox{ }}a{\mbox{ }}{\sigma }_{j,max}^{vm}-$${\sigma }^{project}\leq 0$

( 5)

onde NE é o número total de elementos finitos.

O IP foi proposto por Liang et al. [5] como:

$IP=\left({\frac {{\sigma }_{0,max}^{vm}}{{\sigma }_{i,max}^{vm}}}\right)\cdot {\frac {W_{0}}{W_{i}}}=$$\left({\frac {{\sigma }_{0,max}^{vm}}{{\sigma }_{i,max}^{vm}}}\right)\cdot {\frac {{\rho }_{0}\cdot V_{0}}{{\rho }_{i}\cdot V_{i}}}=$$\left({\frac {{\sigma }_{0,max}^{vm}}{{\sigma }_{i,max}^{vm}}}\right)\cdot {\frac {V_{0}}{V_{i}}}$

( 6)

onde V₀ e V_i são o volume inicial e o volume da i‐ésima iteração, ${\textstyle {\sigma }_{\mbox{0,max}}^{\mbox{vm}}}$ e ${\textstyle {\sigma }_{\mbox{i,max}}^{\mbox{vm}}}$ são a tensão inicial e a tensão máxima de Mises na i ‐ésima iteração, ρ₀ e ρ_i são o peso específico inicial e o peso específico na i ‐ésima iteração, como são iguais para um material incompressível. A suavização gerada devido à equação (4) em termos da matriz constitutiva pode ser escrita em termos de espessura, devido à relação de dependência linear entre elas. Neste contexto, o IP é a equação (7), que leva em conta a expressão (4) em termos de espessura e cada função reguladora do procedimento SESO:

$IP=\left({\frac {{\sigma }_{0,max}^{vm}}{{\sigma }_{i,max}^{vm}}}\right)\cdot {\frac {A_{0}\cdot t_{0}}{\sum _{j=1}^{NE}A_{j}\cdot t_{j}}}=$$\left({\frac {{\sigma }_{0,max}^{vm}}{{\sigma }_{i,max}^{vm}}}\right)\cdot {\frac {A_{0}\cdot t_{0}}{\sum _{j=1}^{NE}A_{j}\cdot t_{j}\cdot \eta (j)}}$

( 7)

onde t₀ é a espessura inicial e t_j é a espessura do elemento na iteração j . O controle do ótimo é feito por este IP, pois ele é um fator de «monitoramento» da região ótima de projeto. O controle da maximização deste parâmetro refere‐se ao controle de minimização do volume, assim, caso este índice caia de forma acentuada é uma forte indicação que se passou por uma configuração de ótima ou um ótimo estacionário é alcançado, no entanto, não se pode garantir que este ótimo é um ótimo global, mas sim uma configuração ótima de projeto de engenharia.

5. Peso próprio

O peso próprio de elementos estruturais correspondente à ação da gravidade sobre os corpos e constitui‐se em um dos carregamentos que devem ser considerados em seu dimensionamento estrutural. Pode‐se observar na figura 3 como se dá a ação do peso próprio.

Figura 3. Ação do peso próprio [19] .

Pode‐se notar que nos elementos estruturais dispostos na posição horizontal o peso próprio constitui‐se em carregamento transversal ao eixo (fig. 3 ), desenvolvendo os esforços de flexão, momento fletor e esforço cortante.

No caso dos elementos dispostos na vertical, o peso próprio (G) atua na direção do eixo longitudinal da peça e provoca esforço axial, que pode ter um efeito diferenciado dependendo da sua vinculação.

Nos elementos suspensos (tirantes, fig. 4 a) o efeito do peso é de tração e nas apoiadas (pilares, fig. 4 b) este efeito é de compressão. O peso próprio de uma peça (G) pode ser calculado multiplicando‐se o volume da mesma pelo peso específico do material, conforme indica a equação 8:

$G=A\rho L$

( 8)

• A ‐ área da seção transversal da peça
• L ‐ comprimento
• ρ ‐ peso específico do material

Figura 4. Ação do peso próprio [19] .

Na tração ou compressão axial a não consideração do peso próprio pode ser eventual e é o caso mais simples.

6. Formulação do elemento finito de alta ordem

O elemento finito empregado é advindo da formulação livre, destacando‐se que toda a formulação foi baseada no trabalho de Bergan e Felippa [20] . Ele é de geometria triangular com 3 graus de liberdade por vértice: 2 deslocamentos translacionais e uma de rotação azimutal, vetor rotacional ortogonal ao plano do elemento. Inicialmente divide‐se o campo de deslocamentos em 2 modos: U_b e U_s , respectivamente modos de deslocamentos básico e superior. O modo básico de deslocamento tem o objetivo de atender as condições de movimento de corpo rígido e de deformação constante com suas funções interpoladoras lineares. O modo de deslocamento superior, ou de alta ordem, são polinômios quadráticos em termos de x e y, levando a campos de tensões de variação lineares. Destarte, a matriz de rigidez final estabelecida como:

$K=K_{b}+K_{s}$

( 9)

A matriz dos modos básicos, K_b , para o caso bidimensional, é dada por:

$K_{b}={\frac {1}{A}}.L.D_{m}.L^{T}$

( 10)

onde L_j é denominada matriz «amontoadora» (lumped) , expressa por:

$L_{j}={\frac {1}{2}}\left[{\begin{array}{ccc}y_{ki}&0&x_{ik}\\0&x_{ik}&y_{ki}\\{\frac {\alpha }{6}}\left(y-ji2y_{kj}^{2}\right)&{\frac {\alpha }{6}}\left(x-ij2x_{jk}^{2}\right)&{\frac {\alpha }{3}}\left(x_{ij}.y_{ji}-x_{jk}.y_{kj}\right)\end{array}}\right]$

( 11)

onde x_ij  = x_i  − x_j , x_ij  = x_i  − x_j e i , j  e k denotam a permutação cíclica de 1,2 e 3. D_m é a matriz constitutiva do elemento. É fundamentalmente a matriz constitutiva da elasticidade, C, integrada ao longo da espessura do elemento. Esta integração reduz o volume V para A , área do elemento.

Para um elemento de espessura t , material com módulo de elasticidade longitudinal E e coeficiente de Poisson ν , a matriz D_m é dada por:

$D_{m}={\frac {E.t}{1-{\nu }^{2}}}\left[{\begin{array}{ccc}1&\nu &0\\\nu &1&0\\0&0&{\frac {1-\nu }{2}}\end{array}}\right]$

( 12)

onde ${\textstyle C={\frac {1}{t}}\cdot D_{m}}$ .

O parâmetro α  ≥ 0 introduzido na formulação da matriz L aparece como um fator de escala nas funções de forma, multiplicando os termos relacionados com os graus de liberdade rotacionais. O valor de α afeta a matriz dos modos básicos e o cálculo das deformações.

Na equação 11, se α  = 0, tem‐se a matriz L para o elemento finito triangular com deformação constante (CST). Para este elemento as forças nodais são relacionadas apenas com translações.

A matriz dos modos de alta ordem é dada por:

$K_{s}=H^{T}.K_{qs}.H_{s}$

( 13)

aqui H  = G⁻¹ , onde G é a matriz que relaciona os coeficientes associados com os deslocamentos nodais e os parâmetros generalizados q, fazendo‐se as devidas substituições das coordenadas nodais nas funções de forma.

Para um elemento bidimensional, a matriz generalizada de alta ordem é:

$K_{qs}={\int }_{A}B_{s}^{T}.D_{m}.B_{s}dA$

( 14)

onde B_S é a matriz de deformação para os modos de alta ordem, definida como:

$B_{s}=\Delta N_{s}=\lambda .\left[{\begin{array}{ccc}2a_{11}\xi +a_{21}\eta &2a_{12}\xi +a_{22}\eta &2a_{13}\xi +a_{23}\eta \\b_{21}\xi +2b_{31}\eta &b_{22}\xi +2b_{32}\eta &b_{23}\xi +2b_{33}\eta \\-4b_{31}\xi -4a_{11}\eta &-4b_{32}\xi -4a_{12}\eta &-4b_{33}\xi -4a_{13}\eta \end{array}}\right]$

( 15)

onde Δ é o operador diferencial para deformações de membrana, expresso por:

$\Delta =\left[{\begin{array}{cc}{\frac {\partial }{\partial x}}&0\\0&{\frac {\partial }{\partial y}}\\{\frac {\partial }{\partial y}}&{\frac {\partial }{\partial x}}\end{array}}\right]$

( 16)

e λ é um fator de escala com dimensões iguais ao inverso do comprimento, dada pela equação:

$\lambda ={\frac {1}{\sqrt {A}}}$

( 17)

com A, a área do triângulo:

$A={\frac {1}{2}}\cdot (x_{j}y_{k})$

( 18)

Os outros elementos da matriz B_s são expressos da seguinte maneira:

$N_{si}=\left\{{\begin{array}{c}u_{si}\\v_{si}\end{array}}\right\}=$ $=\left[{\begin{array}{ccc}a_{1i}&a_{2i}&a_{3i}\\b_{1i}&b_{1i}&b_{1i}\end{array}}\right].\left\{{\begin{array}{c}{\xi }^{2}\\\xi \eta \\{\eta }^{2}\end{array}}\right\},{\mbox{ }}i=1,2,3.$

( 19)

$a_{1i}=-{\frac {1}{2}}s_{i}c_{i}^{2}{\mbox{   }}{\mbox{ }}b_{1i}=$$-s_{i}^{2}c_{i}-{\frac {1}{2}}c_{i}^{3}$ $a_{2i}=c_{i}^{3}{\mbox{   }}{\mbox{ }}b_{2i}=-s_{i}^{3}$ $a_{3i}={\frac {1}{2}}s_{i}^{3}+s_{i}c_{i}^{2}{\mbox{   }}{\mbox{ }}b_{3i}=$${\frac {1}{2}}s_{i}^{2}c_{i}$

( 20)

$s_{i}=sen{\varphi }_{i}={\frac {\left({\eta }_{m}-{\eta }_{i}\right)}{{\mu }_{i}}}=$$-{\frac {3{\eta }_{i}}{2{\mu }_{i}}}$ $c_{i}=sen{\varphi }_{i}={\frac {\left({\xi }_{m}-{\xi }_{i}\right)}{{\mu }_{i}}}=$$-{\frac {3{\xi }_{i}}{2{\mu }_{i}}}$

( 21)

onde o índice m indica o ponto médio entre j e k com ξ_m , η_m  e μ_i são as coordenadas adimensionais desse ponto, definidas pelas expressões:

${\xi }_{m}={\frac {1}{2}}({\xi }_{j}+{\xi }_{k})$ ${\eta }_{m}={\frac {1}{2}}({\eta }_{j}+{\eta }_{k})$ ${\mu }_{i}={\sqrt {{\left({\xi }_{m}-{\xi }_{i}\right)}^{2}+{\left({\eta }_{m}-{\eta }_{i}\right)}^{2}}}=$${\frac {3}{2}}{\sqrt {{\xi }_{i}^{2}+{\eta }_{i}^{2}}}$

( 22)

onde i, j e k representam a permutação cíclica entre 1,2 e 3.

A matriz de alta ordem pode ser multiplicada por um escalar positivo β , logo a matriz de rigidez total do elemento é dada por:

$K=K_{b}+\beta .K_{s}$

( 23)

que satisfaz ao patch‐test e a condição de ser positiva definida para a convergência, conforme provado por [20] . O parâmetro β  ≥ 0 para a matriz de alta ordem é utilizado para calibrar a matriz de rigidez total, K , e com isso aperfeiçoa a convergência. Os valores ótimos sugeridos por [20] para α e β são, respectivamente, 1,5 e 0,5.

7. Exemplos numéricos

De posse da formulação descrita nos itens anteriores, desenvolveu‐se um sistema computacional aplicando o SESO em conjunto com o método dos elementos finitos, usando a formulação livre descrita no item 6 baseada em [20] . Introduziu‐se um algoritmo de discretização automática para malhas não estruturadas baseada na triangularização de Delaunay para agilizar o processo e permitir geometria mais complexas, como a formação de poligonais fechadas não convexas e multiplamente conexas. Além disso, devido à necessidade de resolver o sistema linear centenas de vezes no processo iterativo, programou‐se um conjunto de sub‐rotinas em linguagem Fortran otimizadas para resolver sistemas lineares esparsos e simétricos, indicado por Duff e Reid [21] . Adicionalmente, conjugou‐se a metodologia de cálculo do peso próprio gerando automaticamente uma força concentrada equivalente para cada nó do elemento, já computando seu volume e peso específico. Assim, apresentam‐se resultados numéricos para a avaliação e comparação das configurações obtidas pelos modelos clássicos de bielas e tirantes. Os parâmetros de otimização RR e RE, quando não citados, são iguais a 1%, e a função reguladora é definida como η (Γ ) = 10⁻⁴ .

7.1. Exemplo 1

Uma das questões mais intrigantes na natureza ou em estruturas naturais é: «Porque elas são do jeito que são?». Para revelar a resposta desta pergunta, a técnica de otimização pode ser aplicada ao projeto de uma fruta. Especificamente, a intenção é verificar se o algoritmo SESO poderia produzir a mesma forma criada pela natureza de uma maçã, verificando assim se seu formato foi determinado por uma perspectiva genética ou estrutural.

A geometria e as condições de contorno são indicadas na figura 5 , a espessura considerada para chapa foi de 0,44 mm. O módulo de Young é E  = 70 × 10³MPa e coeficiente de Poisson é ν  = 0,3. Utilizou‐se uma razão de rejeição igual a 0%, uma razão evolucionária de 0,5%, e o volume final obtido, com o controle pelo IP, fora igual a 78% do volume inicial e o volume retirado por iteração igual a 5%.

Figura 5. Domínio de projeto.

Uma malha refinada de (44 x 44) totalizando 3.828 elementos finitos triangulares foi utilizada no processo evolucionário. A estrutura está suspensa por um talo em cujos nós extremos permite‐se deslocamento horizontal e o nó central está fixo. A figura 6 ilustra a topologia ótima obtida com a presente formulação, destacando‐se em azul as tensões positivas como devidas à tração que seu peso próprio sobre a estrutura.

Figura 6. Forma semelhante à de uma maçã – topologia ótima.

Este exemplo foi proposto por Querin [18] e Xie et al. [22] , usando ESO. Sua forma ótima está ilustrada na figura 7 . A diferença básica entre as formulações utilizadas por estes autores para otimizar a forma desta estrutura está no critério de rejeição, o primeiro usa tensão de Von Mises e o segundo usa a tensão principal ou a soma das tensões principais como critério, obtendo a forma com a uniformidade das tensões na superfície da estrutura.

Figura 7. Topologia semelhante a uma maçã –proposta de Querin [18] .

O resultado obtido com a presente formulação mostra que a ação do peso próprio traciona a estrutura fazendo com os elementos próximos ao seu centro de gravidade sejam mais tensionados que os elementos mais próximos ao contorno da estrutura. Nota‐se que o comportamento estrutural é assimétrico, devido a discretização realizada de forma automática por um algoritmo baseado na triangularização de Delaunay, que levou a pequenas assimetrias, mas que não foram significativas para a configuração ótima do problema. Na interface entre os comprimidos e tracionados, respectivamente, vermelho e azul é dentada e descontínua em alguns pontos, características da heurística de remoção dos elementos da malha no procedimento de otimização topológica via SESO.

7.2. Exemplo 2

Este exemplo trata de uma viga baixa biapoiada ilustrada na figura 8 a. Os dados utilizados são espessura de 0,50 m, módulo de Young E  = 20820 MPa e coeficiente de Poisson ν  = 0,3. Aplica‐se a técnica SESO e o critério de parada é realizado em termos do IP da estrutura ou quando o volume desejado é alcançado. Consideraram-se 2 tipos de carregamento, o seu peso próprio e outro de 200 kN uniformemente distribuídos em 0,5 m, atuando simetricamente no topo da viga, como esquematicamente indicado na figura 8 a. Nessa viga baixa utilizou‐se uma razão de rejeição igual a 1%, uma razão evolucionária de 2%. O volume desejado foi de igual 40% do volume inicial e o volume retirado por iteração igual a 1,75%. Durante o procedimento evolucionário o domínio fixo estendido foi definido por uma malha refinada de (160 x 60) totalizando 19.200 elementos finitos triangulares.

Figura 8. a) Viga parede baixa – domínio de projeto; b) Topologia ótima não considerando o peso próprio; c) Topologia ótima considerando o peso próprio.

As configurações ótimas obtidas neste problema foram simuladas com as 2 cargas simultaneamente aplicadas durante o processo de otimização. Assim, o algoritmo é capaz de simular melhor o problema físico real. As figuras 8 b e c mostram a influência do peso próprio para o modelo de bielas e tirantes, neste problema viga baixa.

Por se tratar de um problema de otimização não convexo, isto é, com vários mínimos locais, optou‐se por uma configuração ótima de projeto com volume de 46,1%, que é um ótimo comum às 2 funções objetivo do problema, destacando‐se que o processo ótimo com o peso próprio chegou a um volume final desejado de 33,4%.

O esquema indicado pela figura 9 mostra os elementos do modelo biela‐tirante elaborados segundo as topologias ótimas de maneira direta. Observa‐se que na proposta esquematizada as semirretas em vermelho e azul indicam, respectivamente, os membros comprimidos, bielas, e os tracionados, tirantes.

Figura 9. Proposta de modelo de biela‐tirante via OT (sem influência do peso próprio).

Como o peso próprio é pequeno quando comparado com a carga total distribuída, isto é, o seu peso próprio equivale a 5% da carga total distribuída, então, a proposta para o modelo de biela e tirantes com a ação do peso próprio não tem diferenças perceptíveis do modelo indicado na figura 9 .

Considerando a teoria de treliça clássica aplicada ao concreto, o modelo demonstra‐se consistente, com o arco superior comprimido (C1, C2 e C5) e bielas inclinadas (C3 e C4) compondo a treliça com os tirantes (T1 e T3) inclinados na ordem de 30 graus com a horizontal e o banzo inferior tracionado (T2). Estes tirantes eram representados por armaduras longitudinais denominadas «cavaletes» mas, hoje em dia, eles são representados pelo conjunto de armaduras longitudinais e estribos verticais.

Indicando‐se os membros onde os fluxos de tensões se destacam, podem‐se calcular os esforços resultantes em cada um por meio das tensões médias calculadas no processo. A força resultante em cada membro é obtida pela multiplicação da tensão média pela respectiva área transversal de cada membro. Já a área da seção do membro é obtida pelo produto da espessura da viga (dado do problema) com a largura média da região do fluxo, medidas na direção transversal à direção das tensões principais e calculadas automaticamente pelo programa computacional.

Calculam‐se assim as áreas das armaduras necessárias na região de tirantes (tabela 1 ).

Tabela 1. Dimensionamento das armaduras dos tirantes

Sem consideração do peso próprio			Com consideração do peso próprio
Tirante	Força (kN)	AS,rec (cm2 )	Tirante	Força (kN)	AS,rec (cm2 )
T1	28,18	0,645	T1	28,47	0,655
T2	58,52	1,346	T2	57,26	1,317
T3	26,78	0,616	T3	29,57	0,650

Na tabela 2 apresentam‐se as tensões atuantes no concreto dos membros de bielas. Nota‐se um acréscimo mais significativo de esforços no arco comprimido (banzo superior) de 10‐20%, porém sem alterações importantes nos tirantes.

Tabela 2. Comparação das tensões atuantes nas bielas

Sem consideração do peso próprio		Com consideração do peso próprio
Biela	Força (kN)	Biela	Força (kN)
C1	231,01	C1	252,9
C2	231,01	C2	252,9
C3	25,55	C3	15,1
C4	25,55	C4	15,1
C5	231,42	C5	285,5

A figura 10 mostra o processo evolucionário e desempenho da estrutura durante este processo. Destaca‐se que o IP máximo obtido com a consideração do peso próprio é de 1,81 e sem o peso próprio é de 1,70. Mostrando que o aumento das tensões médias nos centroides dos elementos finitos proporciona um aumento no número de elementos retirados da estrutura determinando um menor custo computacional.

Figura 10. IP por número de iterações.

No procedimento evolucionário considerando o peso próprio é relevante destacar que, mesmo sendo a sua ordem de grandeza menor que o carregamento aplicado à viga, as configurações dos esforços nos elementos do modelo não sofreram alterações significativas, no entanto, a sua consideração proporcionaram maior estabilidade e rapidez ao processo evolucionário.

O procedimento de otimização considerando o peso próprio teve o critério por volume desejado enquanto o IP foi o critério de parada para a viga sem consideração do mesmo. Destaca‐se que o volume desejado no processo evolucionário para ambos os casos foi de 46,6%.

7.3. Exemplo 3

Este exemplo foi apresentado por [13] que é uma viga biapoiada com uma abertura, onde a geometria e a carga de cálculo (F_d ) são apresentadas na figura 11 , a qual é usada como domínio estendido para o processo de otimização. Os parâmetros dos materiais utilizados são o módulo de Young E  = 20820 MPa, coeficiente de Poisson v  = 0,15 e espessura igual a 0,4 m.

Figura 11. Domínio estendido do projeto (mm), [11] .

As resistências de cálculo do aço e do concreto são tomadas com valores, respectivamente, de f_yd  = 434 MPa e f_cd  = 25 MPa . Essa estrutura possui uma região «D» em virtude da descontinuidade geométrica devido à presença de uma abertura imposta por projeto. Nesse caso, essa região deve ser avaliada usando um modelo de bielas e tirantes.

Para a modelagem com a presente formulação, utilizou‐se uma malha de elementos finitos com 13.200 elementos triangulares (150 x 47) e na figura 12 são mostrados os projetos ótimos obtidos pela presente formulação SESO e comparados com os descritos nos trabalhos de [11] e [13] . A topologia ótima do projeto, indicada na figura 12 a, não considera o peso próprio da estrutura.

Figura 12. Topologia ótima obtida: a) sem peso próprio, b) peso próprio aplicado no topo, c) peso próprio aplicado em todo corpo da viga (presente modelo), d) [11] , e) [13] .

As figuras 12 b e c mostram as configurações ótimas considerando o peso próprio da estrutura em 2 casos diferentes: o primeiro considera uma F_R  = 3, 33MN que é obtida por F_R  = F_d  + PP aplicada no topo da viga, como uma força concentrada, na posição indicada na figura 11 , e o segundo, que considera o peso próprio distribuído automaticamente para todos os elementos da estrutura aplicados conjuntamente com a força concentrada F_d .

Neste problema de viga‐parede alta os parâmetros de otimização permaneceram os mesmos para os 3 casos avaliados e a configuração ótima obtida com volume igual a 45,6% sem o peso próprio e, quando se considera peso próprio no topo e distribuído no volume da viga, iguais a 35,1 e 29,9%, respectivamente. Nestas figuras as regiões em vermelho indicam as regiões comprimidas, bielas e as tracionadas, tirantes, em azul. Observa‐se que, mesmo o peso próprio sendo insignificante em relação ao carregamento aplicado à viga, há uma melhor definição dos elementos componentes do modelo e sua configuração.

O gráfico da figura 13 mostra as tensões de Von Mises máxima em cada iteração para os 3 casos estudados neste artigo. As tensões médias quando se considera o peso próprio distribuído na viga foram atenuadas representando melhor o modelo real deste problema uma vez que a função objetivo fica mais bem representada neste caso.

Figura 13. Tensão de Von Mises por número de iterações para o presente procedimento de otimização.

O IP atinge um máximo mais rapidamente quando se considera o peso próprio porque é inversamente proporcional ao número de elementos que permanecem na estrutura. Assim, quanto mais elementos são retirados da estrutura pelo critério de otimalidade imposto no processo de otimização maior será o IP.

Schlaich et al. [6] apresentam uma proposta de reforço para essa estrutura, denominada de «Deep Beam with a large hole», obtida com o emprego de um modelo de bielas e tirantes advindo da combinação do método dos elementos finitos com um procedimento de se obter os fluxos de forças, usando o método denominado de «load path». De forma que na figura 17 é apresentada a disposição das armaduras propostas pelos autores para reforço da viga com cavidade, entretanto eles não consideram o seu peso próprio.

Com a indicação da topologia obtida pela presente formulação é possível definir uma proposta de modelo de biela‐tirante de maneira direta e avaliar a influência de seu peso próprio. Observe‐se a proposta esquematizada na figura 14 onde as semirretas em vermelho e azul indicam, respectivamente, os membros comprimidos, bielas, e os tracionados, tirantes. Indicando‐se os membros onde os fluxos de tensões se destacam, obtendo assim os valores das tensões médias de cada membro e com sua respectiva área, dada pelo produto da espessura da viga com a largura média da região do fluxo.

Figura 14. Proposta de modelo de biela‐tirante via OT.

Deste modo se calculam as áreas das armaduras necessárias na região de tirantes e verifica‐se a resistência do concreto em cada biela.

Na tabela 3 são apresentados os valores das forças médias obtidas nos tirantes. Os tirantes T2 e T4 estão inclinados a 15 e 45 graus, respectivamente, da horizontal.

Tabela 3. Dimensionamento das armaduras dos tirantes sem consideração do peso próprio

Tirante	Força (MN)	AS,rec (cm2 )	Força na armadura (MN)	Armadura proposta	Número de barras
T1	1,50		1,50	A = 34,60	2 × 5 ∅ 20
T2		51,53	0,966 T2	Ac = 49,75	2 × 7 ∅ 20
T4			0,01	A	2 × 2 ∅ 20
T4	2,20	50,70	0,707T4	Ac = 35,85	2 × 5 ∅ 20
			0,707T4	Ad = 35,85	2 × 5 ∅ 20

O tirante T2 é representado pela parcela de força decomposta para as barras longitudinais, obtendo‐se assim a área necessária de armadura As₂ .

O tirante T4 tem sua representação na malha ortogonal, As₃ , que cobre o trecho junto à borda da cavidade e os encontros de bielas C3 e C8 na extremidade esquerda, e na armadura complementar As₄ , a 45 graus, cobrindo os encontros de bielas C2, C5 e C8.

Conforme os procedimentos de cálculo admitidos para as armaduras representativas dos tirantes (tabela 3 e fig. 15 ) apresenta‐se o detalhamento destas armaduras na figura 16 , onde nota‐se a extensão das armaduras As₃ procurando abranger áreas de encontros de bielas sobre a abertura.

Figura 15. Representação esquemática da armadura e do modelo de biela‐tirante advindo da OT.

Figura 16. Proposta de disposição de armaduras pelo presente modelo.

Na tabela 4 apresenta‐se a verificação das tensões atuantes no concreto dos membros de bielas onde nota‐se a não observância do estado último de compressão do concreto nas bielas C3 e C8 (Tabela 5  ;  Tabela 6 .).

Tabela 4. Verificação das tensões atuantes nas bielas

Biela	Força (MN)	σ (MPa)	0,8 fcd  = 20 (MPa)
C3	3,3	–28	Reforço
C5	–1,4	–20	Ok
C6	–1,7	–18,5	Ok
C7	–5,5	–18,5	Ok
C5	–2,6	–22	Reforço

Tabela 5. Dimensionamento das armaduras dos tirantes, considerando o peso próprio distribuído elemento a elemento na viga

Tirante	Força (MN)	AS,rec (cm2 )
T1	1,889	80,3943
T2	1,904	59,0522
T4	1,980	61,3555

Tabela 6. Verificação das tensões atuantes nas bielas, considerando o peso próprio distribuído elemento a elemento na viga

Biela	Força (MN)	AS,rec (cm2 )
C1	–2,640	26,6640
C	–1,190	95,9525
C	–1,394	82,6417
C	–4,485	397,5992
C	–2,080

A proposta para o modelo de bielas e tirantes considerando o peso próprio (de pequena grandeza) da estrutura apresenta configuração semelhante a que não se considera o peso próprio. Nos cálculos para o dimensionamento das armaduras, não houve resultados diferentes das forças resultantes obtidas das tensões médias dos elementos tracionados.

Figura 17. Proposta de disposição de armaduras de para a viga com abertura [6] .

8. Conclusão

O objetivo do artigo é a verificação de um sistema robusto, confiável e automatizado para a obtenção de modelos de bielas e tirantes considerando a influência do peso próprio da estrutura.

Assim, empregou‐se um procedimento de otimização topológica denominado SESO para este fim. Ele, em síntese, emprega o método dos elementos finitos para efetuar a presente análise elástica, de modo a se definir um domínio inicial estendido e, iterativamente, o método busca uma configuração topológica ótima em que, naturalmente, ficam definidos os membros componentes do modelo, bielas e tirantes. Dessa forma, os esforços nos membros podem ser avaliados para possibilitar o dimensionamento e os reforços estruturais necessários para cada trecho.

Os 3 exemplos apresentados demonstraram boa acurácia com as teorias clássicas e os valores descritos por outros autores. Foram propostas também a quantificação e a disposição de armaduras para um exemplo clássico descrito na literatura internacional sobre o assunto e um modelo de bielas e tirantes para uma viga parede baixa com vários casos de carregamento.

Constatou‐se que peso próprio da estrutura influencia positivamente o processo evolucionário e a configuração ótima obtida, mesmo com valores de pequena ordem de grandeza comparados aos carregamentos aplicados aos elementos estruturais. Isto se deve ao fato de que o elemento estrutural passa a ter uma distribuição do peso próprio em cada célula de seu volume proporcionando uma alteração nas tensões internas e no processo suavizado de eliminação que, no entanto, proporciona maior estabilidade e rapidez ao processo evolucionário.

Agradecimentos

Os autores agradecem ao Departamento de Engenharia de Estruturas e Geotécnica da Escola Politécnica da Universidade de São Paulo (EPUSP), à Universidade Federal de Ouro Preto (UFOP) e à Universidade Estadual Paulista (UNESP) pelo suporte financeiro desta pesquisa.

Referências

1. [1] G.I.N. Rozvany, M.P. Bendsøe, U. Kirsch; Layout optimization of structures; Applied Mechanics Review, 48 (1995), pp. 41–119
2. [2] F.V. Stump; Otimização Topológica Aplicada ao Projeto de Estruturas Tradicionais e Estruturas com Gradação Funcional sujeitas a Restrição de Tensão; Dissertação (Mestrado) ‐ Escola Politécnica da USP, São Paulo (2006)
3. [3] M.P. Bendsøe; Optimal Shape design as a material distribution problem; Structural Optimization, 1 (1989), pp. 193–202
4. [4] G.I.N. Rozvany, M. Zhou, T. Birker; Generalized shape optimization without homogenization; Structural Optimization, 4 (1992), pp. 250–252
5. [5] A. Labanowski, E.A. Fancello, A.A. Novotny; SIMP ESO e TSA: uma análise comparativa de métodos de otimização topológica para elasticidade 2d; CILAMCE Recife; 10/11/2004; 12/11/2004; Publicação: Proceedings of the XXV Iberian Latin‐American Congress on Computational Methods in Engineering; CILAMCE (2004)
6. [6] J. Schlaich, K. Schafer, M. Jennewein; Toward a consistent design of structural concrete; PCI‐Journal, 32 (3) (1987), pp. 74–150 May/June
7. [7] M. Ali; Automatic generation of truss model for the optimal design of reinforced concrete structures. Dissertation; Cornell University, Ithaca, NY, USA (1997)
8. [8] Q.Q. Liang, Y.M. Xie, G.P. Steven; On equivalent between stress criterion and stiffness criterion in evolutionary structural optimization; Structural Optimization, 18, Springer – Verlag (1999) 67‐13
9. [9] Q.Q. Liang, Y.M. Xie, G.P. Steven; “Topology optimization of strut‐and‐tie models in reinforced concrete structures using an evolutionary procedure”; ACI Structural Journal, 97 (2) (2000), pp. 322–330
10. [10] Q.Q. Liang, B. Uy, G.P. Steven; Performance‐Based Optimization for Strut‐Tie Modeling of Structural Concrete; Journal of Structural Engineering (2002), pp. 815p–823p
11. [11] 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 (13‐14) (2002), pp. 1471–1489
12. [12] Q.Q. Liang; Performance‐based Optimization of Structures: Theory and applications; Spon Press, London (2005)
13. [13] M. Brugge; Generating strut‐and‐tie patterns for reinforced concrete structures using topology optimization; Computers and Structures, 87 (2009), pp. 1483–1495
14. [14] H.L. Simonetti, V.S. Almeida, F.A. Neves; The influence of geometry on the structural elements of topology via SESO; CILAMCE, Armação de Búzios (2009)
15. [15] H.L. Simonetti, V.S. Almeida, F.A. Neves; Seleção de Topologias Ótimas para Estruturas do Contínuo com minimização de Volume sujeita a restrição de tensão via “Smoothing ESO” (SESO); CILAMCE, Argentina (2010)
16. [16] Y.M. Xie, G.P. Steven; A simple evolutionary procedure for structural optimization; Computers & Structures, 49 (5) (1993), pp. 885–896
17. [17] P. Tanskanen; The evolutionary structural optimization method: Theoretical aspects; Comput Methods Appl Mech Eng, 191 (2002), pp. 5485–5498
18. [18] O.M. Querin; Evolutionary Structural Optimization stress based formulation and implementation; PhD dissertation, University of Sydney (1997)
19. [19] M.R.C. Leggerini, Mecânica dos sólidos; Departamento de Engenharia Civil; PUC‐RJ (2005)
20. [20] P.G. Bergan, C.A. Felippa; A triangular membrane element with rotational degrees of freedom; Comp. Meths. in Appl. Mech. Eng., 50 (1985), pp. 25–69
21. [21] I.S. Duff, J.K. Reid; A set of Fortran subroutines for solving sparse symmetric sets of linear equations; Report R10533, AERE Harwell (1982)
22. [22] Y.M. Xie, X. Huang, J.W. Tang, P. Felicetti; Recent Advances in Evolutionary Structural Optimization; Keynot Lecture for Frontiers of Computational Science Symposium, Nagoya University, Japan (2005)