UMA TÉCNICA PARA OTIMIZAÇÃO ESTRUTURAL MEDIANTE … · Prof. Luiz Eloy Vaz, D.Sc. RIO DE JANEIRO, RJ ... dos elementos finitos é utilizado para discretizar o domínio e obter

UMA TÉCNICA PARA OTIMIZAÇÃO ESTRUTURAL

MEDIANTE A DERIVADA TOPOLÓGICA

Marcelo de França Cordeiro

DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DOS

PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA UNIVERSIDADE

FEDERAL DO RIO DE JANEIRO COMO PARTE DOS REQUISITOS

NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE EM CIÊNCIAS EM

ENGENHARIA MECÂNICA.

Aprovada por:

________________________________________________

Prof. Jose Herskovits Norman, D.Ing.

________________________________________________

Prof. Fernando Pereira Duda, D.Sc.

________________________________________________

Prof. Luiz Eloy Vaz, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

MARÇO DE 2007

ii

CORDEIRO, MARCELO DE FRANÇA

Uma Técnica para Otimização Estrutural

Mediante a Derivada Topológica [Rio de

Janeiro] 2007

XII, 86 p. 29,7 cm (COPPE/UFRJ,

M.Sc., Engenharia Mecânica, 2007)

Dissertação - Universidade Federal do

Rio de Janeiro, COPPE

1. Análise de Sensibilidade Topológica

2. Método dos Elementos Finitos

I. COPPE/UFRJ II. Título (série)

iii

Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos

necessários para a obtenção do grau de Mestre em Ciências (M.Sc.)

UMA TÉCNICA PARA OTIMIZAÇÃO ESTRUTURAL

MEDIANTE A DERIVADA TOPOLÓGICA


Março/2007

Orientador: Jose Herskovits Norman

Programa: Engenharia Mecânica

Esta dissertação aplica o conceito da análise de sensibilidade topológica e da

otimização estrutural evolutiva ao problema da otimização topológica de sólidos.

Estuda-se o caso de minimizar a energia interna de deformação com restrição nas

tensões. É considerado o comportamento do sólido como elástico e linear. O método

dos elementos finitos é utilizado para discretizar o domínio e obter uma solução

aproximada do problema. Os resultados numéricos obtidos são comparados com

resultados existentes na literatura. Por fim, é apresentada uma comparação entre a

presente formulação com formulações diferentes, provando a similaridade entre elas.

iv

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Master of Science (M. Sc.)

A STRUCTURAL OPTIMIZATION TECHNIQUE

USING THE TOPOLOGY SENSITIVITY ANALYSIS


March/2007

Advisor: Jose Herskovits Norman

Department: Mechanical Engineering

This dissertation applies the concepts of topological sensitivity analysis and

evotutionary structural optimization in the problem of topology optimization of solids.

The case studied is to minimize the internal deformation energy with stress constraint.

The comportment of the solid is linear and elastic. The method of finite element is used

to discretize the domain and to obtain an approximate solution of the problem. The

numerical results are compared with existing results in literature. Finally, is presented a

parallel between the present formulation with different ones proving their similarity.

v

Aos meus pais Fernando e Joana, meu irmão Alexandre e Tarsila.

vi

A ciência está em constante evolução. A necessidade do saber é um instinto que

está presente no inconsciente do Homem. Do fogo ao plasma, da célula ao quark. A

ciência empenha papel fundamental na quebra dos paradigmas da Sociedade.

vii

Sumário

1 Introdução 1

1.1 Conceitos Iniciais ..................................................................................1

1.2 Otimização Estrutural ............................................................................3

1.3 Derivada Topológica .............................................................................5

1.4 Histórico Evolutivo da Otimização Topológica ....................................6

1.5 Exemplos de Aplicações........................................................................8

2 Elasticidade Linear de Sólidos 9

2.1 Problema Variacional e de Valor de Contorno da Elasto-Estática

Linear ...............................................................................................................9

3 Derivada Topológica 12

3.1 Formulação da Derivada Topológica ..................................................12

3.2 Parametrização ....................................................................................15

3.3 Análise de Sensibilidade à Mudança de Forma...................................18

3.4 Método de Lagrange............................................................................19

3.5 Cálculo da Derivada Topológica .........................................................21

4 Otimização Topológica de Sólidos 30

4.1 Formulação do Problema de Minimizar a Energia de Deformação com

Restrições nas Tensões ................................................................................................30

4.2 Formulação da Distribuição de Material .............................................31

4.3 Discretização do Domínio ...................................................................32

4.4 Aplicação da Derivada Topológica no Problema de Otimização........33

5 Alivio de Tensões com a Otimização Estrutural

Evolucionária 37

viii

5.1 Formulação de Alívio de Tensões .......................................................38

6 Comparação dos Métodos de Otimização Topológica 42

6.1 Aproximação por Microestrutura ........................................................42

6.2 Aproximação por Macroestrutura........................................................43

6.3 Comparação entre Métodos .................................................................44

6.4 Tensão de Mises x Energia de Deformação Interna ............................45

7 Método dos Elementos finitos 47

7.1 Visão Geral ..........................................................................................47

7.2 Discretização das Equações.................................................................48

7.3 Treliça Tridimensional ........................................................................51

7.4 Triângulo Linear ..................................................................................53

7.5 Quadrado Bi-linear de Malha Regular.................................................55

8 Métodos de Resolução de Sistemas Lineares 58

8.1 Esparsidade..........................................................................................58

8.2 Matriz Simétrica ..................................................................................59

8.3 Matriz Definida Positiva......................................................................60

8.4 Matriz em Banda .................................................................................60

8.5 Decomposição de Cholesky ................................................................62

9 Exemplos 63

9.1 Viga Bi-apoiada com Carga Distribuída na Parte Superior.................64

9.2 Viga Engastada Superiormente com Carga Distribuída Cisalhante ....68

9.3 Viga Bi-apoiada com Carga Distribuída Inferiormente ......................72

10 Implementação Computacional 76

10.1 Programação Otimizada ..................................................................76

10.2 Estrutura do Programa.....................................................................77

ix

10.3 Pré-processamento...........................................................................78

10.4 Processamento .................................................................................78

10.5 Pós-processamento ..........................................................................79

11 Conclusão 80

12 Bibliografia 82

x

Lista de Símbolos

Ω1 Vetor de índices que indica quais pontos do domínio estão ativos. A Área da seção transversal da treliça. a(.,.) ou aτ Forma bi linear contínua e coerciva. b Conjunto de forças de corpo do domínio. Be Esfera de raio e. B Matriz de incidências. D Tensor constitutivo do material, ou tensor de elasticidade. D0 Tensor constitutivo do material na configuração inicial. Dad Conjunto de tensores de elasticidade admissíveis para o resultado

ótimo. D*

Tensor constitutivo do material penalizado. DT(x) ou DT Função vetorial que define a derivada topológica em relação a um

vetor x. E Campo de deformações. E Módulo de Elasticidade, Módulo de Young. e Raio de um furo B. En Energia de deformação interna. f(x) ou f Função objetivo ou função custo. fR(e) Função escalar regularizadora em função do raio e. g(x) ou g Vetor de Restrições de Desigualdade. h(x) ou h Vetor de Restrições de Igualdade. H1 Espaço de Hilbert de ordem 1. I Matriz identidade. ki Matriz de rigidez do elemento. K Matriz de rigidez global. L Comprimento da treliça. l(.) ou lτ Funcional linear contínuo. n Vetor normal a uma superfície pertencente ao domínio. Ni Função de forma de um nó i pertencente a um elemento. N Conjunto dos número naturais. N Dimensão do problema. p Expoente de Penalização na formulação SIMP. q Conjunto de Tensões aplicadas em um domínio. q Número de variáveis do problema. r Número de restrições de desigualdade do problema. Rn Conjunto dos números reais, de dimensão n.

+R Conjunto dos números reais positivos. s Número de restrições de igualdade do problema. t Espessura da malha de elementos finitos. T Campo de tensões. Te Campo de tensões no domínio perturbado. Tτ Campo de tensões parametrizado no domínio perturbado. TVM Tensão de Von Mises. u Conjunto de pontos de engaste do domínio. u Campo de deslocamentos.

xi

ue Campo de deslocamentos no domínio perturbado. uτ Campo de deslocamentos no domínio perturbado e parametrizado,

função implícita de τ. U Conjunto das funções e variações admissíveis do deslocamento u.

eU Conjunto das funções e variações admissíveis do deslocamento ue.

τU Conjunto das funções e variações admissíveis do deslocamento parametrizado uτ.

V Conjunto das funções e variações admissíveis do deslocamento η. eV Conjunto das funções e variações admissíveis do deslocamento ηe.

τV Conjunto das funções e variações admissíveis do deslocamento parametrizado ητ.

ΩV Volume relativo ao domínio Ω. VΓ Volume dos elementos que pertencem a fronteira Γ.

maxV Volume máximo admitido para a solução ótima. v(x) ou v Ação de mudança de forma, ou descrição material de velocidade

que caracteriza a mudança de forma de um corpo. vn Componente normal da Ação de mudança de forma v(x). V Matriz de coeficientes de Von Mises. x Vetor das variáveis de projeto. x Ponto interior ao domínio. x Vetor de solução ótima do projeto. xτ Nomenclatura resumida para ),( τχ x , função que define uma

perturbação parametrizada em um domínio.

x& Derivada parcial de x em relação ao parâmetro τ, xx&=

∂∂τ

.

Xi Vetor que armazena a coordenada nodal de um elemento i Z Matriz de coeficientes de Von Mises com o tensor constitutivo

agregado. α1,...,n Constantes de interpolação do elemento finito.

eB∂ Fronteira do furo Be.

supξ Parâmetro superior para a eliminação de elementos utilizando a Derivada Topológica.

infξ Parâmetro inferior para a eliminação de elementos utilizando a Derivada Topológica.

φ Função regular definida em um domínio não parametrizado.

τφ Função regular parametrizada.

τφΩ Função regular definida em τΩ .

τφΓ Função regular definidas em τΓ .

iφ Função de interpolação dos elementos finitos. μτ Deslocamento arbitrário pertencente ao conjunto τV . η Deslocamento arbitrário pertencente ao conjunto V . ητ Deslocamento arbitrário pertencente ao conjunto τV . ηe Deslocamento arbitrário do domínio perturbado pertencente ao

conjunto eV .

xii

ητ Deslocamento arbitrário do domínio perturbado e parametrizado pertencente ao conjunto τV .

Ω Domínio de projeto, ou de material. Ω0 Domínio inicial de projeto, ou de material. Ωe Domínio de projeto perturbado com um furo Be de raio e. Ωe+δe Domínio de projeto perturbado com a expansão de um furo Be de

raio e para um de raio e+δe. Ωτ Domínio de projeto perturbado e parametrizado. Ωelem Domínio de um elemento pertencente à malha de elementos finitos. τ Parâmetro que define a perturbação no domínio. Г Contorno ou fronteira do domínio. ГD Região da fronteira onde há engastes. Гe Fronteira perturbada com um furo Be de raio e. Гe+δe Fronteira perturbada com a expansão de um furo Be de raio e para

um de raio e+δe. ГN Região da fronteira onde são aplicadas tensões externas. Гτ Fronteira perturbada e parametrizada.

)( τψ Ω Função custo parametrizada em relação ao domínio parametrizado. )( τuΨ ou Ψ Função custo parametrizada em relação ao campo de

deslocamentos parametrizado. λτ Multiplicador de Lagrange.

),( τχ x Função que define uma perturbação parametrizada em um domínio. ),( ⋅⋅L Função Lagrangeano.

Σe Tensor de Eshelby na configuração perturbada do domínio. ρe Vetor de Densidades do domínio. ν Coeficiente de Poisson. σ1,...,i Tensão principal 1,...,i.

e

eρψ∂∂ Sensibilidade da função custo perturbada em relação à densidade

na formulação SIMP e ESO.

1

1 Introdução

1.1 Conceitos Iniciais

O projeto é uma etapa fundamental no desenvolvimento de um produto. Pode-se

dividir esta etapa em duas fases: concepção e dimensionamento. Na concepção, são

definidos os conceitos básicos e tecnologias empregadas para o desenvolvimento do

projeto, obtendo-se uma definição inicial do produto. Em uma segunda etapa, são

definidas as dimensões, forma e topologia através da utilização de técnicas analíticas ou

numéricas que propiciam a obtenção de produtos que atendam a todas as especificações

definidas na etapa da concepção.

A palavra “otimização” é definida como “fazer tão perfeito, efetivo ou

funcional quanto possível”. A otimização tem como meta dimensionar um produto da

melhor forma possível em relação a um objetivo, de modo que este produto atenda às

condições impostas pela sua concepção. Por exemplo: na indústria aeronáutica o peso é

um fator preponderante para se obter um bom projeto de uma peça, e é necessário

reduzir ao máximo este valor. No entanto, esta peça deve resistir a todos os esforços

nela aplicados. Além do peso, outra variável de suma importância para a indústria é o

custo final do produto. Desta forma, minimizando-se o peso, tem-se como conseqüência

o uso de menos matéria prima para se fabricar uma peça, resultando em redução do

custo. Em uma visão empresarial, cada grama de material economizado pode

2

representar um acréscimo considerável nos lucros para uma empresa, fazendo com que

esta se diferencie em relação às concorrentes. Portanto, projetar sem utilizar a

otimização pode acarretar superdimensionamentos desnecessários ao projeto.

Matematicamente, pode-se apresentar o problema de otimização da seguinte

forma:

problema do restrições0)(0)(

:que tal

objetivo funçãoou custo função)(minimize

⇒⎭⎬⎫

⎩⎨⎧

=≤

⇒

xx

xx

hg

f

(1.1)

Em que srq hgf RRRR ∈∈∈∈ e , , ,x , sendo q, r e s escalares que

definem o número de variáveis, restrições de desigualdade e de igualdade.

O problema mostrado em (1.1) consiste em achar, dentre todos os pontos x que

verificam as restrições, aquele que corresponde ao menor valor da função custo.

Definição: x é vetor das variáveis de projeto. f é a função objetivo ou também

chamada de função custo. g é o vetor de restrições de desigualdade e h é o vetor de

restrições de igualdade. Um vetor x satisfazendo todas as restrições impostas para o

problema é chamado de ponto viável do problema. O conjunto de todos os pontos

viáveis do problema é chamado de região viável. A solução ótima do problema é dada

para qR∈x , quando )()( xx ff ≤ , para todo x que pertencer à região viável.

Nas diversas áreas da engenharia, existem técnicas de dimensionamento que,

embora não forneçam um projeto ótimo, resultam em projetos de boa qualidade. Nesta

dissertação serão estudadas algumas técnicas com estas características, aplicadas ao

projeto de estruturas.

3

1.2 Otimização Estrutural

A otimização estrutural vem sendo estudada desde os anos 60 e apresenta

resultados satisfatórios na construção civil, indústria automobilística e indústria

aeroespacial. Estes projetos têm como objetivo básico a redução do custo financeiro e

da quantidade de material empregado, e, além disso, a garantia de que a estrutura

suporte todas as restrições mecânicas impostas, sejam elas estáticas ou dinâmicas.

A função custo mais comumente utilizada na engenharia é o peso total da

estrutura, com restrições locais de tensão. No entanto, nesta dissertação será utilizada

como função custo a energia de deformação interna com a restrição nas tensões. As

restrições mais usuais na engenharia são: tensão máxima de Mises, deslocamentos

máximos, volume máximo ou modos de vibração admitidos.

Existem três classes de otimização, e na concepção do projeto deve-se escolher

qual delas utilizar: otimização dimensional, otimização geométrica e otimização

topológica.

1.2.1 Otimização Dimensional

A otimização dimensional utiliza como variável de projeto um parâmetro de

um elemento estrutural. No caso da treliça, considerado um elemento unidimensional, a

variável de projeto é a área da seção transversal. Para um sólido bidimensional, a

variável de projeto é a sua espessura e necessita-se utilizar técnicas de discretização do

domínio para possibilitar a obtenção de uma solução numérica via método dos

elementos finitos

4

1.2.2 Otimização de Forma

A otimização de forma tem como objetivo definir a melhor fronteira de um

sólido com relação a uma função custo e restrições mecânicas de projeto. O domínio

deste problema necessita ser discretizado para existir um número finito de pontos ao

longo da fronteira e assim possibilitar uma solução numérica via método dos elementos

finitos.

1.2.3 Otimização Topológica

A otimização topológica tem o objetivo de definir da melhor forma possível a

distribuição de material em um domínio pré-determinado considerando-se uma função

custo e as restrições mecânicas. A variável de projeto, neste caso, é o próprio domínio,

no qual se pode eliminar material, ou seja, criar furos.

O domínio do problema da topologia necessita ser discretizado para se

trabalhar apenas com um número finito de pontos para se obter uma solução numérica

aproximada das equações diferenciais parciais do problema. Os elementos pertencentes

a este domínio são interpretados como bits, cujo valor pode ser 0 ou 1, ou seja,

pertencentes ou não ao domínio ótimo. A solução do problema discreto é definida como

a união de todos os elementos que tiverem o valor 1. Nesta dissertação, optou-se por

trabalhar com a formulação discreta das variáveis de projeto porque desta forma se

obtém uma solução final otimizada sem a necessidade da aplicação de filtros e de

penalidades, como, por exemplo, é feita no modelo SIMP (material sólido isotrópico

com penalização). Além disso, utilizando-se a formulação discreta não ocorre o

problema do “tabuleiro de xadrez”.

Como mencionado anteriormente, nesta dissertação será utilizada a energia

interna de deformação como função custo. A restrição mecânica empregada será a

tensão máxima de Mises admitida na estrutura. Apesar de não se trabalhar com o peso

5

como função custo, o resultado final sempre acarretará em uma redução de peso

considerável.

1.3 Derivada Topológica

A derivada topológica é uma função escalar que estuda a sensibilidade da

função custo com relação ao domínio, levando em conta determinadas restrições. A

função custo, relembrando o problema de otimização topológica desta dissertação, será

a energia interna de deformação e a restrição será a tensão máxima de Mises. Em outras

palavras: faz-se um furo infinitesimal no domínio, verifica-se a variação da energia

interna de deformação em relação àquele ponto, e se obtém o valor da derivada

topológica neste ponto.

Considera-se um domínio inicial de material NR⊂Ω , com N sendo um

escalar que representa a dimensão do problema. Um novo domínio perturbado é criado

sendo ee B−Ω=Ω . Be é uma esfera de raio e. O contorno é dado por ee B∂∪Γ=Γ .

Conforme mostrado na Ilustração 1.1.

Em linguagem matemática, pode-se representar a derivada topológica da

seguinte forma:

( ) ( )( )ef

DR

ee

TΩ−Ω

=→

ψψ0

lim)ˆ(x (1.2)

)ˆ(xTD é um escalar em função do ponto x , ψ é a função custo parametrizada

e fR é uma função regularizadora.

Esta é a definição clássica da análise de sensibilidade de forma[1] que será

utilizada para o cálculo da derivada topológica. Desta maneira, podem-se obter os

6

pontos do domínio que são pouco sensíveis à função custo em comparação com outros

pontos, e, portanto, são candidatos a serem eliminados do domínio.

Ilustração 1.1 – Representação da Derivada Topológica.

1.4 Histórico Evolutivo da Otimização Topológica

Focalizando o problema da otimização topológica, será apresentado um sucinto

histórico evolutivo[2] das descobertas e formulações relativas a esta área de pesquisa,

assim como o seu estado da arte.

Em MICHELL[3], no ano de 1904, encontra-se o primeiro registro sobre o

tema da topologia. Michell fundamenta uma teoria sobre a topologia de um domínio

resultando em treliças como resultado ótimo da minimização de peso. Apesar disso, o

embasamento teórico sobre o assunto surgiu algumas décadas depois. CHENG e

OLHOFF[4], em 1981, apresentaram o problema da otimização dimensional de placas e

provaram que a utilização de material compósito aparece naturalmente na formulação

da otimização topológica. O alavancar dos métodos de otimização topológica se deu

com BENDSOE e KIKUCHI[5], em 1988. Neste artigo foi usado o método dos

elementos finitos aliado à teoria da homogenização. Esta teoria resultava em uma bem-

fundada dissertação em que se simulavam furos considerando a rigidez deste ponto do

eΩ xΩ

eβ

Γ

b

qu

e

7

domínio tendendo para zero. Um ano depois, BENDOSE[6], em 1989, formulou a

teoria do Material Sólido Isotrópico com Penalização (SIMP, em Inglês) que realizava a

penalização dos elementos de densidade intermediária para se aproximar da solução

discreta do problema. Estas publicações pioneiras causaram vasto desenvolvimento no

campo da otimização topológica.

Paralelamente ao estudo da otimização topológica, outra área de pesquisa

estava se desenvolvendo: a análise de sensibilidade à mudança de forma. O estudo da

sensibilidade em relação a perturbações na fronteira teve origem nos anos 70 com

MURAT e SIMON[7]. Este trabalho forneceu as bases matemáticas para a análise a

mudança de forma[8]. Em 1994, ESCHENAUER et al.[9] apresentaram uma

formulação que atendia tanto à mudança de forma quanto a topologia ótima usando as

técnicas da análise de sensibilidade à mudança de forma. No entanto, a formulação do

problema continha limitações com relação à função custo. Em 1998, o trabalho de CÉA

et al.[10] calculava a derivada topológica via análise de sensibilidade à mudança de

forma de maneira correta, em comparação com a formulação de NOVOTNY[11].

Entretanto, apenas era válida para a condição de contorno de Newman homogênea nos

furos. Finalmente, em 2003, NOVOTNY[11], em sua tese de doutorado, propõe uma

formulação completa para o problema, provando ser possível utilizar condições de

contorno de Newman, Robin e Dirichilet em diversas situações da engenharia,

estendendo enormemente a gama de problemas nos quais podem ser utilizados os

conceitos da derivada topológica. Atualmente, já existem alguns softwares comerciais

de otimização topológica difundidos pelo mundo. Tem-se como exemplo o NASTRAN

Topological Optimization, um dos mais usados programas comerciais de análise

estrutural, desenvolvido pela National Aeronautics and Space Administration (NASA).

8

1.5 Exemplos de Aplicações

No decorrer de todos estes anos de pesquisa, muitos exemplos práticos da

otimização topológica foram propostos. Como exemplo, cita-se:

1. Condução estacionária de calor em sólidos rígidos, considerando a

função objetivo como sendo a energia potencial total e condições de

contorno de Newman, Dirichilet ou Robin;

2. Flexão elástica linear de placas de Kirchhoff considerando a energia

potencial total como função objetivo;

3. Torsão de barras prismáticas sujeitas à fluência estacionária

considerando a energia de dissipação complementar como função

objetivo;

4. Elasticidade linear de sólidos, com o peso como função objetivo e

restrições locais de tensão;

5. Elasticidade linear de Sólidos, com a Energia de Deformação Interna

como função objetivo, que será abordada nesta dissertação.

9

2 Elasticidade Linear de Sólidos

Para descrever o fenômeno físico da elasticidade linear, e conseqüentemente

otimizar topologicamente o domínio de projeto, é necessário a definição de um modelo

físico de comportamento do sólido, considerando algumas simplificações. Adotar-se-á

neste trabalho um material isotrópico linear sob o estado plano de tensões. As equações

que modelam este fenômeno serão descritas a seguir.

2.1 Problema Variacional e de Valor de Contorno da Elasto-

Estática Linear

Seja um corpo elástico ocupando a região caracterizada por um domínio

regular nR⊂Ω , e limitado pelo contorno Γ, sofrendo tensões externas NΓ∈q e

forças de corpo b , engastado nos pontos DΓ∈u ,e n é um vetor normal à superfície.

Conforme mostrado na Ilustração 2.1.

O problema de valor de contorno consiste em determinar o campo de

deslocamentos u, o campo de deformações E, e o campo de tensões T, existentes no

corpo. Considera-se a hipótese cinemática para pequenas deformações e deslocamentos,

na qual não serão consideradas tensões, deformações iniciais, e força de corpo. Será

adotado um material isotrópico, homogêneo e elástico-linear.

10

Ilustração 2.1 – Exemplo genérico de um problema da mecânica dos sólidos.

As seguintes condições devem ser satisfeitas:

Equação diferencial de equilíbrio:

0=+ bTdiv (2.1)

Equação cinemática:

[ ] ( ) sT uuuE ∇=∇+∇=21 (2.2)

Equação constitutiva:

ET D= (2.3)

D é o tensor constitutivo do material, ou tensor de elasticidade. [E] é o tensor-

deformação. E é o tensor-deformação em forma de vetor, para o caso do estado plano de

tensões. T é o tensor de tensões em forma de vetor para o estado plano de tensões.

As seguintes condições de contorno serão consideradas:

Pontos engastados:

uu = sobre ΓD (2.4)

Pontos de aplicação de tensão na superfície:

qTn = sobre ΓN (2.5)

Resumidamente o problema resulta no seguinte sistema de equações (forma

diferencial ou forte):

Encontre u tal que:

ΓΩ

b

qu

n

11

⎪⎩

⎪⎨

⎧

Γ=Γ=Ω=+

em , em , em , 0

D

N

div

qTnuu

bT (2.6)

Reescrevendo o problema em sua forma fraca, usando o princípio dos trabalhos

virtuais, tem-se:

Encontre U∈u tal que:

V∈∀Γ+=Ω∇ ∫∫∫ ΓΩΩ

ηηqηbsηT ddΩ.d. . (2.7)

O conjunto das funções e das variações admissíveis é definido respectivamente

por [11]:

0))((

))((

Γ21

Γ21

=Ω∈=

=Ω∈=

vη

uuu

H

H

V

U (2.8)

H1 é o espaço de Hilbert de ordem 1.

A equação pode ser representada na forma abstrata a seguir:

Ω),( d.a ∫Ω ∇= sηTηu (2.9)

Γ+= ∫∫ ΓΩ

dd.l ηqηbη .Ω)( (2.10)

Esta equação relaciona os esforços externos com os internos aplicando um

deslocamento arbitrário (virtual) η e fazendo-se um balanço de energia em todo o

domínio Ω.

12

3 Derivada Topológica

3.1 Formulação da Derivada Topológica

Para se obter o valor da derivada topológica, necessita-se especificar o

problema variacional que descreve este fenômeno. As equações que regem seu

comportamento são as mesmas descritas no capítulo 2. Trabalhar-se-á com a forma

fraca destas equações porque isto viabiliza a utilização de espaços topologicamente

mais fracos, e que estabelecem facilmente condições de existência e unicidade da

solução u. O problema é definido da seguinte forma:

Definição: Considera-se um domínio inicial de material Ne R⊂Ω , com um

furo Be de raio e. Expande-se este furo para um raio e+δe. Um novo domínio perturbado

é criado sendo eeeee B δδ ++ −Ω=Ω . O contorno é dado por eeee B δδ ++ ∂∪Γ=Γ . Ilustração

3.1

A equação de Euler-Lagrange, associada ao problema variacional, agora com o

domínio perturbado Ne R⊂Ω , é dado como o seguinte problema de valor de contorno

vetorial de segunda ordem:

Encontre ue tal que:

13

⎪⎪⎩

⎪⎪⎨

⎧

∂=Γ=Γ=Ω=+

em , 0 em ,

em , em , 0

D

B

div

ee

Ne

e

e

nTqnT

uubT

(3.1)

Definindo os domínios inicial e perturbado desta maneira, é possível criar um

homeomorfismo entre eles.

Definição: Dois espaços topológicos dizem-se homeomorfos se existir uma

aplicação entre esses espaços que seja contínua, invertível e a sua inversa seja contínua.

Portanto, como será provado mais adiante, o cálculo da derivada topológica

pode ser obtido usando-se todos os conceitos da análise de sensibilidade à mudança de

forma. É importante salientar que expandir um furo de raio e quando este tende para

zero, nada mais é do que criá-lo.

Ilustração 3.1 – Figura representativa do modelo com os domínios inicial e perturbado.

ee δ+Γee δ+Ω x

eΓeΩ x

ee δβ +

eβ

n n

14

O Domínio inicial eΩ deve satisfazer as condições de contorno estipuladas

pelo problema de valor de contorno. Inclusive o contorno do furo eB∂ , que dependendo

do problema analisado pode conter uma restrição de Newman, Dirichilet ou Robin. No

caso específico da otimização topológica de sólidos, a condição de contorno eB∂ no furo

é a de Newman homogênea. O problema perturbado automaticamente satisfará as

condições do problema inicial, visto que a perturbação eeB δ+ é infinitesimal. O

problema se resume em:

Encontre ee U∈u , tal que:

eeeee la V∈∀= ηηηu )(),( (3.2)

Satisfazendo as condições de contorno eB∂ .

Estando o problema variacional bem definido, define-se a derivada topológica

da seguinte maneira:

( ) ( )( )

( ) ( )( ) ( )

( ) ( )( ) ( )efeef

efeefefD

RR

eee

ee

RR

eeeeeR

ee

T

−+Ω−Ω

=

=⎭⎬⎫

⎩⎨⎧

−+Ω−Ω

=Ω−Ω

=

+

→→

+

→→→

δψψ

δψψψψ

δ

δ

δδ

00

000

lim

limlimlim)ˆ(x

(3.3)

Para usar todo o ferramental da teoria da análise de sensibilidade à mudança de

forma, necessita-se criar um vínculo desta com a derivada topológica. Objetivando isto,

a equação (1.2) será parametrizada, e dessa forma o vínculo será criado.

15

3.2 Parametrização

Para se calcular analiticamente a expressão (3.3), necessita-se fazer a

parametrização do domínio de modo que as condições de homeomorfismo sejam

satisfeitas. Considera-se então uma perturbação dependente do parâmetro τ, denotada

por ),( τχ x com eΩ∈x e R∈τ . Então o domínio perturbado parametrizado Ωτ e o seu

contorno Γτ , são definidos por:

e

N

eeN

eR

eR

Γ=Γ==Γ∈∃∈=Γ

Ω=Ω==Ω∈∃∈=Ω

==

==

00

00

),,(,

),,(,

ττττττττ

τττττττ

τχ

τχ

xxxxxx

xxxxxx (3.4)

Faz-se uma linearização na perturbação ),( τχ x , usando a série de Taylor em

torno de 0τ , e se define:

)())(,(),(),( 000 ττττχτ

τχτχ O+−∂∂

+= xxx (3.5)

Fazendo 00 =τ , e simplificando a notação ),( τχ x por τx , )0,(xχ por x e

)()0,( xvx =∂∂ χτ

, tem-se que todo ponto com τ próximo de zero pode ser escrito por:

)(xvxx ττ += (3.6)

)(xv é uma ação de mudança de forma, ou descrição material de velocidade

que caracteriza a mudança de forma de um corpo.

Definição: A ação de mudança de forma )(xv , como foi provado em

ZOLÉSIO[12], somente tem valor significativo em sua componente normal às fronteiras

Г e B∂ do domínio Ω. Esta componente normal será denominada vn. Esta afirmação

16

baseia-se na idéia de que somente vn produz de maneira efetiva uma mudança de forma

no corpo.

Redefinindo a descrição material de velocidade, e levando em conta que a

perturbação no domínio é devida à expansão de um furo eB∂ , tem-se a definição deste

campo de velocidade:

⎩⎨⎧

∂∪Γ=ΓΓ=∂≤=

e

enn

BB

e que lembrando , sobre 0 sobre 0 com

vvnvv

(3.7)

n é um vetor normal à superfície. Analisando a definição acima, observa-se que

só ocorrerá mudança de forma no ponto onde foi expandido um furo. Nas outras áreas

do contorno a velocidade é nula.

Re-escrevendo (3.6):

en B∂∈∀+= xnvxx ττ (3.8)

Desta maneira, é possível associar o parâmetro τ com a perturbação eδ :

eee BBe δττδ +∂∈∂∈−= xxxx e, com (3.9)

Então:

nne vnv ττδ == (3.10)

Portanto, a relação entre o parâmetro e a perturbação está definida.

Com a equação (3.10), as formulações da análise de sensibilidade à mudança

de forma e a da derivada topológica passam a ser equivalentes, e todo o ferramental

usado para aquela formulação será usado também para a derivada topológica.

17

Este modelo parametrizado de domínio também deve respeitar as equações do

problema de valor de contorno definidas pela equação (3.2).

Encontre uτ τU∈ , tal que:

0,)(),( ≥∀∈∀= ττττττ Vηηηu la (3.11)

τU e τV são os conjuntos das funções e das variações admissíveis.

O conjunto das funções e das variações admissíveis é definido respectivamente

por:

0))((

))((

Γ21

Γ21

=Ω∈=

=Ω∈=

τττ

τττ

vη

uuu

H

H

V

U (3.12)

H1 é o espaço de Hilbert de ordem 1.

Agora, aplica-se esta parametrização na formulação da derivada topológica:

Sejam os domínios definidos anteriormente em (3.4), eΩ e ee δ+Ω , a equação

(1.2), pode ser re-escrita como:

( ) ( )( ) ( )

ee

efeefD

RR

eee

ee

Tδ

δδ

ψψ δ

δ−+

Ω−Ω= +

→→

00

lim)ˆ(x (3.13)

Simplificando a equação com:

( ) ( ) ( )e

efeefefe

R δδ

δ

−+=

→0lim' (3.14)

Chega em:

( )( ) ( )

eefD eee

eReT δ

ψψ δδ

Ω−Ω= +

→→ 00lim

'1lim)ˆ(x (3.15)

Usando a relação (3.10) na equação (3.15), obtém-se:

18

( )( ) ( )

( ) 00

000

)('1lim1

lim'1lim)ˆ(

=→

=

→→

Ω=

=Ω−Ω

=

ττ

δ

ττττδ

ψτ

τψψ

dd

ef

efD

Ren

nReT

v

vx

(3.16)

Considerando-se )(ef R uma função regularizadora escolhida de modo que

∞≤≤ )ˆ(0 xTD .

Finalmente a derivada topológica está definida no domínio parametrizado da

mesma maneira que na análise de sensibilidade à mudança de forma, assim como no

trabalho de SOKOLOWSKI e ZOCHOWSKI[1].

3.3 Análise de Sensibilidade à Mudança de Forma

Após a definição do modelo parametrizado, parte-se para os conceitos da

análise de sensibilidade à mudança de forma.

Determina-se a sensibilidade da função custo )( τψ Ω em relação à perturbação

caracterizada pelo campo de velocidade v governada pelo parâmetroτ :

τψψ

ψτ

τττ

τττ

)()(lim)( 0

00

=

→=

Ω−Ω=Ω

dd (3.17)

Esta derivada é descrita em LUEMBERGER[13].

De uma maneira generalizada, a função custo é definida como:

19

ττττττ

τ

τ

τ

τφφψ Γ+Ω=Ψ=Ω ∫∫ ΓΓ

ΩΩ dd )()()()( uuu (3.18)

τφΩ e

τφΓ são funções regulares definidas em τΩ e τΓ ,e τu é uma função

implícita deτ através da equação de estado(3.11).

Então a derivada (3.17) é definida por:

τψψ

ψτ

τττ

τ)()(

lim)( 000

uu −=Ω

→=dd (3.19)

Observa-se que euuu ===0

0τ

τ é a solução associada ao domínio eΩ do

problema variacional dado pela equação (3.2).

Formalmente, o cálculo da derivada topológica pode ser escrito da seguinte

forma:

Calcule:

0))((:que tal

)(0

=

Ψ

Ω

=

τ

ττ

τφ

τ

u

u

h

dd

(3.20)

3.4 Método de Lagrange

Para se determinar o valor da derivada topológica, deve-se escolher o melhor

método para resolver as equações (3.20)(ver HAUG et al.[14]). Os métodos existentes

para isto são: métodos direto, método adjunto e o método de Lagrange. A escolha do

método varia de problema para problema, e neste trabalho se escolherá o método do

20

lagrangeano devido a sua fácil generalização para variadas funções de restrição e custo.

Este método consiste em relaxar as restrições do problema através de multiplicadores,

chamados de multiplicadores de Lagrange. Assim o lagrangeano da equação na forma

perturbada, pode ser escrito como:

τττττττ λλλ V∈∀+Ψ= ),()(),( uuu hL (3.21)

Se a equação de estado (3.11) é satisfeita para todo τ , basta calcular a derivada

total do lagrangeano, respeitando as condições necessárias de primeira ordem de

Lagrange, em relação aτ para se obter a sensibilidade da função custo:

( ) ( ) ( )λλτ

λλτ

λτ ττ

τ

τττ

τττττ ∂

∂+

∂∂

∂∂

+∂∂

= ),(),(),(),( uu

uu

uu LLLLdd

(3.22)

Substituindo (3.21) em (3.22) e denominando xx&=

∂∂τ

, tem-se:

),(),(),(),( ττττττττ λλλτ

λτ

&& uuuu LLLL ++∂∂

=dd (3.23)

A equação adjunta é representada por:

0),( =∂∂

τττ

ηuu

L (3.24)

Dependendo do problema a ser resolvido, os procedimentos de cálculo serão

diferentes. Cada caso será estudado à parte.

21

3.5 Cálculo da Derivada Topológica

Considera-se o problema formulado em (3.11), de maneira que no domínio

perturbado a equação de equilíbrio bi linear e o funcional linear permanecem válidos:

+∈∀∈∀= Rττττττ ,)(),( Vηηηu la (3.25)

Com:

ττττ

τ

Ω∇= ∫Ω d.a sηTηu ),( (3.26)

ττττ

ττ

Γ+= ∫∫ ΓΩ

dd.l ηqηbη .Ω)( (3.27)

O cálculo da derivada topológica será realizado considerando como função

custo a energia potencial total definida por:

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Γ+−Ω∇=Ψ ∫∫∫ ΓΩΩττττττ

τττ

dd.d. uqubsuTu .Ω21)( (3.28)

Ou ainda:

)(),(21)( ττττ uuuu la −=Ψ (3.29)

Relembrando a definição da derivada topológica (3.20) e o seu lagrangeano

(3.21), basta calcular a derivada total do lagrangeano com relação ao parâmetro τ , para

se obter a derivada topológica.

O problema fica da seguinte maneira:

Calcule:

22

+

=

∈∀∈∀=

⎟⎠⎞

⎜⎝⎛ −

Rτ

τ

ττττ

ττττ

,)(),(

:que tal

)(),(21

0

Vtla

ladd

ηηηu

uuu

(3.30)

O lagrangeano deste problema é:

ττ

ττττττττ

V∈∀

−+−=

μ

μμuuuuμu )(),()(),(21),( lalaL

(3.31)

Da equação (3.23), obtém-se a derivada do lagrangeano, que neste caso, é a

própria função custo:

⎟⎠⎞

⎜⎝⎛ −

∂∂

= )(),(21),( τττττ ττ

uuuμu ladd L (3.32)

Portanto a equação acima, mostra que a derivada parcial é igual à derivada total

com τu e τμ fixos. Então basta calcular a derivada parcial da função custo para se obter a

derivada topológica. Avaliando esta derivada para 0=τ , obtém-se a derivada

topológica.

( )000

)(),(21

=== ∂∂

−⎟⎠⎞

⎜⎝⎛

∂∂

=Ψτ

ττ

τττ τττ

uuu ladd (3.33)

Para se obter o valor analítico desta derivada, as seguintes formulações serão

usadas[11][15][16]:

Teorema do transporte de Reynolds:

23

Γ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛+

∂∂

=Γ∂∂

Ω⎟⎟

⎠

⎞

⎜⎜

⎝

⎛+

∂∂

=Ω∂∂

∫∫∫∫Γ =Γ =

Ω =Ω =

dd

dd

v

v

div

div

00

00

φτφ

φτ

φτφ

φτ

τ

τ

τττ

τ

τ

τττ

τ

τ (3.34)

Teorema do transporte de Reynolds na configuração perturbada:

∫∫∫∫∫

ΩΓ

=Γ =

ΓΩ =Ω =

Γ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛+∇+

∂∂

=Γ∂∂

Γ+Ω∂∂

=Ω∂∂

dd

dndd

vdivv

v

φφτφ

φτ

φτφ

φτ

τ

τ

τττ

τ

τ

τττ

τ

τ

.

.

00

00 (3.35)

Utilizando o teorema do transporte de Reynolds para calcular as derivadas,

obtém-se:

ττ

ττ

τ

τττ

τ

τ

τ

τ

Ω⎥⎦

⎤⎢⎣

⎡∇+∇

∂∂

=

=⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω∇

∂∂

∫

∫

Ω=

=Ω

d.

d.

see vuTsuT

suT

div.

21

0

0 (3.36)

A derivada do gradiente simétrico de um campo vetorial é dada por:

( )See. vuTsuT ∇∇−=∇∂∂

=.2

0ττττ

(3.37)

Substituindo (3.37) em (3.36), tem-se:

24

( )[ ]

( )[ ] τ

τ

τ

τττ

τ

ττ

τ

Ω∇∇−∇=

=Ω∇∇−∇=⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω∇

∂∂

∫

∫∫

Ω

Ω=

Ω

d

dd.

Te

see

See

see

vTuIuT

vuTvuTsuT

. 2.

.2div.21

e

0

(3.38)

A derivada do funcional linear )( τul fica da mesma maneira, analisado no

ponto 0=τ :

( ) ee dl

e

Ω∇=∂∂ ∫Ω=

vIub ..0τ

ττ (3.39)

Substituindo as equações (3.38) e (3.39) na equação (3.36), chega-se na

relação:

ee dl

e

Ω∇Σ=∂∂ ∫Ω=

v.0τ

ττ

(3.40)

Com ( ) eTee

Seee TuIubuT ∇−−∇=Σ .2.

21 sendo tensor momento de energia

de Eshelby. Interpreta-se o tensor de Eshelby[17] como sendo um conjunto de forças

necessárias para provocar uma mudança de forma no domínio Ωe caracteriza por v∇ .

Considerando a relação tensorial uTuTuTdiv ∇+= ..)div()( T e aplicando o

Teorema da Divergência[15], a equação (3.40) se torna:

eeee ddl

ee

ΩΣ+ΓΣ=∂∂ ∫∫ ΩΓ

=vvn .div.

0τ

ττ

(3.41)

25

Proposição: Sejam uτ e τλ , respectivamente, soluções da equação de estado e

adjunta, então:

0div =Σe (3.42)

Prova: Considerando o teorema do transporte de Reynolds na configuração

perturbada e a relação:

( ) vuuu

.00

∇+′=∂∂

==

τττ

ττ

(3.43)

Para ( ) ( ) fixoττττ

´..0=

∂∂

=.

As derivadas da forma bi linear ),( ττ uua e do funcional )( τul , avaliadas em

0=τ , podem ser representadas como:

( ) ( )∫∫ ΓΩ=Γ∇+Ω′∇=

∂∂

ee

eseee

see dda nvuTuT ...2

0τττ

(3.44)

( )

( )

( )[ ]∫∫∫∫∫

∫∫

ΓΓΩ

ΓΓ

ΓΩ=

Γ∇++Γ′+Ω′=

=Γ∇+Γ′+

+Γ+Ω′=∂∂

eNe

NN

ee

eeTeeNeeee

NeNe

eeee

ddd

dd

ddl

vnTunubuqub

vuquq

nvubub

....

..

...0τ

ττ

(3.45)

Substituindo as relações (3.44) e (3.45) na equação (3.33). Chega em

26

( )

( )

∫∫∫∫

∫∫

ΓΓ

ΓΓ

ΩΩ=

ΓΣ=ΓΣ+′−′=

=ΓΣ+Γ∇

−Ω′−Ω′∇=∂∂

ee

eN

ee

eeee

eeNe

eeees

ee

ddla

dd

dbdl

vnvnuuu

vnvuq

uuT

..)(),(

..

..0τ

ττ

(3.46)

Comparando as expressões (3.40) e (3.46), e lembrando que para este caso

particular ee V∈′u , pode-se concluir que o tensor de Eshelby possui divergência nula:

vv ∀=Σ⇒=ΩΣ∫Ω ,0div0.div eee d

e

(3.47)

E a proposição está provada .

A prova desta proposição pode ser considerada uma forma alternativa para o

cálculo das derivadas com relação ao domínio parametrizado. Esta forma se mostra

mais simples e direta uma vez que as derivações são realizadas na configuração

perturbada τΩ .

Substituindo a definição de campo de velocidade (3.7) na equação (3.46), tem-

se:

∫∂=∂Σ=

∂∂

eBeen Bdnnv .

0τ

ττ

L (3.48)

Substituindo o resultado dos cálculos (3.48) na definição de derivada

topológica (3.16), e considerando uma expansão no furo, ou seja, 1)sign( −=nv :

27

( ) ∫∂→∂Σ

′−=

eBee

ReT Bd

efD nnx .1lim)ˆ(

0δ (3.49)

Substituindo a relação cinemática(2.2), desconsiderando as forças de corpo,

0=b , e usando as condições de contorno do problema de elasticidade linear definidas

em (2.4) e (2.5), tem-se 0=nTe . Então a relação acima se torna:

( ) ∫∂→∂⎟

⎠⎞

⎜⎝⎛

′−=

eBeee

ReT Bd

efD ETx .

211lim)ˆ(

0δ (3.50)

Para facilitar os cálculos, considera-se um sistema de coordenadas curvilíneas:

( ) ( ) ( )ttTntTtnTnT tte

tne

nte

nnee ⊗+⊗+⊗+=T (3.51)

n e t são vetores normal e tangencial, ortogonais entre si (n . t=0), e a condição

de contorno do problema (2.6), é escrita em coordenadas curvilíneas:

( ) ( ) etne

nte

nnee BntTtnTnTn ∂=⊗+⊗+= sobre 0.T (3.52)

Então:

ttee T=T (3.53)

Substituindo na equação (3.50):

( ) ( )∫∂∈→∂⎟⎟

⎠

⎞⎜⎜⎝

⎛′

−=

eBe

tte

RT BdT

tefD

2

0 211lim)ˆ(Dδ

x (3.54)

Para se obter a expressão final da derivada topológica, necessita-se escolher

adequadamente a função ( )ef R′ , com a restrição de atender o teorema (3.16). Esta

28

escolha é feita fazendo uma análise assintótica[18] da solução eu ou realizando

experimentos numéricos.

Seja um sistema de coordenadas cujos eixos estão alinhados com as tensões

principais do tensor tensão T associado ao domínio original, e com sua origem

posicionada no centro x da esfera B . Adotando um sistema polar de coordenadas

),( θR , a componente tteT definida sobre eB∂ é dada por:

θτ

2cos)(2)(0

RDRST tte −=

= (3.55)

Com 0→R , quando 0→e e eR >> . R é o raio da expansão.

Com:

)()()()()()(

21

21

RRRDRRRS

σσσσ−=+=

(3.56)

Substituindo a equação (3.55) em (3.54), conclui-se que a função

regularizadora )(efR deve ter um valor igual a:

2)()(2)()( eBáreaefeBperímetroef eReR ππ −=−=⇒−=∂−=′ (3.57)

Considerando este último resultado, calcula-se o limite da equação (3.54), que

resulta em:

( ) ( )[ ]221

221 2

21)ˆ( σσσσρ

−++=E

DT x (3.58)

Sendo 1σ e 2σ as tensões principais do tensor tensão xT ˆ avaliado no ponto

Ω∈x .

29

Reescrevendo a relação com um pouco de algebrismo, resulta na expressão

final da derivada topológica para o estado plano de tensões:

ETETx tr.tr)1(2

13.1

2)ˆ( 2νν

ν −

−+

+=TD (3.59)

Para o estado plano de deformações[11], tem-se:

ETETx tr.tr)1(2

)14)(1(.)1(2)ˆ( 2νννν

−

−−+−=TD (3.60)

Para o caso geral tridimensional, tem-se[19]:

⎥⎦⎤

⎢⎣⎡

−−

−−−

= ETETx tr.tr2151.10

)57(4)1(3)ˆ(

νν

νν

TD (3.61)

Para o caso especial em que 31

=ν e 41

=ν , para as equações (3.59) e (3.60),

respectivamente, tem-se a derivada topológica:

ETx .23)ˆ( =TD (3.62)

Para o caso tridimensional, quando 51

=ν , tem-se:

ETx .)ˆ( =TD (3.63)

Portando, obteve-se o valor analítico da derivada topológica apenas

conhecendo o campo de tensões T e deformações E, para qualquer ponto x pertencente

ao domínio.

30

4 Otimização Topológica de Sólidos

Neste capítulo, formula-se o problema da otimização topológica considerando

a energia de deformação interna como função objetivo e restrição na tensão máxima de

Mises permitida no domínio, levando em conta o equilíbrio do sólido. A derivada

topológica é aplicada ao algoritmo de otimização com o objetivo de se obter o ponto

ótimo do problema de otimização topológica proposto.

4.1 Formulação do Problema de Minimizar a Energia de

Deformação com Restrições nas Tensões

O modelo utilizado para a formulação do problema é o mesmo proposto por

BENDSOE[20], porém considera-se a formulação direta sem a aplicação de um modelo

de penalização.

Definição: Considere um elemento mecânico como um corpo ocupando um

domínio Ω, que faz parte do domínio total Ω em R2 ou R3. O domínio de referência é

escolhido de modo que seja possível aplicar restrições de carregamentos e condições de

contorno. Referindo-se ao domínio Ω, define-se o problema de otimização topológica

como a tarefa de achar o tensor ótimo de elasticidade D, que é variável ao longo do

31

domínio. Relembrando o capítulo anterior: a forma de energia bi linear (2.9),

deformações linearizadas (2.2) e (2.3), carregamento na forma linear (2.10). O problema

da mínima flexibilidade (máxima rigidez global) toma a forma:

⎩⎨⎧

≤∈∀=

−

eVM

,

,Tla

la

σ21))(max(),(),(

que tal

)(),(21min

uηηηu

uuuu

U

D

(4.1)

TVM é a tensão de Mises, e σe é a tensão de escoamento do material. max é

uma função que fornece o valor máximo.

As equações deste modelo são todas escritas na sua forma “fraca”, ou seja,

forma variacional com U denotando o espaço cinematicamente admissível de

deslocamentos. Este problema converge para um resultado ótimo, mas é extremamente

dependente do modo como foi feita a discretização do domínio[20].

4.2 Formulação da Distribuição de Material

Para se obter a formulação do problema, ou seja, trabalhar com o domínio do

problema de modo que se possa acrescentar ou remover material, é necessário

formalizar como será a definição de distribuição de material representada pelo tensor D,

ao longo do domínio do problema Ω. Em outras palavras, precisa-se de um modelo de

distribuição de material. Em BENDSOE[20], encontra-se a formulação do

comportamento do material.

Definição: O valor do tensor de elasticidade D será nulo onde não existir

material, com relação ao domínio inicial Ω0.

32

Esta afirmação é representada por:

⎩⎨⎧

ΩΩ∈Ω∈

== ΩΩ \01

1,)(10

0 xx

xsese

DD (4.2)

D é o tensor de elasticidade para um dado material isotrópico, pertencente ao

domínio Ω. 1Ω é uma função que indica quais pontos do domínio inicial, Ω0, pertencem

ao domínio ótimo, Ω.

4.3 Discretização do Domínio

Para se obter uma solução para o problema da otimização topológica proposto

será necessário utilizar uma técnica de discretização do domínio, que no caso será o

método dos elementos finitos. Com o domínio discretizado, as formulações propostas

anteriormente sofrem alterações. Porém, as equações diferenciais de Euler para o

problema da elasticidade (2.6) devem permanecer válidas.

O domínio discretizado pode ser interpretado como uma estrutura composta

Ilustração 4.1 – Discretização do domínio.

por diversas células menores, os elementos, e estas células podem ter um valor binário 0

ou 1. Ou também, para melhor ilustrar, pode-se pensar em uma foto digitalizada em 2

Ω

Ωelem

33

bits, preto e branco, na qual os pixels representam cada elemento e estes podem apenas

ter a cor preta ou branca. Ver Ilustração 4.1.

O modelo de distribuição de material e a formulação da derivada topológica

são modificados de modo que:

qelemTDV R∈ΩΩΩ x),(,,1

0 (4.3)

q é o número de elementos do domínio Ω.

Para o elemento elemΩ , tem-se:

⎩⎨⎧

ΩΩ∈ΩΩ∈Ω

=Ω= ΩΩ \01

1,)(10

0elem

elemelem se

seDD (4.4)

D é o tensor de elasticidade individual para cada elemento da malha.

4.4 Aplicação da Derivada Topológica no Problema de

Otimização

A derivada topológica fornece os pontos do domínio menos sensíveis à função

custo, que no caso estudado é a energia interna de deformação. O valor da derivada

topológica será calculado no centro de cada elemento Ωelem pertencente ao domínio Ω.

Eliminando os elementos menos sensíveis, ou seja, que tenham a derivada topológica

menor que um parâmetro supξ , é garantido que no domínio perturbado se tem a melhor

distribuição de material para determinado volume fixo VV ∂− . É importante observar

que o parâmetro supξ deve ser escolhido de modo criterioso[11], de maneira que haja

um decréscimo infinitesimal no volume. Devido à discretização do domínio pelo

método dos elementos finitos este valor infinitesimal V∂ é aproximado para um VΔ .

Através de testes experimentais realizados em [11] se chegou a conclusão de

34

que VV %1=Δ é um valor admissível. Portanto, no algoritmo proposto a seguir, a

redução no volume admitida a cada iteração será sempre de 1%. O algoritmo descrito a

seguir cria furos através da eliminação dos nós da malha de elementos finitos porque

este procedimento reduz instabilidades numéricas[11].

4.4.1 Eliminação de elementos no domínio

A derivada topológica é calculada no baricentro de cada elemento pertencente

ao domínio. Porém, para facilitar a implementação do algoritmo de otimização e evitar

instabilidades numéricas, será calculada a derivada topológica em função dos nós e não

em função dos elementos. A nova definição da derivada topológica em função do nó é:

))(max()(jelemTnóT DD Ω=Ω (4.5)

Em que jelemΩ são os elementos que contém o referido nó e max é uma função

que retorna o valor máximo. Ωnó é o domínio que representa a união de todos os

elementos que contém este nó. Ver Ilustração 4.2.

Ilustração 4.2 – Definição de Ωelem e Ωnó.

O critério de atualização do vetor de índices e de eliminação dos nós é definido

por:

⎪⎩

⎪⎨⎧

≥Ω

≤Ω=ΩΩ

jnóTjnóT

nój

j

Dse

Dse

sup

sup

)(,1

)(,0)(1 ξ

ξ (4.6)

Ωelem

Ωelem Ωelem

Ωelem Ωnó

35

Este critério é calculado para todos os nós pertencentes ao domínio. Ver

Ilustração 4.3.

Ilustração 4.3 – Representação da eliminação de elementos.

4.4.2 Cálculo do parâmetro de eliminação ξ

O parâmetro de eliminação é calculado arranjando o vetor de derivadas

topológicas em função dos nós em ordem crescente e em seguida seleciona-se o

elemento que corresponde a 1% do tamanho total do vetor. Este elemento é definido

como ξsup. Ver Ilustração 4.4.

Ilustração 4.4 – Definição do parâmetro de corte ξsup.

nó que atende ao critério

Ω Ωe

furo no domínio

sort(DT(Ωnó) =

1% ξsup

36

4.4.3 Algoritmo para a Criação dos Furos

Este algoritmo é conhecido como método “hard kill” devido à sua

característica de se eliminar os elementos da malha.

Algoritmo:

Seja a seqüência jΩ dos domínios de projeto, em que j é a j-ésima iteração:

1. Fornecer:

a. Domínio inicial Ωj=0 discretizado pelo método dos elementos

finitos;

b. As condições de contorno ( carregamentos externos e pontos de

apoio);

c. O critério de parada nas tensões: eVMT σ2,1≤ ;

2. Calcular a derivada topológica jTD , considerando o vetor de índices

Ω1 ;

3. Calcular jsupξ de modo que jj VV %1=Δ ;

4. Atualizar o vetor de índices Ω1 (Criar os furos):

⎪⎩

⎪⎨⎧

≥Ω

≤Ω=ΩΩ

jnóTjnóT

nój

j

Dse

Dse

sup

sup

)(,1

)(,0)(1 ξ

ξ;

5. Definir o novo domínio Ωj=j+1;

6. Voltar para o passo 1 enquanto não atingir o critério de parada;

7. Determinar geometria final de projeto e realizar o pós-processamento.

37

5 Alivio de Tensões com a Otimização Estrutural

Evolucionária

No ano de 1993, XIE e STEVEN[21] desenvolveram uma metodologia para

resolver problemas de otimização estrutural chamada otimização estrutural

evolucionária. Esta metodologia é baseada na idéia de que uma estrutura ótima pode ser

gerada através da eliminação gradual do material pouco atuante no domínio de projeto.

Há muitos atrativos para esta formulação, como por exemplo, a facilidade de

implementação do algoritmo e a velocidade de cálculo das iterações. Porém, para

algumas funções custo, não existe prova analítica de que as iterações convergem para

um mínimo local do problema de minimização.

Muitos trabalhos estão sendo publicados nesta área de pesquisa. Por exemplo,

em TANSKANEN [22]. Neste artigo ele enfatiza o aspecto teórico da otimização

estrutural evolutiva. Já em RONG et. Al [23], é apresentada uma aplicação em sólidos

tridimensionais. Para problemas em larga escala, JANG e KWAK[24] desenvolveram

uma metodologia de refinamento de malha a cada iteração evolutiva do algoritmo.

38

5.1 Formulação de Alívio de Tensões

A proposta deste capítulo é utilizar os conceitos de otimização evolutiva e

derivada topológica juntos. Após calcular uma geometria ótima utilizando o algoritmo

de otimização topológica, algumas vezes a tensão de Mises excede a tensão permitida

de projeto. A idéia é trabalhar com a fronteira do domínio e adicionar gradualmente

elementos onde a derivada topológica tiver maior valor. Ver Ilustração 5.1.


Este processo iterativo ameniza concentradores de tensão ao longo do contorno

do domínio devido ao aumento de rigidez na estrutura naquele ponto. Para garantir a

eficiência desta formulação o crescimento do domínio dever ser infinitesimal.

Consequentemente a variação do volume ΓΓ ∂+ VV também deve ser. Ao discretizar o

domínio pelo método dos elementos finitos, esta variação se torna ΓΓ =Δ VV %1 . Como

é mostrado em [11] , este valor de acréscimo de material é uma boa aproximação para

se obter uma distribuição uniforme e viável de material ao longo do domínio ativo.

eΩ xΩ

eΓ

eΓ

39

Similarmente ao algoritmo de otimização topológica, necessita-se calcular um valor

infξ no qual o material adicionado seja ΓΓ =Δ VV %1 .

5.1.1 Adição de elementos na fronteira

A derivada topológica é calculada no baricentro do elemento. Porém, para

facilitar a implementação computacional e evitar instabilidades numéricas, trabalhar-se-

á com a derivada topológica em função dos nós. Assim como no algoritmo de

otimização topológica, a seguinte expressão é definida:

))(max()(jelemTnóT DD Ω=Ω (5.1)

Em que jelemΩ são os elementos que contém o referido nó e max é uma função

que retorna o valor máximo. Ωnó é o domínio que representa a união de todos os

elementos que contém este nó.

O critério de atualização do vetor de índices e de adição dos elementos na

fronteira é definido por:

inf)(,1)(1 ξ≥Ω=ΩΩ nóTnó jDse (5.2)

Este critério é calculado para todos os nós pertencentes à fronteira, e se o

critério de adição for válido, os elementos pertencentes a este nó da fronteira são

adicionados ao domínio. Ver Ilustração 5.2.

Ilustração 5.2 – Figura representativa da adição de elementos na fronteira no nó que atende ao

critério de adição.

Г Гe

nó que atende ao critério

40

5.1.2 Cálculo do parâmetro de adição

O parâmetro de adição, assim como na otimização topológica, é calculado

arranjando o vetor de derivadas topológicas em função dos nós em ordem crescente e

em seguida seleciona-se o elemento que corresponde a 1% do tamanho total do vetor.

Este elemento é definido como ξinf. Ver Ilustração 5.3.


Algoritmo:

Seja a seqüência jΩ dos domínios de projeto, em que j é a j-ésima iteração,

então:

1. Fornecer:

a. Domínio inicial Ω0 discretizado;

b. As condições de contorno ( carregamentos externos e pontos de

apoio);

c. O critério de parada nas tensões: eVMT σ< ;

2. Calcular a Derivada Topológica DTj utilizando a análise estrutural

através do método dos elementos finitos considerando o vetor de

índices Ω1 ;

sort(DT(Ωnó) =

1% ξinf

41

3. Calcular infξ de modo que ΓΓ =Δ VV %1

4. Atualizar o vetor de índices jΓ

1 (criar as incrustações):

jnóTnó jDse inf)(,1)(1 ξ≥Ω=ΩΩ ;

5. Definir o novo domínio Ωj=j+1;

6. Voltar para o passo 1 enquanto não atingir o critério de parada;

7. Determinar geometria final de projeto e realizar o pós-processamento.

42

6 Comparação dos Métodos de Otimização

Topológica

Devido à otimização topológica ter um grande desenvolvimento nas últimas

décadas, surgiram dois grupos de classificação para os métodos de solução para a

otimização topológica: aproximação por microestrutura ou aproximação por

macroestrutura[25].

6.1 Aproximação por Microestrutura

A aproximação por microestrutura usa a estratégia de variar a densidade dos

pontos do domínio admitindo existir um material compósito capaz de ter a sua rigidez

variando entre zero (sem material) até um (com material). O método da homogenização

e o método SIMP[26](Solid Isotropic Material with Penalization) são métodos eficazes

para solucionar o problema de minimizar a energia interna do problema com restrição

de volume. Já o método da relaxação por campo de fase provou ser uma ferramenta útil

para se resolver problemas nos quais existem muitas restrições de desigualdade , como

por exemplo, tensões locais[2]. Em BURGER e STAINKO[27] é demonstrado um

método eficaz para se resolver o problema de minimizar o volume com restrições locais

43

de tensão. Este método é baseado na técnica da relaxação-ε como proposto por

CHENGO e GOU[28], mas adaptada para muitas restrições.

6.2 Aproximação por Macroestrutura

Nesta aproximação, a topologia do domínio é modificada através da inserção

de furos em lugares onde a sensibilidade com relação à função custo tiver valor menor.

Podem-se classificar os métodos desta categoria em métodos baseados em sistemas

contínuos ou sistemas discretos. Um exemplo da aplicação de aproximação por

macroestrutura em sistemas discretos é encontrada em XIE e STEVEN[29], a chamada

ESO, Evolutionary Structural Optimization, ou Otimização Estrutural Evolutiva.

Para resolver o problema da Otimização Topológica, usar-se-á neste trabalho

uma aproximação por Macroestrutura baseada em Sistemas Contínuos, conhecida como

Derivada Topológica, ou Análise de Sensibilidade Topológica.

No trabalho de NOVOTNY[11], foi proposta uma nova relação entre a

Sensibilidade à Mudança de Forma, comumente usada na Otimização de Forma[1] e a

Derivada Topológica, resultando em um teorema que conduz a uma vasta classe de

problemas de Engenharia. O ferramental matemático utilizado para as deduções,

inicialmente foi desenvolvido para a mecânica do contínuo e se adaptou perfeitamente

para a Análise de Sensibilidade à mudança de forma. Esta nova formulação não

apresenta restrições com relação à função custo, função de restrições, condições de

contorno e ao fenômeno estudado. Restrições estas, que eram encontradas em trabalhos

anteriores, como por exemplo, em [30], que ainda utilizavam uma formulação mais

simples do problema.

44

6.3 Comparação entre Métodos

A característica comum a estes dois grupos de classificação, microestrutura e

macroestrutura é a utilização da sensibilidade da função objetivo para a obtenção do

resultado ótimo do problema. No caso particular de se usar a energia interna como

função objetivo, estas formulações geram sensibilidades semelhantes, diferindo apenas

de constantes. Esta sensibilidade é definida por e

eρψ∂∂ para a formulação SIMP e ESO e

definida como )ˆ(xTD na derivada topológica. Portanto todas as expressões acarretam

soluções semelhantes. Como mostrado em LABANOWSKI et al.[31]. Na formulação

SIMP, a sensibilidade e

eρψ∂∂ é calculada levando em conta a variação da densidade de

um ponto do domínio:

ess

e

e de

Ω∇∇=∂∂ ∫

Ω

uu .21 *D

ρψ

(6.1)

com DD 1* −= pepρ , correspondendo ao tensor de elasticidade do material com

densidade intermediária, p é o parâmetro de penalização e ρ é a densidade de um

elemento que pode variar entre 0 e 1.

Na formulação da otimização estrutural evolutiva (ESO)[32] a sensibilidade

e

eρψ∂∂ é calculada ao se retirar um elemento finito do espaço de aproximação do

domínio, e resulta em:

ess

e

e de

Ω∇∇=∂∂ ∫

Ω

uu .21

Dρψ

(6.2)

A expressão da derivada topológica, para ν = 1/3 , é representada por:

45

sseTD uux ∇∇=

∂∂

==

.23)ˆ(

0D

ττψ

(6.3)

Pode-se interpretar o resultado (6.3) da derivada topológica como sendo a

formulação local da equação ESO, equação (6.2), que por sua vez é igual a expressão da

formulação SIMP, equação (6.1), no ponto ótimo do problema. Ou seja, quando

DD =* .

6.4 Tensão de Mises x Energia de Deformação Interna

Com o objetivo de se utilizar a sensibilidade à mudança de topologia no

algoritmo de alívio de tensões, utilizou-se um trabalho de STEVEN et al.[33] [34], no

qual é comparado o critério falha de Von Mises avaliado em um ponto do domínio com

a energia interna de deformação avaliada neste mesmo ponto, para o caso do estado

plano de tensões.

Levando em conta as relações da teoria da elasticidade linear (2.1), (2.2) e

(2.3), o vetor de tensões, no estado plano de tensão, é representado por:

suT ∇= D (6.4)

Para o estado plano de tensões, a tensão de Mises é definida como:

TVTT TVM =2 (6.5)

Com:

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡−

−=

300015,005,01

V (6.6)

Substituindo (6.4) em (6.5), obtém-se:

46

( ) sTTsVM uVuT ∇∇= DD2 (6.7)

Chamando ZDD =VT , tem-se:

( ) sTsVM uuT ∇∇= Z2 (6.8)

Analisando a energia de deformação interna, tem-se:

( ) sTsTEn uuET ∇∇== D21.

21 (6.9)

T e E estão expressos em forma vetorial considerando o estado plano de

tensões.

Comparando as expressões (6.8) e (6.9), conclui-se que os valores calculados

diferem apenas dos valores das constantes Z e D . Para o estado plano de tensões, por

exemplo, que é o caso desta monografia, tem-se:

( ) ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

=2/100

0101

1 2

νν

ν

νE

D (6.10)

Calculando o valor de Z :

( ) ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−+−−+−−+−+−

⎟⎠

⎞⎜⎝

⎛−

=2

22

222

2175,000015,025,005,025,01

1 ννννννννν

νE

Z (6.11)

Portanto, as duas formulações, da tensão de Von Mises e da energia de

deformação, são parecidas, diferindo apenas dos valores das constantes.

47

7 Método dos Elementos finitos

7.1 Visão Geral

O método dos elementos finitos é uma ferramenta essencial para a engenharia.

A sua simplicidade e robustez o torna muito atraente para solucionar diversos

problemas matemáticos que até então não tinham sequer uma solução aproximada.

Além disso, este método é generalizado, ou seja, sempre o mesmo procedimento é

adotado para diferentes problemas, facilitando assim a confecção de algoritmos e a

implementação em linguagens de programação.

Muitas das vezes, uma solução analítica para um problema se torna muito

trabalhosa ou até mesmo impossível de se obter devido a uma série de fatores:

geometrias complexas, estados de carga complexos, restrições complexas. O conceito

do método dos elementos finitos tenta contornar esta dificuldade assumindo o domínio

contínuo do problema ser um domínio discreto.

Essa discretização permite obter soluções aproximadas ou até mesmo iguais à

solução analítica do problema. Em [35], encontra-se uma excelente fonte de pesquisa

para o método dos elementos finitos.

48

7.2 Discretização das Equações

O domínio Ω, que ora era contínuo, agora passa a ser discreto e contendo um

número finito de elementos. Apesar disso, cada elemento é contínuo individualmente. O

comportamento individual destes elementos com relação ao fenômeno estudado é de

fácil entendimento e será mostrado agora.

Sob o domínio Ω, será criado um número finito de elementos. Cada elemento

tem o seu próprio domínio contínuo Ωelem. Estes elementos são interconectados através

dos nós, que são pontos de interseção entre elementos. Ao conjunto de todos os

elementos e nós, chama-se malha.

Ilustração 7.1 – Esquema do problema discretizado.

Dado o problema de valor de contorno (2.7) cuja solução seja u. A

discretização do problema se baseia na aproximação desta solução u , baseada em

pontos nodais, para uma solução un, usando funções de interpolação iφ com certas

restrições que permitem determinar as constantes nααα ,...,, 21 , com n finito:

∑=

=

n

i

iin α

1

φu (7.1)

De maneira que:

Γ Ω im

j

n

Ωelem

49

uu →∞→

n

n (7.2)

Esta interpolação se dá nos deslocamentos dos nós para se obter um campo de

deslocamento diferenciável pertencente ao elemento. Para cada nó de um elemento,

resolve-se o seguinte sistema e determinam-se as constantes nααα ,...,, 21 :

⎩⎨⎧

==

0)(1)(

jii

iii

αα

XX

φφ

(7.3)

Sendo Xn um vetor que representa a coordenada nodal n do elemento.

Após o cálculo dos coeficientes, determina-se a função de forma do nó i, Ni:

iii φα=N (7.4)

A matriz de forma é obtida:

[ ]nNININ K1= (7.5)

I é a matriz identidade e n é o número de nós de um elemento.

Re-escrevendo a equação (7.1) , tem-se:

uNu =n (7.6)

A função de forma varia de elemento para elemento e geralmente é um

polinômio. Casos especiais desta função de forma serão descritos nos itens seguintes.

A equação (2.2) discretizada, torna-se:

uNE s∇= (7.7)

Simplificando a operação, define-se a matriz de derivada de forma B, que

relaciona diretamente as deformações com os deslocamentos nodais para cada elemento.

50

nós de número],,[

0

0

1

==

⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

∂∂

∂∂

∂∂

∂∂

=∇=

N

xN

yN

yN

xN

N

ii

i

i

ii

BBB

NB

K

(7.8)

Nesta definição, apresenta-se o caso bidimensional para a matriz B.

E a relação fica:

uE B= (7.9)

Usando a equação (2.7), e discretizando-a para um número limitado de

elementos no domínio, obtém-se:

Fu

FuuU

Fuuu

∑ ∫

∑ ∫

∑ ∫

=⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω

=⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω

=⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω

Ω

Ω

= Ω

N

iiTi

N

Tii

Ti

T

T

N

i

iiT

i

i

i

i

d

d

d

1

1

1

)(

)(

)(

DBB

DBB

BDB

(7.10)

F é um vetor de forças externas, u é vetor de deslocamentos nodais.

Define-se a rigidez de um elemento por:

∫Ω Ω=

i

iiTii d)( DBBk (7.11)

Observa-se que sempre a matriz de rigidez será simétrica, pelo fato da matriz

constitutiva D ser simétrica e haver uma dupla multiplicação da matriz de derivada de

forma B, como será mostrado mais adiante.

A matriz de rigidez global é definida por:

51

∑=N

i

1

kK (7.12)

Tem-se a expressão final da formulação:

KuF = (7.13)

Definida a teoria geral dos elementos finitos, o próximo passo é definir a

matriz constitutiva D e de a derivada de forma B, que variam de acordo com o tipo de

elemento em questão. Neste trabalho optou-se por usar três tipos de elementos: treliça

tridimensional, triângulo linear, e quadrado bi-linear. Todos estes elementos são de fácil

formulação e não necessitam da integração numérica para calcular a matriz de derivada

de forma.

7.3 Treliça Tridimensional

A treliça é um elemento estrutural com apenas uma dimensão e é formado pela

ligação de dois nós. É esbelto e só há esforço axial ao longo de sua estrutura. As

ligações entre um elemento e outro são em forma de rótula, portanto não transmitem

momento. Somente pode ser aplicada carga nos nós. Ilustração 7.2

Ilustração 7.2 – Representação do elemento treliça.

j

x z

i

dy

Ωel

, iii yx=d

52

A matriz de rigidez local do elemento é definida como:

⎥⎦

⎤⎢⎣

⎡−

−=

1111

LEA

D (7.14)

E é o modulo de elasticidade, A é a área da seção transversal da treliça e L é o

comprimento da barra. Observa-se que neste caso a matriz D é apenas uma matriz de

segunda ordem. Devido à simplicidade deste elemento e por conter apenas tensões no

sentido axial. O vetor de tensões neste caso é:

[ ]xσ=T (7.15)

A função de interpolação generalizada para a treliça é:

xNi 21 αα += (7.16)

O sistema a ser resolvido para se determinar os coeficientes α da função de

interpolação Ni, fica:

⎩⎨⎧

=+==+=

01

21

21

jj

ii

xNxN

αααα

(7.17)

E assim sucessivamente para Nj.

Calculam-se todos os coeficientes e se determina a função de forma para o

elemento. Usando as relações (7.5) e (7.7) , chega-se na matriz B:

⎥⎦

⎤⎢⎣

⎡=

cbacba

L 0000001

B (7.18)

( ) ( ) ( )212

212

212

121212 ;; com

zzyyxxL

zzayybxxa

−+−+−=

−=−=−=

A matriz de rigidez, definida pela equação (7.11), neste caso, se torna um

produto matricial, lembrando que as matrizes D e B são constantes:

iTi DBB=ik (7.19)

Calculando este produto, resulta em uma matriz simétrica:

53

⎥⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

−−−−−−

−

=

2

2

2

22

22

22

3

ccbbcabaa

ccbcaccbbbacbb

caba-ac-abaa

LEAk s

i (7.20)

121212 ;;a com zzcyybxx −=−=−=

7.4 Triângulo Linear

O triângulo linear é o elemento bidimensional mais simples de todos devido a

sua facilidade de cálculo e sua versatilidade. É constituído por três nós e é realizada uma

interpolação linear para se obter o campo de deslocamentos pertencentes ao elemento. O

estado plano de tensões será adotado para este elemento. Ilustração 7.3:

Ilustração 7.3 – Representação do elemento triângulo linear.

A matriz constitutiva do material para o estado plano de tensões é definida

como:

j

x

, iii yx=d

y i m

n

di

Ωele

54

( ) ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

=2/100

0101

1 2

νν

ν

νE

D (7.21)

E é o modulo de elasticidade, ν é o coeficiente de Poisson. Devido ao estado

plano de tensões, a matriz D nesse caso é de terceira ordem. Só existem três

componentes de tensão. São elas σx, σy, τxy, e formam o tensor de tensões de segunda

ordem:

⎥⎦

⎤⎢⎣

⎡=

yxy

xyx

σττσ

T (7.22)

A função de interpolação generalizada para o triângulo é:

yxNi 321 ααα ++= (7.23)



⎪⎩

⎪⎨

⎧

=+=

=+==+=

00

1

21

21

21

mm

jj

ii

xNxNxN

αα

αααα

(7.24)

E assim sucessivamente para Nj e Nm.



⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

Δ=

mmjji

mji

mji

bcbcbcccc

bbb000

000

21

B (7.25)

mm

jj

ii

jmi

mj

yxyxyx

ijm

yxc

yy

111

det) triângulodo (área 22

bcom

i

==Δ

−=

−=

55

Os outros coeficientes são obtidos por uma permutação cíclica dos subscritos

na ordem i, j, m.

Observa-se que a matriz de derivada de forma B é constante para qualquer

ponto dentro do elemento. Então a equação (7.11) se torna um produto matricial não

necessitando o uso da integração numérica.

iTit DBB=ik (7.26)

t é a espessura da malha que é considerada constante.

Calculando-se este produto, obtém-se a matriz de rigidez do elemento:

( )

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( ) ( )

( ) ( )

( ) ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

−+

−+−+

−+−+−+

−+−+−+−+

−+−+−+−+−+

−+−+−+−+−+−+

Δ−=

ν

ννν

νννν

νννννν

ννννννν

ννννννννν

ν

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

12

14k

22

22

22

22

22

222

22i

mm

mmmm

mm

mjmj

mjmj

jj

mjmj

mjmj

jjjj

jj

mimi

mimi

jiji

jiji

ii

mimi

mimi

jiji

jiji

iii

ii

bc

cbcbcbSim

bbcc

cbbc

bc

bccb

ccbb

bccb

cb

bbcccbbcbb

cccb

bcbc

bccbccbbbc

cbcc

bbccbcb

E

(7.27)

7.5 Quadrado Bi-linear de Malha Regular

Este elemento é composto de quatro nós e a função de interpolação usada para

o seu campo de deslocamento é linear. Devido à malha ser regular, este elemento

permite simplificações em sua formulação que facilitam o cálculo de sua matriz de

rigidez. Ou seja: 31 xx = e 31 yy = . Ilustração 7.4:

56

Ilustração 7.4 – Representação do elemento quadrado bi linear.

A matriz constitutiva do material é dada por:

( ) ⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

=2/100

0101

1 2

νν

ν

νE

D (7.28)

E é o modulo de elasticidade, ν é o módulo de Poisson.

O tensor de tensões é:

⎥⎦

⎤⎢⎣

⎡=

yxy

xyx

σττσ

T (7.29)

Observa-se que a matriz constitutiva e o tensor de tensões são o mesmo do

elemento triângulo linear devido ao estado plano de tensões ser adotado nos dois casos.

A função de interpolação generalizada para o quadrado é:

yxNi 321 ααα ++= (7.30)



⎪⎪

⎩

⎪⎪

⎨

⎧

=+=

=+=

=+==+=

00

01

21

21

21

21

mm

jk

jj

ii

xNxN

xNxN

αα

αα

αααα

(7.31)

E assim sucessivamente para Nj , Nk e Nm.

y

x

i m

j

n

di

, iii yx=d

k

Ωelem

57



⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−−−−

−−

Δ=

cbcdadabbddb

ccaa0000

00001

B (7.32)

)).(()abcd quadrado do (área ;;;

com

1133

bdacxxdyycxxbyya−−==Δ

−=−=−=−=

Neste caso a matriz de derivada de forma B não é constante. Necessita-se

realizar a integração da equação (7.11) para se calcular a matriz de rigidez do elemento.

Não é necessária a utilização da integração numérica, pois a malha é regular e isto

facilita a obtenção da matriz de rigidez de uma maneira direta.

dxdytx

x

y

yi

Tii ∫ ∫=

3

1

3

1

)( DBBk (7.33)

t é a espessura da malha que é considerada constante.

Devido a complexidade do cálculo da equação (7.33), não será apresentada

aqui a expressão literal final da matriz de rigidez. No entanto, o calculo é obtido

facilmente utilizando uma ferramenta de cálculo simbólico como, por exemplo,

Mathematica e Maple. O resultado obtido do cálculo é então exportado para um

programa de elementos finitos para a utilização prática.

58

8 Métodos de Resolução de Sistemas Lineares

O método dos elementos finitos sempre resulta na solução de um sistema linear

com muitas equações. Relembrando a equação (7.13), necessita-se desenvolver técnicas

para se resolver de uma maneira satisfatória esse sistema para que o cálculo seja rápido

e eficiente. Para se determinar qual o melhor método para se utilizar, relacionam-se as

características desta matriz K:

1. Esparsa

2. Simétrica

3. Definida positiva

4. De banda

8.1 Esparsidade

Uma matriz é esparsa quando o número de elementos nulos é muito maior que

o número de elementos não nulos. Ilustração 8.1.

59

0 500 1000 1500

0

200

400

600

800

1000

1200

1400

1600

nz = 21650

Ilustração 8.1 – Exemplo de uma matriz esparsa de ordem 1722. Os elementos não nulos são os

pontos pretos neste diagrama totalizando nz = 21650.

A organização dos elementos desta matriz é otimizado da seguinte maneira: o

elemento jxin não nulo é armazenado na matriz esparsa junto com as coordenadas de

sua localização na matriz não esparsa. Deste modo, os elementos nulos não necessitam

serem armazenados de modo que esparsojxi Kn ∉= 0 . O vetor I e o vetor J armazenam

os índices de localização do elemento; na terceira coluna N, é armazenado o valor.

[ ]NJIK =esparso (8.1)

8.2 Matriz Simétrica

Uma matriz é simétrica quando ela é igual a sua transposta:

TAA = (8.2)

A matriz de rigidez para o caso da elasticidade linear, também é simétrica. A

equação (7.11) prova isso: observa-se que sempre a matriz de rigidez será simétrica,

pelo fato da matriz constitutiva D ser simétrica e haver uma dupla multiplicação da

matriz de derivada de forma B.

60

iTi

iiTii

TTi

TTii

Ti

Ti

Ti

T

iiTi

Ti

iii

i

ddd

d

kk

k

k

=

Ω=Ω=Ω=

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Ω=

∫∫∫∫

ΩΩΩ

Ω

)()()(

)(

DBBBDBDBB

DBB

(8.3)

8.3 Matriz Definida Positiva

Uma matriz é definida positiva quando:

0>vMvT (8.4)

Para qualquer nXnn RR ∈∈≠ Mv ,0

Re-escrevendo esta relação para a matriz de rigidez, tem-se:

0>= UKUTdefE (8.5)

Esta expressão é válida considerando que o resultado deste produto matricial é

a energia interna de deformação do sólido que sempre será positiva.

8.4 Matriz em Banda

Matriz em banda é aquela em que todos os elementos estão concentrados

próximo à diagonal principal da matriz.

Definição: A banda da matriz é a maior distância que um elemento não nulo

tem com a diagonal principal. Entende-se distância como a quantidade de diagonais

61

existentes com relação a diagonal principal. Uma matriz em banda sempre será esparsa,

mas nem sempre o oposto é verdadeiro. Ilustração 8.2

0 500 1000 1500

0

200

400

600

800

1000

1200

1400

1600

nz = 21650

Ilustração 8.2 – Representação gráfica dos elementos de uma matriz,destacando-se a banda

superior e inferior.

Uma vantagem de a matriz ser em banda, é reduzir a memória de

armazenamento para a matriz. Este método somente guarda os elementos que pertencem

à banda da matriz. Se a matriz for simétrica, ao mesmo tempo em que for de banda,

reduz-se mais ainda o tamanho da matriz armazenada, guardando somente uma das

metades simétricas da banda em questão.

O esquema de armazenamento para uma matriz em banda e simétrica é descrito

da seguinte forma:

),(),()1(

jjildbajianXldbnXn

−+⇒

⇒ +KK (8.6)

Com

simétrica. banda de largura a é ; , e ],1[];),1,1[max( ldbjijjjldbji N∈=+−=

62

8.5 Decomposição de Cholesky

Considerando todas as características da matriz de rigidez do problema -

simétrica, definida positiva e em banda - o melhor método para solucionar este sistema

de equações lineares é o método da decomposição de Cholesky, com o caso especial de

a matriz ser em banda. Esta decomposição resulta numa simplificação do tipo:

TLLK = (8.7)

L é uma matriz diagonal superior.

Esta simplificação facilita muitas operações que envolvem o cálculo da inversa

da matriz, ou quando a mesma matriz deve ser usada para diferentes vetores.

( ) ( )( ) FLLFKU

LLLLK111

1111

−−−

−−−−

==

==T

TT

(8.8)

Um algoritmo especial para a solução deste sistema é encontrado no LAPACK,

Linear Algebra Package[36]. Uma coletânea de algoritmos em FORTRAN 77 para

solucionar diversos tipos de problemas da álgebra linear. A grande vantagem desta

biblioteca virtual é a de que todas as funções que a integram suportam a plataforma de

multiprocessamento. Esta plataforma é possível através da linguagem OPEN MP[37],

que é uma série de diretivas de compilação e de variáveis que determinam o paralelismo

de processamento.

63

9 Exemplos

Todas as ferramentas necessárias para aplicar a derivada topológica e a

otimização evolucionária foram descritas nos capítulos anteriores. Com o objetivo de

validar toda a teoria proposta nesta dissertação, formularam-se alguns exemplos

existentes na literatura específica[38]. Os algoritmos citados nesta dissertação, de

otimização topológica e de alívio de tensões serão utilizados. Após executar o algoritmo

de otimização topológica, será executado o algoritmo de alívio de tensões partindo do

ponto ótimo obtido pelo algoritmo de topologia.

Em todos os exemplos será usado o mesmo material: aço carbono com módulo

de elasticidade GPaE 210= , tensão de escoamento MPae 250=σ , coeficiente de

Poisson 3,0=ν . O critério de parada para o algoritmo de otimização topológica é:

eVMT σ2,1≤ . O critério de parada para o algoritmo de alívio de tensões é:

eVMT σ< .

O elemento finito utilizado em todos os exemplos é o quadrado bi linear

isoparamétrico de malha regular.

64

9.1 Viga Bi-apoiada com Carga Distribuída na Parte Superior

Neste exemplo a espessura do material mt 01,0= é considerada constante.

Devido à simetria do problema, será discretizada apenas uma metade do domínio. Foi

utilizada uma malha de 3750 elementos. Ilustração 9.1:

Ilustração 9.1 – a) Geometria da viga (cotas em metro) b) Distribuição da Tensão de Von

Mises no domínio inicial (MPa).

10

20

30

40

50

60

70

80

90

100

Ω

0,6 0,2

1,2 0,6

1,0

3,0

ΩmkN /400

65

Definidas a geometria e condições de contorno do problema, executa-se o

algoritmo de otimização, cujo resultado final é mostrado na Ilustração 9.2, e o histórico

de iterações é mostrado na Ilustração 9.3:

Ilustração 9.2 – a) Resultado do Algoritmo de Otimização b) Distribuição da Tensão de Von


Ilustração 9.3 – Histórico de iterações do Volume, Energia e Tensão Máxima de Mises.

Ω

40

60

80

100

120

140

160

180

200

220

240

Ω

0 50 100 150 200 2500

0.005

0.01

0.015

0.02

0.025

0.03Iterações do Algoritmo - Volume

Vol

ume(

m3)

Iterações0 50 100 150 200 250

20

30

40

50

60

70

80

90

100Iterações do Algoritmo - Energia Interna

Ene

rgia

(J)

Iterações0 50 100 150 200 250

100

120

140

160

180

200

220

240

260

280

300Tensão de Mises Máxima

TMax

(MP

a)

Iterações

66

Partindo do ponto de parada do algoritmo de otimização, executa-se o alívio de

tensões, cujo resultado final é mostrado a seguir na Ilustração 9.4, e o histórico de

iterações é mostrado na Ilustração 9.5:

Ilustração 9.4 – a) Resultado do Alívio de Tensões b) Distribuição da Tensão de Von

Mises(MPa).


Ω

20

40

60

80

100

120

140

160

180

200

Ω

0 5 10 15 20 254.1

4.15

4.2

4.25

4.3

4.35

4.4

4.45x 10-3 Iterações do Algoritmo - Volume

Vol

ume(

m3)

Iterações0 5 10 15 20 25

84

85

86

87

88

89

90

91

92


Ene

rgia

(J)

Iterações0 5 10 15 20 25

245

250

255

260

265

270

275

280

285


TMax

(MP

a)

Iterações

67

Resumo do exemplo 1 na Tabela 9.1:

Viga bi apoiada com carga distribuída superior

Algoritmo de Otimização

Iteração

Descrição 0 175 208

Volume (m3) 0,0300 0,0056 0,0040

Energia Interna de Deformação (J) 26,003 65,28 92,470

Tensão Máxima de Von Mises (MPa) 118,9 248,08 288,46

Algoritmo de Alívio de Tensões

Iteração

Descrição 22

Volume final (m3) 0,0044

Energia Interna de Deformação final (J) 84,258

Tensão Máxima de Von Mises final (MPa) 248,84

Tabela 9.1

68

9.2 Viga Engastada Superiormente com Carga Distribuída

Cisalhante


Foi utilizada uma malha de 7500 elementos. Ilustração 9.6:

Ilustração 9.6 – a) Viga engastada com carga cisalhante (cotas em metro). b) Distribuição da

Tensão de Von Mises (MPa).

1,2 0,6

1,0

3,0

Ω

mkN /400

10

20

30

40

50

60

70

80

90

Ω

69


algoritmo de otimização, cujo resultado final é mostrado na Ilustração 9.7, e o histórico

de iterações é mostrado na Ilustração 9.8:

Ilustração 9.7 – Resultado após o algoritmo de Otimização. a) Distribuição de tensão (MPa).


Ω

40

60

80

100

120

140

160

180

200

220

240

Ω

0 50 100 150 200 2500

0.005

0.01

0.015

0.02

0.025


Vol

ume(

m3)

Iterações0 50 100 150 200 250

0

50

100

150

200


Ene

rgia

(J)

Iterações0 50 100 150 200 250

50

100

150

200

250


TMax

(MP

a)

Iterações

70




Ilustração 9.9 – Resultado após o algoritmo de Alívio de Tensão. a) Distribuição de tensão (MPa).

Ilustração 9.10 – Histórico do alívio de Tensões.

Ω

40

60

80

100

120

140

160

180

200

Ω

1 2 3 44.1

4.11

4.12

4.13

4.14

4.15

4.16

4.17x 10-3 Iterações do Algoritmo - Volume

Vol

ume(

m3)

Iterações1 2 3 4

234

235

236

237

238

239

240


Ene

rgia

(J)

Iterações1 2 3 4

240

245

250

255

260

265

270

275


TMax

(MP

a)

Iterações

71


Viga engastada com carga cisalhante distribuída


Iteração


Volume (m3) 0,0300 0,0048 0,0040




Iteração

Descrição 4




Tabela 9.2

72

9.3 Viga Bi-apoiada com Carga Distribuída Inferiormente


Devido à simetria do problema, será discretizada apenas uma metade do domínio. Foi

utilizada uma malha de 3750 elementos. Ilustração 9.11:

Ilustração 9.11 – a) Geometria da viga (cotas em metro) b) Distribuição da Tensão de Von


0,4 0,9

1,0

3,0

Ω

mkN /400

0,4

10

20

30

40

50

60

70

80

90

100

Ω

73


algoritmo de otimização, cujo resultado final é mostrado na Ilustração 9.12, e o

histórico de iterações é mostrado na Ilustração 9.13:

Ilustração 9.12 - Resultado após o algoritmo de Otimização. a) Distribuição de tensão (MPa).


Ω

40

60

80

100

120

140

160

180

200

220

240

Ω

0 50 100 150 2000

0.005

0.01

0.015

0.02

0.025


Vol

ume(

m3)

Iterações0 50 100 150 200

20

40

60

80

100

120

140


Ene

rgia

(J)

Iterações0 50 100 150 200

100

120

140

160

180

200

220

240

260


TMax

(MP

a)

Iterações

74




Ilustração 9.14 – Resultado após o algoritmo de Otimização. a) Distribuição da tensão de Mises (MPa).

Ilustração 9.15 – Histórico do alívio de Tensões.

Ω

20

40

60

80

100

120

140

160

180

200

Ω

0 5 10 15 204

4.2

4.4

4.6

4.8

5

5.2x 10-3Iterações do Algoritmo - Volume

Vol

ume(

m3)

Iterações0 5 10 15 20

125

130

135

140

145

150


Ene

rgia

(J)

Iterações0 5 10 15 20

240

250

260

270

280

290


TMax

(MP

a)

Iterações

75


Viga bi apoiada com carga distribuída inferior


Iteração


Volume (m3) 0,0300 0,0055 0,0040




Iteração

Descrição 19




Tabela 9.3

76

10 Implementação Computacional

Neste capitulo será mostrado como está organizada a estrutura das funções

computacionais que compõem este trabalho. A plataforma utilizada para a programação

foi o MatLab 6.5. A facilidade de programação com a linguagem MatLab facilitou

muito a implementação e a alteração do código fonte já existente de elementos finitos.

Este código necessitou apenas de algumas alterações para trabalhar em conjunto com a

formulação do algoritmo de otimização. Não houve nenhuma utilização de programas

comerciais de elementos finitos e nem de pós e pré-processamento. Com esta atitude,

conseguiu-se uma maior versatilidade nas funções programadas e também facilitou

alterações e modificações no programa de otimização.

10.1 Programação Otimizada

Procurou-se neste trabalho desenvolver um código fonte que seja eficiente em

velocidade e economia de memória. Para atingir esse objetivo alguns preceitos foram

considerados durante a fase de programação.

1. Trabalhar sempre com a formulação da matriz de rigidez esparsa para

economizar memória;

2. Evitar criar loops desnecessários, principalmente um loop dentro do

outro;

77

3. Trabalhar com os tipos de dados que são acelerados pelo MatLab:

logical, char, int8, int16, int32, uint8, uint16, uint32, double;

4. Evitar chamar funções durante um loop, pois diminuem a performance

do loop;

5. Evitar usar mais de uma operação por linha;

6. Evitar mudar o tipo de uma variável ou o tamanho de uma matriz ou

vetor;

7. Compilar funções gargalo em um arquivo tipo MEX, para melhorar a

velocidade do programa;

8. Trabalhar preferencialmente com funções e não com scripts;

9. Enquanto o programa estiver sendo executado, evitar carregar outros

aplicativos que comprometam a memória;

10. Pré-alocar valores iniciais de vetores e matrizes e definir seu tamanho

final;

11. Trabalhar preferencialmente com loops vetorizados. Isto significa

transformar o loop for e while em uma operação equivalente

vetorizada.

10.2 Estrutura do Programa

A estrutura do programa pode ser dividida em pré-processamento, programa

principal e pós-processamento. Cada etapa desta contém funções com finalidades

distintas, como mostrado na Ilustração 10.1.

78

Ilustração 10.1 – Diagrama representativo da estrutura do programa.

10.3 Pré-processamento

Nesta fase do programa, é definida a geometria do problema, condições de

contorno, carregamentos e o domínio discretizado. Todas estas informações são

armazenadas em um arquivo do tipo .mat (arquivo que armazena matrizes e vetores no

MatLab). Este arquivo será carregado posteriormente na fase de processamento pelas

funções probderivadatopol.m e probaliviotensao.m. Cada exemplo é gerado por

um script diferente.

10.4 Processamento

Nesta fase é iniciado o algoritmo de Otimização Topológica

(probderivadatopol.m) e em seguida o algoritmo de Alívio de Tensões

Pré-processamento

probderivadatopol.m

Processamento Pós-processamento

exemplo1.m

exemplo2.m

exemplo3.m

derivadatopologica.m

probaliviotensao.m


plotadadostopol.m

plotadadosalivio.m

Pré-processamento

probderivadatopol.m

Processamento Pós-processamento

exemplo1.m

exemplo2.m

exemplo3.m


probaliviotensao.m


plotadadostopol.m

plotadadosalivio.m

79

(probaliviotensao.m). A informação contida em cada iteração dos algoritmos é

armazenada em uma matriz denominada histórico. Esta matriz será utilizada no pós-

processamento para recuperar informações de iterações intermediárias. A função

derivadatopologica.m é utilizada por ambos os programas desta fase de

processamento. Daí surgiu a necessidade de encapsular este procedimento em uma

função à parte. Dentro deste procedimento se encontra o código dos elementos finitos

adaptado para a aplicação da derivada topológica e alívio de tensões.

10.5 Pós-processamento

Nesta fase é interpretado o resultado gerado pelas funções do processamento.

Uma análise gráfica do domínio final, da distribuição de tensão de Mises e da

convergência do algoritmo são apresentadas na tela. Isto permite um estudo detalhado

da variação de todos os parâmetros envolvidos no problema de Otimização e no alívio

de tensões. Há ainda a possibilidade de se recuperar iterações intermediárias através da

matriz histórico .

80

11 Conclusão

A derivada topológica é um tema em pleno desenvolvimento e esta dissertação

ajudou a difundir os conceitos desta nova teoria. A inovação proposta nesta dissertação

foi formulação da análise de sensibilidade topológica utilizada em conjunto com a teoria

da otimização evolutiva. Provou-se que a presente formulação é mais eficiente que a

formulação tradicional da análise de sensibilidade topológica, ou seja, fornece pontos

ótimos melhores do ponto de vista da engenharia.

Aplicando os conceitos apresentados aqui, possibilitou-se uma discretização

refinada da malha de elementos finitos, 7500 elementos, com um tempo de

processamento igual a 20 minutos na obtenção do ponto ótimo, para um computador

caseiro. O método dos elementos finitos se adeqüou perfeitamente à metodologia da

derivada topológica e denotou importante papel na discretização do domínio, o que

permitiu trabalhar com geometrias e condições de contorno complexas, de maneira que

não seria possível uma solução analítica dos exemplos demonstrados.

Apresentou-se uma comparação entre formulações distintas da sensibilidade da

função custo em relação às variáveis de projeto para a solução do problema da

minimização de energia de deformação com restrição no volume. Todas as formulações

chegaram a um resultado semelhante, apesar de serem originárias de diferentes teorias.

Além do fenômeno da elasticidade linear, muitos outros fenômenos podem ser

modelados de acordo com o conceito aplicado nesta dissertação. Como, por exemplo,

81

no fenômeno da condução e convecção de calor, na torsão de barras prismáticas e em

placas de Kirchhoff.

Uma sugestão para trabalhos futuros é a utilização de novas formulações de

função custo e restrições aplicando os conceitos da derivada topológica, e como

sugestão de estudo, tem-se o problema de minimizar o peso com restrições locais de

tensão. Por fim, objetiva-se aplicar esta teoria em fenômenos ainda não explorados com

o conceito de derivada topológica.

82

12 Bibliografia

[1] SOKOŁOWSKI, J., ZOCHOWSKI, A., On topological derivative in shape

optimization, Relatório de pesquisa 3170, INRIA - Lorraine,França, Maio,

1997.

[2] STAINKO, R., Advanced Multilevel Techniques to Topology Optimization, Tese de

Doutorado, Johannes Kepler Universit¨at, Linz, February 2006

[3] MICHELL, A.G.M. “The limit of economy of material in frame structures”.

Philosophical Magazine, v.8(6), pp. 589–597, 1904.

[4] CHENG ,G.D., OLHOFF ,N . “An investigation concerning optimal design of solid

elastic plates”. International Journal of Solids and Structures, v.17, pp.

305–323, 1981.

[5] BENDSØE, M.P., KIKUCHI, N., “Generating optimal topologies in structural

design using a homogenization method”, Computer Methods in Applied

Mechanics and Engineering, v.71(2), pp.197–224, 1988.

[6] BENDSØE, M.P., “Optimal shape design as a material distribution problem”,

Structural Optimization, v.1, pp.193–202, 1989.

[7] MURAT, F., SIMON, J. ,Sur le côntrole par un domaine géométrique, Tese de

Doutorado, Université Pierre et Marie Curie, Paris , França, 1976.

83

[8] CÉA, E.J., “Problems of Shape Optimal Design”. In: Proceedings of Optimization of

Distributed Parameters Structures. Iowa, Estados Unidos da América,

1981.

[9] ESCHENAUER, H.A., KOBELEV, V.V., SCHUMACHER, S., “Buble Method for

Topology and Shape Optimization of Structures”, Structural Optimization,

v.8, pp.42-51, 1994.

[10] CÉA, J., GARREOU, S., GUILLAUME, Ph., MASMOUDI, M., The Shape and

Topological Optimizations Connection, Relatório de Pesquisa, UFR MIG,

Université Paul Sabatier, Toulouse, França, 1998, também publicado na

revista Computer Methods in Applied Mechanics and Engineering, v.188,

pp.713-726, 2000.

[11] NOVOTNY, A.A., Análise de Sensibilidade Topológica, Tese de D.Sc,

Laboratório Nacional de Computação Científica, Petrópolis, RJ , Brasil,

Fevereiro de 2003.

[12] ZOLÉZIO, J.P., “The Material Derivative (or Speed) Method for Shape

Optimization”, In: Optimization of Distributed Parameters Structures, Iowa,

EUA, 1981.

[13] LUEMBERGER, D.G., Optimization by Vector Space Methods, John Wiley &

Sons, 1969.

[14] HAUG, E.J., CHOI, K.K., KOMBOV, V., Design of Sensitivity Analysis of

Structural Systems. Academic Press, 1986.

[15] GONZALES, O., STUART, A.M., Introduction to Continuum Mechanics,

Relatório de Pesquisa, University of Warwick , Conventry, Inglaterra, 2005.

[16] GURTIN, M.E., “An Introduction to Continuum Mechanics”, Mathematics in

Science and Engineering, v.158, Academic Press, 1981.

84

[17] LI, S. WANG, G., SAUER, R.A., “The Eshelby Tensors in a Finite Spherical

Domain: II. Applications in Homogenization”, Accepted for Publication in

ASME Journal of Applied Mechanics, Department of Civil and

Environmental Engineering, University of California, Berkeley, CA94720,

U.S.A.

[18] GÖTZ, T., Asymptotic Analysis, Relatório de Pesquisa, Technomathematics Group

Department of Mathematics University of Kaiserslautern, Kaiserslautern,

Germany, e-mail:[email protected] , 22 de Julho, 2005.

[19] NOVOTNY, A.A., FEIJÓ, R.A., TAROCO, E. , PADRA, C., “Topological

Sensitivity Analysis for Three-dimensional Linear Elasticity Problem” ,6th

World Congress on Structural and Multidisciplinary Optimization, Rio de

Janeiro, Brazil, 30 de Maio à 03 Junho 2005.

[20] BENDSØE, M.P., SIGMUND, O, Topology Optimization: Theory, Methods and

applications. Springer, Berlin, 2003.

[21] XIE, Y.M., STEVEN, G.P., “A simple evolutionary procedure for structural

optimization”, Comput. Struct., 49 (1993) 885–896.

[22] TANSKANEN, P., “The evolutionary structural optimization method: theoretical

aspects”, Computer Methods in Applied Mechanical Engineering, v. 191,

pp.5485–5498, 2002.

[23] RONG, J, TANG, G., LANG, Q.Q., YANG, Y., “A Topology Optimization

Method for Three-dimensional Continuum Structures”, 6th World Congress

on Structural and Multidisciplinary Optimization, Rio de Janeiro, Brazil, 30

de Maio à 03 Junho 2005.

[24] JANG, I.G., KWAK, B.M., “Evolutionary topology optimization for large-scale

problems using design space adjustment and refinement”, 6th World

85

Congress on Structural and Multidisciplinary Optimization, Rio de Janeiro,

Brazil, 30 de Maio à 03 Junho 2005.

[25] ESCHENAUER, H.A., OLHOFF , N., “Topology of continuum structures: a

review”, Applied Mechanics Review, v.54, pp.331-390, 2001.

[26] MARTINEZ, J.M., On the Theoretical Convergence Properties of the SIMP

method, Instituto de Matemática da Universidade de Campinas, CP 6065,

13081-970 Campinas, SP, Brasil, 4 de Março, 2004.

[27] BURGER, R.S., STAINKO, R., Phase-Field Relaxation of Topology Optimization

with Local Stress Constraints, Relatório de Pesquisa, Institut fur

Industriemathematik, Johannes Kepler Universitat, Altenbergerstr, 14 de

Janeiro, 2005.

[28] CHENG, G.D., GOU, X, “ǫ-relaxed approach in topology optimization”,

Structural and Multidisciplinary Optimization, 13:258–267, 1997.

[29] XIE, Y.M., STEVEN, G.P., “A simple evolutionary procedure for structural

optimization”, Computers & Structures, v.49(5), pp.885-896, 1993.

[30] GARREAU, S., GUILLAUME, P., MASMOUDI, M., The Topological Gradient,

Relatório de Pesquisa, UFRJ MIG, Université Paul Sabatier, Toulouse,

França, 1998.

[31] LABANOWSKI, JR.A., FANCELLO, E.A., Simp, Eso e Tsa: Uma Análise

Comparativa Demétodos de Otimização Topológica para Elasticidade 2d,

Relatório de Pesquisa, Grupo de Análise e Projetos Mecânicos,

GRANTE/EMC/UFSC, Laboratório Nacional de Computação Científica.

[32] TANSKANEN, P. “The evolutionary structural optimization method:theoretical

aspects”, Computer Methods in Applied Mechanics and Engineering, v.

191(2002), pp. 5485-5498, 5 de Agosto, 2002.

86

[33] ZIENKIEWICZ, The Finite Element Method in Engineering Science, McGraw-

Hill, 2a. ed., Londres, Inglaterra, 1971.

[34] STEVEN, G., LI, Q., QUERIN, O. “Some Thoughts on the Physics And

Mechanics of the Evolutionary Structural Optimization Process”, 6th

ASMO-UK/ISSMO International Conference on Engineering Design

Optimization, Oxford University, United Kingdon, 3 - 4 July 2006.

[35] HUGHES, T.J.R., The Finite Element Method, Englewood Cliffs, New Jersey,

Prentice-Hall, 1987.

[36] ANDERSON, E., BAI, Z., BISCHOF, C., BLACKFORD, S., DEMMEL, J.,

DONGARRA, J., DU CROZ, J., GREENBAUM, A., HAMMARLING, S.,

MCKENNEY ,A., SORENSEN D., LAPACK User's Guide,

(http://www.netlib.org/lapack/lug/lapack_lug.html), Third Edition, SIAM,

Philadelphia, 1999.

[37] OpenMP Application Program Interface, site:http://www.openmp.org/drupal/,

Version 2.5, May 2005.

[38] MASSERA, J.M.A., Uma nova técnica para otimização de estruturas de grande

porte, Tese de M.Sc., COPPE/UFRJ, Rio de Janeiro, RJ, Brasil, 2005.

Documents

UMA TÉCNICA PARA OTIMIZAÇÃO ESTRUTURAL MEDIANTE … · Prof. Luiz Eloy Vaz, D.Sc. RIO DE JANEIRO, RJ ... dos elementos finitos é utilizado para discretizar o domínio e obter