Anderson Pereira Projeto ótimo de pórticos planos com ... otimização de estruturas, programação matemática, análise não-linear e elementos finitos. ... tais como as que utilizam

Anderson Pereira

Projeto ótimo de pórticos planos com restrição à flambagem

Dissertação apresentada como requisito parcial para obtenção do grau de Mestre pelo Programa de Pós-graduação em Engenharia Civil do Departamento de Engenharia Civil da PUC-Rio. Área de Concentração: Estruturas.

Orientador: Luiz Eloy Vaz

Co-orientador: Paulo Batista Gonçalves

Rio de Janeiro

Agosto de 2002

DBD

PUC-Rio - Certificação Digital Nº 0025018/CA

Anderson Pereira

Projeto ótimo de pórticos planos com restrição à flambagem

Dissertação apresentada como requisito parcial para obtenção do grau de Mestre pelo Programa de Pós-graduação em Engenharia Civil do Departamento de Engenharia Civil do Centro Técnico Científico da PUC-Rio. Aprovada pela Comissão Examinadora abaixo assinada.

Prof. Luiz Eloy Vaz Orientador

Departamento de Engenharia Civil – PUC-Rio

Prof. Paulo Batista Gonçalves Co-orientador

Departamento de Engenharia Civil – PUC-Rio

Prof. Raul Rosas e Silva Departamento de Engenharia Civil – PUC-Rio

Profa. Cláudia Ribeiro Eboli Universidade Federal do Rio de Janeiro

Prof. Ney Augusto Dumont Departamento de Engenharia Civil – PUC-Rio

Rio de Janeiro, 29 agosto de 2002

DBD


Todos os direitos reservados. É proibida a reprodução total ou parcial do trabalho sem autorização da universidade, do autor e do orientador.

Anderson Pereira Graduou-se em Engenharia Civil na UDESC/Joinville (Universidade do Estado de Santa Catarina) em 2000. Interesses acadêmicos em áreas de pesquisa que envolvam otimização de estruturas, programação matemática, análise não-linear e elementos finitos. Atualmente é aluno de doutorado no departamento de engenharia civil da PUC-Rio e pesquisador do Tecgraf – Grupo de tecnologia em computação gráfica.

Ficha Catalográfica

Pereira, Anderson Projeto ótimo de pórticos planos com restrição à flambagem / Anderson Pereira; orientador: Luiz Eloy Vaz; co-orientador: Paulo Batista Gonçalves. – Rio de Janeiro : PUC, Departamento de Engenharia Civil, 2002. v., 99 f.: il. ; 29,7 cm 1. Dissertação (mestrado) – Pontifícia Universidade Católica do Rio de Janeiro, Departamento de Engenharia Civil. Inclui referências bibliográficas 1. Engenharia Civil – Teses. 2. Otimização. 3. Análise Não Linear. 4. Análise de Sensibilidade. 5. Pórticos Planos. 6. Programação Matemática. 7. Instabilidade. 8. Projeto Ótimo. I. Vaz, L. E. (Luiz Eloy); Gonçalves, P. B. (Paulo Batista). II Pontifícia Universidade Católica do Rio de Janeiro. Departamento de Engenharia Civil. III. Título.

DBD


Dedico este trabalho, que se torna pe-queno diante do sofrimento que ela tem passado, a minha querida mãe, Idete de Souza Pereira, com muito amor.

DBD


Agradecimentos Aos meus orientadores Luiz Eloy Vaz e Paulo Batista Gonçalves pelo estímulo e

parceria para a realização deste trabalho.

Aos professores do Departamento de Engenharia Civil da PUC-Rio, por todos os

conhecimentos transmitidos durante a pós-graduação.

A Janaína, pelo carinho, compreensão, paciência e constante apoio que fizeram

possível a realização deste trabalho.

Ao amigo Sidiclei, a quem responsabilizo pelo inicio da minha vida acadêmica,

graças ao seu incentivo e exemplo.

A toda minha família, em especial meus pais, Osvy Manoel Pereira e Idete de

Souza Pereira, por todos os gestos de carinho e otimismo.

Aos colegas de república, Alan Wilter S. da Silva, Antonio Miranda e Ricardo

Alexandre de Oliveira Passos pelo convívio saudável e pelo ambiente de estudo.

Aos colegas “ótimos”, Ivy e Sandoval, pelas discussões ao longo deste trabalho.

A Claudia Eboli, pelas aulas sobre o algoritmo de Han-Powell.

A todos os amigos da PUC, em especial Galvão, Joabson, Antonio Sérgio, Chan,

Salete, Walter e Jaqueline.

Aos funcionários da PUC-Rio, em particular a Ana Roxo, Lenilson, Cristiano,

Euclides, Haroldo, José Nilson e Evandro.

Ao CAPES, pelo apoio financeiro.

DBD


Resumo

Pereira, Anderson; Vaz, Luiz Eloy; Gonçalves, Paulo Batista. Projeto ótimo de pórticos planos com restrição à flambagem. Rio de Janeiro, 2002. 99p. Dissertação de Mestrado - Departamento de Engenharia. Pontifícia Universidade Católica do Rio de Janeiro.

O objetivo deste trabalho é apresentar uma formulação e uma

correspondente implementação computacional para otimização de dimensões de

estruturas evitando os problemas de instabilidade apresentados pela formulação

convencional. Para atingir este objetivo, a formulação utilizada considera os

efeitos da não-linearidade geométrica no comportamento da estrutura e inclui uma

restrição sobre a carga de colapso. Elementos finitos reticulados planos e a

formulação Lagrangiana Atualizada forma utilizados para análise de estruturas

com comportamento geometricamente não-linear. As varáveis de projeto são as

alturas das seções transversais dos elementos. O método de Newton-Raphson é

utilizado acoplado a diferentes estratégias de incremento de carga e de iteração,

tais como as que utilizam a restrição do comprimento de arco e as baseadas no

controle dos deslocamentos generalizados, que permitem a ultrapassagem de

pontos críticos que possam existir ao longo da trajetória de equilíbrio. Os

algoritmos de programação matemática utilizados neste trabalho empregam os

gradientes da função objetivo e das restrições, que são calculados com base nos

gradientes das respostas da estrutura. Partindo-se das equações gerais de equilíbrio

válidas para qualquer elemento, foram desenvolvidas expressões analíticas

aproximadas que permitem o cálculo das sensibilidades em relação as variáveis de

projeto aproveitando as características da análise.

Palavras-chave Otimização; análise de sensibilidade; programação matemática; projeto

ótimo; pórticos planos; análise não-linear; instabilidade.

DBD


Abstract

Pereira, Anderson; Vaz, Luiz Eloy; Gonçalves, Paulo Batista. Optimal design of planar frames with stability constraints. Rio de Janeiro, 2002. 99p. Msc. Dissertation - Departamento de Engenharia. Pontifícia Universidade Católica do Rio de Janeiro.

The aim of this work is to present a formulation and corresponding

computational implementation for the sizing optimization of structures. To

achieve this goal, the formulation considers the structural geometric nonlinear

behavior and include a constraint related to the collapse load. Plane frame finite

elements and Updated Lagrangian approach are used for the geometric nonlinear

analysis. The standard Newton-Raphson method, in connection with different load

increment strategies and iteration, such as use the arch length method and

strategies based on the control of generalized displacements, which allow the

algorithm to transpose the critical points that happen to appear along the

equilibrium path. The mathematical programming algorithms applied in this work

make use of the gradients of the objective function and of the constraints, which

depend on the gradients of the structural response. Starting from general

equilibrium equations for the Update Lagrangian approach, valid for any finite

element, approximate analytical expressions for the sensitivity analysis whit

respect of design variables were developed taking advantage of the structural

characteristics.

Keywords Optimization; sensitivity analysis; mathematical programming; optimal design; planar frames; nonlinear analysis; stability.

DBD


Sumário

1 Introdução..........................................................................................16

1.1 Escopo do trabalho ...................................................................................18

2 Programação Matemática – Princípios Básicos.................................19

2.1 Considerações Gerais ..............................................................................19

2.2 Condições de Ótimo .................................................................................20

2.3 Forma Geral dos Algoritmos de Otimização .............................................21

2.4 Método de Newton para Problemas de Otimização sem Restrição .........22

2.5 Busca Linear .............................................................................................23

2.6 Programação Quadrática..........................................................................24

2.7 Algoritmo de Han-Powell - Programação Quadrática Seqüencial ............25 2.7.1 Etapas do Algoritmo Não-Linear Han-Powell (PQS).....................................27

2.8 Método dos Pontos Interiores ...................................................................29 2.8.1 Etapas do Algoritmo de Pontos Interiores (PI)..............................................32

2.9 Implementação .........................................................................................33

3 Análise Não-Linear Geométrica.........................................................34

3.1 Comentários Iniciais..................................................................................34

3.2 Comportamento Não-Linear, Análise e Projeto ........................................34 3.2.1 Fontes de Não-Linearidade...........................................................................35

3.3 Formulação para a Análise Não-Linear Geométrica de Estruturas

Reticuladas..........................................................................................................36 3.3.1 Descrição do Problema .................................................................................36 3.3.2 Principio dos Deslocamentos Virtuais ...........................................................38

DBD


3.3.3 Aplicação a Elementos de Pórtico Plano ......................................................41 3.3.4 Montagem das Equações da Estrutura.........................................................46

3.4 Estratégias de Solução para Problemas Não-Lineares ............................47 3.4.1 Análise Incremental-Iterativa.........................................................................48 3.4.2 Estratégias de Iteração .................................................................................52 3.4.3 Critérios de Convergência.............................................................................53 3.4.4 Incremento Automático de Carga..................................................................55

3.5 Determinação dos Pontos Críticos ...........................................................57

3.6 Exemplos ..................................................................................................58 3.6.1 Pórtico de Lee ...............................................................................................59 3.6.2 Pórtico de Williams........................................................................................61

4 Análise de Sensibilidade....................................................................64


4.2 Método Analítico .......................................................................................64 4.2.1 Sensibilidade dos Deslocamentos ................................................................65 4.2.2 Método Analítico Aproximado (MAA) ............................................................67 4.2.3 Sensibilidade da Carga Limite.......................................................................69

4.3 Método das Diferenças Finitas (MDF) ......................................................71

4.4 Exemplos ..................................................................................................71 4.4.1 Pórtico de Lee ...............................................................................................72 4.4.2 Pórtico de Williams........................................................................................73

5 Otimização de Dimensões .................................................................75


5.2 Metodologia de Otimização ......................................................................76

5.3 Formulação do problema de Otimização ..................................................77 5.3.1 Fatores de escala ..........................................................................................78 5.3.2 Cálculo dos gradientes..................................................................................79

5.4 Implementação da Formulação ................................................................81 5.4.1 Programa Principal ........................................................................................81

DBD


5.4.2 Processo de Otimização ...............................................................................82

5.5 Exemplos ..................................................................................................85 5.5.1 Pórtico de Williams........................................................................................86 5.5.2 Pórtico de Lee ...............................................................................................88 5.5.3 Pórtico com três barras .................................................................................90 5.5.4 Pórtico com 7 barras .....................................................................................93

6 Comentários Finais............................................................................95

6.1 Sugestões.................................................................................................96

7 Referências Bibliográficas .................................................................98

DBD


Lista de figuras Figura 3.1. Deformações do elemento. ...................................................36

Figura 3.2. Elemento de pórtico. .............................................................43

Figura 3.3. Elemento de pórtico. .............................................................44

Figura 3.4. Curva carga-deslocamento. ..................................................48

Figura 3.5. Solução incremental-iterativa: sistema com um grau de liberdade...................................................................................................51

Figura 3.6. Variação do sinal do parâmetro de rigidez generalizado (GSP). ......................................................................................................57

Figura 3.7. Pontos críticos de uma estrutura...........................................58

Figura 3.8. Pórtico de Lee. ......................................................................59

Figura 3.9. Pórtico de Lee – malhas utilizadas........................................59

Figura 3.10. Pórtico de Lee – curvas de equilíbrio. .................................60

Figura 3.11. Pórtico de Lee – configuração deformada...........................60

Figura 3.12. Pórtico de Williams..............................................................61

Figura 3.13. Pórtico de Williams – curvas de equilíbrio...........................62


Figura 4.2. Pórtico de Williams................................................................73

Figura 5.1. Situações encontradas na análise estrutural.........................76

Figura 5.2. Fluxograma do programa principal........................................82

Figura 5.3. Forma geral dos algoritmos de PM. ......................................83

Figura 5.4. Fluxograma da função Análise/Sensibilidade........................84

Figura 5.5. Pórtico de Williams................................................................86


Figura 5.7. Alturas das seções transversais / momento fletor. ................89

Figura 5.8. Pórtico com três barras. ........................................................90

Figura 5.9. Deslocamento horizontal do ponto a para caso linear e não-linear. .......................................................................................................92

Figura 5.10. Pórtico com 7 barras. ..........................................................93

DBD


Lista de tabelas Tabela 2.1 Divisão dos problemas de Programação Matemática ............21 Tabela 3.1 Pórtico de Lee – cargas críticas .............................................61 Tabela 3.2 Pórtico de Williams – cargas críticas ......................................63 Tabela 4.1 Pórtico de Lee – deslocamentos e sensibilidades para h1 .....73 Tabela 4.2 Pórtico de Lee – deslocamentos e sensibilidades para h2 .....73 Tabela 4.3 Pórtico de Lee – cargas críticas e sensibilidades...................73 Tabela 4.4 Pórtico de Williams – deslocamentos e sensibilidades para h1.................................................................................................................74 Tabela 4.5 Pórtico de Williams – cargas críticas e sensibilidades............74 Tabela 5.1 Valores usuais dos parâmetros. .............................................85

Tabela 5.2 Pórtico de Williams – resumo dos resultados.........................86

Tabela 5.3 Pórtico de Williams – dimensões finais ..................................87

Tabela 5.4 Pórtico de Williams – resumo dos resultados.........................89

Tabela 5.5 Pórtico de Lee – dimensões finais - pilar................................89

Tabela 5.6 Pórtico de Lee – dimensões finais - viga ................................89

Tabela 5.7 Pórtico com três barras – nós restritos ...................................91

Tabela 5.8 Pórtico com três barras – resumo dos resultados ..................91

Tabela 5.9 Pórtico com três barras – dimensões finais - pilares ..............91

Tabela 5.10 Pórtico com três barras – dimensões finais - vigas ..............91

Tabela 5.11 Pórtico com 7 barras – resumo dos resultados ....................93

Tabela 5.12 Pórtico com 7 barras – dimensões finais..............................94

DBD


Lista de Símbolos, Siglas e Abreviaturas LETRAS ROMANAS a,b,...z - nós do elemento ai - gradientes da restrição ci A - matriz dos gradientes das restrições A - área da seção transversal do elemento b0 - valor inicial da diagonal da aproximação da Hessiana bz , gz - parâmetros do algoritmo SQP b - vetor dos termos independentes das restrições, largura de uma

seção retangular B - aproximação da Hessiana B0 - aproximação inicial da Hessiana ci - restrição do problema de otimização c - restrição do problema de otimização adimensionalizada c - vetor das restrições C - tensor incremental de tensão-deformação C - matriz diagonal que contem os valores das restrições d - direção de busca d0, d1 - direções de busca intermediárias do algoritmo PI e - parcela linear da deformação E - módulo de elasticidade f - função objetivo f - função objetivo adimensionalizada

f - vetor de forças internas do elemento F - vetor de forças internas da estrutura F - vetor de pseudo-forças g - gradiente da função objetivo h - altura de uma seção transversal retangular h - variável adimensional relativa à altura da seção H - Hessiana da função objetivo I - matriz identidade k - contador de iterações do algoritmo de PM ka, ke, kf - parâmetros do algoritmo de pontos interiores ke - matriz de rigidez elástica do elemento kg - matriz de rigidez geométrica do elemento K(T) - matriz de rigidez tangente da estrutura l - número de restrições de igualdade, comprimento do elemento L - comprimento do elemento L(x, λ) - função Lagrangiana m - número total de restrições M - momento fletor n - número de variáveis de projeto nr - número de iterações para reinício da matriz B N1, N3 - funções de forma ne - número de elementos

DBD


ndr - número de deslocamentos restritos nsecs - número de grupos de seções transversais p - ponto material p(t) - função unidimensional utilizada na busca linear P - força axial P, P1, P2 - forças aplicadas nos elementos P - vetor de carregamento externo

refP - vetor de carregamento externo de referência W - trabalho W - Hessiana da função Lagrangiana w - contador de iterações da análise não-linear q - intensidade da carga distribuída q - vetor da equação de PQ Q - força cisalhante Q - matriz da equação de PQ y - coordenada cartesiana ri - fator de penalidade da busca linear R - trabalho virtual externo R - vetor de forças residuais ℜ - conjunto de números reais t - tamanho do passo (busca linear) tol1, tol2 - tolerâncias do algoritmo PQS tol - tolerância do algoritmo PI T - matriz de rotação T - tensor de tensões de Cauchy T - tensor de tensões Piola-Kirchhoff II

,u v - componentes de deslocamento u,v - componentes de deslocamento u - vetor de deslocamentos nodais ug - parcela de u referente às forças residuais g ur - parcela de u referente às forças de referência Pref. ut - vetor dos deslocamentos nodais tangenciais uj,lim - valor absoluto admissível para o deslocamento uj V - volume x - vetor das variáveis de projeto, coordenada cartesiana x, y, z - coordenadas cartesianas x0 - vetor inicial das variáveis de projeto x - vetor das variáveis de projeto adimensionalizada LETRAS GREGAS α - seqüência de valores utilizados na busca linear

( )δ ⋅ - correção iterativa, indicador variacional ( )∆ ⋅ - incremento

ε ,ε ij - deformações de Green-Lagrange γ - parâmetro de controle da busca linear η - perturbação relativa (diferenças finitas), parcela não-linear da

deformação

DBD


ϕ - autovetor correspondente ao autovalor nulo da matriz de rigidez tangente

λ - parâmetro de carga responsável pelo escalonamento de Pref λi - multiplicador de Lagrange da restrição ci λ - vetor dos multiplicadores de Lagrange λ0, λ1, λ - estimativa do vetor dos multiplicadores de Lagrange no algoritmo

de pontos interiores iλ - valor admissível para a carga crítica λ - estimativa do multiplicador de Lagrange no algoritmo de pontos

interiores Λ - matriz diagonal contendo os multiplicadores de Lagrange θ - ângulo (coordenada polar) ρ - coeficiente de deflexão da direção de busca τ - componentes da tensão de Cauchy ζ - tolerância ao resíduo requerida no processo de convergência ζ1 - fator de convergência baseado em relações de força ζ2 - fator de convergência baseado em relações de deslocamentos SOBRESCRITOS E SUBSCRITOS

( )t ⋅ - função avaliada no instante t *( )⋅ - função avaliada no ponto crítico T( )⋅ - transposta do vetor ou matriz

( )ij⋅ - notação indicial ( )l⋅ - limite inferior da variável de projeto, restrições da carga crítica ( )u⋅ - limite superior da variável de projeto ( )d⋅ - restrições de deslocamento SIGLAS E ABREVIATURAS BFGS - Broyden - Fletcher - Goldfarb - Shanno GSP - Parâmetro de rigidez generalizado MAA - Método Analítico Aproximado MDF - Método das Diferenças Finitas PM - Programação Matemática PQ - Programação Quadrática PQS - Programação Quadrática Seqüencial RLT - Referencial Lagrangiano Total RLA - Referencial Lagrangiano Atualizado

DBD


1 Introdução

Ao longo das últimas décadas, a otimização estrutural vem se destacando

dentro da engenharia. O crescente desenvolvimento dos microcomputadores

aliado à automação dos procedimentos de análise viabilizaram a sua aplicação.

Desta forma, as técnicas de otimização numérica se tornaram valiosas na busca

pelo projeto ótimo.

A formulação clássica dos problemas de otimização de dimensões busca

minimizar o peso de uma estrutura, onde é adotada a hipótese do comportamento

linear físico e geométrico, o que simplifica bastante as etapas de análise estrutural

e de avaliação da sensibilidade, aumentando a eficiência e contribuindo para a

popularização da otimização.

Em geral, a formulação clássica gera projetos seguros e econômicos. No

entanto, freqüentemente o modelo linear não permite a avaliação correta da

capacidade de carga da estrutura.

Na busca pela estrutura ideal, utilizando-se a formulação clássica, muitas

vezes se aumenta a esbeltez de um dado elemento estrutural, o que pode gerar

significativas mudanças no seu mecanismo de colapso. Desta forma, a não-

linearidade geométrica se torna cada vez mais importante e dá origem a vários

fenômenos que não são encontrados em sistemas lineares como à existência de

múltiplas configurações de equilíbrio (estáveis e instáveis) e de pontos críticos

(limite e bifurcação) ao longo do caminho não-linear de equilíbrio.

Assim, durante muito tempo a otimização foi criticada por gerar estruturas

com sérios problemas de instabilidade. Contudo, estes problemas não são

inerentes à otimização, podendo ser evitados através da utilização de restrições e

procedimentos de análise apropriados. A formulação do problema é de suma

importância na solução encontrada, pois nenhuma solução será mais “ótima” do

que a sua formulação permite.

O objetivo principal deste trabalho é apresentar uma metodologia de

otimização de estruturas reticuladas planas com segurança em relação à

instabilidade. Para atingir este objetivo, a formulação proposta considera os

DBD


Introdução

17

efeitos da não-linearidade geométrica no comportamento da estrutura e inclui uma

restrição sobre a carga crítica

Esta dissertação é parte integrante de algumas linhas de pesquisa do

DEC/PUC-Rio, em particular da linha de Aplicação de Técnicas de Otimização,

Instabilidade e Dinâmica das Estruturas.

As técnicas de otimização vêm sendo largamente aplicadas a problemas de

engenharia estrutural no DEC/PUC-Rio e diversos trabalhos vêm sendo

publicados nesta área. O trabalho de Eboli (1989) foi o precursor desta linha e traz

uma descrição detalhada do algoritmo de Han-Powell de programação não-linear.

Mais recentemente, Parente (2000) estudou a otimização de forma de estruturas

geometricamente não-lineares.

Na parte de instabilidade, Silveira (1995) forneceu uma metodologia geral

de solução de sistemas de equações algébricas não-lineares. Utilizando a mesma

metodologia, Rocha (2000), realizou um estudo comparativo de diversas

estratégias de iteração e incremento de carga através da análise de vários

exemplos numéricos de sistemas estruturais e Galvão (2000) implementou

diversas formulações de elementos finitos geometricamente não-lineares.

Os algoritmos de programação matemática utilizados neste trabalho

empregam os gradientes da função objetivo e das restrições em relação às

variáveis de projeto para determinar a direção de busca do processo de

otimização. Esses gradientes são calculados com base nos gradientes das respostas

da estrutura (sensibilidades). Portanto, a convergência do processo de otimização

é fortemente influenciada pela qualidade das sensibilidades calculadas. A análise

de sensibilidade de estruturas geometricamente não-lineares tem apresentado

grandes progressos nos últimos anos e as equações básicas já são bem conhecidas

(Kleiber, 1997; Haftka, 1992).

As expressões necessárias ao cálculo analítico das sensibilidades são

apresentadas no presente trabalho. É importante ressaltar que os resultados obtidos

através da aplicação dessas expressões representam as sensibilidades exatas de

uma dada malha de elementos finitos.

DBD


Introdução

18

1.1 Escopo do trabalho

Para facilitar o entendimento, a dissertação foi dividida em diversos

capítulos, cujo conteúdo é apresentado a seguir.

No Capítulo 2 são apresentados os conceitos gerais de programação

matemática e os algoritmos utilizados neste trabalho. Este capítulo mostra quais as

informações necessárias para os algoritmos de otimização, com o objetivo de

facilitar o entendimento da organização do sistema computacional implementado.

No Capítulo 3 estuda-se a análise não-linear de estruturas através do Método

dos Elementos Finitos. A formulação dos elementos de pórtico de acordo com o

procedimento Lagrangiano Atualizado é discutida em detalhes. Os métodos de

determinação do caminho de equilíbrio de estruturas não-lineares são

apresentados, bem como o método utilizado para a determinação dos pontos

críticos. Por fim, são apresentados exemplos numéricos que mostram a qualidade

dos elementos finitos e procedimentos numéricos implementados.

No Capítulo 4 é feita uma explanação sobre o cálculo dos gradientes dos

deslocamentos e da carga crítica (limite) de um modelo de elementos finitos com

comportamento geometricamente não-linear onde, considerando-se os

procedimentos de análise da análise não-linear, algumas aproximações são feitas

tendo em vista uma melhor eficiência computacional. São apresentados exemplos

numéricos com o objetivo de validar as implementações realizadas e de comparar

a precisão do método.

No Capítulo 5 é apresentada a formulação do modelo de otimização de

dimensões de estruturas bidimensionais. O modelo proposto inclui restrições

sobre deslocamentos e carga crítica, podendo ser utilizado tanto para estruturas

lineares quanto não-lineares. Também são apresentados os programas que fazem

parte do sistema computacional para otimização de dimensões e exemplos de

estruturas otimizadas por este sistema.

Finalmente, no Capítulo 6 são apresentadas as conclusões obtidas neste

trabalho e as sugestões para trabalhos futuros.

DBD


2 Programação Matemática – Princípios Básicos

2.1 Considerações Gerais

Os objetivos deste capítulo são apresentar os conceitos de Programação

Matemática (PM) necessários à compreensão do processo de otimização de

dimensões e descrever os algoritmos de otimização utilizados.

Em problemas típicos de engenharia, podem ser obtidas várias, ou

possivelmente infinitas, soluções. Em um problema de otimização deseja-se obter

um projeto ótimo, maximizando ou minimizando uma função a qual

denominamos função objetivo. Isto deve ser realizado através da determinação dos

parâmetros que definem o sistema. Estes parâmetros são chamados de variáveis

de projeto. Na maioria dos problemas encontraremos restrições impostas para que

o projeto seja admissível ou viável, devido às leis físicas da natureza, leis

políticas, limitações de orçamento, etc.

A Programação Matemática é a disciplina que estuda a minimização de

funções em problemas com ou sem restrições. Matematicamente, estes problemas

são enunciados como:

Minimizar ( )f x nx∈ℜ sujeito a ( ) = 0ic x 1...i l= ( ) 0ic x ≤ 1...i l m= + l u

i i ix x x≤ ≤ 1...i n=

(2.1)

onde x é um ponto do nℜ sobre o qual são impostos os limites mínimos e

máximos (restrições laterais), ( )f x é a função a ser minimizada e as funções

( )ic x representam as restrições de igualdade e desigualdade. Assume-se que tanto

a função objetivo quanto as restrições são funções contínuas no nℜ . Em geral,

elas são funções não-lineares e implícitas das variáveis ( )x que definem o

problema.

DBD


Programação Matemática - Princípios Básicos

20

Um ponto que satisfaça todas as restrições é denominado um ponto viável e

o conjunto de todos os pontos que satisfaçam todas as restrições é conhecido

como região viável. Uma restrição de desigualdade define uma fronteira que

divide o nℜ em uma região viável e outra inviável. Quando um ponto está sobre

esta fronteira, a restrição é dita ativa; quando um ponto está no interior da região

viável, a restrição está inativa e, quando um ponto está fora desta região, à

restrição está violada.

2.2 Condições de Ótimo

A solução *x do problema enunciado em (2.1) tem que necessariamente

atender as condições de Kuhn-Tucker enunciadas por:

* *( , ) 0xL x∇ =λ *( ) 0ic x = 1...i l= *( ) 0ic x ≤ 1...i l m= + * 0iλ ≥ 1...i l m= + * *( ) 0i ic xλ = i∀

(2.2)

onde * *( , )L x λ é a função Lagrangiana dada pela expressão a seguir:

* * * * *

1( , ) ( ) ( )

l

i ii

L x f x c xλ=

= +∑λ (2.3)

onde *iλ são os multiplicadores de Lagrange associados às restrições no ponto *x ,

solução do problema.

As condições de Kuhn-Tucker são também conhecidas como condições de

primeira ordem. Para determinadas classes de problemas de programação

matemática as condições de Kuhn-Tucker são suficientes para a determinação de

uma solução ótima global. São incluídos nessas classes os problemas de

programação convexa, tais como os de programação linear e quadrática. O

problema de programação convexa é caracterizado por função objetivo e

restrições convexas.

DBD



21

Porém, se o problema não é de programação convexa, o que é mais comum,

as condições de primeira ordem não são mais suficientes para a determinação da

solução ótima global, devendo ser verificada a condição de segunda ordem,

expressa na equação (2.4) a seguir

* 0t ≥d W d , *0 tal que 0t

i∀ ≠ =d d a (2.4)

onde * *( )i ic x= ∇a para todas as restrições ativas e * 2 *( )L x= ∇W é a Hessiana da

função Lagrangiana. O que significa que *W em *x é positiva definida no ponto

ótimo para qualquer direção estacionária d .

2.3 Forma Geral dos Algoritmos de Otimização

Para resolver um problema de otimização, além dos algoritmos ditos

evolucionários, existem diversos algoritmos de programação matemática que são

definidos de acordo com as características da função-objetivo e das restrições.

Assim, os problemas de otimização podem se dividir em diferentes formas, como

mostra a Tabela 2.1.

TABELA 2.1 Divisão dos problemas de Programação Matemática Tipos de Otimização (x)f (x)ic Programação Linear linear linear Programação Quadrática quadrática linear Programação Não-Linear linear / não-linear não-linear / linear

Algoritmos de otimização para problema de programação linear e

programação quadrática têm solução em um número finito de passos, já os

algoritmos de programação não-linear podem não ter solução em um número

finito de passos, mas espera-se que a seqüência gerada convirja (no limite) para

um mínimo local. Portanto, um problema adicional no processo de otimização

ocorre quando a função objetivo e as restrições são funções não-lineares do vetor

de variáveis de projeto, nx∈ℜ .

Os algoritmos de programação não-linear, restrita e irrestrita, são

procedimentos iterativos em que novos pontos x são gerados a partir do ponto

corrente 0x através da expressão:

DBD



22

0x x t= + d (2.5)

Assim, os algoritmos podem ser divididos em duas etapas principais: a

primeira etapa é a determinação da direção de busca d e a segunda é a avaliação

do parâmetro escalar t, que representa o tamanho do passo a ser dado ao longo da

direção de busca. A partir da expressão (2.5) diversos algoritmos podem ser

construídos utilizando diferentes técnicas para a determinação da direção de busca

e do tamanho do passo.

Os algoritmos de PM podem ser classificados de acordo com a ordem da

derivação da função objetivo e das restrições utilizadas para a determinação da

direção de busca. Desta forma, um algoritmo é dito de primeira ordem se utilizar

apenas os gradientes da função objetivo e das restrições para calcular a direção de

busca. Por outro lado, se o algoritmo utiliza informações sobre as Hessianas

destas funções, então ele é dito de segunda ordem.

2.4 Método de Newton para Problemas de Otimização sem Restrição

O método de Newton utiliza a informação de segunda ordem da função a

otimizar. Assim, a sua convergência é quadrática. Com tal propósito a função

( )f x é expandida até a segunda ordem, ou seja, a expansão de Taylor em torno

do ponto 0x será:

21

0 0 0 0 0 02( ) ( ) ( )( - ) ( - ) ( )( - )tf x f x f x x x x x f x x x= +∇ + ∇ (2.6) se

0 0( - )x x x x x= ∆ = → = +d d (2.7) e

0( )f x= ∇g e 20( )f x= ∇H (2.8)

Substituindo-se (2.7) e (2.8) em (2.6), tem-se

t t1

0 0 2( ) ( )f x f x+ = + +d d g d Hd (2.9)

onde d é o incremento de 0x , g é vetor gradiente de f e H, uma matriz simétrica

positiva definida, é a hessiana da função f no ponto 0x . A equação (2.9) é uma

DBD



23

equação quadrática cuja variável é d. Portanto, o algoritmo de otimização procura

determinar um d tal que 0 0( ) ( )f x f x+ <d em cada passo, ou seja, uma direção de

decréscimo em f, assim:

t t1

0 2min ( x ) min( H )f + = +d d g d d (2.10)

Escrevendo a condição de otimalidade de (2.10) ( d 0( ) 0f x∇ + =d ), obtém-

se:

-1= −d H g (2.11)

Assim, (2.11) fornece um mínimo global único para a função aproximadora

de f. A única desvantagem deste método é que os cálculos para a montagem da

matriz H solicitam um grande esforço computacional, sobretudo em problemas

com grande número de variáveis.

Os métodos Quase-Newton surgiram para resolver esse problema sem

perder as boas propriedades de convergência do método de Newton. Nesses

métodos, uma aproximação da Hessiana (ou de sua inversa) é construída a partir

dos valores dos gradientes ao longo das iterações. Esses métodos, dos quais o

BFGS (Broyden - Fletcher - Goldfarb - Shanno) é o mais popular, possuem

convergência superlinear e são amplamente utilizados em problemas de

otimização.

2.5 Busca Linear

A busca linear é um procedimento adotado tanto nos algoritmos sem

restrição como nos com restrição. Após a determinação da direção de busca d é

necessário calcular o tamanho do passo a ser dado nessa direção, a fim de se obter

o novo vetor das variáveis de projeto em (2.5). O tamanho do passo é calculado

fazendo-se uma minimização da função unidimensional p definida através da

expressão:

0( ) (x )p t f t= + d (2.12)

DBD



24

A partir desta definição, pode-se verificar que:

0(0) (x )p f= (2.13) e

t

0

(x) (x)(0)t

f dfpx dt =

′ =∂

(2.14)

onde p´ indica a derivada em relação à t.

A busca linear pode ser exata ou aproximada, dependendo do método

utilizado para a minimização. A busca aproximada é uma forma mais moderna, na

qual o objetivo é determinar t de forma que f apresente um certo nível de

decréscimo, segundo um critério preestabelecido, como:

0( ) ( ) + tp t f x tγ=≤ d g , (0,1)γ ∈ (2.15)

De acordo com esta equação, o parâmetro γ controla o tamanho do passo.

Assim, um γ pequeno permite a utilização de passos maiores e a utilização de um

γ grande força a utilização de passos pequenos.

Uma forma bastante popular de busca linear é fazer uma aproximação

quadrática de p e calcular t como o mínimo desta aproximação, verificando se a

equação (2.15) é satisfeita. Se isto não ocorrer, então a aproximação é atualizada

utilizando o novo ponto e o processo é repetido. Uma forma ainda mais simples é

o método de Armijo, no qual t é igual ao primeiro número da seqüência {1,α, α2,

α3, ...}, α ∈ (0; 1), para o qual p(t) satisfaz a condição (2.15).

2.6 Programação Quadrática

A Programação Quadrática (PQ) tem como objetivo determinar o vetor

solução *x do problema colocado na seguinte forma::

minimizar t t12x x x+q Q

sujeito a t =i ix ba 1...i l= t

i ix b≤a 1...i l m= + (2.16)

DBD



25

onde a é uma matriz que contem os coeficientes dos gradientes das restrições, b é

o vetor dos termos independentes das restrições.

Sendo Q uma matriz positiva definida, o problema quadrático é convexo e

pode-se garantir a existência de um único mínimo local.

A solução deste problema pode ser obtida em três etapas bem definidas

(Eboli, 1989 e Parente, 2000):

1. As l restrições de igualdade são eliminadas do problema diminuindo-se o

número das variáveis independentes para n - l, obtendo-se um problema de

programação quadrática (reduzida), chamado problema padrão de PQ, só com as

restrições de desigualdade.

2. O problema quadrático reduzido é transformado em um Problema Linear

Complementar (PLC), que pode ser resolvido através de métodos de pivoteamento

como o de Lemke.

3. Recupera-se a solução para o espaço original com o cálculo das variáveis

eliminadas na primeira etapa, obtendo-se os valores de x e λ .

2.7 Algoritmo de Han-Powell - Programação Quadrática Seqüencial

O algoritmo de otimização de Han-Powell proposto por Han em 1976 e

1977 e por Powell em 1978 (Eboli, 1989), foi implementado e aplicado a

problemas de Engenharia Estrutural no DEC/PUC-Rio por Eboli (1989), Parente

(2000) e Farfán (2000). Este algoritmo utiliza a técnica de Programação

Quadrática Seqüencial (PQS) através da resolução de um subproblema quadrático

(PQ).

O método de PQS pode ser considerado como o resultado da aplicação do

método de Newton à minimização de uma aproximação quadrática da função

Lagrangiana do problema. Este método fornece a cada iteração os vetores d

(correção de x) e ∆λ (correção dos multiplicadores de Lagrange λ ), os quais

atualizados são aproximadores da solução *x e. *λ Este fato pode ser

demonstrado considerando o problema:

minimizar ( )f x sujeito a ( ) = 0ic x (2.17)

DBD



26

cuja função Lagrangiana é dada por:

( , ) ( ) ( )i ii

L x f x c xλ= +∑λ (2.18)

Desenvolvendo ( , )L x∇ λ em séries de Taylor em torno de ( , )k kx λ até a

primeira ordem, obtém-se

1

1 1 21

( , ) ( , ) ( , )k

k k k k k k k kkL x L x L x+

+ ++

∇ + + ∆ = ∇ + ∇ ∆

dd λ λ λ λ

λ (2.19)

Considerando 1 1k k kx x+ += −d e 1 1k k k+ +∆ = −λ λ λ e aplicando a condição de

estacionariedade a (2.19) no ponto 1 1, )k k k kx + ++ + ∆( d λ λ , resulta:

1

21

( , ) ( , )k

k k k kkL x L x+

+

∇ = −∇ ∆

dλ λ

λ (2.20)

ou, expresso matricialmente, como

1

10

t k k k kk k

k kk

+

+

= −

∆

d g + AW AcA

λλ

(2.21)

Substituindo 1k+λ por 1k k++ ∆λ λ , tem-se:

1

10

t k kk k

k kk

+

+

= −

d gW AcA λ

(2.22)

onde, kA é a matriz dos gradientes das restrições, kW é a Hessiana da

Lagrangiana, e kg é o gradiente de f(x) sendo todos avaliados no ponto xk. A

solução de (2.22) equivale à solução do subproblema de PQ (Eboli,1989):

minimizar t12

tk k+g d d W d

sujeito a 0tk k+ =c A d (2.23)

DBD



27

Ou seja, cada iteração k da solução do problema original é aproximada pela

solução do PQ obtido pela linearização das restrições e pela expansão quadrática

de f em torno de 0x .

Em problemas em que todas as restrições são de igualdade, a direção de

busca e os multiplicadores de Lagrange podem ser obtidos pela solução do

sistema de equações lineares gerado pelo método de Newton aplicado a

Lagrangiana do problema, como mostrado em (2.22).

Para considerar o caso de restrições de desigualdade, pode-se resolver o

problema geral de PM da seguinte forma (Eboli, 1989):

minimizar ( )f x sujeito a ( ) = 0ic x 1...i l=

( ) 0ic x ≤ 1...i l m= + (2.24)

definindo uma direção de busca d e uma nova estimativa dos multiplicadores de

Lagrange λ através da solução do PQ:

Minimizar t12

tk k+g d d W d

sujeito a 0tk k

i ic + =a d 1...i l=

0tk k

i ic + ≤a d 1...i l m= + (2.25)

cujo método de solução foi visto na seção anterior.

2.7.1 Etapas do Algoritmo Não-Linear Han-Powell (PQS) As etapas que formam o algoritmo Han-Powell são (Parente, 2000):

1. Dado um ponto inicial 0x e uma aproximação da Hessiana da função

Lagrangiana 0B , fazer k = 0. 0B é dada pela seguinte função:

0 ob=B I (2.26)

onde ob é um parâmetro definido pelo usuário do algoritmo. O número

de reinícios da matriz B é controlado pelo parâmetro nr definido pelo

DBD



28

usuário. Segundo Parente (2000), o reinício de B serve para descartar a

influência de pontos muito distantes do ponto corrente.

2. Para k = k+1, montar e resolver o problema de programação quadrática

definido pela equação (2.25) determinando os vetores dk e λk:

Minimizar 1 t 112

tk k− −+g d d B d n∈ℜd sujeito a 1 1 0

tk ki ic − −+ =a d 1...i l=

1 1 0tk k

i ic − −+ ≤a d 1...i l m= + (2.27)

onde 1kic − é o vetor com as restrições, 1 tk

i−a é uma matriz com o gradiente

das restrições e 1k−B é uma aproximação da Hessiana no ponto 1kx − .

3. Verificar os critérios de convergência do algoritmo:

11

2max( )

tk k

ki

tol

c tol

− ≤

≤

g d (2.28)

onde o primeiro critério representa a variação da função objetivo na

direção dk e o segundo critério verifica explicitamente o valor da

restrição mais violada.

Verificar também os critérios de parada tais como: número de avaliações

da função objetivo e número de iterações.

4. Se os critérios de convergência e/ou os de parada não são atendidos, faz-

se então uma busca linear unidimensional para determinar o tamanho do

passo tk , na direção dk de forma que o novo estimador da solução 1k k k kx x t−= + d seja um ponto que contribua para o decréscimo da

função objetivo. A busca é feita sobre a função de penalidade (p),

construída no intuito de impor um alto custo à violação das restrições.

Esta função é definida pela expressão:

[ ]1 1

( ) ( ) ( ) ( ) max ( ),0l m

i i i ii i l

p t p x t f x r c x r c x= = +

= + = + +∑ ∑d (2.29)

DBD



29

onde os ir são os fatores de penalidades. A busca é aproximada, isto é a

solução *t não é o mínimo de ( )p t , mas atende a um certo decréscimo

pré-estipulado em ( )p t considerado satisfatório. O coeficiente de

decréscimo da função é dado pelo parâmetro γ definido pelo usuário.

5. Atualização da matriz Bk do subproblema quadrático através do método

BFGS.

6. Retorno à etapa 2.

2.8 Método dos Pontos Interiores

O algoritmo de Pontos Interiores (PI) foi implementado e aplicado a

problemas de Engenharia Estrutural no DEC/PUC-Rio por Parente (2000).

O algoritmo utilizado neste trabalho baseia-se na aplicação do método de

Newton para a solução do sistema de equações não-lineares obtidas a partir da

aplicação das condições de Kuhn-Tucker do problema de otimização (Herskovitz,

1995). Neste trabalho, apenas o algoritmo para restrições de desigualdade será

discutido, uma vez que os problemas de projeto ótimo a serem resolvidos não

possuem restrições de igualdade. No entanto, as mesmas idéias aqui apresentadas

também são válidas para os problemas que possuem simultaneamente restrições

de igualdade e de desigualdade e podem ser vistas em mais detalhes em

(Herskovitz, 1995; Herskovitz & Santos, 1997).

O método de Pontos Interiores tem como característica gerar uma seqüência

de pontos no interior da região viável que converge para a solução do problema.

Outra propriedade importante deste algoritmo é que cada um dos pontos

intermediários possui valores decrescentes da função objetivo, ou seja, se por

algum motivo a convergência não for alcançada o ponto final é sempre viável.

Considere o problema de otimização:

minimizar ( )f x sujeito a ( ) 0ic x ≤ 1...i m= (2.30)

cujas condições de Kuhn-Tucker são:

DBD



30

10

m

i iiλ

=

+ =∑g a

* *( ) 0i ic xλ = *( ) 0ic x ≤ * 0iλ ≥

(2.31)

Sendo A a matriz dos gradientes das restrições e C uma matriz diagonal

contendo os valores das restrições, as duas primeiras equações podem ser escritas

como:

0t+ =g A λ 0=Cλ

(2.32)

Aplicando o método de Newton para resolver o problema acima, obtém-se o

sistema:

0

0 0

t = −

d gW AA C λΛ

(2.33)

Na equação acima, Λ é uma matriz diagonal para a qual Λii = λi, d0 é a

direção de busca e λ0 é a estimativa dos multiplicadores de Lagrange. Pode-se

demonstrar que d0 é uma direção de decréscimo de f e que d0 = 0 se x for um

ponto estacionário (Parente, 2000).

A direção de busca fornecida por (2.33) nem sempre é uma direção viável.

Expandindo-se uma equação da parte inferior do sistema (2.33), chega-se a:

0 0 0i

ti i icλ λ+ =a d (2.34)

Esta equação implica que 0 0ti =a d para todo i tal que ci = 0.

Geometricamente, isto significa que d0 é tangente às restrições ativas, indicando

que a direção aponta para fora da região viável.

Uma solução para evitar este efeito é adicionar uma constante negativa do

lado direito da equação acima:

DBD



31

ti i i i icλ λ ρλ+ = −a d (2.35)

onde iλ é a nova estimativa de iλ .

Este procedimento faz com que a direção original seja defletida, de um valor

proporcional a ρ, para o interior da região viável. Como a deflexão é proporcional

a ρ e d0 é uma direção de decréscimo de f, é possível encontrar limites em ρ para

que d também seja uma direção de decréscimo. Este objetivo pode ser atingido

impondo-se que:

0t t

ak≤g d g d (2.36)

para ka ∈ (0; 1). Em geral, a taxa de decréscimo de f ao longo de d é menor que ao

longo de 0d . No entanto, este é o preço a ser pago para se obter uma direção de

decréscimo viável.

Considerando o sistema auxiliar:

1

1

t = −

d gW AA C λ λΛ

(2.37)

é fácil mostrar que:

0 1ρ+d = d d (2.38) e

0 1ρ+=λ λ λ (2.39)

Substituindo (2.38) em (2.36) chega-se a:

0

1

( -1)t

a tkρ ≤g dg d

(2.40)

Definida a direção de busca d, é necessário realizar uma busca linear restrita

ao longo dessa direção, de forma a garantir que o ponto gerado esteja no interior

da região viável. Além disso, é necessário atualizar os valores dos multiplicadores

de Lagrange de maneira a assegurar a convergência para a solução correta.

DBD



32

2.8.1 Etapas do Algoritmo de Pontos Interiores (PI) O algoritmo de Pontos Interiores para problemas de restrições de

desigualdade necessita de um ponto inicial viável 0x , uma estimativa para os

multiplicadores de Lagrange de forma que λi > 0 e uma matriz B simétrica e

positiva definida, que é uma aproximação de W. O algoritmo pode ser dividido

nos seguintes passos (Herskovits & Santos, 1997):

1. Obter a direção de busca d:

a) Determinar os vetores 0 0( , )d λ através da solução do sistema linear

definido em (2.33).

b) Verificar o critério de convergência:

tol≤d (2.41)

c) Determinar os vetores 1 1( , )d λ através da solução do sistema linear

definido em (2.37).

d) Calcular o valor de ρ:

2

1 0 0 1

21 0

0, , ( 1) /

0,

t t tf a

tf

se entao = min k k

se entao = k

ρ

ρ

> − ≤

g d d g d g d

g d d(2.42)

sendo 0fk > .

e) Calcular a direção de busca d:

0 1ρ+d = d d (2.43) e

0 1ρ+=λ λ λ (2.44)

2. Fazer uma busca linear sobre d, determinando o tamanho do passo t que

satisfaça um critério sobre o decréscimo da função objetivo e para o

qual:

DBD



33

( ) 0, s 0( ) ( ), s 0

i i

i i i

c x t e c x t c x e

λλ

+ ≤ ≥

+ ≤ <

dd

(2.45)

e o novo ponto x:

0x x t= + d (2.46)

3. Atualizar a matriz B, que é uma aproximação da Hessiana da função

Lagrangiana, através do método BFGS.

4. Definir uma nova estimativa para os multiplicadores de Lagrange:

2

0 0max ,ii ekλ λ = d (2.47)

sendo ke > 0.

5. Fazer x igual a 0x e retornar ao passo 1.

A aproximação inicial e o reinício da Hessiana da função Lagrangiana são

controlados pelos mesmos parâmetros utilizados pelo algoritmo de Programação

Quadrática Seqüencial.

2.9 Implementação

Os algoritmos foram implementados em linguagem FORTRAN-90 a partir

dos códigos utilizados por Parente (2000), sendo a maioria das rotinas transcritas

de C e outras modificadas de acordo com a adaptação do programa aos demais

módulos.

DBD


3 Análise Não-Linear Geométrica

3.1 Comentários Iniciais

Este capítulo começa com uma breve discussão sobre o comportamento não

linear, o objetivo da análise não linear, e o seu lugar na engenharia estrutural. As

fontes de não linearidade mais importantes no projeto de pórticos são listadas e a

formulação do problema será definida, onde são apresentados os referenciais

Lagrangianos. Estratégias para a solução numérica de equações não lineares são

mostradas, bem como as estratégias de incremento de carga, as estratégias de

iteração e os critérios de convergência. O capítulo termina mostrando soluções

clássicas existentes de alguns problemas elementares.

3.2 Comportamento Não-Linear, Análise e Projeto

O objetivo da análise estrutural é determinar o comportamento da estrutura

quando submetida a ações externas, ou seja, obter tensões, deformações e

deslocamentos. É importante notar que, a cada passo do processo de otimização, é

necessário se determinar a resposta da estrutura de forma a verificar se as

restrições de projeto são atendidas. Conseqüentemente, o processo de análise deve

ser o mais eficiente possível.

A maioria das estruturas de engenharia exibem um comportamento linear

elástico sob cargas de serviços. Existem exceções como arcos e edifícios altos, e

estruturas sujeitas a um escoamento localizado prematuro ou fissuração, por

exemplo, que apresentam um comportamento não-linear. Antes de alcançar o seu

limite de resistência, quase todas essas estruturas vão apresentar uma resposta

não-linear significante.

Na análise não-linear tenta-se melhorar a simulação do comportamento de

uma estrutura em alguns aspectos. O objetivo fundamental é se obter para fins de

projeto uma previsão segura do comportamento do sistema. Como conseqüência,

DBD


Análise Não-Linear Geométrica

35

tem-se um aumento da complexidade do problema e conseqüente aumento do

custo computacional.

3.2.1 Fontes de Não-Linearidade O comportamento não-linear de uma estrutura, sob ação de um

carregamento qualquer, pode ser classificado de acordo com seus efeitos. Dentre

as várias fontes de não linearidade, destacam-se:

Não-linearidade Física. Decorre do fato do material não apresentar uma relação

tensão-deformação linear (não segue a lei de Hooke), isto é, o comportamento do

material não é elástico linear. Os efeitos não lineares são descritos por formas

mais complexas de equações constitutivas (matrizes constitutivas não-lineares

e/ou equações constitutivas em termos de “taxas” ou “incrementos”). Pode-se ter

também não linearidade física nas relações momento-rotação de conexões semi-

rígidas ou flexíveis, ou de rótulas inelásticas oriundas de mecanismos de colapso

localizados (flambagem, plastificação ou fissuração localizadas em componentes

estruturais).

Não-Linearidade Geométrica. Uma estrutura pode ter um comportamento não-

linear, ainda que constituída de um material que obedeça à lei de Hooke. Para

valores relativamente grandes de deslocamentos, a deflexão lateral de um membro

pode trazer como conseqüência, o aparecimento de momentos fletores adicionais

(denominadas de segunda ordem), em virtude da presença de um esforço normal.

A esse tipo de comportamento não-linear, dá-se o nome de não-linearidade

geométrica. Neste caso os efeitos não lineares estão associados às equações de

equilíbrio, que consideram a configuração deformada, e as relações deformação-

deslocamento.

No presente trabalho será considerada somente a não-linearidade

geométrica.

DBD



36

3.3 Formulação para a Análise Não-Linear Geométrica de Estruturas Reticuladas

A formulação para a análise não-linear geométrica de estruturas tem seus

fundamentos teóricos na teria da elasticidade não-linear, que faz parte da

mecânica dos sólidos. A não-linearidade geométrica aparece, na teoria da

elasticidade tanto nas equações de equilíbrio, que são escritas utilizando-se as

configurações deformadas do corpo, quanto nas relações deformação-

deslocamento, que incluem termos não lineares nos deslocamentos e suas

derivadas.

3.3.1 Descrição do Problema Neste trabalho, um procedimento incremental-iterativo será utilizado para

traçar o caminho de equilíbrio da estrutura ao longo do tempo. O problema

fundamental é ilustrado na figura 3.1 onde se examina a resposta de um membro

típico de pórtico (figura 3.1a), o tramo ab da estrutura. Sua configuração inicial,

V V

(a)

H

∆

(b)

H

∆

tt+∆t

y

0

t+∆t

(c) x

ba

a

ta

a

0

t+∆t

b

tb

b

0p

tp

t+∆tp

0x'

x'

x'

t

t+∆t

0y'

y't

y't+∆t

Configuração 0

Configuração t

Configuração t+∆t

Figura 3.1. Deformações do elemento.

DBD



37

descarregada, indeformada (configuração 0), pode ser definida em termos do

sistema de coordenadas global fixo, x, y, ou em termos do sistema de coordenadas

locais, 0 0,x y′ ′ , no qual 0 x′ corresponde ao eixo da barra, no sentido a0, b0 (figura

3.1c). Considerando-se que, após a aplicação gradual do carregamento o sistema

muda da configuração 0 para a configuração t e que todas as variáveis do

problema já tenham sido determinadas nesta última configuração, estando o

sistema em equilíbrio. Pode-se agora tomar como referência o elemento ab na

configuração t ou no sistema de coordenadas globais ou no novo sistema local

atualizado, ,t tx y′ ′ , com t x′ determinado pelos extremos do elemento na nova

configuração. O elemento terá mudado a sua forma e dimensão neste processo,

mas as equações de equilíbrio formuladas desta forma estão satisfeitas e a posição

de qualquer ponto material, na posição inicial, p0, pode ser mapeada na nova

posição pt.

Tendo-se por base um estado de equilíbrio conhecido, em uma configuração

t, os procedimentos incrementais iterativos procuram determinar o próximo estado

de equilíbrio, em uma nova configuração +∆t t . As equações incrementais de

equilíbrio são obtidas a partir de aproximações lineares para os incrementos de

deslocamentos e deformações. Portanto, o equilíbrio em +∆t t não é satisfeito

exatamente e é necessário utilizar um procedimento iterativo em cada passo de

carga. Estes procedimentos serão estudados posteriormente. A relação correta

entre carga-deslocamento está indicada na figura. 3.1b.

Existem duas formas de descrição do movimento de um ponto material, p, a

descrição Lagrangiana e a Euleriana (McGuire, 2000). Para a análise de estruturas

a formulação Lagrangiana é mais natural, sendo aqui empregada. Na formulação

Lagrangiana usa-se as coordenadas de pontos materiais referidas à estrutura

indeformada (configuração 0) ou a uma estrutura de referência temporária

(configuração t). No referencial Lagrangiano Total (RLT), todas as variáveis

estáticas e cinemáticas no tempo +∆t t são referidas à configuração inicial

(indeformada) da estrutura, ou seja, o membro ab é referido a 0 0,x y′ ′ . Por outro

lado, no referencial Lagrangiano Atualizado (RLA), todas as variáveis estáticas e

cinemáticas são referidas à última configuração de equilíbrio da estrutura, ou seja,

o membro ab é referido a ,t tx y′ ′ .

DBD



38

Comumente, as formulações RLT e RLA têm sido usadas na análise

incremental não-linear de pórticos. Quando desenvolvidas consistentemente, as

duas formulações geram matrizes de rigidez global e vetor de forças idênticos. Ao

longo deste trabalho, o RLA será adotado por sua eficiência computacional em

resolver problemas do tipo viga (Bathe, 1996).

3.3.2 Principio dos Deslocamentos Virtuais O princípio dos deslocamentos virtuais para corpos deformáveis é dado por

int extW Wδ δ= , para o equilíbrio em +∆t t , nós temos

t

t t t t t tt ij t ij t

V

T dV Rδ ε+∆ +∆ +∆=∫ (3.1)

onde ijT é o tensor de tensões Piola-Kirchhoff II, ε ij é o tensor de deformações de

Green-Lagrange e t tR+∆ é o trabalho virtual das forças externas. Aqui, o

sobrescrito +∆t t refere-se à configuração final e o subscrito t à configuração de

referência.

As tensões +∆t tt ijT podem ser decompostas em

+∆ ∆ ∆= + = +t t t t t t

t ij t ij t ij ij t ijT T T T T (3.2)

onde tt ijT , tensões Piola-Kirchhoff II na configuração t , são idênticas às tensões de

Cauchy, tijT , na mesma configuração, e ∆t

t ijT é o incremento das tensões de Piola-

Kirchhoff II entre [ , ]+ ∆t t t .

Para as deformações de Green-Lagrange e os incrementos de deslocamento

expressos em termos da configuração de referência, ε ε+∆t tt ij t ij= , tem-se,

1

, , , ,2 ( )ε = + +t ij t i j t j i t k i t k ju u u u (3.3)

Decompondo (3.3) em parcelas lineares e não-lineares, t ij t ij t ijeε η= + , tem-

se

DBD



39

1, ,2 ( )= +t ij t i j t j ie u u e 1

, ,2 ( )η =t ij t k i t k ju u (3.4)

Substituindo (3.2) e (3.4) em (3.1), obtém-se

δ ε δ δ η∆ +∆+ + =∫ ∫ ∫t t t

t t t t tt ij t ij t ij t ij t ij t ij t

V V V

T dV T e dV T dV R (3.5)

A solução para a equação (3.5) não pode ser obtida diretamente, uma vez

que ela é não-linear nos incrementos dos deslocamentos. Para casos onde o

incremento de deformação t ijε dentro de cada passo da análise incremental pode

ser considerado pequeno, as seguintes suposições podem ser feitas (Yang & Kuo

1994)

∆ =t

t ij t ijrs t rsT C e e t ij t ijeδ ε δ= (3.6)

onde t ijrsC é o tensor da relação incremental de tensão-deformação. O resultado é

a equação linearizada

δ δ δ η +∆+ + =∫ ∫ ∫t t t

t t t tt ijrs t rs t ij t ij t ij t ij t ij t

V V V

C e e dV T e dV T dV R (3.7)

onde a segunda integral do lado esquerdo é igual ao trabalho virtual das forças

externas atuando no elemento na configuração t de equilíbrio

t

t tij t ij t

V

R T e dVδ= ∫ (3.8)

Então (3.7) se torna

t t

t t t tt ijrs t rs t ij t ij t ij t

V V

C e e dV T dV R Rδ δ η +∆+ = −∫ ∫ (3.9)

Na literatura, esta equação ou sua forma equivalente em termos de energia

são geralmente adotadas como base para o desenvolvimento de equações de

DBD



40

elementos finitos para análise não-linear de elementos do tipo barra (treliça,

pórtico), placas e sólidos.

Com a discretização do campo de deslocamentos pelo método dos

elementos finitos, os termos da equação (3.9) podem ser representados por

produtos de matrizes e vetores como se segue:

{ } [ ]{ }T

t ijrs t rs t ij tV

C e e dVδ δ= ∆∫ eu k u (3.10)

e

{ } { }Ttij t ij t

V

T dVδ η δ = ∆ ∫ gu k u (3.11)

onde [ ]ek é a matriz de rigidez elástica, gk é a matriz de rigidez geométrica,

u∆ são os incrementos de deslocamento gerados durante o intervalo [ ],t t t+ ∆ .

Assumindo que somente cargas concentradas nodais são aplicadas nos

elementos, tem-se que:

{ } { }Tt tR δ= u f (3.12) e

{ } { }Tt t t tR δ+∆ +∆= u f (3.13)

onde { }t f são as forças atuando no elemento em t e { }t t+∆ f são as forças atuando

no elemento em t t+∆ .

Substituindo (3.10)-(3.13) em (3.9) e tomando um campo de deslocamento

arbitrário { }uδ , tem-se

{ } { } { }t t t+∆ + ∆ + = e gk k u f f (3.14)

que representa o incremento de forças no intervalo [ ],t t t+ ∆ . Uma interpretação

física de (3.14) pode ser dada como se segue : o incremento de forças

{ } { }t t t+∆ −f f será resistido não somente palas ações elásticas geradas por [ ]ek ,

DBD



41

mas também pelas forças devidas à mudança da geometria representadas pela

matriz gk .

3.3.3 Aplicação a Elementos de Pórtico Plano Em elementos de pórtico plano existem dois componentes de tensões

independentes associados com dois componentes de deformação. Sejam as

tensões de Cauchy tτxx e tτxy e o vetor de deformações incrementais de Green de

t xxε e t xyε , onde as deformações são decompostas em parcelas lineares e não-

lineares da seguinte forma:

t xx t xx t xxeε η= + e t xy t xy t xyeε η= + (3.15)

Substituindo (3.15) em (3.9) tem-se

( ) ( )4 2t t t t txx xx xy xy xx xx xy xy

V V

Ee e Ge e dV dV R Rδ δ τ δη τ δη +∆+ + + = −∫ ∫ (3.16)

onde os fatores 4 e 2, foram adicionados para se levar em conta a simetria das

tensões cisalhantes e deformações, e o subscrito direito t para V e o subscrito

esquerdo t para t ije , t ijη , t ijrsC foram omitidos para melhor visualização. Tem-se

então que todas as variáveis da equação estão referidas à configuração t.

As deformações em um ponto arbitrário da seção x podem ser relacionadas

com os aos deslocamentos longitudinais e transversais, xu e yu , por

,xx x xe u= e 1, ,2 ( )xy x y y xe u u= + (3.17)

2 21, ,2 ( )xx x x y xu uη = + e 1

, , , ,2 ( )xy x y x x y y y xu u u uη = + (3.18)

Com base nas hipóteses de Euler-Bernouli, onde as seções transversais

inicialmente planas de uma viga permanecem planas e normais ao eixo da viga

após a deformação, os deslocamentos xu e yu em um ponto arbitrário podem ser

referenciados aos deslocamentos u e v do eixo da viga da seguinte forma

DBD



42

xu u yv′= − e yu v= (3.19)

onde ( )´ denota a diferenciação em relação a x. Substituindo (3.19) em (3.17),

tem-se os termos lineares:

xxe u yv′ ′′= − e 0xye = (3.20)

e os não-lineares

2 2 2 21

2 ( 2 )xx u v y v yu vη ′ ′ ′′ ′ ′′= + + − e 12 ( )xy u v yv vη ′ ′ ′ ′′= − + (3.21)

As forças inicias são obtidas através de integrações das tensões de Cauchy

τxx e τxy, ou seja:

xx t t

A

P dAτ= ∫ ; t txy

A

Q dAτ= ∫ ; e t txx

A

M y dAτ= ∫ (3.22)

onde tP são as forças axiais, tQ o cisalhamento transversal, tM o momento fletor e

A é a área de seção transversal da viga.

Substituindo as equações (3.20) - (3.22) em (3.16), tem-se

( ) 2 2 2

0 0

1 ( ) ( )2

L Lt t z

zIEAu u EI v v dx P u v P v dxA

δ δ δ δ ′ ′ ′′ ′′ ′ ′ ′′+ + + + ∫ ∫

{ } { }0

( ) ( )L

Tt t t t tM u v Q u v dxδ δ δ +∆′ ′′ ′ ′ + − = − ∫ u f f

(3.23)

O elemento de viga-coluna adotado é o esquematizado na figura 3.2. Trata-

se de um segmento reto, limitado pelos nós 1 e 2, que se deforma no plano de

definição da estrutura. As funções que relacionam os deslocamentos em um ponto

qualquer do elemento com os deslocamentos nodais são

{ }{ }1u u= N e { }{ }3v v= N (3.24)

DBD



43

y

x

x',u

y',v

F ,ux2 2

F ,vy1 1

F ,vy2 2

M ,θz2 z2

M ,θz1 z1

1

2

L

β

F ,ux1 1

Figura 3.2. Elemento de pórtico.

onde N1 e N3 são funções de interpolação lineares dadas por

1 1 x xL L

= −

N (3.25)

2 3 2 3 2 3 2 3

3 2 3 2 2 3 2

3 2 2 3 21 x x x x x x x xxL L L L L L L L

= − + − + − − +

N (3.26)

e os vetores de deslocamentos { }u e { }v são definidos como

{ }1 2u u u= e { }1 1 2 2v v vθ θ= (3.27)

O primeiro termo de (3.23), de acordo com (3.10), representa a matriz de

rigidez elástica. Fazendo uso de (3.24) ele pode ser escrito como

{ } { }{ } { } { } { }{ } { }T T T T1 1

0 0

L L

zu EA dx u v EI dx vδ δ′ ′ ′′ ′′+∫ ∫ 3 3N N N N

(3.28)

Fazendo uso das funções de forma dadas em (3.25) e (3.26), integrando e

combinando-se os termos, a matriz elástica fica da seguinte forma

DBD



44

3 2 3 2

2

3 2

0 0 0 0

12 6 12 60

4 6 20

0 0

12 6

4Simétrica

z z z z

z z z

z z

z

EA EAL L

EI EI EI EIL L L L

EI EI EIL L L

EAL

EI EIL L

EIL

− − − =

−

ek

(3.29)

Os componentes do segundo termo de (3.23), relacionados com a

deformação longitudinal e a interação entre a força axial e flexão, podem ser

escritos da seguinte maneira

{ } { }{ } { } { } { }{ } { }T T T T1 1

0 0

L Lt tu P dx u v P dx vδ δ′ ′ ′ ′+∫ ∫ 3 33N N N N

{ } { }{ } { }T T

0

Lt zIv P dx v

Aδ ′′ ′′+ ∫ 3 33N N

(3.30)

De acordo com a figura 3.3 e considerando a força cisalhante constante ao

longo do elemento, tM e tQ podem ser obtidos por

1 21

t tt t M MM M x

L+

= − + e 1 2 - t t

t M MQL+

= (3.31)

P

PM 1

M 2t

t

t

t

Figura 3.3. Elemento de pórtico.

O último termo de (3.23) representa as fontes de interação entre a flexão e a

deformação longitudinal. Fazendo uso de (3.24) e (3.31) tem-se

DBD



45

{ } { }{ } { }T T1 21 1

0

L t tt M Mv M x dx u

Lδ

+ ′′ ′+ ∫ 3N N

{ } { }{ } { }T T1 21 1

0

L t tt M Mu M x dx v

Lδ

+ ′ ′′+ + ∫ 3N N

{ } { }{ } { }T T1 23 1

0

L t tM Mv dx uL

δ + ′ ′+ ∫ N N

{ } { }{ } { }T T1 21 3

0

L t tM Mu dx vL

δ + ′ ′+ ∫ N N

(3.32)

Fazendo uso das funções de forma dadas em (3.25) e (3.26), integrando e

combinando os termos das equações (3.30) e (3.32), a matriz geométrica fica da

seguinte forma

1 2

2 3

T

=

g g

g

g g

k kk

k k (3.33)

onde

1

3 21

0

12 665 10

42Simétrica15

tt

t tt tz z

ttz

MPL L

I P I PP PL AL AL

I PPLAL

−

= + + +

gk (3.34)

2

3 22

12 2

0

12 6605 10

6 610 10

tt

t tt tz z

t t tt tz z

MPL L

I P I PP PL AL AL

M I P I PP PL AL AL

− − = − − + − − − −

gk (3.35)

DBD



46

2

33

0

12 265 30

42Simétrica15

tt

t tt tz z

ttz

MPL L

I P I PP PLL AL AL

I PPLAL

= + − + +

gk (3.36)

3.3.4 Montagem das Equações da Estrutura No referencial Lagrangiano Atualizado, a última configuração calculada é

tomada como referência para descrever o movimento da estrutura no passo

incremental de t para t+∆t. Com base nesta formulação, foi desenvolvida a

equação incremental de equilíbrio, equação (3.14). Nos métodos utilizados neste

trabalho a matriz + e gk = k k é chamada de matriz de rigidez tangente.

Uma vez que a equação incremental de equilíbrio foi desenvolvida para

cada elemento, o próximo passo é transformar essa relação para a estrutura em

estudo. Desta forma a equação (3.14) fica da seguinte forma

t t t+∆∆K u = P - P (3.37)

onde P é o vetor das forças externas e K é a matriz de rigidez da estrutura, ou de

maneira equivalente

∆ ∆ ∆F = K u = P (3.38)

onde ∆P é o incremento das forças externas e ∆F é o incremento das forças

internas.

O vetor de forças internas e os deslocamentos totais da estrutura são obtidos

de maneira incremental, ou seja, a cada passo de carga o acréscimo nas forças

internas e deslocamentos devem ser calculados. Para isso consideram-se válidas as

seguintes relações

t t t+∆ = ∆F F + F (3.39)

t t t+∆ = ∆u u + u (3.40)

DBD



47

O problema estrutural não-linear a ser resolvido pode ser expresso da

seguinte forma:

0t t t t t t+∆ +∆ +∆= =R P - F (3.41)

onde R é o vetor de forças residuais e o sobrescrito à esquerda representa a

configuração considerada.

Pode-se ainda expressar P da seguinte forma

t t t tλ+∆ +∆

refP = P (3.42)

onde λ é o parâmetro de carga, refP é o vetor de forças nodais tomado como

referência.

O vetor de forças internas é obtido no sistema local e transformado para o

sistema global através de uma matriz de rotação T da seguinte forma

1

net t t t t t

i

+∆ +∆ +∆

=

=∑F T f (3.43)

com t t+∆ f dado pela equação (3.14).

Da mesma maneira, a matriz de rigidez global da estrutura é obtida a partir

das matrizes de rigidez do elemento, (3.29) e (3.33), e transformada para o sistema

global através de uma matriz de rotação T da seguinte forma.

1

neT

i=

= + ∑ e gK T k k T (3.44)

3.4 Estratégias de Solução para Problemas Não-Lineares

No esquema tradicional do método de Newton-Raphson, o parâmetro de

carga λ é mantido constante durante os ciclos iterativos, funcionando bem na parte

ascendente do caminho de equilíbrio (trecho OA), figura 3.4, mas falha ao

descrever esta curva após o primeiro ponto limite (ponta A), o que levaria a uma

incorreta avaliação da capacidade pois o equilíbrio será atingido no ponto C.

DBD



48

Para se traçar à curva carga-deslocamento completa (trecho OABC), com

possíveis passagens pelos pontos limites, é necessário que seja permitida a

variação de λ a cada iteração.

λ

A

O

B

C

u Figura 3.4. Curva carga-deslocamento.

Basicamente, as dificuldades para o traçado da curva completa se devem ao

mau condicionamento da matriz de rigidez tangente nos pontos limites, onde ela é

singular,e o algoritmo apresentará um erro de over-flow na fatorização da matriz.

Felizmente esse não é um problema muito sério, pois é praticamente impossível

chegar precisamente em um ponto crítico.

3.4.1 Análise Incremental-Iterativa Considerando um instante t+∆t, que representa as diferentes etapas de

aplicação do carregamento e as correspondentes configurações de equilíbrio da

estrutura, tem-se que o vetor de forças residuais ( ) , 0,1,...t t w w+∆ =R , computado

após a w-ésima iteração de Newton-Raphson é dado por:

( ) ( ) ( ) 0t t w t t w t t wλ+∆ +∆ +∆= =refR P - F (3.45)

e ( ) ( ) )t t w t w+∆ ∆F = F( F, u (3.46)

A fim de se obter o próximo ponto de equilíbrio (w+1), as estratégias

incrementais para o tratamento de efeitos não-lineares consideram que em torno

de uma configuração deformada t t+∆ u , o problema é localmente linear. Desta

forma é feita uma expansão em série de Taylor da equação (3.45), sendo esta

DBD



49

aproximada por termos lineares obtidos a partir do truncamento dos termos de

ordem superior da série:

( ) ( )

( 1) ( ) ( 1) ( 1)( ) ( ) 0

t t w t t wt t w t t w w w

w t t wδ δλλ

+∆ +∆+∆ + +∆ + +

+∆

∂ ∂≅ + + ≅

∂∆ ∂R RR R u

u (3.47)

onde

( ) ( )

( )( )( ) ( )

t t w t t wt t T w

w w

+∆ +∆+∆∂ ∂

= − = −∂∆ ∂∆

R F Ku u (3.48)

sendo ( )TK a matriz de rigidez tangente, e

( )

( )

t t w

t t wλ

+∆

+∆

∂=

∂R P (3.49)

Substituindo as equações (3.48) e (3.49) em (3.47) e reorganizando os

termos se obtém a equação de equilíbrio

( )( ) ( 1) ( 1) ( )t t T w w w t t wδ δλ+∆ + + +∆K u = P + R (3.50)

onde δλ e δu são as correções do parâmetro de carga e dos deslocamentos nodais

obtidas durante o processo iterativo. De (3.50) tem-se que os deslocamentos

nodais iterativos (δu) podem ser decompostos em duas parcelas:

( 1) ( 1) ( 1) ( 1)w w w wδ δ δλ δ+ + + += +g ru u u (3.51)

onde:

( )( ) ( 1) ( )t t T w w t t wδ+∆ + +∆=gK u R (3.52)

( )( ) ( 1)t t T w wδ+∆ + =rK u P (3.53)

A correção do parâmetro de carga, ( 1)wδλ + , única incógnita da equação

(3.51), é determinada seguindo uma das estratégias de iteração fornecidas a

DBD



50

seguir, onde será introduzida uma equação de restrição que deve ser respeitada a

cada iteração. Para se considerar esses métodos de iteração tem-se inicialmente

que se supor que, para w=0, ∆λ0 tenha um valor prescrito dado pelo usuário ou

calculado automaticamente como será visto na seção 3.4.4.

Após a seleção de ∆λ0, determina-se o incremento inicial dos deslocamentos

nodais ∆u0. As aproximações ∆λ0 e ∆u0 caracterizam a chamada solução

incremental predita.

O primeiro passo para a obtenção da solução incremental inicial tangente

( 0λ∆ e 0u∆ ) consiste na montagem, usando informações da última configuração

de equilíbrio da estrutura, da matriz de rigidez tangente K. Após a definição de K,

resolve-se o sistema de equações:

( )t T δ =tK u P (3.54)

para determinar os deslocamentos nodais tangente, δut.

Com a definição de ∆λ0 e δut tem-se:

tuu δλ∆=∆ 00 (3.55)

Nesse estágio o parâmetro de carga e os deslocamentos nodais totais são

atualizados, ou seja:

0t t tλ λ λ+∆ = + ∆ (3.56)

0t t t+∆ = + ∆u u u (3.57)

onde λt e ut caracterizam o ponto de equilíbrio obtido no último passo de carga.

A Figura 3.5 fornece um esquema de solução incremental-iterativa para

sistemas com um grau de liberdade, onde os parâmetros de carga e o

deslocamento são atualizados seguindo a restrição de comprimento de arco

cilíndrico (Crisfield, 1991).

DBD



51

δu1 δu2

∆u 0

∆u1

∆u2

ut

∆λ

∆λ0

2∆λ1

δλ2

tλ

∆l

λ

u

solução predita

restrição

δλ1

Figura 3.5. Solução incremental-iterativa: sistema com um grau de liberdade.

Com a obtenção da solução iterativa ( 1) ( 1)( )w weδλ δ+ +u , faz-se a atualização

das variáveis incrementais do problema:

( 1) ( ) ( 1)w w wλ λ δλ+ +∆ = ∆ + (3.58)

( 1) ( ) ( 1)w w wδ+ +∆ = ∆ +u u u (3.59)

Para o parâmetro de carga e os deslocamentos nodais totais tem-se que:

( 1) ( 1)t t w t wλ λ λ+∆ + += + ∆ (3.60)

( 1) ( 1)t t w t w+∆ + += + ∆u u u (3.61)

Os procedimentos descritos nessa seção são repetidos até que um dado

critério de convergência seja atendido (veja a seção 3.4.3).

DBD



52

3.4.2 Estratégias de Iteração A determinação do parâmetro de carga iterativo, δλ(w+1) é função de uma

dada estratégia de iteração, ou equação de restrição imposta ao problema, que tem

a função de otimizar a convergência do processo iterativo. A seguir são

apresentadas duas estratégias bastante eficientes que serão utilizadas nos capítulos

seguintes.

3.4.2.1 Carga Constante Essa estratégia de iteração caracteriza o método tradicional de controle de

carga constante, no qual o parâmetro de carga é mantido constante durante o ciclo

iterativo. Para esse caso, tem-se que a equação de restrição se reduz à expressão

trivial:

( 1) 0wδλ + = w ≥ 0 (3.62)

Dessa forma a Equação (3.51) é reduzida aos deslocamentos fornecidos pelo

método convencional de Newton-Raphson.

3.4.2.2 Comprimento de Arco Cilíndrico Segundo Crisfield (1981), a cada iteração, a seguinte equação de restrição

deve ser satisfeita:

2T l∆=∆∆ uu (3.63)

Substituindo (3.51) em (3.59) e o resultado desta operação em (3.63), chega-

se a uma equação quadrática em δλ, a saber:

A B Cδλ δλ2 0+ + = (3.64)

onde, os coeficientes A, B e C têm a seguinte forma:

DBD



53

( 1) ( 1)Tw wA δ δ+ += r ru u ( 1) ( ) ( 1)2 ( )

Tw w wB δ δ+ += ∆ +r gu u u (3.65) ( ) ( 1) ( ) ( 1) 2( ) ( )w w T w wC lδ δ+ += ∆ + ∆ + −∆g gu u u u

Com a resolução de (3.64), chega-se aos dois valores δλ1 e δλ2, de forma

que se deve escolher entre as soluções:

( ) ( 1) ( 1)

1w w wδ δλ δ+ +∆ = ∆ + +1 g ru u u u

( ) ( 1) ( 1)2

w w wδ δλ δ+ +∆ = ∆ + +2 g ru u u u (3.66)

aquela que mais se aproxima da solução incremental da iteração anterior, ∆u(w).

Essa escolha deve prevenir um possível retorno, o que faria a solução regredir ao

longo do caminho já calculado. Um procedimento utilizado, consiste em se achar

o menor ângulo entre ∆u(w+1) e ∆u(w). Isto equivale a achar o máximo co-seno do

ângulo:

( )T T( ) ( ) ( 1)T T( ) ( 1) ( ) ( 1)

1,2 1,22 2 2cosw w w

w w w w

l l l

δ δθ δλ+

+ +∆ ∆ +∆ ∆ ∆= = +

∆ ∆ ∆

gr

u u uu u u u

(3.67)

Como (3.64) é uma equação quadrática, ela poderá ter raízes imaginárias se 2( 4 )B AC− for menor que zero. Essa situação pode existir quando o incremento

inicial do parâmetro de carga for muito grande, ou se a estrutura exibir múltiplos

caminhos de equilíbrio em torno de um ponto.

3.4.3 Critérios de Convergência O processo iterativo descrito termina indicando uma nova posição de

equilíbrio para a estrutura em análise quando um dos dois, ou os dois critérios de

convergência apresentados abaixo forem atendidos:

(1) o primeiro critério de convergência é baseado em relações de forças e é

calculado no início da iteração corrente utilizando parâmetros da iteração anterior.

Ele é definido como segue:

DBD



54

( )

1 ( )

w

wζ ζ

λ= ≤

∆

R

P (3.68)

onde ( )wR é igual à norma euclidiana do vetor das forças desequilibradas, que é

calculada usando-se o parâmetro de carga e os deslocamentos nodais totais da

iteração anterior, ( )wλ∆ P é a norma euclidiana do vetor de incremento de

carregamento externo e ζ é um fator de tolerância fornecido pelo usuário do

programa como dado de entrada.

(2) o segundo critério de convergência obedece a relações de deslocamentos e é

sempre verificado no final da iteração corrente. Ele é definido por:

2 ( 1)w

δζ ζ

+= ≤

∆

uu

(3.69)

onde δu é a norma euclidiana dos deslocamentos iterativos (residuais), ( 1)w+∆u

é a norma Euclidiana dos deslocamentos incrementais, que são obtidos após a

correção do processo iterativo, e ζ segue a mesma definição do critério anterior.

(3) o terceiro critério de convergência consiste em obedecer a ambas as relações

(forças e deslocamentos) dadas em (3.68) e (3.69), assim este critério é verificado

se:

ζ≤ζζ≤ζ 21 e (3.70)

onde 21 e , ζζζ são definidos nos itens (1) e (2).

DBD



55

3.4.4 Incremento Automático de Carga A obtenção da solução incremental inicial tem como passo fundamental a

definição do parâmetro de carga inicial ∆λ0. A seleção automática do tamanho do

incremento desse parâmetro é importante, e deve refletir o grau de não-linearidade

corrente do sistema estrutural em estudo. Em outras palavras, uma estratégia

eficiente de incremento automático de carga deve satisfazer basicamente os

seguintes requisitos: (i) produzir grandes incrementos quando a resposta da

estrutura for aproximadamente linear; (ii) gerar pequenos incrementos quando a

resposta da estrutura for fortemente não-linear; (iii) ser capaz de escolher o sinal

correto para o incremento, introduzindo medidas capazes de detectar quando os

pontos limites são ultrapassados.

A seguir são apresentadas algumas estratégias de incremento de carga que

satisfazem esses requerimentos.

3.4.4.1 Incremento do Comprimento de Arco Como proposto por Crisfield (1991), o incremento do comprimento de arco

a ser adotado como parâmetro de controle no passo de carga corrente pode ser

definido como:

1/ 2

t dt

Nl lN

∆ = ∆ (3.71)

onde t l e l∆ ∆ representam os incrementos do comprimento de arco no passos de

carga anterior (valor conhecido) e no passo de carga corrente (incógnita),

respectivamente; Nd é o numero de iterações desejadas para o processo iterativo

corrente, especificado pelo usuário, e t N é o numero de iterações que foram

necessárias para convergir no passo de carga anterior.

Através da Equação (3.71) e da condição de restrição escrita para a solução

incremental inicial:

0 0 2T

l∆ ∆ = ∆u u (3.72)

DBD



56

chega-se facilmente, usando-se (3.55) e (3.72), à expressão do incremento inicial

do parâmetro de carga:

0

T

lλδ δ∆

∆ = ±t tu u (3.73)

O critério utilizado para escolher o sinal correto na expressão (3.73) é o

sugerido por Yang e Kuo (1994), baseando-se no sinal do parâmetro GSP, que

será apresentado na seção seguinte.

No programa desenvolvido nesse trabalho, o usuário deve especificar 01λ∆ como

dado de entrada, sendo este valor usado em seguida para calcular 01∆u através de

(3.55). Substituindo-se 0∆u na Equação (3.72), chega-se a 1l∆ . Para os passos de

carga seguintes, calcula-se automaticamente ∆l através de (3.71).

3.4.4.2 Incremento Baseado no Parâmetro GSP Uma estratégia baseada na introdução de um parâmetro de rigidez

generalizado foi adotada por Yang e Kuo (1994) para limitar o incremento inicial

do parâmetro de carga. O método de solução é denominado de estratégia de

controle de deslocamento generalizado. Yang e Kuo (1994) propuseram a seguinte

equação para avaliar o parâmetro de carga:

1 1

0 01

T

t T

δ δλ λδ δ

∆ = ± ∆ t t

t t

u uu u (3.74)

considerando-se o parâmetro de rigidez generalizado (GSP) do sistema como se

segue:

1 1 GSP

T

t T

δ δδ δ

= t t

t t

u uu u (3.75)

Pode-se, portanto, reescrever (3.74) da seguinte forma:

DBD



57

0 01 GSPλ λ∆ = ± ∆ (3.76)

Observa-se que o sinal do incremento inicial de carga pode ser positivo ou

negativo. A escolha do sinal correto é de suma importância na definição de

seqüências de soluções (u, λ) que permitam um avanço contínuo na resposta

carga-deslocamento. De acordo com Yang e Kuo (1994), o sinal do parâmetro de

rigidez corrente depende exclusivamente dos vetores tδ tu (passo de carga

anterior) e δ tu (passo de carga corrente), conforme (3.75).

´+´

´+´´+´

´+´

´-´

´-´

û = δurt

ê = δur

û

û = δur

û

û

û

û

û

ê

ê

ê

ê

ê

ê

λ

u

Figura 3.6. Variação do sinal do parâmetro de rigidez generalizado (GSP).

O GSP torna-se negativo somente para o passo imediatamente após o ponto

limite, para os demais passos este permanecerá sempre positivo, o que pode ser

visto na figura 3.6.

3.5 Determinação dos Pontos Críticos

Os pontos críticos são aqueles em que um caminho de equilíbrio atinge um

valor extremo ou aqueles onde diferentes caminhos de equilíbrio se encontram.

Na figura 3.7 podem ser observados três pontos críticos (A,B e C), onde os pontos

(B e C) são chamados de pontos limite e o ponto (A) ponto de bifurcação. No

presente trabalho serão considerados somente pontos limite. Caso a estrutura

DBD



58

apresente ponto de bifurcação, uma pequena imperfeição será considerada para se

eliminar a bifurcação.

Matematicamente, um ponto de equilíbrio é um ponto crítico se a matriz de

rigidez do modelo de elementos finitos for singular ( )det 0T =K . Além disso,

sabe-se que o equilíbrio é estável quando todos os autovalores são positivos o que

leva a ( )det 0T >K e torna-se instável quando o menor autovalor se torna negativo

e, portanto, ( )det 0T <K .

A

E

B

C

D

λ

u Figura 3.7. Pontos críticos de uma estrutura.

Conforme visto na seção acima, as estratégias de incremento automático de

carga têm como objetivo gerar pequenos incrementos quando a resposta da

estrutura for fortemente não-linear além de detectar pontos limites. A

determinação dos pontos críticos se faz de forma aproximada, através do sinal do

GSP, que controla os incrementos de carga, ou seja, a medida que o problema se

aproxima de um ponto limite, os incremento de carga dados pela equação (3.76)

se tornam muito pequenos, sendo tomado como ponto limite o valor

imediatamente após a inversão de sinal do GSP.

3.6 Exemplos

Nesta seção são apresentadas as soluções de alguns problemas estruturais

encontrados freqüentemente na literatura em função da acentuada não-linearidade

da relação carga-deslocamento. Pretende-se, assim, verificar a eficiência da

formulação de elementos finitos não-linear aqui apresentada, bem como

determinar a melhor malha a ser utilizada no processo de otimização. As

estruturas escolhidas apresentam grandes deslocamentos e pontos limite.

DBD



59

Detalhes adicionais e outros problemas podem ser encontrados nas

referências ao longo do texto.

3.6.1 Pórtico de Lee Neste exemplo é feita a análise do chamado Pórtico de Lee. Esta estrutura é

interessante porque apresenta instabilidade por ponto limite depois de sofrer

grandes deslocamentos e rotações. É usada freqüentemente por pesquisadores para

validar estratégias de solução não-linear (Galvão, 2000; Parente, 2000). A

geometria do pórtico é mostrada na figura 3.8. Os valores numéricos empregados

são L = 120, H = 120, Lp = 24 e P = 1. A seção transversal do pórtico é retangular

de dimensões b = 3 e h = 2 e o módulo de elasticidade do material é 720.

Lp P

uθw

L

H

b

hA

Figura 3.8. Pórtico de Lee.

A estrutura foi analisada utilizando-se 5 malhas distintas, sendo 3 delas

ilustradas na figura 3.9. As malhas são formadas por elementos de pórtico. As

malhas com 20 e 100 elementos, não mostradas na figura, são subdivisões dos

modelos de 10 elementos

40

24 48 48

8 elem.

24 32 3232

10 elem.

24

30

6 elem.

Figura 3.9. Pórtico de Lee – malhas utilizadas.

Adotou-se a estratégia de iteração comprimento de arco cilíndrico

juntamente com o método de Newton-Raphson, com incremento automático do

DBD



60

comprimento de arco controlando o valor inicial do parâmetro de carga, 0λ∆ . No

início do processo adotou-se: 01 0.01λ∆ = . Para controlar a convergência foi

adotado 310ζ −= .

Foram efetuadas análises para modelos estruturais com: 6, 8, 10 e 20

elementos finitos, mostradas na figura 3.10. A análise considerando um modelo

com 100 elementos foi feita para servir de referência e avaliar a precisão dos

resultados analisados.

0 20 40 60 80 100

-1

0

1

2

fato

r de

carg

a

deslocamento w

6 elementos 8 elementos 10 elementos 20 elementos 100 elementos

Figura 3.10. Pórtico de Lee – curvas de equilíbrio.

L1

L2

Figura 3.11. Pórtico de Lee – configuração deformada.

Utilizando os resultados apresentados na figura 3.10, pode-se verificar que a

estrutura apresenta dois pontos limite. Estes pontos, encontrados ao longo da

curva, foram detectados através da mudança de sinal do parâmetro GSP.

DBD



61

TABELA 3.1 Pórtico de Lee – cargas críticas

Elementos da malha L1 L1 – erro (%) L2 L2 - erro (%)

06 1.97095 6.21 -1.24505 32.23 08 1.90181 2.48 -1.09332 16.11 10 1.88497 1.58 -1.04141 10.60 20 1.86291 0.39 -0.96560 2.55

100 1.85570 - -0.94161 -

A configuração deformada correspondente a cada um dos dois pontos limite

(L1 e L2) é ilustrada na figura 3.11. Pode-se verificar que esses pontos ocorrem

depois da estrutura sofrer grandes deslocamentos e rotações. O valor do fator de

carga crítica para as cinco malhas é mostrado na tabela 3.1. Estes resultados

mostram que o fator de carga em L1, obtido para a malha com 20 elementos, é

praticamente o mesmo que para a malha de 100. Portanto, como o interesse na

otimização é o ramo ascendente, utilizando-se apenas 20 elementos ao longo das

barras, o comportamento do pórtico é modelado com a precisão desejada. Os

valores obtidos para L1 e L2 concordam com os valores encontrados na literatura:

1.855cλ = e 0.9298cλ = − (Parente, 2000), respectivamente.

3.6.2 Pórtico de Williams Neste exemplo é feita a análise não-linear do Pórtico de Williams. Este

problema possui resultados analíticos e experimentais apresentados por (Williams

,1964), sendo freqüentemente utilizados para validar modelos numéricos (Yang e

Kuo, 1994). Esta estrutura apresenta um caminho de equilíbrio acentuadamente

P

uθw

L

H

b

hA

Figura 3.12. Pórtico de Williams.

não linear com perda de estabilidade por ponto limite. A geometria do pórtico é

mostrada na figura 3.12. Os valores numéricos empregados são L = 12.943, H =

DBD



62

0.386, e P = 1. A seção transversal do pórtico é retangular de dimensões b = 0.753

e h = 0.243 e o módulo de elasticidade é 1.03 x107.

A estrutura foi analisada utilizando 5 malhas distintas formadas por 6, 8, 10,

20 e 100 elementos, sendo a subdivisão feita de forma uniforme para cada barra.

Adotou-se a estratégia de iteração com comprimento de arco cilíndrico

juntamente com o método de Newton-Raphson Padrão, com incremento

automático do comprimento de arco controlando o valor inicial do parâmetro de

carga, 0λ∆ . No início do processo adotou-se: 01 0.01λ∆ = . Para controlar a

convergência foi adotado 310ζ −= .

Os resultados das análises para as malhas selecionadas são apresentados na

figura 3.13. A análise considerando um modelo com 100 elementos foi feita para

servir de referência e avaliar a precisão dos resultados obtidos.

0,00 0,15 0,30 0,45 0,60 0,750

10

20

30

40

50

60

70

6 elementos 8 elementos 10 elementos 20 elementos 100 elementos

fato

r de

carg

a

deslocamento w

Figura 3.13. Pórtico de Williams – curvas de equilíbrio.

Analisando as curvas da figura 3.13, pode-se verificar que a estrutura

apresenta dois pontos limite. Estes pontos foram detectados através da mudança

de sinal do parâmetro GSP.

O valor do fator de carga crítica para as cinco malhas é mostrado na tabela

3.2. Estes resultados mostram que o fator de carga em L1 obtido para a malha com

20 elementos é praticamente o mesmo que o obtido com a malha de 100. Portanto,

DBD



63

TABELA 3.2 Pórtico de Williams – cargas críticas

Elementos da malha L1 L1 – erro

(%) L2 L2 – erro (%)

06 35.46690 4.41 30.40257 2.81 08 34.75060 2.60 30.89583 1.24 10 34.42747 1.65 31.06543 0.69 20 34.00402 0.40 31.23855 0.14

100 33.87000 - 31.28282 -

como o interesse na otimização é apenas ramo ascendente, utilizando-se apenas

20 elementos ao longo das barras, o comportamento do pórtico é modelado com a

precisão desejada.

DBD


4 Análise de Sensibilidade


Conforme visto no Capítulo 2, os algoritmos utilizados neste trabalho

necessitam das derivadas da função objetivo e das restrições em relação às

variáveis de projeto para determinar a direção de busca do processo de

otimização.

De forma geral, os gradientes da função objetivo e das restrições são

calculados a partir dos gradientes das respostas da estrutura determinadas na etapa

de análise. Dependendo do problema, as respostas de interesse podem ser

deslocamentos, tensões, freqüências naturais e cargas críticas. Outra grandeza de

interesse, muito utilizada como função objetivo, é o peso (volume) da estrutura.

A análise de sensibilidade, também chamada de gradientes das respostas da

estrutura, desempenha um papel central no processo de otimização pois é avaliada

a cada passo do algoritmo.

Neste capítulo, são desenvolvidas analiticamente as sensibilidades dos

deslocamentos e a sensibilidade da carga limite para o elemento de pórtico com

comportamento geometricamente não-linear, bem como uma aproximação

utilizada neste trabalho a fim de se adaptar as características da análise não-linear

vista no capítulo anterior. É apresentado também o método das diferenças finitas a

fim de se validar as formulações analíticas. Por fim são analisados alguns

exemplos para validar as implementações.

4.2 Método Analítico

O método analítico consiste na diferenciação direta das equações de

equilíbrio lineares e não-lineares. Para facilitar o desenvolvimento das equações

para o cálculo da sensibilidade, considera-se uma estrutura descrita por uma única

variável de projeto x.

DBD


Análise de Sensibilidade

65

4.2.1 Sensibilidade dos Deslocamentos A equação de equilíbrio na análise linear é dada por:

( ) ( ) ( )x x xK u = P (4.1)

derivando-se (4.1) em relação a variável de projeto x, obtém-se

( )d d dxdx dx dx

+K u Pu K = (4.2)

Reorganizando (4.2) temos

( )d d d xdx dx dx

−u P KK = u (4.3)

onde d dxu representa a sensibilidade dos deslocamentos com relação às

variáveis de projeto x. Para calcular d dxu , os seguintes passos devem ser

seguidos:

(i) resolver a equação (4.1) para u,

(ii) computar o lado direito da equação (4.3) (chamado de pseudo-forças),

(iii) resolver a equação (4.3) para d dxu usando a já decomposta matriz de

rigidez, realizando somente retrosubstituições com o vetor de pseudo-

forças.

Na análise não-linear, a equação de equilíbrio incremental é dada da

seguinte forma:

( ), ( ); ) ( ) ( ), ( ); ) 0t tx x x x x x x∆ ∆ = ∆ ∆ ∆ =R( F u P - F( F u (4.4)

onde P é suposto, no caso geral, dependente das variáveis de projeto.

Derivando-se (4.4) em relação a variável de projeto x

0t

t

d d ddx dx dx x∆ ∂∆ ∂∆ ∆ ∂∆

− − − =∂ ∂∆ ∂

P F F F u FF u (4.5)

DBD



66

Reorganizando os termos de (4.5) e fazendo uso da seguinte relação

( )t t T+∆ ∂∆=

∂∆FKu (4.6)

tem-se

( )t

t t Tt

d d ddx dx dx x

+∆ ∆ ∆ ∂∆ ∂∆= − −

∂ ∂u P F F FK

F (4.7)

onde d dx∆u representa a sensibilidade dos incrementos de deslocamentos com

relação às variáveis de projeto x. A equação (4.7) pode ser escrita usando a

seguinte notação

( )

( ) ( )

t tt t T

x x

d d d ddx dx dx dx

+∆+∆

∆ ≠∆ ∆ ≠∆

∆ ∆ ∆= = −

u u u u

u R P FK (4.8)

onde o lado direito de (4.7) é diferenciado de forma total desconsiderando-se a

dependência implícita de ∆u sobre x.

Observa-se agora que, devido a formulação utilizar o RLA, tem-se um

problema dependente do caminho (path-dependent), ou seja, a sensibilidade das

forças internas td dxF computada na etapa anterior entra no lado direito da

equação (4.7), e sua atualização para o próximo passo é dada da seguinte forma:

( )

( )

t tt t T

x

d d ddx dx dx

+∆+∆

∆ ≠∆

∆ ∆= +

u u

F F uK (4.9)

A equação (4.8) é a equação básica para a solução da sensibilidade

incremental dos deslocamentos e algumas observações podem ser feitas ao seu

respeito:

• A matriz do lado esquerdo da equação (4.8) é a mesma matriz de rigidez

tangente obtida após a convergência para um ponto de equilíbrio.

• Na etapa de referência t, o lado direito da equação (4.8) é conhecido.

DBD



67

• A equação de sensibilidade em (4.8) é linear em d dx∆u . Assim nenhuma

iteração é necessária para resolvê-la.

• A sensibilidade é computada de forma incremental no final do passo, uma

vez conhecida a sensibilidade de td dxF no início do passo.

Conseqüentemente as sensibilidades incrementais têm que ser computadas

no final do passo, atualizadas e armazenadas para serem usadas no

próximo passo.

4.2.2 Método Analítico Aproximado (MAA) A equação (4.8), descrita acima, necessita de um grande esforço

computacional para o cálculo das sensibilidades incrementais. Em cada passo de

carga é necessário avaliar a sensibilidade dos deslocamentos incrementais, o que

envolve a resolução de um sistema de equações lineares para cada variável de

projeto x.

A determinação dos pontos críticos feita no presente trabalho (seção 3.5)

gera, com a diminuição do GSP, muitos incrementos de carga antes do ponto

crítico, o que tornaria o processo de otimização bastante lento se para cada ponto

de equilíbrio fossem atualizadas as sensibilidades incrementais. Para adaptar a

análise de sensibilidade ao presente trabalho, foram feitas algumas considerações

apresentadas a seguir.

A sensibilidade dos deslocamentos incrementais pode ser escrita da seguinte

forma:

t t td d d

dx dx dx

+∆∆=

u u u- (4.10)

Pré-multiplicando (4.10) por ( )t t T+∆ K e substituindo o resultado em (4.8),

tem-se

( ) ( )

( )

t t tt t T t t T

x

d d d ddx dx dx dx

+∆+∆ +∆

∆ ≠∆

∆ ∆= − +

u u

u P F uK K (4.11)

Tomando-se derivada de x da relação de equilíbrio no passo anterior dada

abaixo

DBD



68

( )( ( ); ) 0t t t

t t t Td d dx xdx dx dx x

∂ = − − = ∂P u FP -F u K (4.12)

e somando-se (4.12) a (4.11), tem-se a seguinte relação

( ) ( ) ( )

( )(

t t t t t tt t T t t T t T

x

d d d ddx dx x dx dx

+∆ +∆+∆ +∆

∆ ≠∆

∂ ∆= − − +

∂ u u

u P F F uK K - K ) (4.13)

onde

t t td d ddx dx dx

+∆ ∆= +

P P P (4.14)

Conforme visto na seção 3.44, as estratégias de incremento automático de

carga têm como objetivo gerar pequenos incrementos quando a resposta da

estrutura for fortemente não-linear e produzir grandes incrementos quando a

resposta da estrutura for aproximadamente linear, o que faz com que a variação da

rigidez da estrutura no intervalo [t, t +∆t] seja pequena, então

( ) ( ) 0t t T t T+∆ ≅K - K (4.15)

e o termo td dxu é eliminado de (4.13).

Finalmente, quando as forças externas P independem de x,

( )

( )

t t tt t T t t

x

d ddx x dx

+∆+∆ +∆

∆ ≠∆

∂ ∆= − − =

∂ u u

u F FK F (4.16)

onde t t+∆ F seria o vetor de pseudo-forças.

A equação (4.16) é uma equação aproximada para o cálculo da sensibilidade

total dos deslocamentos e algumas observações podem ser feitas a seu respeito:

• A matriz do lado esquerdo da equação (4.16) é a mesma matriz de rigidez

tangente obtida após a convergência para equilíbrio em t+∆t.

• Na etapa de referência t, o lado direito da equação (4.16) é conhecido.

DBD



69

• A equação de sensibilidade em (4.16) é linear em t td dx+∆ u . Assim

nenhuma iteração é necessária para resolvê-la.

• A sensibilidade t td dx+∆ u é computada a qualquer instante desejado, uma

vez conhecida a sensibilidade de t x∂ ∂F . Conseqüentemente, somente a

sensibilidade da força interna tem que ser computada no final do passo,

atualizada e armazenada para ser usada no próximo passo.

4.2.3 Sensibilidade da Carga Limite Para se calcular a sensibilidade da carga limite basta considerar o fator de

carga λ como dependente implicitamente de x, ou seja, o vetor de carregamento

externo fica definido como :

( ) ( ) ( ) ( ) ( )t t t t tx x x x xλ λ λ+∆ +∆= = + ∆ref ref refP P P P (4.17)

onde refP é um vetor de referência de magnitude arbitrária.

Derivando-se (4.17) em relação a x, tem-se:

t t t t

t t dd ddx dx dx

λ λ+∆ +∆

+∆= + refref

PP P (4.18)

Assim a equação (4.13) fica da seguinte forma

( )t t t t

t t T t t t tdd ddx dx dx

λ λ+∆ +∆

+∆ +∆ +∆= + −refref

PuK P F (4.19)

No nível da carga crítica ( *t tλ λ+∆ = ), um asterisco é adicionado nos termos

de (4.19) para identificar o ponto onde eles estão sendo avaliados. Então se

escreve

* *

( )* *T dd ddx dx dx

λ λ= + − *refref

PuK P F (4.20)

DBD



70

A matriz de rigidez tangente ( )*TK é singular e ( )* ( )* 0T Tϕ ϕ= =TK K , onde

ϕ é o autovetor associado com o autovalor nulo da matriz ( )*TK . Pré-

multiplicando a equação anterior por ϕ T , tem-se

* *

( )* *T dd ddx dx dx

λϕ ϕ ϕ λ ϕ= + −T T T T *refref

PuK P F (4.21)

Eliminando os termos nulos ( )* ( )* 0T Tϕ ϕ= =TK K e isolando o termo com *d dxλ , tem-se

*

* dddx dxλϕ ϕ ϕ λ= −T T * T ref

refPP F (4.22)

Finalmente, tem-se a expressão para o cálculo da sensibilidade do ponto

limite:

**

dd dxdx

ϕ λλ

ϕ

− =

T * ref

Tref

PF

P (4.23)

ou, quando as forças externas refP independem de x

*d

dxλ ϕ

ϕ=

T *

Tref

FP (4.24)

A equação (4.24) é uma equação aproximada para o cálculo da sensibilidade

da carga crítica e algumas observações podem ser feitas a seu respeito:

• A sensibilidade é calculada de maneira direta, uma vez que o vetor das

pseudo-forças é calculado de maneira incremental e conhecido no ponto

crítico.

• O autovetor associado com o autovalor nulo da matriz é calculado através

do método das inversas generalizadas para a matriz ( )*TK .

DBD



71

4.3 Método das Diferenças Finitas (MDF)

A mais simples técnica para cálculo da sensibilidade com respeito a variável

de projeto é a aproximação por diferenças finitas. Esta técnica é geralmente cara

computacionalmente, mas é de fácil implementação e por isto muito usada.

A mais simples aproximação por diferenças finitas é a de primeira ordem,

chamada de diferença a frente. Seja a função ( )f x , onde x é a variável de projeto.

A aproximação de primeira ordem, f x∆ ∆ , para a derivada, df dx , é dada por:

( ) ( )f f x x f xx x

∆ + ∆ −=

∆ ∆ (4.25)

onde x∆ é uma perturbação absoluta, pequena o suficiente para produzir

resultados satisfatórios. Geralmente essa perturbação é definida através da

seguinte expressão:

x xη∆ = (4.26)

onde η é o valor da perturbação relativa.

A maior dificuldade no MDF é selecionar o valor da perturbação η , se for

selecionada uma pequena perturbação, para reduzir o erro de truncamento (pois a

derivada só é exata quando x∆ tende a zero), pode-se ter um excessivo erro de

arredondamento causado pela forma como os números reais são representados nos

computadores. Perturbações relativas entre 10-4 a 10-8 geralmente levam a bons

resultados, o suficiente para aplicações práticas.

4.4 Exemplos

Nesta seção são apresentados exemplos de análise de sensibilidade de

estruturas geometricamente não-lineares. Os exemplos tratam tanto da

sensibilidade dos deslocamentos nodais como da sensibilidade das cargas críticas

das estruturas. Os objetivos são verificar a validade das expressões apresentadas

neste capítulo, testar a implementação computacional e comparar a precisão dos

métodos discutidos anteriormente.

DBD



72

As estruturas foram analisadas através do MEF utilizando a estratégia de

iteração comprimento de arco cilíndrico juntamente com o método de Newton-

Raphson Padrão, com incremento automático do comprimento de arco

controlando o valor inicial do parâmetro de carga, 0λ∆ . No início do processo

adotou-se: 01 0.01λ∆ = . Para controlar a convergência foi adotado 310ζ −= . A

malha menos refinada possui 20 elementos e a malha mais refinada 100 elementos

sendo as mesmas utilizadas na seção 3.6.

O MDF foi utilizado para verificar a consistência das sensibilidades

calculadas pelo MAA, utilizou-se η variando de 10-4 a 10-8, sendo considerado

como válido o que melhor representa a tendência dos resultados. Geralmente a

perturbação de 10-5 apresentou bons resultados.

4.4.1 Pórtico de Lee A geometria do pórtico é mostrada na figura 4.1. O comportamento desta

estrutura foi estudado no capítulo anterior. Os valores numéricos empregados são

L = 120, H = 120, Lp = 24 e P = 1. A seção transversal do pórtico é retangular de

dimensões b = 3 e h1 = h2 = 2 e o material utilizado apresenta E = 720. As

variáveis de projeto são as alturas das barras (h).

H

L

Lp P

h1

h2uθw


O deslocamento vertical (w) do ponto de aplicação da carga e a sua

sensibilidade calculada pelos MAA e MDF, para diversos fatores de carga, são

apresentados nas tabelas 4.1 e 4.2.

DBD



73

TABELA 4.1 Pórtico de Lee – deslocamentos e sensibilidades para h1

dw/dh1 λ n º de Elementos w MAA MDF

20 4.8014 -3.8981 -3.8949 0,6 100 4.8105 -3.9101 -3.9103 20 25.7312 -47.1949 -47.0603 1,5 100 25.8404 -47.1719 -47.1741 20 41.0602 -76.3053 -76.0002 1,8 100 41.3742 -77.7288 -77.7172

TABELA 4.2 Pórtico de Lee – deslocamentos e sensibilidades para h2

dw/dh2 λ n º de Elementos w MAA MDF

20 4.8014 -5.8162 -5.8122 0,6 100 4.8105 -5.8314 -5.8314 20 25.7312 -45.6404 -45.5511 1,5 100 25.8404 -45.8695 -45.8706 20 41.0602 -119.6585 -119.5125 1,8 100 41.3784 -126.3417 -126.3255

A carga crítica e a sua sensibilidade calculada pelos MAA e MDF são

mostradas na tabela 4.3.

TABELA 4.3 Pórtico de Lee – cargas críticas e sensibilidades

d crλ /dh1 d crλ /dh2 n º de elementos crλ

MAA MDF MAA MDF 20 1.86291 0.5665 0.5635 2.2323 2.2323 100 1.85570 0.5591 0.5623 2.2240 2.2280

4.4.2 Pórtico de Williams A geometria do pórtico é mostrada na figura 4.2. O comportamento desta

estrutura foi estudado no capítulo anterior. Os valores numéricos empregados são

P

uθw

L

Hh2h1


DBD



74

L = 12.943, H = 0.386, e P = 1. A seção transversal do pórtico é retangular de

dimensões b = 0.753 e h1 = h2 = 0.243 e o módulo de elasticidade é 1.03 x 107. As

variáveis de projeto são as alturas das barras (h).

TABELA 4.4 Pórtico de Williams – deslocamentos e sensibilidades para h1

dw/dh1 λ n º de elementos W

MAA MDF 20 3.0740 x 10-2 -0.1154 -0.1155 10 100 3.0760 x 10-2 -0.1156 -0.1156 20 9.8397 x 10-2 -0.5899 -0.5908 25 100 9.8694 x 10-2 -0.5956 -0.5954 20 18.1613 x 10-2 -3.4035 -3.4066 33 100 18.4330 x 10-2 -3.7314 -3.7355

O deslocamento vertical (w) do ponto de aplicação da carga e a sua

sensibilidade calculada pelos MAA e MDF, para diversos fatores de carga, são

apresentados na tabela 4.4. A carga crítica e a sua sensibilidade calculada pelos

MAA e MDF são mostradas na tabela 4.5.

TABELA 4.5 Pórtico de Williams – cargas críticas e sensibilidades

d crλ /dh1 n º de elementos crλ

MAA MDF 20 34.005 17.1285 17.1001

100 33.870 17.0632 17.0823

DBD


5 Otimização de Dimensões


O desejo de se obter o projeto ideal, considerando aspectos relacionados

com o consumo, desempenho ou eficiência, tais como quantidades mínimas de

peso, volume, massa, sempre foi um dos principais objetivos da engenharia

estrutural. Nos projetos atuais as estruturas são cada vez mais esbeltas e

complexas, no caso de pórticos as seções são reduzidas para se ganhar espaço

interno e diminuir o custo. As técnicas de otimização são ferramentas valiosas na

busca desse objetivo.

Um fator impulsionador no campo de otimização foi o desenvolvimento da

indústria aeroespacial onde o objetivo principal é se obter estruturas mais leves.

Em outras áreas, como a engenharia civil, mecânica e automotiva, o objetivo

principal é minimizar o custo da estrutura, embora o peso afete o custo e o

comportamento final do projeto.

A otimização de estruturas pode ser dividida em otimização de dimensões,

otimização de forma e otimização topológica. O objetivo da otimização topológica

é determinar a topologia ótima de uma estrutura através da eliminação de

elementos desnecessários e da criação de vazios. Na otimização de forma, busca-

se determinar a geometria ótima dos contornos externos e internos de estruturas

contínuas e das coordenadas nodais de estruturas reticuladas cujas dimensões e

topologia são fixas. A otimização de dimensões tem por objetivo determinar a

dimensão (seções transversais, espessuras, etc) de cada componente de uma

estrutura para a qual a forma e a topologia são fixas.

Durante o processo de otimização busca-se promover o aumento da

eficiência na utilização do material da estrutura. Uma otimização com base no

comportamento linear (físico e geométrico) gera, para um grande número de casos

práticos, uma estrutura com sérios problemas de instabilidade, pois a análise

fornece uma incorreta capacidade de carga da estrutura. Esses problemas podem

DBD


Otimização de Dimensões

76

ser evitados desde que o modelo de otimização seja formulado com o objetivo de

evitar problemas de instabilidade. Este é interesse principal deste trabalho.

5.2 Metodologia de Otimização

O problema de otimização a ser considerado neste trabalho consiste em

minimizar o peso da estrutura, considerando restrições sobre a carga crítica e

sobre os deslocamentos. Esta metodologia de projeto ótimo tem a vantagem de

tratar tanto o caso de estruturas que apresentam pontos críticos (ex: pilares) como

o de estruturas que apresentam grandes deslocamentos sem problemas de

estabilidade (ex: vigas).

Nesta metodologia, dois fatores de carga distintos devem ser considerados.

O primeiro corresponde à carga total aplicada (λa = 1), onde os deslocamentos

devem ser avaliados, e o segundo corresponde à carga crítica (λ*) da estrutura.

Assim, existe a possibilidade da ocorrência de duas situações, conforme ilustrado

na figura 5.1.

λλ*

λa

λ

λa

(a) λa < λ* (b) λa > λ*

u u

λ*

Figura 5.1. Situações encontradas na análise estrutural.

Na primeira situação, a carga aplicada é menor que a carga crítica (figura

5.1a). Neste caso, a análise pode ser completada sem problemas. Ao final, as

respostas da estrutura e os respectivos gradientes são avaliados.

A outra situação ocorre quando a carga aplicada é superior à carga limite da

estrutura (figura 5.1b). Neste caso, a análise não pode ser levada ao fim, pois a

carga aplicada não corresponde a nenhum ponto do caminho de equilíbrio.

Portanto, é impossível determinar os deslocamentos correspondentes ao nível de

DBD



77

carga desejado. Como conseqüência, é impossível avaliar as sensibilidades destas

grandezas, e o algoritmo de programação matemática é incapaz de prosseguir.

Diferentes procedimentos podem ser utilizados para contornar este problema

(Parente, 2000; Melo, 2000). A idéia é empregar um problema modificado, de

maneira a sempre utilizar os deslocamentos avaliados em um nível de carga

inferior à carga crítica. Neste trabalho optou-se por uma solução mais simples.

Sempre que houver risco da carga crítica ficar abaixo da carga aplicada, utiliza-se

o algoritmo de pontos interiores.

Conforme discutido no Capítulo 2, o algoritmo de pontos interiores gera

uma seqüência de pontos no interior da região viável com valores decrescentes da

função objetivo. Portanto, a carga crítica de todos os projetos gerados é

obrigatoriamente maior que a carga aplicada. A única desvantagem deste

algoritmo é a necessidade de gerar um projeto inicial viável. Contudo, na maioria

das vezes, este não é um problema sério, pois é simples criar uma estrutura com

resistência maior que a necessária.

5.3 Formulação do problema de Otimização

O problema de otimização de estruturas geometricamente não lineares,

submetidos a carregamento estático, pode ser formulado da seguinte forma:

minimizar ( )f x nx ∈ℜ sujeito a ( ) 0jc x ≤ 1...j m= l u

i i ix x x≤ ≤ 1...i nsecs= (5.1)

O modelo considera dois tipos de restrições: estruturais e geométricas. As

restrições estruturais visam garantir que o projeto atenda aos limites especificados

para os deslocamentos e a carga crítica. As restrições geométricas têm por

objetivo garantir que a geometria do modelo seja válida.

As variáveis de projeto serão as alturas das seções transversais, hi. Para cada

grupo de seções transversais, nsecs, tem-se uma variável de projeto. A altura de

um elemento qualquer é a medida da seção na direção paralela ao plano da

estrutura.

Para tornar a apresentação mais geral será mantido x como notação, ficando

assim válida para qualquer outra variável de projeto.

DBD



78

A função objetivo adotada será o volume da estrutura, no caso de estruturas

compostas por mais de um material, o volume deverá ser substituído diretamente

pelo peso da estrutura.

i i i1

ne

if V b h l

=

= = ∑ (5.2)

onde bi é a largura da seção, hi é a altura da seção, li é o comprimento do elemento

e o índice i se refere ao elemento.

As restrições sobre os deslocamentos nodais são calculadas com as ações

externas de serviço Q ( )1λ = e escritas como:

,lim= sign( ) 0dj j j jc u u u− ≤ , 1...j ndr= (5.3)

onde ,limju é o deslocamento máximo permitido e ndr é o número de

deslocamentos nodais restritos.

O valor do fator de carga crítico é controlado pela restrição:

*= - 0l

jc λ λ ≤ (5.4)

onde λ é o fator de carga limite e *λ é o fator de carga crítico.

Restrições laterais, na forma de limites diretamente impostos à variável de

projeto, são também incluídas:

l ui i ix x x≤ ≤ (5.5)

onde lix é o limite inferior e u

ix é o limite superior da variável de projeto ix .

5.3.1 Fatores de escala A diversidade de medidas presente no problema de otimização pode

acarretar diferenças significativas entre as suas magnitudes e causar problemas na

estabilidade numérica do algoritmo de solução. Desta forma a variável de projeto,

DBD



79

x, é definida como a razão entre o valor corrente do parâmetro e o correspondente

valor inicial, i. e.,

= ii 0

i

xxx

(5.6)

A derivada em relação a x é obtida na forma

(.) (.)0d d= xdx dx

(5.7)

A função objetivo e as restrições ficam definidas em termos da variável de

projeto adimensional, onde a função objetivo é agora a relação entre o volume

corrente da estrutura (V) e o volume inicial (V0):

0

VfV

= (5.8)

As restrições são reescritas de forma adimensionalizada, a saber

,lim

= sign( ) 1 0jdj j

j

uc u

u− ≤

(5.9)

*

= 1- 0lc λλ

≤ (5.10)

e as restrições laterais

l ui i ix x x≤ ≤ (5.11)

5.3.2 Cálculo dos gradientes Os algoritmos de programação matemática utilizados neste trabalho

precisam dos gradientes da função objetivo e das restrições para a determinação

DBD



80

da direção de busca. Assim, derivando-se a equação (5.8), obtém-se o gradiente da

função objetivo:

1

i 0 i

df dVdx V dx

= (5.12)

onde, conforme (5.7), tem-se

0

i i

dV dVxdx dx

= (5.13)

Para a variável de projeto x igual a h, a equação (5.12) toma a seguinte

forma:

1

1 ne0

i i iii 0

df b h ldh V =

= ∑ (5.14)

Verifica-se que esta expressão é uma função dos valores iniciais de h, desta

forma ela é constante durante o processo de otimização.

Partindo-se da equação (5.9), os gradientes das restrições sobre os

deslocamentos nodais podem ser escritos como:

,lim

sign( )=

d d dj j j j j j

i i i j i

dc c c du u dudx x u dx u dx

∂ ∂+ =

∂ ∂ (5.15)

onde a derivada do deslocamento resultante é dada pela expressão:

j j0

i i

du dux

dx dx= (5.16)

que é função das sensibilidades dos deslocamentos j idu dx .

Derivando-se a equação (5.10), tem-se o gradiente da restrição de

instabilidade, que pode ser escrito como:

DBD



81

*1= -l

i i

dc ddx dx

λλ

(5.17)

onde

* *

0

i i

d dxdx dxλ λ

= (5.18)

que é função da sensibilidade do ponto limite *id dxλ .

5.4 Implementação da Formulação

A formulação apresentada foi desenvolvida em linguagem FORTRAN 90, a

escolha da linguagem é justificada pelo fato do programa de análise não-linear já

estar implementado em FORTRAN 77/90. O programa de análise foi inicialmente

implementado por Silveira (1995) para problemas de contato e posteriormente

expandido por Rocha (2000), que implementou com sucesso algumas estratégias

de solução não-linear encontradas recentemente na literatura e Galvão (2000) que

implementou formulações de elementos finitos geometricamente não-lineares. O

programa foi adaptado ao presente trabalho para incorporar as rotinas de análise

de sensibilidade como será visto a seguir.

O sistema computacional desenvolvido neste trabalho é composto

basicamente por dois módulos: análise (estrutural e sensibilidade) e otimização. O

gerenciamento dos módulos se faz através de um programa principal, conforme

será visto abaixo.

5.4.1 Programa Principal O fluxograma do programa principal é mostrado na figura 5.2. Nele é feita a

leitura de dados (1) que define o tipo de análise a ser utilizada, linear ou não-

linear. Neste mesmo arquivo encontram-se as definições do modelo de elemento

finito adotado (número de pontos nodais, número de elementos, número de

materiais, número de seções, dimensão do problema, número de pontos nodais,

coordenadas nodais; incidência nodal, propriedades físicas e geométricas dos

elementos, condições de contorno, etc).

DBD



82

Alocar memória

Retornar

Início

Pontos Interiores Han-Powell

Define os parâmetros

do algoritmo

Algoritmo?

Define os parâmetros

do algoritmo

Leitura de dados (1)

Leitura das restrições

Imprimir saida da otimização

Alocar memória

Alocar memóriaDesalocar memória

Algoritmo de Pontos Interiores

Algoritmo de Pontos Interiores

Figura 5.2. Fluxograma do programa principal.

Após esta fase, o controle é passado para o algoritmo de programação

matemática. Quando necessário este algoritmo vai requisitar os valores da função

objetivo, das restrições e dos respectivos gradientes.

5.4.2 Processo de Otimização O problema de otimização descrito na seção 5.3 tem a forma padrão dos

problemas de programação matemática. Portanto, qualquer algoritmo de

otimização com restrições de desigualdade pode ser utilizado para sua solução.

Contudo, no caso de estruturas não-lineares com problemas de instabilidade, pode

ser mais conveniente utilizar o algoritmo de pontos interiores.

A forma geral dos algoritmos de programação matemática é ilustrada na

figura 5.3. Estes algoritmos são compostos por duas etapas principais: o cálculo

da direção de busca e a determinação do tamanho do passo (busca linear). Para

DBD



83

determinar a direção de busca, os algoritmos utilizam os valores correntes da

função objetivo e das restrições, bem como os gradientes destas grandezas em

relação às variáveis de projeto. Na etapa de busca linear, apenas os valores da

função objetivo e das restrições são utilizados.

Alocar memória

Retornar

Calcular a direção de busca

Definir xo

Fazer busca linear para obter o novo x

Convergiu ?

não

sim

Início

Figura 5.3. Forma geral dos algoritmos de PM.

A comunicação entre o algoritmo de otimização e o restante do sistema

computacional se dá através de duas funções, denominadas Análise/Sensibilidade

e Gradiente. A primeira é responsável pelo cálculo do valor da função objetivo,

das restrições e dos gradientes das respostas da estrutura, enquanto que a segunda

é responsável pelo cálculo dos gradientes da função objetivo e gradientes das

restrições.

Quando o algoritmo de programação matemática necessita do valor da

função objetivo e das restrições, a função Análise/Sensibilidade é chamada. As

tarefas realizadas por esta função são descritas na figura 5.4. Na análise estrutural

dispõe-se de outro arquivo de dados com as informações que controlam a solução

incremental-iterativa (número de incrementos de carga, número máximo de

iterações, incremento de carga inicial, tolerância para convergência, etc). Detalhes

adicionais sobre a implementação da análise estrutural podem ser encontrados em

Silveira (1995), Rocha (2000) e Galvão (2000). Finalmente, com base nos

resultados da análise, determinam-se os valores da função objetivo e das

restrições.

DBD



84

Conforme visto no capítulo 4, a análise de sensibilidade é feita em conjunto

com a análise estrutural, a cada passo de carga é atualizada a sensibilidade das

forças internas que será utilizada no nível de carga prescrito para cálculo da

sensibilidade dos deslocamentos e no ponto limite para cálculo da sensibilidade

do ponto limite.

Linear Não-linear

Montagem da matriz de rigidez K

Cálculo de ∆λ e ∆u0 0

Processo iterativo Newton-Rapshon

Cálculo de parâmetros do próximo incremento

Proc

esso

incr

emen

tal i

tera

tivo

A

Montagem do vetor de forças externas Q

Montagem da matriz de rigidez linear KL

Cálculo dos deslocamentos

nodais u

Montagem do vetor de pseudo-load p

Cálculo da sensibilidade dos deslocamentos

nodais du

Análise?

Montagem do vetor de cargas de referencia P

Retornar

λ (?)

λ≠[1,λcr]

λ =λcr

Cálculo da sensibilidade da carga crítica dλ

Cálculo da sensibilidade dos deslocamentos

nodais du

Leitura de dados (2)

λ =1

Início

Avaliar função objetivo e restrições

Figura 5.4. Fluxograma da função Análise/Sensibilidade.

A função Gradiente é mais simples, pois ela é chamada após uma chamada

à função Análise/Sensibilidade. Assim, na função Gradiente, é necessário apenas

avaliar o gradiente da função objetivo e o gradiente das restrições em relação às

variáveis de projeto, observando que o gradiente das restrições, equações (5.15)-

(5.18), são funções das sensibilidades da estrutura que já retornam da função

Análise/Sensibilidade.

DBD



85

5.5 Exemplos

Nesta seção serão apresentados exemplos de otimização de dimensões

utilizando análises lineares e não-lineares, empregando o sistema computacional

desenvolvido neste trabalho. As estruturas foram otimizadas utilizando-se os

algoritmos de Programação Quadrática Seqüencial (PQS) desenvolvido por

Wilson, Han e Powell e o algoritmo de Pontos Interiores (PI) desenvolvido por

Herskovits.

Os dois algoritmos foram apresentados no Capítulo 2, onde foram

apresentados os diversos parâmetros que influenciam o desempenho de cada

algoritmo. Entre os parâmetros numéricos do algoritmo de PQS estão a tolerância

para convergência (tol1) e a tolerância para violação de restrições (tol2). O

algoritmo possui ainda dois parâmetros adicionais, bz e gz. Os parâmetros

numéricos do algoritmo de PI são a tolerância para convergência (tol), os

coeficientes de deflexão da direção de busca (ka e kf) e o coeficiente para a

atualização dos multiplicadores de Lagrange (ke).

Os dois algoritmos utilizam os mesmos procedimentos para a atualização da

Hessiana da função Lagrangiana e para a busca linear. Os parâmetros de controle

destas etapas são o número de iterações para o reinício da matriz B (nr), o valor

inicial dos elementos desta matriz (b0), o coeficiente de decréscimo da função

objetivo (γ) e o valor base para a definição da seqüência de valores do tamanho do

passo (α).

Os valores usuais dos diversos parâmetros são mostrados na tabela 5.1. Nos

exemplos apresentados neste capítulo, alguns destes parâmetros foram variados de

maneira a melhorar o desempenho dos algoritmos. Quando valores diferentes dos

contidos na tabela 5.1 forem utilizados, eles serão explicitamente indicados.

TABELA 5.1 Valores usuais dos parâmetros.

rn 0b γ α 1tol 2tol zb zg tol ak ek fk

5 0.1 0.1 0.5 10-4 10-4 10-4 102 10-4 0.7 1.0 1.0

Os dois algoritmos requerem que os gradientes sejam calculados uma vez a

cada iteração, portanto o número de vezes que as sensibilidades são avaliadas é

igual ao número de iterações. Por outro lado, o número de análises é maior ou

igual ao número de iterações, pois durante a busca linear pode ser necessário

DBD



86

avaliar a função objetivo e as restrições várias vezes. Portanto, o desempenho de

cada algoritmo será avaliado de acordo com o número de iterações para

convergência (nit) e o número de chamadas à função de análise (nf).

5.5.1 Pórtico de Williams O primeiro exemplo a ser otimizado é o pórtico de Williams. O

comportamento desta estrutura foi estudado nos capítulos 3 e 4. O objetivo é

determinar as dimensões das barras de maneira a minimizar o volume, impondo-

se restrições sobre os deslocamentos e a carga crítica. A geometria, as condições

de contorno e o carregamento da estrutura são mostrados na figura 5.5. Os demais

dados do problema são L = 12.943, H = 0.386 e P = 10. A seção transversal do

pórtico é retangular com b = 0.753, e o módulo de elasticidade do material (E) é

1.03 x 107.

P

uθw

L

H

5 x hi


O modelo possui 10 variáveis de projeto, sendo cinco por barra, devido a

simetria esse número se reduz à metade. O valor inicial das alturas é 0.2430ih = ,

i= 1..5, divididas igualmente da esquerda para a direita. São utilizadas restrições

geométricas de maneira que a altura de cada grupo de elementos não seja inferior

a 0.05. O volume inicial do pórtico é igual a 4.7387.

TABELA 5.2 Pórtico de Williams – resumo dos resultados

caso Vol ,maxju *λ nit nf

Linear (PI, b0 =1.0) 1.19062 0.14998 0.1029 8 11 Não Linear (PI, b0 =1.0) 2.84012 0.08240 1.2000 11 21

O volume da estrutura foi minimizado considerando dois procedimentos de

análise, linear e não linear. Uma restrição estrutural foi imposta no nó central de

forma que o deslocamento vertical seja inferior a 0.15. A malha utilizada para a

DBD



87

análise estrutural é formada por 20 elementos de pórtico distribuídos igualmente

nas barras. O resumo dos resultados obtidos nos dois casos é mostrado na tabela

5.2.

Na minimização, supondo um comportamento linear, o volume obtido é

bastante inferior ao volume inicial, demonstrando o sucesso da otimização, mas a

estrutura apresenta instabilidade por ponto limite para um nível de carga muito

inferior à carga aplicada quando é feita a análise não-linear da mesma.

Com o objetivo de evitar os problemas encontrados na análise linear, a

estrutura foi otimizada utilizando uma análise não linear, incluindo-se a restrição

de estabilidade. Esta restrição impõe que a carga crítica seja pelo menos 20%

maior que a carga aplicada ( 1.2λ = ).

Após a otimização, verifica-se que o volume não diminuiu tão

sensivelmente quanto no caso anterior, mas a estrutura otimizada apresenta uma

carga crítica superior à carga aplicada, o que garante a sua estabilidade. Pode-se

notar também que o deslocamento obtido está bem abaixo do limite máximo, ou

seja, a restrição estrutural ativa é a restrição sobre a carga crítica.

TABELA 5.3 Pórtico de Williams – dimensões finais

caso 1h 2h 3h 4h 5h Linear 0.063 0.060 0.059 0.060 0.063 Não Linear 0.150 0.128 0.139 0.162 0.149

A descrição da situação final das dimensões pode ser vista na tabela 5.3.

Com base nestes resultados, verifica-se que a altura da seção transversal é

praticamente constante em todo o comprimento do pórtico na análise linear. Com

a inclusão da restrição de estabilidade e a consideração da não-linearidade

geométrica, há uma variação maior nas dimensões.

Para a análise não-linear utilizou-se a estratégia de iteração comprimento de

arco cilíndrico juntamente com o método de Newton-Raphson Padrão ( 310ζ −= ),

com incremento automático do comprimento de arco ( 01 0.01λ∆ = ).

Os segundo caso descrito foi otimizado através do algoritmo de PI, pois o

algoritmo PQS apresentou sérios problemas de convergência.

DBD



88

5.5.2 Pórtico de Lee A estrutura a ser otimizado agora é o pórtico de Lee. O comportamento

desta estrutura foi estudado nos capítulos 3 e 4. O objetivo é determinar as

dimensões das barras de maneira a minimizar o volume, impondo-se somente uma

restrição na carga crítica. A geometria do pórtico é mostrada na figura 5.6. Os

valores numéricos empregados são L = 120, H = 120, Lp = 24 e P = 1. A seção

transversal do pórtico é retangular com b = 3, e o módulo de elasticidade do

material (E) é 720.

L

H

LpP

a

cb d

5 x hi

5 x hi


A malha utilizada para a análise estrutural é formada por 20 elementos

distribuídos igualmente nas barras. O modelo possui 10 variáveis de projeto,

sendo cinco por barra. O valor inicial das alturas é 20ih = , i= 1..10, divididas

igualmente de baixo para cima no pilar e da esquerda para a direita na viga.

Como a estrutura apresenta grandes deslocamentos e problemas de

instabilidade, as respostas da estrutura são calculadas através de uma análise não-

linear. Para a análise não-linear utilizou-se a estratégia de iteração comprimento

de arco cilíndrico juntamente com o método de Newton-Raphson Padrão

( 310ζ −= ), com incremento automático do comprimento de arco ( 01 0.01λ∆ = ).

São utilizadas restrições geométricas de maneira que a altura de cada grupo

de elementos não seja inferior a 1. O volume inicial do pórtico é igual a 1440.

No capítulo 3 foi determinada a carga crítica desta estrutura, *λ =1.86291.

Pretende-se aqui minimizar o volume considerando uma carga crítica 50% e 100%

maior que a da estrutura inicial, ou seja λ =2.79 e 3.72, valores estes que tornam o

DBD



89

projeto inicial inviável, portanto, apenas o algoritmo de PQS será utilizado. A

descrição dos resultados se encontra na tabela 5.4.

TABELA 5.4 Pórtico de Lee – resumo dos resultados

λ Vol ,maxju *λ nit nf 2.79 1187.8 41.404 2.789 15 33 3.72 1280.4 37.645 3.721 14 28

Como não foram impostas restrições de deslocamento, o algoritmo gerou

uma estrutura com grandes deslocamentos e, conforme o desejado, com carga

crítica superior aos limites impostos, conforme visto na tabela 5.4. As descrições

finais das dimensões se encontram nas tabelas 5.5 e 5.6. Uma representação

esquemática das dimensões finais e do gráfico do momento fletor é mostrada na

figura 5.7.

TABELA 5.5 Pórtico de Lee – dimensões finais - pilar

λ 1h 2h 3h 4h 5h 2.79 1.381 1.531 1.000 2.073 2.124 3.72 1.512 1.751 1.000 2.299 2.317

TABELA 5.6 Pórtico de Lee – dimensões finais - viga

λ 6h 7h 8h 9h 10h 2.79 2.328 2.458 1.603 1.000 1.000 3.72 2.527 2.678 1.699 1.000 1.000

h2

h1

h3

h6

h5

h4

h7 h8 h9 h10

h2

h1

h3

h5

h4

h6 h7 h8h9 h10

Estrutura Inicial λ = 2.79h2

h1

h3

h5

h4

h6 h7 h8h9 h10

λ = 3.72

Figura 5.7. Alturas das seções transversais / momento fletor.

DBD



90

Neste exemplo, observa-se que, mesmo partindo de uma estrutura com carga

crítica bastante inferior àquela desejada, o algoritmo de otimização de dimensões

foi capaz de reduzir o volume.

5.5.3 Pórtico com três barras A geometria do pórtico é mostrada na figura 5.8. Uma carga distribuída

horizontal (h=κq) aplicada na mesma região que a carga vertical (q) foi

considerada para forçar a flambagem lateral do pórtico, eliminando assim a

flambagem simétrica. Os valores numéricos empregados são L = 1000, q = 100,

κ=0.001, e E = 21000000. A seção transversal do pórtico é retangular com a base

igual a 5.

A geometria da estrutura é modelada com 10 elementos por barra divididos

uniformemente. O modelo de otimização é formado por 20 alturas (h), sendo 10

para a viga e 5 para cada pilar, para garantir a simetria considerou-se as alturas

indicadas na figura 5.8. Assim, o número de variáveis de projeto caiu de 20 para

10. O valor inicial das alturas é 200ih = , i= 1..5 nos pilares, distribuídas

igualmente de baixo para cima, e 600ih = i= 6..10 nas vigas, distribuídas

igualmente da esquerda para a direita em ½ vão.

q

ac

ed

b

L

L

L/2

h

fg5 x hi

5 x hi

Figura 5.8. Pórtico com três barras.

São utilizadas restrições geométricas de maneira que a altura de cada grupo

de elementos não seja inferior a 10. O volume inicial do pórtico é igual a 500000.

Os deslocamentos, verticais e horizontais, nos nós a, b, c, d, e e são limitados ao

DBD



91

valor máximo de 7 em módulo. As coordenadas dos nós restritos estão na tabela

5.7.

TABELA 5.7 Pórtico com três barras – nós restritos

Nó A b c d e f g x 500 1000 1000 1000 1000 1000 1000 y 500 1000 900 800 700 600 500

Os resultados obtidos são apresentados nas tabelas 5.8 a 5.10:

TABELA 5.8 Pórtico com três barras – resumo dos resultados Caso Vol ,maxju nit nf Linear 220459.66 6.9936 18 25 Não Linear 264968.36 7.0000 16 36

TABELA 5.9 Pórtico com três barras – dimensões finais - pilares Caso 1h 2h 3h 4h 5h Linear 10.002 10.002 10.002 10.002 10.002 Não-Linear 10.000 12.208 14.942 15.861 15.085

TABELA 5.10 Pórtico com três barras – dimensões finais - vigas Caso 6h 7h 8h 9h 10h Linear 14.968 21.466 25.388 28.312 30.309 Não-Linear 21.690 24.877 26.606 27.579 28.025

Observando as dimensões finais do pórtico nota-se que o volume no caso

linear poderia ser menor ainda, pois todas as restrições geométricas dos pilares

estavam ativas. Verifica-se que a otimização fez com que a altura da viga seja

maior no centro, onde o momento fletor é maior.

A estrutura perfeita apresenta bifurcação associada a um modo de

flambagem lateral. Com a inclusão da imperfeição o caminho de equilíbrio se

afasta do caminho da estrutura perfeita e a mesma deixa de apresentar pontos

críticos. A inclusão da carga horizontal, κq, torna os deslocamentos horizontais

significativos mesmo para níveis baixos de carregamento.

Uma análise não-linear foi feita para as dimensões finais da otimização

considerando o comportamento linear, observa-se na figura 5.9 que o pórtico

DBD



92

apresenta grandes deslocamentos laterais para um nível de carga muito inferior a

aplicada. Isto se deve principalmente ao aumento da esbeltez dos pilares e da

diminuição da altura da viga. As restrições de deslocamento são violadas para um

nível de carga aproximadamente 36% menor que a aplicada.

0 10 20 30 400,0

0,2

0,4

0,6

0,8

1,0

1,2

1,4

fato

r de

carg

a

deslocamento

ua - Não Linear ua - Linear

Figura 5.9. Deslocamento horizontal do ponto a para caso linear e não-linear.

A curva do deslocamento horizontal no ponto a para as dimensões obtidas,

considerando o comportamento não linear é apresentada na figura 5.9. Verifica-se

que, a utilização da análise não linear provocou um aumento de 177% na

resistência da estrutura e um aumento de apenas 20 % no volume da estrutura. O

aumento da resistência da estrutura é obtido principalmente com o aumento da

largura dos pilares.

Neste exemplo utilizou-se, para a análise não-linear, a estratégia de iteração

comprimento de arco cilíndrico juntamente com o método de Newton-Raphson

Padrão ( 310ζ −= ), com incremento automático do comprimento de arco

( 01 0.01λ∆ = ).

Os dois algoritmos foram utilizados com os parâmetros usuais e forneceram

resultados semelhantes. Na tabela 5.8 são apresentados os resultados do PQS que

apresentou melhores desempenhos neste exemplo. O algoritmo de PI, para o caso

não-linear precisou de 44 iterações e 71 análises não-lineares.

DBD



93

5.5.4 Pórtico com 7 barras A estrutura a ser otimizada agora é um pórtico com 7 barras mostrada na

figura 5.10. Os valores numéricos utilizados são L1 = 1000, L2 = 250, L3 = 600,

P1 = 30000, P2 = 40000 e E = 21000000. A seção transversal do pórtico é

retangular com b = 20.

P2

L1

P2

P1

L1

L2

L3

a

c

e

db

h

g

f

h2h1 h3

h6h5h4 h7

Figura 5.10. Pórtico com 7 barras.

A geometria da estrutura é modelada com 22 elementos, sendo 2 por pilares

e 4 por viga. O modelo de otimização é formado por 7 alturas (hi) indicadas na

figura 5.10. O valor inicial das alturas é 400ih = , i= 1..3 (pilares) e 200

ih = i= 4..7

(vigas). São utilizadas restrições geométricas de maneira que a altura de cada

grupo de elementos não seja inferior a 5. O volume inicial do pórtico é igual a

3089242. Os deslocamentos, verticais e horizontais, nos nós b, c, d, f e g são

limitados ao valor máximo de 2 em módulo.

Novamente foram considerados dois procedimentos de análise, linear e não

linear. Com o objetivo de evitar os problemas encontrados na análise linear, a

estrutura foi otimizada utilizando uma análise não linear, incluindo-se a restrição

de estabilidade. Esta restrição impõe que a carga crítica seja pelo menos 20%

maior que a carga aplicada ( 1.2λ = ). No presente problema, foi adotado um

procedimento para abortar a análise caso λ seja maior do que 1.5.

O resumo dos resultados obtidos nos dois casos é mostrado na tabela 5.11.

TABELA 5.11 Pórtico com 7 barras – resumo dos resultados Caso Vol ,maxju *λ nit nf Linear 1723523.1 2.000 0.606 7 9 Não-Linear 2174406.1 2.000 1.848 16 47

DBD



94

A descrição da situação final das dimensões pode ser vista na tabela 5.12.

Com base nestes resultados, verifica-se que, para a análise, linear o volume

poderia ser menor ainda, pois algumas restrições geométricas estavam ativas. No

caso linear observa-se o aumento nas seções dos pilares das extremidades, já as

demais seções, tanto das vigas como o pilar central, foram reduzidas ao valor

mínimo. No caso não-linear houve também um aumento nas seções das

extremidades, mas as principais diferenças foram nas seções das vigas.

Verifica-se que, a utilização da análise não linear triplicou a carga crítica da

estrutura com apenas um aumento de 14 % no volume da estrutura.

TABELA 5.12 Pórtico com 7 barras – dimensões finais Caso 1h 2h 3h 4h 5h 6h 7h Linear 48.806 5.000 55.450 5.000 5.000 5.000 5.007 Não-Linear 47.065 16.985 64.088 12.168 5.160 5.448 8.185

Para a análise não-linear utilizou-se a estratégia de iteração comprimento de

arco cilíndrico juntamente com o método de Newton-Raphson Padrão ( 310ζ −= ),

com incremento automático do comprimento de arco ( 01 0.01λ∆ = ).

Os segundo caso descrito foi otimizado através do algoritmo de PI (b0 =0.1,

nr = 5) pois o algoritmo PQS, apresentou sérios problemas de convergência, já no

caso linear o PQS (b0 =0.01, nr = 10) teve melhor desempenho.

DBD


6 Comentários Finais

Neste trabalho foi desenvolvido um sistema computacional para a

otimização de dimensões de estruturas geometricamente não-lineares. Este

desenvolvimento exigiu a implementação e a integração dos módulos de

programação matemática, análise estrutural e análise de sensibilidade. Cada uma

destas etapas foi discutida nos capítulos anteriores e vários aspectos merecem ser

destacados:

As equações gerais para o cálculo das sensibilidades incrementais de

estruturas geometricamente não-lineares foram apresentadas. Em seguida, foram

desenvolvidas as expressões analíticas aproximadas necessárias ao cálculo das

sensibilidades em elementos de pórtico formulados com base no procedimento

Lagrangiano Atualizado. Utilizando estas expressões é possível calcular as

sensibilidades (deslocamento e carga crítica) aproximadas de uma dada malha de

elementos finitos em relação às variáveis de projeto de maneira eficiente,

armazenando-se apenas a sensibilidade do vetor de forças internas.

Com base nos exemplos analisados, pode-se dizer que o algoritmo de pontos

interiores é mais robusto, pois, partindo de um ponto interior, ele sempre converge

para a solução. Além disso, mesmo que uma solução ótima não seja obtida, o

projeto final sempre é viável.

O algoritmo de programação quadrática seqüencial, embora tenha

apresentado sérios problemas de convergência, foi mais eficiente que o algoritmo

de pontos interiores, além de poder partir de qualquer ponto.

Os métodos utilizados para a determinação de caminhos de equilíbrio

funcionaram de acordo com o esperado. Nenhum dos métodos é o mais indicado

para todas as situações, por isso é interessante dispor de diversas alternativas. No

processo de otimização foi utilizada a estratégia de iteração a carga constante,

DBD


Comentários Finais

96

uma vez que é necessário calcular as respostas da estrutura para um nível de carga

prescrito.

O método aproximado para a determinação dos pontos críticos apresentou

bons resultados. Mas o grande número de pontos gerados próximo ao ponto

crítico pode inviabilizar a sua utilização em estruturas de grande porte.

O sistema foi capaz de tratar estruturas lineares e não-lineares de uma

maneira uniforme. No caso de estruturas não-lineares, foram otimizadas com

sucesso estruturas apresentando grandes deslocamentos e pontos limite. Apesar

dos efeitos das imperfeições sobre a capacidade de carga da estrutura não terem

sido incluídos de maneira explícita na formulação, a otimização de estruturas com

imperfeições prescritas foi realizada com sucesso.

6.1 Sugestões

A seguir são apresentadas algumas sugestões para o desenvolvimento de

trabalhos futuros:

• Implementar um ambiente gráfico e iterativo para controlar todo o

processo de otimização, uma vez que a convergência nem sempre é atingida e,

mesmo que o processo convirja, os resultados nem sempre são satisfatórios,

podendo o usuário intervir sempre que julgar necessário.

• Incluir restrições sobre os esforços internos, além de outras restrições,

tendo em vista uma adequação a aspectos práticos.

• Estudar novos métodos de análise de sensibilidade como o método semi-

analítico e o método semi-analítico refinado, bem como aprofundar o estudo de

métodos analíticos para o Referencial Lagrangiano Atualizado.

• Implementar a sensibilidade da carga crítica linearizada, devido ao alto

custo computacional para a determinação dos pontos críticos e sua sensibilidade

através de uma análise não-linear.

DBD


Comentários Finais

97

• Estudar métodos exatos para a determinação dos pontos críticos (limite e

bifurcação) e para o cálculo das suas sensibilidades.

• Estudar novas estratégias de carga e iteração uma vez que, na

determinação do caminho de equilíbrio, nem sempre se obteve sucesso com as

formulações utilizadas.

• Utilizar novos algoritmos de otimização, incluindo pacotes comerciais.

Comparar o desempenho dos diversos algoritmos.

DBD


7 Referências Bibliográficas

BATHE, K. J. Finite Element Procedures. 3rd ed. Prentice Hall, 1996. CRISFIELD, M. A. Non-linear Finite Element Analysis of Solids and Structures, Chichester: John Wiley & Sons, Inc., 1991 EBOLI, C. R. Dimensionamento Ótimo de Seções de Concreto Armado à Flexão Composta Oblíqua, Dissertação de Mestrado, PUC, Rio de Janeiro, Brasil, 1989. FARFÁN, A. D. Aplicação da Análise Limite a Problemas Geotécnicos Modelados como Meios Contínuos Convencionais e Meios de Cosserat. Tese de Doutorado., PUC, Rio de Janeiro, Brasil, 2000. GALVÃO, A. S. Formulações Não-Lineares de Elementos Finitos para Análise de Sistemas Estruturais Metálicos Reticulados Planos. Dissertação de Mestrado., UFOP, Ouro Preto, Brasil, 2000. HAFTKA, R. T. & GÜRDAL, Z. Elements os Structural Optimization. 3rd ed. Dordrecht. Kluwer Academic Publishers, 1992. HERSKOVITS, J. A View on Nonlinear Optimization. Em J. Herskovits, editor, Advances in Structural Optimization, pp 71-117. Kluwer Academic Publishers, 1995. HERSKOVITS, J. & SANTOS, G. On the Computer Implementation of Feasible Direction Interior Point Algorithms for Nonlinear Optimization., Structural Optimization, 14, pp 165-172. 1997. KLEIBER, M. & et all. Parameter Sensitivity in Nonlinear Mechanics. Chichester: John Wiley & Sons, Inc., 1997 McGUIRE, W. & et all. Matrix Structural Analysis. 2nd ed. New York: John Wiley & Sons, Inc., 2000 MELO, A. M. C. Projeto Ótimo de Pórticos Planos de Concreto Armado. Tese de Doutorado, UFRJ/COPPE, Rio de Janeiro, Brasil, 2000. PARENTE JR, E. C. Análise de Sensibilidade e Otimização de Forma de Estruturas Geometricamente Não-Lineares. Tese de Doutorado., PUC, Rio de Janeiro, Brasil, 2000.

DBD


Referências Bibliográficas

99

ROCHA, G. Estratégias de Incremento de Carga e de Iteração para Análise Não-Linear de Sistemas Estruturais. Dissertação de Mestrado, UFOP, Ouro Preto, Brasil, 2000. SILVEIRA R, A. M. Análise de Elementos Estruturais Esbeltos com Restrições Unilaterais de Contato. Tese de Doutorado., PUC, Rio de Janeiro, Brasil, 1995. WILLIAMS, F.W. An Approach to the Nolinear Behavior of the Members of a Rigid Jointed Plane Framework with Finite Deflections, Quart. J. Mech. Appl. Math., 17(4), p. 451-469, 1964. YANG, Y. B. & KUO, S. B. Theory & Analysis of Nonlinear Framed Structures. Prentice Hall, 1994.

DBD


Documents

Anderson Pereira Projeto ótimo de pórticos planos com ... otimização de estruturas, programação matemática, análise não-linear e elementos finitos. ... tais como as que utilizam