CINTHIA ANDREIA GARCIA SOUSA
UM ESTUDO SOBRE MÉTODOS DE
CONTINUAÇÃO PARA ANÁLISE DE
ESTRUTURAS NÃO LINEARES
São Paulo
2012
CINTHIA ANDREIA GARCIA SOUSA
UM ESTUDO SOBRE MÉTODOS DE
CONTINUAÇÃO PARA ANÁLISE DE
ESTRUTURAS NÃO LINEARES
Dissertação de Mestrado apresentada à Escola Politécnica da Universidade de São Paulo para obtenção
de título de mestre.
Área de Concentração:
Engenharia de Estruturas
Orientador:
Prof. Dr. Paulo de Mattos Pimenta
São Paulo
2012
À minha mãe, meu irmão e ao meu noivo, Rodrigo.
ii
Agradecimentos
Ao orientador, Prof. Dr. Paulo de Mattos Pimenta, pelo apoio, incentivo, orientação e
sobretudo pela amizade e confiança em mim depositadas.
Ao Prof. Gabriel, por ter me apresentado a essa área fascinante, que é a Engenharia de
Estruturas. Agradeço imensamente por ter me apoiado e inspirado nas minhas decisões
profissionais.
À minha querida mãe, Maria das Graças, por ser meu grande exemplo de vida, amor,
força, esperança e fé. Ao meu irmão Demis, pela amizade e companheirismo. Agradeço
ainda por sempre acreditarem na minha capacidade e me darem suporte para superar
todas as etapas e dificuldades durante todo esse período.
Ao Rodrigo Rocha, meu grande amor e inspiração, por toda dedicação incondicional,
amor e carinho. Agradeço por sempre estar disposto a me ajudar e por tudo que eu
aprendi com ele até hoje. Obrigada por me dar o privilégio de viver ao seu lado.
A todos os colegas do JAC/LMC/AADP pela amizade. Em especial, ao Fernando Gonçalves, Leonardo Lago, Marcelo Teixeira, Jorge Costa, Henrique Campelo, Paulo Nigro
e a Gabriela Lins, pelo apoio, incentivo e companheirismo.
Às minhas grandes amigas Ludmila Figueiredo e Flávia Evangelista, por sempre estarem ao meu lado, pelo amor fraterno e amizade eterna. Ao meu grande amigo Ronaldo,
por sempre estar disposto a me ajudar e pela ótima companhia.
À FAPESP pelo apoio financeiro.
"O inimigo mais perigoso que você poderá encontrar será sempre você mesmo."
Friedrich Nietzsche.
iv
Resumo
Nas últimas décadas, considerável evidência tem sido direcionada no desenvolvimento de métodos computacionais que analisam os comportamentos não lineares das
estruturas. O método da corda é considerado o tipo de método de seguimento de trajetórias que possui a técnica mais comum e versátil para analisar tais comportamentos
não lineares. No entanto, testes realizados pelos autores concluem que tais métodos
ainda apresentam algumas dificuldades quando as estruturas possuem formas complexas e comportamentos fortemente não lineares. Devido a isso, o objetivo do trabalho em
questão é acrescentar duas modificações ao método da corda, a fim de desenvolver um
modelo mais seguro e estável. A esses métodos, inclui-se (i) uma equação alternativa
para o passo previsor, que tornará o método mais estável ao ultrapassar pontos críticos;
e (ii) um parâmetro de controle, que através do conceito de soma dos comprimentos de
corda, é chamado de método da corda acumulado. Ao final, por consequência destas
alterações, um procedimento mais robusto para o método da corda é obtido. Sua validação é apresentada, pelos autores, através de exemplos utilizando estruturas treliçadas
(bidimensionais e tridimensionais) em um programa de simulação numérica desenvolvida.
Palavras-chave: Métodos de Continuação, Análise Não Linear de Estruturas, Algoritmo(Otimização).
Abstract
In the last decades, a considerable effort has been directed to develop of computational methods to analyze the nonlinear behavior of structures. The arc-length method
is considered the type of path following method which have the most common and versatile techniques for analyzing such nonlinear behaviors. Nevertheless, tests performed
by the authors conclude that such methods present some difficulties when the structures
possess complex shapes and strongly nonlinear behaviors. Due to this, the purpose of
the present work is to add two amendments to the arc-length method, in order to develop a model safer and more stable. In such methods, it is included (i) an alternative
equation for the predictor step, that will make the method most stable to overcoming
critical points, and (ii) a control parameter, through the concept of sum of the arc length
up to now, referred to as accumulated arc-length method. By the end, as a result of these
modifications, a more robust procedure for the method of the string is obtained. And
its validation is presented by the authors through examples using truss structures in a
two-dimensional and three-dimensional numerical simulation program developed.
Keywords: Arc-Length Method, Nonlinear Structural Analysis, Algorithm (Optimization).
Sumário
Lista de Figuras
ix
Nomenclatura
xiii
1
Introdução
1
2
Teoria Geometricamente Exata para Treliças
6
2.1
Considerações Sobre Tensões e Deformações . . . . . . . . . . . . . .
6
2.2
Equações Constitutivas . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.3
Equações de Equilíbrio . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3
Métodos de Continuação
15
3.1
Formulação Inicial para os Métodos de Continuação . . . . . . . . . . .
21
3.2
Métodos Corretores . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.2.1
Método da Corda Esférico . . . . . . . . . . . . . . . . . . . .
24
3.2.2
Método da Corda Cilíndrico . . . . . . . . . . . . . . . . . . .
26
3.2.3
Método da Corda de Schweizerhof-Wriggers . . . . . . . . . .
30
Métodos Previsores . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
3.3.1
34
3.3
Sinal do Determinante da Matriz de Rigidez . . . . . . . . . . .
3.4
4
38
4.1
Método da Corda Acumulado . . . . . . . . . . . . . . . . . . . . . . .
40
4.1.1
O Comprimento da Corda Acumulado . . . . . . . . . . . . . .
41
4.1.2
Aproximação Quadrática . . . . . . . . . . . . . . . . . . . . .
43
4.1.3
Aproximação Linear . . . . . . . . . . . . . . . . . . . . . . .
45
4.1.4
Aproximação Quadrática Alternativa . . . . . . . . . . . . . .
46
4.1.5
Parâmetros de Controle . . . . . . . . . . . . . . . . . . . . . .
48
Novo Método Previsor . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.2.1
50
Critério do Produto Interno . . . . . . . . . . . . . . . . . . . .
Simulação Numérica
53
5.1
Métodos Previsores . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
5.1.1
Treliça de duas barras geometricamente simétricas . . . . . . .
53
5.1.2
Domo em forma de estrela . . . . . . . . . . . . . . . . . . . .
56
Métodos Corretores . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58
5.2.1
Domo em forma de estrela . . . . . . . . . . . . . . . . . . . .
59
5.2.2
Arco Circular bidimensional . . . . . . . . . . . . . . . . . . .
63
Método Acumulado . . . . . . . . . . . . . . . . . . . . . . . . . . . .
71
5.2
5.3
6
35
Otimização do Método da Corda
4.2
5
Comprimento de Corda Incremental . . . . . . . . . . . . . . . . . . .
Conclusões Finais
Referências
82
84
Lista de Figuras
2.1
Descrição da deformação de uma barra. . . . . . . . . . . . . . . . . .
6
3.1
Diagrama carga-deformação. . . . . . . . . . . . . . . . . . . . . . . .
16
3.2
Trajetórias primária e secundária. . . . . . . . . . . . . . . . . . . . . .
16
3.3
Comportamento de carga limite (snap-through). . . . . . . . . . . . . .
17
3.4
Comportamento de deslocamento limite (snap-back). . . . . . . . . . .
17
3.5
Método de Controle de Carga. . . . . . . . . . . . . . . . . . . . . . .
18
3.6
Método de Controle de Deslocamento. . . . . . . . . . . . . . . . . . .
18
3.7
Método da Corda proposto por Riks (1979). . . . . . . . . . . . . . . .
19
3.8
Método da Corda proposto por Crisfield (1991). . . . . . . . . . . . . .
20
3.9
Método da Corda de Schweizerhof-Wriggers (SCHWEIZERHOF; WRIGGERS, 1986).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
5.1
Treliça com duas barras simétricas e carga aplicada . . . . . . . . . . .
54
5.2
Diagrama de carga-deformação da treliça utilizando o previsor det(Kt ).
54
5.3
Diagrama de carga-deformação da treliça utilizando o Critério do Produto Interno. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4
55
Comparação dos Diagramas carga-deformação da treliça utilizando os
dois métodos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
5.5
Domo tridimensional em forma de estrela. . . . . . . . . . . . . . . . .
57
5.6
Diagrama de carga-deformação do Domo utilizando o previsor det(Kt ).
57
5.7
Diagrama de carga-deformação utilizando o Critério do Produto Interno.
5.8
Trajetória de equilíbrio com pontos convergidos do Domo utilizando o
5.9
58
Método da Corda Cilíndrico. . . . . . . . . . . . . . . . . . . . . . . .
59
Diagrama de convergência do Domo no Método da Corda Cilíndrico.
.
60
Método de Schweizerhof-Wriggers. . . . . . . . . . . . . . . . . . . .
60
5.10 Trajetória de equilíbrio com pontos convergidos do Domo utilizando o
5.11 Diagrama de convergência do Domo no Método de Schweizerhof-Wriggers. 61
5.12 Trajetória de equilíbrio do domo com pontos de deformação. . . . . . .
61
5.13 Estado deformado do domo no ponto A. . . . . . . . . . . . . . . . . .
62
5.14 Estado deformado do domo no ponto B. . . . . . . . . . . . . . . . . .
62
5.15 Estado deformado do domo no ponto C. . . . . . . . . . . . . . . . . .
62
5.16 Estado deformado do domo no ponto D. . . . . . . . . . . . . . . . . .
63
5.17 Estado deformado do domo no final da trajetória de equilíbrio. . . . . .
63
5.18 Arco circular treliçado. . . . . . . . . . . . . . . . . . . . . . . . . . .
64
5.19 Trajetória de equilíbrio do Arco utilizando o Método da Corda Cilíndrico. 64
5.20 Trajetória de equilíbrio do Arco, com pontos convergidos, utilizando o
Método da Corda Cilíndrico. . . . . . . . . . . . . . . . . . . . . . . .
65
5.21 Diagrama de convergência do Arco no Método da Corda Cilíndrico. . .
65
5.22 Diagrama de carga-deformação do Arco utilizando o Método de SchweizerhofWriggers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
5.23 Diagrama de carga-deformação do Arco, com pontos convergidos, utilizando o Método de Schweizerhof-Wriggers. . . . . . . . . . . . . . .
66
5.24 Diagrama de convergência do Arco no Método de Schweizerhof-Wriggers. 67
5.25 Trajetória de equilíbrio do Arco com pontos de deformação. . . . . . .
67
5.26 Estado deformado do arco no ponto A. . . . . . . . . . . . . . . . . . .
68
5.27 Estado deformado do arco no ponto B. . . . . . . . . . . . . . . . . . .
68
5.28 Estado deformado do arco no ponto C. . . . . . . . . . . . . . . . . . .
69
5.29 Estado deformado do arco no ponto D. . . . . . . . . . . . . . . . . . .
69
5.30 Estado deformado do arco no ponto E. . . . . . . . . . . . . . . . . . .
70
5.31 Estado deformado do arco no ponto F. . . . . . . . . . . . . . . . . . .
70
5.32 Estado deformado do arco no final do processo. . . . . . . . . . . . . .
71
5.33 Propriedades geométricas do cilindro (CAMPELLO, 2005). . . . . . . . .
72
5.34 Barras de treliça do cilindro. . . . . . . . . . . . . . . . . . . . . . . .
72
5.35 Diagrama de carga-deformação do cilindro utilizando o Método da Corda
Cilíndrico. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
5.36 Diagrama de carga-deformação do cilindro utilizando o Método da Corda
Acumulado (λd = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . .
74
5.37 Diagrama de comparação de convergência. . . . . . . . . . . . . . . . .
74
5.38 Diagrama de carga-deformação do cilindro para λd = 100. . . . . . . .
75
5.39 Configuração deformada do cilindro em λd = 100. . . . . . . . . . . .
76
5.40 Diagrama de carga-deformação do cilindro para λd = 1000. . . . . . . .
76
5.41 Configuração deformada do cilindro em λd = 1000. . . . . . . . . . . .
77
5.42 Diagrama de carga-deformação do cilindro para λd = 5000. . . . . . . .
77
5.43 Configuração deformada do cilindro em λd = 5000. . . . . . . . . . . .
78
5.44 Diagrama de carga-deformação do cilindro para λd = 10000. . . . . . .
78
5.45 Configuração deformada do cilindro em λd = 10000. . . . . . . . . . .
79
5.46 Configuração atual e deformada do cilindro. . . . . . . . . . . . . . . .
80
5.47 Configuração deformada do cilindro, visto pelo ângulo 1. . . . . . . . .
80
5.48 Configuração deformada do cilindro, visto pelo ângulo 2. . . . . . . . .
81
5.49 Configuração deformada do cilindro, visto pelo ângulo 3. . . . . . . . .
81
Nomenclatura
Letras Latinas
x Vetor posição de um ponto no espaço em relação à sua origem
u Vetor deslocamento do ponto no espaço
ℓr
Comprimento da barra de treliça na configuração de referência é
ℓ Comprimento da barra de treliça na configuração deformada
N
Força normal
Ar
Área da seção transversal na configuração de referência
E
Módulo de Elasticidade
U
Energia de deformação da barras de treliça
Vr
Volume da barra na configuração de referência
D
Módulo tangente de rigidez elástica
L
Matriz composta por matrizes identidade
p Vetor de deslocamentos nodais
lr
Vetor de comprimento da barra na configuração de referência
I Matriz identidade
l
Vetor de comprimento da barra na configuração deformada
P
Vetor dos esforços nodais internos
kt
Matriz de rigidez tangente de uma barra
ke
Matriz de rigidez elástica da barra de treliça
kg
Matriz de rigidez geométrica da barra de treliça
Rint
Vetor dos esforços internos da estrutura
Kt
A
Matriz de rigidez tangente da estrutura
Matriz de conectividade do elemento de barra e
r Vetor dos deslocamentos nodais da estrutura no sistema global
g Vetor de forças residuais
p Vetor de deslocamento nodal
qef
Vetor de força externa aplicada
qi
Vetor de forças internas
s
Comprimento de corda
D
Matriz diagonal
I0
Número de iterações requeridas
Id
Número de iterações desejadas
Li
Comprimento do arco acumulado
Ld
Comprimento de corda acumulado desejado
Letras Gregas
ρ
Estiramento
ε
Deformação linear
σ
Tensão nominal
σ0
Tensão inicial
ϕ
Energia de deformação específica por unidade de volume
ξ
Vetor de coordenadas nodais
δ
Representa as grandezas virtuais
δε
Deformação virtual
λ
Fator de carregamento
ψ
Parâmetro de dimensionamento
∆l
Comprimento de corda incremental
∆p Vetor de deslocamento incremental
∆λ
Fator de carregamento incremental
δ p Vetor de deslocamento iterativo
δλ
Fator de carregamento iterativo
δ p̄ Vetor de deslocamento que se origina do método de Newton-Raphson
δ pt
Vetor de deslocamento tangencial
∆qef
Vetor forças externas incremental
λd
∆ld
γ
Fator de carregamento predefinido
Comprimento de corda incremental desejado
Parâmetro de comprimento acumulado
1
1
Introdução
Com o crescente avanço tecnológico e a consolidação do Método dos Elementos
Finitos, houve uma grande necessidade de analisar estruturas com características complexas, levando em consideração as não linearidades geométricas e materiais. Desse
modo, o desenvolvimento e a aplicação de procedimentos numéricos para seguir trajetórias de solução em problemas estruturais não lineares têm recebido uma atenção
significativa nas últimas décadas.
As análises do comportamento de sistemas estruturais requerem a solução das equações de equilíbrio mediante a variação de parâmetros do sistema, tais como: intensidade
de carga, amplitudes de deslocamento, variáveis de projeto, entre outros. Estas análises são de considerável interesse para a localização e caracterização dos pontos críticos
(pontos limites e de bifurcação) nestes conjuntos de soluções (RIKS, 1984).
Um método de análise estrutural ideal deve ser capaz de traçar inteiramente a trajetória de equilíbrio de uma estrutura, podendo apresentar um comportamento suave
quanto um fortemente não linear, mesmo com a presença de pontos limites de carga e
de deslocamento e, possivelmente, pontos de bifurcação. Isto é possível através da utilização de procedimentos incrementais-iterativos como o de Método de Newton-Raphson.
Para tanto, foram desenvolvidas várias técnicas de solução para obter a trajetória
de equilíbrio. O Método de de Newton-Raphson com controle de Carga foi o primeiro
método desenvolvido com este propósito, mas ele falha na vizinhança do ponto limite.
Para superar as dificuldades em ultrapassar os pontos limites, foram introduzidas as
técnicas de controle de deslocamento (BATOZ; DHATT, 1979). No entanto, para sistemas
1 Introdução
2
estruturais que exibem comportamento de snap-through ou snap-back, estas técnicas
conduzem ao erro.
Como alternativa, ou era necessário alternar entre os controles de carga e de deslocamento (SABIR; LOCK, 1972) utilizando as molas artificiais de Wrigth e Gaylord (1968),
ou lançar mão de técnicas que optam por abandonar as iterações de equilíbrio nas imediações do ponto limite (BERGAN; SOREIDE, 1978; BERGAN et al., 1978). Contudo, estas
técnicas são desenvolvidas especificamente para uma determinada estrutura.
A motivação dos Métodos de Continuação consistia em encontrar respostas para
perguntas sobre o comportamento de uma estrutura como as citadas por Crisfield (1991)
(com tradução nossa):
• Realmente era um ponto limite ou foi o procedimento de solução iterativa que entrou em colapso?
• Isso foi apenas um máximo local e a estrutura ainda pode ser carregado, estando longe de sofrer um colapso?
• Como é o processo de colapso, dúctil ou frágil?
Os Métodos de Continuação são, basicamente, métodos previsores-corretores. Um
número considerável de diferentes Métodos de Continuação são propostos na literatura,
entre os quais, o mais popular é o método chamado de Método da Corda originalmente
desenvolvido por Riks (1972, 1979) e Wempner (1971).
O Método da Corda é composto pelo seguintes procedimentos de solução: formulação inicial, onde uma equação de restrição é formulada, método previsor, método corretor e controle do comprimento de corda. Além disso, os Métodos da Corda diferenciamse uns dos outros através das variações de ao menos um destes procedimentos.
Diferentes versões do Métodos da Corda desenvolvido por Riks (1972, 1979) e
Wempner (1971) foram propostas por Ramm (1981), Crisfield (1981, 1991), Riks (1984),
Forde e Stieme (1987), Al-Rasby (1991), Fafard e Massicotte (1991), Wagner (1991),
Wriggers P. e Miehe (1987) com o propósito de aprimorar o método original modifi-
1 Introdução
3
cando a equação de restrição ou algum outro conceito pertencente ao passo corretor.
Do mesmo modo, Crisfield (1991), Feng et al. (1996), Neto e Feng (1999), Bergan
e Soreide (1978), Bergan et al. (1978) propuseram diferentes métodos para aprimorar a
determinação do sinal correto para o carregamento incremental inicial, pertencente ao
passo previsor.
Além disso, Teng e Lou (1997), Teng e Luo (1998) propuseram a adição um parâmetro de comprimento de corda acumulado ao Método da Corda desenvolvido por Crisfield (1991), com o objetivo de obter um comprimento de corda melhor dimensionado.
Este método ainda proporciona a possibilidade de convergência em carregamentos predefinidos pelo usuário.
Devido ao grande volume de publicações que propõem diferentes técnicas e procedimentos para aprimorar algum ponto do Método da Corda, foi percebida a necessidade
de verificar e testar estas modificações e, reunir as que apresentarem um melhor desempenho, desenvolvendo, deste modo, um Método da Corda mais robusto.
Com isso, a motivação principal que levou a realização deste trabalho consiste em
contribuir, consistentemente, para o aprimoramento do Método da Corda.
O objetivo principal deste trabalho é testar e propor um procedimento mais eficiente
para o Método da Corda por meio da análise de alguns Métodos da Corda distintos,
tornando-o mais robusto.
Primeiramente, será considerado o Método da Corda Cilíndrico, desenvolvido por
Crisfield (1991), que utiliza uma equação de restrição quadrática. Este método fornecerá
duas raízes possíveis, havendo possibilidade de ocorrerem raízes complexas. Com isso,
uma versão linearizada deste método foram propostas por Schweizerhof e Wriggers
(1986), na qual apenas uma raiz apropriada para o passo corretor é encontrada. Esta
versão também será considerada neste trabalho para fins comparativos.
Além disso, métodos previsores mais eficientes também serão analisados, levando
em consideração os trabalhos desenvolvidos por Crisfield (1991, 1997) e Feng et al.
(1995, 1996). A fim de aumentar a eficiência do Método da Corda, evitando que a
1 Introdução
4
trajetória de equilíbrio siga pela direção errada, ocorrendo o problema de doubling back
e , consequentemente a a falência do método.
Com o intuito de obter um maior controle sobre a execução do Método da Corda,
este será modificado através da introdução de um novo parâmetro fundamentado no
conceito do comprimento de corda acumulado. Contudo, este novo método é nomeado
por Método da Corda Acumulado e foi proposto por Teng e Lou (1997), Teng e Luo
(1998).
Através desta modificação, serão acrescentados ao Método da Corda dois parâmetros de controle. O primeiro diz respeito à convergência do método em uma carga
definida no início do processo. Este controle se faz necessário para analisar o comportamento das estruturas para determinados níveis de carga. O segundo parâmetro refere-se
ao controle do comprimento da corda, este controle é essencial para que seja obtido uma
trajetória de equilíbrio de um modo eficiente através da escolha apropriada do comprimento de corda. Assim, o problema de retorno da trajetória para passos já convergidos
(doubling back) é totalmente resolvido.
Os procedimentos estudados serão testados através de simulações numéricas utilizando estruturas treliçadas. As simulações serão desenvolvidas computacionalmente
utilizando o programa Matlab.
Dentro de todo esse contexto, o conteúdo dessa dissertação de mestrado foi dividido
em seis capítulos.
O Capítulo 2 apresenta a formulação da Teoria Geometricamente Exata utilizada
para estruturas em forma de treliças.
O Capítulo 3 refere-se aos Métodos de Continuação. Neste capítulo serão descritos
os Métodos da Corda considerados para este trabalho, apresentando suas respectivas
formulações.
No Capítulo 4, serão apresentadas as formulações referentes aos métodos propostos
para otimização do Método da Corda. Este capítulo diz respeito à implementação do
Método da Corda Acumulado e do Método previsor baseado no critério do produto
1 Introdução
5
interno.
No Capítulo 5, toda a teoria desenvolvida nos capítulos anteriores serão aplicadas
a diferentes estruturas treliçadas para testar o desempenho e a eficiência dos Métodos
descritos.
Finalmente, o Capítulo 6 apresenta as conclusões e ponderações referente aos métodos testados no capítulo anterior.
6
2
Teoria Geometricamente Exata
para Treliças
O objetivo deste capítulo é apresentar a formulação da teoria estrutural geometricamente exata para a análise não linear de treliças espaciais. As formulações apresentadas
foram propostas por Pimenta (1986, 1996) e serão restritas apenas aos materiais elásticos.
Para tanto, supor-se-á que as barras das treliças sejam prismáticas – na qual sua
seção transversal é considerada constante ao longo do eixo – e que o carregamento
externo seja aplicado nos nós da estrutura.
2.1
Considerações Sobre Tensões e Deformações
O sistema de coordenadas globais que será considerado para as estruturas aqui estudadas encontra-se representado pela Figura 2.1, na qual são apresentadas as configurações de referência e deformada (atual) de uma barra de treliça através do espaço.
br
Conf. de
referência
lr
ar
x ar
e3
e1
b
rb
l
ra
Conf.
atual
a
xa
e2
Figura 2.1: Descrição da deformação de uma barra.
2.1 Considerações Sobre Tensões e Deformações
7
Neste sistema, o vetor x é denominado de vetor posição de um ponto no espaço em
relação à sua origem e o vetor u é chamado de vetor deslocamento do ponto no espaço.
O comprimento da barra de treliça na configuração de referência é representado por ℓr
e na configuração deformada por ℓ .
Define-se como medida de deformação qualquer grandeza que compare os comprimentos da barra de treliça nas configurações de referência e atual. Uma medida básica
de deformação é o estiramento, representado pela letra grega ρ , e é dado por
ρ=
ℓ
ℓr
(2.1)
A família de medidas de deformação é definida através da seguinte relação
{
ε=
(ρ m − 1) /m , m ̸= 0
ln ρ ,
m=0
(2.2)
No entanto, neste trabalho, será considerada apenas a medida de deformação linear,
definida na Equação 2.3, como medida de deformação da barra de treliça.
ε = ρ −1
(2.3)
Ao examinar a Figura 2.1 que apresenta as configurações de referência e deformada
de uma barra de treliça, nota-se que atua sobre a barra na configuração deformada uma
força normal denominada por N. A tensão nominal (σ ) atuante na barra na configuração
deformada é dada por
σ=
N
Ar
(2.4)
onde Ar é a área da seção transversal na configuração de referência. Sendo assim, a
partir da Equação 2.4 observa-se que
2.2 Equações Constitutivas
8
N = σ Ar
(2.5)
2.2 Equações Constitutivas
Nesta seção, tratar-se-á, de forma concisa, as relações constitutivas elásticas para as
barras de treliças.
Uma barra está em regime elástico se existir uma relação em que cada deformação
associa a uma só tensão. A equação constitutiva elástica pode ser expressa pela função
σ = σ (ε )
(2.6)
A Equação 2.7 se refere à equação constitutiva linear dada pela lei de Hooke e
trata-se de um caso particular da Equação 2.6.
σ = E ε + σ0
(2.7)
onde E é o módulo de Elasticidade e σ0 é a tensão inicial.
A seguir, pode-se definir que
ϕ (ε ) =
∫ ε
0
σ (η )d η
(2.8)
onde ϕ é a energia de deformação específica por unidade de volume da configuração de
referência.
Com isso, a energia de deformação (U) da barras de treliça é dada por
U = Ar ℓr ϕ
(2.9)
2.3 Equações de Equilíbrio
9
onde o volume da barra na configuração de referência é dado por V r = Ar ℓr .
Para a lei de Hooke definida pela Equação 2.7, tem-se
1
ϕ = E ε 2 + σ0 ε
2
(2.10)
Considerando a definição da Equação 2.8, conclui-se que
σ=
dϕ
dε
(2.11)
onde, observa-se que, ϕ pode ser utilizado como potencial para as tensões.
O módulo tangente de rigidez elástica, representado pela letra D, é definido por
dσ
dε
(2.12)
d2σ
D= 2
dε
(2.13)
D = D (ε )
(2.14)
D=
ou
Em geral, pode-se dizer que
Contudo, para a lei de Hooke, D é considerada uma constante dada por
D=E
2.3
(2.15)
Equações de Equilíbrio
Considerando novamente a Figura 2.1, nota-se que a e b correspondem às extremidades de uma barra de treliça na configuração deformada. Com isso, é possível deter-
2.3 Equações de Equilíbrio
10
minar os vetores x e u destas extremidades e assim obter os vetores de coordenadas (ξ )
e deslocamentos (p) nodais da barra representados, respectivamente, por
[
ξ=
]
(2.16)
xb
[
e
xa
p=
ua
6x1
]
ub
(2.17)
6x1
O vetor de comprimento da barra na configuração de referência (lr ) é dado pelo
vetor que une a extremidade a com a extremidade b. Considera-se a sua origem na
extremidade a, da seguinte forma
lr = L ξ
(2.18)
considerando que L é uma matriz composta por matrizes identidade (I), como
[
−I3 I3
L=
]
(2.19)
Da mesma forma, o vetor de comprimento da barra na configuração deformada (l)
é definido por
l = L (ξ + p)
(2.20)
Os comprimentos da barra nas configurações de referência e deformada são calculados, respectivamente, do seguinte modo
√
ℓr =
lrT lr
(2.21)
lT l
(2.22)
√
e
ℓ=
2.3 Equações de Equilíbrio
11
Por conseguinte, pode-se definir a deformação linear ε em função dos deslocamentos nodais como se observa na Equação 2.23
ε = ε (p)
(2.23)
ou, explicitamente,
√
ε=
(ξ + p)T LT L (ξ + p)
ℓr
−1
(2.24)
Através do Princípio dos Trabalhos Virtuais, torna-se possível obter o vetor dos
esforços nodais internos, representado pela letra P, para a barra de treliça da seguinte
forma
PT δ p = Ar ℓr σ δ p
(2.25)
onde δ representa as grandezas virtuais. No entanto, é preciso lembrar que a deformação
virtual (δ ε ) é definida como
(
δε =
∂ε
∂p
)T
δp
(2.26)
Para tanto, define-se o novo vetor
b=
∂ε
∂p
(2.27)
Ao passo que a Equação 2.26 torna-se
δ ε = bT δ p
Introduzindo a expressão acima na Equação 2.25, obtém-se
(2.28)
2.3 Equações de Equilíbrio
12
P = Ar ℓr σ b
(2.29)
Ao aplicar as técnicas de diferenciação para a Equação 2.24, considerando a definição dada pela Equação 2.27, obtém-se a seguinte expressão para o vetor b
b=
1 T
L l
ℓr ℓ
(2.30)
Substituindo o vetor b (Equação 2.30) na Equação 2.29, obtém-se que o vetor de
forças nodais internas da barra é dado por
1
P = N LT l
ℓ
(2.31)
Do mesmo modo, a matriz de rigidez tangente (kt ) de uma barra é definida por
kt =
∂P
∂p
(2.32)
a qual, por diferenciação da Equação 2.29, é possível obter sua definição da seguinte
forma
(
∂σ
kt = A ℓ b
∂p
r r
)T
(
+A ℓ σ
r r
∂b
∂p
)
(2.33)
Para resolver a expressão acima, deve-se considerar, utilizando a regra da cadeia,
que
∂σ
dσ ∂ ε
=
∂p
dε ∂ p
(2.34)
Nota-se que, ao considerar as Equações 2.12 e 2.15 fornecidas na seção anterior e
as Equações 2.28 e 2.30, obtém-se
2.3 Equações de Equilíbrio
13
∂σ
D
= Db = r LT l
∂p
ℓℓ
(2.35)
Além disso, também por diferenciação da Equação 2.30, tem-se
∂b
1
1
= r LT L − r 3 LT l lT L
∂p ℓ ℓ
ℓℓ
(2.36)
Finalmente, a matriz de rigidez tangente da barra é obtida introduzindo as Equações
2.35 e 2.36 na Equação 2.33
DAr
kt = r LT
ℓ
(
)
)
(
1 T
N T
1 T
l l L + L I3 − 2 l l L
ℓ2
ℓ
ℓ
(2.37)
A matriz de rigidez apresentada acima pode ser dividida do seguinte modo
kt = ke + kg
(2.38)
As parcelas da Equação 2.38 são apresentadas a seguir pelas Equações 2.39 e 2.40.
DAr
ke = r LT
ℓ
(
)
1 T
ll L
ℓ2
(
)
N T
1 T
kg = L I3 − 2 l l L
ℓ
ℓ
(2.39)
(2.40)
onde ke é a matriz de rigidez elástica e kg a matriz de rigidez geométrica da barra de
treliça.
A partir de então, é possível calcular o vetor dos esforços internos (Rint ) e a matriz
de rigidez tangente (Kt ) de uma dada estrutura através de P e kt para todas as barras
desta estrutura, como se observa nas Equações 2.41 e 2.42.
2.3 Equações de Equilíbrio
14
Rint = ∑ Ae T P
(2.41)
Kt = ∑ Ae T kt Ae
(2.42)
e
e
onde A é a matriz de conectividade do elemento de barra e definida por
p=Ar
(2.43)
O vetor r, definido pela Equação 2.43, consiste no vetor dos deslocamentos nodais
da estrutura no sistema global.
15
3
Métodos de Continuação
A análise não linear de estabilidade estrutural considera os comportamentos não lineares tanto geométricos quanto materiais de uma dada estrutura. Quando há uma maior
complexidade estrutural, os sistemas de equações não lineares que regem seu comportamento necessitam de técnicas que descrevam por completo suas trajetórias de equilíbrio. Consequentemente, os métodos de seguimento de trajetórias foram desenvolvidos,
como o próprio nome diz, para traçar as trajetórias de soluções não lineares arbitrárias,
identificando os pontos críticos e comportamentos de equilíbrio dessas estruturas.
O produto final de um dado sistema de equações não lineares de uma estrutura
é transferido para um sistema de coordenadas retangulares, como o apresentado pela
Figura 3.1, cujos eixos representam o fator de carga e o deslocamento de um dado
grau de liberdade e recebe o nome de diagrama de carga-deslocamento. A trajetória
adquirida neste diagrama representa, ponto a ponto, o estado de equilíbrio da estrutura
e, dessa forma, é conhecida como trajetória de equilíbrio da estrutura.
Convém observar que, quando o estado de equilíbrio inicial de uma estrutura é considerado com deformações e tensões iniciais nulas (como adotado no presente trabalho),
sua trajetória de equilíbrio terá como ponto de partida a origem do diagrama de cargadeslocamento, apresentado na Figura 3.1.
Os diagramas de carga-deslocamento são compostos por trajetórias de equilíbrio e
pontos críticos que representam o comportamento de uma estrutura. Geralmente, são
encontradas mais de uma trajetória de equilíbrio para uma determinada estrutura, como
mostra a Figura 3.2. A trajetória que cruza o estado inicial é designada por trajetó-
3 Métodos de Continuação
16
Figura 3.1: Diagrama carga-deformação.
ria fundamental ou primária e as demais são chamadas de trajetórias secundárias e
terciárias.
Figura 3.2: Trajetórias primária e secundária.
Por sua vez, os pontos críticos são divididos em pontos limites e pontos de bifurcação. Estes últimos caracterizam-se por localizar o cruzamento (intersecção) entre duas
ou mais trajetórias de equilíbrio. Os pontos limites encontram-se localizados tangen-
3 Métodos de Continuação
17
cialmente à trajetória de equilíbrio, isto é, são paralelos ou ao eixo dos deslocamentos
(pontos limite de deslocamento) ou ao eixo do fator de carga (pontos limite de carga).
Caso, durante o processo do Método de Continuação, estes pontos forem ultrapassados ou por um incremento de carga ou por um deslocamento adicional, o equilíbrio
estático da estrutura é perdido e esta passará dinamicamente para a próxima posição de
equilíbrio pós-crítico. Estes comportamentos dinâmicos são, respectivamente, chamados de snap-through, como se vê na Figura 3.3, e de snap-back, apresentado na Figura
3.4.
Figura 3.3: Comportamento de carga limite (snap-through).
Figura 3.4: Comportamento de deslocamento limite (snap-back).
As várias técnicas de seguimento de trajetórias diferenciam-se umas das outras pelas
condições de restrição adotadas para a resolução do sistema de equações não lineares.
Portanto, nos próximos parágrafos, os métodos que representam maior relevância para
este trabalho serão sucintamente apresentados.
O Método de Controle de Carga caracteriza-se por apresentar uma equação de restrição fundamentada na equação de equilíbrio da estrutura, definida como a diferença
entre as forças internas e o carregamento externo prescrito. Esta condição de restrição resulta exatamente no método de Newton-Raphson padrão e, consequentemente,
apresenta comportamento instável ao ultrapassar os pontos de carga limites, como re-
3 Métodos de Continuação
18
presenta a Figura 3.5. Dessa forma, torna-se incapaz de traçar por completo a trajetória
de equilíbrio de uma estrutura com comportamento não linear.
Figura 3.5: Método de Controle de
Carga.
Figura 3.6: Método de Controle de
Deslocamento.
O Método de Controle de Deslocamento difere-se do Método de Controle de Carga
por ser capaz de ultrapassar os pontos limites de carga. Sua equação de restrição é uma
função dependente dos incrementos dos deslocamentos e pode ser formulada para um
único componente do campo de deslocamentos que seja decisivo para o processo. Este
método mostra-se sensível na presença dos pontos limites de deslocamento, tornando-se
instável ao ultrapassá-los, como ilustra a Figura 3.6. Consequentemente, a análise não
é completada.
O Método da Corda é, portanto, o único tipo de Método de Seguimento de Trajetórias capaz de obter uma trajetória de equilíbrio completa, ultrapassando os pontos
limites e de bifurcação e ainda descrevendo os comportamentos de snap-through e snapback de uma estrutura.
Como já mencionado, Riks (1972, 1979) e Wempner (1971) foram os pioneiros no
desenvolvimento do Método da Corda, cuja característica principal consiste em encontrar a intersecção entre a equação de equilíbrio da estrutura e um dado comprimento de
corda.
3 Métodos de Continuação
19
Para tanto, estes autores propuseram uma equação de restrição linear em relação
ao fator de carga e ao campo de deslocamentos que fosse capaz de descrever um plano
normal tangente à curva calculada para o último passo convergido (Figura 3.7).
Figura 3.7: Método da Corda proposto por Riks (1979).
Poucos anos depois, Crisfield (1981) propôs que o Método da Corda fosse utilizado
com uma equação de restrição não linear capaz de descrever uma superfície esférica em
torno do último passo convergido, denominado de Método da Corda Esférico.
Em um trabalho posterior, Crisfield (1991) propôs que os termos de carga da equação de restrição não linear fossem desprezados com o objetivo de tornar o método mais
estável. Com isso, a superfície descrita por esta equação passou a ser cilíndrica ao invés
de esférica. Este método foi chamado de Método da Corda Cilíndrico e é ilustrado pela
Figura 3.8.
Posteriormente, uma linearização de forma consistente para o Método de NewtonRaphson da equação de restrição quadrática (CRISFIELD, 1991) foi desenvolvida por
Schweizerhof e Wriggers (1986), como mostra a Figura 3.9.
O propósito desta versão linearizada do Método da Corda Cilíndrico é obter apenas
uma raiz como solução para o método, eliminando o surgimento de raízes complexas.
Consequentemente, este método necessita de um menor esforço computacional para ser
3 Métodos de Continuação
20
Figura 3.8: Método da Corda proposto por Crisfield (1991).
executado, convergindo de forma mais rápida e eficiente. No presente trabalho, este
método será denominado por Método da Corda de Schweizerhof-Wriggers.
Figura 3.9: Método da Corda de Schweizerhof-Wriggers (SCHWEIZERHOF; WRIGGERS,
1986).
O Método de Continuação é composto por 4 procedimentos de solução: (1) Formulação inicial, onde uma equação de restrição é formulada 1 ; (2) Corretor; (3) Previsor;
1 Esta
etapa é comumente chamada de Estratégia de Parametrização.
3.1 Formulação Inicial para os Métodos de Continuação
21
e (4) Controle do comprimento de corda. Estes procedimentos serão apresentados a
seguir.
3.1
Formulação Inicial para os Métodos de Continuação
A formulação que será apresentada abaixo, desenvolvida por Riks (1972, 1979) e
Wempner (1971), trata-se da formulação geral para o desenvolvimento de todos os tipos
de Método da Corda existentes.
O objetivo do Método de Continuação é encontrar a intersecção entre o sistema de
equações não lineares de uma dada estrutura e um dado comprimento de corda através
de uma equação de restrição. Para que isso se torne possível, é preciso criar um conjunto
de equações envolvendo a equação de equilíbrio da estrutura e a equação de restrição
escolhida.
Inicialmente, considera-se que o sistema de equações não lineares da estrutura seja
apresentado na forma da seguinte equação de equilíbrio
g (p, λ ) = qi (p) − λ qef = 0
(3.1)
no qual g é o vetor de forças residuais, p é o vetor de deslocamento nodal, λ é um
escalar que representa o fator de carregamento, qef é o vetor de força externa aplicada e
qi é o vetor de forças internas.
Para completar este conjunto de equações, adiciona-se uma equação de restrição
através do conceito do comprimento de corda, representado pela letra s e definido pela
Equação 3.2.
∫
s=
ds
(3.2)
3.1 Formulação Inicial para os Métodos de Continuação
22
O termo diferencial ds é dado pela Equação 3.3.
ds =
√
dpT dp + d λ 2 ψ 2 qef 2 qef
(3.3)
onde ψ é um parâmetro de dimensionamento e sua função é determinar a influência dos
termos de carga da equação no comportamento do Método da Corda.
A interseção entre o comprimento de corda s, representado pela Equação de restrição 3.2, e a Equação de equilíbrio 3.1 só será possível se a equação de equilíbrio se
tornar uma função de s, como mostrado na Equação 3.4.
g (s) = qi (p (s)) − λ (s) qef = 0
(3.4)
Contudo, ao adotar esta aproximação, torna-se muito difícil limitar com êxito a
tendência ao equilíbrio e, neste caso, utiliza-se Métodos Previsores-Corretores
2
para
solucionar a Equação 3.4.
Dessa forma, modifica-se a forma diferencial da Equação 3.3 por sua forma incremental correspondente, que resulta em
(
)
a = ∆pT ∆p + ∆λ 2 ψ 2 qef T qef − ∆l 2 = 0
(3.5)
onde a é a forma geral da equação de restrição que é escrita em termos do vetor de deslocamento e do fator de carga, ambos desconhecidos. Ainda na Equação 3.5, define-se
que o ∆l é o comprimento de corda incremental e representa o raio fixo para a interseção desejada com a trajetória de equilíbrio. O vetor de deslocamento ∆p e o fator de
carregamento ∆λ são incrementais e relativos ao último estado convergido da equação.
No Método da Corda, o fator de carregamento (λ ) é considerado uma variável desconhecida. Dessa forma, o número de variáveis a serem determinadas será n + 1, onde
n é o número de variáveis desconhecidas do vetor de deslocamento (p) e 1 é a variável
2 Original
do inglês: Predictor-Corrector Method.
3.1 Formulação Inicial para os Métodos de Continuação
23
que corresponde ao fator de carregamento (λ ).
As variáveis mencionadas acima são determinadas através da aplicação do método
iterativo de Newton-Raphson para as Equações 3.1 e 3.5. Com isso, são obtidas as
mudanças iterativas referentes ao vetor de deslocamentos e ao fator de carga.
Quando isso se faz necessário, o Método de Newton-Raphson é melhor introduzido
através da linearização das Equações 3.1 e 3.5. Essa linearização é obtida através da
Série de Taylor truncada, conforme as Equações 3.6 e 3.7.
gi+1 = gi +
ai+1 = ai +
∂g i ∂g i
δp +
δ λ = gi + Kt δ pi − qef δ λ i = 0
∂p
∂λ
( )T
∂a i ∂a i
δp +
δ λ = ai + 2 ∆pi δ pi + 2∆λ i δ λ i ψ 2 qef T qef = 0
∂p
∂λ
(3.6)
(3.7)
Nota-se que os índices sobrescritos i + 1 representam o novo passo iterativo e i o
passo iterativo imediatamente anterior a este3 . Ainda nas Equações 3.6 e 3.7 , Kt é
a matriz de rigidez tangente da estrutura e δ pi e δ λ i são, respectivamente, o vetor de
deslocamento iterativo e o fator de carregamento iterativo.
Primeiramente, combinam-se as Equações 3.6 e 3.7 e, num segundo momento,
considera-se gi+1 = ai+1 = 0 para calcular δ pi e δ λ i na forma de um sistema de matrizes, conforme a Equação 3.8.
(
δ pi
δλi
)
[
=−
Kt
−qef
( )T
2 ∆pi
2∆λ i ψ 2 qef T qef
]−1 (
gi
ai
)
(3.8)
É importante observar que a matriz estendida (Equação 3.8) não é simétrica para
valores incrementais desconhecidos do vetor de deslocamento e do fator de carga. Esta
matriz não se torna singular nos pontos limites, embora a matriz de rigidez tangente
3A
partir daqui, estes índices sobrescritos serão frequentes nas equações. Portanto, considera-se sempre que i refere-se ao passo iterativo atual e (i + 1) ao passo iterativo sucessor.
3.2 Métodos Corretores
24
(Kt ), que a compõe, seja singular nestes pontos. Porém, nos pontos de bifurcação esta
matriz se torna singular.
O sistema de equações não-simétricas supracitado é normalmente resolvido através
de técnicas de particionamento capazes de utilizar a estrutura simétrica da matriz de
rigidez Kt e, com isso, se torna um algoritmo eficiente. No entanto, a propriedade de
não singularidade nos pontos limites é perdida (WRIGGERS, 2008).
A partir deste ponto, existem várias maneiras diferentes de se obter a solução do
sistema apresentado na Equação 3.8. Porém, neste trabalho, serão abordados apenas
dois métodos formulados nas seções a seguir.
3.2
Métodos Corretores
Nesta fase de solução do Método da Corda, são aplicados os procedimentos para
resolução do sistema de equações estendido. São desenvolvidas as diferentes equações
de restrição que define o definirá o tipo de Método da Corda.
Neste seção são apresentados dois tipos de Métodos Corretores, um com equação
de restrição quadrática e o outro com equação de restrição linearizada.
3.2.1
Método da Corda Esférico
Nesta seção, será apresentada a formulação desenvolvida por Crisfield (1981, 1991)
através das equações gerais apresentadas anteriormente.
Ao invés de solucionar o sistema definido pela Equação 3.8, Crisfield (1981) propôs
a resolução direta da equação de restrição (Equação 3.6), levando em consideração o
procedimento desenvolvido por Batoz e Dhatt (1979) para o Método de Controle de
Deslocamento 4 . O diferencial deste método advém da divisão do vetor de deslocamento
iterativo, δ p, em duas partes.
4 Original
do inglês: Displacement Control Method.
3.2 Métodos Corretores
25
Para tanto, o Método de Newton-Raphson é alterado ao se consider um novo fator
de carregamento incremental desconhecido (λ i+1 = λ i + δ λ i ). Com isso, o vetor de
deslocamento iterativo (δ p) se torna
(
)
δ pi = −(Kt )−1 g pi , λ = −(Kt )−1 gi + δ λ i (Kt )−1 qef
(3.9)
Considerando a última forma apresentada na Equação 3.9, a subdivisão de δ p pode
ser expressa como
δ pi = δ p̄i + δ λ i δ pt i
(3.10)
As duas parcelas que compõem a Equação 3.10 são dadas por
δ p̄i = −(Kt )−1 gi
(3.11)
δ pt i = (Kt )−1 qef
(3.12)
onde δ p̄i é a mudança iterativa que se origina do método de Newton-Raphson para o
carregamento controlado (com fator de carga fixo λ i ) e δ pt i é o vetor de deslocamento
iterativo correspondente ao vetor de força externa qef .
Ao aplicar as mudanças acima mencionadas, o novo vetor de deslocamento incremental é dado pela Equação 3.13.
∆pi+1 = ∆pi + δ pi+1 = ∆pi + δ p̄i + δ λ i δ pt i
onde δ λ i é a única variável desconhecida.
Por consequência, a equação de restrição (Equação 3.5) é reescrita como
(3.13)
3.2 Métodos Corretores
((
∆p
)
i T
(
∆p + ∆λ
i
)
i 2
26
)
ψ qef qef =
2
T
((
)
i+1 T
∆p
∆p
i+1
( i+1 )2 2 T )
+ ∆λ
ψ qef qef = ∆l 2
(3.14)
O valor de ai foi considerado como igual a 0 pelo simples fato de que o raio da
interseção desejada é sempre mantido constante durante o processo iterativo.
A substituição de ∆pi+1 (Equação 3.13) na Equação 3.14 resulta em uma equação
quadrática, como mostra a Equação 3.15.
a1 δ λ 2 + a2 δ λ + a3 = 0
(3.15)
Percebe-se que os coeficientes desta equação quadrática são iguais a
(
)T
a1 = δ pt i δ pt i + ψ 2 qef T qef
)
(
a2 = 2δ pt i ∆pi + δ p̄i + 2∆λ i ψ 2 qef T qef
)T ( i
)
(
(
)2
a3 = ∆pi + δ p̄i
∆p + δ p̄i − ∆l 2 + ∆λ i ψ 2 qef T qef
(3.16)
(3.17)
(3.18)
A Equação 3.15 é nomeada de equação de restrição esférica, pois uma esfera é
introduzida próxima ao último passo convergido e, consequentemente, este Método da
Corda ficou conhecido como Método da Corda Esférico5 .
No presente trabalho, será considerado o Método da Corda Cilíndrico, que será
discorrido a seguir.
3.2.2
Método da Corda Cilíndrico
O Método da Corda Cilíndrico é uma variação do Método da Corda Esférico, onde o
parâmetro de dimensionamento de carga (ψ ) presente na equação de restrição esférica é
desconsiderado. Sua equação de restrição é chamada de equação de restrição cilíndrica.
5 Original
do inglês: Spherical Arc-length Method.
3.2 Métodos Corretores
27
O parâmetro de dimensionamento ψ determina a influência do controle de carga e
de deslocamentos no comportamento do método da corda. Quando este é desconsiderado, os termos referentes ao controle de carga da equação de restrição (Equação 3.15)
são desprezados e, com isso, o Método da Corda torna-se mais eficiente em ultrapassar
os pontos limites e os de bifurcação (AL-RASBY, 1991).
Isto se deve ao fato de que se o parâmetro de dimensionamento tender a um (ψ → 1),
o termo da equação de restrição (Equação 3.5), que contém forças externas e carregamentos, tenderá ao infinito. Por consequência, o Método da Corda se comportará de
forma similar ao Método de Controle de Carga e perderá a sua característica de ultrapassar os pontos críticos.
Não obstante, se o parâmetro de dimensionamento tende a se tornar igual a zero
(ψ → 0), o termo que contém cargas e forças externas da Equação 3.15 também tenderá a zero. Assim, sua capacidade de ultrapassar os pontos limites e os de bifurcação
continuará vigente.
Portanto, ao se considerar ψ = 0, a Equação 3.15 se torna igual a
( i )T i
∆p ∆p = ∆l 2
(3.19)
Como consequência da Equação 3.15, a equação quadrática Equação 3.15 é apresentada como mostra a Equação 3.20
a1 δ λ 2 + a2 δ λ + a3 = 0
(3.20)
os coeficientes da Equação 3.20 são dados pelas Equações 3.21, 3.22 e 3.23.
(
)T
a1 = δ pt i δ pt i
)
(
)T ( i
a2 = 2 δ pt i
∆p + δ p̄i
(
)T ( i
)
a3 = ∆pi + δ p̄i
∆p + δ p̄i − ∆l 2
(3.21)
(3.22)
(3.23)
3.2 Métodos Corretores
28
No Método da Corda desenvolvido por Crisfield (1991), por haver uma equação
de restrição quadrática, três situações distintas podem ocorrer na sua resolução. Essa
característica é muito importante neste método, pois o êxito do seu desenvolvimento
depende da escolha apropriada do fator incremental de carga δ λ .
Dentre as possíveis situações para a escolha da raiz apropriada na Equação 3.20,
estão:
1. Duas raízes reais δ λ1 e δ λ2 :
Apenas uma delas será a raiz que levará a trajetória pela direção correta. A outra raiz, provavelmente, terá o sentido contrário. A escolha da raiz que levará à
direção correta é feita através da comparação entre as duas soluções encontradas para o vetor de deslocamento incremental atualizado ∆pi+1 . Isto quer dizer
que serão encontrados ∆p1 i+1 e ∆p2 i+1 , ambos referentes às raízes δ λ 1 e δ λ 2 ,
respectivamente.
A raiz apropriada será aquela que possuir a direção mais próxima do vetor de
deslocamento incremental anterior ∆pi , considerando o menor ângulo entre este
e ∆pi+1 . Isso é obtido ao analisar o maior cosseno do ângulo, através da seguinte
equação
(
cos θ =
∆pi
)T
∆pi+1
∆l 2
)
(
)T
( i )T ( i
δ λ i δ pt i δ pt i
∆p
∆p + δ p̄i
=
+
∆l 2
∆l 2
(3.24)
De um modo mais simplificado, considera-se que
cos θ =
a4 + a5 δ λ
∆l 2
(3.25)
Os coeficiente da Equação 3.25 são dados por
) ( )T
(
a4 = ∆pi + δ p̄i + ∆pi ∆pi
( )T
a5 = ∆pi δ pt i
(3.26)
(3.27)
3.2 Métodos Corretores
29
2. Apenas uma raiz real δ λlin :
Haverá apenas uma raiz real se o valor de a1 da Equação 3.21 for menor ou igual
a zero (a1 ≤ 0). Neste caso, a Equação 3.20 se tornará linear. Com isso, δ λ é
dado por
δ λlin = −
a3
a2
(3.28)
Os coeficientes a2 e a3 são expressos nas Equações 3.22 e 3.23.
3. Raízes não reais (raízes complexas):
A presença de raízes complexas podem ocorrer na vizinhança de pontos de bifurcação quando o comprimento de corda incremental (∆l) atual é razoavelmente
grande. Quando isso ocorre, o processo iterativo é cessado imediatamente e o
comprimento de corda incremental (∆l) é recalculado para o início do novo passo
incremental.
Encontrar raízes complexas é uma das maiores dificuldades deste método e não
há, na literatura, registro de um método que resolva, por definitivo, este problema.
Ao final de todo esse processo δ pi e δ λ i são calculados, atualizando o vetor de deslocamento incremental ∆pi+1 e fator incremental de carga ∆λ i+1 , conforme as Equações
3.52 e 3.30.
∆pi+1 = ∆pi + δ pi
(3.29)
∆λ i+1 = ∆λ i + δ λ i
(3.30)
O fato de se trabalhar com três possibilidades de escolha das raízes torna o Método
da Corda Cilíndrico desvantajoso, pois um esforço computacional extra é necessário
para o cálculo e escolha das raízes apropriadas.
3.2 Métodos Corretores
3.2.3
30
Método da Corda de Schweizerhof-Wriggers
Nesta seção, será apresentada uma equação de restrição para o Método da Corda
proposto por Schweizerhof e Wriggers (1986), derivado de um processo de linearização
consistente.
Para o desenvolvimento deste método, Schweizerhof e Wriggers (1986) utilizam
com base os Métodos da Corda propostos por Crisfield (1981) e Ramm (1981).
Quanto ao método desenvolvido por Ramm (1981), os autores retiram a ideia de
utilizar uma equação de restrição linear. Devido ao fato de que a equação de restrição
linear, por possuir apenas uma raiz, torna-se exatamente satisfeita em toda iteração,
além do beneficio de não ser preciso escolher a raiz correta.
Do mesmo modo, consideram o método desenvolvido por Crisfield (1981) por possui a característica de ser mais estável perante fortes não linearidades. Devido ao fato
deste método utilizar uma equação de restrição quadrática, o que a torna apenas iterativamente satisfeita.
Deste modo, Schweizerhof e Wriggers (1986) sugerem a utilização de uma equação de restrição linearizada de uma maneira consistente utilizando um método do tipo
Newton 6 para que os problemas mencionados acima sejam o resolvidos.
Para isso, ao invés de resolver diretamente a Equação 3.8, opta-se por resolver as
Equações 3.6 e 3.7, considerando gi+1 = ai+1 = 0, como mostram as Equações 3.31 e
3.32.
−gi = Kt δ pi − qef δ λ i
( )T
−ai = 2 ∆pi δ pi + 2∆λ i δ λ i ψ 2 qef T qef
(3.31)
(3.32)
O vetor de deslocamento iterativo δ pi é isolado na Equação 3.31, do seguinte modo
δ pi = −(Kt )−1 gi + δ λ i (Kt )−1 qef
6A
Newton-type scheme.
(3.33)
3.2 Métodos Corretores
31
Nota-se que a Equação 3.33 é idêntica à Equação 3.9 proposta por Batoz e Dhatt
(1979), no qual o vetor de deslocamento incremental é dividido em dois termos, conforme a Equação 3.34.
δ pi = δ p̄i + δ λ i δ pt i
(3.34)
Os dois termos do vetor de deslocamento são dados pelas Equação 3.35 e 3.36.
δ p̄i = −(Kt )−1 gi
(3.35)
δ pt i = (Kt )−1 qef
(3.36)
Ao introduzir diretamente a Equação 3.34 na Equação 3.32, obtém-se
( i)
( i )T ( i
)
a
i
i
i
i 2
T
∆p
δ p̄ + δ λ δ pt + ∆λ δ λ ψ qef qef = −
2
(3.37)
onde, a é dado pela Equação 3.5 e representa a forma geral da equação de restrição
escrita em termos do vetor de deslocamentos e do fator de carga.
A Equação 3.37 é satisfeita para a quando a convergência é atingida. Esta equação
de restrição não precisa ser satisfeita em todo o conjunto de iterações.
O fator de carga é obtido isolando δ λ na Equação 3.37, como mostra a Equação
3.38.
δλ =
( / ) ( )T
− ai 2 − ∆pi δ p̄i
i
(∆pi )T δ pt i + ∆λ i ψ 2 qef T qef
(3.38)
A equação acima é chamada de equação de restrição do Método de Schweizerhof e
de Wriggers.
Os termos de carga da equação de restrição (Equação 3.38) podem ser negligenciados, como já proposto anteriormente na Seção 3.2.2, ao se considerar ψ = 0. Deste
modo, a Equação 3.38 vai influenciar apenas os termos de deslocamentos, de acordo
com a Equação 3.39.
δλ =
i
( / ) ( )T
− ai 2 − ∆pi δ p̄i
(∆pi )T δ pt i
(3.39)
3.3 Métodos Previsores
32
Como já mencionado, em comparação com o Método da Corda criado por Crisfield
(1981), este método possui duas vantagens: a equação de restrição não é apenas iterativamente satisfeita; e sua forma linear apresenta apenas uma única raiz à ser determinada
e, devido a esse fato apresenta uma maior eficiência computacional ao elimina o problema de cálculo e escolha das raízes apropriadas e de ocorrência de raízes complexas.
Deste modo, Schweizerhof e Wriggers (1986) sugerem que a versão consistentemente linearizada seja utilizada sempre que as raízes do Método de Crisfield se tornarem
complexas.
Após calcular o vetor de deslocamento iterativo, δ p, e o fator de carga iterativo, δ λ ,
o novo vetor de deslocamento e o fator de carga incremental são atualizados mediante
as Equação 4.1 e 4.2.
3.3
∆pi+1 = ∆pi + δ pi
(3.40)
∆λ i+1 = ∆λ i + δ λ i
(3.41)
Métodos Previsores
Os Métodos de Continuação são compostos de duas fases: uma previsora e a outra
corretora, descrita na Seção XXX. Esta seção será dedicada à apresentação do Passo
Previsor 7 , cuj objetivo é indicar a direção que será tomada pelo passo iterativo inicial
do Método da Corda.
A fase previsora é de grande importância para todo o processo iterativo, pois é
através desta fase que se constrói um valor de partida adequado para a fase de correção,
utilizando a informação pertencente ao último ponto calculado.
Caso a escolha do passo previsor não seja apropriada, o corretor poderá convergir
para um ponto indesejado, por conta do problema de não linearidade em questão. Desse
modo, o conjunto de pontos de solução encontrado não seguirá uma sequência correta
e, consequentemente, poderá ser encontrada uma trajetória de solução inadequada.
7 Tradução
livre do termo em inglês Predictor Step.
3.3 Métodos Previsores
33
Quando este fenômeno indesejável ocorre, os métodos de continuação podem convergir em pontos previamente calculados, fazendo com que a trajetória de equilíbrio
retorne a estes pontos. Este comportamento é comumente conhecido por double back
ou track back, e caracteriza uma falha na execução do método, pois o comportamento
da estrutura não será conhecido.
A expressão para o previsor origina-se da equação de estimativa desenvolvida por
Crisfield (1981). Nesta equação é adotado o Método Explícito de Euler 8 para previsores
tangenciais 9 , que são obtidos apenas com base na última solução calculada. Para tanto,
considera-se a Equação 3.42 para o vetor de deslocamento incremental.
∆p1 = Kt −1 ∆qef
(3.42)
onde Kt é a matriz de rigidez tangente no início do incremento e ∆qef o vetor incremental de forças externas.
É importante salientar que o número sobreposto ao vetor de deslocamento incremental, ∆p, representa o número de iterações do passo previsor, que sempre tem o
incremento i igual a 1. Isto significa que o passo previsor será calculado apenas para a
primeira iteração.
O vetor incremental de forças externas dado pela Equação 3.42 pode ser definido
como
∆qef = ∆λ 1 qef
(3.43)
no qual ∆λ 1 é o carregamento incremental na primeira iteração e qef o vetor de forças
externas da estrutura.
Substituindo a Equação 3.43 na Equação 3.42, obtém-se
∆p1 = ∆λ 1 Kt −1 qef
(3.44)
Considerando a definição do vetor δ pt dada pela Equação ?? e substituindo-a na
8 Na
nomenclatura inglesa: forward-Euler Scheme.
do inglês: Tangencial predictor.
9 Original
3.3 Métodos Previsores
34
Equação 3.44, tem-se que
∆p1 = ∆λ 1 δ pt 0
(3.45)
Neste caso, δ pt 0 é o vetor iterativo de deslocamento tangencial no início do processo
incremental.
A Equação 3.45 tem que ser restringida pelo comprimento de corda incremental ∆l.
Para tanto, deve-se substitui-la na Equação 3.19, resultando a Equação 3.46.
∆l
∆l
∆λ 1 = ± √
= s√
T
T
(δ pt 0 ) δ pt 0
(δ pt 0 ) δ pt 0
(3.46)
onde, s representa o sinal (+ ou −) que acompanhará a equação.
É importante observar na Equação 3.46 a possibilidade de se terum sinal positivo ou
negativo. Isto significa que existem duas soluções que possuem o mesmo valor, porém
as suas direções são opostas.
Apenas uma dessas soluções seguirá a trajetória de equilíbrio correta. A outra solução levará à falência do método. Desta forma, é extremamente importante escolher
corretamente o sinal desta equação.
3.3.1
Sinal do Determinante da Matriz de Rigidez
Diferentes critérios são encontrados na literatura para determinar a direção correta
da trajetória de equilíbrio. O método de escolha do sinal apropriado mais utilizado foi
originalmente introduzido por Crisfield (1981).
Neste método, a escolha do sinal s para a Equação 3.46 é fundamentada na definição
da matriz de rigidez tangente, do seguinte modo
s(∆λ 1 ) = s(det(Kt ))
(3.47)
onde Kt é a matriz de rigidez tangente da estrutura no início do incremento, quando ∆λ 1
é calculado.
3.4 Comprimento de Corda Incremental
35
A Equação 3.47 determina que o sinal de ∆λ 1 será sempre positivo quando o determinante de Kt for positivo, o que significa que a matriz de rigidez é positiva definida
(estrutura está em equilíbrio estável). Do mesmo modo, quando o determinante de Kt
for negativo (estrutura está em equilíbrio instável) o sinal de ∆λ 1 será negativo.
Além do mais, neste critério, o sinal do determinante poderá ser facilmente calculado através da análise do sinal dos termos da matriz diagonal D, originada do processo
de fatorização LDLT da matriz de rigidez Kt . Esta analise é feita do seguinte modo
{
todos D > 0 → s = +1
(3.48)
ao menos um D < 0 → s = −1
Este método responde bem na presença de pontos limites. No entanto, na presença
de pontos de bifurcação o sinal geralmente fica oscilando entre positivo e negativo em
sua vizinhança.
3.4
Comprimento de Corda Incremental
A determinação do tamanho apropriado para o comprimento de corda incremental
∆l ainda é um procedimento subjetivo. Por conta disso, seu tamanho ideal depende de
alguns fatores do comportamento da trajetória de equilíbrio, tais como o grau de não
linearidade da curvatura da trajetória de equilíbrio e a posição dos pontos limites e de
bifurcação.
Como ainda não existe uma regra geral para definir comprimento de corda ideal, a
primeira estimativa depende de um processo de tentativa e erro.
O comprimento de corda pode ser considerado fixo, onde o valor é constante ao
longo de todo o processo do Método da Corda, ou variável, no qual seu valor varia
durante o processo incremental.
A magnitude do primeiro comprimento de corda incremental (∆l) pode ser um valor
arbitrado pelo analista antes do início do processo ou esta magnitude apropriada pode
3.4 Comprimento de Corda Incremental
ser calculada no início do processo de acordo com a Equação 3.49.
√
T
∆l = (δ pt 0 ) δ pt 0
36
(3.49)
onde δ pt 0 é o vetor de deslocamento iterativo calculado no início do processo do Método da Corda, o que justifica o índice sobreposto 0.
Para os demais passos incrementais o comprimento de corda ∆l, caso seja variável,
pode ser determinado de forma automática durante o processo incremental através do
procedimento chamado de passo de dimensionamento automático 10 .
Este procedimento tem como finalidade mensurar o grau de não linearidade da trajetória de equilíbrio. Esta técnica leva ao fornecimento de pequenos incrementos quando
a resposta é mais não-linear e de grandes incrementos quando a resposta é mais linear.
O algoritmo deste procedimento ajusta o tamanho do próximo incremento através
de uma relação entre o comprimento de corda incremental no passo atual e número de
iterações necessários para alcançar a convergência. Esta técnica foi desenvolvida por
Crisfield (1981) e é dada pela Equação 3.50.
( )n
( )1/2
Id
Id
∆li = ∆l0
= ∆l0
I0
I0
(3.50)
onde ∆li é o novo comprimento de corda incremental, ∆l0 é o comprimento de corda incremental no passo atual, I0 é número de iterações requeridas e Id o número de iterações
desejadas. Nota-se que estas duas últimas variáveis são definidas pelo usuário no início
do processo. Ainda na Equação 3.50, observa-se que o parâmetro n foi considerado
como 0, 5 por Ramm (1981), com o intuito de evitar que ocorram rápidas oscilações no
comprimento de corda incremental ∆l.
Caso raízes complexas sejam encontradas o comprimento de corda deve ser reduzido e o processo do Método da Corda reiniciado. O procedimento utilizado para reduzir
este comprimento de corda é chamado de corte automático de incremento 11 , e é dado
10 Termo
traduzido livremente do original em inglês: Automatic step sizing.
livremente do inglês: automatic increment cutting.
11 Traduzido
3.4 Comprimento de Corda Incremental
37
pela Equação 3.51.
∆l i+1 = β ∆l i
(3.51)
onde, β é um escalar com valor definido entre
0, 1 ≤ β ≥ 0, 5
(3.52)
O corte automático de incremento também é utilizado quando o sistema não alcança
a convergência após atingir o número máximo de iterações permitidas, cujo valor é
adotado no início do processo.
38
4
Otimização do Método da Corda
Nos Capítulos anteriores foram descritas e formuladas todas as etapas referentes
ao Método da Corda. Foi visto que existem alguns pontos que podem dificultar o bom
funcionamento deste método. Por conta disso, e o objetivo deste Capítulo consiste
sistematicamente em apresentar soluções para tais imperfeições.
Como já foi supracitado na Seção 3.4, ainda não existe na literatura um método
específico com a finalidade de determinar e controlar o tamanho apropriado para o comprimento de corda ao longo do passo incremental no Método da Corda.
Este é um critério muito importante para o desenvolvimento dos métodos de continuação, pois o tamanho do comprimento de corda irá afetar todos os passos necessários para circunscrever a trajetória de solução considerada. Além disso, a cada passo
incremental, este também afetará o número de iterações de equilíbrio requeridas para
alcançar convergência no passo corretor e, ainda, poderá afetar o custo total para traçar
toda a trajetória.
As técnicas apresentadas na Seção 3.4 para determinar e controlar o tamanho do
comprimento de corda levam a procedimentos que trabalham razoavelmente bem para
problemas com suaves não linearidades. Todavia, apresentam algumas dificuldades em
estruturas que possuem mudanças repentinas e fortes não linearidades em suas trajetórias de equilíbrio. Também é importante ressaltar que, caso o comprimento de corda
seja mal dimensionado, há possibilidade de raízes complexas aparecerem.
Um dos pontos chave em um algoritmo de Método da Corda diz respeito à escolha correta do sinal do parâmetro de carregamento incremental na primeira iteração em
4 Otimização do Método da Corda
39
cada passo incremental-iterativo. Entretanto, conforme a Seção 3.3, existem vários critérios propostos na literatura para tal determinação, porém nenhum com desempenho
satisfatório.
A escolha incorreta deste sinal inevitavelmente levará a uma resposta indesejada e,
consequentemente, causará a falência do método ou por divergência no passo corretor
ou pelo retorno da trajetória (doubling back) em pontos previamente calculados. Isto
normalmente acontece na presença de pontos limites e de bifurcação, principalmente
quando as estruturas possuem comportamentos fortemente não lineares.
Finalmente, quanto aos problemas de retorno da trajetória (doubling back) em pontos já convergidos, estes podem ocorrer em qualquer estágio de um método de continuação e caracterizam-se como o "fenômeno"mais indesejável que pode ocorrer em um
método de continuação.
Este fenômeno ocorre quando, por alguma razão, uma solução atual (representada
por esta letra k), converge para uma solução já calculada anteriormente (k − 1), ao invés
de convergir para a próxima solução desejada (k + 1). Quando isto ocorre não é possível
traçar a trajetória de equilíbrio completa da estrutura em questão.
Em geral, o doubling back pode ocorrer apenas se uma estimativa inicial ruim é
construída no estágio previsor do método, onde o critério para a escolha do incremento
de carregamento inicial correto não funciona apropriadamente. Dessa forma, o doubling
back pode ser evitado se um critério correto for válido.
Na realidade, a ocorrência dos doubling backs se dão através da deficiência entre
o parâmetro de trajetória usado e o verdadeiro comprimento de corda, o qual aumenta
quando as características não lineares se tornam mais intensas.
A única maneira para diminuir esta deficiência é adaptar uma estratégia própria de
controle do comprimento de corda, tal que a mudança entre quaisquer duas soluções
adjacentes não sejam muito acentuadas. A escolha de um preditor tangente próximo à
solução desejada garantirá encontrar um corretor adequado.
Para que os problemas descritos acima sejam resolvidos, serão propostos, nas pró-
4.1 Método da Corda Acumulado
40
ximas seções, dois métodos que serão acrescentados ao Método da Corda. Um método
define um comprimento de corda incremental e o outro propõe um novo critério de
determinação do sinal do previsor.
Primeiramente, será proposto o Método da Corda Acumulado, desenvolvido por
Teng e Luo (1998). Este método possui parâmetros de controle do comprimento de
corda e é muito estável, mesmo em estruturas com comportamentos estruturais fortemente não lineares.
O segundo método foi desenvolvido por Feng et al. (1995) e diferencia-se dos demais por determinar o sinal apropriado (positivo ou negativo), levando em consideração
o histórico do comportamento da trajetória de equilíbrio em questão.
Ao modificar o Método da Corda utilizando estes dois métodos, os problemas acima
citados serão resolvidos. Este novo método possui maior estabilidade de execução e garante que a trajetória de equilíbrio de uma estrutura seja sempre obtida, não importando
qual o seu grau de não linearidade.
4.1
Método da Corda Acumulado
O Método da Corda Acumulado é uma forma melhorada do método da Corda, de
forma que ele possa alcançar a convergência para fatores de carga predefinidos, pontos de bifurcação previstos ou valores predefinidos de qualquer outros parâmetros que
estejam em função do carregamento.
Com o objetivo de aumentar o controle do usuário sobre o Método da Corda, Teng e
Lou (1997) propuseram que um novo parâmetro fosse acrescentado ao método da corda
convencional. Este novo parâmetro foi chamado de comprimento de corda acumulado
1.
Ao agregar este parâmetro ao Método da Corda, este passou a se denominar Método
da Corda Acumulado 2 .
1 Na
2 Na
nomenclatura inglesa, diz-se accumulated arc length.
nomenclatura inglesa: Accumulated Arc-Length Method.
4.1 Método da Corda Acumulado
41
A teoria utilizada por Teng e Lou (1997) foi primeiramente proposta por Riks (1984)
para o cálculo de pontos críticos em uma trajetória de equilíbrio.
A principal motivação para utilizar o comprimento de corda acumulado neste trabalho, se baseia no ganho de eficiência para a convergência em pontos críticos. Além
disso, com o advento deste parâmetro, é eliminada a possibilidade de retorno da trajetória em pontos que já foram convergidos anteriormente.
Quando não existe a necessidade de monitorar o comportamento de uma determinada estrutura para um determinado fator de carga, seja para ter acesso às características
de uma estrutura para um determinado nível de carga ou para que o Método da Corda
Acumulado inicie desde o primeiro passo iterativo, o valor para o fator de carregamento
predefinido pode ser considerado nulo.
4.1.1
O Comprimento da Corda Acumulado
A idéia essencial do Método da Corda Acumulado é, inicialmente, utilizar o Método
da Corda Convencional para traçar a trajetória de equilíbrio de uma determinada estrutura. Durante esta execução, um monitoramento contínuo do processo é realizado para
observar o fator de carga convergido. Este monitoramento é uma poderosa ferramenta
de controle do usuário ao processo. Quando este fator convergido se aproximar do fator
de carga predefinido, o Método da Corda Acumulado é imediatamente acionado.
Para incluir este processo no método da corda convencional, ele deve ser modificado. Isso se torna possível através do conceito do comprimento de corda acumulado.
O comprimento da corda acumulado é definido como a soma de todos os incrementos de comprimento do corda ∆l calculados até o momento presente, incluindo o fator
de carga atual. Sua representação é ilustrada pela Equação 4.1.
i
Li =
∑ ∆lk
(4.1)
k=1
O parâmetro Li é o comprimento do arco acumulado no passo de carga atual i e ∆lk
4.1 Método da Corda Acumulado
42
é o comprimento da corda incremental do método da corda convencional no passo de
carga k.
O parâmetro L representa o estado atual da estrutura e depende não apenas das
características da estrutura e seus carregamentos, mas também do processo incremental
de carga durante a análise.
Para que seja possível modificar o método da corda convencional, um novo parâmetro γ é definido como
γ (Li ) = λi − λd
(4.2)
onde λi é o fator de carregamento convergido (no método da corda convencional) no
passo de carregamento incremental i e λd é o fator de carregamento predefinido
3
que
será introduzido pelo analista no início do processo.
Nota-se que, por simplicidade de representação na formulação do Método da Corda
Acumulado proposto por Teng e Luo (1998), a convergência é requerida apenas para
um fator de carga predefinido. Em uma única análise, vários fatores de carga poderão
ser predefinidos.
O processo relativo ao Método da Corda Acumulado é iniciado quando, através da
escolha de um critério de aproximação racional, o fator de carregamento convergido
(λi ) aproxima-se do fator de carregamento predefinido (λd ). Contudo, neste trabalho,
pode se assumir que o valor para o fator de carregamento predefinido seja negligenciado
(λd = 0), fazendo com que o Método Acumulado se inicie junto com o Método da Corda
Convencional.
Para que o Método da Corda convencional seja modificado, calcula-se um comprimento da corda desejado, representado por ∆ld , para o passo de carga incremental
seguinte (i + 1). Para isso, o comprimento de corda acumulado desejado, Ld , deve satisfazer a Equação 4.3.
γ (Ld ) = 0
3O
índice d deriva-se da palavra inglesa "desired", porém aqui será utilizado como predefinido.
(4.3)
4.1 Método da Corda Acumulado
43
A Equação 4.4 define o comprimento de corda acumulado desejado.
Ld = Li+1 = Li + ∆li+1 = Li + ∆ld
(4.4)
onde, ∆ld representa o comprimento de corda desejado para o próximo passo incremental.
Utilizando a expressão de expansão da Série de Taylor de forma truncada para a
Equação 4.3, temos
γ (Ld ) = γ (Li + ∆ld ) = γ (Li ) +
d γ (Li )
1 d 2 γ (Li ) 2
∆ld +
∆ld + ... = 0
dL
2 dL2
(4.5)
A Equação 4.5 permite encontrar uma aproximação para o comprimento do corda
incremental desejado (∆ld ). Neste trabalho serão apresentadas três formas distintas de
aproximação para este comprimento, todas estas foram também propostas por Teng e
Lou (1997).
4.1.2
Aproximação Quadrática
Nesta aproximação, apenas os termos truncados até a segunda ordem da Equação
4.5 são considerados, sendo que os termos de terceira ordem e superiores são omitidos.
Portanto, a equação quadrática é, deste modo, ilustrada pela Equação 4.6.
γ (Li ) +
d γ (Li )
1 d 2 γ (Li ) 2
∆ld +
∆ld = 0
dL
2 dL2
(4.6)
As derivadas de γ em relação à L, da Equação 4.6, são obtidas utilizando aproximação por diferenças finitas regressivas.
d γ (Li ) γ (Li ) − γ (Li−1 )
=
dL
Li − Li−1
(4.7)
4.1 Método da Corda Acumulado
44
/
/
d 2 γ (Li ) d γ (Li ) dL − d γ (Li−1 ) dL
=
dL2
Li − Li−1
(4.8)
d γ (Li−1 ) γ (Li−1 ) − γ (Li−2 )
=
dL
Li−1 − Li−2
(4.9)
As derivadas referem-se ao passo de carregamento incremental atual (i) e ao passo
de carga incremental anterior a este (i − 1).
A substituição de 4.2 em 4.7 e 4.9, respectivamente, leva às Equações 4.10 e 4.11.
d γ (Li ) λi − λi−1
=
dL
Li − Li−1
(4.10)
d γ (Li−1 ) λi−1 − λi−2
=
dL
Li−1 − Li−2
(4.11)
Para o desenvolvimento da Equação 4.8, considera-se os resultados das derivadas
de primeira ordem obtidos pelas Equações 4.10 e 4.11, que são substituídas na Equação
4.8, resultando na Equação 4.12.
d 2 γ (Li )
λi − λi−1
λi−1 − λi−2
=
−
2
2
dL
(Li−1 − Li−2 ) (Li − Li−1 )
(Li − Li−1 )
(4.12)
A expressão final da aproximação quadrática é obtida introduzindo os resultados
das derivadas de primeira e segunda ordem dadas pelas Equações 4.10, 4.11 e 4.12,
respectivamente, na Equação 4.6. Com isso,
b1 ∆ld 2 + b2 ∆ld + b3 = 0
(4.13)
Os coeficientes b1 , b2 e b3 são definidos pelas Equações 4.14, 4.15 e 4.16, respectiviamente.
(
1
b1 =
2
λi − λi−1
λi−1 − λi−2
−
2
(Li−1 − Li−2 ) (Li − Li−1 )
(Li − Li−1 )
)
(4.14)
4.1 Método da Corda Acumulado
45
(
b2 =
λi − λi−1
Li − Li−1
)
b3 = (λi − λd )
(4.15)
(4.16)
A Equação 4.13 permite encontrar o novo comprimento de corda incremental do
passo (i + 1) para alcançar o nível de carga predefinido. Contudo, por ser uma equação
quadrática, existem duas raízes possíveis para ∆ld . A raiz apropriada é aquela que mais
se aproxima da solução linear. Para tanto, comparam-se as duas raízes, em termos de
valores absolutos, e a raiz apropriada é, em geral, muito menor do que a outra raiz.
Observa-se que os valores absolutos são utilizados, pois o sinal de ∆ld pode ser positivo
ou negativo dependendo da localização do ponto desejado.
Raízes imaginárias (ou complexas), teoricamente, podem ser obtidas em alguns
/
casos especiais, quando o valor de d γ (Li ) dL se aproxima de zero. Não obstante, conforme os autores Teng e Luo (1998), esta possibilidade raramente ocorre em cálculos
numéricos. Em suas simulações numéricas este tipo de problema nunca ocorreu.
4.1.3 Aproximação Linear
A aproximação linear do comprimento do arco incremental desejado, ∆ld , é encontrada de forma análoga à aproximação quadrática (Seção 4.1.2). No entanto, consideramse apenas os termos lineares da Equação 4.5. Com isso, tem-se que
γ (Li ) +
d γ (Li )
∆ld = 0
dL
(4.17)
Substituindo as Equações 4.10 e 4.2 na Equação 4.17, obtém-se a seguinte expressão
)
λi − λi−1
(λi − λd ) +
∆ld = 0
Li − Li−1
(
(4.18)
4.1 Método da Corda Acumulado
46
Afim de simplificar a Equação, isola-se o ∆ld na Equação 4.18, resultando em
∆ld =
λd − λi
λd − λi
(Li − Li−1 ) =
∆li
λi − λi−1
λi − λi−1
(4.19)
Finalmente, a Equação 4.19 ilustra é a forma final da equação de aproximação linear
de ∆ld . Em consequência de sua linearidade, apenas um valor para ∆ld será encontrado.
A aproximação linear é uma opção vantajosa, devido à sua maior estabilidade por
determinar apenas uma raiz. Com isso, menos esforço computacional será demandado.
A aproximação linear tem como principal vantagem a não geração de raízes complexas
na solução do problema.
4.1.4
Aproximação Quadrática Alternativa
Teng e Luo (1998) propuseram uma formulação alternativa à aproximação quadrática desenvolvida na Seção 4.1.2, considerando-a mais vantajosa, pois há escolha da raiz
apropriada e não há possibilidade de gerar raízes complexas.
Nesta formulação, o comprimento de corda acumulado (L) é definido em função do
fator de carga, como mostra a Equação 4.20.
L = L (λ )
(4.20)
Após a convergência do fator de carga λi , no final do passo de carregamento i, o
comprimento de corda acumulado é igual a Li .
Adicionalmente, expande-se o comprimento de corda acumulado desejado (Ld ) em
função do comprimento de corda atual (Li ) através da expressão de expansão da Série
de Taylor de forma truncada.
Ld = L (λd ) = L (λi + ∆λd ) = Li +
dLi
1 d 2 Li
∆λd +
∆λd 2 + ...
dλ
2 dλ 2
(4.21)
onde
∆λd = λd + λi
(4.22)
4.1 Método da Corda Acumulado
47
Considerando apenas os termos de primeira e segunda ordem da Equação 4.21, a expressão para o comprimento da corda incremental desejado (∆λd ) é dada pela Equação
4.23.
∆ld = Ld − Li =
dLi
1 d 2 Li
∆λd +
∆λd 2
dλ
2 dλ 2
(4.23)
Novamente, para resolver as derivadas de L com relação a λ , utiliza-se a expressão
de diferenças finitas regressivas.
dLi Li − Li−1
=
dλ
λi − λi−1
(4.24)
/
/
d 2 Li dLi d λ − dLi−1 d λ
=
dλ 2
λi − λi−1
(4.25)
dLi−1 Li−1 − Li−2
=
dλ
λi−1 − λi−2
(4.26)
No entanto, substituindo as Equações 4.24 e 4.26 na Equação 4.25, tem-se
Li−1 − Li−2
d 2 Li
Li − Li−1
=
−
2
2
dλ
(λi − λi−1 ) (λi−1 − λi−2 )
(λi − λi−1 )
(4.27)
Finalmente, a Equação para a aproximação alternativa é obtida substituindo as
Equações 4.24 e 4.27 na Equação 4.23.
)
(
Li − Li−1
Li−1 − Li−2
Li − Li−1
1
∆λd 2
∆λd +
∆ld =
−
λi − λi−1
2 (λi − λi−1 )2 (λi − λi−1 ) (λi−1 − λi−2 )
(4.28)
Com isso, conclui-se que, ao utilizar o comprimento da corda acumulado (L) em
função do fator de carga (λ ), a equação final não se torna quadrática, pois o valor do
fator de carga (λd ) é conhecido.
Nota-se que, se apenas os termos lineares da Equação 4.23 forem considerados,
encontra-se uma aproximação linear para o comprimento do corda incremental idêntica
à Equação 4.19 da Seção 4.1.3.
4.1 Método da Corda Acumulado
4.1.5
48
Parâmetros de Controle
O Método da Corda Acumulado completa-se com a adição de dois novos parâmetros
que controlam o seu funcionamento. Estes, que são nomeados de parâmetros de controle, são de essencial importância para o método e desempenham as seguintes funções:
controlar o comprimento da corda incremental e controlar o fator de carga convergido.
Como já mencionado anteriormente, o Método da Corda Acumulado só é acionado
através da análise da magnitude entre o fator de carga convergido e o fator de carga
predefinido. Para tanto, adiciona-se um critério de aproximação (α ) entre estes fatores,
critério este que deve ser imposto pelo usuário no início do método.
O parâmetro de controle de carregamento determina quando o Método da Corda
Acumulado será ativado, através do monitoramento, a cada passo incremental, do fator
de carga convergido.
λd − λi ≤ α
(4.29)
A análise da diferença entre o fator de carga convergido (λi ) e o fator de carga
predefinido (λd ) impõe que, quando esta diferença for menor ou igual ao critério de
aproximação α , o Método da Corda Acumulado deve ser acionado.
Quando este for acionado, será calculado o comprimento de corda desejado através
de uma das aproximações mencionadas nas seções anteriores. Consequentemente, haverá dois comprimentos de corda diferentes, um dado pelo método convencional (∆l) e
o outro pelo método acumulado (∆ld ).
Portanto, faz-se necessário adicionar outro parâmetro de controle ao processo para
determinar qual o comprimento de corda mais apropriado para ser utilizado na continuação do método.
E este é designado por parâmetro de controle de comprimento de corda, que por
sua vez, analisa a magnitude de ∆ld e ∆l e determina qual o comprimento de corda que
4.2 Novo Método Previsor
49
será utilizado para o próximo passo, do seguinte modo,
∆l < ∆ld → ∆l ou ∆l > ∆ld → ∆ld
(4.30)
O comprimento de corda ideal será, sempre, o de menor magnitude entre os dois
comprimentos encontrados. Isso se deve ao fato de que grandes comprimentos de corda
podem levar o método mais facilmente à divergência.
A cada passo incremental, o parâmetro de controle de carga será acionado e, caso o
Método da Corda Acumulado seja solicitado, o parâmetro de controle de comprimento
de corda também será utilizado. Deste modo, observa-se que os parâmetros de controle
são essenciais para o funcionamento do Método da Corda Acumulado.
O Método da Corda Acumulado também pode ser utilizado para alcançar convergência de valores predefinidos de outros parâmetros, como por exemplo, deslocamentos
ou tensões. Logo, troca-se apenas as variáveis de carregamento (λ ) pelas variáveis correspondentes ao parâmetro que será analisado.
4.2
Novo Método Previsor
Nesta seção, será apresentado o critério proposto por Feng et al. (1995) para determinar o sinal correto para o fator de carregamento incremental. Este critério produz
uma previsão muito confiável para a próxima direção que a trajetória de equilíbrio deverá seguir.
Apenas por conveniência, a Equação 3.46, que calcula o incremento de carga para
o próximo passo, será reescrita aqui.
∆l
∆λ 1 = s √
t
(δ pt 0 ) δ pt 0
(4.31)
A letra s representa o sinal (+ ou −) que acompanhará a equação. Portanto, esta
equação terá duas soluções possíveis que, por sua vez, definem duas direções tangentes
4.2 Novo Método Previsor
50
e opostas uma à outra e, com uma delas, o procedimento retornará à solução anterior.
Como já mencionado, o critério mais utilizado para determinar o sinal do fator de
carga é dado pelo determinante da matriz de rigidez da estrutura.
s(∆λ 1 ) = s(det(Kt ))
(4.32)
Contudo, o método proposto pela Equação 4.32 só funciona adequadamente no
Método da Corda a partir do início do seu processamento até o instante em que o ponto
de bifurcação é encontrado.
Por conta disso, um novo critério foi desenvolvido por Feng et al. (1995, 1996),
para a escolha do sinal adequando para o fator de carga incremental.
4.2.1
Critério do Produto Interno
Segundo Neto e Feng (1999), a maioria dos critérios desenvolvidos para determinar
o sinal do fator de carga incremental dependem exclusivamente de informações relativas
ao ponto de equilíbrio atual (no início do incremento). Consequentemente, a escolha do
sinal é feita sem considerar o histórico da trajetória de equilíbrio atual.
Em situações nas quais as estruturas possuem comportamentos fortemente não lineares, a escolha de tais critérios resultam na predição errada da direção das trajetórias
de equilíbrio para o primeiro passo incremental. Diante disso, Feng et al. (1995, 1996)
propuseram um critério diferente para a determinação do sinal de ∆λ ao considerar a
história da trajetória de equilíbrio atual da estrutura, evitando o erro de predição na
direção das trajetórias de equilíbrio.
Neste critério, o sinal do fator de carga previsto deve coincidir com o sinal do produto interno 4 entre o vetor de deslocamento incremental convergido no passo anterior,
∆p, e o vetor de deslocamento tangencial iterativo do passo atual, δ pt , como mostra a
4 Produto
Interno é uma operação algébrica entre dois vetores.
4.2 Novo Método Previsor
51
Equação 4.33.
sinal(∆λ 1 ) = sinal(∆pT δ pt )
(4.33)
Desta forma, este método foi nomeado de critério do produto interno. O ponto fundamental relativo a este critério se deve ao fato de que o ∆p carrega consigo informações
sobre a história da trajetória de equilíbrio atual.
A fim de continuar traçando a trajetória atual sem retornar aos pontos previamente
convergidos, a ideia básica deste critério requer que o vetor de deslocamento do incremento inicial esteja sempre em direção crescente na trajetória, obedecendo a condição
mostrada pela Equação 4.34.
∆λ ∆pT δ pt > 0
(4.34)
Além disso, o produto positivo entre ∆pT e δ pt indica que o fator de carga está,
a cada nova iteração, aumentando de valor, de forma que o sinal positivo deva ser escolhido para o previsor ∆λ 1 continuar seguindo a trajetória de equilíbrio atual. Similarmente, caso o produto entre ∆pT e δ pt for negativo, o valor do fator de carga irá
decrescer e, assim, o sinal negativo deve ser escolhido para ∆λ 1 .
A importância de se considerar o histórico da trajetória de equilíbrio pode ser provada através de argumentos geométricos. Para tanto, sabe-se que o vetor de deslocamento tangente δ pt já definido pela Equação ??b, será repetido aqui, é dado pela
Equação 4.35.
δ pt = (Kt )−1 qef
(4.35)
O vetor δ pt é um vetor tangente à trajetória de equilíbrio no eixo dos deslocamentos.
Este vetor aponta para a direção do gradiente positivo de λ , fornecendo a informação
da direção ao longo da qual o fator de carga aumenta. Isso mostra que esta direção está
associada com a escolha positiva de λ . Desta forma, o vetor δ pt indica a nova direção
da trajetória de equilíbrio atual.
Por outro lado, o vetor de deslocamento incremental ∆p é um vetor secante à trajetória de equilíbrio e, quando suficientemente pequeno, dá uma aproximação satisfatória
4.2 Novo Método Previsor
52
da curva de solução dentro do intervalo. O valor de ∆p, deste modo, se aproxima da
história da trajetória de equilíbrio dentro de um dado intervalo.
O vetor δ pt precisa ser suficientemente pequeno para determinar precisamente a
próxima direção da trajetória, esta necessidade sempre será mantida, pois uma limitação
para o tamanho máximo deste vetor é imposta pelo comprimento de corda incremental.
Além do mais, observa-se que critério desenvolvido por Feng torna-se insensível à
presença de pontos de bifurcação e de pontos limites, assim como as dificuldades associadas ao critério baseado no sinal do determinante de Kt e as trajetórias com acentuados
snap-backs.
53
5
Simulação Numérica
Neste capítulo, serão apresentadas simulações numéricas computacionais necessárias para validar os métodos propostos e formulados neste trabalho. Para tais simulações
utilizam-se estruturas em formato de treliça discretizadas através da Teoria Geometricamente Exata, proposta no Capítulo 2.
As estruturas aqui utilizadas foram desenhadas utilizando um programa de pré e pós
processamento em elementos finitos chamado GID. Este programa é ideal para gerar
toda a informação requerida através de uma análise numérica estrutural porém, neste
trabalho, utiliza-se apenas as ferramentas de geração de malha para variadas estruturas.
Em seguida, foi criada uma interface para a transferência de dados entre o GID
e o programa do Método da Corda desenvolvido em ambiente Matlab, utilizando um
problem type chamado MAT-fem.
5.1
Métodos Previsores
O desempenho dos métodos previsores serão provados por meio de uma análise
comparativa dos seus comportamentos utilizando o Método da Corda Cilíndrico.
5.1.1
Treliça de duas barras geometricamente simétricas
A primeira estrutura utilizada é a treliça simples de duas barras geometricamente
simétricas, como mostra a Figura 5.1. Esta treliça foi discretizada em dois elementos
5.1 Métodos Previsores
54
de barra e possui material de comportamento elástico linear, cujo o Módulo de Elasticidade é igual a 2.1 × 105 . A seção transversal da barra é igual a 20 e as demais
características geométricas são mostradas na Figura 5.1. Esta estrutura também possui
uma carga unitária, aplicada ao seu quarto grau de liberdade global, que será aumentada
incrementalmente.
Figura 5.1: Treliça com duas barras simétricas e carga aplicada
Na primeira simulação realizada, utilizou-se o método previsor baseado no sinal
do determinante da matriz de rigidez da estrutura Kt (Seção 3.3.1). A trajetória de
equilíbrio desta estrutura encontra-se na Figura 5.2.
Figura 5.2: Diagrama de carga-deformação da treliça utilizando o previsor det(Kt ).
5.1 Métodos Previsores
55
Nesta figura, também foi plotado o determinante de Kt para uma avaliação do seu
comportamento ao longo do processo.
Em seguida, utilizou-se o método previsor chamado de Critério do Produto Interno,
apresentado na Seção 4.2. A Figura 5.3 ilustra o diagrama de carga-deformação utilizando este previsor e também o determinante de Kt .
Figura 5.3: Diagrama de carga-deformação da treliça utilizando o Critério do Produto
Interno.
Observa-se na Figura 5.4, que a estrutura comportou-se de forma muito uniforme
ao utilizar os dois métodos, de modo que a trajetória de equilíbrio convergiu exatamente
nos mesmos pontos.
Deste modo, pode-se concluir que este exemplo não é indicado para a avaliação
de desempenho dos Métodos Previsores por não apresentar situações fortemente não
lineares.
5.1 Métodos Previsores
56
Figura 5.4: Comparação dos Diagramas carga-deformação da treliça utilizando os dois
métodos.
5.1.2
Domo em forma de estrela
A segunda estrutura se refere a um domo em treliça espacial com o formato de uma
estrela, como mostra a Figura 5.5. Este domo já foi proposto na literatura por Crisfield
(1997) e Wriggers (2008) e possui, como característica, um comportamento fortemente
não linear.
Este domo possui material elástico não linear de St. Venant com Módulo de Young
igual a 1079.6. Quanto a sua geometria, os nós externos do domo estão localizados
em um círculo de raio igual a 50 e os nós internos em um círculo de raio igual a 25.
Os nós internos estão localizados a uma altura de 6.216, enquanto os nós do centro do
domo estão a uma altura de 8.216. Nesta estrutura, os nós externos estão simplesmente
apoiados e o ponto de carga unitária aplicada situa-se no topo do Domo, conforme a
Figura (5.5).
Assim como no exemplo anterior, foi aplicado, primeiramente, o Método da Corda
5.1 Métodos Previsores
57
Figura 5.5: Domo tridimensional em forma de estrela.
Cilíndrico, utilizando o previsor baseado no sinal do determinante de Kt . O comportamento desta estrutura encontra-se na Figura 5.6.
Figura 5.6: Diagrama de carga-deformação do Domo utilizando o previsor det(Kt ).
Observa-se que, circunstancialmente, não foi possível completar a trajetória de
equilíbrio completa desta estrutura. Após o incremento de número 83 o Método da
Corda retorna ao passo incremental anterior e fica oscilando entre os passos 82 e 83 até
o término do Método.
5.2 Métodos Corretores
58
Em seguida, foi aplicado o Método da Corda Cilíndrico utilizando o Critério do
Produto Interno proposto por Feng et al. (1995). A trajetória de equilíbrio desta estrutura
é mostrada pela Figura 5.7.
Figura 5.7: Diagrama de carga-deformação utilizando o Critério do Produto Interno.
Desta vez, a trajetória de equilíbrio foi completamente traçada, utilizando para isso
um total de 315 passos incrementais. A Figura 5.7 também mostra a trajetória de comportamento do determinante de Kt .
Nota-se que a trajetória de equilíbrio do domo possui um comportamento fortemente não linear com vários snap-through e snap-backs. Devido a este fato, o método
previsor que utiliza o sinal o determinante de Kt não consegue traçar completamente a
trajetória.
5.2
Métodos Corretores
Nesta seção, os Métodos Corretores descritos na Seção 3.2 terão seus desempenhos
testados e comparados. Duas estruturas com comportamentos diferentes serão utilizadas. Uma delas diz respeito a um Domo, em formato de estrela (já apresentado na Seção
5.2.1), e a outra corresponde a um arco circular bidimensional.
5.2 Métodos Corretores
5.2.1
59
Domo em forma de estrela
As propriedades materiais e geométricas do domo espacial, em formato de estrela
(Figura 5.5), foram apresentadas na Seção 5.2.1.
Primeiramente, foi aplicado a esta estrutura o Método da Corda com equação de
restrição cilíndrica desenvolvido por Crisfield (1981). O diagrama de carga-deformação
encontra-se na Figura 5.8.
Figura 5.8: Trajetória de equilíbrio com pontos convergidos do Domo utilizando o Método da Corda Cilíndrico.
Ao utilizar esta equação de restrição cilíndrica, a trajetória de equilíbrio da estrutura é traçada utilizando 315 passos incrementais. A cada passo incremental, são necessários, em média, 2 passos iterativos para alcançar a convergência. Esta relação de
incrementos-iterações é ilustrada pela Figura 5.9.
Do mesmo modo, aplicou-se nesta estrutura o Método Corretor desenvolvido por
Schweizerhof e Wriggers (1986). A trajetória de equilíbrio do domo é mostrada através
do diagrama de carga-deformação da Figura 5.10.
5.2 Métodos Corretores
60
Figura 5.9: Diagrama de convergência do Domo no Método da Corda Cilíndrico.
Figura 5.10: Trajetória de equilíbrio com pontos convergidos do Domo utilizando o
Método de Schweizerhof-Wriggers.
Quanto ao desempenho, para se obter a mesma trajetória de equilíbrio do Método
Cilíndrico, foram necessários 465 passos incrementais. A taxa de convergência por
passo incremental ficou com uma média de 3 passos iterativos, como ilustra a Figura
5.11.
5.2 Métodos Corretores
61
Figura 5.11: Diagrama de convergência do Domo no Método de SchweizerhofWriggers.
Na Figura 5.12, são assinalados alguns pontos convergidos da trajetória de equilíbrio. As configurações deformadas do domo, referente a esses pontos, são apresentadas
em vermelho nas Figuras de 5.13 a 5.17.
Figura 5.12: Trajetória de equilíbrio do domo com pontos de deformação.
5.2 Métodos Corretores
Figura 5.13: Estado deformado do domo no ponto A.
Figura 5.14: Estado deformado do domo no ponto B.
Figura 5.15: Estado deformado do domo no ponto C.
62
5.2 Métodos Corretores
63
Figura 5.16: Estado deformado do domo no ponto D.
Figura 5.17: Estado deformado do domo no final da trajetória de equilíbrio.
5.2.2
Arco Circular bidimensional
A segunda análise retrata uma estrutura de arco circular bidimensional em treliça
(Figura 5.18), já proposto na literatura por Crisfield (1997).
Este arco possui 101 barras, totalizando 42 nós. O Módulo de Young considerado
é igual a 1 × 107 . Quanto às características geométricas, as barras externas estão localizadas em um raio com valor igual a 50 e as barras internas em um raio igual a 48.
Este arco está engastado em suas duas extremidades e o ponto de carga unitária aplicada
5.2 Métodos Corretores
64
situa-se no centro do arco, conforme a Figura 5.18.
2
45°
R=48
Figura 5.18: Arco circular treliçado.
Do mesmo modo, aplicou-se a esta estrutura o Método da Corda com equação de
restrição cilíndrica. Consequentemente, foi obtido o diagrama de carga-deformação,
como ilustra a Figura 5.19.
Figura 5.19: Trajetória de equilíbrio do Arco utilizando o Método da Corda Cilíndrico.
Para traçar completamente a trajetória de equilíbrio da estrutura, foram necessários
670 passos incrementais. O diagrama carga-deformação com os pontos convergidos é
5.2 Métodos Corretores
65
mostrado na Figura 5.20. Além disso, foram utilizados, em média, 3 passos iterativos
para alcançar a convergência em cada incremento, como ilustra a Figura 5.21.
Figura 5.20: Trajetória de equilíbrio do Arco, com pontos convergidos, utilizando o
Método da Corda Cilíndrico.
Figura 5.21: Diagrama de convergência do Arco no Método da Corda Cilíndrico.
Ao se utilizar o Método Corretor através da equação de restrição de SchweizerhofWriggers, a trajetória de equilíbrio do arco é mostrada pela Figura 5.22.
5.2 Métodos Corretores
66
Figura 5.22: Diagrama de carga-deformação do Arco utilizando o Método de
Schweizerhof-Wriggers.
Quanto ao desempenho, foram necessários 692 passos incrementais para traçar totalmente a trajetória de equilíbrio, conforme a Figura 5.23.
Figura 5.23: Diagrama de carga-deformação do Arco, com pontos convergidos, utilizando o Método de Schweizerhof-Wriggers.
5.2 Métodos Corretores
67
A Figura 5.24 apresenta um diagrama com a relação entre os incrementos e as iterações convergidas. Foram necessários uma média de 4 passos iterativos para alcançar
a convergência em cada passo incremental.
Figura 5.24: Diagrama de convergência do Arco no Método de Schweizerhof-Wriggers.
Novamente, a configuração deformada da estrutura em alguns pontos convergidos
da trajetória de equilíbrio (Figura 5.25), são apresentadas pelas Figuras de 5.26 a 5.32.
Figura 5.25: Trajetória de equilíbrio do Arco com pontos de deformação.
5.2 Métodos Corretores
Figura 5.26: Estado deformado do arco no ponto A.
Figura 5.27: Estado deformado do arco no ponto B.
68
5.2 Métodos Corretores
Figura 5.28: Estado deformado do arco no ponto C.
Figura 5.29: Estado deformado do arco no ponto D.
69
5.2 Métodos Corretores
Figura 5.30: Estado deformado do arco no ponto E.
Figura 5.31: Estado deformado do arco no ponto F.
70
5.3 Método Acumulado
71
Figura 5.32: Estado deformado do arco no final do processo.
5.3
Método Acumulado
Nesta seção, será implementado o Método da Corda Acumulado descrito na Seção
4.1. Para tanto, é utilizada uma estrutura em treliça no formato de um cilindro. Esta
estrutura será pressionada por duas forças concentradas, de forma que esteja em condições que simulam o amassamento de uma lata de cerveja pelos dedos de uma pessoa
(Figura 5.33).
Este tipo de estrutura já foi utilizado por Campello (2005) e por Tiago (2007), com
a ressalva de ter sido empregado elementos de casca na modelagem. Este exemplo é
conhecido na literatura como pinched cylinder. Não obstante, foi levado ao contexto de
uma lata de cerveja por Campello (2005).
As propriedades geométricas do cilindro são apresentadas na Figura 5.33. A seção
transversal considerada é igual a 1 e o Módulo de Young adotado é igual a 3.0 × 104 .
Foram aplicadas duas forças opostas e concentradas no meio do cilindro e os graus de
liberdade das extremidades são restringidos. Além disso, o cilindro foi discretizado em
5.3 Método Acumulado
72
Figura 5.33: Propriedades geométricas do cilindro (CAMPELLO, 2005).
596 barras de treliça e possui 152 nós, como ilustra a Figura 5.34.
Figura 5.34: Barras de treliça do cilindro.
Em primeiro lugar, será comparado o desempenho do Método da Corda Acumulado
em relação ao Método da Corda Cilíndrico. Para tanto, a estrutura cilíndrica foi submetida à análise pelo Método da Corda Cilíndrico, cujo diagrama de carga-deformação é
mostrado pela Figura 5.35.
Do mesmo modo, a estrutura foi submetida à análise pelo Método da Corda Acumulado sem se considerar um carregamento predefinido, ou seja, foi adotado como λd = 0.
Com isso, o Método Acumulado se inicia logo no primeiro incremento. A Figura 5.36
5.3 Método Acumulado
73
Figura 5.35: Diagrama de carga-deformação do cilindro utilizando o Método da Corda
Cilíndrico.
ilustra a trajetória de equilíbrio da estrutura utilizando este método.
Quanto à análise de desempenho, o Método Cilíndrico necessita de menos passos
incrementais para traçar a trajetória de equilíbrio completa, totalizando 113 passos incrementais. Fazendo uma comparação, o Método Acumulado necessita de 125 passos
incrementais para traçar a trajetória até o mesmo ponto. Este último método necessita
de alguns passos incrementais a mais pelo fato do comprimento de corda ∆l calculado
ser ligeiramente menor do que o comprimento de corda calculado no Método Cilindrico.
Em contrapartida, o Método Acumulado possui uma taxa de convergência mais
eficiente do que o Método Cilíndrico, por necessitar de menos passos iterativos para
alcançar a convergência. Esta comparação pode ser vista através da Figura 5.37.
Deste modo, mesmo utilizando mais passos incrementais, o Método da Corda Acumulado continua sendo uma boa opção para aprimorar o Método da Corda Cilíndrico,
uma vez que, a diferença de incrementos não é expressiva para comprometer o desempenho do Método da Corda. Mesmo assim, proporciona uma maior estabilidade ao
5.3 Método Acumulado
74
Figura 5.36: Diagrama de carga-deformação do cilindro utilizando o Método da Corda
Acumulado (λd = 0).
Figura 5.37: Diagrama de comparação de convergência.
5.3 Método Acumulado
75
Método, melhorando a sua taxa de conversão.
A partir de então, será empregado apenas o Método da Corda Acumulado, pois será
testada a sua capacidade de convergir para um determinado carregamento predefinido.
Serão considerados 4 valores diferentes de carregamento desejado λd . Estes valores são:
100, 1000, 5000 e 10000. O critério de aproximação adotado é igual a α = 1.0 × 10−5 .
As trajetórias de equilíbrio da estrutura com convergência em cada um dos carregamentos predefinidos λd são apresentadas pelos diagramas de carga-deformação dados
pelas Figuras de 5.38 a 5.44, respectivamente.
As configurações deformadas da estrutura são apresentadas concomitantemente em
relação aos seus diagramas de carga-deformação correspondentes.
Figura 5.38: Diagrama de carga-deformação do cilindro para λd = 100.
5.3 Método Acumulado
Figura 5.39: Configuração deformada do cilindro em λd = 100.
Figura 5.40: Diagrama de carga-deformação do cilindro para λd = 1000.
76
5.3 Método Acumulado
Figura 5.41: Configuração deformada do cilindro em λd = 1000.
Figura 5.42: Diagrama de carga-deformação do cilindro para λd = 5000.
77
5.3 Método Acumulado
Figura 5.43: Configuração deformada do cilindro em λd = 5000.
Figura 5.44: Diagrama de carga-deformação do cilindro para λd = 10000.
78
5.3 Método Acumulado
79
Figura 5.45: Configuração deformada do cilindro em λd = 10000.
Observa-se que independente do tamanho do carregamento predefinido, esse método sempre irá alcançar a convergência neste.
Quanto aos diferentes métodos de aproximação desenvolvidos, diferentemente do
que Teng e Lou (1997) relataram, não apresentaram nenhuma diferença no desempenho do Método Acumulado. Deste modo, recomenda-se que o método de aproximação
linear seja utilizado ao invés dos demais para evitar a possibilidade de ocorrência de
futuros problemas na escolha das raízes.
A configuração deformada do cilindro é mostrada, por vários ângulos, através das
Figuras de 5.46 a 5.49. Essa configuração deformada ocorre tanto no final do processo
do Método da Corda Cilíndrico quanto do Acumulado.
5.3 Método Acumulado
Figura 5.46: Configuração atual e deformada do cilindro.
Figura 5.47: Configuração deformada do cilindro, visto pelo ângulo 1.
80
5.3 Método Acumulado
Figura 5.48: Configuração deformada do cilindro, visto pelo ângulo 2.
Figura 5.49: Configuração deformada do cilindro, visto pelo ângulo 3.
81
82
6
Conclusões Finais
Nessa dissertação de mestrado, os Métodos de Continuação denominados Métodos
da Corda foram minuciosamente estudados, com o objetivo de otimizar o seu funcionamento.
Foram formulados dois tipos de Métodos Corretores, o Método da Corda Cilíndrico,
desenvolvido por Crisfield (1991), e o Método da Corda de Schweizerhof-Wriggers
(SCHWEIZERHOF; WRIGGERS, 1986), cujos desempenhos foram comparados.
Além disso, também foram propostos dois Métodos Previsores, o Método do determinante da matriz de rigidez e o Critério do produto interno, desenvolvidos, respectivamente, por Crisfield (1981) e por Feng et al. (1996), Neto e Feng (1999).
Finalmente, o Método da Corda Cilíndrico foi modificado através da introdução
de um parâmetro de comprimento de corda acumulado desenvolvido por Teng e Lou
(1997), Teng e Luo (1998), sendo assim nomeado de Método da Corda Acumulado.
Estes métodos supracitados foram aplicados a estruturas treliças com comportamento não linear através de programação computacional e , deste modo, foi possível
obter as conclusões descritas abaixo.
Primeiramente, foi observado que, é de grande importância dimensionar cuidadosamente o tamanho do comprimento de corda que será considerado na análise. Pois,
foi observado, neste trabalho, que um comprimento de corda mal dimensionado pode
ocasionar um número maior de divergências nas iterações e aumentar a possibilidade
de ocorrência de raízes complexas. Quanto maior o comprimento de corda escolhido,
6 Conclusões Finais
83
maior será a probabilidade dos problemas acima ocorrerem. Deste modo conclui-se
que, a determinação correta do comprimento de corda incremental inicial possui uma
importância significativa para o sucesso dos Métodos da Corda.
Com relação aos dois métodos corretores apresentados, as simulações numéricas
mostraram que o Método da Corda Cilíndrico apresenta um melhor desempenho em
comparação com o Método da Corda de Schweizerhof e de Wriggers. Isso é evidenciado com um menor número de passos incrementais para traçar a trajetória de equilíbrio,
a taxa de convergência apresentada é melhor. Devido a necessidade de se satisfazer, a
cada iteração, a equação de restrição este método possui um maior domínio de atração
sobre os pontos da trajetória de equilíbrio. Como proposto por Schweizerhof e Wriggers (1986), o Método consistentemente linearizado pode ser adicionado ao Método da
Corda Cilíndrico para ser usado apenas quando as raízes complexas aparecerem.
O método previsor, baseado no sinal do determinante da matriz de rigidez, pode
ser utilizado apenas em estruturas simples que não possuam comportamento de snapback. Por outro lado, as simulações numéricas utilizando o método previsor do Produto
interno evidenciaram a sua grande eficácia e estabilidade, mesmo quando a estrutura
apresenta um comportamento fortemente não linear.
Em primeiro lugar, pode-se concluir que o Método da Corda Acumulado é um procedimento de fácil implementação no Método da Corda. Além do mais, neste trabalho,
consolidou-se como uma boa alternativa para o aprimoramento do Método da Corda,
uma vez que, apresentou uma taxa de convergência incremental melhor quando comparada com a taxa de convergência do Método da Corda Cilíndrico. Isto se deve ao
conceito de comprimento de corda introduzido através do parâmetro de comprimento de
corda acumulado proporcionou uma otimização no tamanho do comprimento de corda
incremental. Esta característica também é útil para a para prevenir que não ocorra de
retorno da trajetória a pontos já convergidos.
Finalmente, como opção de otimização do Método da Corda, sugere-se fortemente
utilizar o método previsor do Critério de Produto Interno e como método corretor o
Método da Corda Acumulado com aproximação linear em análises estruturais.
84
Referências
AL-RASBY, S. N. Solution techniques in nonlinear structural analysis. Computer and
Structures, v. 40, p. 985–993, 1991.
BATOZ, J. L.; DHATT, G. Incremental displacement algorithms for nonlinear problems.
Int. J. Num. Methods in Engineering, v. 14, p. 1262–1266, 1979.
BERGAN, P. G.; HORRIGMORE, G.; KRAKELAND, B.; SOREIDE, T. H. Solution
techniques for non-linear finite element problems. Int. J. Num. Methods in Engineering,
v. 12, p. 1677–1696, 1978.
BERGAN, P. G.; SOREIDE, T. H. Solution of large displacement and instability problems using the current stiffness parameter. In: BERGAN, P. G. (Ed.). Finite Element
in Nonlinear Mechanics. Trondheim: Tapir Press, 1978.
CAMPELLO, E. M. B. Modelos não-lineares de casca em elasticidade e elastoplasticidade com grandes deformações - teoria e implementação em elementos finitos. Tese
(Tese de Doutorado) — Universidade de São Paulo, 2005.
CRISFIELD, M. A. A fast incremental-iterative solution procedure that handles ’snapthrough’. Computers and Structures, v. 13, p. 55–62, 1981.
CRISFIELD, M. A. Non-Linear Finite Element Analysis Of Solids And Structures. London: John Wiley and Sons, 1991.
CRISFIELD, M. A. Non-Linear Finite Element Analysis Of Solids And Structures. London: John Wiley and Sons, 1997.
FAFARD, M.; MASSICOTTE, B. Geometrical interpretation of the arc-length method.
Computers and Structures, v. 46, p. 603–615, 1991.
FENG, Y.; PERI’C, D.; OWEN, D. Determination of travel direction in path-following
methods. Mathl. Comput. Modelling, v. 21(7), p. 43–59, 1995.
FENG, Y. T.; PERI’C, D.; OWEN, D. R. J. A new criterion for determination of initial
loading parameter in arc-length methods. Computers and Structures, v. 58, p. 479–485,
1996.
Referências
85
FORDE, B.; STIEME, S. F. Improved arc - length orthogonality methods for nonlinear
finite element analysis. Computers and Structures66705-709, v. 66, p. 705–709, 1987.
NETO, E. A. de S.; FENG, Y. On the determination of the path direction for arc length
methods in the presence of bifurcations and snap through. Comput. Methods Appl.
Mech. Engrg., v. 179, p. 81–89, 1999.
PIMENTA, P. M. Análise Não Linear de Treliças Espaciais. São Paulo, 1986. Boletim
Técnico.
PIMENTA, P. M. Análise Não Linear de Treliças Espaciais: teoria Exata vs. Teoria de
Segunda Ordem. São Paulo, 1996. Boletim Técnico.
RAMM, E. Strategies for tracing the nonlinear response near limit points. In: WUNDERLICH, W. (Ed.). Nonlinear finite element Analysis in Structural Mechanics. Berlin:
Spriger-Velag, 1981.
RIKS, E. The application of newton’s method to the problem of elastic stability. J. of
Applied Mechanics, v. 39, p. 1060–1066, 1972.
RIKS, E. An incremental approach to the solution of snapping and buckling problems.
Int. J. Solids Structures, v. 15, p. 529–551, 1979.
RIKS, E. Some computational aspects of the stability analysis of nonlinear structures.
Comput. Methods Appl. Mech. Engrg., v. 47, p. 219–259, 1984.
SABIR, B.; LOCK, A. C. The application of finite elements to the large-deflection geometrically nonlinear behavior of cylindrical shells. Proceedings of International Conference on Variational Mechanics, 1972.
SCHWEIZERHOF, K. H.; WRIGGERS, P. Consistent linearization for path following
methods in nonlinear f.e. analysis. Comput. Methods Appl. Mech. Engrg., v. 59, p. 261–
279, 1986.
TENG, J. G.; LOU, Y. F. Post-collapse bifurcation analysis of shells of revolution by
the accumulated arc-length method. Int. J. Numer. Meth. Engng., v. 40, p. 2369–2383,
1997.
TENG, J. G.; LUO, Y. F. A user-controlled arc-length method for convergence to predefined deformation states. Comm. In Numer. Meth. Engng., v. 14, p. 51–58, 1998.
TIAGO, M. C. Meshless methods-Extending the linear formulation and its generalization to geometrically exact structural analysis. Tese (Tese de Doutorado) — Instituto
Superior Técnico - Lisboa, 2007.
Referências
86
WAGNER, W. A path-following algorithm with quadratic predictor. Computers and
Structures, v. 39, p. 339–348, 1991.
WEMPNER, G. A. Discrete approximations related to nonlinear theories of solids. Int.
J. Solids Structures, v. 7, p. 1581–1599, 1971.
WRIGGERS, P. Nonlinear Finite Element Methods. Berlin: Springer, 2008.
WRIGGERS P., W. W.; MIEHE, C. A quadratically convergent procedure for the calculation of stability points in finite element analysis. Comput. Methods Appl. Mech.
Engrg., v. 70, p. 329–347, 1987.
WRIGTH, E. W.; GAYLORD, E. H. Analysis of unbraced multistory steel rigid frames.
Int. J. of Structural Division, v. 94, p. 1143–1163, 1968.