( PDF ) Estudo de algoritmos de integração elemento por elemento para análise dinâmica não linear de estruturas

Download PDF

ads:

LUCIANA CORREIA LAURINDO MARTINS VIEIRA

ESTUDO DE ALGORITMOS DE

INTEGRAÇÃO ELEMENTO POR

ELEMENTO PARA ANÁLISE DINÂMICA

NÃO LINEAR DE ESTRUTURAS

Dissertação apresentada ao Programa de Pós-

Graduação em Engenharia Civil da Universidade

Federal de Alagoas como requisito parcial para

obtenção do título de Mestre em Engenharia Civil

MACEIÓ

Fevereiro/2004

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

LUCIANA CORREIA LAURINDO MARTINS VIEIRA

ESTUDO DE ALGORITMOS DE

INTEGRAÇÃO ELEMENTO POR

ELEMENTO PARA ANÁLISE DINÂMICA

NÃO LINEAR DE ESTRUTURAS

Dissertação apresentada ao Programa de Pós-

Graduação em Engenharia Civil da Universidade

Federal de Alagoas como requisito parcial para

obtenção do título de Mestre em Engenharia Civil

Área de concentração: Estruturas

Orientador: Prof. Dr. Adeildo Soares Ramos Júnior

MACEIÓ

Fevereiro/2004

ads:

iii

Vieira, Luciana Correia Laurindo Martins

Estudo de algoritmos de integração elemento por elemento para análise dinâmica

não linear de estruturas. Maceió, 2004. 101p.

Dissertação (Mestrado) – Universidade Federal de Alagoas. Programa de Pós-

Graduação em Engenharia Civil.

1. Análise dinâmica 2. Análise não linear 3. Algoritmos de integração 4.

Elemento por elemento. I. Universidade Federal de Alagoas. Centro de

Tecnologia. Programa de Pós-Graduação em Engenharia Civil. II-Título

A meu pai, Ricardo,

minha mãe, Nininha,

e meu irmão, Guilherme.

Agradecimentos

A minha família, por tudo. Amor, incentivo, paciência, dedicação, solidariedade,

vida.

Ao Prof. Dr. Adeildo Soares Ramos Júnior, pela orientação, confiança e amizade

desenvolvidas nesses anos.

Ao Prof. Dr. Eduardo Nobre Lages, pelas suas contribuições a este trabalho.

Aos Engenheiros Isaías Quaresma Masetti e Rogério Diniz Machado pelo apoio

e incentivo durante o período da graduação e do mestrado.

A todos os meus amigos, professores e colegas, pois com cada um de vocês

ganhei ensinamentos que permitiram esta conquista.

À Universidade Federal de Alagoas – UFAL, pela minha formação e à

Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – CAPES, pelo auxílio

financeiro.

Finalmente, a Deus, que nos permite viver para alcançar conquistas como esta.

vii

Sumário

Agradecimentos ..........................................................................................vi

Lista de Figuras...........................................................................................x

Lista de Tabelas........................................................................................xiii

Lista de Símbolos......................................................................................xiv

Lista de Abreviaturas ............................................................................xviii

Resumo ......................................................................................................xix

Abstract...................................................................................................... xx

Capítulo 1..................................................................................................... 1

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

1.1. RELEVÂNCIA DO TEMA ..................................................................... 1

1.2. OBJETIVOS........................................................................................ 3

1.3. ESTRUTURA DA DISSERTAÇÃO ..........................................................3

1.4. NOTAÇÃO UTILIZADA ....................................................................... 4

Capítulo 2..................................................................................................... 5

2. Introdução à Dinâmica de Estruturas................................................5

2.1. INTRODUÇÃO .................................................................................... 5

2.2. ABORDAGEM GERAL SOBRE DINÂMICA DE ESTRUTURAS...................5

2.3. DISCRETIZAÇÃO ESPACIAL................................................................ 6

2.3.1. Formulação em deslocamentos do Método dos Elementos Finitos.............. 8

2.3.2. Considerações sobre análises lineares e não lineares ............................... 11

2.4. DISCRETIZAÇÃO TEMPORAL............................................................12

2.4.1. Métodos Modais.......................................................................................... 12

2.4.2. Métodos Diretos.......................................................................................... 12

Capítulo 3................................................................................................... 14

3. Algoritmos de integração para dinâmica não linear ......................14

3.1. INTRODUÇÃO .................................................................................. 14

viii

3.2. ALGORITMOS DE INTEGRAÇÃO NO TEMPO.......................................14

3.3. ALGORITMOS EXPLÍCITOS ............................................................... 16

3.3.1. Algoritmo da diferença central................................................................... 17

3.4. ALGORITMOS IMPLÍCITOS ...............................................................19

3.4.1. Técnicas de solução de sistemas não lineares............................................ 20

3.4.2. Algoritmos da família Newmark ................................................................. 27

3.5. TÉCNICA ELEMENTO POR ELEMENTO ..............................................31

3.5.1. Formulação elemento por elemento em algoritmos explícitos ................... 33

3.5.2. Formulação elemento por elemento em algoritmos implícitos................... 33

3.6. MEDIDAS DE ERRO EM ANÁLISE DINÂMICA .....................................39

Capítulo 4................................................................................................... 43

4. Exemplos e Aplicações .......................................................................43

4.1. INTRODUÇÃO .................................................................................. 43

4.2. APLICAÇÃO LINEAR: SISTEMA MASSA-MOLA .................................. 44

4.3. VIGA ENGASTADA........................................................................... 51

4.4. CABO TRACIONADO ........................................................................57

4.4.1. Obtenção da configuração estática ............................................................ 58

4.4.2. Influência do refinamento da malha de elementos finitos .......................... 64

4.5. PÊNDULOS ......................................................................................66

4.5.1. Pêndulo Triplo ............................................................................................ 67

4.5.2. Pêndulo Quíntuplo ...................................................................................... 75

5. Conclusões...........................................................................................84

5.1. CONCLUSÕES, COMENTÁRIOS E SUGESTÕES....................................84

Referências bibliográficas ........................................................................87

Apêndice A................................................................................................. 89

A Exemplo numérico de multiplicação de matriz-vetor e vetor-

matriz-vetor elemento por elemento .......................................................89

A.1 DESCRIÇÃO DO EXEMPLO................................................................89

A.1.1 Multiplicação matriz-vetor elemento por elemento.................................... 91

A.1.2 Multiplicação vetor-matriz-vetor elemento por elemento .......................... 93

Apêndice B .................................................................................................95

B Formulação utilizada do elemento de barra....................................95

B.1 ELEMENTO DE BARRA ..................................................................... 95

Apêndice C................................................................................................. 99

C Solução estática para a viga engastada sujeita a momento

concentrado (exemplo 4.3)........................................................................ 99

C.1 DESLOCAMENTO VERTICAL DA EXTREMIDADE DIREITA ..................99

C.2 MOMENTO DE REFERÊNCIA...........................................................100

Lista de Figuras

Figura 1.1 – (a) Lançamento de estacas torpedo (âncoras para estrutura offshore); (b)

Estrutura flutuante submetida aos efeitos da natureza.................................. 1

Figura 2.1 – Idéia básica do Método dos Elementos Finitos............................................ 7

Figura 3.1 – Algoritmo do Método das Diferenças Centrais.......................................... 18

Figura 3.2 – Estratégia do Método dos Gradientes Conjugados..................................... 24

Figura 3.3 – Algoritmo do Método dos Gradientes Conjugados Não Linear................. 27

Figura 3.4 – Algoritmo de Newmark versão deslocamento. .......................................... 31

Figura 3.5 – Algoritmo do Método dos Gradientes Conjugados Linear ........................ 35

Figura 4.1 – Sistema massa-mola. .................................................................................. 44

Figura 4.2 – História de deslocamento do grau de liberdade 1 para ∆t = 0,7s . ............. 46

Figura 4.3 – História de deslocamento do grau de liberdade 1 para ∆t = 0,6s . ............. 46

Figura 4.4 – História de deslocamento do grau de liberdade 1 para ∆t = 0,4s . ............. 47

Figura 4.5 – História de deslocamento do grau de liberdade 1 para ∆t = 0,01s . ........... 47

Figura 4.6 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado. ......................................................................................... 48

Figura 4.7 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado. ......................................................................................... 49

Figura 4.8 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado. ............................................................. 50

Figura 4.9 – (a) Desenho esquemático e (b) discretização da viga analisada................. 51

Figura 4.10 – Curva de aplicação da força momento. .................................................... 52

Figura 4.11 – Máximas freqüências do sistema ao longo da análise.............................. 52

Figura 4.12 – Deslocamento vertical da extremidade direita da viga para

101t

−

⋅=∆ s.

..................................................................................................................... 53

Figura 4.13 – Deslocamento vertical da extremidade esquerda da viga para

(a)

105t

−

⋅=∆ e (b)

1045,1t

−

⋅=∆ s; Detalhe do deslocamento vertical da

extremidade esquerda da viga para (c)

105t

−

⋅=∆ e (d)

1045,1t

−

⋅=∆

s.. 54

Figura 4.14 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado. ......................................................................................... 55

Figura 4.15 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado. ......................................................................................... 56

Figura 4.16 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado. ............................................................. 57

Figura 4.17 – Desenho esquemático da malha de elementos finitos. ............................. 58

Figura 4.18 – Curva de aplicação do peso próprio. ........................................................ 59

Figura 4.19 – História de deslocamentos verticais do nó central. .................................. 60

Figura 4.20 – Detalhe dos deslocamentos verticais do nó central.................................. 60

Figura 4.21 – Deformada do cabo no tempo final de análise. ........................................ 61

Figura 4.22 – Detalhe da deformada no nó central......................................................... 62

Figura 4.23 – Etapas da análise. ..................................................................................... 62

Figura 4.24 – (a) Tempos de simulação; (b) Média de iterações por intervalo de tempo.

..................................................................................................................... 63

Figura 4.25 – (a) Deslocamentos verticais do nó central para o algoritmo NR – PCG; (b)

Detalhe dos deslocamentos verticais do nó central..................................... 65

Figura 4.26 – Desenho esquemático da malha de elementos finitos. ............................. 67

Figura 4.27 – História das máximas freqüências para o pêndulo triplo. ........................ 67

Figura 4.28 – História das mínimas freqüências para o pêndulo triplo. ......................... 68

Figura 4.29 – Etapas da análise do pêndulo triplo.......................................................... 69

Figura 4.30 – Deslocamento horizontal na extremidade direita para s01,0t =∆ . .......... 70

Figura 4.31 – Deslocamento vertical na extremidade direita para s01,0t =

∆

. .............. 70

Figura 4.32 – Comparação entre as respostas obtidas com o algoritmo NR - CG para

dois intervalos de tempo. ............................................................................ 71

Figura 4.33 – Detalhe da comparação entre as respostas obtidas com o algoritmo NR -

CG para dois intervalos de tempo............................................................... 71

Figura 4.34 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado. ......................................................................................... 72

Figura 4.35 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado. ......................................................................................... 73

Figura 4.36 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado. ............................................................. 74

Figura 4.37 – Etapas da análise do pêndulo quíntuplo. .................................................. 75

Figura 4.38 – (a) e (b) Deslocamentos horizontal e vertical na extremidade direita para

s01,0t =∆ ; (c) e (d) Deslocamentos horizontal e vertical na extremidade

direita para s001,0t =∆ . ............................................................................. 76

Figura 4.39 – (a) e (b) Deslocamentos horizontal e vertical na extremidade direita para

s0001,0t =∆ ................................................................................................ 77

Figura 4.40 – História no tempo da máximas freqüências. ............................................ 77

Figura 4.41 – História no tempo da mínimas freqüências. ............................................. 78

Figura 4.42 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado. ......................................................................................... 79

Figura 4.43 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado. ............................................................. 80

Figura 4.44 – (a), (c) e (e) Erros em energia para s01,0t

∆

,s001,0t

∆

e s0001,0t =∆ ;

(b), (d) e (f) Detalhes dos gráficos (a) (c) e (e)........................................... 81

Figura 4.45 – Deslocamento horizontal na extremidade direita para s01,0t =∆ , usando o

algoritmo NR - Método direto. ................................................................... 83

Figura 4.46 – Deslocamento vertical na extremidade direita para s01,0t =

∆

, usando o

algoritmo NR - Método direto. ................................................................... 83

Figura A. 1 – Estrutura de molas conectadas em série................................................... 89

Figura A. 2 – Numeração dos elementos........................................................................ 89

Figura B. 1 – Elemento de barra..................................................................................... 95

xii

Figura B. 2 – Geração do vetor de forças internas de matriz de rigidez do elemento de

barra. ........................................................................................................... 98

Figura C. 1 – Desenho esquemático do modelo. ............................................................ 99

Figura C. 2 – Configuração genérica............................................................................ 100

Figura C. 3 – Configuração de referência..................................................................... 101

xiii

Lista de Tabelas

Tabela 3.1 – Métodos da família Newmark.................................................................... 28

Tabela 4.1 – Nomenclatura simplificada para os algoritmos estudados......................... 44

Tabela 4.2 – Propriedades do sistema massa-mola......................................................... 44

Tabela 4.3 – Condições iniciais do sistema massa-mola................................................ 45

Tabela 4.4 – Valores dos intervalos de tempo................................................................ 45

Tabela 4.5 – Valores dos erros relativos em cada intervalo de tempo analisado. .......... 48

Tabela 4.6 – Valores dos tempos de simulação (s) para cada intervalo de tempo

analisado. .................................................................................................... 49

Tabela 4.7 – Valores das médias de iterações por intervalo de tempo para cada intervalo

de tempo analisado...................................................................................... 50

Tabela 4.8 – Propriedades físicas e geométricas da viga................................................ 51

Tabela 4.9 – Valores dos intervalos de tempo................................................................ 53

Tabela 4.10 – Valores dos erros relativos em cada intervalo de tempo analisado. ........ 55

Tabela 4.11 – Valores dos tempos de simulação (s) para cada intervalo de tempo

analisado. .................................................................................................... 56

Tabela 4.12 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado. ..................................................................... 57

Tabela 4.13 – Propriedades físicas e geométricas do cabo............................................. 58

Tabela 4.14 – Valores dos intervalos de tempo.............................................................. 59

Tabela 4.15 – Tempos de simulação e médias de iterações por intervalo de tempo. ..... 63

Tabela 4.16 – Freqüências máximas............................................................................... 64

Tabela 4.17 – Valores críticos de intervalo de tempo..................................................... 65

Tabela 4.18 – Análise geral do desempenho dos algoritmos iterativos.......................... 65

Tabela 4.19 – Máximo número de condição da matriz efetiva....................................... 66

Tabela 4.20 – Valores dos erros relativos em cada intervalo de tempo analisado. ........ 72

Tabela 4.21 – Valores dos tempos de análise (s) para cada intervalo de tempo analisado.

..................................................................................................................... 73

Tabela 4.22 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado. ..................................................................... 74

Tabela 4.23 – Comparação entre implementações para o algoritmo NR-CG. ............... 74

Tabela 4.24 – Valores dos tempos de análise para cada intervalo de tempo analisado.. 79

Tabela 4.25 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado. ..................................................................... 80

Tabela 4.26 – Valores das médias de iterações por intervalo de tempo para cada......... 82

Tabela A. 1 – Valores das rigidezes dos elementos........................................................ 90

xiv

Lista de Símbolos

A Área da seção transversal

Matriz que relaciona deslocamentos nodais e deformações

c Velocidade de propagação da onda

Matriz de amortecimento

Matriz de amortecimento tangente

energ

Erro para o critério em energia

força

Valor do erro para o critério de força

rel

Erro relativo

E Módulo de elasticidade longitudinal

)(f x

Função contínua

)(

Vetor de forças de volume

)(

Vetor de cargas concentradas

)(

Vetor de forças de superfície

Vetor de forças efetivo

ext

Vetor de forças externas

int

Vetor de forças internas

ref

Vetor de forças de referência para cálculo do erro usando o critério

de força

G Módulo de elasticidade transversal

)(xg

Gradiente de )(f x

Hessiana da função

)(f x

I Momento de inércia

J Momento de inércia polar

Matriz de rigidez

Matriz de rigidez efetiva na versão deslocamento da formulação

Matriz de rigidez efetiva na versão aceleração da formulação

Matriz de rigidez do elemento

Matriz de rigidez tangente

k Curvatura

L Comprimento do elemento

Matriz boleana de conectividade

Matriz de massa da estrutura

Momento de referência

glbn

Número de graus de liberdade do modelo

Matriz das funções de interpolação

Direção da busca linear

p Coordenadas não nulas do vetor

Vetor p sendo anuladas as posições dos graus de liberdade que não

estão conectados ao elemento

Matriz de pré-condicionamento

Vetor de forças desequilibradas ( ),(

intext

uuFF

−

)

Vetor de forças desequilibradas local do elemento

xvi

Vetor de forças desequilibradas do elemento, expandida para o

sistema global

t Valor do tempo corrente

Menor período do sistema de elementos finitos com

graus de

liberdade

Energia cinética

Vetor de deslocamentos nodais

Vetor de velocidades nodais

Vetor de acelerações nodais

Campo de deslocamentos no interior do elemento

Vetor de velocidades no interior do elemento

Vetor de acelerações no interior do elemento

Preditor de deslocamento do algoritmo de Newmark

Preditor de velocidade do algoritmo de Newmark

Vetor com as condições iniciais para deslocamento

Vetor com as condições iniciais para velocidades

ext

Trabalho das forças externas

int

Trabalho das forças internas

Escalar que minimiza )(f pu

Escalar que define a nova direção de busca / parâmetro do algoritmo

de Newmark

Parâmetro do algoritmo de Newmark

uδ

Vetor de deslocamentos virtuais

externo

Wδ

Trabalho virtual das forças externas

interno

Wδ

Trabalho virtual das forças internas

xvii

εδ

Vetor de deformações virtuais

t∆

Valor do intervalo de tempo

t∆

Valor do intervalo de tempo crítico para o MDC

u∆

Incremento de deslocamento

Campo de deformações

Densidade mássica do material

()

εεσ

Tensor de tensões

Variável

χ no tempo tn

∆

⋅

Variável χ no tempo t)1n(

∆

, na iteração k

max

Máxima freqüência da estrutura

xviii

Lista de Abreviaturas

CG Método dos Gradientes Conjugados

ExE Elemento por Elemento

MEF Método dos Elementos Finitos

MDC Método das Diferenças Centrais

MDF Método das Diferenças Finitas

NLCG Método dos Gradientes Conjugados Não Linear

NR Método de Newton-Raphson

nelem

Número de elementos do modelo de elementos finitos

PCG Método dos Gradientes Conjugados Pré-condicionados

PTV Princípio dos Trabalhos Virtuais

xix

Resumo

A dinâmica de estruturas compreende uma vasta área de atuação visto que as

ações da natureza são geralmente variáveis com o tempo. Análises de aeronaves,

estruturas

offshore, estruturas sujeitas a terremotos, entre outras, devem ser realizadas

considerando também os efeitos inerciais, que têm papel importante na performance e

segurança dessas estruturas. A análise do comportamento dinâmico dessas estruturas é

bastante complexa em virtude das diversas não linearidades que podem estar presentes.

Atualmente, essas análises são realizadas utilizando-se algoritmos numéricos

computacionais para discretização espacial e temporal, exigindo, em geral, grandes

esforços computacionais. Com o advento dos

clusters de computadores na década de

90, o emprego de algoritmos numéricos paralelos utilizando a técnica elemento por

elemento tem se difundido e se mostrado eficiente na solução desses problemas. Nesse

contexto, este trabalho está direcionado ao estudo de algoritmos de integração no tempo

para problemas de dinâmica não linear de estruturas pelo método dos elementos finitos,

que possibilitem seu uso futuro em arquiteturas computacionais paralelas. Os principais

algoritmos para dinâmica não linear de estruturas são revisados, assim como as

alternativas de implementação através da utilização da técnica elemento por elemento.

Em seguida, utilizando-se o programa computacional MATLAB, implementam-se

alguns desses algoritmos e avaliam-se suas características através de exemplos

ilustrativos. Observam-se, nesses exemplos, características como custo computacional,

número de iterações e qualidade das respostas geradas para problemas quase estáticos e

dinâmicos não lineares considerando-se diversos níveis de discretização espacial e

temporal.

Palavras-Chave: Análise dinâmica, Análise não linear, Algoritmos de integração,

Elemento por elemento

Abstract

The topic of structural dynamics involves an ample actuation field since the most

of nature actions are time dependent. The analyses of structures in aerospace

engineering, offshore engineering, civil engineering and others might be realized

including the inertial effects since they have important role in reliability and safety of

them. Due to the various nonlinearities that can be present in the dynamic analysis, the

preview of the dynamic behavior of a structure can become a complex task. That

analysis is generally executed using computational numeric algorithms improved for

spatial and temporal discretization where computational requirements are excessive.

Techniques based on element-by-element concept have gained recognition due to its

increasing efficiency caused by the recent development of parallel computers. This

work is focused on time integration algorithms for nonlinear dynamic problems

formulated by the finite element approach. Besides, the considered algorithms have

been prepared to be used, in the future, on parallel computational architectures. The

main algorithms for nonlinear structural dynamics and alternatives of implementation of

the element-by-element techniques are studied. The software MATLAB has been used

for implementation of some algorithms, whose features are analyzed by examples.

Some characteristics – as computational cost, number of iterations and solutions quality

– are observed for nonlinear quasistatic and dynamic problems, adopting several spatial

and temporal discretization levels.

Key-words: Dynamic analyses, Nonlinear analyses, Integration methods, Element-by-

element.

Capítulo 1

1. Introdução

1.1. Relevância do tema

Os modelos físicos e matemáticos utilizados na engenharia buscam sempre

representar, de forma mais fiel possível, os eventos que ocorrem na natureza. O

interesse pela dinâmica de estruturas sempre esteve presente no contexto da engenharia,

visto que na natureza as ações aplicadas às estruturas são geralmente variáveis com o

tempo. Quando possível, utilizam-se modelos independentes do tempo, ditos estáticos,

para representar algum dado fenômeno. Porém, em casos como aeronaves, estruturas

offshore, terremotos e muitos outros, as ações variantes com o tempo têm papel

principal na performance e segurança da estrutura. A Figura 1.1 ilustra alguns dos casos

acima citados.

(a) (b)

Figura 1.1 – (a) Lançamento de estacas torpedo (âncoras para estrutura offshore); (b)

Estrutura flutuante submetida aos efeitos da natureza.

Especialmente a partir da última metade do século passado, os recursos

computacionais passaram a integrar, de forma mais viável, o rol de ferramentas para a

solução de problemas de dinâmica estrutural, alavancando os desenvolvimentos nessa

área da engenharia. Mesmo passadas algumas décadas e considerando a incrível

velocidade com que a tecnologia da informação tem se desenvolvido, as soluções de

muitos problemas de dinâmica de estruturas são limitadas pelo custo computacional

exigido.

Devido à variedade de fenômenos incluídos no estudo de dinâmica de estruturas,

usa-se classificá-los em problemas quase estáticos, de dinâmica estrutural e de

propagação de onda. Essas classes são definidas de acordo com a participação maior ou

menor de altas e baixas freqüências na resposta da estrutura. Nesta dissertação o foco

dos estudos está direcionado para os problemas quase estáticos e de dinâmica estrutural,

problemas onde as freqüências de excitação são da ordem de grandeza das freqüências

naturais mais baixas da estrutura.

Para obtenção das respostas da estrutura, especialmente em problemas não

lineares, faz-se necessário o uso de técnicas computacionais para integração no tempo.

Várias opções estão disponíveis na literatura, possuindo, cada uma, particularidades que

as tornam mais adequadas a certos tipos de problemas. De um modo geral, esses

algoritmos são classificados em algoritmos implícitos e explícitos.

A formulação dos algoritmos explícitos é, em geral, mais adequada aos

problemas de propagação de onda, pois apresentam baixo custo computacional por

intervalo de tempo. Porém, essa formulação apresenta uma restrição para o intervalo de

tempo, baseada nas freqüências naturais mais altas da estrutura, penalizando o seu

desempenho em problemas regidos pelas baixas freqüências. Em geral, para que esses

algoritmos apresentem a performance desejada, exige-se a solução de um sistema

diagonal de equações (desacoplado), em cada passo de tempo. Para os problemas quase

estáticos e de dinâmica estrutural, os algoritmos ditos implícitos apresentam melhor

performance que os explícitos (COOK e outros, 1989). Essa constatação deve-se

especialmente ao fato de que nesses casos há maior interesse na integração das

freqüências mais baixas, permitindo, para os algoritmos implícitos, o uso de intervalos

de tempo maiores. Porém, tradicionalmente, o uso de algoritmos implícitos implica na

necessidade da montagem das matrizes globais do sistema estrutural, demandando mais

memória para armazenamento de dados durante a análise. Outra característica comum

nesses algoritmos é a construção da matriz de rigidez tangente da estrutura, que para

problemas não lineares podem apresentar dificuldades para geração, a exemplo dos

problemas de contato e dos que envolvem o atrito.

As situações que envolvem análises numéricas de grande porte podem ter seu

tempo de análise reduzido através de aplicações realizadas em paralelo. Em especial, ao

se tratar de algoritmos de integração, a utilização de técnicas elemento por elemento se

destaca para esse tipo aplicação. Essas técnicas tiram partido da independência dos

cálculos realizados no nível do elemento, o que permite a sua paralelização e economia

de memória.

Justifica-se assim a proposta de se estudar estratégias para alterações em

algoritmos de integração, embutindo neles algumas características particulares das

técnicas elemento por elemento, como a economia de memória e uma estrutura que

facilite a futura paralelização do código computacional.

1.2. Objetivos

Esta dissertação tem como objetivos estudar estratégias elemento por elemento

para algoritmos de integração no tempo; implementar algumas estratégias em

algoritmos não lineares de integração; e avaliação do comportamento desses algoritmos

diante de casos de análises dinâmicas e quase estáticas de estruturas.

Busca-se gerar conhecimento acerca do comportamento dos algoritmos

combinados com as estratégias elemento por elemento, visando o futuro uso destes em

ambientes de computação paralela.

1.3. Estrutura da dissertação

Esta dissertação está organizada em cinco capítulos mais dois apêndices. O

primeiro capítulo, de introdução, contém a relevância do tema abordado, onde justifica-

se o interesse pela análise dinâmica de estruturas, em especial pelo estudo de algoritmos

de integração. Ainda nesse capítulo, apresentam-se os objetivos da dissertação e uma

breve descrição da sua estrutura e conteúdo.

No segundo capítulo é introduzido o tema de dinâmica de estruturas. Apresenta-

se a discretização espacial, através da formulação em deslocamentos do Método dos

Elementos Finitos, e de forma sucinta aborda-se a discretização temporal, introduzindo-

se o tema principal do trabalho.

O terceiro capítulo possui todo o conteúdo explorado acerca de algoritmos de

integração. Inicia-se revisando as classes de algoritmos, apresentando-se os mais

tradicionais, e descrevem-se brevemente algumas técnicas de solução de sistemas não

lineares, usadas nos algoritmos estudados. Aborda-se a técnica elemento por elemento,

descrevendo suas características e sua formulação, e introduzindo-a nos algoritmos

antes apresentados. Finalizando esse capítulo, encontra-se uma descrição de algumas

medidas de erros utilizadas para critérios de convergência e quantificação da qualidade

das respostas obtidas pelos métodos implementados.

No quarto capítulo, apresentam-se exemplos onde são utilizadas as

implementações dos algoritmos estudados. São abordados casos lineares e não lineares,

e de dinâmica estrutural e quase estáticos. Avalia-se a qualidade das respostas geradas,

o tempo gasto pelos algoritmos, número de iterações, comentando-se os resultados

encontrados.

As conclusões e sugestões para trabalhos futuros estão descritas no quinto

capítulo deste trabalho. Esse capítulo é seguido pela lista de referências utilizadas no

decorrer do trabalho, organizadas em ordem alfabética.

Para complementar o conteúdo do trabalho, apresentam-se dois apêndices. O

apêndice A aborda as operações passíveis de serem realizadas elemento por elemento,

executando-as passo a passo para esclarecer em quais etapas a técnica é aplicada. Por

fim, o apêndice B descreve a geração das matrizes de massa e rigidez e do vetor de

forças internas para o elemento de barra bidimensional, usado nos exemplos explorados.

1.4. Notação utilizada

Nesta dissertação adota-se a notação matricial, utilizando-se o negrito para

representar vetores, tensores e matrizes. As letras maiúsculas simbolizam as matrizes e

as minúsculas os vetores e tensores.

Para representar o valor de uma variável em um determinado intervalo de tempo,

utilizam-se índices subscritos que indicam o intervalo de tempo corrente. No caso de

processos iterativos, o valor da variável em cada iteração é indicado por índices

sobrescritos contendo o valor da iteração corrente.

Capítulo 2

2. Introdução à Dinâmica de Estruturas

2.1. Introdução

Este é um capítulo introdutório onde se apresenta a dinâmica de estruturas como

o estudo da formulação e da solução de equações diferenciais parciais onde o tempo e o

espaço são variáveis.

Adota-se o processo de semidiscretização, realizando-se inicialmente a

discretização espacial através do Método dos Elementos Finitos. Descreve-se a

formulação em deslocamentos desse método, gerando um novo sistema de equações,

porém equações diferenciais ordinárias no tempo.

Ao final, de forma sucinta, trata-se do tema da discretização temporal,

descrevendo os principais tipos de métodos de integração no tempo, introduzindo o

tema principal do trabalho que é abordado no próximo capítulo.

2.2. Abordagem geral sobre dinâmica de estruturas

A formulação matemática utilizada em dinâmica de estruturas é baseada em

equações diferenciais parciais hiperbólicas, tendo como variáveis o tempo e o espaço. A

obtenção dessas equações diferenciais depende de uma modelagem física e matemática

do fenômeno de interesse. O modelo é, portanto, uma tentativa de se simular, através de

equações ou sistemas de equações, o fenômeno observado.

Existem procedimentos formalizados para sistematizar a obtenção das equações

diferenciais que regem o comportamento do sistema. Dois deles são bastante usados, o

equilíbrio direto de forças e os métodos energéticos, cujos procedimentos podem ser

encontrados na literatura (COOK e outros, 1989).

Nos itens 2.3 e 2.4 estão apresentados métodos para solução das equações

diferenciais que regem o fenômeno em estudo. Porém, a escolha dos métodos mais

apropriados depende da natureza do problema, incluindo os tipos das ações e as

propriedades do sistema estrutural. Dentro do contexto da análise dinâmica, a depender

das freqüências de excitação e da correspondente resposta da estrutura, classificam-se as

análises em problema de propagação de onda, de dinâmica estrutural ou quase estático

(COOK e outros, 1989). É importante observar que essa é uma divisão criada para

facilitar os estudos dos diversos fenômenos físicos existentes em dinâmica das

estruturas. Na natureza, o que se observa é uma combinação desses casos gerando a

resposta do sistema. A observação e o bom senso do engenheiro é que indicará qual dos

casos é mais predominante na análise a ser realizada.

Nos casos em que o carregamento e a resposta da estrutura são ricos em altas

freqüências, os efeitos das ondas de tensão são do interesse da engenharia. Esses são os

chamados problemas de propagação de onda que geralmente ocorrem nas situações de

impacto, cargas explosivas, mudanças bruscas de condição de contorno e outros. O

tempo de análise é usualmente curto, se restringido ao tempo que uma onda transversal

leva para percorrer a estrutura (COOK e outros, 1989).

Os problemas de dinâmica estrutural, também denominados inerciais, têm as

respostas governadas por modos de baixas freqüências. São os problemas onde os

efeitos da inércia são importantes, mas o foco da análise está direcionado para as

conseqüências de cargas tipo harmônicas, rampas, ou seja, cujas freqüências estejam na

faixa das mais baixas freqüências naturais da estrutura. O tempo de análise, nesses

casos, é bem maior do que nos problemas de propagação de onda (COOK outros, 1989).

Os efeitos inerciais no caso de problemas quase estáticos são bastante reduzidos

devido às baixas freqüências de excitação. Se conveniente, pode-se desprezar esses

efeitos e realizar uma análise estática nesse caso, simplificando o processo de solução.

COOK e outros (1989) limitam essas análises aos casos em que as freqüências de

excitação são menores que um terço das menores freqüências naturais da estrutura.

2.3. Discretização espacial

Para solucionar as equações diferenciais parciais geradas pela modelagem e obter

a resposta do sistema, usa-se aproximar os campos contínuos no espaço por variáveis

discretas, através de técnicas de aproximação como, por exemplo, o Método dos

Elementos Finitos ou o Método das Diferenças Finitas. Esse processo de discretização

no espaço independente do tempo é dito semidiscretização (HUGHES &

BELYTSCHKO, 1983; ZIENKIEWICZ & MORGAN, 1983). Neste trabalho adota-se o

Método dos Elementos Finitos como a alternativa para discretização espacial.

O Método dos Elementos Finitos (MEF) é um dos Métodos dos Resíduos

Ponderados, tendo como base a forma fraca do Método de Galerkin, que do ponto de

vista da mecânica corresponde ao princípio dos deslocamentos virtuais. Consiste

basicamente em dividir o domínio em sub-regiões, ditas elementos, cujo

comportamento é conhecido, conectadas por alguns pontos, os nós, através dos quais

interagem entre si, como apresentado na Figura 2.1.

A solução de um problema, através da formulação em deslocamentos do Método

dos Elementos Finitos, pode ser sistematizada através das etapas a seguir:

Divide-se o domínio do problema em subdomínios denominados elementos

finitos, conectados entre si através de um número finito de pontos denominados

pontos nodais ou nós.

Aproxima-se a distribuição da variável cuja solução é procurada por uma função

particular, função de interpolação.

Relaciona-se o valor da variável nos nós de cada elemento com as suas

propriedades e geometria, dando origem ao sistema de equações do elemento.

Associam-se as equações dos elementos considerando suas conexões através dos

pontos nodais. Tem-se um sistema global de equações para o problema.

Induz-se no contorno os valores conhecidos da variável do problema (introdução

das condições de contorno).

Resolve-se o sistema de equações global, obtendo-se os valores da variável do

problema nos nós.

Figura 2.1 – Idéia básica do Método dos Elementos Finitos.

● ●

●●

Elemento

Nó

Domínio

Contorno

2.3.1. Formulação em deslocamentos do Método dos Elementos Finitos

Como comentado anteriormente, o domínio do problema é subdividido em uma

série de elementos finitos. O comportamento de cada elemento é aproximado por

funções de interpolação. Tais funções devem ser escolhidas de modo a garantir a

continuidade de todo meio discretizado. A qualidade da solução aproximada depende do

grau de refinamento das funções de interpolação e da malha de elementos finitos.

A formulação em deslocamentos do método dos elementos finitos supõe que no

interior do elemento o valor dos deslocamentos e das deformações seja dado de uma

forma pré-estabelecida, em função dos deslocamentos nodais. O campo de

deslocamentos em cada elemento é aproximado por:

Nuu = ,

( 2.1 )

onde

é o campo de deslocamentos no interior do elemento que é função do tempo e

do espaço;

N é a matriz das funções de interpolação, função apenas do espaço e u é o

vetor de deslocamentos nodais, função apenas do tempo.

No caso de análise dinâmica, a discretização estende-se às derivadas do campo

de deslocamento. Têm-se:

uNu

= e

( 2.2 )

uNu

&&&&

( 2.3 )

onde

é o vetor de velocidades no interior do elemento; u

é o vetor de acelerações no

interior do elemento;

é o vetor de velocidades nodais e u

é o vetor de acelerações

nodais.

O campo de deformações em cada elemento está relacionado com os

deslocamentos nodais através da relação:

Buε = ,

( 2.4 )

onde

ε é o campo de deformações e B é a matriz que relaciona deslocamentos nodais e

deformações.

Para um corpo deformável em equilíbrio, o Princípio dos Trabalhos Virtuais

estabelece, em cada instante de tempo, a igualdade entre o trabalho das forças internas e

das forças externas durante o desenvolvimento de um campo de deslocamentos virtuais,

compatível com os vínculos externos e a continuidade do corpo (COOK e outros, 1989).

As forças internas de um corpo em equilíbrio podem ser obtidas a partir do

estado de tensões e das forças inerciais presentes em cada ponto. Ao campo de

deslocamentos virtuais corresponde um estado de deformação virtual associado. O

trabalho das forças internas é dado por:

()

∫∫

ρδ+δ=δ

interno

dVdV,W uuεεσε

&&&

( 2.5 )

sendo

interno

Wδ o trabalho virtual das forças internas; ε

o vetor de deformações

virtuais;

uδ o vetor de deslocamentos virtuais;

(

)

εεσ

o tensor de tensões, função das

deformações e das taxas de deformações e

a densidade mássica do material.

O trabalho das forças externas é dado por:

)(dS)(dV)(W

externo

ufuufuufu

∑

∫∫

δ+δ+δ=δ ,

( 2.6 )

sendo

externo

Wδ o trabalho virtual das forças externas; )(

uf o vetor de forças de

volume; o )(

uf o vetor de forças de superfície e )(

uf o vetor de cargas

concentradas.

No Método dos Elementos Finitos, para atender a condição de equilíbrio em

todo o domínio, faz-se o somatório da contribuição de cada elemento. Portanto,

considerando que as forças externas sejam representadas pelas forças nodais resultantes

dos carregamentos e fazendo o trabalho das forças externas igual ao trabalho das forças

internas, tem-se a equação obtida pelo Princípio dos Trabalhos Virtuais:

()

nelem

SddVdVdV, ffNfNuNNεεσB

∑∑

∫

∑

∫

∑

∫

∑

∫

=====

++=ρ+

&&&

( 2.7 )

onde

nelem é o número de elementos do modelo de elementos finitos e ρ é a densidade

mássica do material.

A equação ( 2.7 ) pode ser resumida na forma apresentada na equação ( 2.8 ):

)(),(

extint

ufuufuM =+

&&&

( 2.8 )

As ações externas estão apresentadas como dependentes dos deslocamentos,

porém neste trabalho essa dependência não é considerada. No item 2.3.2 são feitas

algumas observações a respeito de fontes de não linearidades, incluindo o caso acima

citado. A equação ( 2.8 ) torna-se então:

extint

),( fuufuM =+

&&&

( 2.9 )

Na equação ( 2.10 ) apresenta-se a matriz de massa do modelo, que representa

uma discretização da distribuição contínua de massa da estrutura.

∑

∫

ρ=

nelem

dVNNM

( 2.10 )

O vetor de forças internas, função dos deslocamentos e velocidades nodais, é

representado por:

()

∑

∫

nelem

int

dV,),( εεσBuuf

( 2.11 )

Em particular, se a linearidade é assumida com relação aos deslocamentos e

velocidades, a equação ( 2.11 ) toma a forma clássica:

∑

nelem

eeint

),( uCuKuuf

( 2.12 )

onde

K é a matriz de rigidez do elemento e

C a matriz de amortecimento.

Na equação ( 2.13 ) apresenta-se a parcela correspondente às ações externas,

função apenas do tempo.

nelem

ext

SddV ffNfNf

∑∑

∫

∑

∫

===

++=

( 2.13 )

A equação ( 2.8 ) ou ( 2.9 ), conhecida por equação de movimento, foi resultado

de uma discretização no espaço, porém, para completar o processo de solução, esse

novo sistema deve ser resolvido. No item 2.4, algumas técnicas de integração no tempo

são apresentadas.

2.3.2. Considerações sobre análises lineares e não lineares

Na natureza, os fenômenos são em geral não lineares. Há casos, porém, em que a

influência das fontes de não linearidades não é relevante para a análise. Nesses casos, a

análise linear é geralmente utilizada devido a sua simplicidade para formulação e

implementação.

Neste trabalho, objetiva-se abordar problemas de natureza não linear, portanto, a

formulação utilizada permite a inclusão de não linearidades, porém, restritas ao vetor de

forças internas. Fenômenos como não linearidade física (do material), geométrica,

geradas por problemas de contato, atrito, interação fluido-estrutura e outros podem ser

abordados pela formulação apresentada.

Cargas externas dependentes da geometria, também denominadas cargas

seguidoras, não estão no âmbito deste trabalho. Esse tipo de não linearidade é causa de

alterações no sistema de equações como a perda da simetria das matrizes do sistema

(BELYTSCHKO & HUGHES, 1983), condição não tratada pelos algoritmos

apresentados.

2.4. Discretização temporal

Após ter partido de um sistema de equações diferenciais parciais e ter sido

utilizada a técnica de discretização no espaço por elementos finitos, obteve-se um novo

sistema de equações, porém, um sistema de equações diferenciais ordinárias. O

problema a ser abordado tornou-se, portanto, um problema de valor inicial apenas,

representado pela solução da equação ( 2.8 ).

Há vários métodos para obtenção da resposta dinâmica após um processo de

semidiscretização, como o realizado no item anterior. Esses métodos podem ser

divididos em duas principais categorias, os métodos modais e os métodos diretos.

2.4.1. Métodos Modais

O Método da Superposição Modal consiste em transformar as equações de

equilíbrio expressas em coordenadas físicas em novas equações expressas em

coordenadas generalizadas (BATHE, 1996; COOK e outros, 1989), usualmente os

modos naturais de vibração da estrutura. Nessas coordenadas as equações tornam-se

desacopladas, podendo ser resolvidas separadamente, se não houver amortecimento, ou

se este for ortogonal, como ocorre com o amortecimento de Rayleigh. Esse método se

restringe aos casos de dinâmica linear de estruturas e se todos os modos forem

integrados, a resposta obtida não possui aproximações.

2.4.2. Métodos Diretos

Os métodos diretos são aplicados a problemas lineares ou não lineares e

consistem na obtenção da solução aproximada da equação de movimento em

determinados instantes de tempo discretos. WEAVER e JOHNSTON (1987) visualizam

as técnicas numéricas utilizadas para tal objetivo como um tipo de aproximação por

diferenças finitas. COOK e outros (1989) também apresentam essa abordagem para

definição dos métodos diretos.

Uma abordagem mais moderna sobre esse tema é dada por ZIENKIEWICZ e

MORGAN (1983) que apresentam o procedimento para discretização no tempo como

uma aplicação da técnica de Resíduos Ponderados, onde, a depender da escolha das

funções de ponderação, recai-se no método das diferenças finitas ou no método dos

elementos finitos. No último caso, procede-se de forma semelhante à realizada na

discretização espacial, usando-se, neste caso, elementos finitos lineares para representar

o domínio do tempo.

Em ambas abordagens, conhecendo-se as equações governantes e as condições

iniciais, obtém-se a solução para o final do primeiro incremento de tempo, ou elemento

de tempo. Repete-se o procedimento para os demais incrementos, adotando-se a solução

anterior como condição inicial do passo seguinte.

Capítulo 3

3. Algoritmos de integração para dinâmica não

linear

3.1. Introdução

Apresenta-se nos itens que seguem uma revisão dos principais tipos de

algoritmos, suas classificações, formulações e estratégias de solução para os sistemas de

equações gerados pelas formulações. Dentre os métodos explícitos, apresenta-se o

algoritmo do Método das Diferenças Centrais, enquanto que na classe de métodos

implícitos, apresenta-se a Regra Trapezoidal, ambos membros da família Newmark de

algoritmos de integração.

Após a revisão dos algoritmos acima citados, abordam-se os conceitos da técnica

elemento por elemento, aplicando-a as alternativas anteriormente apresentadas.

Ao final, apresentam-se algumas alternativas para medidas de erros, utilizadas

como critério de convergência ou para verificação da qualidade da resposta gerada pelos

algoritmos.

3.2. Algoritmos de integração no tempo

Como abordado nos capítulos anteriores, a análise dinâmica é regida por

equações diferenciais parciais. O uso da semidiscretização por elementos finitos permite

transformar o sistema de equações diferenciais parciais em um problema de valor inicial

apenas. O sistema de equações para análise não linear de estruturas resultante de uma

discretização no espaço por elementos finitos está apresentado na equação ( 2.9 ), cujas

condições iniciais estão expressas em:

)0( uu = e

( 3.1 )

)0( vu =

( 3.2 )

Os algoritmos de integração, aqui abordados, são métodos diretos de

discretização no tempo, consistindo em obter de forma aproximada o valor das variáveis

em questão (os deslocamentos e suas derivadas) em determinados intervalos de tempo.

A equação de movimento ( 3.3 ) expressa o equilíbrio estático para um determinado

instante de tempo:

extnnintn

),( fuufuM =+

&&&

( 3.3 )

Nos métodos diretos o espaço de tempo foi dividido em intervalos ditos

∆

e o

subscrito

n existente na equação ( 3.3 ) identifica o tempo tn

∆

⋅

Vários métodos diretos para integração das equações de movimento foram

desenvolvidos, porém a escolha do método mais adequado depende principalmente do

tipo de análise dinâmica que se deseja realizar.

Comumente classificam-se os diversos algoritmos em dois principais grupos, os

algoritmos explícitos e os algoritmos implícitos. Neste trabalho, utiliza-se a definição de

COOK e outros (1989) e BATHE (1996) para a classificação dos algoritmos de

integração no tempo.

Os algoritmos explícitos são aqueles em que as variáveis no intervalo de tempo

seguinte são determinadas apenas em função das variáveis obtidas nos intervalos de

tempo passados. Essa definição pode ser expressa pela relação funcional:

()

&&&

,,,,

1nnnn1n −+

= uuuufu

( 3.4 )

Nos algoritmos implícitos o valor da incógnita base no intervalo de tempo 1n

é dependente do seu próprio valor, além da história ao longo dos tempos passados. Essa

definição pode ser resumida na relação:

()

&&&

,,,,

n1n1n1n1n

uuuufu

++++

( 3.5 )

Nos itens a seguir apresenta-se uma descrição mais detalhada das características

de cada classe de algoritmos de integração.

3.3. Algoritmos explícitos

Como apresentado no item 3.2, a principal característica das técnicas de

integração explícitas é o fato da solução no tempo

∆

ser obtida considerando-se as

condições de equilíbrio do tempo t. Outras características importantes dessa formulação

foram apresentadas por COOK e outros (1989) e BATHE (1996) e estão citadas nos

parágrafos seguintes.

Tradicionalmente os algoritmos explícitos apresentam baixo custo

computacional por intervalo de tempo, pois o sistema de equações gerado pela

formulação é geralmente desacoplado. Não há portanto a necessidade de processos

iterativos para obtenção da solução no passo de tempo t + ∆t. Porém esse

desacoplamento está condicionado ao fato das matrizes de massa e de amortecimento

serem diagonais. Essa restrição, em alguns tipos de problemas, representa um sério

obstáculo para uso dos algoritmos explícitos, já que a sua não observância implica em

custos maiores por intervalo de tempo.

O desacoplamento anteriormente citado permite que os cálculos sejam feitos no

nível do elemento, sem a montagem de matrizes globais do sistema. Essa característica

além de ser bastante vantajosa nos casos de problemas de grande porte, já que

representa uma importante economia de memória computacional, torna a

implementação desses algoritmos mais simples do que os implícitos.

Quanto à estabilidade, os algoritmos explícitos apresentam uma restrição ao

valor de incremento de tempo máximo usado para cada análise. Essa estabilidade

condicional, garantida para problemas lineares, pode levar a pequenos valores para o

intervalo de tempo. O valor crítico para

∆

é inversamente proporcional à máxima

freqüência natural do sistema discreto, independente do tipo de solicitação aplicada à

estrutura. Portanto, em problemas quase estáticos ou de dinâmica estrutural, os métodos

explícitos tornam-se pouco eficazes já que o

∆

crítico é bem menor que o necessário

para uma integração razoavelmente precisa dos modos solicitados. Em problemas de

propagação de onda, o

t∆ necessário para uma precisa integração das freqüências

existentes na resposta do sistema é equivalente (ou menor) ao necessário para a garantia

da estabilidade. Nesses casos os algoritmos explícitos tornam-se mais adequados.

Em problemas não lineares, o valor crítico para o intervalo de tempo é variável

com o tempo, já que é atrelado às freqüências naturais do sistema que se modificam ao

longo da análise (ZIENKIEWICZ e MORGAN, 1983).

3.3.1.

Algoritmo da diferença central

O Método da Diferença Central está entre os mais populares algoritmos

explícitos usados em mecânica computacional. Sua formulação é baseada em

aproximações por diferenças centrais para as velocidades e acelerações (COOK e

outros, 1989):

()

1n1nn

−+

−

∆⋅

= uuu

( 3.6 )

()

1nn1n

−+

+⋅−

∆

= uuuu

( 3.7 )

Essas expressões são geradas por uma interpolação quadrática no tempo para os

vetores de deslocamento entre os instantes de tempo t e t +

∆

A manipulação das relações ( 3.6 ) e ( 3.7 ) em conjunto com a equação de

movimento gera o algoritmo apresentado na Figura 3.1. Alguns autores, a depender de

como se procede a manipulação das equações, apresentam diferentes variantes do

Método da Diferença Central. O algoritmo implementado neste trabalho (Figura 3.1)

tem a forma da variante 1 apresentada por KRYSL e BELYTSCHKO (1998).

A. INICIALIZAÇÕES:

1. Define o intervalo de tempo:

tt ∆<∆

2. Inicializa os vetores

u e

(condições iniciais)

3. Inicializa a matriz de massa (

M )

4. Calcula o vetor de forças desequilibradas:

),(

00int

ext0

uuffr

−=

5. Calcula:

rMu ⋅=

−

PARA CADA INCREMENTO DE TEMPO:

1. Calcula um preditor de velocidade:

tuuu

&&&&

∆+=

2. Calcula o vetor deslocamentos:

t uuuu

&&&

∆

+∆+=

3. Calcula o vetor de forças desequilibradas:

),()t(

1nintext

−= uuffr

4. Calcula o vetor de acelerações:

rMu

−

5. Calcula o vetor de velocidades:

()

1nnn

∆

+= uuuu

&&&&&&

6. Volta para B.1

C. FIM

Figura 3.1 – Algoritmo do Método das Diferenças Centrais.

No algoritmo acima apresentado, observa-se que se a matriz de massa for

diagonal, não é preciso uma técnica especializada para solução do sistema de equações

lineares, já que este se torna desacoplado. É tradicional o uso de matrizes de massa

diagonais em algoritmos explícitos, evitando-se a solução do sistema de equações em

cada incremento de tempo. A diagonalização dessas matrizes torna o custo

computacional por intervalo de tempo, para esses algoritmos, consideravelmente menor

que para os algoritmos implícitos.

Como mencionado anteriormente, o Método das Diferenças Centrais é

condicionalmente estável, exigindo valores de intervalos de tempo relativamente

pequenos. BATHE (1996) mostra que para o Método das Diferenças Centrais o valor

crítico para o intervalo de tempo (

∆

), o qual não deve ser ultrapassado, é dado por:

=∆

( 3.8 )

onde

T é o menor período do sistema de elementos finitos com m graus de liberdade.

Em problemas de propagação de onda, casos de melhor performance dos

métodos explícitos, é comum o aparecimento de altas freqüências que não têm relação

com a natureza física do problema, sendo causadas por erros provenientes do uso de

ferramentas numéricas. O Método das Diferenças Centrais aqui apresentado não possui

mecanismos para dissipação dessas freqüências.

3.4. Algoritmos implícitos

A principal característica dos algoritmos implícitos, já citada anteriormente,

reside no fato do valor da incógnita base no intervalo de tempo 1n

ser dependente do

seu próprio valor, além da história ao longo dos tempos passados (COOK e outros,

1989). Essa dependência gera um sistema não linear de equações representado por:

()

0xΦ =

+1n

( 3.9 )

onde

1n+

x é a variável independente.

A solução do sistema ( 3.9 ) para cada passo de tempo representa a principal

desvantagem dos algoritmos implícitos. Esse processo implica em um alto custo

computacional por instante de tempo em comparação com os algoritmos explícitos.

Outra conseqüência é a necessidade de montagem das matrizes globais do problema,

gerando um maior uso de memória computacional, o que pode provocar a inviabilidade

do algoritmo em problemas de grande porte.

Embora existam essas desvantagens, os algoritmos implícitos tornaram-se

populares em análises dinâmicas. Características como a estabilidade incondicional,

presente na maioria dos algoritmos implícitos, os tornam competitivos e, para certos

tipos de problemas, mais adequados que os algoritmos explícitos.

Em problemas lineares a estabilidade incondicional permite a livre escolha do

valor do incremento de tempo, ficando este apenas limitado pelos critérios de precisão

(COOK e outros, 1989). Nos casos em que o valor de

∆

máximo para garantia da

estabilidade é bem menor que o necessário para precisão desejada, os algoritmos

implícitos apresentam-se vantajosos em relação aos explícitos. Essas são geralmente

situações de análises quase estáticas ou de dinâmica estrutural, onde as baixas

freqüências envolvidas no problema permitem valores de incremento de tempo

relativamente altos para uma razoável precisão. Porém, mesmo nesses casos, quando a

análise envolve não linearidades, a estabilidade dos algoritmos implícitos não fica

garantida

a priori, mas estes têm a tendência de serem mais estáveis que os explícitos,

gerando menos dificuldades para a obtenção da resposta do sistema.

Quanto ao uso e à implementação, nos casos lineares, os algoritmos implícitos

não oferecem maiores dificuldades que os explícitos, porém, ao se tratar de problemas

não lineares, a busca do equilíbrio em cada instante de tempo passa a exigir técnicas

mais elaboradas para solução de sistemas não lineares (COOK e outros, 1989 e

SUBBARAJ e DOKAINISH, 1989).

No item a seguir, apresentam-se os métodos de solução de sistemas não lineares

utilizados neste trabalho.

3.4.1.

Técnicas de solução de sistemas não lineares

Na formulação em elementos finitos para dinâmica das estruturas, as variáveis

em questão no sistema de equações ( 3.9 ) são aceleração, velocidade e deslocamento

nodais no tempo n + 1. Porém apenas uma dessas variáveis é independente. Dessa

forma, manipulando-se as relações de dependência, que serão tratadas posteriormente,

tanto o deslocamento como a aceleração podem ser escolhidos. O sistema de equações

não lineares ( 3.9 ) torna-se então:

(

)

ext1n1n1nint1n1n1n

)(,)()(

++++++

−

+= fuuufuuMuΦ

&&&

( 3.10 )

se a formulação for na versão deslocamento, ou, analogamente,

()

ext1n1n1n1nint1n1n

)(),()(

++++++

−

+= fuuuufuMuΦ

&&&&&&&&&

( 3.11 )

se a formulação for na versão aceleração.

Neste trabalho o Método de Newton-Raphson e Gradientes Conjugados Não

Linear (PAPADRAKAKIS e GHIONIS, 1986) são utilizados na obtenção da solução do

sistema não linear ( 3.9 ).

•

Método de Newton-Raphson

Inicialmente aborda-se a formulação do Método de Newton-Raphson, cuja idéia

básica é aproximar o sistema não linear em série de Taylor e o utilizar um processo

iterativo para obter uma solução dentro de uma tolerância admitida.

Aproximando-se o sistema

)(

1n+

uΦ , para a primeira ordem, tem-se:

(

)

0uu

uΦ

uΦuΦ

=−

∂

+≈

1n1n

)(

)()(

( 3.12 )

onde

ext

1nint

),()(

++++

−+= fuufuMuΦ

&&&

( 3.13 )

int

1n1n

)(

∂

uΦ

( 3.14 )

Nessas equações, o índice inferior indica o instante de tempo

1n + , enquanto que

o índice superior indica a k-ésima iteração.

Neste caso, utiliza-se como variável independente o deslocamento.

Analogamente, para aceleração, obtém-se:

(

)

0uu

uΦ

uΦuΦ

=−

∂

+≈

1n1n

)(

)()(

&&&&

( 3.15 )

onde

ext

1nint

),()(

++++

−+= fuufuMuΦ

&&&&&

( 3.16 )

int

1n1n

)(

∂

uΦ

&&&&

( 3.17 )

O sistema linearizado é tradicionalmente denominado sistema efetivo e

representado por:

1n +++

=∆ fuK ou

( 3.18 )

1n +++

=∆ fuK

( 3.19 )

sendo

1n1n

)(

∂

uΦ

( 3.20 )

1n1n

)(

∂

uΦ

&&&&

( 3.21 )

),(

1nint

ext

1n +++

−−= uufuMff

&&&

( 3.22 )

onde,

1n+

K é a matriz de rigidez efetiva na versão deslocamento da formulação,

1n+

é a matriz de rigidez efetiva na versão aceleração e

1n+

é o vetor de forças efetivo.

Nas equações ( 3.14 ) e ( 3.17 ), o termo

int

∂

representa a variação das forças

internas com os deslocamentos, sendo geralmente denominado de matriz de rigidez

tangente (

1n+

). Da mesma forma, o termo

int

∂

representa a variação das forças

internas com as velocidades, sendo portanto denominado matriz de amortecimento

tangente (

1n+

Ct ). No caso de análise linear, as matrizes de rigidez e amortecimento

tangentes permanecem constantes. A denominação “tangente” para essas matrizes perde

o sentido, passando assim a serem representadas respectivamente por

. A matriz

de rigidez efetiva passa a ser expressa na seguinte forma:

MtK +

∂

&&&

( 3.23 )

∂

CMtK

&&&&

( 3.24 )

•

Método dos Gradientes Conjugados Não Linear

O Método dos Gradientes Conjugados pode ser visto como um problema de

minimização de uma função contínua

)(f x cujo gradiente possa ser calculado

(PAPADRAKAKIS e GHIONIS, 1986). O vetor

x que minimiza essa função equivale

ao vetor solução do sistema não linear de equações:

0xg =)( ,

( 3.25 )

onde

)(xg é o gradiente de )(f x , que, considerando a formulação em elementos finitos

para dinâmica das estruturas, toma a forma:

(

)

ext1n1n1nint1n1n1n

)(,)()(

++++++

−

+= FuuuFuuMug

&&&

( 3.26 )

Nas equações seguintes deste item omitem-se os subscritos que indicam o tempo

1n +

para tornar a notação mais leve.

A estratégia para minimização da função escalar consiste na busca sucessiva de

mínimos em dadas direções. A Figura 3.2 ilustra as etapas desse processo.

Figura 3.2 – Estratégia do Método dos Gradientes Conjugados.

A trajetória na qual se realiza a busca pelo mínimo é dada pela linha que passa

por

u e percorre a direção

p :

kkk1k

puu α+=

( 3.27 )

sendo

α o escalar que leva o ponto

u ao mínimo de )(f x na direção

p , ou seja,

que minimiza

)(f

kkk

pu α+

No Método dos Gradientes Conjugados, o vetor direção de busca

1k +

p é

expresso por:

kk1k1k

pgp β+−=

( 3.28 )

Ponto de partida

Direção de busca

Mínimo na

direção

Mínimo na

direção

sendo

β um escalar que define a nova direção de busca e )(

1k1k ++

= ugg .

SHEWCHUK (1994) apresenta em seu trabalho uma descrição detalhada das

características do Método dos Gradientes Conjugados.

Alguns autores apresentam sugestões para valores de

β , cuja melhor escolha

permanece como alvo de investigações. PAPADRAKAKIS e GHIONIS (1986) fazem

comparações entre as expressões de Fletcher-Reeves, Polak-Ribière e Soreson. Neste

trabalho opta-se pela expressão sugerida por SORENSON (1969), dada por:

()

k1k

ggp

ggg

−

=β

( 3.29 )

Definida a direção de busca, o mínimo na trajetória ( 3.27 ) é obtido buscando-se

anular a derivada direcional da função escalar

(

)

xf , ou seja:

)(f)(f

1k1k

α∂

∂

α∂

∂

( 3.30 )

0)(

)(f

kT1k

α∂

∂

pug

( 3.31 )

Substituindo o valor de

1k+

u em ( 3.31 ), tem-se a seguinte equação:

0)(

kkk

=α+ pugp .

( 3.32 )

Observa-se na equação ( 3.32 ) que não há uma expressão explícita para

α , ou

seja, uma solução analítica da equação não linear acima descrita. Neste trabalho, opta-se

inicialmente por utilizar uma ferramenta numérica para cálculo de zero de função e,

como uma segunda alternativa, por adotar uma aproximação linear da função gradiente,

obtém-se uma expressão explícita para

α .

No primeiro caso, utilizam-se técnicas de ordem zero para determinação de zero

de função. A técnica utilizada consiste em uma combinação dos métodos da bisseção,

secante e interpolação quadrática inversa (MATLAB, 2000), encontrando o zero de uma

função de uma variável, dentro de um intervalo no qual a função muda de sinal. Caso o

intervalo não seja conhecido, o próprio algoritmo se encarrega de encontrá-lo. Neste

caso, conhece-se um ponto de partida próximo do qual um intervalo válido é

investigado.

Observa-se que, ao se calcular o zero da função ( 3.32 ), não é necessária a

montagem da matriz de rigidez do problema, já que essa equação envolve apenas o

cálculo do gradiente, dado na equação de movimento ( 3.26 ).

No segundo caso, a função gradiente na direção

p é aproximada por:

+≈α

=α 0

)(d

)0()(

( 3.33 )

⋅

∂

α∂

+≈α

=α

d)(

)0()(

( 3.34 )

o que resulta em:

pHp

−=α

( 3.35 )

onde

é a hessiana da função )(f x avaliada em

u , ou seja,

)(

∂

( 3.36 )

Observa-se que nesse segundo caso é necessária a determinação da hessiana de

)(f

x , que em termos de dinâmica estrutural significa a montagem da matriz de rigidez

efetiva, apresentada na equação ( 3.20 ).

Na Figura 3.3 apresenta-se o algoritmo do Método dos Gradientes Conjugados

Não Linear.

A. INICIALIZAÇÕES:

1. Inicializa o vetor

u (valor inicial de partida do processo iterativo)

2. Inicializa a direção de busca:

)(

ugp −=

INICIAM-SE AS ITERAÇÕES: ( k = 1, 2, ... )

1. Calcula

α , solução de:

0)(

kkk

=α+ pugp

2. Calcula novo ponto

kkk1k

puu α+=

3. Calcula novo gradiente

)(

1k1k ++

= ugg

4. Calcula nova direção de busca

(

)

()

k1k

ggp

ggg

−

=β

kk1k1k

pgp β+−=

5. Se não convergiu, volta para B.1

FIM

Figura 3.3 – Algoritmo do Método dos Gradientes Conjugados Não Linear.

Observa-se que o passo

B.5 na Figura 3.3 consiste em realizar um teste de

convergência, cujo critério é abordado no item 3.6 que trata das medias de erros

utilizadas neste trabalho.

3.4.2. Algoritmos da família Newmark

Dentre as famílias de algoritmos de integração, a família Newmark

(NEWMARK, 1959) possui os mais populares algoritmos de integração. Consiste

basicamente em expressar as velocidades e deslocamentos segundo aproximações por

diferenças finitas no domínio do tempo, dadas por:

()

[]

1nn

nn1n

221

β+β−

∆

+∆+= uuuuu

&&&&&

( 3.37 )

(

)

[]

1nnn1n

−∆+= uuuu

&&&&&&

( 3.38 )

Os parâmetros

determinam propriedades de estabilidade e precisão dos

métodos. HUGHES (1987) apresenta uma análise de estabilidade para os métodos da

família Newmark, para problemas lineares, colocando que se

2 ≥γ≥β ,

( 3.39 )

os algoritmos são incondicionalmente estáveis.

Muitos métodos conhecidos podem ser gerados a partir de particularizações das

expressões ( 3.37 ) e ( 3.38 ). A Tabela 3.1 apresenta alguns desses métodos e suas

características. Outras informações sobre cada método em particular podem ser

encontradas em literaturas tradicionais como HUGHES (1987), COOK e outros (1989),

BATHE (1996) entre outros. Neste trabalho adota-se a Regra Trapezoidal devido à sua

estabilidade incondicional, garantida para análises lineares.

Tabela 3.1 – Métodos da família Newmark.

MÉTODO TIPO

ONDIÇÃO DE

ESTABILIDADE

ORDEM DE

PRECISÃO

Aceleração Média

(Regra Trapezoidal)

implícito

Incondicional 2

Aceleração Linear implícito

condicional 2

Fox-Goodwin implícito

condicional 2

Diferença Central explícito 0

condicional 2

As equações ( 3.37 ) e ( 3.38 ) em conjunto com a equação de equilíbrio ( 2.9 ) ,

representada também pelo sistema não linear ( 3.9 ), formam um sistema com três

equações e três incógnitas (

1n +

u ,

1n+

1n +

), considerando que os estados nos tempos

passados são conhecidos. Várias são as estratégias para a solução desse sistema, sendo

as versões aceleração e deslocamento as mais tradicionais.

Na versão aceleração, como esta é considerada a variável independente, utilizam-

se as equações ( 3.37 ) e ( 3.38 ) tal como foram apresentadas:

() ()

[]

1nn

nn1n1n

221

+++

β+β−

∆

+∆+= uuuuuu

&&&&&&&

( 3.40 )

() ()

[]

1nnn1n1n

+++

−∆

= uuuuu

&&&&&&&&

( 3.41 )

Utilizando-se a versão da linearização do sistema não linear ( 3.9 ) que considera

a aceleração como variável independente, equação ( 3.15 ), substituem-se as equações

( 3.40 ) e ( 3.41 ) na expressão da matriz efetiva ( 3.17 ), obtendo-se:

MCtKtK +⋅∆⋅γ+⋅∆⋅β=

+++

tt .

( 3.42 )

Observa-se que se for utilizado o Método das Diferenças Centrais (

0=β

) e as

matrizes de massa e amortecimento forem diagonais, o sistema efetivo perde o seu

acoplamento.

Quando a estratégia considerada admite o deslocamento como variável

independente, obtêm-se expressões de velocidade e aceleração em função do

deslocamento através da manipulação das equações ( 3.40 ) e ( 3.41 ), resultando em:

() ()

1n1n1n1n1n

+++++

−

∆⋅β

+= uuuuu

( 3.43 )

() ()

1n1n

++++

−

∆⋅β

= uuuu

( 3.44 )

onde

()

nn1n

21t

uuuu

&&&

⋅β⋅−⋅∆+⋅∆+=

( 3.45 )

(

)

nn1n

uuu

&&&&

⋅γ−⋅∆+=

( 3.46 )

Substituindo-se as equações ( 3.43 ) e ( 3.44 ) na expressão da matriz efetiva

( 3.14 ), obtém-se:

MCtKtK

∆β

+++

( 3.47 )

O algoritmo implementado neste trabalho utiliza esta última versão e está

apresentado na Figura 3.4.

A. INICIALIZAÇÕES:

1. Define o intervalo de tempo

∆

2. Define o os valores para as constantes

e γ .

3. Inicializa os vetores

u e

v (condições iniciais)

4. Inicializa a matriz de massa (

5. Calcula o vetor de forças desequilibradas:

()

00int

ext0

,uuffr

−

6. Calcula:

rMu

−

PARA CADA INCREMENTO DE TEMPO:

1. Incrementa o tempo:

ttt

∆+=

2. Calcula o vetor de forças externas:

ext

3. Calcula o preditor de deslocamento:

()

nn1n

21t

uuuu

&&&

β−∆+∆+=

4. Calcula o preditor de velocidade:

()

nn1n

uuu

&&&&

γ−∆+=

5. Soluciona o sistema não linear:

()

0uΦ =

+1n

6. Calcula os valores de velocidade e aceleração:

()

−

∆β

+= uuuu

()

−

∆β

= uuu

7. Volta para B.1

FIM

Figura 3.4 – Algoritmo de Newmark versão deslocamento.

3.5. Técnica elemento por elemento

Problemas formulados através do Método dos Elementos Finitos recaem na

solução de um sistema de equações ordinárias, linear ou não linear, a depender do

problema tratado. Em análise dinâmica de estruturas, o uso de algoritmos explícitos,

com matrizes de massa e amortecimento diagonais, gera um desacoplamento do sistema

de equações, evitando-se a montagem das matrizes e dos vetores globais da estrutura.

Essa é a principal característica da técnica elemento por elemento, a obtenção da

solução do sistema de equações sem que tenham sido montados os vetores e matrizes

globais da estrutura. Observa-se que essa característica proporciona uma razoável

economia de memória computacional em problemas de larga escala, pois evita

redundância de armazenamento de informações, tornando-se uma das principais

vantagens da técnica (JOUGLARD e COUTINHO, 1998).

Os métodos elemento por elemento foram introduzidos por HUGHES e outros

(1983) e ORTIZ e outros (1983) e posteriormente vários outros autores fizeram

contribuições ao tema como BELYTSCHKO (1984), WINGET e HUGHES (1985),

WATHEN (1989), entre outros. HUGHES e outros (1983) elaboraram uma técnica que,

através de aproximações, evita a inversão das matrizes globais da solução do sistema de

equações efetivo em problemas lineares onde a matriz de massa é diagonal.

À medida que mais pesquisadores passaram a utilizar essa técnica, abordagens

diferentes foram dadas, mas todas visando a obtenção da solução do sistema de

equações (lineares ou não lineares) gerado pela formulação em elementos finitos sem a

montagem das matrizes e vetores globais. Em geral, tem-se por objetivo transformar

operações que envolvam as soluções desses sistemas em produtos de matrizes por

vetores e produtos escalares entre vetores comuns nos cálculos locais dos elementos.

Nos métodos implícitos, a solução do sistema de equações não lineares acoplado

apresenta-se como a principal barreira à técnica elemento por elemento. A escolha do

método de solução desse sistema passa a ser um importante fator, facilitando ou

dificultado o uso da técnica. Quando o sistema de equações não lineares é linearizado

(solução pelo método de Newton-Raphson), recai-se na solução de um sistema de

equações lineares em cada passo do processo iterativo. Uma das alternativas para esse

caso é utilizar um algoritmo iterativo para solução do sistema de equações lineares,

como, por exemplo, o método dos gradientes conjugados, já que nas suas operações

dominam os produtos de matrizes por vetores e produtos escalares entre vetores. Uma

outra alternativa à essa é a otimização do método dos gradientes conjugados não linear.

Cada uma dessas técnicas será detalhada a seguir.

3.5.1.

Formulação elemento por elemento em algoritmos explícitos

Nos algoritmos explícitos, poucas alterações são necessárias para a aplicação da

técnica elemento por elemento. Nesses algoritmos necessita-se apenas da montagem do

vetor de forças desequilibradas:

),(

1n1nint

ext1n ++

−= uuffr

( 3.48 )

que pode ser gerado a partir de contribuições individuais de cada elemento,

∑

nelem

e1n

( 3.49 )

onde

r é o vetor de forças desequilibradas do elemento, expandida para o sistema

global, sendo dado por:

rLr =

( 3.50 )

onde

L é a matriz de incidência cinemática do elemento e

r o vetor de forças

desequilibradas local do elemento, com dimensão

local

glbn x 1. O símbolo glbn

significa o número de graus de liberdade, ora global, da estrutura, ora local, do elemento

individualmente.

Essas alterações são realizadas no itens A.4 e B.3 do algoritmo do Método das

Diferenças Centrais apresentado na Figura 3.1.

3.5.2.

Formulação elemento por elemento em algoritmos implícitos

Devido à solução do sistema de equações não lineares em cada passo de tempo, o

uso da técnica elemento por elemento requer algumas alterações na estrutura

tradicionalmente utilizada no caso dos algoritmos implícitos. Apresentam-se neste item

algumas estratégias aplicadas ao algoritmo de Newmark, sem que haja perda de

generalidade.

Para solução do sistema não linear, gerado pela formulação em elementos finitos,

consideram-se o Método de Newton-Raphson e o Método dos Gradientes Conjugados

Não Linear, ambos abordados no item 3.4.1. Como cada método tem suas

características particulares, estratégias diferentes são usadas em cada um.

•

Método de Newton-Raphson

No Método de Newton-Raphson a principal dificuldade para o uso direto da

técnica elemento por elemento é a solução do sistema linearizado acoplado em cada

incremento de tempo, ou seja, a solução do sistema efetivo:

1n +++

=∆ fuK .

( 3.51 )

Torna-se necessário o uso de técnicas de solução de sistemas lineares que sejam

adaptáveis à técnica elemento por elemento. O uso do Método dos Gradientes

Conjugados para sistemas lineares apresenta-se como uma alternativa viável para esse

fim. A Figura 3.5 apresenta o algoritmo desse método. Nesta figura foram omitidos os

índices que indicam o tempo

com o objetivo de simplificar a notação.

A. INICIALIZAÇÕES:

1. Inicializa o vetor

u∆

(valor inicial do processo iterativo)

2. Inicializa o gradiente:

1111

fuKg −∆=

3. Inicializa a direção de busca:

gp −=

INICIAM-SE AS ITERAÇÕES: ( k = 1, 2, ... )

1. Realiza a busca linear

pKp

pg−

=α

2. Calcula novo ponto de mínimo

kkk1k

puu α+=∆

3. Calcula novo gradiente

1k1k1k1k ++++

−∆= fuKg

4. Calcula nova direção de busca

=β

kk1k1k

pgp β+−=

5. Se não convergiu, volta para B.1

FIM

Figura 3.5 – Algoritmo do Método dos Gradientes Conjugados Linear

Assim como na Figura 3.3, o passo B.5 da Figura 3.5 consiste em realizar um

teste de convergência, cujo critério é abordado no item 3.6 que trata das medias de erros

utilizadas neste trabalho.

A expressão para cálculo da direção de busca é a sugerida por FLETCHER AND

REEVES (1964).

Observa-se que as operações realizadas durante a solução iterativa do sistema

envolve apenas produtos matriz-vetor ou vetor-vetor. Os itens A.2 e B.3 da Figura 3.5

apresentam a formação semelhante ao apresentado para os métodos explícitos. O item

B.1 apresenta uma multiplicação vetor-matriz-vetor, dada por:

aux pKp= ,

( 3.52 )

sendo

aux um escalar auxiliar, correspondente ao denominador de

α .

A matriz efetiva

K pode ser expressa como um somatório das contribuições de

cada elemento, ou seja:

∑

nelem

( 3.53 )

onde

K é a matriz efetiva do elemento, expandida para o sistema global, sendo dada

por:

LKLK =

( 3.54 )

onde

K é a matriz efetiva local do elemento, com dimensão

local

glbnx

local

glbn.

O produto ( 3.52 ), omitindo-se os sobrescritos k para simplificar a notação, pode

ser então escrito como:

pKp

∑

nelem

aux

( 3.55 )

∑

nelem

aux pKp

( 3.56 )

onde

p corresponde ao vetor

tendo anuladas as posições dos graus de liberdade que

não estão conectados ao elemento

Não é necessário armazenar a matriz

que tem a mesma dimensão da matriz

global. Uma alternativa ao produto ( 3.56 ) é:

∑

nelem

aux pKp

( 3.57 )

onde

p corresponde apenas aos elementos não nulos do vetor

p , tendo a dimensão

local

glbn x 1, e

K a matriz de rigidez do elemento.

O Apêndice A apresenta um exemplo numérico de multiplicação matriz-vetor e

vetor-matriz-vetor elemento por elemento, aplicado a um sistema de molas.

Uma forma alternativa bastante popular do Método dos Gradientes Conjugados é

o uso de pré-condicionadores, que aumentam significativamente a eficiência do método

(PAPADRAKAKIS e DRACOPOULOS, 1991).

Basicamente o Método dos Gradientes Conjugados Pré-Condicionados (PCG)

consiste em resolver o sistema de equações ( 3.18 ) transformado em:

−

=∆ FPuKP

( 3.58 )

sendo P a matriz de pré-condicionamento. Essa matriz deve ser tal que o número de

condição do produto

−

⋅KP seja o mais próximo possível da unidade, assim, quanto

mais semelhante a matriz P for de

1n +

, mais próximo da unidade será o número de

condição do produto acima citado. Outra característica necessária é que um sistema

linear cuja matriz de coeficientes é

deve precisar de poucas iterações para ser

resolvido. Sendo assim, a escolha da matriz de pré-condicionamento deve satisfazer, na

medida do possível, essas condições. Ao contrário da matriz efetiva, a montagem da

matriz de pré-condicionamento como um somatório de contribuições de cada elemento

não é tão natural.

Vários autores, como PAPADRAKAKIS e DRACOPOULOS (1991), abordam o

tema pré-condicionamento e técnicas elemento por elemento. Neste trabalho essa

estratégia não é tratada. Utiliza-se o método PCG, na sua formatação tradicional, ou

seja, montando-se as matrizes globais do sistema, assim como o método direto de

solução de sistemas lineares, apenas para fins de comparação com o Método dos

Gradientes Conjugados sem pré-condicionamento. A matriz de precondicionamento

utilizada é uma matriz diagonal, formada pela diagonal da matriz

1n+

K .

•

Método dos Gradientes Conjugados Não Linear

Uma segunda alternativa à solução do sistema de equações não lineares

resultante da formulação é o uso do Método dos Gradientes Conjugados Não Linear

(NLCG), como abordado no item 3.4.1.

Neste caso, as alterações são bastante semelhantes às realizadas no tópico

anterior, quando se utiliza a forma linear do Método dos Gradientes Conjugados.

Considerando o conteúdo da Figura 3.3, que descreve o algoritmo do NLCG,

observa-se que os itens A.2, B.1 e B.3 podem ser calculados elemento por elemento.

Nos itens A.2 e B.3 são realizadas avaliações do vetor gradiente, que são

montados de forma semelhante ao vetor resíduo nos algoritmos explícitos.

A depender de como é realizada a busca linear, o item B.1 pode apresentar o

cálculo do gradiente (como feito nos itens A.2 e B.3) e a multiplicação deste pelo vetor

direção de busca, ou seja, avaliação da expressão:

0)(

kkk

=α+ pugp

( 3.59 )

ou pode apresentar as mesmas operações da versão linear do Método dos Gradientes

Conjugados:

pHp

−=α ,

( 3.60 )

na qual o denominador pode ser calculado elemento por elemento, tal como mostrado

no tópico anterior.

A realização da busca linear através da solução da equação ( 3.59 ) tem a

vantagem de ser necessário apenas calcular o gradiente da função )(f x durante todo o

processo.

As demais operações do método não envolvem matrizes e podem ser realizadas

com os vetores no sistema global, sem grande perda de memória.

3.6. Medidas de erro em análise dinâmica

Devido ao uso de um processo iterativo para solução do sistema não linear de

equações, torna-se necessária a definição de um critério de convergência. Este critério

tem a finalidade de julgar se os erros gerados pela resposta estão dentro de uma faixa

aceitável ou não, ou seja, dentro de uma tolerância admitida. Para tanto, é necessária

uma medida dos erros gerados pela resposta após cada iteração.

BELYTSCHKO e HUGHES (1983) sugerem um critério baseado no vetor de

forças desequilibradas, no qual a norma desse vetor deve ser menor que uma dada

fração da norma de um vetor de forças de referência:

ref

força

int

ext

⋅≤−− ffuMf

( 3.61 )

onde

força

é o valor da tolerância admitida para esse critério e

ref

é a norma do

vetor de forças de referência.

COOK e outros (1989) definem que o vetor de forças de referência é dado por:

ext

ref

= ff .

( 3.62 )

Porém, da forma como acima apresentado, é necessário garantir que o vetor de

forças externas não seja nulo. Assim, para considerar os casos em que não há aplicação

de cargas externas, neste trabalho utiliza-se a seguinte adaptação do vetor de forças de

referência sugerido por BELYTSCHKO e HUGHES (1983):

ext

int

ref

+= fff .

( 3.63 )

Além de necessitar de um critério de convergência, é preciso também avaliar a

qualidade da resposta gerada pelo algoritmo, já que, em sua maioria, as análises

dinâmicas não lineares não têm resposta analítica.

Uma alternativa é o uso da equação ( 3.61 ), sendo o erro dado pelo valor de

força

e . Porém, como o critério de convergência utiliza essa condição, os erros obtidos

assim seriam sempre da ordem de grandeza da tolerância admitida, mascarando a real

qualidade da resposta.

COOK e outros (1989) sugerem, como forma de medida de qualidade da

resposta em problemas não lineares, a checagem do balanço energético em cada instante

de tempo, dado por:

ext1n

int

WTW

( 3.64 )

onde

int

é o trabalho das forças internas,

é a energia cinética e

ext

e o

trabalho das forças externas.

Essa igualdade expressa que o trabalho das forças externas é convertido parte em

energia cinética e parte em trabalho das forças internas que pode conter energia

potencial elástica armazenada e energia dissipada por algum agente plástico ou viscoso.

O erro, nesse caso, seria dado por:

extn

int

energ

extn

int

WTWeWTW ++⋅≤−+ ,

( 3.65 )

adaptação da expressão sugerida por COOK e outros (1989), representando a relação

entre o erro em energia (lado esquerdo da equação) e a energia total do sistema (lado

direito da equação).

O trabalho interno é dado por:

∫

∆+

∆

t)1n(

int

dtWWW

( 3.66 )

sendo

int

W fu

= .

( 3.67 )

Aproximando-se a integral através da regra trapezoidal, tem-se:

()

int

∆

+= fufu

( 3.68 )

O trabalho externo é obtido de forma análoga ao trabalho interno, sendo,

portanto:

()

ext

∆

+= fufu

( 3.69 )

A energia cinética e calculada através da seguinte expressão:

()

T uMu

= .

( 3.70 )

Para algoritmos explícitos o erro

energ

e não deve ultrapassar o valor 05,0 , já que

isso gera instabilidades no algoritmo. No caso dos implícitos essa medida de erro

apresenta a mesma dificuldade, para medir a qualidade da resposta, que a medida

força

e ,

em força, apresenta. Opta-se então por uma medida relativa, para a qual é necessária

uma resposta considerada precisa o suficiente para que seja usada como base de

comparação.

O erro relativo é então definido como:

()

ref

rel

max

−

( 3.71 )

onde

(

)

ref

max u representa a chamada norma do máximo de um vetor coluna

(BOLDRINI et al, 1980).

A resposta analítica é tomada como base nos casos lineares, quando possível.

Nos casos não lineares, avalia-se uma faixa de valores provável para o

t∆ do método

das Diferenças Centrais, então obtém-se uma resposta com um

∆

bem menor que o

crítico. Nos exemplos deste trabalho são indicados os intervalos de tempo usados nas

simulações base.

Usa-se

rel

e como uma medida global da qualidade da resposta de um dado

algoritmo para fins de comparação.

Capítulo 4

4. Exemplos e Aplicações

4.1. Introdução

Neste capítulo apresentam-se quatro exemplos que ilustram algumas

características dos algoritmos estudados. O primeiro utiliza um sistema massa-mola

para validação das implementações realizadas, comparando-se com a solução por

superposição modal do sistema. O segundo exemplo apresenta o desempenho dos

algoritmos no caso em que se utiliza um elemento de pórtico tridimensional, para

simular uma viga engastada sujeita a grandes deslocamentos. O terceiro apresenta uma

análise quase estática de um cabo horizontal, apoiado nas extremidades, onde se observa

o comportamento dos algoritmos na obtenção de uma configuração estática e diante de

três discretizações da malha de elementos finitos. O último exemplo trata de uma

análise dinâmica, não linear, onde se avaliam as respostas dos algoritmos para um cabo

que pendula por alguns segundos. São usadas duas discretizações, três e cinco

elementos, analisando-se as respostas para cada uma.

Nesses exemplos, avaliam-se duas estratégias: solução do sistema não linear

usando o método de Newton-Raphson em conjunto com o método dos gradientes

conjugados para solução do sistema linearizado e solução do sistema não linear usando

o método dos gradientes conjugados não linear, sendo a busca linear ora realizada

analiticamente, ora iterativamente. Utilizam-se algumas formatações tradicionais, como

o uso do método direto e gradientes conjugados pré-condicionados para solução do

sistema linearizado, para fins de comparação com as três anteriormente citadas. Além

desses, utiliza-se o Método das Diferenças Centrais, explícito, como alternativa aos

implícitos.

Nesses exemplos, os algoritmos estudados são identificados através de uma

nomenclatura simplificada que está apresentada na Tabela 4.1.

Tabela 4.1 – Nomenclatura simplificada para os algoritmos estudados.

Solução do sistema não linear usando Newton-Raphson com solução

do sistema linearizado pelo método direto.

NR - Método direto

Solução do sistema não linear usando Newton-Raphson com solução

do sistema linearizado pelo método dos gradientes conjugados pré-

condicionados.

NR - PCG

Solução do sistema não linear usando Newton-Raphson com solução

do sistema linearizado pelo método dos gradientes conjugados

(convencional).

NR - CG

Solução do sistema não linear usando gradientes conjugados com

busca linear realizada iterativamente.

NLCG - Busca linear iterativa

Solução do sistema não linear usando gradientes conjugados com

busca linear realizada analiticamente.

NLCG - Busca linear analítica

Método das diferenças centrais Diferença central

Em todos os exemplos aqui apresentados utiliza-se um computador com

processador Pentium IV, 2 Gigahertz com 256 Megabytes de memória RAM.

4.2. Aplicação linear: sistema massa-mola

Trata-se de um caso de vibração livre de um sistema massa-mola com dois graus

de liberdade. O seu objetivo principal é a aferição dos algoritmos estudados, para um

caso simples de análise linear.

O sistema discreto está apresentado na Figura 4.1.

Figura 4.1 – Sistema massa-mola.

A Tabela 4.2 apresenta as propriedades do sistema utilizado neste exemplo.

Tabela 4.2 – Propriedades do sistema massa-mola.

RIGIDEZES MASSAS

= 1 m

= 1

= 10 m

= 1

O tempo de análise é de 84 segundos, sujeito às condições iniciais indicadas na

Tabela 4.3.

Tabela 4.3 – Condições iniciais do sistema massa-mola.

)0(u

A máxima freqüência do sistema é de 4,58 rad/s, o que limita o valor do

incremento de tempo para o Método das Diferenças Centrais em: 0,4364s. Opta-se por

observar a resposta obtida em quatro intervalos de tempo, sendo dois maiores que o

crítico para o MDC e dois menores que este. Os valores adotados estão apresentados na

Tabela 4.4.

Tabela 4.4 – Valores dos intervalos de tempo.

∆

0,7 s

∆

0,6 s

∆

0,4 s

∆

0,01 s

As Figuras 4.2, 4.3, 4.4 e 4.5 apresentam a história de deslocamentos do grau de

liberdade 1 obtida pela integração numérica através dos algoritmos estudados, para os

valores de intervalo de tempo acima mencionados, em comparação com a solução

gerada pelo método da superposição modal, que neste exemplo, por se ter utilizado

todos os modos de vibração do sistema, é referenciado como sendo a resposta analítica.

0 10 20 30 40 50 60 70 80 90

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Tempo

Deslocamento u1

Resposta Analítica

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.2 – História de deslocamento do grau de liberdade 1 para ∆t = 0,7s .

0 10 20 30 40 50 60 70 80 90

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Tempo

Deslocamento u1

Resposta Analítica

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.3 – História de deslocamento do grau de liberdade 1 para ∆t = 0,6s .

0 10 20 30 40 50 60 70 80 90

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Tempo

Deslocamento u1

Resposta Analítica

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.4 – História de deslocamento do grau de liberdade 1 para ∆t = 0,4s .

0 10 20 30 40 50 60 70 80 90

-0.8

-0.6

-0.4

-0.2

0.2

0.4

0.6

0.8

Tempo

Deslocamento u1

Resposta Analítica

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.5 – História de deslocamento do grau de liberdade 1 para ∆t = 0,01s .

Observando-se estas quatro figuras, não é possível distinguir as soluções geradas

pelos algoritmos implícitos estudados. De fato, neste exemplo, esses algoritmos

apresentaram a mesma solução.

As Figuras 4.2 e 4.3 não apresentam a solução gerada pelo Método das

Diferenças Centrais pois este divergiu durante o processo de solução, comportamento

esperado para algoritmos explícitos quando o intervalo de tempo não obedece à

condição de estabilidade (

∆

≤∆ ).

Para facilitar a observação da qualidade da resposta gerada pelos algoritmos,

apresenta-se graficamente, na Figura 4.6, os erros relativos para os quatro intervalos de

tempo analisados. A Tabela 4.5 apresenta os valores numéricos correspondentes à

Figura 4.6.

0,0

2,0

4,0

6,0

8,0

10,0

12,0

14,0

1234

Variação do intervalo de tempo

Erro relativo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

Figura 4.6 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado.

Tabela 4.5 – Valores dos erros relativos em cada intervalo de tempo analisado.

ERRO RELATIVO À SOLUÇÃO POR SUPERPOSIÇÃO MODAL

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto 13,26 12,07 8,16 0,64

NR - PCG 13,26 12,07 8,16 0,64

NR - CG 13,26 12,07 8,16 0,64

NLCG - Busca linear analítica 13,26 12,07 8,16 0,64

NLCG - Busca linear iterativa 13,26 12,07 8,16 0,64

Diferença central - - 8,00 0,32

- O algoritmo divergiu.

Observa-se o crescimento do erro a medida que o valor do

t∆

é maior. Esse

comportamento constata o fato de que apesar da estabilidade dos algoritmos implícitos,

nestes casos o valor do intervalo de tempo fica a critério da precisão necessária para os

objetivos da análise realizada.

Um fator importante na análise dos algoritmos é o tempo gasto para simulação.

Na Figura 4.7 e na Tabela 4.6 apresentam-se, graficamente e numericamente, os valores

obtidos em cada algoritmo para cada intervalo de tempo.

1234

Variação do intervalo de tempo

Tempo de simulação (s)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

211,78

Figura 4.7 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado.

Tabela 4.6 – Valores dos tempos de simulação (s) para cada intervalo de tempo

analisado.

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto 0,98 1,37 1,42 34,61

NR - PCG 1,16 1,22 1,58 42,24

NR - CG 1,03 1,12 1,43 38,77

NLCG - Busca linear analítica 1,14 1,30 1,58 42,95

NLCG - Busca linear iterativa 4,10 4,84 7,04 211,78

Diferença central - - 1,04 23,42

- O algoritmo divergiu.

Observa-se inicialmente o aumento do tempo de análise com a redução do valor

do incremento de tempo nos algoritmos estudados, em especial para o NLCG - Busca

linear iterativa. Constata-se também que quando a condição de estabilidade do MDC é

satisfeita, seu desempenho, com relação ao tempo de análise, é superior aos demais

algoritmos.

Uma outra característica comumente observada na avaliação de algoritmos de

integração, que possuem processos iterativos, é a média de iterações por intervalo de

tempo. A Figura 4.8 e a Tabela 4.7 apresentam esses valores.

0,5

1,5

1234

Variação do intervalo de tempo

Número de iterações por intervalo de tempo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.8 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado.

Tabela 4.7 – Valores das médias de iterações por intervalo de tempo para cada intervalo

de tempo analisado.

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto 1,03 1,03 1,02 1,00

NR - PCG 1,03 1,03 1,02 1,00

NR - CG 1,03 1,03 1,02 1,00

NLCG - Busca linear analítica 1,99 1,99 2,00 2,00

NLCG - Busca linear iterativa 1,99 1,99 2,00 2,00

Observa-se que devido à linearidade deste exemplo, o número de iterações é

muito próximo da unidade, para os algoritmos que usam o Método de Newton-Raphson,

pois o sistema linearizado para este exemplo não tem aproximações. Os algoritmos que

utilizam o Método dos Gradientes Conjugados Não Linear realizam duas iterações

devido às suas aproximações.

4.3. Viga engastada

Este exemplo tem por objetivo ilustrar o uso dos algoritmos estudados em

problemas não lineares quase estáticos. Trata-se de uma viga engastada na extremidade

esquerda, sujeita a uma força momento aplicada na extremidade direita. A Figura 4.9

apresenta o desenho esquemático da viga acima citada e a discretização utilizada.

(a) (b)

Figura 4.9 – (a) Desenho esquemático e (b) discretização da viga analisada.

O modelo é composto por dez elementos de pórtico tridimensionais,

desenvolvidos por PACOSTE e ERIKSSON (1997). As propriedades físicas e

geométricas estão apresentadas na Tabela 4.8.

Tabela 4.8 – Propriedades físicas e geométricas da viga.

Comprimento da viga (L) 10 m

Módulo de elasticidade longitudinal (E)

2,05 ⋅ 10

kN/m

Módulo de elasticidade transversal (G) E/2

Momento de inércia (I) 0,083 m

Momento de inércia polar (J) I

Área da seção transversal (circular) (A) 1 m

Densidade mássica (ρ)

7,83 Mg/m

A matriz de massa aqui utilizada é diagonalizada, possibilitando as comparações

dos resultados com o algoritmo explícito da Diferença Central.

O momento na extremidade é aplicado linearmente com o tempo como

apresentado na Figura 4.10. O tempo de aplicação é de 0,5s, o mesmo de duração da

análise.

M(t)

...

M(t)

Figura 4.10 – Curva de aplicação da força momento.

O momento de referência

M apresentado na Figura 4.10 foi obtido a partir da

solução estática do caso abordado neste exemplo. O cálculo desse valor pode ser

encontrado no Apêndice C deste trabalho.

Utilizou-se o intervalo de tempo de 1⋅10

-5

s e o algoritmo NR – Método Direto

para obter a solução de referência e se avaliar as máximas freqüências do sistema a fim

de se obter o intervalo de tempo crítico para o método da Diferença Central. A Figura

4.11 apresenta a variação das máximas freqüências do sistema ao longo do tempo de

análise. O intervalo de tempo crítico encontrado é de 1,49·10

-4

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

1.331

1.332

1.333

1.334

1.335

1.336

1.337

1.338

1.339

1.34

1.341

x 10

Máximas freqüências

Tempo

Figura 4.11 – Máximas freqüências do sistema ao longo da análise.

Os intervalos de tempo utilizados nas análises estão apresentados na Tabela 4.9.

A tolerância permitida para a solução do sistema não linear é de 1⋅10

-8

⋅

⋅π⋅=













M)t(M

M(t)

Tabela 4.9 – Valores dos intervalos de tempo.

∆

1⋅10

-2

∆

5⋅10

-4

∆

1,45⋅10

-4

Os resultados das análises são apresentados em gráficos tempo

versus

deslocamento vertical. Apresentam-se, além dos resultados dos algoritmos estudados, a

solução de referência e a solução estática para fins de comparação.

Para o primeiro intervalo de tempo os algoritmos NR - Método direto, NLCG -

Busca linear iterativa, NLCG - Busca linear analítica e Diferença central divergiram.

Apresenta-se na Figura 4.12 a solução obtida para

101t

−

⋅=∆ s.

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

-1

Tempo

Deslocamento vertical na extremidade direita (m)

Referência

Solução estática

NR-PCG

NR-CG

Figura 4.12 – Deslocamento vertical da extremidade direita da viga para

101t

−

⋅=∆ s.

A Figura 4.13 apresenta os resultados para os demais intervalos de tempo.

Observa-se no detalhe a proximidade das respostas geradas pelos algoritmos, tornando-

se difícil distingui-los visualmente na Figura 4.13 (a) e (b).

(a) (b)

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

-1

Tempo

Deslocamento vertical na extremidade direita (m)

Referência

Solução estática

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

-1

Tempo

Deslocamento vertical na extremidade direita (m)

Referência

Solução estática

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

0.185 0.1855 0.186 0.1865 0.187 0.1875 0.188 0.1885 0.189 0.1895 0.19

7.33

7.331

7.332

7.333

7.334

7.335

7.336

7.337

7.338

7.339

7.34

Tempo

Deslocamento vertical na extremidade direita (m)

Referência

Solução estática

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

0.185 0.1855 0.186 0.1865 0.187 0.1875 0.188 0.1885 0.189 0.1895 0.19

7.33

7.331

7.332

7.333

7.334

7.335

7.336

7.337

7.338

7.339

7.34

Tempo

Deslocamento vertical na extremidade direita (m)

Referência

Solução estática

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.13 – Deslocamento vertical da extremidade esquerda da viga para

(a)

105t

−

⋅=∆ e (b)

1045,1t

−

⋅=∆ s; Detalhe do deslocamento vertical da extremidade

esquerda da viga para (c)

105t

−

⋅=∆ e (d)

1045,1t

−

⋅=∆ s.

Para melhor análise da qualidade das respostas obtidas, apresentam-se na Figura

4.14 e na Tabela 4.10 os erros relativos à resposta gerada com o intervalo de tempo de

1⋅10

-5

Área do detalhe

0,000

0,002

0,004

0,006

0,008

0,010

0,012

0,014

0,016

0,018

0,020

0,022

0,024

0,026

0,028

0,030

123

Variação do intervalo de tempo

Erro relativo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

~ 0,40

Figura 4.14 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado.

Tabela 4.10 – Valores dos erros relativos em cada intervalo de tempo analisado.

ERRO RELATIVO À SOLUÇÃO COM

101t

−

⋅=∆

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto * 0,0068 0,0027

NR - PCG 0,3975 0,0068 0,0027

NR - CG 0,4071 0,0068 0,0027

NLCG - Busca linear analítica * 0,0261 0,0042

NLCG - Busca linear iterativa * 0,0266 0,0037

Diferença central * * 0,0016

* O algoritmo divergiu.

O tempo de simulação de cada algoritmo está apresentado na Figura 4.15 e na

Tabela 4.11. Observa-se que o grupo de algoritmos que utilizam o Método de Newton-

Raphson para solução do sistema não linear apresenta menor tempo de simulação que o

grupo que utiliza o Método dos Gradientes Conjugados Não Linear.

100

200

300

400

500

600

700

800

900

1000

123

Variação do intervalo de tempo

Tempo de computação (s)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

2.081,35 17.230,00

7.168,84 70.571,0

Figura 4.15 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado.

Tabela 4.11 – Valores dos tempos de simulação (s) para cada intervalo de tempo

analisado.

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto * 250,89 437,97

NR - PCG 93,71 254,90 441,14

NR - CG 105,67 264,30 468,96

NLCG - Busca linear analítica * 2081,35 7168,84

NLCG - Busca linear iterativa * 17200,00 70570,98

Diferença central * * 192,15

* O algoritmo divergiu.

Comparando-se a média de iterações por intervalo de tempo, observa-se na

Figura 4.16 e na Tabela 4.8 que esse número cai para os algoritmos NR - Método direto,

NR - PCG e NR - CG. Porém, para os algoritmos NLCG - Busca linear analítica e

NLCG - Busca linear iterativa, a média permanece alta.

123

Variação do intervalo de tempo

Número de iterações por intervalo de tempo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.16 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado.

Tabela 4.12 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado.

ALGORITMO UTILIZADO

∆

t∆

NR - Método direto * 2,00 1,02

NR - PCG 13,90 2,00 1,02

NR - CG 14,02 2,00 1,03

NLCG - Busca linear analítica * 15,00 15,00

NLCG - Busca linear iterativa * 15,00 15,00

* O algoritmo divergiu.

4.4. Cabo tracionado

Este exemplo trata do caso de um cabo, inicialmente tracionado e na posição

horizontal, tendo as extremidades esquerda e direita presas e seu peso próprio aplicado

gradualmente de forma que o carregamento seja quase estático. A matriz de massa

utilizada é diagonal, possibilitando comparações das respostas obtidas pelos algoritmos

implícitos e o explícito (método da diferença central). O objetivo deste exemplo é

verificar a performance desses métodos na obtenção da configuração estática não linear

de cabos para 3 níveis de discretização da malha.

A Figura 4.17 apresenta a malha de elementos finitos utilizada e a Tabela 4.13,

as propriedades físicas e geométricas do cabo.

Figura 4.17 – Desenho esquemático da malha de elementos finitos.

Tabela 4.13 – Propriedades físicas e geométricas do cabo.

Comprimento do cabo 500 m

Diâmetro* 0,214 m

Área da seção transversal (circular)*

⋅π

Peso específico linear* 3,1977 kN/m

Peso específico linear submerso* 2,7742 kN/m

Módulo de elasticidade longitudinal*

10907,3 ⋅

kN/m

Tensão inicial 10

kN/m

* Dados obtidos em SILVEIRA (2001).

Observa-se que neste exemplo as propriedades físicas e geométricas utilizadas se

referem a um tipo de linha de ancoragem utilizada nos exemplos de SILVEIRA (2001).

Por se tratar de um caso de cabo que geralmente está submerso quando solicitado,

utilizou-se o valor do seu peso submerso, ao invés do seu peso ao ar.

4.4.1.

Obtenção da configuração estática

Nesta análise utiliza-se a malha de 12 elementos finitos de treliça, cuja

configuração inicial está apresentada na Figura 4.17. O peso próprio foi aplicado

linearmente até um dado instante de tempo, a partir do qual seu valor permanece

constante. Inicialmente, avalia-se qual o tempo necessário de aplicação do carregamento

para que se obtenha uma configuração quase estática. A menor freqüência da estrutura

tem valor de 0,2082 rad/s, levando a um período de 30,2s. Adota-se um tempo de 1000s

para aplicação do peso próprio, bem superior ao maior período natural da estrutura, e

mais 200s para observação da configuração da estrutura sem aumento do carregamento,

totalizando 1200s de análise. A Figura 4.10 ilustra a aplicação do peso próprio.

c d

...

Figura 4.18 – Curva de aplicação do peso próprio.

Para execução do algoritmo da Diferença Central, é necessária uma avaliação das

máximas freqüências do modelo. A maior freqüência registrada é de 9,88 rad/s, levando

a um valor crítico de intervalo de tempo de 0,2024s. Opta-se por executar uma análise

de referência, além das análises com cada algoritmo, com um intervalo de tempo no

valor de 0,01s. Devido ao pequeno valor do intervalo de tempo para análise de

referência, utiliza-se o Método da Diferença Central na sua execução.

Os valores de intervalo de tempo utilizados em cada algoritmo estão

apresentados na Tabela 4.14.

Tabela 4.14 – Valores dos intervalos de tempo.

ALGORITMO UTILIZADO

∆

(s)

NR – Método direto 8,0

NR – PCG 8,0

NR – CG 8,0

NLCG – Busca linear analítica 8,0

NLCG – Busca linear iterativa 8,0

Diferença Central 0,2

Diferença Central - referência 0,01

A história de deslocamentos verticais do nó central, obtida pelos algoritmos

estudados, está apresentada na Figura 4.19.

sub

W(t)

1000s

1200s

0 200 400 600 800 1000 1200

-80

-70

-60

-50

-40

-30

-20

-10

Tempo (s)

Deslocamento vertical no nó central (m)

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.19 – História de deslocamentos verticais do nó central.

Para melhor visualização dos resultados, apresenta-se na Figura 4.20, um detalhe

dos deslocamentos em um intervalo de tempo após a aplicação total do peso próprio.

1000 1010 1020 1030 1040 1050 1060

-77

-76.5

-76

-75.5

-75

-74.5

-74

-73.5

-73

Tempo (s)

Deslocamento vertical no nó central (m)

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.20 – Detalhe dos deslocamentos verticais do nó central.

Área do detalhe

Apesar do carregamento ser aplicado lentamente, a estrutura apresenta oscilações

ao final da análise (Figura 4.20). Este fato é decorrente da falta de amortecimento no

modelo estudado.

Observa-se (Figura 4.20) que os algoritmos implícitos não integram

corretamente as maiores freqüências da estrutura, porém, obtêm resposta satisfatória

considerando que o objetivo é encontrar uma configuração estática para o cabo.

Apresenta-se na Figura 4.21 a deformada do cabo para todas as análises

realizadas.

0 50 100 150 200 250 300 350 400 450 500

-200

-150

-100

-50

100

150

ormada do cabo para pequenos deslocamentos.

metros

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

Figura 4.21 – Deformada do cabo no tempo final de análise.

A proximidade das soluções dificulta a visualização de cada uma

individualmente. Portanto, apresenta-se na Figura 4.22 um detalhe da deformada no nó

central, onde pode-se distinguir as soluções dos algoritmos NR - CG, NLCG - Busca

linear iterativa, NLCG - Busca linear analítica e Diferença Central, além da análise de

referência. As soluções dos algoritmos NR - Método direto e NR - PCG estão

encobertas pela solução do algoritmo NR – CG. Observa-se que os deslocamentos estão

em torno de 75,5m enquanto que as diferenças entre as deformadas são menores que

50cm, ou seja, menores que 0,66% do deslocamento total nesse nó.

249.6 249.8 250 250.2 250.4 250.6 250.8

-76

-75.9

-75.8

-75.7

-75.6

-75.5

-75.4

-75.3

-75.2

-75.1

-75

Deformada do cabo para pequenos deslocamentos.

metros

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

Figura 4.22 – Detalhe da deformada no nó central.

A Figura 4.23 apresenta algumas deformadas do decorrer da análise e a

deformada final.

Figura 4.23 – Etapas da análise.

Os tempos de simulação e as médias de iterações por intervalo de tempo de cada

algoritmo estão apresentados na Figura 4.24.

(a) (b)

100

110

120

Tempo de simulação (s)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

783,77

Média de iterações por intervalo de tempo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.24 – (a) Tempos de simulação; (b) Média de iterações por intervalo de tempo.

Observa-se que os algoritmos implícitos que usam o método de Newton-

Raphson para solução do sistema não linear apresentam tempo de simulação bem

inferior que os demais algoritmos, incluindo-se o Método da Diferença Central. Este

último só foi ultrapassado pelo algoritmo NLCG - Busca linear iterativa.

Observando-se a Figura 4.24 (b), a média de iterações por intervalo de tempo dos

algoritmos que usam Gradientes Conjugados Não Linear apresenta-se bem superior aos

demais algoritmos implícitos, justificando o maior tempo gasto por esses métodos.

Os valores numéricos dos tempos de simulação e médias de iterações,

representados graficamente na Figura 4.24, estão apresentados na Tabela 4.15.

Tabela 4.15 – Tempos de simulação e médias de iterações por intervalo de tempo.

ALGORITMO UTILIZADO TEMPOS (S) MÉDIAS DAS ITERAÇÕES

NR – Método direto 16,46 2,61

NR – PCG 18,89 2,61

NR – CG 18,46 2,61

NLCG – Busca linear analítica 75,61 14,91

NLCG – Busca linear iterativa 783,77 14,91

Diferença Central 94,58 -

4.4.2.

Influência do refinamento da malha de elementos finitos

Após analisar uma malha com 12 elementos, avalia-se qual o desempenho dos

algoritmos NR - PCG e NR - CG, em relação ao NR - Método direto, diante do aumento

do número de elementos do modelo. Comparam-se os resultados para as malhas de 12,

120 e 480 elementos.

Com a alteração da malha as freqüências máximas do modelo foram alteradas.

COOK e outros (1989) apresentam uma forma de cálculo estimado da máxima

freqüência da estrutura, baseado na velocidade de propagação da onda no elemento

linear, como apresentado na equação:

max

=ω ,

( 4.1 )

onde

( 4.2 )

sendo c a velocidade de propagação da onda, L o comprimento do elemento, E o

módulo de elasticidade longitudinal e ρ a densidade do material. COOK e outros (1989)

mostram que a máxima freqüência para um modelo formado por elementos finitos é

limitada pela máxima freqüência individual dos elementos que a constituem.

Apresentam-se na Tabela 4.16 as freqüências estimadas da forma acima citada, além

dos valores obtidos através do cálculo numérico das freqüências naturais. Observam-se

poucas diferenças entre as duas formas utilizadas.

Tabela 4.16 – Freqüências máximas.

Método 12 elem. 120 elem. 480 elem.

Autovalores 9,88 99,66 398,66

Estimativa analítica 9,97 99,65 398,58

O refinamento da malha de elementos finitos gera mudanças na ordem de

grandeza das freqüências, o que altera o valor máximo para o incremento de tempo do

algoritmo da diferença central. A Tabela 4.17 apresenta os valores críticos para o

intervalo de tempo do Método da Diferença Central.

Tabela 4.17 – Valores críticos de intervalo de tempo.

12 elem. 120 elem. 480 elem.

t∆ (s)

2,02⋅10

-1

2,0⋅10

-2

5,02⋅10

-3

O valor do intervalo de tempo utilizado nos algoritmos implícitos para a malha

de 12 elementos permanece o mesmo, 8s. A Figura 4.25 apresenta os deslocamentos

verticais do nó central gerados pelo algoritmo NR - PCG para as três malhas avaliadas.

(a) (b)

0 200 400 600 800 1000 1200

-80

-70

-60

-50

-40

-30

-20

-10

Tempo (s)

Deslocamento vertical no nó central (m)

12 elementos

120 elementos

480 elementos

1000 1020 1040 1060 1080 1100 1120 1140

-76.5

-76

-75.5

-75

-74.5

12 elementos

120 elementos

480 elementos

Figura 4.25 – (a) Deslocamentos verticais do nó central para o algoritmo NR – PCG; (b)

Detalhe dos deslocamentos verticais do nó central.

A Tabela 4.18 apresenta os tempos de simulação e médias de iterações por

intervalo de tempo para as algoritmos iterativos.

Tabela 4.18 – Análise geral do desempenho dos algoritmos iterativos.

TEMPOS (S) MÉDIAS

ALGORITMO UTILIZADO

120 elem. 480 elem. 120 elem. 480 elem.

NR – Método direto 169,33 731,30 2,92 3,23

NR – PCG 234,74 1774,1 2,92 3,23

NR – CG 284,87 * 3,77 *

* O algoritmo divergiu.

JOUGLARD e COUTINHO (1998) classificam como problemas de grande porte

aqueles com mais de 100000 elementos. Assim sendo, apesar do aumento do sistema,

Área do detalhe

de 22 para 958 equações, ou de 12 para 480 elementos, o problema não é caracterizado

ainda como de grande porte, o que justifica o fato do método direto apresentar-se mais

eficiente que os iterativos neste exemplo.

Um outro fato que se observa é a divergência do algoritmo dos Gradientes

Conjugados sem precondicionamento. Apresenta-se na Tabela 4.19 o máximo número

de condição das matrizes efetivas para cada malha, indicando o grau de

condicionamento em cada caso.

Tabela 4.19 – Máximo número de condição da matriz efetiva.

12 elem. 120 elem. 480 elem.

cond

415,05 4204,0 670000,0

Quanto mais próximo da unidade estiver o valor do número de condição, melhor

é o seu condicionamento. Observa-se que o alto valor de

cond

N aliado à falta do uso do

precondicionador, para 480 elementos, levou à divergência do algoritmo NR - CG.

Nessas condições, o algoritmo dos Gradientes Conjugados Pré-Condicionados

apresenta-se mais eficiente e robusto que o método tradicional.

4.5. Pêndulos

Este é um exemplo no qual um cabo, cujo único carregamento é seu próprio

peso, é posicionado inicialmente em uma configuração parabólica tendo sua

extremidade direita solta no instante inicial da análise.

A matriz de massa utilizada é diagonal, fazendo o cabo comportar-se como um

pêndulo composto por n hastes. O objetivo deste exemplo é avaliar o comportamento

dos algoritmos estudados diante de problemas dinâmicos com elevado grau de não

linearidade.

A Figura 4.17 apresenta as malhas de elementos finitos utilizadas, com a

indicação da numeração dos elementos e dos graus de liberdade livres. As propriedades

físicas e geométricas do cabo estão apresentadas na Tabela 4.13, sendo neste caso o

valor do comprimento L igual a 10m.

(a) (b)

Figura 4.26 – Desenho esquemático da malha de elementos finitos.

4.5.1. Pêndulo Triplo

Inicialmente analisa-se a malha com três elementos finitos, utilizando-se o

algoritmo NR - Método direto para avaliar as máximas freqüências do sistema,

verificando-se qual o valor crítico de

∆

para o Método das Diferenças Centrais. O

histórico de máximas freqüências assim obtido está apresentado na Figura 4.27, para

s01,0t =∆ . A tolerância admitida nos algoritmos que utilizam um processo iterativo

para solução do sistema não linear é de 10

-6

0 5 10 15 20 25

100

102

104

106

108

Tempo (s)

Máximas freqüências (rad/s)

Figura 4.27 – História das máximas freqüências para o pêndulo triplo.

O maior valor entre as máximas freqüências foi de 107,2 rad/s, levando a um

valor máximo de 0,0187s para o incremento de tempo do Método das Diferenças

Centrais. Avaliam-se então os valores de intervalos de tempo iguais a 0,01s e 0,001s.

A fim de se avaliar o grau de não linearidade do problema, apresenta-se também,

na Figura 4.28 a variação das mínimas freqüências do sistema com o tempo.

0 5 10 15 20 25

0.5

1.5

2.5

Tempo (s)

Mínimas freqüências (rad/s)

Figura 4.28 – História das mínimas freqüências para o pêndulo triplo.

Observa-se grande variação de freqüências em curtos intervalos de tempo,

indicando que este é um problema com alto grau de não linearidade.

As respostas geradas com os valores de incremento de tempo, antes

apresentados, são comparadas com a solução obtida pelo MDC para um intervalo de

tempo de 0,0001s, considerada como referência.

A Figura 4.29 apresenta a resposta da análise realizada com o algoritmo NR -

Método direto para s01,0t =∆ em alguns instantes de tempo iniciais.

Figura 4.29 – Etapas da análise do pêndulo triplo.

Esta resposta pode ser melhor avaliada através da observação do histórico de

deslocamentos nos nós no modelo. Observa-se nas Figuras 4.28 e 4.29 os

deslocamentos horizontal e vertical do nó extremo da direita.

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tempo (s)

Deslocamento horizontal (m)

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.30 – Deslocamento horizontal na extremidade direita para s01,0t =∆ .

0 5 10 15 20 25

-12

-10

-8

-6

-4

-2

Tempo (s)

Deslocamento vertical (m)

Referência

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.31 – Deslocamento vertical na extremidade direita para s01,0t =∆ .

Observa-se a coerência entre os resultados, de forma que não se distingue

facilmente as respostas de cada algoritmo.

Executam-se novamente as simulações com os diversos algoritmos estudados,

porém com o segundo intervalo de tempo, de valor 0,001s. Comparando-se os

resultados do algoritmo NR - CG para os dois valores de intervalo de tempo, observa-se

na Figura 4.32 e no detalhe apresentado na Figura 4.33 que a redução do intervalo de

tempo aproximou a resposta obtida da solução de referência.

Figura 4.32 – Comparação entre as respostas obtidas com o algoritmo NR - CG para

dois intervalos de tempo.

23.7 23.75 23.8 23.85 23.9 23.95 24

-19.62

-19.6

-19.58

-19.56

-19.54

-19.52

-19.5

-19.48

Tempo (s)

Deslocamento horizontal (m)

Referência

dt=0,01s

dt=0,001s

Figura 4.33 – Detalhe da comparação entre as respostas obtidas com o algoritmo NR -

CG para dois intervalos de tempo.

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tem

(

)

Deslocamento horizontal (m)

Referência

dt=0,01s

dt=0,001s

Área do detalhe

A Figura 4.34 e a Tabela 4.5 apresentam os erros relativos para os algoritmos

estudados, para cada intervalo de tempo utilizado. Observa-se a razoável redução nos

erros com a diminuição do

t∆ , demonstrando a convergência dos vários algoritmos

para uma única resposta.

0,0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

1,1

1,2

1,3

1,4

1,5

1,6

1,7

1,8

Variação do intervalo de tempo

Erro relativo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

Figura 4.34 – Representação gráfica dos valores dos erros relativos em cada intervalo de

tempo analisado.

Tabela 4.20 – Valores dos erros relativos em cada intervalo de tempo analisado.

ERRO RELATIVO À SOLUÇÃO COM

0001,0t

∆

ALGORITMO UTILIZADO

∆

NR - Método direto 1,19 0,09

NR - PCG 1,19 0,09

NR - CG 1,19 0,09

NLCG - Busca linear analítica 1,19 0,09

NLCG - Busca linear iterativa 1,19 0,09

Diferença Central 1,70 0,05

Um importante fator observado é o tempo de análise gasto pelos algoritmos, que

são apresentados graficamente na Figura 4.35 e numericamente na Tabela 4.21.

Observa-se o natural aumento do tempo de análise com a redução do intervalo de

tempo, para cada algoritmo. Entre os algoritmos, o grupo que utiliza o método de

Newton-Raphson para solução do sistema não linear permanece mais eficiente que o

grupo que utiliza o Método dos Gradientes Conjugados Não Linear.

120

150

180

210

240

270

300

330

360

390

420

450

480

510

540

570

600

Variação do intervalo de tempo

Tempo de computação (s)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

2.011,71 1.700,98 12.029,2

Figura 4.35 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado.

Tabela 4.21 – Valores dos tempos de análise (s) para cada intervalo de tempo analisado.

ALGORITMO UTILIZADO

∆

NR - Método direto 38,61 395,13

NR - PCG 43,04 398,46

NR - CG 42,07 396,78

NLCG - Busca linear analítica 254,64 1700,98

NLCG - Busca linear iterativa 2011,71 12029,20

Diferença Central 15,15 148,74

- O algoritmo não foi executado.

A Figura 4.36 e a Tabela 4.22 apresentam esses resultados, onde observa-se a

redução da média com a diminuição do intervalo de tempo.

0,0

1,0

2,0

3,0

4,0

5,0

6,0

7,0

8,0

9,0

10,0

11,0

Variação do intervalo de tempo

Número de iterações por intervalo de tempo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.36 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado.

Tabela 4.22 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado.

ALGORITMO UTILIZADO

∆

NR - Método direto 1,08 1,00

NR - PCG 1,08 1,00

NR - CG 1,08 1,00

NLCG - Busca linear analítica 10,37 6,61

NLCG - Busca linear iterativa 9,30 5,97

- O algoritmo não foi executado.

Utilizando-se este modelo de cabo com três elementos finitos, e o algoritmo NR-

CG, realiza-se uma comparação entre a implementação do Método dos Gradientes

Conjugados que utiliza as matrizes globais e a implementação elemento por elemento

desse método.

A Tabela 4.23 apresenta uma comparação dos tempos de análise que cada

implementação utiliza.

Tabela 4.23 – Comparação entre implementações para o algoritmo NR-CG.

IMPLEMENTAÇÃO UTILIZADA TEMPO

NR – CG (global) 41,64

NR – CG (ExE) 56,58

Por se tratarem do mesmo algoritmo, as implementações apresentam soluções

idênticas, ou seja, o erro relativo entre elas é nulo. As diferenças surgem apenas na

quantidade de memória utilizada durante o processamento e no tempo de execução. A

implementação ExE executa

loops a mais do que a implementação que utiliza as

matrizes globais. Nesses

loops, necessários sempre que forem executadas operações do

tipo multiplicação vetor-matriz, realizam-se serialmente processos que poderiam ser

executados paralelamente.

4.5.2. Pêndulo Quíntuplo

Avalia-se também um cabo discretizado com cinco elementos, formando um

pêndulo quíntuplo. Neste caso, a estrutura tem as configurações apresentadas na Figura

4.37 nos instantes iniciais da simulação.

Figura 4.37 – Etapas da análise do pêndulo quíntuplo.

Na Figura 4.38 apresentam-se os deslocamentos vertical e horizontal na

extremidade direita do cabo, para os dois intervalos avaliados. Observa-se que, para

essa discretização e esses valores de incrementos de tempo, os algoritmos não

convergem para uma única solução, resultado diferente do obtido para três elementos.

Neste caso, admite-se a mesma tolerância (10

-6

) do caso com três elementos para os

algoritmos que utilizam um processo iterativo para solução do sistema não linear.

(a) (b)

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tempo (s)

Deslocamento horizontal (m)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

0 5 10 15 20 25

-12

-10

-8

-6

-4

-2

Tempo (s )

Deslocamento horizontal (m)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tempo (s )

Deslocamento horizontal (m)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

0 5 10 15 20 25

-12

-10

-8

-6

-4

-2

Tempo (s)

Deslocamento horizontal (m)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença Central

Figura 4.38 – (a) e (b) Deslocamentos horizontal e vertical na extremidade direita para

s01,0t =∆ ; (c) e (d) Deslocamentos horizontal e vertical na extremidade direita para

s001,0t

∆

Realiza-se um teste com os algoritmos NR - Método direto e Diferença central

com um intervalo de tempo menor que os dois anteriores e de valor s0001,0t =∆ . O

resultado obtido está apresentado na Figura 4.39.

(a) (b)

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tempo (s)

Deslocamento horizontal (m)

NR-Método direto

Diferença Central

0 5 10 15 20 25

-12

-10

-8

-6

-4

-2

Tempo (s )

Deslocamento vertical (m)

NR-Método direto

Diferença Central

Figura 4.39 – (a) e (b) Deslocamentos horizontal e vertical na extremidade direita para

s0001,0t

∆

Constata-se que, mesmo utilizando-se o intervalo de tempo pequeno, as respostas

continuam apresentando discordância, principalmente nos últimos 10s de simulação.

Observando o histórico de máximas freqüências para esta última discretização,

Figura 4.40, tem-se que a maior freqüência registrada tem valor de 187,5 rad/s. Esse

histórico registra maiores variações das máximas freqüências e um aumento de cerca de

75% na maior freqüência em relação à malha de três elementos.

0 5 10 15 20 25

160

165

170

175

180

185

190

Tempo (s)

Máximas freqüências (rad/s)

Figura 4.40 – História no tempo da máximas freqüências.

Apresenta-se também o histórico de mínimas freqüências, onde observa-se o

aumento da variabilidade dos valores em relação ao caso do pêndulo triplo, indicando

um maior nível de não linearidade deste caso.

0 5 10 15 20 25

0.5

1.5

2.5

3.5

4.5

Tempo (s)

Mínimas freqüências (rad/s)

Figura 4.41 – História no tempo da mínimas freqüências.

Os tempos de análise neste segundo caso foram maiores que no primeiro e estão

apresentados na Figura 4.42 e na Tabela 4.24.

120

150

180

210

240

270

300

330

360

390

420

450

480

510

540

570

600

Variação do intervalo de tempo

Tempo de simulação (s)

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Diferença central

4.302,72

3.368,64

24.728,60

Figura 4.42 – Representação gráfica dos tempos de simulação para cada intervalo de

tempo analisado.

Tabela 4.24 – Valores dos tempos de análise para cada intervalo de tempo analisado.

ALGORITMO UTILIZADO

∆

NR – Método direto 85,75 563,47

NR – PCG 95,07 569,57

NR – CG 93,75 573,24

NLCG – Busca linear analítica 514,50 3368,64

NLCG – Busca linear iterativa 4302,72 24728,60

Diferença central 21,42 205,91

Observa-se também um aumento na média de iterações por intervalo de tempo

em comparação com a malha de três elementos. Os resultados estão apresentados na

Figura 4.43 e na Tabela 4.25.

0,0

1,5

3,0

4,5

6,0

7,5

9,0

10,5

12,0

13,5

15,0

Variação do intervalo de tempo

Número de iterações por intervalo de tempo

NR-Método direto

NR-PCG

NR-CG

NLCG-Busca linear analítica

NLCG-Busca linear iterativa

Figura 4.43 – Representação gráfica das médias de iterações por intervalo de tempo para

cada intervalo de tempo analisado.

Tabela 4.25 – Valores das médias de iterações por intervalo de tempo para cada

intervalo de tempo analisado.

ALGORITMO UTILIZADO

∆

NR - Método direto 1,76 1,00

NR - PCG 1,75 1,00

NR - CG 1,76 1,00

NLCG - Busca linear analítica 13,71 8,63

NLCG - Busca linear iterativa 13,03 8,00

Com relação às medidas de erros, para essa malha, não é possível encontrar uma

configuração de referência, pois mesmo usando s0001,0t

∆

, as respostas dos

algoritmos NR - Método direto e Diferença central não geram respostas coincidentes

(Figura 4.39). Como os valores dos intervalos de tempo são bastante pequenos, avaliam-

se os erros do algoritmo explícito utilizado, Diferença Central, através da medida do

balanço energético. Esses valores estão apresentados na Figura 4.44.

(a) (b)

7.4 7.6 7.8 8 8.2 8.4 8.6 8.8

-2

-1

x 10

-3

Tempo

Resíduo - Energia

dt=1e-2s

0 5 10 15 20 25

0.2

0.4

0.6

0.8

1.2

x 10

-3

Tempo

Resíduo - Energia

dt=1e-3s

7.4 7.6 7.8 8 8.2 8.4 8.6 8.8

-2

-1

x 10

-5

Tempo

Resíduo - Energia

dt=1e-3s

(e) (f)

0 5 10 15 20 25

0.2

0.4

0.6

0.8

1.2

x 10

-5

Tempo (s)

Resíduo - Energia

dt=1e-4s

7.4 7.6 7.8 8 8.2 8.4 8.6 8.8

-2

-1

x 10

-7

Tempo (s)

Resíduo - Energia

dt=1e-4s

Figura 4.44 – (a), (c) e (e) Erros em energia para s01,0t

∆

,s001,0t

∆

e s0001,0t =∆ ;

(b), (d) e (f) Detalhes dos gráficos (a) (c) e (e).

Observa-se que no caso de s01,0t

∆

, após o oitavo segundo o erro passa a

aumentar chegando a ultrapassar o limite de 0,05 citado por COOK e outros (1989)

como um limite para estabilidade. O mesmo ocorre para s001,0t

∆

e s0001,0t =∆ ,

0 5 10 15 20 2

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

Tempo

íd

uo -

nerg

dt=1e-2s

Limite de 0,05

Área do detalhe

porém a ordem de grandeza do resíduo para esses casos é duas e quatro vezes

respectivamente menor, não atingindo o limite de erro de 0,05.

O fato de não se obter uma concordância entre as respostas dos diversos

algoritmos, para o pêndulo quíntuplo, demonstra que este exemplo possui características

particulares de não linearidades que devem ser analisadas. Uma possibilidade levantada

é este caso apresentar natureza caótica, caracterizando-se pela alta sensibilidade da

resposta obtida para mínimas diferenças nas condições iniciais do problema. Essa

possibilidade é investigada ao se comparar as respostas obtidas para um mesmo

algoritmo, porém com pequenas perturbações nas condições iniciais. As Figuras 4.45 e

4.46 apresentam esta comparação para o caso do algoritmo NR - Método direto, onde

observa-se o mesmo fenômeno de discordância de respostas apresentado quando

comparados os diferentes algoritmos estudados para condições iniciais idênticas. As

condições iniciais utilizadas nessa comparação são apresentadas na Tabela 4.26.

Este trabalho, porém, objetiva o estudo dos algoritmos, não sendo exploradas

essas características, ficando esse aspecto como sugestão para outros trabalhos que

abordem temas afins.

Tabela 4.26 – Valores das médias de iterações por intervalo de tempo para cada.

CONDIÇÕES INICIAIS NO GRAU DE LIBERDADE

VERTICAL NA EXTREMIDADE DIREITA

= 0 m v

= 0

= 1·10

-4

m v

= 0

= 1·10

-6

m v

= 0

0 5 10 15 20 25

-25

-20

-15

-10

-5

Tempo (s)

Deslocamento horizontal (m)

= 0

= 1e-4

= 1e-6

Figura 4.45 – Deslocamento horizontal na extremidade direita para s01,0t =∆ , usando o

algoritmo NR - Método direto.

0 5 10 15 20 25

-12

-10

-8

-6

-4

-2

Tempo (s)

Deslocamento vertical (m)

= 0

= 1e-4

= 1e-6

Figura 4.46 – Deslocamento vertical na extremidade direita para s01,0t =

∆

, usando o

algoritmo NR - Método direto.

5. Conclusões

5.1. Conclusões, comentários e sugestões

Este trabalho tem como objetivo principal realizar um estudo acerca de

algoritmos não lineares de integração no tempo e estratégias elemento por elemento

para esses algoritmos, implementando-os e avaliando seus comportamentos diante de

casos de análises dinâmicas e quase estáticas de estruturas. No decorrer deste trabalho,

apresentam-se resultados demonstrando que esse objetivo é alcançado.

A principal contribuição deste trabalho está na geração de conhecimento acerca

desses algoritmos e na sua disponibilização, estando preparados para implementações

futuras em ambientes de computação paralela. Os algoritmos implementados

apresentam, nos exemplos, comportamento compatível com o descrito pela literatura.

Os exemplos iniciam-se com a apresentação de um caso linear simples, sistema

massa-mola, onde as implementações demonstram estarem aferidas para os casos

lineares. Em seguida, aborda-se uma simulação não linear de uma viga engastada, cuja

modelagem utiliza um elemento de pórtico tridimensional. A matriz de massa aqui

utilizada é diagonal, permitindo comparações com os resultados obtidos entre

algoritmos implícitos e explícitos. Neste exemplo tem-se um caso de divergência de

algoritmos implícitos, fato que pode ocorrer por se tratar de analises não lineares. As

respostas geradas apresentam-se coerentes com o problema avaliado.

O terceiro exemplo apresenta um cabo tracionado, modelado por elementos de

barra, com matriz de massa concentrada, sujeito a seu peso próprio. Neste exemplo o

carregamento é aplicado gradualmente, tratando-se de uma análise quase estática.

Constata-se que, para essas classes de análise, os algoritmos implícitos obtém boa

performance em relação aos explícitos. Três discretizações são analisadas, obtendo-se

sucesso na definição da configuração estática do cabo. Observa-se que, apesar do

aumento do número de equações com o refinamento da malha até 480 elementos, o

método direto apresenta-se mais eficiente que os iterativos, neste exemplo. Esse

resultado é justificado pelo fato desse problema pertencer ainda à classe dos problemas

de pequeno porte, onde os métodos diretos são mais eficientes. Uma outra característica

observada nesse exemplo é a influência do pré-condicionamento na estabilidade do

Método dos Gradientes Conjugados, pois ao se analisar a malha com 480 elementos o

algoritmo sem precondicionador apresenta divergência. Isso mostra a necessidade da

utilização de precondicionadores nas futuras implementações dos algoritmos elemento

por elemento.

O último exemplo apresentado consta de uma análise de um cabo que pendula

durante alguns segundos. Esta é uma análise dinâmica com alta não linearidade sendo

apresentada para dois níveis de discretização.

Usando-se três elementos, os algoritmos apresentam o comportamento esperado,

convergindo para uma mesma solução à medida que o incremento de tempo era

reduzido. Usando este exemplo, comparam-se duas implementações do algoritmo NR-

CG, uma que monta as matrizes globais do sistema e outra elemento por elemento.

Constata-se que, apesar de serem o mesmo algoritmo, a implementação seqüencial do

algoritmo elemento por elemento no programa MATLAB apresenta maior custo

computacional. Nessa implementação necessita-se de mais estruturas do tipo

loops do

que quando se utiliza a matriz global. Como todas as implementações são realizadas no

programa MATLAB, é importante observar que esse programa apresenta melhor

performance para operações matriciais do que para estruturas de repetição. Ao se

substituir multiplicações de matrizes globais por uma série de multiplicações de

pequenas matrizes (dos elementos), a performance do programa é comprometida.

Sugere-se testar a eficiência dessa implementação em um ambiente de clusters de

computadores, realizando paralelamente essas operações elemento por elemento e

criação de uma classe de algoritmos de integração para dinâmica de estruturas, contendo

as implementações tradicionais e elemento por elemento, nas versões seqüencial e

paralela.

No segundo nível de discretização, utilizam-se cinco elementos, formando-se um

pêndulo quíntuplo. Assim como no primeiro caso, a matriz de massa utilizada foi

diagonal. São testados três valores para intervalos de tempo, porém, em todos os casos,

as respostas geradas por cada algoritmo, a partir de determinados instantes de tempo,

apresentam-se diferentes entre si. A dificuldade em encontrar uma única solução indica

que esse problema não linear possui características de um sistema caótico, gerando-se

diferentes soluções a partir de determinado instante em virtude do diferente acúmulo de

erros de cada algoritmo.

Por fim, sugere-se o estudo de problemas com outros tipos de não linearidades

(fisíca, contato, atrito) e restrições cinemáticas generalizadas, associados a algoritmos

elemento por elemento.

Referências bibliográficas

BATHE, K.J. (1996) Finite Element Procedures. Prentice-Hall International Editions.

BELYTSCHKO, T., (1984) An Improved Element-by-Element Semi-Implicit Scheme

for Dynamic Problems. Nuclear Engineering and Design, 80, 127-139.

BELYTSCHKO, T. e HUGHES, T., (1983) Computational Methods in Mechanics –

Computational Methods for Transient Analysis, V-1, Elsevier Science Publishers B. V.

BOLDRINI, J.L.; COSTA, S.I.R.; FIGUEIREDO, V.L. e WETZLER, H.G. (1980)

Álgebra Linear, Terceira Edição, Harper & Row do Brasil.

COOK, R.D.; MALKUS, D.S. e PLESHA, M.E. (1989) Concepts and Applications of

Finite Element Analysis, Third Edition, Jonh Wiley & Sons.

FLETCHER, R. e REEVES, C.M. (1964) Function minimization by conjugate

gradients. Computer J., Vol 6, 149-154.

HUGHES, T.J.R. (1987) The Finite Element Method – Linear static and dynamic finite

element analysis. Prentice-Hall International Editions.

HUGHES, T.J.R. e BELYTSCHKO, T. (1983) A Précis of Developments in

Computational Methods for Transient Analysis. Journal of Applied Mechanics, Vol. 50,

1033-1041.

HUGHES, T.J.R., LEVIT, I. e WINGET, J. (1983) An Element-by-Element Solution

Algorithm for Problems of Structural and Solid Mechanics. Computer Methods in

Applied Mechanics and Engineering, 36, 241-254.

JOUGLARD, C.E. e COUTINHO, A.L.G.A. (1998) A comparison of iterative multi-

level finite element solvers, Computers & Structures, 69, 655-670.

KRYSL, P. e BELYTSCHKO, T. (1998) Object-oriented parallelization of explicit

structural dynamics with PVM, Computers & Structures, 66, 259-273.

MATLAB (2000) User’s Guide. The MathWorks Inc., Natick, MA, USA.

NEWMARK, N. M. (1959) A Method of Computation for Structural Dynamics, ASCE

Journal of Engineering Mechanics Division, Vol. 85, 67-94.

ORTIZ, M., PINSKY, P.M. e TAYLOR, R.L. (1983) Unconditionally Stable Element-

by-Element Algorithms for Dynamics Problems. Computer Methods in Applied

Mechanics and Engineering, 36, 223-239.

PACOSTE, C. e ERIKSSON, A. (1997) Beam elements in instability problems.

Computer Methods in Applied Mechanics and Engineering, 144, 163-197.

PAPADRAKAKIS, M. e DRACOPOULOS, M.C. (1991) Improving the Efficiency of

Preconditioning for Iterative Methods. Computers & Structures, 41, 1263-1272.

PAPADRAKAKIS, M. e GHIONIS, P. (1986) Conjugate Gradient Algorithms in

Nonlinear Structural Analysis Problems. Computer Methods in Applied Mechanics and

Engineering, 59, 11-27.

SHEWCHUK, J.R. (1994) Na Introduction to the Conjugate Gradient Method Without

the Agonizing Pain. Edition 1¼, School of Computer Science Carnegie Mellon

University Pittsburgh, PA 15213.

SILVEIRA, E.S.S. (2001) Análise Dinâmica de Linhas de Ancoragem com Adaptação

no Tempo e Subciclagem. Tese de Doutorado , PUC/Rio de Janeiro, Programa de Pós-

Graduação em Engenharia Civil / Estruturas.

SORENSON, H.W. (1969) Comparison of some conjugate directions procedures for

function minimization, J. Franklin Inst. 288 421-442.

SUBBARAJ, K. e DOKAINISH, M.A. (1989) A Survey of Direct Time-Integration

Methods in Computational Structural Dynamics – II. Implicit Methods. Computers &

Structures, Vol. 32 (6), 1387-1401.

WATHEN, A.J., (1989) An Analysis of Some Element-by-Element Techniques.

Computer Methods in Applied Mechanics and Engineering, 74, 271-287.

WEAVER, W.Jr. e JOHNSTON, P.R. (1987) Structural Dynamics by Finite Elements,

Prentice-Hall, Inc.

WINGET, J.M. e HUGHES, T.J.R. (1985) Solution Algorithms for Nonlinear Transient

Heat Conduction Analysis Employing Element-by-Element Iterative Strategies.

Computer Methods in Applied Mechanics and Engineering, 52, 711-815.

ZIENKIEWICZ, O.C e MORGAN, K. (1983) Finite Elements and Aproximations, Jonh

Wiley & Sons.

Apêndice A

A Exemplo numérico de multiplicação de matriz-

vetor e vetor-matriz-vetor elemento por

elemento

Para facilitar a compreensão das operações matriciais elemento por elemento

exploradas neste trabalho, apresenta-se neste apêndice um exemplo de uma estrutura

composta por barras de rigidez

k , com o objetivo de realizar as suas operações matriz-

vetor e vetor-matriz-vetor explicitamente.

A.1 Descrição do exemplo

A estrutura utilizada como exemplo está apresentada na Figura A. 1. Trata-se de

três molas conectadas serialmente, sendo a primeira presa na extremidade esquerda.

Figura A. 1 – Estrutura de molas conectadas em série.

Os elementos, suas rigidezes e deslocamentos apresentam-se na Figura A. 2.

Figura A. 2 – Numeração dos elementos.

A matriz de rigidez de cada elemento é dada por:



















−

k ,

( A. 1 )

onde

k é o valor numérico da cada rigidez e está apresentado na Tabela A. 1 que

segue.

Tabela A. 1 – Valores das rigidezes dos elementos.

Em cada nó de cada elemento, aplica-se as seguintes cargas nodais:













F ,

( A. 2 )













−

F e

( A. 3 )













F .

( A. 4 )

Para mapear os graus de liberdade locais nos graus de liberdade globais, utiliza-

se as seguintes matrizes de incidência para cada elemento:













001

000

L ,

( A. 5 )













010

001

L e

( A. 6 )













100

010

L .

( A. 7 )

A matriz global do sistema é dada segundo uma montagem das matrizes dos

elementos mapeadas para as coordenadas globais. Tem-se:

∑

nelem

LkLK

, ou seja

( A. 8 )













−

−−

−

440

462

025

( A. 9 )

Da mesma forma, o vetor de forças externas é montado a partir da contribuição

de cada elemento:

∑

nelem

FLF , ou seja

( A. 10 )













−=

( A. 11 )

O vetor deslocamento solução desse sistema é dado por:

FKu

1−

= ou

( A. 12 )













58,1

83,0

33,1

( A. 13 )

A.1.1 Multiplicação matriz-vetor elemento por elemento

As parcelas envolvidas na multiplicação elemento por elemento apresentada são

o deslocamento e a matriz de rigidez. Portanto é preciso montar o vetor dos

deslocamentos correspondentes a cada elemento. Para tal objetivo, definem-se matrizes

booleanas que fazem o mapeamento inverso ao apresentado anteriormente, ou seja,

matrizes que mapeiam os graus de liberdade globais nos graus de liberdade locais. Tem-

se:













( A. 14 )













( A. 15 )













Li .

( A. 16 )

A fim de utilizar as matrizes e vetores dos elementos mapeadas tanto nas

coordenadas locais como nas globais, montam-se os vetores de deslocamentos dos

elementos nos dois sistemas coordenadas.

Utilizando-se das matrizes de conectividade

Li , obtêm-se os vetores de

deslocamentos nodais dos elementos, nas coordenadas locais:













33,1

uLiu ,

( A. 17 )













83,0

33,1

uLiu e

( A. 18 )













58,1

83,0

uLiu .

( A. 19 )

Mapeando-se esses vetores para as coordenadas globais, tem-se:













33,1

uLu

( A. 20 )













83,0

33,1

uLiu

( A. 21 )













58,1

83,0

uLiu .

( A. 22 )

A primeira alternativa apresentada é realizar o somatório utilizando as matrizes e

vetores nas coordenadas globais:

FuKKu =













−==

∑

( A. 23 )

A segunda alternativa é realizar a operação nas coordenadas locais, necessitando-

se fazer o mapeamento em cada parcela do somatório:

FukLKu =













−==

∑

( A. 24 )

Observa-se que neste último caso não é necessário o armazenamento das

matrizes, nem dos vetores de deslocamentos dos elementos, nas coordenadas globais.

A.1.2 Multiplicação vetor-matriz-vetor elemento por elemento

Uma outra operação passível de ser realizada elemento por elemento é a

multiplicação vetor-matriz-vetor, realizada durante as iterações do Método dos

Gradientes Conjugados entre a direção de busca, vetor

, e a hessiana, matriz

Neste exemplo, apenas para ilustrar o procedimento da operação, são

multiplicados o vetor

u e a matriz K . O resultado esperado é então:

083,8=uKu

( A. 25 )

Assim como na operação de multiplicação vetor-matriz, apresenta-se primeiro o

caso em que o vetor e a matriz estão mapeados para as coordenadas globais:

083,825,25,0333,5

ˆˆ

eee

=++==

∑

uKuuKu .

( A. 26 )

Neste caso, as parcelas do somatório são escalares, portanto, para utilizar o vetor

e a matriz de rigidez dos elementos nas coordenadas locais, não é necessário fazer o

mapeamento:

083,825,25,0333,5

eee

=++==

∑

ukuuKu

( A. 27 )

Apêndice B

B Formulação utilizada do elemento de barra

Neste apêndice apresenta-se sucintamente a formulação das matrizes de massa e

rigidez e vetor de forças internas do elemento de barra utilizado nos exemplos que

envolvem simulações com cabos. Maiores informações podem ser encontradas em

literaturas tradicionais que tratam do Método dos Elementos Finitos, como por exemplo

COOK e outros (1989).

B.1 Elemento de barra

O elemento finito utilizado é um elemento linear de barra (treliça) com quatro

graus de liberdade de translação, conforme exposto na Figura A. 1. Este é um elemento

simples que apresenta apenas esforço normal.

Figura B. 1 – Elemento de barra.

Para a geração da matriz de massa, opta-se por utilizar a matriz concentrada, na

qual a massa do elemento é considerada concentrada em partículas nos nós. Sendo

assim, a matriz de massa do elemento é diagonal e dada por:

L, A,

, E













⋅⋅=

1000

0100

0010

0001

M ,

( B. 1 )

onde m é a massa total do elemento.

Para obtenção da matriz de rigidez e vetor de forças internas utiliza-se o

software

Maple, gerando-se a seguinte seqüência de cálculos:

> restart:

> with(linalg): with(codegen):

Warning, the protected names norm and trace have been redefined and

unprotected

Warning, the protected name MathML has been redefined and

unprotected

Comprimento do elemento na configuração final:

> Le := sqrt(((xj+u3)-(xi+u1))^2+((yj+u4)-(yi+u2))^2):

Medida de deformação de engenharia:

> eps := (Le-Lo)/Lo:

Energia potencial elástica:

> U := int((sigmao*(eps)+E*((eps)^2)/2)*A,L=0..Lo);

U sigmao sqrt xj

2 xj u3 2 xj ξ 2 xj u1 u3

2 u3 ξ 2 u3 u1 ξ

2 ξ u1 + − − + − − + + ((







2 yj u4 2 yj yi 2 yj u2 u4

2 u4 yi 2 u4 u2 yi

2 yi u2 + + + − − + − − + +

+ ) Lo − ) Lo/

E sqrt xj

2 xj u3 2 xj ξ 2 xj u1 u3

2 u3 ξ 2 u3 u1 + − − + − − (( +

2 ξ u1 u1

2 yj u4 2 yj yi 2 yj u2 u4

2 u4 yi 2 u4 u2 yi

+ + + + +

−

− − +

2 yi u2 u2

+ + ) Lo − )







ALo

Vetor das forças internas:

> Fint := grad(U,[u1,u2,u3,u4]):

Matriz de rigidez:

> Kt := hessian(U,[u1,u2,u3,u4]):

Geração do código computacional:

> codegen[C](Fint,optimized):

t1 = xj*xj;

t8 = u3*u3;

t13 = xi*xi;

t16 = u1*u1;

t17 = yj*yj;

t24 = u4*u4;

t29 = yi*yi;

t32 = u2*u2;

t33 = t1+2.0*xj*u3-2.0*xj*xi-2.0*xj*u1+t8-2.0*u3*xi-2.0*u3*u1

+t13+2.0*xi*u1+t16+t17+2.0*yj*u4-2.0*yj*yi-2.0*yj*u2+t24-2.0*u4*yi-

2.0*u4*u2+t29+2.0*yi*u2+t32;

t34 = sqrt(t33);

t35 = 1/t34;

t36 = sigmao*t35;

t37 = -xj-u3+xi+u1;

t38 = 1/Lo;

t42 = E*(t34-Lo);

t43 = Lo*Lo;

t45 = 1/t43*t35;

t51 = -yj-u4+yi+u2;

Fint[0] = (2.0*t36*t37*t38+2.0*t42*t45*t37)*A*Lo/2.0;

Fint[1] = (2.0*t36*t51*t38+2.0*t42*t45*t51)*A*Lo/2.0;

Fint[2] = (-2.0*t36*t37*t38-2.0*t42*t45*t37)*A*Lo/2.0;

Fint[3] = (-2.0*t36*t51*t38-2.0*t42*t45*t51)*A*Lo/2.0;

> codegen[C](Kt,optimized);

t1 = xj*xj;

t8 = u3*u3;

t13 = xi*xi;

t16 = u1*u1;

t17 = yj*yj;

t24 = u4*u4;

t29 = yi*yi;

t32 = u2*u2;

t33 = t1+2.0*xj*u3-2.0*xj*xi-2.0*xj*u1+t8-2.0*u3*xi-2.0*u3*

u1+t13+2.0*xi*u1+t16+t17+2.0*yj*u4-2.0*yj*yi-2.0*yj*u2+t24-

2.0*u4*yi-2.0*u4*u2+t29+2.0*yi*u2+t32;

t34 = sqrt(t33);

t36 = 1/t34/t33;

t37 = sigmao*t36;

t38 = -xj-u3+xi+u1;

t39 = 4.0*t38*t38;

t40 = 1/Lo;

t44 = 1/t34;

t46 = sigmao*t44*t40;

t48 = E/t33;

t49 = Lo*Lo;

t50 = 1/t49;

t55 = E*(t34-Lo);

t56 = t50*t36;

t61 = t55*t50*t44;

t64 = (-t37*t39*t40/4.0+t46+t48*t39*t50/4.0-t55*t56*t39/4.0+

t61)*A*Lo;

t65 = -yj-u4+yi+u2;

t69 = 2.0*t38*t50;

t72 = t55*t50;

t78 = (-4.0*t37*t65*t40*t38+2.0*t48*t69*t65-

4.0*t72*t36*t65*t38) *A*Lo/4.0;

t79 = -2.0*t38*t40;

t86 = -2.0*t36*t38;

t92 = (-t37*t79*t38/2.0-t46-t48*t69*t38/2.0-t72*t86*t38/2.0-

t61)*A*Lo;

t93 = -2.0*t65*t40;

t98 = -2.0*t36*t65;

t103 = (-2.0*t37*t93*t38-2.0*t48*t69*t65-2.0*t72*t98*t38)*A*

Lo/4.0;

t104 = 4.0*t65*t65;

t116 = (-t37*t104*t40/4.0+t46+t48*t104*t50/4.0-

t55*t56*t104/4.0+t61)*A*Lo;

t119 = 2.0*t65*t50;

t126 = (-2.0*t37*t79*t65-2.0*t48*t119*t38-

2.0*t72*t86*t65)*A*Lo/4.0;

t138 = (-t37*t93*t65/2.0-t46-t48*t119*t65/2.0-t72*t98*t65/2.0-

t61)*A*Lo;

t148 = (2.0*t37*t93*t38+4.0*t48*t38*t50*t65+2.0*t72*t98*t38)

*A*Lo/4.0;

Kt[0][0] = t64;

Kt[0][1] = t78;

Kt[0][2] = t92;

Kt[0][3] = t103;

Kt[1][0] = t78;

Kt[1][1] = t116;

Kt[1][2] = t126;

Kt[1][3] = t138;

Kt[2][0] = t92;

Kt[2][1] = t126;

Kt[2][2] = t64;

Kt[2][3] = t148;

Kt[3][0] = t103;

Kt[3][1] = t138;

Kt[3][2] = t148;

Kt[3][3] = t116;

Figura B. 2 – Geração do vetor de forças internas de matriz de rigidez do elemento de

barra.

Apêndice C

C Solução estática para a viga engastada sujeita

a momento concentrado (exemplo 4.3)

Neste apêndice descreve-se a formulação para obtenção da solução estática do

caso abordado no exemplo do item 4.3.

C.1 Deslocamento vertical da extremidade direita

Como apresentado no exemplo do item 4.3, trata-se de uma viga engastada na

extremidade esquerda, sujeita a uma força momento aplicada na extremidade direita

(Figura C. 1).

Figura C. 1 – Desenho esquemático do modelo.

Para uma viga sujeita a flexão devido a uma força momento aplicada, sua

curvatura pode ser expressa como:

k == ,

( C. 1 )

sendo k a curvatura; R o raio de curvatura; M o momento fletor; E o módulo de

elasticidade e I o momento de inércia.

Observando-se a Figura C. 2, o deslocamento vertical (

vert

u ) na extremidade

direita é dado por:

E, I, L

100

−= cosRRu

vert

( C. 2 )

Figura C. 2 – Configuração genérica.

Sabe-se da geometria de figuras planas que a relação entre o raio, o ângulo e o

comprimento de um dado arco de círculo é dada por:

θ= RL .

( C. 3 )

Substituindo-se as equações ( C. 1 ) e ( C. 3 ) em ( C. 2 ) e simplificando-se, tem-

se que o deslocamento vertical é expresso por:

























−=

cos1

vert

( C. 4 )

Essa é a expressão utilizada no exemplo do item 4.3 para obtenção da resposta

estática apresentada graficamente.

C.2 Momento de referência

Para obtenção do momento de referência M

utilizado no exemplo do item 4.3,

considera-se para a viga a configuração apresentada na Figura C. 3.

vert

101

Figura C. 3 – Configuração de referência.

Nessa configuração o comprimento da viga é dado por:

R2L π= .

( C. 5 )

Assim, substituindo-se o comprimento expresso por ( C. 5 ) em ( C. 1 ), tem-se o

momento de referência:

EI2

= .

( C. 6 )

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 