TECNICAS DE OTIMIZAÇ´ AO APLICADAS A PROBLEMAS DE … · elásticos, na hipótese de elasticidade infinitesimal sem atrito. Trata-se de um problema não linear devido à presença

TECNICAS DE OTIMIZACAO APLICADAS A PROBLEMAS DE CONTATO

ELASTICO

Gabriel Mario Guerra Bernada

DISSERTACAO SUBMETIDA AO CORPO DOCENTE DA COORDENACAO

DOS PROGRAMAS DE POS-GRADUACAO DE ENGENHARIA DA

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS

REQUISITOS NECESSARIOS PARA A OBTENCAO DO GRAU DE MESTRE

EM CIENCIAS EM ENGENHARIA MECANICA.

Aprovada por:

Prof. Jose Herskovits Norman, Dr. Ing.

Prof. Marcelo Amorim Savi, D.Sc.

Prof. Angela Cristina Cardoso de Souza, D.Sc.

RIO DE JANEIRO, RJ - BRASIL

OUTUBRO DE 2006

BERNADA, GABRIEL MARIO GUERRA

Tecnicas de otimizacao aplicadas a

problemas de contato elastico. [Rio de

Janeiro] 2006

IX, 79p. 29,7 cm (COPPE/UFRJ, M.Sc.,

Engenharia Mecanica, 2006)

Dissertacao - Universidade Federal do Rio

de Janeiro, COPPE

1. Otimizacao nao linear.

1. Mecanica de solidos.

I. COPPE/UFRJ II. Tıtulo (serie)

ii

A Mario Gaston Bernada

iii

Agradecimentos

Ao professor Herskovits, pela orientacao.

A Sandro, Evandro e Carmen Nilda por os valiosos aportes e observacoes.

Ao os colegas e amigos do laboratorio Optimize: Alfredo, Paulo, Henry, Moises,

Miguel, Veranisse, Vania, Denise, Rossane e Silvio.

Ao corpo docente e pessoal administrativo do Programa de Engenharia Mecanica

da COPPE.

Ao programa PEC-PG da CAPES pelo apoio financeiro.

Aos professores do Instituto de Ingenierıa Mecanica y Produccion Industrial da

Faculdade de Engenharia do Uruguai pelo incentivo.

A minha avo, mae e irma pelo grande carinho e apoio constante a distancia.

iv

Resumo da dissertacao apresentada COPPE/UFRJ como parte dos requisitos

necessarios para a obtencao do grau de Mestre em Ciencias (M.Sc.)

TECNICAS DE OTIMIZACAO APLICADAS A PROBLEMAS DE CONTATO

ELASTICO


Outubro/2006

Orientador: Jose Herskovits Norman

Programa: Engenharia Mecanica

Este trabalho apresenta a formulacao do problema de contato entre corpos

elasticos, na hipotese de elasticidade infinitesimal sem atrito. Trata-se de

um problema nao linear devido a presenca de restricoes unilaterais de nao

interpenetracao. Por conseguinte, utilizando conceitos de otimizacao, pode-se

formular como um problema de minimizacao com restricoes de desigualdade, ou

de forma equivalente como um problema de complementaridade na variavel dual.

Devido ao carater localizado do contato, e conveniente adotar tecnicas que

trabalhem somente com a regiao de contato potencial como variaveis de otimizacao.

Com essa finalidade, foi desenvolvida uma interface entre os algoritmos e um

codigo comercial de elementos finitos, a qual permite o desacoplamento da analise

estrutural.

A presente pesquisa compara o desempenho dos algoritmos de ambas

formulacoes, e explora as vantagens da utilizacao de um codigo comercial como

ferramenta auxiliar na analise de problemas de grande porte. Em razao disso,

sao apresentados alguns exemplos de contato entre corpos bidimensionais e

tridimensionais no capıtulo final.

v

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

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

OPTIMIZATION TECHNIQUES APPLIED FOR ELASTIC CONTACT

PROBLEMS


October/2006

Advisor: Jose Herskovits Norman

Department: Mechanical Engineering

This work presents the formulation of the contact problem between elastic bodies,

in the hypothesis of infinitesimal elasticity without friction. This is a non linear

problem due to the presence of unilateral constraints of non interpenetrability.

Consequently, using concepts of optimization, a minimization problem with

inequality constraints, or in equivalent way as a complementarity problem in the

dual variable can be formulated.

Due to the local character of the contact, it is convenient to adopt techniques

that only work with this potential contact area as optimization variables. Thus, an

interface between the algorithms and a commercial code of finite elements analysis,

which allows the uncoupling of structural analysis, was developed.

The present research compares the acting of the algorithms of both formulations,

and explores the advantages of the commercial code as an auxiliary tool of the

optimization algorithms for the analysis of large problems. Then, examples are

presented among two-dimensional and three-dimensional bodies in the final chapter.

vi

Sumario

1 Introducao 1

1.1 Conceitos preliminares . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Objetivos e organizacao deste trabalho . . . . . . . . . . . . . . . . . 5

2 Equilıbrio em elasticidade linear 7

2.1 Formulacao forte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Formulacao fraca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.3 Princıpio de mınima energia potencial total . . . . . . . . . . . . . . . 10

2.4 Solucao aproximada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Apresentacao do problema de contato 14

3.1 Problema de Signorini . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2 Contato entre corpos elasticos . . . . . . . . . . . . . . . . . . . . . . 18

3.2.1 Geometria do contato . . . . . . . . . . . . . . . . . . . . . . . 18

4 Metodos de resolucao 21

4.1 Metodo 1 - Problema de otimizacao . . . . . . . . . . . . . . . . . . . 21

4.2 Descricao do algoritmo (FDIPA) . . . . . . . . . . . . . . . . . . . . . 23

4.2.1 Aplicacao ao problema de contato . . . . . . . . . . . . . . . . 25

4.2.2 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.3 Metodo 2 - Problema de complementaridade . . . . . . . . . . . . . . 31

4.4 Descricao do algoritmo (ALG-NCP) . . . . . . . . . . . . . . . . . . . 33

4.4.1 Aplicacao ao problema de contato . . . . . . . . . . . . . . . . 36

4.4.2 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

vii

5 Interface com ABAQUS 40

5.1 Descricao do programa . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.1.1 Acessibilidade e interacao . . . . . . . . . . . . . . . . . . . . 42

5.2 Descricao da interface . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

5.3 Funcionamento da interface . . . . . . . . . . . . . . . . . . . . . . . 43

6 Resultados numericos 45

6.1 Exemplo 1 - Semicilindro (2D) . . . . . . . . . . . . . . . . . . . . . . 46

6.2 Exemplo 2 - Mola (3D) . . . . . . . . . . . . . . . . . . . . . . . . . . 50

6.3 Exemplo 3 - Pinca (3D) . . . . . . . . . . . . . . . . . . . . . . . . . 54

6.4 Exemplo 4 - Pinhao e cremalheira (2D) . . . . . . . . . . . . . . . . . 58

7 Consideracoes finais 62

7.1 Conclusoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

7.2 Sugestoes para trabalhos futuros . . . . . . . . . . . . . . . . . . . . 63

Referencias Bibliograficas 65

A Algumas consideracoes sobre otimizacao 67

B Metodo de Newton Raphson 70

C Codigos da interface 72

viii

Lista de sımbolos

IRn Vetores reais de dimensao n

IRm×n Matrizes reais de dimensoes m× n

x Vetor de variaveis primais ou variaveis de projeto

u Vetor de deslocamentos

σ Tensor de tensoes

ε Tensor de deformacoes infinitesimais

b Campo de forcas de volume

f Campo de forcas de superfıcie

F Campo de forcas externas

F Funcao de complementaridade

Φ Funcao potencial energıa elastica, Φ : IRn → IR

f Funcao objetivo. f : IRn → IR

g Funcao das restricoes de desigualdade. g : IRn → IRm

∇f Gradiente da funcao f(x), ∇f : IRn → IRn

∇g Gradiente da funcao g(x), ∇g : IRn → IRm×n

Ω Regiao viavel, conjunto de pontos que satisfazem as restricoes

λ Multiplicadores de Lagrange das restricoes de desigualdade, λ ∈ IRm

L Funcao lagrangeana. L(u, λ) : IRn × IRm × IRp → IR

Π Funcao energıa total de deformacao Π(u) : IRn → IR

H Matriz hessiana. H : IRn × IRm × IRp → IRn×n

K Matriz de rigidez do solido k

t Passo, t ∈ IR

d, d0, d1 Direcoes de busca. d ∈ IRn, d0 ∈ IRm, d1 ∈ IRp

Λk, Gk Matriz diagonal, diag(λk), diag(gk)

• Produto de Hadamard

ix

Lista de Figuras

2.1 Solido em equilibro . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3.1 Condicoes geometricas de contato . . . . . . . . . . . . . . . . . . . . 14

3.2 Problema de Signorini . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.3 Malhas incompatıveis . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.4 Situacao de contato no - segmento . . . . . . . . . . . . . . . . . . . . 19

3.5 Situacao de contato no - no . . . . . . . . . . . . . . . . . . . . . . . 20

3.6 Campo de deslocamentos incompatıveis na configuracao no - no . . . 20

4.1 Direcao de busca do algoritmo FDIPA . . . . . . . . . . . . . . . . . 24

4.2 Direcao de busca do algoritmo ALG-NCP . . . . . . . . . . . . . . . 34

5.1 Esquema da Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.2 Diagrama da interacao Matlab - ABAQUS . . . . . . . . . . . . . . . 44

6.1 Modelo 1 condicoes de contorno . . . . . . . . . . . . . . . . . . . . . 46

6.2 Configuracao deformada com tensoes de Von-Mises . . . . . . . . . . 47

6.3 Malha deformada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

6.4 Curva pressao na regiao de contato . . . . . . . . . . . . . . . . . . . 47

6.5 Criterio de parada . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48


6.7 Configuracao deformada com tensoes de Von-Mises . . . . . . . . . . 51

6.8 Nos da face do modelo em contato na superfıcie analıtica . . . . . . . 51



6.11 Geometria da malha e tensoes de Von-Mises . . . . . . . . . . . . . . 55

x

6.12 Ampliacao da zona de contato . . . . . . . . . . . . . . . . . . . . . . 55



6.15 Geometria da malha . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

6.16 Tensoes de Von-Mises nos elementos . . . . . . . . . . . . . . . . . . . 59


B.1 Direcao de busca do algoritmo de Newton . . . . . . . . . . . . . . . 71

xi

Capıtulo 1

Introducao

O estudo do contato entre corpos deformaveis apresenta grande interesse na

engenharia devido a que os esforcos sao transmitidos nesta condicao na maioria dos

componentes mecanicos. Quando o interesse na distribuicao de tensoes encontra-se

em regioes distantes a regiao de contato, a aplicacao do principio de Saint-Venant

e suficiente. Nele, afirma-se que a distribuicao longe do ponto de aplicacao das

forcas depende somente da resultante. Deste modo, as forcas de contato podem

ser simplificadas a forcas de superfıcie, constituindo-se num dado do problema de

equilıbrio. Porem, existem casos onde devido as grandes tensoes geradas nas areas de

contato, torna-se indispensavel conhecer sua distribuicao na regiao. Como exemplos

podem-se citar o desenho de engrenagens e eixos, onde tensoes de contato cıclicas sao

responsaveis pelas falhas por fadiga ou fenomenos mais complexos como os processos

de conformacao metalica ou colisoes. A partir destes exemplos, pode ser feita uma

classificacao para a analise do contato, a saber:

- Contato intencional: previamente planejado, procura maximizar a eficiencia do

contato. Este e o caso do contato entre a matriz de conformacao mecanica e a peca

a ser conformada.

- Contato acidental: ocorre em situacoes nao previstas e procura minimizar os

efeitos, como por exemplo o impacto de veıculos.

Como pode-se observar, sao inumeras as aplicacoes onde intervem o fenomeno

de contato, assim como as dificuldades que apresentam-se no tratamento analıtico,

constituindose portanto num tema ativo de pesquisa.

1

Os primeiros estudos do problema de contato aparecem na literatura por volta

do seculo XIX, com o objetivo de desenvolver uma teoria para o impacto. Entre

os autores encontram-se nomes como POISSON, VOIGHT e SAINT-VENANT

responsaveis pelas primeiras formulacoes para esta teoria.

Em (1882) Heinrich HERTZ considerou o equilıbrio de dois corpos elasticos

em contato superficial. Seu interesse surgiu da analise da deformacao das

lentes submetidas a acao de forcas de contato e como estas poderiam gerar

interferencia otica. Assim, por analogia com a distribuicao de pressoes produzida por

deslocamentos compatıveis com uma area elıptica, ele obteve com sucesso formulas

analıticas para serem aplicadas em alguns problemas especiais, as quais ate hoje sao

utilizadas, TIMOSHENKO [17], JOHNSON [11]. Posteriormente, SIGNORINI em

(1933) apresento seu primer estudo do equilıbrio entre corpos linearmente elasticos

com contato no qual trabalho ate sua versao definitiva em (1959). Durante a primeira

metade do seculo XX a escola russa trouxe importantes contribuicoes, BELIAEV

e DINNIK verificaram experimentalmente os resultados obtidos por HERTZ, e

GALIN e SHTERMAN propuseram novos metodos para analise do contato.

A primeira analise rigorosa do problema de SIGNORINI, foi apresentada por

FICHERA em (1963), realizando valiosas contribuicoes teoricas sobre existencia e

unicidade das solucoes alcancadas a partir da formulacao variacional, este e outros

detalhes sobre a evolucao historica podem-se ver em KIKUCHI [12]. Daı em diante

muitos avancos foram efetuados, tanto na formulacao quanto na resolucao. Devido

principalmente a evolucao dos metodos numericos, e em particular, ao metodo dos

elementos finitos. Nesta tecnica, as variaveis e funcoes que definem o problema

derivam da discretizacao do corpo, obtendo desta forma solucoes aproximadas. No

entanto, com estes avancos, novas dificuldades apresentarem-se, tanto do ponto

de vista matematico quanto de sua resolucao computacional. Constituindo-se na

principal motivacao para que durante os ultimos anos um esforco consideravel tenha

sido direcionado a analise e solucao destes problemas.

Hoje em dia, como consequencia do crescente poder de processamento dos

computadores, e possıvel atingir maior sofisticacao e detalhe nos modelos de

contato. Em virtude disto, encontram-se disponıveis numerosos pacotes comerciais

2

de analise por elementos finitos (FEA), que incluem em seu codigo ferramentas

para o tratamento de alguns problemas de contato. Neste sentido, este trabalho

procura avaliar novas ferramentas aplicadas a este tipo de problemas, utilizando

como ferramenta auxiliar um software (FEA) comercial para analise estrutural.

1.1 Conceitos preliminares

Quando solidos estao em contato, existe uma restricao geometrica de nao

penetrabilidade que deve ser satisfeita na interface. Como nao se conhece a regiao

de contato a priori, nao e possıvel determinar deslocamentos nem tensoes, ou seja,

nao se podem estabelecer as condicoes de contorno. Como consequencia, mesmo

em elasticidade linear na hipotese de deformacoes infinitesimais, o problema de

valor de contorno deixa de ser linear. No caso geral, ha dois fenomenos fısicos

a serem considerados nesta interface: o comportamento na direcao normal e o

comportamento na direcao tangencial. Na direcao normal, a formulacao da condicao

de nao penetrabilidade como restricao geometrica e o procedimento adotado em

muitos trabalhos como KIKUCHI, SIMO, WRIGGERS [12] [4] [14]. Na direcao

tangencial, o comportamento e mais complexo devido a existencia de atrito. Para

modelar este fenomeno existem muitas aproximacoes, a mais frequente e baseada

na lei classica de Coulomb. Nesta, o movimento na direcao tangencial e dividido

em duas acoes. A primeira acao se denomina condicao de aderencia perfeita

- Stick - onde nao ocorre deslocamento relativo tangencial. A segunda acao e

denominada de deslizamento - Slip -, e existe um movimento tangencial relativo

na interface de contato. Outra forma de modelar o atrito e atraves de analogias com

o fenomeno de plasticidade. Nele, as equacoes constitutivas tem como ideia basica

dividir o deslocamento tangencial em uma parte elastica reversıvel e outra plastica

irreversıvel, constituindo assim um problema nao conservativo BARBOZA [7].

Neste trabalho nao sera considerado o efeito do atrito, portanto, o problema

de contato e conservativo e, como sera visto mas adiante pode-se considerar como

um problema de minimizacao de um funcional sujeito a restricoes de desigualdade.

Portanto, recai-se num problema de programacao matematica que podera ser tratado

3

atraves das tecnicas de otimizacao usuais. Entre estas tecnicas pode-se citar: o

metodo de programacao quadratica ou o metodo de penalidade, este ultimo consiste

em transformar o problema num outro de minimizacao sem restricoes atraves da

inclusao das restricoes de impenetrabilidade ao funcional. Outra alternativa e

descrita no trabalho de BJORKMAN, KLARBRIG e SJODIN [6], onde mostra como

um problema quadratico de minimizacao com restricoes de desigualdade, pode ser

reduzido a um problema de complementaridade linear (PLC) nas variaveis duais,

atraves do estabelecimento das condicoes de otimalidade (condicoes de Karush-

Kuhn-Tucker) e posterior eliminacao das variaveis primais. Para resolver problemas

de complementaridade existem uma grande variedade de tecnicas as quais serao

vistas em capıtulos posteriores.

Finalmente, outro ponto importante a ser considerado e quando apresenta-se

contato entre dois ou mais solidos deformaveis. Neste caso, o calculo das restricoes

geometricas nao e trivial, pois e necessario estabelecer uma estrategia para medir

a distancia entre corpos discretizados. Este ponto e mais importante em caso de

trabalhar com grandes deformacoes, levando em conta que um ponto da superfıcie de

um corpo pode entrar em contato com qualquer regiao da superfıcie de outro corpo.

Quando o contato e entre um corpo deformavel e um corpo rıgido, o problema e mais

simples ja que o corpo rıgido pode ser descrito por funcoes analıticas, o que facilita

o calculo das restricoes, KLARBRING [2]. Uma revisao detalhada do estado da

arte do problema de contato, assim como das vantagens e desvantagens de distintas

estrategias para o calculo das restricoes de impenetrabilidade pode-se ver no trabalho

de ARORA [3].

A seguir, sao apresentados os objetivos da presente dissertacao junto a

consideracoes adotadas para abordagem do problema de contato.

4

1.2 Objetivos e organizacao deste trabalho

Este trabalho tem como objetivo apresentar o problema de contato sem atrito

na hipotese de pequenas deformacoes. Para sua resolucao, serao adotadas duas

formulacoes alternativas. A primeira formulacao define o problema de minimizacao

na variavel primal, e a segunda define um problema de complementaridade linear

na variavel dual.

Para a resolucao do problema de minimizacao sera adotado uma variante do

metodo de Newton da famılia dos metodos de Lagrange. Esta variante proposta por

HERSKOVITS [8] chama-se ’FDIPA’ (Feasible direction interior point algorithm).

A tecnica baseia-se no calculo da direcao de descida a partir de uma iteracao

de Newton nas condicoes de Karush-Kuhn-Tucker (KKT) de primeira ordem,

seguida de uma pequena deflexao para garantir uma direcao viavel no passo.

Este procedimento iterativo gera, a partir de um ponto viavel u0, uma sequencia

de pontos viaveis convergentes para a solucao do problema, onde as restricoes

sempre sao satisfeitas, o que representa uma vantagem se comparado ao metodo de

penalidade. Para a resolucao do problema de complementaridade, sera utilizado um

novo algoritmo chamado(ALG-NCP) baseado na ideia do (FDIPA). Este algoritmo

gera uma sequencia de pontos interiores a regiao viavel, minimizando uma funcao

potencial associada ao problema.

Devido ao carater localizado do contato, as restricoes envolvem apenas os nos de

contato potencial, o que sugere trabalhar somente com o nos da regiao de contato

potencial como variaveis de otimizacao. Com esse objetivo, em geral adotam-se

tecnicas de sub-estruturacao ou condensacao estatica da matriz de rigidez. Neste

trabalho foi possıvel atingir esse objetivo atraves do desacoplamento do calculo

estrutural, permitindo resolver o problema sem a necessidade de ter acesso a matriz

de rigidez.

Em razao disto, sera desenvolvida uma interface entre os algoritmos apresentados

e um aplicativo comercial de elementos finitos ABAQUS. Isto permitira aos

algoritmos trabalhar somente com os nos da regiao de contato potencial como

variaveis de otimizacao, diminuindo assim o custo computacional. Outra vantagem

5

desta estrategia e permitir atingir problemas de maior porte com menores recursos

computacionais, assim como utilizar plataformas computacionais diferenciadas

segundo a necessidade. Por conseguinte, procura-se explorar as vantagens do uso

desta interface como ferramenta auxilar dos algoritmos, nao so no processamento,

mas tambem no pre-processamento e pos-processamento.

Com esta finalidade, sao apresentados no final deste trabalho alguns exemplos

bidimensionais e tridimensionais de contato com os quais a interface foi calibrada.

A partir destes exemplos tambem sera feita a comparacao do desempenho dos

algoritmos apresentados em termos de tempo e resultados numericos obtidos.

Este trabalho encontra-se organizado da seguinte forma:

• Capıtulo 1. Faz uma introducao do problema de contato e descreve o objetivo

do trabalho.

• Capıtulo 2. Apresenta conceitos de equilıbrio em elasticidade linear junto as

formulacoes forte e fraca.

• Capıtulo 3. Trata o problema geral de contato apresentando duas formulacoes

do mesmo.

• Capıtulo 4. Mostra esquemas numericos para a solucao das formulacoes e um

resumo dos algoritmos.

• Capıtulo 5. Apresenta a interface desenvolvida com o ABAQUS e sua

aplicacao.

• Capıtulo 6. Mostra os resultados numericos obtidos nos exemplos.

• Capıtulo 7. Finaliza a pesquisa, expondo as conclusoes e sugestoes para

trabalhos futuros.

6

Capıtulo 2

Equilıbrio em elasticidade linear

Previamente a analise do problema de contato, efetuamos algumas consideracoes a

respeito do problema de equilıbrio geral num solido elastico. Com este objetivo, sao

apresentadas as formulacoes forte (diferencial) e fraca (variacional) do problema

de valor de contorno (PVC) associado a hipotese de pequenos deslocamentos e com

condicoes de contorno bilaterais, ou seja, igualdades.

2.1 Formulacao forte

Seja um solido contınuo que na sua configuracao indeformada ocupa um subconjunto

do espaco Ω , aberto, limitado no contorno Γ, e simplesmente conexo como na

(Figura 2.1).

Figura 2.1: Solido em equilibro

7

Assume-se que o contorno Γ de Ω e dividido em regioes,

- Γu onde os deslocamentos sao prescritos,

- Γt onde atuam forcas de superfıcie f ,

que verificam,

Γt ∪ Γu = Γ, (2.1)

Γt ∩ Γu = ∅. (2.2)

Para simplificar considera-se o corpo fixo em Γu,

u(x) = 0 ∀ x ∈ Γu. (2.3)

O campo de tensoes internas σ do corpo estara em equilıbro com o carregamento

externo (b, f) se sao verificadas as seguintes equacoes,

div(σ) + b = 0 em Ω (2.4)

σ.n = f em Γt (2.5)

onde b representa o campo das forcas de volume em Ω, f o campo das forcas de

superfıcie em Γt , n e um vetor unitario normal exterior a esta superfıcie e σ o tensor

de tensoes.

Considerando-se as tensoes como uma funcao do tensor de deformacoes

infinitesimais ε, pode-se admitir a existencia de uma relacao constitutiva da forma:

σ = σ(ε(u)) (2.6)

Portanto, o problema de equilıbrio de um solido, Fig. 2.1, consiste em achar o

par (u, σ) tal que,

u satisfaca as condicoes de contorno geometricas (Eq. 2.3),

σ esteja em equilıbrio com (b, f) (Eqs. 2.4 e 2.5),

(σ, u) estejam ligados pela relacao constitutiva (Eq. 2.6).

8

2.2 Formulacao fraca

Considere agora um espaco vetorial V formado pelo campo de deslocamentos com

regularidade suficiente para que as operacoes a seguir facam sentido. Nele, considere

o subespaco K formado pelos deslocamentos cinematicamente admissıveis,

K = v ∈ V : v(x) = 0 ∀ x ∈ Γu. (2.7)

Multiplicando-se escalarmente a Eq. 2.4 pelo campo de deslocamentos v ∈ V

tem-se,

div(σ) · v + b · v = 0, (2.8)

levando em conta a seguinte propriedade,

div(σ · v) = div(σ) · v + σ · ∇v, (2.9)

e integrando no domınio Ω, obtem-se a seguinte expressao,

−

∫

Ω

σ · ∇v dΩ +

∫

Ω

div(σ · v) dΩ +

∫

Ω

b · v dΩ = 0, (2.10)

logo, aplicando o teorema da divergencia no segundo termo da equacao integral,

obtem-se,∫

Ω

σ · ∇v dΩ =

∫

Γt

f · v dΓ +

∫

Ω

b · v dΩ. (2.11)

Devido a hipotese de deformacoes infinitesimais, o tensor de deformacoes pode ser

dado pela parte simetrica do gradiente do campo de deslocamentos,

ε(v) = (∇v)s =1

2(∇v + ∇vt). (2.12)

Assim, substituindo a (Eq. 2.12) na (Eq. 2.11), obtem-se a seguinte expressao

conhecida como o princıpio dos trabalhos virtuais,

∫

Ω

σ · ε(v) dΩ =

∫

Ω

b · v dΩ +

∫

Γt

f · v dΓ, (2.13)

expressao que e satisfeita para todo campo de deslocamentos v ∈ K, desde que σ

esteja em equilıbrio com as forcas externas (b, f).

Por conseguinte, esta equacao representa o equilıbrio global e constitui a

caracterizacao fraca ou integral do equilıbrio. Levando em conta a relacao

9

constitutiva do material (Eq. 2.6) e o princıpio dos trabalhos virtuais (Eq. 2.13);

aplicando v ∈ V e a solucao do problema de valor de contorno (u = v), obtem-se a

expressao,

∫

Ω

σ(ε(u)) ·ε(v−u) dΩ =

∫

Ω

b · (v−u) dΩ+

∫

Γt

f · (v−u) dΓ ∀v ∈ K, (2.14)

que corresponde a equacao variacional do equilıbrio do solido da (Fig.2.1) com

solucao em u.

2.3 Princıpio de mınima energia potencial total

Dada uma equacao constitutiva do tipo elastica linear, pode-se introduzir uma

funcao densidade de energia de deformacao Φ quadratica nas componentes do tensor

de deformacoes infinitesimais ε,

Φ(ε) :=1

2ε

t · C · ε =1

2εij · Cijkl · εkl, (2.15)

onde C e o tensor constitutivo de elasticidade linear de quarta ordem, que satisfaz

as propriedades usuais de simetria e elipticidade,

C(u) A ·B = A · C(u)B, (2.16)

C(u) A ·B ≥ α A · A, (2.17)

sendo α positivo , u ∈ Ω e (A e B) tensores simetricos de segunda ordem.

Portanto, escrevendo a constitutiva como,

σ(u) =∂Φ(ε)

∂ε= C · ε(u) , σij = Cijkl · εkl (2.18)

a equacao variacional (Eq. 2.14) fica,

∫

Ω

C · ε(u) · ε(v − u)dΩ =

∫

Ω

b · (v − u)dΩ +

∫

Γf

f · (v − u)dΓ ∀v ∈ K. (2.19)

assim, levando em conta que Φ e convexa e ε e linear obtem-se,

Φ(ε(v)) − Φ(ε(u)) ≥ Cε(u) · ε(v − u). (2.20)

10

Finalmente, introduzindo o funcional da energia potencial total do corpo,

Π(v) :=1

2

∫

Ω

Cε(v) · ε(v) dΩ −

∫

Ω

b · v dΩ −

∫

Γt

f · v dΓ, (2.21)

e substituindo a (Eq. 2.20) na (Eq. 2.19), chega-se a formulacao equivalente do

problema de equilıbrio conhecida como ”Principio de mınima energia potencial

total”,

u ∈ K : Π(u) ≤ Π(v) ∀v ∈ K. (2.22)

que expressa a condicao de equilıbrio de um corpo elastico como aquela que minimiza

o funcional da energia potencial total do sistema Π(u).

De uma forma geral, o problema do equilıbrio de uma estrutura elastica linear

com restricoes de igualdade ou bilaterais pode ser escrito em sua formulacao fraca

como,

a(u, u− v) = l(v − u) ∀ v ∈ V (2.23)

onde, a(·, ·) a forma bilinear simetrica, positiva e contınua que representa o trabalho

das forcas internas associadas a um deslocamento virtual,

a(u, v) =

∫

Ω

C ε(u) · ε(v) dΩ, (2.24)

e l(·) e o funcional linear que representa o trabalho virtual das forcas externas,

l(v) =

∫

Γt

f · v +

∫

Ω

b · v dΩ. (2.25)

Em virtude disto, a energia potencial total do sistema pode ser expressa em sua

forma variacional como,

Π(v) =1

2a(v, v) − l(v). (2.26)

11

2.4 Solucao aproximada

Como foi visto, o problema de equilıbrio formulado a partir da equacao variacional

requer a minimizacao de um funcional em dimensao infinita. A fim de obter uma

solucao aproximada deste problema, reformula-se num espaco de dimensao finita

atraves do ”Metodo dos Elementos Finitos” (MEF). Nesta tecnica, o domınio e

discretizado com elementos finitos onde sao geradas aproximacoes da solucao a

partir de um conjunto finito de ”funcoes-base”denominadas funcoes de interpolacao.

Assim, a primeira aproximacao e obtida a partir da substituicao do domınio original

Ω por outro Ωh formado por uma geometria de elementos finitos ZIENCKIEWICS

[13].

Denotando as funcoes base como φh, o espaco de aproximacao gerado fica definido

por,

Vh = vh : vh =n

∑

i=1

βiφhi (x), βi ∈ R, x ∈ Ωh, (2.27)

onde a aproximacao Vh, se constitui num subespaco interior ao espaco original do

problema ou seja, Vh ⊂ V.

Neste espaco, tanto forma a bilinear a(·, ·) quanto o funcional linear l(·) se

aproximam por,

ah(u, v) =

∫

Ωh

C ε(u) · ε(v) dΩ, (2.28)

lh(v) =

∫

Γht

f · v dΓ +

∫

Ωh

b. · v dΩ (2.29)

Considerando estas formas em Vh resulta,

ah(uh, vh) =

n∑

i=1

n∑

j=1

αiβj

∫

Ωh

Cε(φhi ) · ε(φh

j ) dΩ, (2.30)

lh(vh) =

n∑

i=1

βi[

∫

Γht

f · φhj dΓ +

∫

Ωh

b · φhj dΩ], (2.31)

onde os termos,

Kij =

∫ h

Ω

C ε(φhi ) · ε(φh

j )dΩ, (2.32)

12

estao associados a matriz de rigidez, enquanto os elementos fj, representam o vetor

de forcas nodais equivalentes ao carregamento (b, f),

fj =

∫

Γhf

f · φhj dΓ +

∫

Ωh

b · φhj dΩ. (2.33)

Como resultado desta aproximacao, o funcional da energia potencial total Π em sua

forma discreta se representa como,

Πh(vh) =1

2

n∑

i=1

n∑

j=1

βiβjKij −n

∑

j=1

βjfj, (2.34)

ou em notacao compacta,

Π(u) =1

2utKu− f tu, (2.35)

onde K e a matriz de rigidez global do sistema, u o vetor de deslocamentos e f e

o vetor de esforcos nodais equivalentes, incluindo esforcos de superfıcie e de corpo.

Portanto, a condicao necessaria de mınimo da energia potencial e,

∂Π(u)

∂u= 0 =⇒ Ku− f = 0. (2.36)

Esta condicao conduz a um sistema linear de equacoes, onde se os movimentos rıgidos

do corpo sao eliminados a positividade da matriz K do sistema e garantida. Logo

como Π(u) e uma funcao quadratica, a funcao e estritamente convexa e a solucao

achada unica representando o mınimo global do funcional, isto pode-se ver com mais

detalhe em LUENBERGER [5].

13

Capıtulo 3

Apresentacao do problema de

contato

No capıtulo anterior, foram apresentadas as formulacoes forte e fraca do equilıbro

em solidos elasticos sob a hipotese de deformacoes infinitesimais. Agora, apresenta-

se a formulacao geral do problema de contato. Neste sentido, novas restricoes

aos deslocamentos do problema de equilıbrio sao adicionadas como desigualdades

representando as condicoes de impenetrabilidade na fronteira.

Assim sendo, consideram-se dois corpos A e B em contato eventual nas regioes

ΓAc e ΓB

c e suas normais externas nA e nB, como na Fig.3.1.

Figura 3.1: Condicoes geometricas de contato

14

Para simplificar o desenvolvimento posterior e a notacao, supoe-se nA = nB = n

portanto, se as respectivas normais se confundem, pode-se escrever a condicao

linearizada de nao interpenetrabilidade entre os corpos A e B como:

s(x) − u(x) · n(x) ≥ 0 ∀ x ∈ Γc onde Γc = ΓAc = ΓB

c (3.1)

Assim s(x) e a funcao folga na direcao normal n as superfıcies, entre os pontos

de ΓAc e ΓB

c das regioes de contato potencial. Denotando com o sub-ındice (n) a

normal do campo vetorial de deslocamentos, tem-se:

un − s ≤ 0 e un = uAn + uB

n (3.2)

sendo rn a reacao no ponto de contato da regiao Γc e, portanto, as seguintes condicoes

devem ser satisfeitas:

(un − s) < 0 =⇒ rn = 0 (3.3)

(un − s) = 0 =⇒ rn < 0 (3.4)

As (Eq. 3.3 e Eq. 3.4) podem ser reescritas como:

(un − s) ≤ 0 (3.5)

rn ≤ 0 (3.6)

(un − s) · rn = 0 (3.7)

sendo a (Eq. 3.7) conhecida como relacao de complementaridade entre (un − s) e rn

do problema de contato.

15

3.1 Problema de Signorini

Este problema descreve o contato sem atrito de um corpo elastico linear com um

corpo rıgido, sendo base para muitas generalizacoes. Considerando a configuracao

indeformada que ocupa o domınio Ω, o problema de Signorini divide em tres areas

o contorno Γ, como mostra a Figura 3.2, a saber:

Γu onde os deslocamentos prescritos sao iguais a zero;

Γt onde as forcas de superfıcie sao prescritas;

Γc onde os nos podem ter contato;

Figura 3.2: Problema de Signorini

Γu ∪ Γt ∪ Γc = Γ (3.8)

Γu ∩ Γt = Γu ∩ Γc = Γt ∩ Γc = ∅ (3.9)

A partir das novas condicoes de fronteira, o equilıbrio pode ser formulado na sua

forma forte, adicionando-se as novas condicoes descritas nas (Eqs. 3.6, 3.7 e 3.7).

16

Assim temos:

divσ(u) + b = 0 em Ω (3.10)

σ(u) · n = f em Γt (3.11)

u = 0 em Γu (3.12)

(un − s) ≤ 0

rn = [σ(u) · n] · n ≤ 0

(un − s) · rn = 0

em Γc (3.13)

Definindo agora o conjunto de deslocamentos cinematicamente admissıveis para

este problema como,

Definicao 3.1.1 Campo de deslocamentos admissıveis

Σ = u ∈ U | un ≤ s em Γc (3.14)

onde U e um espaco de deslocamentos suficientemente regulares.

Considerando o contato como um caso particular do problema de equilıbrio,

temos que sua solucao u minimiza a energia potencial total Π(u) (Eq. 2.35) num

conjunto de deslocamentos cinematicamente admissıveis ver KIKUCHI [12].

Partindo da formulacao variacional do equilıbrio, chega-se de igualmente a

expressao da energia total de deformacao do solido, como a diferenca entre a energia

elastica de deformacao e a energia potencial das forcas externas. Portanto, o

problema de contato resume-se a minimizacao da energia de deformacao no campo

de deslocamentos admissıveis Σ.

Ou seja, se expressa como,

u ∈ Σ : Π(u) ≤ Π(v) ∀v ∈ Σ, (3.15)

onde a energia total de deformacao e representada na sua forma discreta por (Eq.

2.35). Para escrever o campo de deslocamentos admisıveis como restricoes de

desigualdade define-se a funcao gap como,

17

Definicao 3.1.2 Funcao ”gap”

g(u) = un − s (3.16)

onde s e a distancia inicial entre os pares de nos correspondentes na direcao n.

onde g(u) ≤ 0 no problema de otimizacao associado.

3.2 Contato entre corpos elasticos

Generalizando os resultados vistos anteriormente, em caso de contato entre dois

corpos A e B, minimiza-se a soma das energias potenciais totais,

ΠA(u) =1

2ut

AKAuA − f tAuA, (3.17)

ΠB(u) =1

2ut

BKBuB − f tBuB. (3.18)

Entao, pode-se escrever a energia potencial total dos corpos A e B numa

unica equacao, obtendo desta forma grandes vantagens do ponto de vista da

implementacao computacional.

ΠAB(u) = ΠA(u) + ΠB(u) =1

2

uA

uB

t[KA 0

0 KA

]

uA

uB

−

fA

fB

t uA

uB

,

(3.19)

ou de forma compacta,

ΠAB(u) =1

2ut

ABKABuAB − f tABuAB, (3.20)

onde por simplicidade na notacao, os sub-ındices AB podem ser omitidos.

3.2.1 Geometria do contato

Um aspecto importante a ser considerado em problemas de contato e a geometria da

malha. Devido ao fato de que os corpos estao discretizados, e necessario estabelecer-

se uma estrategia na geometria das malhas que permita o calculo da s(x) ∀x ∈ Γc

na funcao ”gap” (g(u) = uAn + uB

n − s ≤ 0). Como apresenta-se na (Fig.3.3) podem

18

existir diferentes graus de refinamentos nas malhas ou elementos o que dificulta esse

calculo.

Figura 3.3: Malhas incompatıveis

Uma primeira estrategia pode ser calcular a distancia entre os nos de uma malha

e os segmentos superficies opostas conforme mostrado na (Figura 3.4). Nota-se que

os campos de deslocamentos sao pontualmente compatıveis em apenas um ponto

por segmento, onde a distancia dp = 0.

Figura 3.4: Situacao de contato no - segmento

Portanto, embora as restricoes g(u) ≤ 0 estejam satisfeitas pontualmente, pode

existir interpenetracao dp. Uma forma de evitar isto e refinando a malha, o que

minimiza os valores de dp, mas acarreta um maior custo computacional. Outra

possibilidade e utilizar a estrategia ”no - no”conforme mostrado na (Figura

3.5).

Neste caso, devido ao posicionamento dos nos, a distancia de penetracao dp

tende a ser nula mas tem com a desvantagem de que e necessario conhecer os nos

associados nas superficies A e B.

Como observa-se na (Figura 3.6), se os nos nao tem um posicionamento correto

19

Figura 3.5: Situacao de contato no - no

Figura 3.6: Campo de deslocamentos incompatıveis na configuracao no - no

apos a deformacao, como resultado de deslocamentos laterais, sera necessario refinar

a malha nesta area.

Tambem existem outras estrategias mais complexas, como a de utilizar elementos

finitos de maior ordem, que contribuem para uma melhor compatibilidade das

malhas em vista de que utiliza os nos intermediarios para definir novas restricoes. Ou

definir superficies de nos ”mestres” e superfıcies de nos ”escravos” com o objetivo

de gerar funcoes analıticas que permitam alcancar maior precisao na solucao. Todas

elas com a desvantagem de um aumento grande do custo computacional. Neste

trabalho foi adotada como estrategia para obtencao de s(x) a configuracao de malhas

compatıveis (no - no). Devido a adocao de um software comercial, a geracao

de malhas compatıveis foi feita utilizando malhas estruturadas e ajuste manual

dos nos para atingir a melhor compatibilidade entre pares de nos. Finalmente,

estas consideracoes podem ser estendidas a geometrias em tres dimensoes de forma

natural.

20

Capıtulo 4

Metodos de resolucao

4.1 Metodo 1 - Problema de otimizacao

Como foi mostrado na (Secao 3.1), o problema de contato sem atrito se resume

a um problema de minimizacao com restricoes de desigualdade. Desta forma a

funcao objetivo a ser minimizada Π(u) representa a energia potencial total dos

corpos e as restricoes de desigualdade correspondem as condicoes cinematicas de

nao penetrabilidade da fronteira na regiao de contato g(u). Como a funcao objetivo

e quadratica, as condicoes de primeira ordem de KKT sao suficientes e pode-se

afirmar que achada a solucao, ela representa um mınimo global e e unica, ver

LUENBERGER [5].

Sendo assim, o problema pode ser represando por,

(P2)

minimizar Π(u)

sujeito a: gj(u) ≤ 0 j = 1, .....,m

(4.1)

onde,

Π(u) =1

2uKut − Fu, (4.2)

g(u) = Au− s, (4.3)

sendo A a matriz (m × n) gradiente das restricoes, onde as filas representam

21

a derivada da restricao j-esima, portanto, tem o valor ”1”na coordenada

correspondente.

Aplicando as condicoes de primeira ordem (KKT) apresentadas no (Anexo B)

obtem-se,

Ku− F + Atλ = 0, (4.4)

G(u)λ = 0, (4.5)

λ ≥ 0, (4.6)

g(u) ≤ 0. (4.7)

onde,

∇f(u) = Ku− F gradiente da funcao objetivo ∈ Rn×1

∇2f(u) = K matriz de rigidez ∈ Rn×n

∇g(u) = A gradiente das restricoes ∈ Rm×n

λ parametros de Lagrange, (forca de contato) ∈ Rm

Gii = gi matriz diagonal das restricoes ∈ Rm×m

Analisando as equacoes, nota-se que a (Eq. 4.4) representa a equacao de equilıbrio

onde os multiplicadores de Lagrange tem o significado fısico de representar as forcas

de contato nos nos. Portanto, o gradiente das restricoes atua como um ponderador,

distribuindo os esforcos de contato entre os corpos na equacao de equilıbrio.

Para resolver este problema de otimizacao, existe uma grande variedade de

tecnicas. A modo de exemplo, pode-se citar o metodo de penalidade. Este,

consiste em transformar o problema de minimizacao com restricoes num problema

sem restricoes atraves da inclusao das restricoes impenetrabilidade ao funcional do

problema. Neste caso, as forcas de contato atuam expulsando a parcela do corpo

que penetrou no outro, portanto, as restricoes nem sempre sao satisfeitas. Para

evitar isto, neste trabalho adotou-se uma tecnica de ponto interior e direcoes viaveis

FDIPA - Feasible Direction Interior Point Algorithm, proposta por HERSKOVITS

[8] [9].

22

4.2 Descricao do algoritmo (FDIPA)

Esta tecnica baseia-se no calculo da direcao de descida da funcao objetivo, a partir de

uma iteracao de Newton nas condicoes de Karush-Kuhn-Tucker (KKT) de primeira

ordem, seguida de uma pequena deflexao para garantir uma direcao viavel no passo.

Para o calculo dessa direcao, resolvem-se dois sistemas lineares nas variaveis

primais e duais x ∈ IRn e λ ∈ IRm, com a matriz chamada primal-dual. Achada a

direcao, efetua-se uma busca linear inexata na funcao objetivo para determinar o

novo ponto. Desta forma, repetindo este processo estabelece-se um procedimento

iterativo que a partir de um valor viavel u0, gera uma sequencia de pontos viaveis

convergentes a solucao do problema, onde a funcao objetivo e reduzida a cada

iteracao.

Portanto, aplicando as condicoes de otimalidade de (KKT) ao problema (Eq.

4.1), obtem-se as (Eqs. 4.4 e 4.5) que determinam o sistema de equacoes primal-

dual. Aplicando o metodo de Newton-Raphson descrito no (Anexo B), obtem-se,

H(xk, λk) ∇g(xk)

Λk∇gt(xk) G(xk)

xk+1 − xk

λk+10 − λk

= −

∇f(xk) + ∇g(xk)λk

G(xk)λk,

, (4.8)

onde (xk, λk) e o ponto de partida das iteracoes, e (xk+1, λk+10 ) o novo ponto estimado.

H(xk, λk) e a matriz Hessiana da funcao Lagrangeana do problema, e Λ a matriz

diagonal com Λii ≡ λi. Alem disto, define-se o sistema, (xk0, λ

k+10 ).

H(xk, λk)dk0 + ∇g(xk)λk+1

0 = −∇f(xk) (4.9)

Λk∇gt(xk)dk0 +G(xk)λk+1

0 = 0, (4.10)

onde dk0 e a direcao de descida no espaco primal definida por dk

0 = xk+10 − xk, e λk+1

0

e a nova estimativa do multiplicador de Lagrange λ.

Pode-se demostrar que dk0 e uma direcao de descida para a funcao f(x), porem, dk

0

nao e uma direcao de busca factıvel, devido ao fato de que se algumas restricoes sao

ativas, forcam dk0 a uma direcao tangente ao conjunto viavel. O seja, se gi(x

k) = 0,

a direcao calculada dk0 sera tangente a restricao Λk∇gt(xk)dk

0 = 0, gerando uma

23

direcao inviavel. Para evitar isto, define-se um novo sistema linear,


0 = −∇f(xk), (4.11)


0 = −ρkλk. (4.12)

Este sistema adiciona um vetor negativo na (Eq. 4.10), onde ρ ∈ IR+. Como

resultado desta perturbacao, se produz uma deflexao em dk0 proporcional a ρk.

Como dk0 e uma direcao de descida, e possıvel encontrar uma estimativa para ρk

que assegure que dk seja uma direcao viavel e de descida.

Figura 4.1: Direcao de busca do algoritmo FDIPA

Como dkt

0 ∇f(xk) ≤ 0, pode-se estimar ρk,

dkt

∇f(xk) ≤ αdkt

0 ∇f(xk), (4.13)

o que implica, dkt

∇f(xk) < 0.

Em geral, a taxa de descida ao longo de dk e menor que ao longo de dk0

constituindo-se no custo a pagar para que a direcao calculada seja factıvel.

24

Considerando agora o sistema auxiliar,


0 = 0 (4.14)


0 = −λk (4.15)

deduz-se que, dk = dk0 + ρdk

1 assim, substituindo esta equacao na (Eq. 4.13) chega-se

a uma expressao para ρ.

ρ ≤(α− 1)dkt

0 ∇f(xk)

dkt∇f(xk)(4.16)

Com a direcao calculada, procede-se a busca linear inexata na funcao objetivo f(x)

com alguma tecnica como Armijo, Wolf ou Goldstein descritas em LUENBERGER

[5]. No caso da funcao objetivo ser quadratica, existem metodos diretos para achar

este mınimo.

Finalmente um novo ponto interior a regiao viavel e calculado, a partir do qual,

atualizando as variaveis, inicializa-se o processo iterativo na procura da solucao, ate

que a condicao de parada seja verificada. Uma aplicacao do FDIPA a problemas de

contato em elasticidade linear pode ser visto no trabalho de AUATT [18].

4.2.1 Aplicacao ao problema de contato

Com o proposito de desacoplar o calculo estrutural do FDIPA, sao efetuadas algumas

adaptacoes no algoritmo. Isto evita a analise da matriz de rigidez K, mas como

consequencia impede o calculo da funcao objetivo no problema de contato (Eq. 4.1).

Utilizando u como a variavel primal, o sistema de equacoes (Eq. 4.8) aplicado ao

problema (Eq. 4.1) e dado por,

K At

ΛA G(uk)

uk+10 − uk

λk+10 − λk

= −

Kuk − F + Atλk

G(uk)λk,

(4.17)

onde define-se dk0 = (uk+1

0 −uk). Substituindo esta expressao no sistema obtem-se

as seguintes equacoes,

Kdk0 + Atλk+1

0 = −(Kuk − F ), (4.18)

ΛAdk0 +G(uk)λk+1

0 = 0, (4.19)

25

chegando a,

λk+10 = (−Λ−1G(uk) + AK−1At)−1(−A(uk −K−1F )) (4.20)

dk0 = −K−1Atλ0 − (uk −K−1F ). (4.21)

De forma analoga, opera-se no sistema auxiliar,

Kdk1 + Atλk+1

1 = 0, (4.22)

ΛAdk1 +G(uk)λk+1

1 = −λk, (4.23)

chegando a,

λk+11 = (−ΛAK−1At +G)−1(−λ), (4.24)

dk1 = −K−1Atλ1. (4.25)

Observa-se que no calculo de d0 e d1 sao necessarias analises de sistemas S2 = K−1F

e S1 = K−1At, que podem ser executados num software comercial, obtendo

dk = dk0 + ρd1, onde ρ e estimado devido a impossibilidade de seu calculo.

Uma vez calculada a direcao dk, estima-se o passo para determinar o novo ponto

da variavel primal. Por conseguinte, o passo t e calculado de forma a garantir a

viabilidade e o mınimo da funcao objetivo na direcao dk. Como a funcao objetivo

e uma funcao quadratica, nao e necessario fazer uma busca inexata e o passo pode

ser determinado pelas as seguintes condicoes,

Primeira condicao

Como a funcao e quadratica o passo que minimiza a funcao pode ser calculado

conforme a seguir,∂

∂tf(u+ td) = 0 (4.26)

onde,

f(u+ td) =1

2(u+ td)K((u+ td)t − F (u+ td)

26

simplificando,

f(u+ td) = f(u) + uKtd+1

2t2dtKd− Ftd

na Eq. 4.26 temos que,

uKd+ tdtKd− Fd = 0 t1 =F td− utKd

dtKd. (4.27)

Como neste trabalho nao e conhecida a matriz de rigidez, esta condicao nao pode

ser aplicada.

Segunda condicao

A segunda condicao impoe que o passo seja factıvel, entao tem que cumprir,

g(u+ tid) ≤ γg(u), (4.28)

onde γ = min(γ, ‖d0‖2), alem disso,

A(u+ t2d) − s ≤ γ(Au− s)

Das restricoes, obtem-se um vetor [m× 1], de onde o passo mınimo e,

t2 = min[ti2], i− 1, ....,m, (4.29)

Em virtude de que o passo deve cumprir ambas condicoes, ele e calculado como

mınimo das duas condicoes definidas nas equacoes (Eq. 4.27) e (Eq. 4.29). Como t1

nao pode ser calculado, a escolha e feita entre 1 e t2, onde,

t = min[1, t2].

Finalmente, com t e d calculados, sao atualizados os parametros deslocamento u e

o coeficiente de Lagrange λ para a proxima iteracao.

u = u+ td,

λ = max(λ0,min[1, ‖d0‖2]).

Devido a comparacao com o problema de complementaridade, o criterio a ser

utilizado e,

‖ ut.F(u) ‖< ǫtol (4.30)

27

onde F(u) e a funcao associada ao problema de complentaridade e ε e a tolerancia

prescrita. Alem disto, monitora-se o resıduo da equacao de equilıbrio atraves das

condicoes,

Condicao 1

‖ u−K−1F −K−1Atλ ‖≤ ε (4.31)

Condicao 2

‖ G(u)λ0 ‖< ε (4.32)

A positividade de λ e garantida pela regra de atualizacao desta variavel ao final de

cada iteracao.

28

4.2.2 Algoritmo

Parametros:. α ∈ (0, 1)

Passo 0: Inicializacao de variaveis e parametros.

k = 0, xk ∈ Ω, λ, ϕ > 0

Passo 1:Direcao de busca

λ0 = (AS1 − Λ−1G)−1A(S2 − u)

d0 = −(S1λ0 + (u− S2))

λ1 = (G− ΛAS1)−1λ

d1 = −(S1λ1)

Calculo de ρ

se dt1(u− S2) > 0

entao ρ = min(ϕ‖d0‖2,

(α−1).dt0

.dt1

)

se nao ρ = min(ϕ‖d0‖2)

fim

λ = λ0 + ρλ1

d = d0 + ρd1

Passo 2: Passo

γ = min(γ, ‖d0‖2)

t = min(1,(1 − γ)λ

λ0

)

29

Passo 3: Atualizacao

u = u+ td

λi = sup(λ0,min(1, ‖d0‖2))

g = Au− s

G = diag(g)

Λ = diag(λ)

Passo 4: Criterio de parada

‖ λk+1.F(λk+1) ‖< ǫtol

Passo 5: Voltar ao passo 1

30

4.3 Metodo 2 - Problema de complementaridade

Como descrito no (Capıtulo 1), no trabalho de BJORKMAN [6] mostra-se que todo

o problema quadratico de minimizacao com restricoes de desigualdade pode ser

reduzido a um problema de complementaridade linear (PCL) nas variaveis duais.

Em virtude disto, o problema de contato pode-se escrever como um problema de

complementaridade atraves do estabelecimento das condicoes de otimalidade (KKT)

e posterior eliminacao das variaveis primais. Para isto, define-se a formulacao geral

de um problema de complementaridade.

Definicao 4.3.1 Seja F : D ⊆ IRn → IRn uma funcao vetorial nao linear, problema

de complementaridade consiste em encontrar x ∈ IRn tal que;

x ≥ 0, F (x) ≥ 0 e x.F (x) = 0, (4.33)

com,

x ≥ 0 ⇔ xi ≥ 0 ∀ 1 ≤ i ≤ n,

F (x) ≥ 0 ⇔ Fi(x) ≥ 0 ∀ 1 ≤ i ≤ n,

xF (x) = 0 ⇔ xiFi(x) = 0 ∀ 1 ≤, i ≤ n

onde, se F e uma funcao ( F (x) = Ax+b , A ∈ Rn×n e b ∈ R

n ), temos um problema

de complementaridade linear (PCL).

Definicao 4.3.2 O conjunto de pontos viaveis Ω do problema e definido como:

Ω := x ∈ IRn | x ≥ 0, F (x) ≥ 0 (4.34)

Portanto, a partir das condicoes de KKT definidas nas (Eqs.4.4 a 4.7) reescrevem-se

as (Eqs. 4.4 e Eq.4.7) em funcao de u,

u(λ) = K−1F −K−1Atλ,

g(u(λ)) = −Au(λ) + s ≥ 0,

31

substituindo, podem-se reescrever como um problema de complementaridade na

variavel dual,

G(λ) = (AK−1At)λ− AK−1F + s ≥ 0,

G(λ)λ = 0, (4.35)

λ ≥ 0,

representando as forcas de contato.

Para resolver problemas desta classe, existem varias tecnicas que basicamente

podem-se dividir em dois grandes grupos.

- No primeiro grupo encontram-se aquelas tecnicas baseadas em metodos de

minimizacao sem restricoes, para o qual sao utilizadas funcoes de merito associadas

ao problema de complementaridade.

Definicao 4.3.3 Uma funcao φ : IRn → IR e denominada de funcao de merito para

problemas de complementaridade se e so se satisfaz: a) φ(x) ≥ 0 para todo x ∈ IRn

b) φ(x) = 0 se e somente se x e uma solucao do problema de complementaridade

Desta forma consegue-se transformar o problema de complementaridade num

problema de minimizacao sem restricoes, que pode ser resolvido por metodos usuais

de programacao matematica.

- No segundo grupo encontram-se aquelas tecnicas baseadas em sistemas de

equacoes, que reformulam o problema de complementaridade como um sistema de

equacoes atraves de NCP-funcoes, ou seja,

Definicao 4.3.4 Uma funcao ψ : IR2 → IR e chamada de NCP-funcao para

problema de complementaridade se satisfaz a condicao:

ψ(a, b) = 0 ⇔ a ≥ 0, b ≥ 0 e a.b = 0. (4.36)

32

Assim, dada uma NCP-funcao ψ, define-se o operador Ψ : IRn → IRn como o sistema

de equacoes associado, da seguinte forma,

Ψ(x) :=

ψ(x1, F1(x))...

ψ(xn, Fn(x))

,Ψ(x) = 0. (4.37)

Portanto, seja ψ uma NCP-funcao e Ψ seu operador de equacoes correspondente, x

e solucao do problema de complementaridade (Eqs. 4.3), se e somente se e solucao

do sistema de equacoes (Eqs. 4.37). Em virtude disto, se F e ψ sao continuamente

diferenciaveis, o operador Ψ tambem sera, e pode-se aplicar o metodo de Newton

para achar o x que faz Ψ = 0

xk+1 = xk − [∇Ψ(xk)]−1Ψ(xk), k = 0, 1, 2 . . . (4.38)

Neste trabalho, sera adotada uma nova tecnica para resolver problemas de

complementaridade proposta em MAZORCHE [16] ou [15].

4.4 Descricao do algoritmo (ALG-NCP)

Como foi visto no capıtulo anterior, um problema de complentaridade consiste em

determinar x ∈ IRn tal que x ≥ 0, F(x) ≥ 0 e xt.F(x) = 0 onde F : IRn → IRn; ou

seja na regiao viavel do problema Ω := x ∈ IRn/x ≥ 0,F(x) ≥ 0

Para resolver este problema e utilizado um metodo baseado no metodo de

sistemas de equacoes, x e solucao do problema de complementaridade (Eq. 4.3),

se e somente se e solucao do sistema de equacoes (Eqs. 4.37).

Baseada nesta propriedade, a tecnica desenvolvida por MAZORCHE [16] se

constitui num novo metodo para problemas de complementaridade. Este metodo

utiliza ideias de FDIPA, gerando uma sequencia de pontos interiores a Ω, que

convergem a uma solucao do problema de complementaridade mediante a solucao

de dois sistemas lineares e uma busca linear na funcao potencial φ = xt.F(x)

respeitando as restricoes x ≥ 0 e F(x) ≥ 0.

33

Figura 4.2: Direcao de busca do algoritmo ALG-NCP

A proposta e resolver o seguinte sistema de equacoes dentro da regiao Ω ⊂ IRn

utilizando uma iteracao de Newton para construir a sequencia de pontos viaveis que

convergira a solucao do problema.

H(x) ≡ x • F(x) ≡

x1.F1(x)...

xn.Fn(x)

=

0...

0

(4.39)

A direcao de busca d1 pode ser calculada a partir do seguinte sistema nao linear,

∇H(x)d1 = −x • F(x) (4.40)

onde,

∇H(x) =

F1(x). . .

Fn(x)

+

xk+11

. . .

xk+1n

∇F1(x)...

∇Fn(x)

, (4.41)

com,

∇Fi(x) = (∂Fi(x)

∂x1

,∂Fi(x)

∂x2

. . .∂Fi(x)

∂xn

). (4.42)

Como dk1, em geral, nao e uma direcao viavel, de forma semelhante ao FDIPA,

calcula-se uma nova direcao utilizando a mesma matriz (Eq. 4.41) com a finalidade

34

de gerar uma restauracao, que torne a direcao viavel. Para isto define-se o novo

sistema como,

∇H(x)d2 = µE (4.43)

onde,

E = [1, 1, ...., 1]t ∈ IRn e µ =xt.F(x)

n

como dk2 nao e uma direcao de descida, a nova direcao d e calculada como uma

combinacao linear entre a direcao de Newton dk1, e a direcao dk

2 viavel para todo ponto

de Ω. Esta combinacao linear possui as duas propriedades o que e demonstrado em

MAZORCHE [16].

Assim, atraves da resolucao destes dois sistemas lineares, calcula-se a direcao

de busca d como uma combinacao da direcao de Newton dk1 e de um desvio dk

2 que

restaura a viabilidade do passo, portanto e interior a Ω e de descida para a funcao

potencial.

d = d1 + ρd2 com ρ ∈ (0, 1) (4.44)

Apos o calculo da direcao, o novo ponto e calculado fazendo uma busca linear

inexata na funcao potencial e nas restricoes da funcao potencial. Para isto, adota-se

a tecnica de Armijo LUENBERGER [5]. Nesta busca inexata, calcula-se o passo t

como o primeiro termo da sequencia,

1, ν, ν2, ν3, .....,

que satisfaz a seguinte inequacao,

φ(xk) + tν1∇φ(xk)tdk ≥ φ(xk + tdk) (4.45)

xk1 + tdk ≥ 0 (4.46)

F(xk + tdk) ≥ 0 (4.47)

onde η1 ∈ (0, 1) e ν1 ∈ (0, 1).

35

Portanto, logo procedese ao calculo do novo ponto x,

xk+1 = xk + dxk.

Finalmente, efetuadas as atualizacoes, repete-se o procedimento verificando o

seguinte criterio de parada,

‖ xt.F(x) ‖< ǫtol, (4.48)

gerando assim uma sequencia de pontos em Ω que e decrescente para a funcao

potencial, e converge a solucao do problema de complementaridade.

4.4.1 Aplicacao ao problema de contato

A partir do problema de complentaridade da forma (Eqs. 4.3) tem-se que a funcao

que define o problema de contato e sua derivada e escrito por,

F(λ) = A(K−1Atλ−K−1F ) + s (4.49)

∇F(λ) = AK−1At (4.50)

Como os termos S2 = K−1At e S1 = K−1F permanecem constantes, sao

calculados antes da analise do algoritmo, permitindo assim desacoplar a analise

estrutural do problema de complementaridade.

Definindo os seguintes sistemas:

S1 = K−1At (4.51)

S2 = K−1F (4.52)

Em virtude destes sistemas compartilharem a matrizK, eles podem ser resolvidos

de forma simultanea e eficiente por um software comercial. Assim, uma vez obtida a

funcao e a sua derivada, procede-se a resolucao do problema de complementaridade

com menor custo computacional para o algoritmo (ALG-NCP) .

Em caso de haver contato entre dois corpos A e B, a partir da equacao

generalizada, (Eq. 3.2) chega-se a estrutura das matrizes S1 e S2 como,

36

S1 =

SA1

SB1

; S2 =

SA2

SB2

; A =

AA 0

0 AB

(4.53)

de onde resulta que e necessario uma analise estrutural por corpo, por exemplo uma

para calcular SA1 , SA

2 , e outra para SB1 , SB

2 .

Por conseguinte, pode-se generalizar esta estrutura matricial em caso de contato

entre mais corpos, realizando desta forma analises simultaneas e distribuıdas, se for

necessario com o objetivo de acelerar o processo da geracao das matrices S1 e S2.

37

4.4.2 Algoritmo

Parametros: α ∈ (0, 1), α1 ∈ (0, 1), e ν1 ∈ (0, 1)

Passo 0: Inicializacao de variaveis.

k = 0, λk ∈ Ω, E ∈ IRn, F (λk), H(λk), φ(λk), µk = φ(λk)n

e ∇F (λk)

Passo 1: Direcao de busca

M = ∇H(λ) (4.54)

Mdk1 = −H(λk) (4.55)

Mdk2 = µk.E (4.56)

Calculo de ρ

se ∇φ(λk)dk2 < 0

entao ρ = minα1; ‖d1‖, (α− 1)∇dt

1φ(λk)

∇dt2φ(λk)

se nao ρ = minα1; ‖d1‖

fim

dk = dk1 + ρdk

2

Passo 2 : Busca linear

Define-se o tamanho do passo t como o primeiro valor

da sequencia 1, ν, ν1, ν2, ν3, ν4, ... que satisfaz:

λk1 + tdk ≥ 0 (4.57)

F (λk + tdk) ≥ 0 (4.58)

φ(λk) + tν1∇φ(λk)tdk ≥ φ(λk + tdk) (4.59)

38


λk+1 = λk + tdk

k = k + 1 (4.60)

H = H(λk) (4.61)

M = ∇H(λk) (4.62)

µ =f(λk)

n(4.63)

Passo 4: Criterio de Parada

‖ λk+1 · F(λk+1) ‖< ǫtol

Passo 5: Voltar passo 1

39

Capıtulo 5

Interface com ABAQUS

Como foi visto no capıtulo anterior, o desacoplamento da analise estrutural

apresenta-se como uma estrategia valida para resolver problemas de contato

utilizando os algoritmos descritos. Isto, em razao de que analises do tipo (K−1f = u)

sao varias vezes repetidas para diferentes vetores f com a mesma matriz de rigidezK.

Visando explorar as vantagens deste desacoplamento, foi desenvolvida uma interface

entre os algoritmos e um codigo comercial de elementos finitos.

A estrategia permite analisar problemas de grande porte, ja que a maior parte do

processamento e executada numa plataforma comercial atuando como ferramenta

auxiliar dos algoritmos. Alem da eficiencia na analise, esta estrategia apresenta

outras vantagens importantes, como por exemplo, facilidades no pre-processamento

e pos-processamento. Tambem possibilita realizar analises de forma remota em

plataformas computacionais diferenciadas segundo a necessidade, ou de forma

simultanea no caso de trabalhar com muitos corpos em contato. Para a construcao

da interface, optou-se por utilizar o ABAQUS, mas a presente estrategia e extensiva

a outros programas comerciais.

40

5.1 Descricao do programa

O ABAQUS e um conjunto de programas destinado a modelagem de sistemas

estruturais atraves do metodo de elementos finitos. Nele, podem-se realizar analises

numa grande variedade de problemas do tipo linear, nao linear, estatico, dinamico,

estrutural e termico. Este programa dispoe de ferramentas especıficas para resolucao

de grandes sistemas com estrutura esparsa, assim como tecnicas evoluıdas de

armazenamento que se traduzem numa eficiente e rapida resolucao destes problemas.

Os principais programas que integram ABAQUS sao:

- ABAQUS/Standard e utilizado na analise linear estatica, dinamica e termica,

com controle do passo de tempo.

- ABAQUS/Explicit e utilizado na analise nao linear, transitoria, assim como em

aplicacoes quasi-estaticas que envolvem comportamento nao linear descontınuo.

- ABAQUS/CAE e o ambiente grafico onde se integram todos os programas

ABAQUS. Apresentando uma interface simples para criar, submeter, monitorar e

avaliar os resultados das simulacoes. Sua estrutura encontra-se dividida em modulos

que correspondem a sucessao logica da criacao de um modelo. Evoluindo neles,

define-se a geometria, as propriedades dos materiais, o tipo de analise e a geometria

da malha. O resultado deste processo gera um arquivo de extensao (*.inp) onde

e armazenada esta informacao do pre-processador. O ultimo modulo e onde se

submete a analise e monitora sua evolucao. Como resultado, gera-se uma base de

dados (*.odb) com os resultados.

- ABAQUS/Viewer e ambiente de pos-processamento onde se faz a visualizacao

dos resultados do banco de dados (*.odb). Sua capacidade inclui o tracado de

graficos, animacao de resultados e geracao de informes da evolucao das variaveis.

- ABAQUS/Design e o modulo que faz analise de sensibilidade (DSA).

41

5.1.1 Acessibilidade e interacao

Neste ponto destacam-se tres caracterısticas do programa:

- Acesso aos resultados da base de dados (*.odb) mediante scripts desenvolvidos

em Python ou C++, assim como disponibilidade de outras fontes acessaveis para

pos-processamento, como se mostra na (Figura 5.1) cujas caracterısticas encontram-

se detalhadas na referencia [1].

- Capacidade de importar modelos vetoriais em diversos formatos, criados em

software CAD como Parasolid, SolidWorks, Pro-Engineer, Autocad, assim como

malhas orfas criadas em outros programas de (MEF ).

- Possibilidade de interagir facilmente com outros pacotes comerciais de analise

por elementos finitos como MSC/PATRAN, MSC/NASTRAN ou CATIA atraves

de ferramentas especıficas de traducao dos arquivos de entradas do pre-processador

e saıdas pos-processador.

5.2 Descricao da interface

A interface e composta de tres ferramentas desenvolvidas em Fortran 90 e

Python, integradas aos algoritmos em Matlab de onde sao executados de forma

automatica. A comunicacao efetua-se atraves de arquivos texto (*.txt) ASCII, onde a

compatibilidade com algoritmos desenvolvidos em outras linguagens de programacao

e garantida.

A seguir, faz-se uma descricao das funcoes da cada ferramenta da interface:

• ReadINP.exe - sua funcao e ler o arquivo do pre-processador (*.inp) e

escrever num arquivo texto (Pre.txt) acessıvel pelo algoritmo Matlab/ALG-

PC.

• WriteINP.exe - sua funcao e escrever no arquivo do pre-processador

(*.inp) as condicoes de carregamento, escritas no arquivo (F.txt). Quando

42

armazenadas varias condicoes de carregamento, elas sao analisadas si-

multaneamente.

• ReadODB.py - sua funcao e ler da base de dados (*.odb) os resultados da

analise e escrever-los num arquivo de texto (U.txt) acessıvel pelo algoritmo em

forma de vetores sequenciais, segundo o numero de analises executadas.

O codigo destas ferramentas esta disponıvel no (Anexo C).

Figura 5.1: Esquema da Interface

O codigo destas ferramentas esta disponıvel no (Anexo C).

5.3 Funcionamento da interface

Em linhas gerais, um processo de analise efetua os seguintes passos:

- Primeiro passo, se constroi a geometria do modelo em um formato vetorial, num

software CAD ou no ABAQUS/CAE. Apos disto, executam-se os passos do pre-

processo, definindo materiais, condicoes de contorno, geometria da malha e outros

detalhes. Dados que sao armazenados num arquivo (*.inp) do pre-processador.

- Segundo passo, identificam-se os nos potenciais de contato no algoritmo e

executa-se a analise. Atraves das ferramentas descritas, realiza-se a interacao

43

automatica com ABAQUS. Como resultado da pre-analise, obtem-se o arquivo de

deslocamentos U.txt a partir da qual geram-se as matrizes S1 e S2. Apos disto,

os algoritmos estao em condicoes de realizar a busca da solucao, obtendo como

resultado o vetor forcas nodais de contato.

- Terceiro passo, com as forcas nodais obtidas, e realizada uma nova analise,

chamada de pos-analise com o objetivo de pos-processar os resultados obtidos no

algoritmo. Uma vez finalizada esta etapa, os resultados tornam-se acessaveis no

modulo ABAQUS/Viewer, o que possibilita o acesso a uma grande variedade de

resultados disponıveis na base de dados (*.odb).

Proseguindo, se apresenta-se o diagrama da interacao entre o Matlab e o

ABAQUS atraves da interface.

Figura 5.2: Diagrama da interacao Matlab - ABAQUS

44

Capıtulo 6

Resultados numericos

Neste capıtulo, apresentam-se alguns exemplos que mostram o comportamento

numerico dos algoritmos estudados, efetuando-se a comparacao com os resultados

obtidos em ABAQUS.

Com esta finalidade, apresenta-se o tempo de parede gasto pelo ABAQUS, e

pelos algoritmos junto a interface. Sendo os tempos da interface, os que envolvem

o tempo de pre-analise, analise do algoritmo e pos-analise. Apresenta-se tambem o

numero de iteracoes feitas em cada analise, e o grafico da evolucao do criterio de

parada para cada exemplo.

A precisao dos resultados e calculada com base na norma infinito da diferenca

das pressoes de contato e dos deslocamentos obtidas do pos-processamento.

Utiliza-se a precisao simples (1 × 10−6) do ABAQUS, e tolerancia (1 × 10−8)

no criterio de parada dos algoritmos. As malhas de elementos finitos utilizadas

detalham-se em cada exemplo e correspondem a biblioteca do ABAQUS. Os testes

sao executados num computador AMD Athlon 1800Mhz com 512Mb de memoria

RAM e sistema operacional Windows XP.

45

6.1 Exemplo 1 - Semicilindro (2D)

Este exemplo e utilizado com a finalidade de verificar a interface, pois possui solucao

analıtica. Assim, tem-se um corpo semicilındrico de material homogeneo e isotropico

em contato com uma fundacao rıgida e submetido a um estado plano de tensoes.

Devido a simetria do exemplo, analisa-se a 1/4 do cilindro, impondo restricoes de

deslocamento na vertical como se ilustra na (Figura 6.1). A geometria da malha e

composta de elementos triangulares lineares (CPS3).

RP

1

2

3

Figura 6.1: Modelo 1 condicoes de contorno

Apresentam-se as curvas de pressao analıticas e numericas obtidas com as

diferentes estrategias. Como nos potenciais de contato sao considerados todos os

nos da face curva, totalizando 60 nos.

Dados do modelo.

—————————————————————————

R = 1m E = 207GPa ν = 0.3 q = 5.5GN/m

—————————————————————————

A seguir apresentam-se as equacoes da solucao analıtica das pressoes de contato:

p(x) =2F

πb

√

1 − (x

b)2 b =

√

4FR(1 − ν2)

πE

pmax =2F

πbLB =

E

2(1 − ν2)F = 2Rq

onde b e a largura da regiao de contato, q e a forca por unidade de longitude e R e

o raio do cilindro.

46

Geometria.

Elementos tipo = CPS3 No. de elementos=1062 No. de nos=587

(Ave. Crit.: 75%)S, Mises

+4.958e+00+6.641e+00+8.323e+00+1.001e+01+1.169e+01+1.337e+01+1.505e+01+1.673e+01+1.842e+01+2.010e+01+2.178e+01+2.346e+01+2.515e+01

Step: Step-1Increment 8: Step Time = 1.000Primary Var: S, MisesDeformed Var: U Deformation Scale Factor: +1.000e+00

1

2

3

Figura 6.2: Configuracao deformada com tensoes de Von-Mises

1

2

3

Figura 6.3: Malha deformada

0 0.05 0.1 0.15 0.2 0.25 0.3 0.350

0.05

0.1

0.15

0.2

0.25

0.3

0.35Solução obtida com elementos Tri

x/r

P(x

)/B

AnalíticaABAQUSAlg−PCFDIPA

Figura 6.4: Curva pressao na regiao de contato

47

Tabela 6.1: Tempos das analisesABAQUS Analise 9.00s

FDIPA Pre analise 22.00s ALG-NCP Pre analise 22.00s

Algoritmo 1.70s Algoritmo 0.92s

Pos analise 1.00s Pos analise 1.00s

Tabela 6.2: Diferenca entre solucoes‖PABAQUS − PFDIPA ‖∞=3.5890e-3 ‖uABAQUS − uFDIPA ‖∞=6.1962e-8

‖PABAQUS − PALG−NCP ‖∞=2.6700e-4 ‖uABAQUS − uALG−NCP ‖∞=1.5000e-8

‖PFDIPA − PALG−NCP ‖∞=3.4580e-3 ‖uFDIPA − uALG−NCP ‖∞=6.2981e-8

0 5 10 15 20

0

20

40

60

80

100

120

140

160

180

Nro de iterações

xn*F

n

ALG−NCPFDIPA

Figura 6.5: Criterio de parada

48

Tabela 6.3: Numero de iteracoes com FDIPAiter t norm(d) xn*Fn cond1 cond2

0 1.000000 0.000e+000 4.432e+002 1.116e+001 0.000e+000

1 0.686853 4.208e+001 1.246e+002 3.232e+000 7.342e-001

2 0.563609 2.387e+001 7.505e+001 2.576e+000 5.885e-001

3 0.629837 1.860e+001 4.602e+001 3.320e+000 3.872e-001

4 0.750049 1.375e+001 2.738e+001 3.564e+000 1.890e-001

5 0.916895 9.382e+000 1.503e+001 2.965e+000 4.111e-002

6 1.000000 5.775e+000 7.213e+000 1.904e+000 5.551e-016

7 1.000000 3.533e+000 3.227e+000 1.054e+000 1.146e-003

8 1.000000 2.203e+000 1.412e+000 5.294e-001 1.428e-003

9 1.000000 1.293e+000 5.753e-001 2.379e-001 5.533e-004

10 0.820487 7.414e-001 1.683e-001 9.105e-002 5.612e-003

11 0.680114 3.371e-001 4.326e-002 3.332e-002 6.300e-003

12 0.772360 1.571e-001 1.497e-002 1.092e-002 3.005e-003

13 0.835462 6.825e-002 4.356e-003 3.290e-003 1.165e-003

14 0.871233 2.507e-002 1.154e-003 9.445e-004 4.455e-004

15 0.881758 8.668e-003 3.005e-004 2.675e-004 2.005e-004

16 0.874682 3.087e-003 7.643e-005 7.622e-005 1.052e-004

17 0.876931 1.135e-003 1.906e-005 2.145e-005 5.145e-005

18 0.861946 4.076e-004 4.552e-006 6.024e-006 2.752e-005

19 0.876957 1.505e-004 1.023e-006 1.521e-006 1.125e-005

20 0.940706 5.255e-005 1.613e-007 2.310e-007 2.070e-006

Tabela 6.4: Numero de iteracoes com ALG-NCPIter IterBL t xn*Fn ro norm(d)

0 0 0.000000 4.432e+002 0.000e+000 0.000e+000

1 0 1.000000 3.963e+001 6.396e-001 8.940e-002

2 0 1.000000 2.139e+000 7.639e-002 5.399e-002

3 1 0.800000 3.818e-001 4.510e-003 1.784e-001

4 0 1.000000 5.825e-002 6.413e-001 1.526e-001

5 0 1.000000 1.033e-002 4.363e+000 1.774e-001

6 0 1.000000 1.941e-003 8.963e-001 1.878e-001

7 0 1.000000 1.129e-004 6.034e-002 5.819e-002

8 0 1.000000 1.629e-007 1.453e-003 1.442e-003

9 0 1.000000 6.274e-015 3.990e-008 3.852e-008

49

6.2 Exemplo 2 - Mola (3D)

Neste caso, apresenta-se um exemplo 3D, de um solido com geometria similar a

uma mola. A peca encontra-se engastada num extremo e sujeito a forcas externas

no outro, como e mostrado na (Figura 6.2). O contato e estabelecido com uma

fundacao rıgida representada na funcao analıtica y=0. A discretizacao e feita com

uma geometria de malha composta por 1130 elementos finitos hexaedricos lineares

(C3D8R) e conta com 2052 nos. Como nos potenciais de contato sao considerados

100 nos localizados na superfıcie inferior do solido.

Dados do modelo.

————————————————————————————————–

rint=4.0cm rext=5.2cm h=0.5cm E=207GPa ν=0.3 F=-12KN

————————————————————————————————–

RPX

Y

Z

1

2

3


50


+5.723e-02+1.655e-01+2.737e-01+3.819e-01+4.901e-01+5.983e-01+7.065e-01+8.148e-01+9.230e-01+1.031e+00+1.139e+00+1.248e+00+1.356e+00

1

2

3

Figura 6.7: Configuracao deformada com tensoes de Von-Mises

−6 −4 −2 0 2 4 6−6

−4

−2

0

2

4

6−1

−0.5

0

0.5

1

Nós sem contato

Nós com contato

Nós engastados

Figura 6.8: Nos da face do modelo em contato na superfıcie analıtica

51








0 5 10 15 20 25

0

50

100

150

200

Nro de iterações

xn*F

n

ALG−NCPFDIPA


52


0 1.000000 3.837e-014 2.343e+002 2.598e+000 7.105e-015

1 0.977216 5.223e+001 2.308e+002 2.535e+000 4.905e-002

2 1.000000 2.466e+003 9.704e+001 1.024e+000 2.558e-013

3 1.000000 1.586e+003 4.059e+001 4.122e-001 1.705e-013

4 1.000000 1.006e+003 1.680e+001 1.824e-001 1.421e-013

5 1.000000 6.183e+002 6.784e+000 7.871e-002 9.948e-014

6 1.000000 3.567e+002 2.620e+000 3.240e-002 7.105e-014

7 0.964806 1.849e+002 9.524e-001 1.301e-002 2.570e-001

8 0.933910 8.893e+001 3.511e-001 5.663e-003 2.288e-001

9 0.946743 4.006e+001 1.350e-001 2.714e-003 8.197e-002

10 0.859036 1.644e+001 5.089e-002 1.289e-003 8.172e-002

11 0.781114 9.281e+000 2.491e-002 6.089e-004 6.944e-002

12 0.780828 6.811e+000 1.550e-002 3.052e-004 5.443e-002

13 0.816769 4.486e+000 9.623e-003 1.726e-004 3.028e-002

14 0.748119 2.465e+000 5.388e-003 1.282e-004 2.243e-002

15 0.682439 1.450e+000 3.166e-003 9.664e-005 1.641e-002

16 0.697920 9.917e-001 2.076e-003 8.036e-005 1.082e-002

17 0.733246 7.433e-001 1.445e-003 7.307e-005 7.341e-003

18 0.786338 5.569e-001 8.656e-004 5.361e-005 4.401e-003

19 0.775128 3.807e-001 2.435e-003 2.858e-005 1.205e-002

20 0.769295 2.394e-001 6.521e-003 1.274e-005 3.684e-002

21 0.773411 1.318e-001 1.533e-002 4.675e-006 8.468e-002

22 0.766664 6.059e-002 5.681e-003 1.454e-006 3.115e-002

23 0.813122 2.594e-002 1.378e-003 2.973e-007 7.548e-003

24 0.944746 1.036e-002 2.771e-004 1.733e-008 1.519e-003

25 0.972795 3.145e-003 2.965e-005 3.931e-009 1.625e-004

26 0.991104 9.581e-004 3.093e-006 6.398e-010 1.695e-005

27 0.998529 1.430e-004 6.927e-008 3.408e-011 3.796e-007


0 0 0.000000 2.343e+002 0.000e+000 0.000e+000

1 0 1.000000 2.097e+001 5.843e-001 8.951e-002

2 0 1.000000 4.256e-001 1.852e-002 2.030e-002

3 1 0.800000 8.361e-002 2.157e-004 1.965e-001

4 1 0.800000 1.747e-002 1.607e-002 2.090e-001

5 0 1.000000 2.652e-003 1.588e-001 1.518e-001

6 0 1.000000 4.865e-004 2.094e-001 1.835e-001

7 0 1.000000 4.436e-005 1.054e-001 9.117e-002

8 0 1.000000 4.553e-007 9.901e-003 1.026e-002

9 1 0.800000 8.990e-008 1.002e-004 1.974e-001

10 0 1.000000 1.781e-008 3.059e-001 1.981e-001

11 0 1.000000 2.640e-010 1.483e-002 1.482e-002

53

6.3 Exemplo 3 - Pinca (3D)

Neste exemplo, e resolvido um problema de auto-contato de um solido, adotando-se

como modelo uma geometria em 3D similar a uma pinca. A peca encontra-se su-

jeita a esforcos externos e engastada na sua base como observa-se na (Figura 6.3).

A peca foi discretizada com uma malha de geometria que tem 4630 elementos finitos

hexaedricos lineares do tipo (C3D8R) em 6226 nos. Como restricoes sao considera-

dos 44 nos potenciais de contato em cada superfıcie de contato.

Dados do modelo.

———————————————————

E = 207GPa ν = 0.3 P = 0.01GPa

———————————————————

X

Y

Z

1

2

3


54


+4.294e-05+3.296e-03+6.549e-03+9.803e-03+1.306e-02+1.631e-02+1.956e-02+2.282e-02+2.607e-02+2.932e-02+3.258e-02+3.583e-02+3.908e-02

1

2

3

Figura 6.11: Geometria da malha e tensoes de Von-Mises


+4.294e-05+3.296e-03+6.549e-03+9.803e-03+1.306e-02+1.631e-02+1.956e-02+2.282e-02+2.607e-02+2.932e-02+3.258e-02+3.583e-02+3.908e-02

1

2

3

Figura 6.12: Ampliacao da zona de contato

55








0 2 4 6 8 10 12 14

0

0.1

0.2

0.3

0.4

0.5

0.6

Nro de iterações

xn*F

n

ALG−NCPFDIPA


56


0 1.000000 1.110e-016 6.638e-001 3.033e-002 1.110e-016

1 0.998048 1.375e+000 6.514e-001 2.978e-002 6.853e-005

2 1.000000 5.093e+001 2.750e-001 1.254e-002 1.332e-015

3 1.000000 3.303e+001 1.160e-001 5.279e-003 1.110e-015

4 1.000000 2.135e+001 4.879e-002 2.218e-003 6.661e-016

5 1.000000 1.371e+001 2.040e-002 9.289e-004 4.441e-016

6 1.000000 8.655e+000 8.419e-003 3.838e-004 4.441e-016

7 1.000000 5.267e+000 3.377e-003 1.542e-004 3.331e-016

8 1.000000 2.970e+000 1.283e-003 5.902e-005 2.776e-016

9 1.000000 1.462e+000 4.498e-004 2.131e-005 2.776e-016

10 0.914616 5.929e-001 1.438e-004 8.453e-006 1.292e-003

11 0.794872 2.717e-001 2.448e-005 3.079e-006 1.420e-003

12 0.730647 6.435e-002 2.070e-006 6.702e-007 4.363e-004

13 0.864192 8.644e-003 3.120e-007 1.219e-007 3.463e-005

14 0.965766 1.155e-003 2.225e-008 1.360e-008 2.108e-006


0 0 0.000000 6.638e-001 0.000e+000 0.000e+000

1 0 1.000000 3.386e-002 5.534e-002 5.100e-002

2 0 1.000000 8.907e-004 1.133e-003 2.631e-002

3 1 0.800000 1.780e-004 1.044e-003 1.999e-001

4 0 1.000000 3.452e-005 2.432e-001 1.939e-001

5 0 1.000000 6.844e-006 2.922e-001 1.983e-001

6 0 1.000000 1.986e-007 2.932e-002 2.901e-002

7 0 1.000000 5.629e-012 2.837e-005 2.835e-005

57

6.4 Exemplo 4 - Pinhao e cremalheira (2D)

Neste exemplo, e resolvido um caso de contato entre dois solidos. Para isto, adota-

se como modelo um mecanismo de pinhao e cremalheira criado numa ferramenta

CAD. No modelo, o pinhao esta fixo e a cremalheira desloca num sentido, como

resultado da aplicacao de uma forca no extremo, como e observado se na (Figura

6.4). O pinhao e discretizado atraves de uma malha de 987 nos e 1221 elementos

finitos triangulares lineares do tipo (CPS3); a cremalheira foi discretizada atraves de

uma malha de 1276 nos e 2208 elementos finitos triangulares lineares do tipo (CPS3)

Dados do modelo.

———————————————————

E = 207GPa ν = 0.3 F = 130KN

———————————————————

1

2 3


58

1

2 3

Figura 6.15: Geometria da malha


+4.417e-08+1.333e-01+2.666e-01+3.999e-01+5.332e-01+6.665e-01+7.998e-01+9.331e-01+1.066e+00+1.200e+00+1.333e+00+1.466e+00+1.600e+00

1

2 3

Figura 6.16: Tensoes de Von-Mises nos elementos

59





Tabela 6.14: Diferenca entre solucoes‖PABAQUS − PFDIPA ‖∞=0.3470524 ‖uABAQUS − uFDIPA ‖∞=0.0143314

‖PABAQUS − PALG−NCP ‖∞=0.3470565 ‖uABAQUS − uALG−NCP ‖∞=0.0143314


0 2 4 6 8 10 12 140

5

10

15

20

25

30

35Evolução do criterio de parada

Nro de iterações

xn*F

n

ALG−NCPFDIPA


60


0 1.000000 1.442e-018 3.489e+001 2.707e+000 4.337e-019

1 1.000000 1.865e+001 1.474e+001 1.139e+000 6.661e-016

2 1.000000 1.211e+001 6.219e+000 4.792e-001 3.331e-016

3 1.000000 7.845e+000 2.620e+000 2.029e-001 4.441e-016

4 1.000000 5.061e+000 1.100e+000 8.588e-002 2.220e-016

5 1.000000 3.230e+000 4.572e-001 3.614e-002 1.665e-016

6 1.000000 2.009e+000 1.858e-001 1.486e-002 1.388e-016

7 1.000000 1.281e+000 6.336e-002 5.122e-003 8.327e-017

8 0.896826 6.776e-001 1.305e-002 1.655e-003 2.805e-003

9 0.650938 2.073e-001 1.344e-003 5.779e-004 3.037e-003

10 0.697745 6.848e-002 6.403e-004 1.908e-004 9.638e-004

11 0.886519 3.002e-002 1.234e-004 5.113e-005 1.603e-004

12 0.908248 6.302e-003 1.847e-005 9.063e-006 3.005e-005

13 0.983455 1.324e-003 1.072e-006 6.363e-007 1.335e-006


0 0 0.000000 3.489e+001 0.000e+000 0.000e+000

1 0 1.000000 3.131e+000 5.930e-001 8.972e-002

2 0 1.000000 1.433e-001 4.917e-002 4.578e-002

3 0 1.000000 1.751e-002 3.931e-001 1.222e-001

4 0 1.000000 3.178e-003 1.483e+000 1.815e-001

5 0 1.000000 5.887e-004 2.017e+000 1.853e-001

6 0 1.000000 6.336e-005 1.126e-001 1.076e-001

7 0 1.000000 9.894e-008 1.578e-003 1.561e-003

8 0 1.000000 1.542e-016 1.157e-009 1.559e-009

61

Capıtulo 7

Consideracoes finais

7.1 Conclusoes

Neste trabalho, foi abordado o problema de contato entre corpos elasticos, na

hipotese de pequenos deslocamentos sem atrito. Para a resolucao, foram adotadas

duas formulacoes. A saber, como um problema de minimizacao com restricoes de

desigualdade, e como um problema de complementaridade.

Como ambas formulacoes permitiram desacoplar o calculo estrutural, foi

desenvolvida uma interface com ABAQUS para atuar como ferramenta de analise

auxiliar dos algoritmos.

A utilizacao desta estrategia de analise mostrou-se conveniente pois:

- prescinde de conhecer a matriz de rigidez do problema;

- facilita as tarefas de pre-processamento e pos-processamento;

- permite trabalhar com recursos computacionais diferenciados segundo a

necessidade.

Assim, o uso de um software comercial de analise estrutural interagindo com

algoritmos apresenta-se como uma alternativa valida de analise de problemas de

grande porte, auxiliando na reducao de um sistema fısico a um modelo numerico

que o resolve de forma eficiente.

Nos exemplos apresentados, foi avaliado o desempenho da interface, assim como

dos algoritmos interagindo com ABAQUS. Com esse objetivo, foram realizadas

comparacoes entre tempos de processamento, diferenca entre solucoes do pos-

62

processamento e velocidade de convergencia.

Com relacao aos tempos gastos no processamento dos algoritmos utilizando a

interface, estes foram razoaveis, se comparados ao gasto pelo ABAQUS. Levando

em conta que ainda pode-se melhorar muito a comunicacao, tanto otimizando o

codigo da interface, quanto implementando os algoritmos em linguagem de alto

nıvel, pode-se concluir que a estrategia aplicada e valida para a analise deste tipo

de problemas.

Dos resultados numericos obtidos com os algoritmos, estes apresentaram-se

satisfatorios quando comparados a resultados obtidos com o ABAQUS. Somente

nos exemplos onde existe contato entre superfıcies deformaveis, a diferenca entre

as solucoes dos algoritmos e do ABAQUS foi um pouco maior. O que evidencia a

necessidade de adotar uma melhor estrategia para caracterizar a distancia entre os

corpos a ser considerada no calculo da restricao g(u).

No caso do algoritmo de complementaridade (ALG-NCP) o desempenho foi

superior ao algoritmo de minimizacao (FDIPA). Isto deve-se as dificuldades para

o calculo do passo t visto no capıtulo 4, o que pode ter acarretado o atraso na

velocidade de convergencia do algoritmo. Neste ponto, a formulacao como problema

de minimizacao apresentou mais dificuldade para o desacoplamento do calculo

estrutural do que a equivalente de complentaridade, onde tanto a funcao como sua

derivada sao calculadas sem ter a necessidade de acesso a matriz de rigidez.

Da mesma forma, os resultados revelam a robustez e eficiencia dos algoritmos

aplicados, atingindo em poucas iteracoes as solucoes, e permitem vislumbrar uma

otima perspectiva quanto a resolucao de problemas de grande porte, aproveitando

as vantagens do codigo comercial atuando como ferramenta auxiliar de novos

algoritmos.

7.2 Sugestoes para trabalhos futuros

Varios aspectos podem ser considerados como possıveis linhas de estudo para a

continuidade deste trabalho.

- No que diz respeito a geometria do contato entre corpos deformaveis, e

63

necessario melhorar o modelo utilizado para caracterizar a distancia fısica entre

os corpos, visando trabalhar em grandes deformacoes.

- Ambos algoritmos apresentaram um otimo desempenho, em particular o de

complementaridade mostrou-se mais robusto. Em virtude disto, a extensao as

aplicacoes incluıdo-se outros tipos de nao linearidades presentes no fenomeno de

contato como grandes deformacoes, plasticidade e problemas dinamicos, assim como

atrito, onde introducao de um modelo satisfatorio e uma sugestao para futuros

trabalhos explorando sua formulacao como problema de complementaridade.

- Aprofundar o conhecimento do ABAQUS e outros pacotes comerciais de analise

estrutural, explorando a capacidade destas ferramentas na analise de sensibilidade

(DSA), com o objetivo de atingir analises de topologia em problemas de contato e

outros problemas de otimizacao estrutural.

Finalmente, e importante ressaltar que o emprego de novos metodos numericos

interagindo com software comercial de analise estrutural deve ser mais explorado,

pois desta forma facilita-se a rapida disseminacao nos ambientes industriais dos

avancos nesta area e constitui-se num banco de provas para as novas ferramentas

desenvolvidas no ambiente academico.

64

Referencias Bibliograficas

[1] ABAQUS 6.4. Analisys User’s Guide. c©ABAQUS, Inc.., 2003.

[2] Klarbring A. A mathematical programming approach to three-dimensional

contact problems with friction. Comp. Methods Appl. Mech. Eng., 58:175–200,

1986.

[3] Mijar A.R. and Arora J.S. Review of formulations for elastostaic frictional

contact problems. Structural Multidisciplinary Optimization, (20):167–189,

2000.

[4] Simo J. C., Wriggers P., and Taylor R.L. A perturbed lagrangian formulation

for the finite element solution of contact problems. Comp. Methods Appl. Mech.

Eng., (50):163–180, 1985.

[5] Luenberger D.G. Linear and Non Linear Programming. Addison-Wesley

Publishing Company Inc., 1984.

[6] Bjorkman G., Klarbring A., and Sjodin B. Sequential quadratic programming

for nonlinear elastic contact problems. Int. J. Numer. Methods Eng., 38:137–

165, 1995.

[7] Barboza H. Algoritmos Numericos para problemas de Contato em elasticidade.

PhD thesis, Universidade Federal de Rio de Janeiro, COPPE, Programa de

Engenharia Civil, 1986.

[8] Herskovits J. A View on Nonlinear Optimization, Advances in Structural

Optimization. Kluver-Academic Publishers., 1995.

65

[9] Herskovits J., Mappa P., Goulart E., and Mota Soares C.M. Mathematical

programming models and algorithms for engineering design optimization.

Computational Methods in Applied Mechanical Engineering, 194:3244–3268,

2005.

[10] Nocedal J. and S.J. Wright. Numerical Optimization. New York, Springer-

Verlag, 1999.

[11] Johnson K.L. Contact Mechanics. Cambrigde University Press, London, 1985.

[12] Kikuchi N. and Oden J.T. Contact Problems in Elasticity: A Study of

Variational Inequalities on Finite Element Methods.

[13] Zienckiewics O.C. The Finite Element Method. McGraw-Hill, 3rd. ed edition,

1977.

[14] Wriggers P. and Fisher K. Recent new developments in contact mechanics. 4th

European LS-DYNA Users Conference, 2003.

[15] Mazorche S. Um algoritmo de ponto interior viavel para problemas de

complementaridade nao linear. Exame de Qualificacao de Doutorado -

PEM/COPPE-UFRJ., 2004.

[16] Mazorche S. and Herskovits J. A new interior point algorithm for nonlinear

complementarity problems. WCSMO6 - CD Proceedings., 38:13–14, 2005.

[17] Timoshenko S.P. and Goodier J.N. Teoria da Elasticidade. McGraw-Hill, 1980.

[18] Auatt S.S., Borges L.A., and Herskovits J. An interior poin optimization

algorithm for contact problems in linear elasticity. Numerical Methods in

Engineering, pages 855–861, 1996.

66

Apendice A

Algumas consideracoes sobre

otimizacao

Um problema de otimizacao consiste na minimizacao ou maximizacao de uma

determinada funcao objetivo sujeito a restricoes em suas variaveis. Isto significa

que matematicamente o problema geral de otimizacao nao linear pode ser formulado

como:

minimize f(x)

sujeito a: g(x) ≥ 0

h(x) = 0

(A.1)

Onde na Eq. A.1 x ∈ IRn e vetor de variaveis do problema ou variaveis de projeto,

f : IRn → IR e a funcao objetivo, g : IRn → IRm e a funcao que define as restricoes

de desigualdade e h : IRn → IRp e a funcao das restricoes de igualdade. As funcoes

f , g e h devem ser funcoes contınuas com derivadas contınuas.

Definicao A.0.1 (Regiao viavel) A regiao viavel Ω e o conjunto de pontos que

satisfazem as seguintes restricoes: Ω = x | g(x) ≥ 0, h(x) = 0

Definicao A.0.2 (Mınimo local) x∗ e um mınimo local do problema A.1 se existe

uma vizinhanca N de x∗ tal que f(x) ≥ f(x∗) para todo ponto x pertencente a Ω∩N .

67

Definicao A.0.3 (Conjunto de restricoes ativas) Dado um ponto x ∈ Ω, o conjunto

de restricoes de desigualdade ativas e: A(x) = i ∈ 1, ..,m | gi(x) = 0

Definicao A.0.4 (Independencia Linear) Dado um ponto x ∈ Ω, as restricoes

satisfazem o requisito de independencia linear em x se os gradientes das restricoes

de igualdade e de desigualdade ativas sao linearmente independentes, ou seja, o

conjunto ∇gi(x) | i ∈ A(x),∇hi(x) | i ∈ 1, .., p e linearmente independente.

Definicao A.0.5 (Condicoes de Karush - Kuhn - Tucker) Seja x ∈ IRn, se

satisfazem as condicoes de primeira ordem de Karush-Kuhn-Tucker (KKT) em x

se existem vetores λ ∈ IRm e µ ∈ IRp que verificam:

∇f(x) −m

∑

i=1

λi∇gi(x) −

p∑

i=1

µi∇hi(x) = 0 (A.2)

gi(x)λi = 0, i = 1, ..,m (A.3)

h(x) = 0 (A.4)

g(x) ≥ 0 (A.5)

λ ≥ 0 (A.6)

Teorema A.0.1 (Teorema de Karush - Kuhn - Tucker) Seja x um mınimo local

do problema A.1 onde se satisfaz o requisito de Independencia Linear, entao se

satisfazem as condicoes de primeira ordem de KKT no ponto x e vetores λ e µ sao

unicos.

Demonstracao: Ver na referencia [10].

Adicionando a funcao objetivo uma combinacao linear das restricoes, temos a

funcao Lagrangeana L(u, λ):

68

Definicao A.0.6 (Funcao Lagrangeana) A funcao lagrangeana para o problema A.1

e a funcao L : IRn × IRm × IRp → IR definida como:

L(x, λ, µ) = f(x) −m

∑

i=1

λigi(x) −

p∑

i=1

µihi(x) (A.7)

Em notacao vetorial, pode expressar como:

L(u, λ) = Π(u) + λtg(u) + µth(x) (A.8)

Sendo o vetor λ o multiplicador de Lagrange e g(u) o vetor restricoes, satisfazem-

se as seguintes equacoes:

∂L(u, λ)

∂u= 0 (A.9)

∂L(u, λ)

∂λ≤ 0 (A.10)

λtg(u) = 0 (A.11)

λ ≥ 0 (A.12)

Estas equacoes representam, de outra forma, as condicoes necessarias (KKT) de

tal forma que u seja um otimo local de (P) em termos da funcao Lagrangeana

Definicao A.0.7 (Hessiana da funcao lagrangeana) A hessiana da funcao la-

grangeana para o problema A.1 e a funcao H : IRn × IRm × IRp → IRn×n definida

como:

H(x, λ, µ) = ∇2xxL(x, λ, µ) = ∇2f(x) −

m∑

i=1

λi∇2gi(x) −

p∑

i=1

µi∇2hi(x) (A.13)

69

Apendice B

Metodo de Newton Raphson

E um metodo iterativo para a procura de zeros numa funcao diferenciavel f . O zero

da tangente e proximo ao zero da curva, portanto e utilizado como novo ponto para

uma nova tangente. Assim o metodo estima as raızes de uma funcao χ : IRm → IRn

partindo de um ponto qualquer da funcao e repete o processo, ate atingir a solucao.

Eventualmente pode divergir se a estimativa inicial nao e boa, o qual representa o

ponto fraco do metodo.

Entao, se procuras-se a solucao da seguinte equacao,

χ(x) = 0, (B.1)

a partir dos primeiros termos da serie de Taylor obtem-se:

χ(x) ≈ χ(x) + ∇χ(xk)t(x− xk), (B.2)

e uma estimativa de xk+1 pode ser obtida de,

χ(x) + ∇χ(xk)t(xk+1 − xk) = 0. (B.3)

Definido,

dk+1 := (xk+1 − xk), (B.4)

obtem-se,

[∇χ(xk)]tdk = −χ(xk). (B.5)

70

Figura B.1: Direcao de busca do algoritmo de Newton

Um resumo do algoritmo e dado por:

Passo 0: Valor inicial estimado

x = x0

Passo 1: Calculo do novo passo

[∇χ(xk)]tdk = −χ(xk)


xk+1 = xk + dk

Passo 3: Criterio de parada

‖xk+1 − xk‖ ≤ ε

Passo 4: Volta ao passo 1

71

Apendice C

Codigos da interface

writeINP.exe

!################################################################################! Code : Fortran 90! Author : Gabriel Mario Guerra! Modify : 01/01/2006!################################################################################

!-------------------------------------------------------------------------------- program writeINP!--------------------------------------------------------------------------------

implicit none integer :: ierror,n,NN,Num character(100) :: string,str real*8,allocatable :: MatF(:,:) integer :: l,s,s1,s2 integer :: p=0

!---------------------------------------------------- Read size of matrix F.txt -

open(UNIT=1, FILE='C:\ABAQUS\F.txt' , STATUS='OLD' ,IOSTAT=ierror)

do read(1,*, IOSTAT=ierror) Num if (ierror/=0) then;do n=1,NN+2; backspace(UNIT=1);enddo;exit; end if NN=NN+1 enddo

allocate(MatF(NN,4)) close(1)

call readf(NN,MatF) !allocate matrix F to MatF

!------------------------------------------------------------ Modify inp file F -

open(UNIT=1, FILE='C:\ABAQUS\Job-1.inp' , STATUS='OLD' ,IOSTAT=ierror) open(UNIT=2, FILE='C:\ABAQUS\Job-1F.inp', STATUS='REPLACE',IOSTAT=ierror) open(UNIT=12, STATUS='SCRATCH',IOSTAT=ierror)

if (ierror/=0) then; print*,'error na lectura',ierror;stop; endif

do

read(1,'(a100)',IOSTAT=ierror) string if (ierror/=0) then; exit;endif

if (string=='** BOUNDARY CONDITIONS') then

! seek pickedset in boundary and write ! to picked.txt after erase

72

do read (1,'(a30)') string s = scan(string,',') if (s/=0) then; write(12,'(a)') string(2:s-1); endif if (string=='** LOADS') then;exit; endif end do

endif

end do

rewind(unit=12) !rewind picked.txt rewind(unit=1) !rewind Job-1.inp

!--------------------------------------------------------------------------------

do ! principal loop

read(1,'(a100)',IOSTAT=ierror) string if (ierror/=0) then; exit;endif

write(2,'(a)') string

if (string=='*End Instance') then

do read (1,'(a100)') string;

s1 = scan(string,'_')

if (s1/=0) then; s2 = scan(string(s1:s1+16),',')+s1-2;

do ! seek in picked.txt read (12,'(a)',IOSTAT=ierror) str

if (ierror/=0) then; rewind(unit=12) ;exit;endif

if (string(s1+1:s2)==str) then write(2,'(a)') string do read (1,'(a100)') string if (string(1:1)=='*') exit write (2,'(a)') string enddo backspace(unit=1)

endif

end do

endif

if (string=='*End Assembly') then;

! write forces do n=1,size(MatF(:,1)) write(2,'(a16,i6,a29)') '*Nset, nset=node',MatF(n,1), ', internal, instance=PART-1-1' write(2,'(i6)') MatF(n,1) enddo

write(2,'(a)') '*End Assembly' exit

endif

end do

if (string=='*End Assembly') read(1,'(a100)',IOSTAT=ierror) string

endif

73

!-------------------------------------------------------------- Write only forces

if (string=='** LOADS') then do n=1,size(MatF(:,1)) write(2,'(a6)') '*Cload'

if (MatF(n,2)/=0) then;write(2,'(a4,i6,a4,f12.6)') 'node',MatF(n,1),', 1,',MatF(n,2);endif if (MatF(n,3)/=0)then ;write(2,'(a4,i6,a4,f12.6)') 'node',MatF(n,1),', 2,',MatF(n,3);endif if (MatF(n,4)/=0)then ;write(2,'(a4,i6,a4,f12.6)') 'node',MatF(n,1),', 3,',MatF(n,4);endif

enddo

do read(1,'(a100)',IOSTAT=ierror) string if (string=='** OUTPUT REQUESTS') then;write(2,'(a)') string;exit;endif enddo

endif

!---------------------------------------------------- Write lines to generate .dat

if (string=='*Output, field') then; write(2,'(a)') '*Node Print';write(2,'(a)') 'U,'; endif

!---------------------------------------------------------------------------------

end do

close (1) close (2)

!-------------------------------------------------------------- Modify *.inp with N

open(UNIT=1, FILE='C:\ABAQUS\Job-1.inp' , STATUS='OLD' ,IOSTAT=ierror) open(UNIT=2, FILE='C:\ABAQUS\Job-1N.inp', STATUS='REPLACE',IOSTAT=ierror)

do read(1,'(a100)',IOSTAT=ierror) string

if (ierror/=0) then; exit;endif write(2,'(a100)') string

if (string=='*Output,field') then;write(2,'(a)')'*Node Print';write(2,'(a)')'U,' endif end do

close (1) close (2)

!----------------------------------------------------------------------------------

contains

subroutine readf(NN,Mat)

implicit none integer :: ierror,n,NN real*8 :: Mat(:,:)

open(UNIT=1, FILE='C:\ABAQUS\F.txt' , STATUS='OLD' ,IOSTAT=ierror) if (ierror/=0) then; print*,'error na lectura',ierror;stop; endif

do n=1,NN read (1,'(i6,3f18.12)') MatF(n,:) if (ierror/=0) then;exit;end if enddo

close (1)

end subroutine

!---------------------------------------------------------------------------------- end program!----------------------------------------------------------------------------------

74

readINP.exe

!################################################################################! Code : Fortran 90! Author : Gabriel Mario Guerra! Modify : 01/01/2006!################################################################################

!-------------------------------------------------------------------------------- program readINP!--------------------------------------------------------------------------------

implicit none integer :: ierror,n character(10) :: string10 integer :: NumNodes = 0 integer :: NumElements = 0 character(60) :: string60 integer :: l real*8 :: tstr,s character(60) :: stringbound character(60) :: stringforce,stf integer :: nro,nb,nf integer :: band=0 integer :: sfa,sfb,sfc = 1 real*8 :: fa,fb,fc =0 real*8 :: vf(3) integer :: bnd=0 character(100) :: string100 character(100) :: string integer :: NumBoundaries integer :: NumForces

!---------------------------------------------------------------------- Read data

open(UNIT=1, FILE='C:ABAQUS\Job-1.inp' , STATUS='OLD' ,IOSTAT=ierror) open(UNIT=2, STATUS='SCRATCH',IOSTAT=ierror) open(UNIT=12, STATUS='SCRATCH',IOSTAT=ierror)

write(2,'(A)' ) '--------------------------------------------------------------' write(2,'(A)' ) '| Preprocessor Abaqus - Matlab Translator |' write(2,'(A)' ) '| Job-1.inp to Job-2.txt |' write(2,'(A)' ) '| Gabriel Guerra 2006 |' write(2,'(A/)' )'------------------------------ -------------------------------'


!-------------------------------------------------------------------------------- do

read (1,'(a10)') string10 select case (trim(string10))

case ("*Node")

do read(1,*, IOSTAT=ierror) n if (ierror/=0) then; do n=0,NumNodes; backspace(UNIT=1);enddo;exit end if; NumNodes=NumNodes+1 enddo

write(2,'(A)' ) 'NumNodes' write(2,'(1x,i6/)') NumNodes write(2,'(A)' ) 'MatNodes' write(2,'(A)' ) '%Number; Coord X; Coord Y; Coord X;'

do n=1,NumNodes read (1,'(a60)') string60 write(2,'(a60)') string60 enddo

write(2,'(A/)' )'--------------------------------------------------------- '

75

case ("*Element, ")

do read(1,*, IOSTAT=ierror) n if (ierror/=0) then; do n=0,NumElements; backspace(UNIT=1);enddo;exit end if; NumElements=NumElements+1 enddo

write(2,'(A)' ) 'NumElements' write(2,'(1x,i6/)' ) NumElements write(2,'(A)' ) 'MatElements' write(2,'(A)' ) '%Number;N1;N2;N3;N4;N5'

do n=1,NumElements read (1,'(a40)') string60 write(2,*) string60 enddo

write(2,'(A/)' )'--------------------------------------------------------- '

case ("*Elastic ")

read (1,'(a40)') string60 write(2,'(A)' ) 'NumMaterial' write(2,'(A/)') '1' write(2,'(A)' ) 'MatMaterial' write(2,'(A)' ) '% Element Materials; E; v; D;' write(2,'(a6,a10,a2)') ' 1 ',string60,',1' write(2,'(A/)' )'--------------------------------------------------------- ' write(2,'(A)' ) 'MatBoundaries' write(2,'(A)' ) '%N; Node; grfx; grfy; grfz'

case ("*Boundary ")

read (1,'(a30)') stringbound l = len_trim(stringbound) s = scan(stringbound,',')

rewind(unit=12)

do read (12,'(a)',IOSTAT=ierror) string60 if (ierror/=0) then;exit;endif if(trim(string60)==stringbound(2:s-1)) then

do read (12,'(i6)',IOSTAT=ierror) nro if (ierror/=0) then;exit;endif nb=nb+1 if ('PINNED'==stringbound(s+2:l)) then write(2,'(i6,i6,a)') nb,nro, ' 0, 0, 0';endif if ('XASYMM'==stringbound(s+2:l)) then write(2,'(i6,i6,a)') nb,nro, ' 1, 0, 0';endif if ('YASYMM'==stringbound(s+2:l)) then write(2,'(i6,i6,a)') nb,nro, ' 0, 1, 0';endif if ('ZASYMM'==stringbound(s+2:l)) then write(2,'(i6,i6,a)') nb,nro, ' 0, 0, 1';endif enddo endif

enddo

case ("*Cload ")

if (bnd==0) then write(2,'(A/)' ) write(2,'(A)' ) 'NumBoundaries' write(2,'(1x,i6)' ) nb write(2,'(A/)' )'--------------------------------------------------------- ' write(2,'(A)' ) 'MatForces' write(2,'(A)' ) '%N; Node; fx; fy; fz' bnd=1

endif

76

read (1,'(a30)') stringforce;backspace(UNIT=1) l = len_trim(stringforce); s = scan(stringforce,',')

read (1,'(13x,i4,f12.6)') sfa,fa;vf(sfa)=fa; read (1,'(a1)') string10;backspace(UNIT=1)

if (string10=='_') then read (1,'(13x,i4,f12.6)') sfb,fb;vf(sfb)=fb; endif

read (1,'(a1)') string10;backspace(UNIT=1)

if (string10=='_') then read (1,'(13x,i4,f12.6)') sfc,fc;backspace(UNIT=1);vf(sfc)=fc; endif

rewind(unit=12)

do read (12,'(a)',IOSTAT=ierror) string60

if (ierror/=0) then;exit;endif

if(trim(string60)==stringforce(2:s-1)) then do read (12,'(i4)',IOSTAT=ierror) nro if (ierror/=0) then;exit;endif nf=nf+1 write(2,'(i6,i6,3f12.6)') nf,nro,vf enddo endif

enddo

sfa=1;sfb=1;sfc=1;fa=0;fb=0;fc=0;vf=0

case ("*Nset, nse")

backspace(UNIT=1) read (1,'(a100)') string100 write(12,'(a)') string100(14:10+scan(string100(12:len_trim(string100)),','))

do read (1,'(a100)') string100 if (string100(1:1)=='*')then;backspace(UNIT=1);exit;end if l=len_trim(string100) s=scan(string100,',') tstr=nint(l/s) do n=1,tstr write(12,*) string100(s*n-(s-1):s*n-1) enddo enddo

end select

if (string10=="*End Step") exit

enddo

write(2,'(A/)' ) write(2,'(A)' )'NumForces' write(2,'(1x,i6)' ) nf write(2,'(A/)' )'---------------------------------------------------------------- '

close(1); close(12)

rewind(unit=2)

open(UNIT=3, FILE='C:ABAQUS\Job-matlab.txt', STATUS='REPLACE',IOSTAT=ierror)


77

do read (2,'(a80)',IOSTAT=ierror) string100 if (ierror/=0) then;exit;endif

if (string100=='MatBoundaries') then

do read (2,'(a80)',IOSTAT=ierror) string100 if (string100=='NumBoundaries') then ; read (2,'(i10)',IOSTAT=ierror) NumBoundaries; exit endif enddo

write(3,'(a)') 'NumBoundaries' write(3,'(i6)') NumBoundaries

do n=1,NumBoundaries+6; backspace(UNIT=2);enddo

do n=1,NumBoundaries+3

read (2,'(a80)',IOSTAT=ierror) string100 write(3,'(a80)') string100;enddo

do n=1,3;read (2,*) string100;enddo

endif

if (string100=='MatForces') then

do read (2,'(a80)',IOSTAT=ierror) string100 if (string100=='NumForces') then ; read (2,'(i10)',IOSTAT=ierror) NumForces; exit endif enddo

write(3,'(a)') 'NumForces' write(3,'(i6)') NumForces

do n=1,NumForces+6; backspace(UNIT=2);enddo

do n=1,NumForces+3 read (2,'(a80)',IOSTAT=ierror) string100 write(3,'(a80)') string100; enddo

do n=1,3;read (2,*) string100;enddo

endif

write(3,'(a80)') string100

enddo

close(2); close(3)

!-------------------------------------------------------------------------------- end program!--------------------------------------------------------------------------------

78

readODB.py

################################################################################ Code : Python Author : Gabriel Mario Guerra Modify : 21/03/2006################################################################################

from sys import argv

from odbAccess import *

odb = argv[1]

u = argv[2]

precisao = '%.15e'

local = openOdb(odb)

u = open(u,'w')

for stepName in local.steps.keys():

u.write('**Step\n')

disp = local.steps[str(stepName)].frames[-1].fieldOutputs['U'].values

for v in disp:

for data in v.data:

line = str(precisao % data)

u.write(line+'\n')

u.close()

local.close()

################################################################################

79

Documents

TECNICAS DE OTIMIZAÇ´ AO APLICADAS A PROBLEMAS DE … · elásticos, na hipótese de elasticidade infinitesimal sem atrito. Trata-se de um problema não linear devido à presença