Embed Size (px)
Citation preview
UNIVERSIDADE FEDERAL DO PARÁ
CENTRO TECNOLÓGICO
DEPARTAMENTO DE ENGENHARIA MECÂNICA
CURSO DE MESTRADO EM ENGENHARIA MECÂNICA
A FORMULAÇÃO DE UM ELEMENTO FINITO DE BARRA PARA ANÁLISE
DINÂMICA NÃO-LINEAR GEOMÉTRICA, COM APLICAÇÃO A CABOS DE
LINHAS AÉREAS DE TRANSMISSÃO DE ENERGIA ELÉTRICA
DISSERTAÇÃO APRESENTADA COMO REQUISITO PARCIAL À
OBTENÇÃO DO GRAU DE MESTRE EM ENGENHARIA DO
DEPARTAMENTO DE ENGENHARIA MECÂNICA DO CENTRO
TECNOLÓGICO DA UNIVERSIDADE FEDERAL DO PARÁ.
Orientador: Prof. REMO MAGALHÃES DE SOUZA.
Mestrando: GANDHY YEDDO DA ROCHA ARANHA JÚNIOR.
BELÉM-PARÁ
Dezembro - 2003
GANDHY YEDDO DA ROCHA ARANHA JÚNIOR.
FORMULAÇÃO DE UM ELEMENTO FINITO DE BARRA PARA ANÁLISE
DINÂMICA NÃO-LINEAR GEOMÉTRICA, COM APLICAÇÃO A CABOS DE
LINHAS AÉREAS DE TRANSMISSÃO DE ENERGIA ELÉTRICA
Esta Dissertação foi aprovada como requisito parcial para obtenção do grau de
Mestre em Engenharia, com área de concentração em Vibrações e Acústica, no curso de
Mestrado em Engenharia Mecânica da Universidade Federal do Pará pelos professores
_____________________________________Remo Magalhães de Souza, Phd
Depto. de Construção Civil – UFPa
_____________________________________Newton Sure Soeiro, Dr. Engenharia
Depto. de Engenharia Mecânica – UFPa
_____________________________________Manoel José dos Santos Sena, Dr. Engenharia
Depto. de Engenharia Mecânica – UFPa
_____________________________________José Maria Filardo Bassalo, Dr. Física
Depto. de Física – UFPa
_____________________________________José Perilo da Rosa Neto, Msc. Engenharia
Depto. de Construção Civil – UFPa
Belém (Pa), 22 de dezembro de 2003
“Professores nunca morrem. Vivem em sua memória para sempre.Eles estavam lá quando você chegou; eles ficaram lá quando você foi embora.
Como acessórios. Às vezes lhe ensinavam alguma coisa. Mas nem sempre.E você nunca chegava a conhecê-los realmente, e nem eles a você.
Ainda assim, por algum tempo, você acreditava neles.E, se tivesse sorte, talvez um deles acreditasse em você”.
Kevin Arnold em “ANOS INCRÍVEIS”
(Dedicado aos professores Bassalo e Remo)
ii
AGRADECIMENTOS
Inicialmente quero agradecer ao professor, amigo e orientador Remo Magalhães de
Souza pelo incentivo, apoio e dedicação para concluirmos este trabalho.
Ao coordenador do curso de mestrado em Engenharia Mecânica professor Newton
Sure Soeiro, principalmente por ter acreditado em mim.
Ao meu mestre e professor José Maria Filardo Bassalo que despertou em mim o
interesse pela matemática.
Ao meu querido amigo e professor Paulo Moura Barroso por ter incentivado e
ajudado no caminho a ser seguido.
Aos meus sócios Márcia Fernanda Aranha Maranhão e Franco de Miranda Sério
Filho por terem segurado a barra de ter que trabalhar por três para que eu pudesse concluir o
mestrado.
Aos meus pais Gandhy e Maria da Graça pelo amor e por me darem o prazer de ser
seu filho.
Aos irmãos Paulo, Márcia, Marcelo e Marco por tudo que representaram,
representam e representarão em minha vida.
A minha esposa Jacirema, companheira de todas as horas, pelo amor e carinho de
todos esses anos de convívio.
Aos meus filhos Carolina, Plínio e Camila, razão de minha permanente luta.
Aos companheiros de mestrado Rivanílson e Rodrigo, que juntos conseguimos
vencer as dificuldades que naturalmente surgem.
Ao meu querido sogro Ribamar Fonseca e a minha saudosa sogra Amujacy pelo
amor e amizade em todas as horas.
Enfim aos cunhados, cunhadas, sobrinhos (em especial ao saudoso Jerônimo),
sobrinhas e a todos os meus amigos que de uma forma ou de outra são responsáveis por
aquilo que conquistamos.
iii
SUMÁRIO
AGRADECIMENTOS.......................................................................................................... ii
SUMÁRIO ........................................................................................................................... iii
LISTA DE FIGURAS...........................................................................................................vii
LISTA DE TABELAS............................................................................................................x
LISTA DE SÍMBOLOS.........................................................................................................xi
RESUMO............................................................................................................................xiv
Abstract ...............................................................................................................................xv
CAPÍTULO 1
- INTRODUÇÃO...............................................................................................................001
1.1. OBJETIVOS..................................................................................................................003
1.2. JUSTIFICATIVA ...........................................................................................................004
1.3. REVISÃO BIBLIOGRÁFICA ..........................................................................................004
a) Teorias para Barras e Formulação Co-rotacional....................................................005
b) Teorias de Cabos .........................................................................................................008
1.4. ORGANIZAÇÃO DO TEXTO .........................................................................................011
CAPÍTULO 2
- COMPORTAMENTO DOS CABOS..............................................................................012
2.1. COMPORTAMENTO ESTÁTICO ...................................................................................012
a) Relações Gerais ............................................................................................................013
iv
b) Cabo Parabólico ..........................................................................................................015
c) Cabo em Catenária ......................................................................................................016
2.2. A EQUAÇÃO DO MOVIMENTO DE UMA CORDA VIBRANTE.......................................018
2.3. CÁLCULO DO COMPRIMENTO DO CABO ....................................................................020
2.4. VIBRAÇÃO EM CABOS CONDUTORES DE LINHAS DE TRANSMISSÃO ........................021
a) Vibrações Eólicas .........................................................................................................021
b) Galope...........................................................................................................................024
c) Oscilações Devido a Esteira ........................................................................................026
CAPÍTULO 3
- FORMULAÇÃO DO ELEMENTO................................................................................027
3.1. FORMULAÇÃO DO ELEMENTO DE VIGA.....................................................................029
3.1.1. TEORIAS CLÁSSICAS DE VIGA.................................................................................029
3.1.2. EQUAÇÕES QUE GOVERNAM O PROBLEMA DE FLEXÃO ........................................030
3.1.2.1. EQUAÇÕES DE COMPATIBILIDADE.......................................................................031
3.1.2.2. EQUAÇÕES DE EQUILÍBRIO ..................................................................................035
3.1.2.3. EQUAÇÕES CONSTITUTIVAS.................................................................................037
3.2. PRINCIPIO DOS DESLOCAMENTOS VIRTUAIS.............................................................037
3.3. MATRIZ DE RIGIDEZ PARA O CASO LINEAR..............................................................041
3.3.1. SISTEMA BÁSICO .....................................................................................................043
3.3.2. SISTEMA LOCAL ......................................................................................................047
3.3.3. SISTEMA GLOBAL ....................................................................................................049
3.4. MATRIZ DE MASSA CONSISTENTE .............................................................................051
3.5. FORMULAÇÃO CO-ROTACIONAL...............................................................................052
v
3.5.1. CONFIGURAÇÃO INICIAL DA BARRA DE PÓRTICO .................................................052
3.5.2. CONFIGURAÇÃO DEFORMADA DA BARRA DE PÓRTICO .........................................053
3.5.3. TRANSFORMAÇÃO DE DESLOCAMENTOS ENTRE SISTEMAS DE COORDENADAS...055
3.5.4. TRANSFORMAÇÃO DE FORÇAS................................................................................057
3.5.5. MATRIZ DE RIGIDEZ TANGENTE NO SISTEMA GLOBAL ........................................061
3.5.6. CARGA DISTRIBUÍDA EQUIVALENTE .....................................................................063
CAPÍTULO 4
- MÉTODOS COMPUTACIONAIS.................................................................................075
4.1. MÉTODO DA RIGIDEZ DIRETA PARA A ANÁLISE LINEAR .........................................075
4.2. ANÁLISE MODAL ........................................................................................................077
4.3. ANÁLISE ESTÁTICA NÃO-LINEAR..............................................................................078
4.3.1. ALGORITMO DE NEWTON-RAPHSON ......................................................................079
4.4. ANÁLISE DINÂMICA LINEAR ......................................................................................081
4.4.1. MÉTODO DE NEWMARK PARA A ANÁLISE DINÂMICA LINEAR ..............................082
4.5. ANÁLISE DINÂMICA NÃO-LINEAR .............................................................................086
4.5.1. MÉTODO DE NEWMARK PARA A ANÁLISE DINÂMICA NÃO-LINEAR .....................086
4.6. DESCRIÇÃO DO PROGRAMA COMPUTACIONAL.........................................................088
CAPÍTULO 5
- EXEMPLOS....................................................................................................................090
5.1. TESTE DE VALIDAÇÃO PARA MATRIZ DE RIGIDEZ ...................................................090
5.2. VIGA BI-APOIADA COM CARGA DISTRIBUÍDA ..........................................................091
5.3. COMPORTAMENTO ESTÁTICO DO CABO 1.................................................................093
vi
5.4. VIGA ENGASTADA E LIVRE ........................................................................................095
5.5. COMPORTAMENTO ESTÁTICO DO CABO 2.................................................................096
5.6. PÓRTICO DE LEE ........................................................................................................098
5.7. BARRA COM UM GRAU DE LIBERDADE ......................................................................100
5.8. TORRE CARREGADA NO TOPO ...................................................................................101
5.9. VIGA BI-APOIADA ......................................................................................................103
5.10. VIGA BI-ENGASTADA ...............................................................................................105
5.11. ANÁLISE MODAL DO CABO ......................................................................................108
CAPÍTULO 6
- CONCLUSÕES E SUGESTÕES PARA TRABALHOS FUTUROS .............................112
BIBLIOGRAFIA ...............................................................................................................115
ANEXOS............................................................................................................................122
vii
LISTA DE FIGURAS
Figura 2.1 – Cabo com carregamento transversal..............................................................013
Figura 2.2 – Equilíbrio de um elemento de cabo ...............................................................014
Figura 2.3 – Corda vibrante ...............................................................................................018
Figura 2.4 – Vórtices de Karman numa esteira..................................................................022
Figura 2.5 – Direção do vento e força decorrente I ...........................................................024
Figura 2.6 – Direção do vento e força decorrente II ..........................................................024
Figura 2.7 – Vento sobre um condutor ..............................................................................025
Figura 3.1 – Exemplos de elementos finitos tri, bi e unidimensionais ..............................028
Figura 3.2 – Barra na configuração deformada com deslocamentos ( ) ( )0 0 e u x v x ........031
Figura 3.3 – Deslocamentos em um ponto qualquer de uma seção da barra .....................031
Figura 3.4 – Barra com carregamentos distribuídos e x yq q ............................................035
Figura 3.5 – (a) – Elemento com carregamento xq ...........................................................036
Figura 3.5 – (b) – Elemento com carregamento yq ..........................................................036
Figura 3.6 – Sistema básico ...............................................................................................038
Figura 3.7 – Sistema local..................................................................................................042
Figura 3.8 – Sistema local..................................................................................................047
Figura 3.9 – Equilíbrio de forças no sistema básico ..........................................................048
Figura 3.10 – Esforços e deslocamentos no sistema global...............................................050
Figura 3.11 – Barra na conf. indeformada, em relação aos sistemas local e global ..........053
Figura 3.12 – Barra na conf. deformada, em relação aos sistemas básico, local e global .054
Figura 3.13 – Barra na conf. deformada, com seus respectivos deslocamentos................056
viii
Figura 3.14 – Transformação de forças entre os sistemas global e básico ........................057
Figura 3.15 – Barra com carregamentos distribuídos e x yq q no sistema básico .............064
Figura 3.16 – Cargas distribuídas na configuração indeformada ......................................066
Figura 3.17 – Cargas distribuídas na configuração deformada..........................................066
Figura 3.18 – Equilíbrio das cargas distribuídas................................................................067
Figura 3.19 – Equilíbrio de forças no sistema básico ........................................................068
Figura 3.20 – Transformação de forças do sistema básico para o sistema local ...............069
Figura 4.1 – Convergência através do Método de Newton-Raphson ................................079
Figura 4.2 – Método da Aceleração Média........................................................................083
Figura 4.3 – Método da Aceleração Linear .......................................................................084
Figura 4.4 – Método de Newmark Não-Linear..................................................................087
Figura 5.1 – Viga bi-apoiada sujeita a carregamento uniformemente distribuído.............091
Figura 5.2 – Curva Carga x Deslocamento para diferentes níveis de discretização ..........092
Figura 5.3 – Cabo pré-tensionado com carregamento distribuído transversal...................093
Figura 5.4 – Análise de convergência da flecha do cabo...................................................094
Figura 5.5 – Flecha na extremidade do balanço de uma viga engastada e livre ................096
Figura 5.6 – Cabo sujeito a esforço de tração nos apoios..................................................096
Figura 5.7 – Redução da flecha com o aumento da força de tração ..................................097
Figura 5.8 – Pórtico de Lee................................................................................................098
Figura 5.9 – Desl. vertical do ponto de aplicação da carga para o pórtico de Lee ............099
Figura 5.10 – Trajetória completa de equilíbrio do Pórtico de Lee ...................................100
Figura 5.11 – Barra com um grau de liberdade .................................................................100
Figura 5.12 – Deslocamento horizontal de barra sujeita a uma força senoidal .................101
Figura 5.13 – Torre sujeita a carga degrau ........................................................................102
ix
Figura 5.14 – Torre e massas concentradas .......................................................................102
Figura 5.15 – Deslocamento na extremidade da torre .......................................................103
Figura 5.16 – Viga bi-apoiada sujeita a carga degrau........................................................103
Figura 5.17 – Desl. vertical no pt. médio da viga bi-apoiada sujeita a carga degrau ........104
Figura 5.18 – Viga bi-engastada sujeita a carga degrau ....................................................105
Figura 5.19 – Resposta dinâmica da viga bi-engastada .....................................................106
Figura 5.20 – Resposta dinâmica linear da viga bi-engastada ...........................................107
Figura 5.21 – Resposta dinâmica da viga bi-engastada, casos linear e não-linear ............108
Figura 5.22 – Convergência do primeiro modo de vibração .............................................110
Figura 5.23 – Convergência do segundo modo de vibração..............................................110
Figura 5.24 – Convergência do terceiro modo de vibração ...............................................111
x
LISTA DE TABELAS
Tabela 5.1 – Valores das freqüências naturais do cabo....................................................109
xi
LISTA DE SÍMBOLOS
T – Força de tração no cabo.
xT – Componente da força de tração no cabo na direção x.
0T – Componente horizontal da força de tração no cabo.
yT – Componente da força de tração no cabo na direção y.
x, y – Coordenadas do cabo.
w – Carregamento transversal contínuo por unidade de comprimento do cabo.
s – Comprimento do arco de cabo.
ρ - Densidade linear da corda.
c – Velocidade da onda na corda.
L – Comprimento total do cabo/corda e comprimento indeformado da barra.
ApL – Distância entre apoios do cabo.
n – Modo de vibração do cabo/corda.
nω - Freqüência natural do cabo/corda.
D – Diâmetro do cilindro.
U – Velocidade da corrente de ar.
f – Freqüência.
kF - Força de Karman.
E – Módulo de elasticidade.
I – Momento de inércia.
u, v – Deslocamentos longitudinais e transversais.
( )xe - Vetor de deformações.
xii
( )0 xε - Deformação axial no eixo de referência.
( )xϕ - Curvatura da seção.
e x yq q - Carregamentos distribuídos nas direções x e y.
M – Momento fletor.
V – Esforço cortante.
N – Esforço normal.
int, extW Wδ δ - Trabalho virtual externo e virtual interno, respectivamente.
p – Vetor de forças do elemento.
k – Matriz de rigidez do elemento.
d – Vetor de deslocamentos do elemento.
, , B B Bp k d – Vetor de forças, matriz de rigidez e vetor de deslocamentos do elemento no
sistema básico.
, , L L Lp k d – Vetor de forças, matriz de rigidez e vetor de deslocamentos do elemento no
sistema local.
, , G G Gp k d – Vetor de forças, matriz de rigidez e vetor de deslocamentos do elemento no
sistema global.
( )xS - Campo dos esforços reais da barra.
q – Vetor do carregamento distribuído da barra.
( )xN - Matriz das funções de forma.
( )xS - Campo dos esforços reais da barra.
( )xB - Matriz que relaciona as deformações da seção com os deslocamentos nodais.
T - Matriz de transformação entre os sistemas básico e local.
R - Matriz de rotação.
xiii
M - Matriz de massa.
1Le , 2ˆ
Le - Vetores unitários que formam a base do espaço vetorial do sistema local.
1Be , 2ˆ
Be - Vetores unitários que formam a base do espaço vetorial do sistema básico.
, B Lα α - Ângulos entre os sistemas básico e global e, local e global, respectivamente.
l – Comprimento deformado da barra.
, , P K D – Vetor de forças, matriz de rigidez e vetor de deslocamentos da estrutura no
sistema global.
Geok – Matriz de rigidez geométrica.
e B Bx yq q - Carregamentos distribuídos nas direções x e y no sistema básico.
e L Lx yq q - Carregamentos distribuídos nas direções x e y no sistema local.
Beqp - Carregamento equivalente no sistema básico.
Gqk - Matriz de rigidez do elemento no sistema global devido ao carregamento distribuído.
xiv
RESUMO
Neste estudo, apresenta-se a formulação de um elemento finito de barra para análise
estática e dinâmica geometricamente não-linear de pórticos planos e estruturas formadas por
cabos. Dá-se ênfase ao estudo de cabos condutores de linhas aéreas de transmissão de alta
tensão. O tratamento de grandes rotações das barras é realizado utilizando-se a formulação
co-rotacional, a qual é baseada na utilização de um sistema auxiliar de coordenadas que se
move com o elemento a medida que este se deforma. Uma contribuição original deste
trabalho consiste no desenvolvimento de uma nova técnica para o tratamento consistente de
cargas distribuídas através da formulação co-rotacional. Neste contexto, apresenta-se a
dedução de uma nova matriz de rigidez geométrica dependente do carregamento distribuído.
Para validação da metodologia proposta foi desenvolvido um programa na plataforma
Matlab. O programa desenvolvido permite a realização de análise geometricamente não
linear, estática ou dinâmica, de estruturas planas aporticadas e/ou formadas por cabos,
possibilitando também a determinação das freqüências naturais e modos de vibração, para
diferentes níveis de solicitação, de estruturas com comportamento altamente não-linear. O
programa é baseado no Método dos Elementos Finitos e emprega os Métodos de Newmark e
Newton-Raphson para integração da equação do movimento no tempo, para tratar os efeitos
não lineares. Vários testes numéricos são apresentados, comparando-se os resultados da
metodologia proposta com soluções analíticas simplificadas e com resultados de outros
modelos apresentados na literatura. Os resultados obtidos demonstram a eficiência e
precisão da presente formulação na análise de estruturas submetidas a grandes deformações.
Palavras-chave: análise dinâmica, não-linearidade geométrica, grandes rotações, formulação
co-rotacional, vibrações em cabos, linhas de transmissão.
xv
ABSTRACT
This study presents the formulation of a rod finite element for geometrically
nonlinear static and dynamic analysis of plane frames and cable structures. Special emphasis
is given to the study of overhead cables of transmission lines. The treatment of large
rotations is performed with the corotational formulation, which is based on the use of an
auxiliary system of coordinates that moves with the element as it deforms. One original
contribution of the present work is the development of a new technique, based on the
corotational formulation, for the consistent treatment of distributed loads. Within this
formulation, the derivation of a new geometric stiffness matrix, which depends on the
element distributed load, is presented. To validate the proposed methodology, a computer
program was developed in the Matlab platform. The program permits the geometrically non-
linear analysis of planar frames and cable structures, and also allows for the computation of
natural frequencies and mode shapes of structures with high nonlinear behavior, at different
load levels. The program is based on the Finite Element Method, and employs Newmark and
Newton-Raphson methods for time-integration of the equation of motion, considering
nonlinear effects. Several numerical tests are carried out, with the results of the proposed
methodology being compared against simplified analytical solutions and against the results
of other models presented in the literature. The obtained results demonstrate the efficiency
and precision of the presented formulation in the analysis of structures subjected to large
deformations.
Keywords: dynamic analysis, geometric non-linearity, large rotations, corotational
formulation, cable vibrations, transmission lines.
Capítulo 1
INTRODUÇÃO
A análise estrutural após o fim da II Guerra Mundial teve um avanço extraordinário
em razão da massificação dos computadores. O computador, que só era empregado por
empresas de grande porte e órgãos importantes dos governos de países mais desenvolvidos,
a partir da segunda metade da última década do século XX passou a ser até brinquedo para
crianças.
Com o perfeito casamento dos computadores com a álgebra linear e com os cálculos
vetorial e tensorial se tornou possível resolver problemas antes bastante complicados e
trabalhosos de forma eficiente e rápida, atendendo a exigência dos dias de hoje.
Uma das ferramentas bastante empregadas para a solução desses problemas é o
método dos elementos finitos (MEF). Segundo Pavanello (1997) os primeiros princípios do
método remontam a 1906, porém coube a Courant propor o MEF como é conhecido hoje.
Pavanello (1997) destaca que Courant em 1943 publicou um trabalho matemático
onde empregava o princípio da energia potencial estacionária e a interpolação polinomial
por partes sobre sub-regiões triangulares para estudar o problema de torção de Saint-Venant.
Entretanto, ao trabalho de Courant não foi dada a devida atenção até que engenheiros
tivessem desenvolvido independentemente o mesmo método.
A denominação do método foi consolidada em 1960 por Clough, conforme relata
Pavanello (1997). O valor prático do método era óbvio. Com isso novos elementos foram
desenvolvidos para aplicação em diversas áreas das ciências e tecnologia.
Um grande número de programas de elementos finitos em computadores emergiu no
final das décadas de 60 e 70, tais como SAP, ANSYS, ADINA e NASTRAN. Cada um
desses programas inclui diversos tipos de elementos com capacidade para análise estrutural
estática, dinâmica, transferência de calor, etc. A aplicação do MEF na análise estrutural é
atualmente uma atividade usual.
A análise das estruturas aporticadas é um assunto bastante complexo e vasto.
Basicamente, essa análise envolve a adoção de um modelo de viga, cujos mais conhecidos
são: o de Bernoulli-Euler, o de Timoshenko e o de Vlasov. Deve-se também estabelecer de
2
que forma a variação das ações e seus efeitos, em relação ao tempo, serão considerados. As
análises que podem ser realizadas são a análise estática, a análise quasi-estática e a análise
dinâmica.
A análise estática é aquela em que se considera que as ações (cargas, recalques, etc.)
e os respectivos efeitos (deformações, esforços internos, etc.) não variam em relação ao
tempo. Na análise quasi-estática, considera-se que as ações e seus efeitos são variáveis em
relação ao tempo, mas desprezam-se os efeitos inerciais da estrutura em análise. Finalmente,
na análise dinâmica os efeitos inerciais são considerados.
A análise estrutural pode ser efetuada dentro do regime linear ou do regime não
linear. A análise linear assume uma proporcionalidade entre causa e efeitos. Já na análise
não-linear, não se faz esta hipótese de proporcionalidade. Dois são os tipos de não-
linearidade que se pode encontrar: a não-linearidade física (ou do material) e a não-
linearidade geométrica.
A não-linearidade física está relacionada com a ausência de proporcionalidade entre
tensões e deformações ao nível do material.
Já a não-linearidade geométrica está relacionada com a ausência de
proporcionalidade entre deformações e deslocamentos (relação cinemática). Este efeito é
importantíssimo em corpos sujeitos a grandes deformações.
A solução de problemas com não-linearidade geométrica, é baseada em conceitos da
mecânica do contínuo.
A mecânica do contínuo apresenta duas formulações para a descrição do movimento
das partículas de um corpo, a descrição Lagrangeana e a descrição Euleriana.
A formulação Lagrangeana (também chamada Substantiva, Material ou Histórica) é
aquela em que o movimento é descrito acompanhando o deslocamento das partículas no
espaço através de suas trajetórias (Bassalo (1973) e Coimbra (1967)).
A formulação Euleriana (ou Espacial) é aquela em que a descrição do movimento é
feita por um observador fixo no espaço (Bassalo (1973) e Coimbra (1967)).
Quando o MEF é utilizado na solução de problemas com não-linearidade geométrica,
em geral, utilizam-se três formulações, consideradas as mais usuais. Pode-se empregar a
descrição Lagrangeana (que pode ser a atualizada e a total), a Euleriana e a co-rotacional.
3
É oportuno esclarecer que se pode fazer estudos empregando isoladamente um dos
dois tipos de não-linearidade assim como trabalhar conjuntamente com ambos.
Quando se trabalha no regime não-linear é necessário resolver sistemas de equações
não-lineares, os quais podem ser solucionados através do método de Newton-Raphson,
método de Newton modificado, e outros.
Para a análise dinâmica de uma estrutura, pode-se empregar a chamada análise
modal, a qual, conforme Soeiro (2001), descreve uma estrutura em termos de suas
características naturais, que são as freqüências naturais, os fatores de amortecimento e as
formas modais, ou seja, suas propriedades dinâmicas.
Vale a pena ressaltar que embora a análise modal seja baseada no princípio da
linearidade, pode-se fazer uma análise modal para estruturas não-lineares, em diferentes
níveis de carregamento.
Por exemplo, embora um cabo de linha de transmissão tenha um comportamento
altamente não-linear para a determinação de sua configuração deformada, pode-se fazer uma
análise modal do cabo já deformado, considerando-se o efeito da força de tração na rigidez
do cabo, e conseqüentemente nas suas freqüências naturais e modos de vibração.
Existe uma classificação para a análise dinâmica que está relacionada com o tipo de
resposta observada, a qual pode ser no domínio do tempo ou no domínio da freqüência, cuja
passagem de um para o outro é feito através da transformada de Fourier.
Para resolver as equações diferenciais que surgem quando se trabalha no domínio do
tempo se têm os métodos explícitos e os métodos implícitos. Como exemplos dos métodos
explícitos se têm os métodos de Euler, Taylor, Runge-Kutta e Adams-Bashforth (métodos
multipassos). Entre os métodos implícitos se pode citar os métodos de Newmark (aceleração
média e aceleração linear), Galerkin, Wilson e diferença central entre outros.
1.1 OBJETIVOS
Esta dissertação tem por objetivo principal o estudo de uma formulação de elementos
finitos para análise dinâmica não-linear geométrica para aplicação em estruturas constituídas
de barras e para estruturas formada por cabos. É dada ênfase aos cabos das linhas aéreas de
transmissão de energia, os quais são submetidos a grandes deformações. Outro objetivo é o
4
estudo do comportamento de barras sujeitas a cargas uniformemente distribuídas quando as
barras apresentam grandes deslocamentos (comportamento não-linear).
Como objetivos específicos se pode citar a difusão dentro da comunidade técnica e
científica local, de conhecimentos da análise dinâmica por elementos finitos, assim como a
comparação da eficiência e precisão do elemento proposto com outros elementos existentes
na literatura.
Nesta dissertação se busca fazer o estudo estático e dinâmico de estruturas
constituídas de barras e determinar as freqüências e modos naturais em diversas etapas. Para
isso serão empregados a formulação co-rotacional para a análise não-linear geométrica com
a utilização do método de Newton-Raphson para resolução dos sistemas não-lineares e o
método de Newmark para integração no tempo.
1.2 JUSTIFICATIVA
Este trabalho é justificado pela necessidade de se buscar métodos mais precisos que
garantam não somente a segurança das estruturas assim como a sua viabilidade econômica.
Também há a necessidade de atender as empresas de transmissão de energia na
região que encontram dificuldades para o controle e manutenção preventivos dos cabos das
linhas aéreas de transmissão sujeitos aos esforços do vento.
Além disso, é notória a importância do assunto dentro da comunidade científica
mundial, o que é comprovado pela quantidade de artigos publicados recentemente, tratando
tanto das teorias de barras como do estudo de cabos, nas análises estática e dinâmica.
1.3 REVISÃO BIBLIOGRÁFICA
A literatura sobre as teorias de vigas, formulação co-rotacional e teoria de cabos é
muito vasta, e por esta razão alguns poucos artigos serão citados.
Também foi encontrado um número expressivo de trabalhos abordando o
comportamento dos cabos das linhas aéreas de transmissão de energia elétrica. Dentre os
quais, serão citados principalmente os que foram publicados no Brasil, com a finalidade de
mostrar a importância do problema dentro do país.
5
a) Teorias para Barras e Formulação Co-rotacional
As teorias para barras, conforme relata Yojo (1993), foram inicialmente estudadas a
partir do século XVII. Os primeiros estudos referentes ao comportamento geometricamente
não-linear, com o emprego de material elástico linear, foram desenvolvidos já no século
XVIII por Bernoulli e principalmente por Euler.
Leckar e Sampaio (2000) relacionaram como os três modelos básicos para o estudo
de vibrações em vigas os de Euler-Bernoulli, Vlasov e Timoshenko. No modelo de Euler-
Bernoulli o cisalhamento e a inércia de rotação são desprezados, sendo feita a suposição de
que as seções planas permanecem planas e perpendiculares ao eixo longitudinal da peça,
após a deformação. No modelo de Vlasov o cisalhamento continua sendo desprezado, porém
é levada em conta a inércia de rotação. Finalmente no modelo de Timoshenko, supõe-se
também que as seções planas permaneçam planas, porém não necessariamente
perpendiculares ao eixo longitudinal da viga, pois há um giro da seção em relação a este
eixo, devido ao cisalhamento.
É oportuno esclarecer que quando as dimensões da seção transversal da peça são
pequenas quando comparadas com o seu comprimento (vigas longas e esbeltas), o modelo
normalmente empregado é o de Euler-Bernoulli. Quando as dimensões da seção transversal
não são desprezíveis quando comparadas com o comprimento da viga, o modelo mais
adequado é o de Timoshenko. Vale ainda dizer que para altas freqüências o modelo de
Euler-Bernoulli não fornece bons resultados e o modelo de Timoshenko pode ser empregado
independentemente da geometria da barra.
Na formulação de uma teoria não-linear de barra é comum se empregar o Método
dos Elementos Finitos baseado em deslocamentos com a utilização de técnicas baseadas nas
formulações Euleriana, Lagrangeana e co-rotacional, conforme já mencionado
anteriormente.
Em relação à formulação Lagrangeana ainda pode-se dizer que a sua implementação
computacional abrange as técnicas chamadas de Lagrangeana atualizada e Lagrangeana
total. A formulação Lagrangeana atualizada expressa as equações não-lineares de equilíbrio
da estrutura de forma incremental, ou seja, em função da geometria inicial de cada passo
incremental. A formulação Lagrangeana total emprega a geometria do elemento de
6
referência indeformada (fixa) ou inicial para expressar as equações não-lineares de
equilíbrio da estrutura deformada.
A formulação co-rotacional emprega uma base, em geral formada por nós das
extremidades do elemento, que gira continuamente com o elemento.
Argyris et al (1964) apresentaram um dos primeiros trabalhos aplicando a
formulação co-rotacional para vigas, conforme relata De Souza (2000). Outro dos primeiros
trabalhos foi publicado por Jennings (1968), onde o autor fez a incorporação dos efeitos da
mudança de geometria da estrutura na transformação dos deslocamentos em relação aos
sistemas de coordenadas considerados.
Belytschko e Hsieh (1973) apresentaram um elemento finito para análise transiente
de grandes deslocamentos e pequenas deformações, em problemas com não-linearidade do
material. O método empregou a formulação co-rotacional, chamado então de convected
coordinates. Na formulação do elemento foi utilizada a teoria de viga de Bernoulli-Euler, e
os resultados teóricos foram comparados com resultados experimentais e com outras
formulações numéricas, comprovando a eficiência do elemento proposto.
Oran (1973-a) empregou a formulação co-rotacional para a determinação da matriz
de rigidez tangente de pórticos planos. É oportuno informar que também nesse artigo ainda
não era empregada a denominação co-rotacional. Oran (1973-b) realizou um cálculo
semelhante para pórticos espaciais.
Oran e Kassimali (1976) estudaram o problema de grandes deformações e
estabilidade elástica de estruturas de pórticos. Foi proposto um método tipo viga-coluna
computacionalmente eficiente, que era flexível no sentido de se adaptar ao sistema de
coordenadas adotado. Para o caso da formulação dinâmica foi empregado como técnica de
integração no tempo o método de Newmark.
Mondkar e Powell (1977), apresentaram procedimentos teóricos e computacionais
para análise estática e dinâmica não-linear tanto para estruturas formadas por barras como
para treliças. A formulação empregada para o tratamento das equações do movimento foi a
Lagrangeana. O trabalho apresentou uma seqüência passo a passo para a implementação
computacional, além de diversos exemplos cujos resultados são analisados no texto. As
equações do movimento são integradas com o emprego do método de Newmark.
7
Belytschko e Glaum (1978) propuseram um elemento baseado na formulação co-
rotacional para estruturas curvas (arcos) sujeitas a cargas concentradas.
Yang e Saigal (1984) apresentaram um elemento para análise estática e dinâmica de
vigas com não-linearidade física e geométrica. O elemento simples, como foi chamado,
possuía 6 graus de liberdade e, associado com procedimentos computacionais, apresentava
solução rápida e eficiente. Diversos exemplos foram mostrados para validação do elemento
proposto.
Iura (1994) examinou os efeitos do sistema de coordenadas na precisão da
formulação co-rotacional para vigas planas de Bernoulli-Euler.
Crisfield e Shi (1994) propuseram um novo procedimento para análise dinâmica
implícita utilizando o método dos elementos finitos, obtendo soluções estáveis com
significativos movimentos de corpo rígido.
Neuenhofert e Filippou (1997) compararam modelos baseados no método da rigidez
e no método da flexibilidade para análise através do método dos elementos finitos de
pórticos com não-linearidade física.
Reddy (1997) estudou a deformação por cisalhamento através de um elemento finito
de viga formulado com o emprego das teorias de Bernoulli-Euler e Timoshenko.
Crisfield e Moita (1996) descreveram um método empregando a formulação co-
rotacional para análise de sólidos, cascas e vigas.
Crisfield et al (1997) investigaram diferentes tipos de algoritmos de integração
implícita no tempo para análise dinâmica de vigas espaciais. Os algoritmos empregados no
trabalho são baseados na formulação co-rotacional. Neste trabalho os autores propuseram
um novo método para a introdução de amortecimento numérico.
Petrangeli e Ciampi (1997) descreveram os procedimentos para análise de estruturas
compostas por barras com não-linearidade física empregando o método das forças.
Molins et al (1998) apresentaram uma formulação baseada no método da
flexibilidade com o objetivo de determinar uma nova matriz de massa para análise dinâmica
de pórticos espaciais com elementos curvos e seção transversal variável. A teoria de viga
empregada foi a de Timoshenko. A principal característica da formulação era que o
equilíbrio de forças em um ponto interior qualquer do elemento era exato.
8
Apesar do grande número de artigos encontrados pelo autor sobre teoria de vigas,
assim como sobre a formulação co-rotacional, não foram observados na revisão
bibliográficos trabalhos que tratassem de cargas distribuídas através da formulação co-
rotacional.
b) Teorias de Cabos
De acordo com Irvirne (1981), credita-se a Galileu no século XVII os primeiros
estudos da forma curva de um cabo suspenso sob o seu próprio peso. Sua contribuição se
restringiu a apontar a similaridade entre esta curva e uma parábola. Em 1691 um grupo de
geógrafos e matemáticos obteve a solução desta curva que é hoje conhecida como catenária.
Formavam este grupo os irmãos James e John Bernoulli, Leibniz e Huigens. Além disso,
outros problemas envolvendo a solução da catenária foram investigados por James
Bernoulli, inclusive uma primeira tentativa de avaliar os efeitos da deformação do cabo.
Bassalo (1996) apresenta uma revisão histórica do problema de vibração em cordas.
Ele destaca que nos anos de 1732 a 1733 Daniel Bernoulli, filho de John Bernoulli, fez
estudos sobre pequenas vibrações de uma corda de peso desprezível, pendurada em uma
extremidade e carregada com várias massas igualmente espaçadas.
Outro fato citado por Bassalo (1996) é que Daniel também trabalhou com cordas
vibrantes suspensas de espessura não uniforme
D’Alembert em 1746 publicou dois trabalhos relacionados com o problema de
vibrações de cordas, o primeiro deles sobre a forma tomada por uma corda tensa em
vibração. Nos trabalhos ele chegou a equação diferencial parcial da corda vibrante.
Em 1748, Euler apresentou um trabalho sobre a equação diferencial da corda
vibrante, onde demonstrou que se a forma inicial da corda fosse periódica, então a solução
dessa equação também seria periódica.
Bassalo (1996) também relata que Euler não informou se o somatório envolvia um
número finito ou infinito de termos, apesar de ele já considerar a idéia da superposição dos
modos vibracionais da corda.
9
Em 1753 D’Alembert comunicou à Academia de Ciências e Letras de Berlim que o
somatório tomado por Euler como solução da equação diferencial parcial da corda vibrante
envolvia um número infinito de termos.
Entretanto, as expressões encontradas esbarravam na impossibilidade da
superposição das soluções em razão do comportamento não-linear dos cabos.
Com o aparecimento dos computadores, soluções obtidas através de métodos
iterativos começaram a ser apresentadas.
Inicialmente o procedimento adotado era a divisão do cabo em diversos segmentos
finitos. Partindo de valores iniciais arbitrados para a tensão em um dos extremos do cabo,
calculava-se o resíduo para a outra extremidade, empregando as condições de equilíbrio em
um método iterativo, sendo então os valores iniciais ajustados até que o resíduo chegasse a
valores aceitáveis.
Após esses estudos e trabalhos iniciais, várias publicações apareceram em conexão
com o projeto de pontes suspensas nas quais foram apresentadas análises da resposta de
cabos parabólicos a vários tipos de carregamentos.
Outra aplicação importante que resultou em diversos estudos e publicações é o
emprego dos cabos em linhas aéreas de transmissão de energia elétrica.
Segundo Oliveira (2002) Michalos e Birnstiel propuseram em 1962 um método
numérico de tentativa e erro chamado por eles string polygon approach. Neste método, a
geometria curva do cabo é aproximada por vários segmentos retos e o carregamento
distribuído é substituído por cargas concentradas equivalentes aplicadas nas interseções
destes segmentos.
O’Brien e Francis (1964) apontaram falhas no tratamento teórico de Michalos e
Birnstiel, e apresentaram uma formulação numérica baseada nas expressões analíticas da
catenária elástica. É empregado um processo de aproximações sucessivas para a solução das
equações do equilíbrio estático dos segmentos de cabo. Nessas aproximações a flecha do
cabo tem um tratamento exato.
Com o objetivo de acelerar a convergência dos resultados obtidos por O’Brien e
Francis (1964), Jennings (1965) propôs a reformulação de algumas equações usadas por
eles.
10
O’Brien (1968) publicou um estudo sobre o comportamento dos sistemas de cabos
diante dos carregamentos impostos. Nesse trabalho foram apresentados diversos gráficos e
tabelas relacionando às cargas com a configuração assumida pelo cabo.
Antunes (1970) chama a atenção para a importância da vibração nas linhas de
transmissão, ressaltando que muitas vezes ela era esquecida. Nesse trabalho o autor
apresenta um resumo das causas das vibrações e sugere alguns meios de reduzir o problema.
Buchanan (1970) apresentou uma breve análise bidimensional empregando técnicas
de perturbação das equações diferenciais que governam o problema, com o objetivo de
avaliar a resposta não linear do cabo.
Leonard (1973) propôs o emprego de elementos curvilíneos em substituição ao uso
de elementos retos no Método dos Elementos Finitos com o objetivo de reduzir o grande
número de elementos retilíneos, porém com um elevado custo computacional do
processamento das matrizes.
Irvine (1975), com um tratamento analítico investigou a resposta de cabos
parabólicos sujeitos a vários tipos de carregamentos, porém essa teoria é restrita a cabos
inextensíveis com a razão entre a flecha e o vão inferior a 1:8.
Fried (1981) apresentou um elemento finito de cabo extensível não-linear usando
funções de interpolação polinomiais para os deslocamentos, apresentando alguns resultados
para a resposta estática e dinâmica para um cabo suspenso.
Desai et al (1995) idealizaram um elemento finito eficiente computacionalmente para
análise do fenômeno do galope em linhas de transmissão de energia elétrica. Foram
empregadas para interação com os apoios, torres com rigidez equivalente a dos condutores
nos pontos de suspensão. O programa permite tanto a não-linearidade geométrica como a
física, assim como o caso do amortecimento não-linear.
Oliveira (1996) publicou um estudo sobre os aspectos de projeto e manutenção das
linhas de transmissão de energia elétrica com relação ao problema de vibração, onde
relacionam os diversos tipos de vibrações causadas pelo vento, os tipos de condutores
empregados nas LT’s e os mecanismos empregados para o controle do fenômeno, tais como
amortecedores, espaçadores e grampos de suspensão.
11
Teixeira (1997) apresentou um modelo numérico para estudo e controle da vibração
em cabos condutores de linhas aéreas de transmissão e realizou diversas medidas para
avaliação dos neutralizadores viscoelásticos sugeridos no trabalho.
Diana et al (1998) realizaram a análise dinâmica da linha de transmissão sobre o
Lago de Maracaibo, utilizando um modelo de elementos finitos. No trabalho foram
apresentados os diferentes problemas provocados pelas correntes de vento sobre o lago.
Cella (1999) apresentou um método computacional para facilitar a determinação
matemática exata da geometria da catenária, incluindo o caso em que os apoios se
encontram em diferentes alturas. A metodologia proposta foi ilustrada para condutores de
energia elétrica.
Oliveira (2002) realizou uma análise comparativa com os elementos de cabo
existentes na literatura, apresentando o desenvolvimento de um elemento finito com
formulação isoparamétrica de cabo a partir de formulações Lagrangeana, total e atualizada.
1.4 ORGANIZAÇÃO DO TEXTO
A dissertação está dividida em 6 capítulos, que são descritos resumidamente nesta
seção.
O capítulo 2 trata do comportamento dos cabos. Considera-se tanto o aspecto
estático como o problema dinâmico, relacionado com a corda vibrante. É feita também uma
abordagem sobre a vibração dos condutores de linhas de transmissão.
O capítulo 3 apresenta a formulação do elemento proposto, onde são citadas as
teorias de vigas, as equações que governam o problema de flexão, montadas as matrizes de
massa e rigidez, bem como a formulação co-rotacional e o tratamento para a carga
distribuída equivalente.
O capítulo 4 trata dos métodos computacionais empregados no texto, ou seja, o
método da rigidez direta, o método de Newton-Raphson e o método de Newmark, tanto para
o caso linear como para o não-linear.
No capítulo 5 são apresentados os exemplos para validação do elemento proposto.
Finalmente o capítulo 6 apresenta as principais conclusões relacionadas com o texto
e as sugestões para trabalhos futuros.
Capítulo 2
Equation Section 2COMPORTAMENTO DE CABOS
O cabo flexível é um importante elemento estrutural empregado em pontes suspensas
(pênseis ou estaiadas), linhas de transmissão e suporte de teleféricos, como exemplo das
suas diversas aplicações. Conforme Meriam e Kraige (1999), para o projeto dessas
estruturas é necessário o conhecimento das relações entre os esforços, o vão, a flecha e o
comprimento dos cabos. Tais relações podem ser obtidas através da análise do cabo como
um corpo inextensível. Em geral, assume-se que o cabo apresenta resistência à flexão
desprezível, hipótese que implica que a força no cabo está sempre na sua direção
longitudinal.
Os cabos flexíveis podem suportar tanto cargas concentradas distintas como podem
suportar carregamentos que são distribuídos de forma contínua ao longo do cabo. Em
algumas situações, o peso próprio do cabo é desprezível quando comparado com as cargas
suportadas; em outros casos, o peso próprio do cabo pode ser uma carga de valor apreciável;
e em certas ocasiões é a única carga vertical atuante, não podendo, portanto, ser
negligenciada nos dois últimos casos.
O estudo dos cabos flexíveis está intimamente relacionado com o problema de
vibração de uma corda.
Conforme relata Irvine (1981) a vibração de uma corda foi estudada extensivamente
no princípio do século XVIII. A importância desse estudo é observada pelo interesse de
grandes matemáticos como D’Alembert, Euler, Daniel Bernoulli e Lagrange, que com o
estudo do problema, conseguiram expressivos avanços na teoria das equações diferenciais
parciais, já comentado na revisão bibliográfica.
2.1 COMPORTAMENTO ESTÁTICO
A seguir são feitas algumas considerações sobre o comportamento estático dos
cabos.
13
a) Relações Gerais
Seja w o carregamento transversal contínuo por unidade de comprimento aplicado ao
cabo mostrado na figura 2.1.
Tomando um pequeno elemento de comprimento dx, conforme mostra a figura 2.2, a
posição de equilíbrio requer
d dyT wds ds = −
(2.1)
onde T é a força de tração no cabo e dy/ds é o seno do ângulo de inclinação.
A componente horizontal da força de tração no cabo xT , é constante pois no caso em
estudo não há componentes horizontais nas forças externas que agem sobre o cabo. A
componente horizontal e a componente vertical da força de tração no cabo são
constante==dsdxTTx (2.2)
ydyT Tds
= (2.3)
x
h
y
x(0,0)
y
a
la
Figura 2.1 – Cabo com carregamento transversal.
Levando a equação (2.2) até à equação (2.1) e empregando a regra da cadeia se
chega na expressão
14
2
2xd y dsT wdx dx
= − (2.4)
sds ds
(T )Tds
+d dydy
ds
s w s
sdT+T
dxTds
T
d+
dsT (T
dsdsdx dx s)
Tdsdy
y
x
Figura 2.2 – Equilíbrio de um elemento de Cabo.
O comprimento de um segmento infinitesimal do cabo é dado pela expressão
2 2ds dx dy= + (2.5)
ou seja
2
1
+=
dxdy
dxds (2.6)
Substituindo a equação (2.6) na equação (2.4), tem-se
22
2 1xd y dyT wdx dx
= − +
(2.7)
Na equação (2.4), assumindo w(ds/dx) constante, a solução resulta em uma parábola.
Se apenas w for constante a solução resulta na equação da catenária. Os dois resultados são
mostrados a seguir.
15
b) Cabo Parabólico
Quando a intensidade do carregamento vertical w é constante e considerando
/ 1ds dx ≈ (hipótese válida para pequenas flechas), a equação (2.4) toma a forma2
20
d y wdx T
= − (2.8)
onde 0T é a componente horizontal da tração, constante em toda a extensão do cabo na
situação em estudo.
Resolvendo a equação (2.8), obtém-se2
02wxy
T= − (2.9)
com a aplicação das condições de contorno: 0dydx
= e 0y = , ambas para 0x = .
Para Ax l= (metade da projeção horizontal do comprimento total do cabo) e Ay h=
(flecha máxima), considerando a origem do sistema de coordenadas no meio do vão no
ponto de flecha máxima, conforme a figura 2.1, a equação (2.9) fornece o valor de 0T , que é
dado por2
0 2A
A
wlTh
= (2.10)
de modo que2
AA
xy hl
=
(2.11)
A tração T no cabo é obtida através do emprego do diagrama de corpo livre que pode
ser deduzido a partir da figura 2.2, onde se tem que
2 22 2
x ydx dyT T T T Tds ds
= + = +
(2.12)
que com o emprego das equações (2.2) e (2.3), resulta em
2 2 20T T w x= + (2.13)
Substituindo a equação (2.10) na equação (2.13), fica
16
222
2A
A
lT w xh
= +
(2.14)
onde o valor máximo da tração no cabo ocorre para Ax l= e é dado por
2
max 12
AA
A
lT wlh
= +
(2.15)
Segundo Oliveira (2002) a equação da parábola é bastante empregada no projeto de
linhas de transmissão devido à sua simplicidade matemática. Os resultados obtidos com essa
formulação são próximos aos verificados experimentalmente. Em tais estruturas as flechas
que ocorrem raramente ultrapassam a 10% do tamanho dos respectivos vãos.
c) Cabo em Catenária
Considere-se novamente o cabo da figura 2.1.
A equação (2.6) pode ser arrumada tomando a forma2 2
1dx dyds ds
+ =
(2.16)
2 2 2dx dy ds+ = (2.17)
A equação (2.4) pode ser reescrita como2
20
d y w dsdx T dx
= − (2.18)
Empregando a equação (2.17) devidamente arrumada, a expressão (2.18) toma a
forma
22
20
1d y w dydx T dx
= +
(2.19)
cuja solução é facilitada com o emprego da substituição /q dy dx= , que fornece
201
dq w dxTq
=+
(2.20)
integrando a equação (2.20), tem-se
17
( )2
0
ln 1 wq q x CT
+ + = + (2.21)
A constante C é nula, pois / 0dy dx q= = , para 0x = . Substituindo /q dy dx= na
equação (2.21), passando para a forma exponencial e simplificando, obtém-se0 0/ /
0
senh2
wx T wx Tdy e e wxdx T
−−= = (2.22)
onde a função hiperbólica é introduzida por conveniência. A equação (2.22) é integrada e o
resultado é
01
0
coshT wxy Cw T
= + (2.23)
A constante 1C é calculada a partir da condição de contorno: 0x = quando 0y = .
Com essa condição se chega a
01
TCw
= − (2.24)
e finalmente se obtém
0
0
cosh 1T wxyw T
= −
(2.25)
A equação (2.25) é a equação da curva, no caso a catenária, formada pelo cabo
suspenso sob a ação do seu peso próprio.
Das figuras 2.1 e 2.2 se observa que
0
dy wsdx T
= (2.26)
Com o emprego da equação (2.22) a equação (2.26) fica
0
0
senhT wxsw T
=
(2.27)
Com o emprego da equação (2.27) a tração em um ponto qualquer do cabo é dada
por
2 2 2 2 2 2 20 0 0
0 0
1 senh coshwx wxT w s T T TT T
= + = + =
(2.28)
ou seja
18
00
cosh wxT TT
=
(2.29)
Além disso, é possível também expressar a tração em termos de y com auxílio da
equação (2.25) que, quando substituída na equação (2.29), fornece
0T T wy= + (2.30)
A equação (2.30) mostra que o incremento da tração no cabo, a partir do valor do
ponto mais baixo (ou mais alto, dependendo do referencial), depende apenas de wy.
2.2 A EQUAÇÃO DO MOVIMENTO DE UMA CORDA VIBRANTE
O problema da corda vibrante pode ser simplificado considerando que a corda só
vibre no plano vertical e que sua amplitude de vibração seja suficientemente pequena para
que cada ponto situado na corda se mova apenas na vertical, e de maneira que a tração T na
corda não varie apreciavelmente durante a vibração.
A figura 2.3 mostra uma corda vibrante.
0
lx+x'0 x
dTd+d'
T
dx
Figura 2.3 – Corda Vibrante.
Para se obter a equação do movimento da corda, considera-se um segmento de corda
de comprimento 'x entre x e 'x x+ , onde 'x ( )'x dx= representa um incremento
infinitesimal do comprimento da corda. Seja ρ a densidade da corda por unidade de
19
comprimento, de modo que 'xρ é a massa total do segmento. A velocidade da corda em
qualquer ponto será yt
∂∂
e a inclinação yx∂∂
. A componente vertical da tração exercida da
direita para a esquerda através de qualquer ponto da corda é
senyT T θ= (2.31)
onde θ é o ângulo formado pela corda e a horizontal. Supondo-se na figura 2.5, que θ seja
pequeno, a equação (2.31) fica
sen tgyyT T T Tx
θ θ ∂= ≈ =
∂(2.32)
A força resultante para cima, dF , devido a tração, agindo sobre um segmento 'x da
corda, é a diferença das componentes verticais entre as duas extremidades do segmento, ou
seja
'y yx x x
ydF T T T dxx x+
∂ ∂ = − ≈ ∂ ∂ (2.33)
É oportuno esclarecer que quando não se limita a pequenas inclinações yx∂∂
, então
um segmento da corda poderá ter também uma componente horizontal de força resultante da
tração, e o segmento mover-se-á tanto horizontalmente como verticalmente.
Existindo uma força vertical adicional f por unidade de comprimento, agindo ao
longo da corda, a equação do movimento do segmento 'x será, de acordo com a 2ª lei de
Newton,2
2
y ydx T dx fdxt x x
ρ ∂ ∂ ∂ = + ∂ ∂ ∂ (2.34)
Considerando uma corda horizontal, em que forças horizontais atuam somente em
suas extremidades, e para pequenas amplitudes de vibração, a tração é constante e a equação
(2.34) pode ser escrita2 2
2 2
y yT ft x
ρ ∂ ∂= +
∂ ∂(2.35)
20
A força f pode ser a força gravitacional que age sobre a corda, ou ainda uma força
externa aplicada à corda para iniciar a vibração. Considerando o caso de 0f = , a equação
(2.35) pode ser reescrita como2 2
2 2 2
1 0y yx c t∂ ∂
− =∂ ∂
(2.36)
onde
Tcρ
= (2.37)
A equação (2.36) é uma equação diferencial parcial para a função ( , )y x t ,
representando matematicamente a Lei do Movimento de Newton, aplicada à corda vibrante.
Através da solução da equação (2.36) é possível se obter as freqüências naturais da
corda vibrante dadas por
nn TLπω
ρ= (2.38)
onde n = 1, 2, 3, ..., representa o modo de vibração.
2.3 CÁLCULO DO COMPRIMENTO DO CABO
De acordo com a Geometria Analítica o comprimento de uma curva qualquer é dado
por
2
1
2
1x
x
dyL dxdx
= + ∫ (2.39)
Labegalini et al (1992) resolveram a integral (2.39) e chegaram ao valor do
comprimento deformado do cabo dado pela expressão
11
2 senh2
ApLL C
C
=
(2.40)
onde ApL é a distância entre apoios e 1C é dado por
01
TCw
= (2.41)
com 0T é a componente horizontal da tração e w é o peso por metro de cabo.
21
A equação (2.40) pode ser representada pela série de Taylor3 5
11 1 1 1
1 1 12 ...2 3! 2 5! 2 ! 2
nAp Ap Ap ApL L L L
L CC C C n C
= + + + +
(2.42)
2.4 VIBRAÇÃO EM CABOS CONDUTORES DE LINHAS DE TRANSMISSÃO
Segundo Oliveira (1996) e Labegalini et al (1992) já foram identificados três tipos de
vibrações em condutores tendo como fonte a excitação do vento. Essas vibrações são
eólicas, galopping (ou galope) e de esteira.
Diversos são os efeitos das vibrações em cabos condutores de linhas de transmissão,
sendo o mais importante o prejuízo que trazem às ferragens, isoladores, estruturas e aos
cabos condutores. Os cabos mais afetados são aqueles formados por fios de alumínio, em
razão da fragilidade deste material. Os fios isolados sofrem danos mais rápidos que os
cabos. A ruptura dos condutores ocorre por fadiga do material de que são feitos os fios (aço,
cobre ou alumínio). As fraturas ocorrem nos pontos de fixação dos cabos, em razão de que
uma seção vibra e a seguinte é forçada, pela ferragem de fixação, a permanecer em posição
rígida. Neste ponto é que as ondas oscilatórias são refletidas, retornando em sentido
contrário.
a) Vibrações eólicas
As vibrações eólicas são oriundas dos redemoinhos de ar a sotavento do cabo
(causando as vibrações dos cabos no sentido vertical) no instante em que se igualam as
freqüências do vento e do condutor, que por esta razão entra em ressonância. Essas
oscilações ocorrem com ventos de velocidades na faixa de 2 a 35km/h e são verificadas em
terrenos planos ou levemente ondulados, geralmente ao amanhecer ou entardecer. Elas
geram flexões alternadas e de pequenas amplitudes nos pontos de suspensão do condutor.
Como os esforços são alternados, ocorre a ruptura do cabo por fadiga.
A vibração que mais ocorre no Brasil é essa com alta freqüência e baixa amplitude,
que depende tão somente da existência de vento lateral e cuja explicação é encontrada no
chamado trem de vórtices de Karman (ver Blessmann (1998)).
22
De uma maneira bem sucinta os vórtices de Karman ocorrem quando um fluido
escoa em torno de um obstáculo, sendo que a esteira atrás do obstáculo não é regular,
apresentando vórtices em sentidos alternados, conforme mostra a figura 2.4.
Figura 2.4 – Vórtices de Karman numa esteira.
Conforme relatam Labegalini et al (1992) os vórtices, alternadamente no sentido
horário e anti-horário, originam-se no cilindro de uma maneira perfeitamente regular e estão
associados a uma força lateral alternada. Esse fenômeno foi estudado experimentalmente e
verificou-se que há uma relação entre a freqüência f, o diâmetro do cilindro D e a velocidade
da corrente U, dada pela equação
185,0=UfD (2.43)
O número 0,185 da equação (2.43) é conhecido como número de Strouhal. A
formação alternada dos vórtices nos lados do cilindro provoca uma força harmonicamente
variável sobre o mesmo, na direção perpendicular à da corrente, cujo valor, conforme
Labegalini et al (1992),. é dado por
( )tDUCF kk ωρ sen21 2
= (2.44)
Na equação (2.44) kF é a força de Karman, kC é o coeficiente da força de Karman
(adimensional), ρ é a densidade do fluido (no caso o vento), U é a velocidade da corrente,
D é o diâmetro do cilindro, ω é a freqüência circular e t é o tempo.
Conforme Labegalini et al (1992) o mecanismo de formação de vórtices em um
cilindro estacionário é auto-excitado porque não há propriedades alternadas na corrente de
aproximação, e os vórtices são formados na freqüência natural de Strouhal. A vibração auto-
excitada geralmente é perigosa. Como regra, isso tem muitas conseqüências, quando a
23
freqüência auto-excitada de formação dos vórtices, dada pela equação (2.43), coincide com
a freqüência natural do cabo sobre o qual ela age, ocorre então uma ressonância.
Labegalini et al (1992) destaca que um parâmetro fundamental, que pode ser obtido
da equação (2.36), é a velocidade transversal da onda, dada por
trTUρ
= (2.45)
No estudo dos condutores das linhas de transmissão a equação (2.45) pode ser
refinada com a introdução de termos que levem em conta a rigidez, a flexão e o
amortecimento. No entanto, devido a natureza complexa desses parâmetros, o emprego deles
em condutores reais são apenas aproximações.
Labegalini et al (1992), por exemplo, fornecem a velocidade transversal da onda
levando em conta a rigidez a flexão, cuja expressão é2
212tr tr
EIU UT
ω ρ ′ = +
(2.46)
onde ω = 2πf é a freqüência circular (ou angular do cabo) e EI é a rigidez à flexão.
Da expressão (2.46) observa-se que a velocidade transversal da onda aumenta com a
freqüência e com a rigidez à flexão.
O amortecimento no condutor provoca uma redução na velocidade transversal da
onda. Porém é tão pequeno que pode ser desprezado em casos práticos, além do que a
aceleração da gravidade altera a velocidade de propagação da onda em razão da tensão
variar ao longo do vão.
As influências da rigidez, do amortecimento e da gravidade resultam em alterações
nas freqüências naturais de um vão do cabo. Sendo L o comprimento do vão, trU a
velocidade de propagação da onda correspondente à componente horizontal da tensão no
condutor e n = 1, 2, 3, ..., uma expressão aproximada para as freqüências naturais é dada por
Labegalini et al (1992) como sendo igual a
2n trnf UL
≅ (2.47)
onde trU é dado pela equação (2.45), quando não se considera a rigidez a flexão e pela
equação (2.46) quando a rigidez é levada em conta. Deve-se notar que a equação (2.47) recai
na expressão da freqüência para uma corda vibrante inicialmente reta (ver equação (2.38)).
24
b) Galope
O fenômeno de galope segundo Labegalini et al (1992) corresponde a uma vibração
de baixa freqüência (0,1 a 3 Hz) e de grande amplitude (5 a 300 diâmetros do condutor) que
provoca a movimentação do ponto de suspensão no sentido longitudinal dos condutores,
portanto uns contra os outros. Dependendo do comprimento do vão, a amplitude de
oscilação vertical alcança vários metros, podendo dar origem a curto circuito entre fases e
introduzir perigosos esforços nos condutores e suportes, capazes de destruir a linha. O
galope ocorre tão somente nos trechos de linhas com cadeia de suspensão. O galope é
causado por ventos laterais atuando no condutor com depósito assimétrico de gelo.
Ao soprar sobre um cilindro, o vento exerce sobre ele uma força com a mesma
direção de incidência conforme a figura 2.5, em razão da forma de sua superfície simétrica.
No caso de uma barra de seção transversal não circular, esse efeito geralmente não ocorre
dessa maneira conforme mostra a figura 2.6, havendo um ângulo entre a direção do vento e a
força. Como exemplo desse fenômeno se pode citar o vento atuando sobre a asa do avião em
que a força é aproximadamente perpendicular à direção do vento, conforme mostra
Bisplinghoff et al (1996).
Vento Força
Figura 2.5 – Direção do vento e força decorrente I.
Vento Força
Figura 2.6 – Direção do vento e força decorrente II.
25
As oscilações em linhas de transmissão podem ser estudadas como vibrações auto-
excitadas. Tem-se, então, um exemplo de vibração auto-excitada causada pelo vento que age
sobre o fio, o qual, graças ao acumulo de gelo, apresenta uma seção transversal não circular.
Esse fenômeno raramente é descrito como vibração, sendo em geral conhecido como galope,
não sendo observado, portanto, em clima quente.
A figura 2.7 mostra um cabo de uma linha de transmissão em um processo vibratório
descendente. Quando não há vento, o fio sentirá o ar vindo de baixo, devido a seu próprio
movimento para baixo. Caso haja um vento horizontal com velocidade V, o fio que está se
movendo para baixo com velocidade v sentirá um vento soprando de um ângulo arctg (v/V)
ligeiramente de baixo. Como na figura 2.7 o condutor possui uma seção circular, a força
exercida pelo vento apresenta uma pequena componente para cima. No movimento
descendente esta força gera um pequeno amortecimento. Entretanto, para o caso da seção
transversal não circular, pode ocorrer que a força exercida pelo vento tenha uma
componente para baixo e, desse modo, aparece um amortecimento negativo.
O mesmo raciocínio é empregado para o movimento ascendente da vibração.
No entanto, se o gelo acumulado no fio adquire uma seção transversal que favoreça a
uma instabilidade dinâmica, ou seja, faz com que o fio adquira uma pequena velocidade
ascendente, a ação do vento empurra-o ainda mais para cima, até que a ação elástica ou de
mola do condutor pare o movimento. Então essa força elástica move o fio para baixo, em
cujo processo o vento ajuda de novo, e pequenas vibrações transformam-se rapidamente em
grandes.
V
v Vento
Relativo
Figura 2.7 – Vento sobre um condutor.
26
O fenômeno do galope não é comum no Brasil. Porém, segundo técnicos da
Eletronorte, o fenômeno ocorre na região de fronteira do Brasil com a Venezuela.
c) Oscilações devido a esteira
A oscilação devido a esteira de acordo com Labegalini et al (1992) ocorre na região
de ar rarefeito ou de vácuo parcial formado por ventos de grande velocidade na proximidade
da linha. Quando a rarefação do ar é equivalente ao peso do condutor, normalmente os
esforços perdem a componente vertical, ocorrendo conforme a variação do vento rotações
incontroláveis dos cabos. Da mesma forma que o galope as vibrações devido a esteira
poderão provocar curtos-circuitos e introduzir esforços mecânicos consideráveis, tanto nos
condutores como sobre os suportes. Esse tipo de oscilação é peculiar aos feixes de
condutores e ocorre quando um condutor penetra na esteira gerada por um condutor
adjacente. São oscilações de grande amplitude (até 20 diâmetros do condutor) e baixa
freqüência (1 a 4 Hz), ocorrendo com ventos de velocidade na faixa de 4 a 30 m/s, conforme
relata Oliveira (1996).
Capítulo 3
Equation Section 3FORMULAÇÃO DO ELEMENTO
De acordo com Bathe (1996) o Método dos Elementos Finitos (MEF) é atualmente
um procedimento importante e freqüentemente indispensável como parte de projetos e
análises em engenharia. Programas computacionais vêm sendo desenvolvidos e empregados
em diversas áreas tecnológicas, tais como análise de tensões, transferência de calor,
escoamento de fluidos, campos eletromagnéticos, etc.
O método surgiu como uma nova possibilidade para resolver problemas da teoria da
elasticidade, superando as dificuldades dos métodos até então existentes, dentre os quais
tem-se os Métodos de Rayleigh-Ritz, Galerkin, diferenças finitas e resíduos ponderados.
O MEF comumente utilizado é baseado no método de Rayleigh-Ritz. Nele é prevista
a divisão do domínio de integração em um número finito de pequenas regiões denominadas
elementos finitos, tornando-se assim um meio contínuo em meio discreto. Portanto, o MEF
é um método aproximado de cálculo de sistemas contínuos ou reticulados em que a estrutura
é dividida em um número finito de partes (os elementos), conectadas entre si por intermédio
de pontos discretos que são chamados nós. Através da modelagem matemática a estrutura
tem o seu comportamento especificado por um número finito de parâmetros. Em particular,
nos problemas de análise estrutural, os parâmetros são os deslocamentos nodais, os quais
são as incógnitas principais do problema.
Por ser um método numérico aproximado, os resultados obtidos não são geralmente
exatos, porém os erros vão diminuindo à medida que os modelos vão sendo mais refinados.
Ou seja, o número de elementos empregados pode ser majorado até se atingir a precisão
necessária para o problema.
Os elementos finitos podem ser tri, bi e unidimensionais, conforme mostrado na
figura 3.1. Tais elementos possuem nós em suas extremidades e, de acordo com o tipo de
elemento esses nós podem ser livres para translação ou rotação em uma, duas ou três
direções constituindo os graus de liberdade da estrutura.
A opção por se trabalhar com o método dos elementos finitos se deve à sua
versatilidade, eficiência e simplicidade. Como estes elementos finitos podem ser tri, bi ou
28
unidimensionais cada um deles possui uma formulação específica. Tais formulações,
entretanto, são bastante semelhantes; a diferença básica se encontra no número de graus de
liberdade do elemento e no grau dos polinômios interpoladores. Nesta dissertação, os
elementos finitos empregados são os de barra, que constituem elementos unidimensionais.
Figura 3.1– Exemplos de elementos finitos tri, bi e unidimensionais (Pereira/2002).
Um elemento finito de barra pode ser constituído por dois nós. No caso plano, cada
um desses nós possui três graus de liberdade, dois de translação e um de rotação.
Segundo Cook et al (1989) existem três métodos para se obter a matriz característica
de um elemento: os métodos diretos, os métodos variacionais e os métodos residuais (ou
métodos aproximados).
Os métodos diretos são aqueles baseados na imposição direta do equilíbrio físico do
elemento.
Os métodos variacionais são aqueles baseados no cálculo variacional ou cálculo das
variações, onde não se trabalha com funções, mas com funcionais (um funcional é uma
entidade que depende de uma ou mais funções, ou seja, é uma função de funções).
Os métodos aproximados são aqueles onde é empregada uma função aproximadora
formada por uma combinação de funções (método de Rayleigh-Ritz) ou então utiliza
diretamente a equação diferencial (forma forte) que descreve matematicamente o problema a
ser estudado (método de Galerkin).
A matriz característica do elemento possui nome diferente de acordo com o
problema analisado. Na mecânica estrutural ela é chamada matriz de rigidez e relaciona
deslocamentos nodais e forças nodais; já em condução de calor seu nome é matriz de
condutividade, relacionando temperaturas nodais e fluxos nodais; em problemas de
29
eletromagnetismo recebe o nome de matriz de resistividade (também é chamada de matriz
de rigidez) relacionando o vetor dos potenciais nodais com o vetor das correntes nodais do
elemento e assim por diante.
3.1 FORMULAÇÃO DO ELEMENTO DE VIGA
Alves Filho (2000) estabelece que formular um elemento é basicamente determinar
sua matriz de rigidez (isto para o caso estático linear). Para a formulação da matriz de
rigidez (para o caso linear) inicialmente se deve definir os graus de liberdade do elemento.
Para assim proceder é necessário cumprir os cinco passos seguintes:
a) formulação da função de deslocamentos para o elemento finito e escolha da
função interpoladora adequada ao elemento;
b) cálculo dos coeficientes desconhecidos da função de interpolação relacionando
os deslocamentos dentro do elemento com os deslocamentos nodais;
c) cálculo das deformações internas do elemento a partir dos deslocamentos
nodais;
d) cálculo das forças internas no elemento e, como conseqüência, as tensões a partir
dos deslocamentos nodais;
e) determinação da matriz de rigidez do elemento.
Para o caso não-linear, o objetivo principal é o cálculo das forças. A matriz de
rigidez é obtida através da derivada das forças em relação aos deslocamentos.
Neste trabalho, para a determinação da matriz de rigidez (formulação do elemento de
viga), será empregado o princípio dos trabalhos virtuais, que corresponde a um método
variacional.
No MEF, em geral se utiliza o método dos deslocamentos para a obtenção da matriz
de rigidez da estrutura.
3.1.1 TEORIAS CLÁSSICAS DE VIGA
Segundo Timoshenko (1953) as teorias para barras começaram a ser estudadas
provavelmente a partir do século XVII. Em meados do século XVIII, Bernoulli e
30
principalmente Euler apresentaram trabalhos que podem ser considerados como o limiar da
teoria geometricamente não-linear para barras com o material elástico linear. É evidente que
tais teorias não continham o rigor nas conceituações como se tem atualmente. Basicamente
se pode definir a teoria não-linear, seja de barra ou não, como sendo aquela que trata de não-
linearidade geométrica, exata ou aproximada, com pequenas ou grandes deformações e nos
regimes elástico (linear ou não), plástico e viscoelástico.
Entre os mais conhecidos modelos de viga estão os de Euler-Bernoulli, Vlasov e
Timoshenko. No modelo de Euler-Bernoulli o cisalhamento e a inércia de rotação são
desprezados, sendo feita a suposição que as seções transversais planas permanecem sempre
planas e perpendiculares ao eixo longitudinal da viga após a deformação da viga. No modelo
de Vlasov o cisalhamento continua não sendo considerado, porém é levada em conta a
inércia de rotação. Por fim no modelo de Timoshenko, supõe-se também que as seções
transversais planas permanecem planas, mas não necessariamente perpendiculares ao eixo
longitudinal da viga, pois há um giro da seção em relação a essa perpendicular, devido ao
cisalhamento. É oportuno citar que quando as dimensões da seção transversal são pequenas
em comparação com o seu comprimento (vigas longas e esbeltas), o modelo mais
apropriado é o de Euler-Bernoulli. Nos casos em que as dimensões da viga não são
pequenas em comparação com o comprimento da viga (vigas curtas), situação em que o
cisalhamento deve ser considerado, o modelo de Timoshenko é o mais indicado. Para altas
freqüências o modelo de Euler-Bernoulli não fornece bons resultados e o modelo de
Timoshenko deve ser usado, independente da geometria da viga, conforme Leckar e
Sampaio (2000).
Não se faz uso das teorias de Timoshenko e Vlasov neste trabalho.
3.1.2 EQUAÇÕES QUE GOVERNAM O PROBLEMA DE FLEXÃO
As equações que governam o problema de flexão em vigas são as equações de
compatibilidade, as equações de equilíbrio e as equações constitutivas.
31
3.1.2.1 EQUAÇÕES DE COMPATIBILIDADE.
As equações de compatibilidade relacionam os deslocamentos e as deformações
atuantes em um determinado sistema.
Seja a barra representada na figura 3.2.
u0(x)
x
y������������
v0(x)
θ
Figura 3.2 – Barra na conf. deformada com deslocamentos u0(x) e v0(x) (Pereira (2002)).
Tomando uma seção transversal qualquer da barra da figura 3.2, é possível obter a
localização de um certo ponto da barra na configuração deformada conforme mostrado na
figura 3.3.
θ
y (1–cosθ)
yy cosθ
y senθ
y
Figura 3.3 – Deslocamentos em um ponto qualquer de uma seção da barra (Pereira (2002)).
De acordo com a figura 3.3 o campo de deslocamentos do ponto considerado é dado
por
32
)(sen)(),( 0 xyxuyxu θ−= (3.1)
))(cos1()(),( 0 xyxvyxv θ−−= (3.2)
nas expressões (3.1) e (3.2), )(0 xu e )(0 xv são respectivamente os deslocamentos
longitudinais e transversais de um ponto localizado no eixo de referência.
Considerando pequenos deslocamentos (ângulos pequenos) se pode considerar as
aproximações )()(sen xx θθ ≅ e 1)(cos ≅xθ . Assim as equações (3.1) e (3.2) ficam
)()(),( 0 xyxuyxu θ−= (3.3)
)(),( 0 xvyxv = (3.4)
Da teoria da elasticidade é sabido que
xyxuyxxx ∂
∂=
),(),(ε (3.5)
yyxvyxyy ∂
∂=
),(),(ε (3.6)
xyxv
yyxuyx xyxy ∂
∂+
∂∂
==),(),(2),( εγ (3.7)
Empregando a notação matricial para as equações (3.5),(3.6) e (3.7), tem-se
∂∂
∂∂
∂∂
∂∂
=
),(),(
0
0
yxvyxu
xy
y
x
xy
yy
xx
γεε
(3.8)
ou seja
( ) ( ) ( )x x x= Πe u (3.9)
onde
( )xx
yy
xy
xεεγ
=
e (3.10)
33
é o vetor de deformações,
0
( ) 0
x
xy
y x
∂ ∂ ∂
Π = ∂ ∂ ∂ ∂ ∂
(3.11)
é um operador diferencial e
( ) ( )( )
0
0
u xx
v x =
u (3.12)
é o vetor de deslocamentos.
As equações de números (3.5), (3.6) e (3.7) são chamadas de relações deformação
deslocamentos ou equações de compatibilidade, para problemas bidimensionais. Partindo da
equação (3.9) e sendo conhecido o campo de deslocamentos é possível a determinação do
campo de deformações.
De acordo com o exposto em 3.1.1 para a teoria de viga de Bernoulli-Euler, o ângulo
( )xθ da figura 3.3, para qualquer ponto do eixo de referência pode ser determinado por
)()(
)( 00 xvdx
xdvxtg ′==θ (3.13)
A equação (3.13) para ângulos pequenos pode ser aproximada como )()( xxtg θθ ≅ ,
sendo reescrita da forma
)()( 0 xvx ′=θ (3.14)
Com a substituição da aproximação obtida em (3.14) na equação (3.3), obtém-se o
campo de deslocamentos específico para a viga de Bernoulli-Euler
)()(),( 00 xvyxuyxu ′−= (3.15)
)(),( 0 xvyxv = (3.16)
As equações (3.15) e (3.16) podem ser escritas na forma matricial
34
( )'
0 0
0
( , ) ( ) ( ),( , ) ( )
u x y u x yv xx yv x y v x
− = =
u (3.17)
Com isso, é possível determinar o campo de deformações a partir das equações
dadas por (3.5), (3.6) e (3.7), e com o emprego do campo de deslocamentos que é definido
pela equação (3.17)
0 00 0
( ) ( ( ))( , ) ( ) ( )xxu x yv xu x y u x yv x
x x xε
′∂ ∂∂ ′ ′′= = − = −∂ ∂ ∂
(3.18)
0)(),( 0 =∂
∂=
∂∂
=yxv
yyxv
yyε (3.19)
0)()(
)())(()(),(),(2
00
000
=′−′=∂
∂+
∂′∂
−∂
∂=
∂∂
+∂
∂==
xvxvxxv
yxvy
yxu
xyxv
yyxu
xy
xyxy
γ
εγ(3.20)
De acordo com as equações (3.18), (3.19) e (3.20) se observa que as condições
iniciais adotadas para a viga de Bernoulli-Euler são atendidas, ou seja, não há deformação
por cisalhamento, existindo apenas deformações normais provocadas pelo esforço normal e
por momento fletor.
A equação (3.18) pode ser reescrita como
( ) ( ) ( )0,xx x y x y xε ε ϕ= − (3.21)
onde
( ) ( )0 0x u xε ′= (3.22)
é a deformação axial no eixo de referência e
( ) ( )0x v xϕ ′′= (3.23)
é a curvatura da seção.
Reescrevendo as equações (3.22) e (3.23) na forma matricial, tem-se
( )( )
( )( )
0 02
02
0
0
x u xxx v x
x
εϕ
∂ ∂=
∂ ∂
(3.24)
35
De forma mais compacta a equação (3.24) pode ser escrita como
( ) ( ) ( )xxx ue ∇= (3.25)
onde
( ) ( )( )
0 xx
xεϕ
=
e (3.26)
é o vetor de deformações generalizadas da seção,
∂∂
∂∂
=∇
2
2
0
0)(
x
xx (3.27)
é um operador diferencial.
3.1.2.2 EQUAÇÕES DE EQUILÍBRIO.
As equações de equilíbrio relacionam as forças e tensões que atuam em um
determinado sistema.
A barra mostrada na figura 3.4 está carregada nas direções x e y. A figura 3.5 mostra
um elemento infinitesimal de comprimento dx.
q y
q xd x
x
y
Figura 3.4 – Barra com carregamentos distribuídos qx e qy (Pereira/2002).
Separando as cargas distribuídas e impondo o equilíbrio nas duas direções se chega
às expressões
36
dx
qx
N N+dN
dx
qy
V V+dV
M M+dM
(a) (b)
Figura 3.5 – (a). Elemento com o carregamento qx; (b). Elemento com o carregamento qy.
0=+ xqdxdN (3.28)
0=+ yqdxdV (3.29)
0=+Vdx
dM (3.30)
Nas expressões (3.28), (3.29) e (3.30) N é o esforço normal, V é o esforço cortante e
M é o momento fletor.
Para a teoria de Euler-Bernoulli não ocorrem deformações por cisalhamento.
Conseqüentemente, o campo de deformações é expresso somente em termos das
deformações longitudinais da barra, ou seja, deformações devido a esforço normal e
momento fletor. Derivando-se a equação (3.30) uma vez em relação a x e empregando a
expressão (3.29), as equações diferenciais de equilíbrio podem ser expressas como
=+−
=+
0)(
0)(
2
2
y
x
qdx
xMd
qdx
xdN
(3.31)
37
3.1.2.3 EQUAÇÕES CONSTITUTIVAS.
As equações constitutivas relacionam as tensões e deformações atuantes em um
material, caracterizando seu comportamento como elástico (linear ou não), plástico ou
viscoelástico.
Para o material linear elástico as equações constitutivas são conhecidas como Lei de
Hooke. Para o caso em estudo, como apenas 0xxε ≠ , tem-se que
xx xxEσ ε= (3.32)
onde E é o módulo de elasticidade da seção.
3.2 PRINCÍPIO DOS DESLOCAMENTOS VIRTUAIS
Na formulação dos elementos finitos tanto através do método dos deslocamentos
como pelo método das forças é empregado o Princípio dos Trabalhos Virtuais, expresso pela
equação
intextW Wδ δ= (3.33)
De acordo com a expressão (3.33) o trabalho virtual externo realizado por um
sistema é sempre igual ao trabalho virtual interno realizado por esse mesmo sistema. Este
princípio pode ser aplicado de duas maneiras, a primeira como o princípio dos
deslocamentos virtuais, que relaciona forças reais com deslocamentos virtuais e a segunda
como o princípio das forças virtuais que neste caso relaciona deslocamentos reais com
forças virtuais.
Com o emprego do método dos deslocamentos se chega a uma matriz que transforma
deslocamentos nodais em forças nodais, tal matriz é conhecida como matriz de rigidez.
Já com o emprego do método das forças o que se faz é obter uma matriz que
transforma forças nodais em deslocamentos nodais, a qual é conhecida como matriz de
flexibilidade. No entanto é necessário determinar a matriz de rigidez do elemento que é feito
a partir da inversão da matriz de flexibilidade. Cabe ressaltar que às vezes não é possível se
realizar esta operação pois nem sempre a matriz de rigidez possui inversa. O problema é
resolvido com o emprego de um sistema de coordenadas, conhecido como sistema básico,
que elimina os modos de corpo rígido do elemento.
38
Portanto, a matriz de rigidez transforma deslocamentos em forças e na análise linear,
esta transformação é representada por
=p kd (3.34)
Quando se considera o regime não-linear a equação (3.34) fica
δ δ=p k d (3.35)
ou seja, a matriz de rigidez relaciona incrementos de forças e deslocamentos nodais.
A figura 3.6 mostra o sistema de coordenadas denominado sistema básico. Neste
sistema, um elemento de barra possui 3 (três) graus de liberdade, ou seja, um vetor de
deslocamentos nodais com três posições. Da mesma forma, o vetor das forças nodais
também possui três posições. Logo k da equação (3.34) é uma matriz quadrada de ordem
três.
d1B, p1
B
d2B, p2
B
d3B, p3
B
L
Figura 3.6 – Sistema Básico (Pereira/2002).
O princípio dos deslocamentos virtuais decorre das equações diferenciais de
equilíbrio dadas pelas equações (3.31), que colocadas na forma fraca ficam
( )( ) ( ) ( )( ) ( )0 00
0L
x yN x q u x M x q x dxδ δν ′ ′′+ + − + = ∫ (3.36)
onde L é o comprimento da viga.
39
A equação (3.36) implica na equação (3.31), porém esta equação só implica na
equação (3.36) se os deslocamentos ( )0u xδ e ( )0 xδν forem virtuais, ou seja, forem
totalmente arbitrários.
Arrumando a equação (3.36), tem-se
( ) ( ) ( ) ( ) ( ) ( )0 0 0 00 0 0 0
0L L L L
x yN x u x dx q u x dx M x x dx q x dxδ δ δν δν′ ′′+ − + =∫ ∫ ∫ ∫ (3.37)
Empregando-se a integração por partes para o primeiro e terceiro termos da equação
(3.37), obtém-se
( ) ( ) ( ) ( ) ( ) ( )0 0 00 00
LL L
N x u x dx N x u x N x u x dxδ δ δ ′′ = −∫ ∫ (3.38)
( ) ( ) ( ) ( ) ( ) ( )0 0 000 0
L LL
M x x dx M x x M x x dxδν δν δν ′′′ ′ ′= −∫ ∫ (3.39)
( ) ( ) ( ) ( ) ( ) ( )0 0 000 0
L LL
M x x dx M x x M x x dxδν δν δν′ ′ ′′′ = −∫ ∫ (3.40)
e em seguida levando a equação (3.40) até a equação (3.39), tem-se
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )0 0 0 00 00 0
L LLLM x x dx M x x M x x M x x dxδν δν δν δν′ ′′′′ ′= − +∫ ∫ (3.41)
Aplicando as equações (3.38) e (3.41) na equação (3.37), chega-se a
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
0 0 0 000 0 0
0 0 00 0 0
0
L L LL
x y
LLL
N x u x N x u x dx q u x dx q x dx
M x x M x x M x x dx
δ δ δ δν
δν δν δν
′− + + +
′ ′′′− + − =
∫ ∫ ∫
∫(3.42)
Na teoria de Berrnoulli-Euler a equação (3.30) estabelece que ( ) ( )M x V x′ = − e da
equação de compatibilidade (3.26) tem-se que ( ) ( )0 0x u xε ′= e ( ) ( )0x v xϕ ′′= , que valem
para deslocamentos virtuais.
Empregando as considerações feitas no parágrafo anterior a equação (3.42) fica
40
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
0 0 00 0 0 0
0 0 00 0 0
L L L L
x y
LLL
N x x dx M x x dx q u x dx q x dx
N x u x V x x M x x
δε δϕ δ δν
δ δν δν
+ − − =
′= + +
∫ ∫ ∫ ∫(3.43)
Utilizando os limites de integração para o segundo termo da equação (3.43), obtém-
se
( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
0 0 00 0 0 0
0 0 0 0
0 0
0 0 0 0
0 0
L L L L
x yN x x dx M x x dx q u x dx q x dx
N L u L N u V L L V
M L L M
δε δϕ δ δν
δ δ δν δν
δν δν
+ − − =
= − + − +
′ ′+ −
∫ ∫ ∫ ∫(3.44)
O sistema básico mostrado na figura 3.6 será empregado para a determinação das
condições de contorno. As condições de contorno em termos de força são
( )( )( )
1
2
3
0
B
B
B
N L p
M p
M L p
=
= −
=
(3.45)
Os deslocamentos virtuais 0uδ e 0vδ podem ser expressos em termos dos
deslocamentos nodais como
( ) ( )( ) ( )( ) ( )
0 0 1
0 0
0 2 0 3
0 00 0 0
0
B
B B
u u L dL
d L d
δ δ δδν δν
δν δ δν δ
= == =
′ ′= =
(3.46)
Levando as condições de contorno (3.45) e (3.46) até a equação (3.44), tem-se
( ) ( ) ( ) ( ) ( ) ( )0 0 00 0 0 0
3
1 1 2 2 3 31
L L L L
x y
B B B B B B B Bi i
i
N x x dx M x x dx q u x dx q x dx
p d p d p d p d
δε δϕ δ δν
δ δ δ δ=
+ − − =
= + + =
∫ ∫ ∫ ∫
∑(3.47)
Passando a equação (3.47) para a forma matricial, tem-se
41
( ) ( ) ( )( )
( )( )
0 0
00 0
1
1 2 3 2
3
L L
x y
B
B B B B
B
x u xN x M x dx q q dx
x x
dp p p d
d
δε δδϕ δν
δδδ
− =
=
∫ ∫(3.48)
ou seja
( ) ( ) ( ) ( )0 0
L LTT T B Bx x dx x dxδ δ δ= +∫ ∫S e q u p d (3.49)
onde ( ) ( )( )
N xx
M x =
S é o campo dos esforços internos reais da barra, ( )xδ e é o campo de
deformações virtuais, x
y
=
q é o vetor de carregamento distribuído externo atuante na
barra, ( ) ( )( )
0
0
u xx
v xδ
δδ =
u é o campo dos deslocamentos virtuais medido no eixo de
referência da barra, Bp é o vetor de forças nodais reais no sistema de coordenadas básico eBδd é o vetor de deslocamentos nodais virtuais no sistema de coordenadas básico.
A equação (3.49) verifica o princípio dos deslocamentos virtuais, ou seja
int extW Wδ δ= (3.50)
onde
( ) ( )int0
LTW x x dxδ δ= ∫ S e (3.51)
e
( ) ( )0
LTB B T
extW x dxδ δ δ= + ∫p d q u (3.52)
3.3 MATRIZ DE RIGIDEZ PARA O CASO LINEAR
A matriz de rigidez de um elemento é a matriz que transforma o vetor de
deslocamentos em vetor de forças
42
Quando se emprega o método das forças, a operação realizada é a transformação do
vetor de forças em vetor de deslocamentos, onde a matriz que realiza esta operação é
chamada matriz de flexibilidade. Como o objetivo é transformar deslocamentos em forças é
necessária a determinação da matriz de rigidez que é obtida através da inversão da matriz de
flexibilidade. Porém nem sempre é possível a inversão da matriz de rigidez em razão dos
modos de corpo rígido apresentados pela estrutura.
Este problema é resolvido eliminando estes modos através do emprego de um
sistema de coordenadas, conhecido como sistema básico, que elimina tais modos de corpo
rígido.
O sistema básico tomado nesta dissertação se encontra representado na figura 3.6,
onde se vê os deslocamentos e as forças atuantes no sistema.
O sistema de coordenadas local se encontra representado na figura 3.7 com
referência aos eixos cartesianos globais.
d1L, p1
L
d2L, p2
L
d3L, p3
L
d4L, p4
L
d5L, p5
L
d6L, p6
L
θ
Gy
Gx
Figura 3.7 – Sistema local.
43
3.3.1 SISTEMA BÁSICO
A matriz de rigidez no sistema básico pode ser determinada através da imposição
direta do equilíbrio, ou seja, impõe-se um deslocamento virtual unitário em cada uma das
três direções, restringindo-se os deslocamentos nas outras duas direções, para em seguida
calcular as forças correspondentes a este deslocamento em todos os graus de liberdade do
sistema. Assim procedendo para todos os três graus de liberdade e considerando a teoria de
viga de Bernoulli-Euler, chega-se a matriz de rigidez.
Uma outra maneira de se obter a matriz de rigidez do elemento no sistema básico
consiste em empregar o princípio dos trabalhos virtuais, como será mostrado a seguir.
Considerando o material linear elástico e com seção constante, tem-se que
( )( )
( )( )
000
N x xEAM x xEI
εϕ
=
(3.53)
onde A é a área e I é o momento de inércia (segundo momento de inércia) da seção.
Da equação (3.53) se obtém
( ) ( )( ) ( )
0N x EA x
M x EI x
ε
ϕ
=
=(3.54)
Levando os resultados obtidos em (3.54) até a equação (3.31)
( )( ) ( )
( )( ) ( )
0
2
2
0
0
x
y
d EA x q xdxd EI x q xdx
ε
ϕ
+ =
− + =(3.55)
Considerando barras de seção constante,
( ) ( )
( ) ( )
0
2
2
0
0
x
y
d xEA q x
dxd x
EI q xdx
ε
ϕ
+ =
− + =
(3.56)
e com emprego das expressões (3.22) e (3.23) o grupo de equações (3.56) fica
44
( ) ( )
( ) ( )
20
2
40
4
0
0
x
y
d u xEA q x
dxd v x
EI q xdx
+ =
− + =
(3.57)
Considerando o caso particular em que ( ) ( ) 0x yq x q x= = , as equações (3.57) ficam
( )
( )
20
2
40
4
0
0
d u xEA
dxd v x
EIdx
=
=
(3.58)
ou
( )
( )
20
2
40
4
0
0
d u xdx
d v xdx
=
=
(3.59)
As soluções para as equações (3.59) são do tipo
( )( )
0 0 1
2 30 0 1 2 3
u x a a x
v x b b x b x b x
= +
= + + +(3.60)
cujas condições de contorno são
( )( )( )( )( )( )
0
0 1
0
0
0 2
0 3
0 0
0 0
0
0
B
B
B
u
u L d
v
v L
v d
v L d
=
=
=
=
′ =
′ =
(3.61)
Aplicando as condições de contorno chega-se a
( )
( ) ( )
10
2 3 2 32 30 2 2
2
B
B B B BB
du x xL
d d d dv x d x x xL L
=
+ += − +
(3.62)
As equações (3.62) podem ser reescritas como
45
( )( )
10
22 3 2 30
32 2
0 0
20
B
B
B
x du x L dv x x x x xx d
L L L L
= − + − +
(3.63)
ou seja
( ) ( ) Bx x=u N d (3.64)
onde
( ) 2 3 2 3
2 2
0 0
20
xLx
x x x xxL L L L
= − + − +
N (3.65)
é a matriz das funções de forma no sistema básico.
Considerando as equações (3.25) e (3.64), as deformações generalizadas da seção
podem ser obtidas como
( ) ( ) ( ) ( ) ( ) ( )B Bx x x x x x= ∇ ∇ =e u = N d B d (3.66)
onde
( ) ( ) ( ) 2 2 3 2 3
2 2 2
0 0 0
20 0
xx Lx x x
x x x xxx L L L L
∂ ∂= ∇ =
∂ − + − + ∂
B N (3.67)
( )2 2
1 0 0
4 6 2 60
Lxx x
L L L L
= − + − +
B (3.68)
é a matriz que relaciona as deformações da seção com os deslocamentos nodais.
Através do PTV tem-se que int extW Wδ δ= , onde (desconsiderando o carregamento
distribuído)
( ) ( )int0
LTW x x dxδ δ= ∫ S e (3.69)
e
( )TB BextWδ δ= p d (3.70)
Portanto, igualando as equações (3.69) e (3.70) tem-se
46
( ) ( ) ( )0
LTB B T x x dxδ δ= ∫p d S e (3.71)
Com o emprego da equação (3.66), obtém-se
( ) ( ) ( )0
LTB B T Bx x dxδ δ= ∫p d S B d (3.72)
que arrumada fica
( ) ( ) ( )0
0L
TB B T x x dxδ
− =
∫d p S B (3.73)
como Bδd é arbitrário
( ) ( ) ( )0
LTB T x x dx= ∫p S B (3.74)
ou
( ) ( )0
LB T x x dx= ∫p B S (3.75)
A equação (3.53) pode ser reescrita como
( ) ( )( )
( )( ) ( ) ( )00
0 S
N x xEAx x x
M x xEIεϕ
= = =
S k e (3.76)
onde
( )0
0S
EAx
EI
=
k (3.77)
é a matriz de rigidez da seção.
Substituindo a equação (3.76) na equação (3.75) obtém-se
( ) ( ) ( ) ( ) ( ) ( )0 0
L LB T T B
S Sx x x dx x x x dx= =∫ ∫p B k e B k B d (3.78)
ou seja,B B B=p k d (3.79)
onde
( ) ( ) ( )0
LB T
Sx x x dx= ∫k B k B (3.80)
47
é a matriz de rigidez do elemento no sistema básico.
Resolvendo a integral (3.80) chega-se a
0 0
4 20
2 40
B
EAL
EI EIL LEI EIL L
=
k (3.81)
3.3.2 SISTEMA LOCAL
A matriz de rigidez no sistema local pode ser obtida através da transformação de
coordenadas feita com o emprego da matriz T. A matriz T é chamada de matriz de
transformação do sistema básico para o sistema local. A seguir é mostrada sua dedução para
o caso linear, por questões didáticas. O caso não-linear será estudado mais adiante com o
emprego da formulação co-rotacional.
A figura 3.8 mostra as forças atuantes no elemento em relação ao sistema de
coordenadas locais.
d1L, p1
L
d2L, p2
L
d3L, p3
L
d4L, p4
L
d5L, p5
L
d6L, p6
L
Figura 3.8 – Sistema local (Pereira/2002).
Para a transformação de coordenadas se deve inicialmente fazer o equilíbrio de
forças no sistema básico, conforme mostra a figura 3.9.
48
p1B
p2B
p3B
p1B p2
B+ p3B
L
p2B+ p3
B
L
Figura 3.9 – Equilíbrio de forças no sistema básico (Pereira/2002).
É possível agora fazer a associação dos esforços nodais e reações mostrados na
figura 3.9 com os esforços nodais do sistema local, de acordo com a figura 3.7. Os
resultados são
( )
( )
2 31 1 2 3 2
2 34 1 5 6 3
B BL B L L B
B BL B L L B
p pp p p p p
Lp p
p p p p pL
+= − = =
+= = − =
(3.82)
As equações (3.82) podem ser colocadas na forma matricial como
1
21
32
43
5
6
1 0 01 10
0 1 01 0 0
1 10
0 0 1
L
LB
LB
LB
L
L
pp L L pp
pp
pp
L Lp
− = − − 1442443
TT
(3.83)
ou seja
=L T BP T P (3.84)
Empregando o princípio dos trabalhos virtuais, que implica que o trabalho virtual é o
mesmo em ambos os sistemas, tem-se que
49T TL L B Bδ δ=p d p d (3.85)
Levando a equação (3.84) até a equação (3.85), tem-seB Lδ δ=d T d (3.86)
De posse da equação (3.84) e empregando as equações (3.79) e (3.86), chega-se aL T B T B B T B L= = =p T p T k d T k Td (3.87)
ou sejaL L L=p k d (3.88)
ondeL T B=k T k T (3.89)
é a matriz de rigidez do elemento no sistema local. A título de ilustração, a matriz de rigidez
no sistema local para uma barra no regime linear, é
3 2 3 2
2 2
3 2 3 2
2 2
0 0 0 0
12 6 12 60 0
6 4 6 20 0
0 0 0 0
12 6 12 60 0
6 2 6 40 0
L
EA EAL L
EI EI EI EIL L L LEI EI EI EIL L L L
EA EAL L
EI EI EI EIL L L LEI EI EI EIL L L L
− − − = −
− − − −
k (3.90)
3.3.3 SISTEMA GLOBAL
Após a determinação da matriz de rigidez no sistema local é possível se encontrar
facilmente a rigidez no sistema global de coordenadas.
Na figura 3.7 são mostrados os sistemas local e global com θ representando o
ângulo entre esses sistemas de coordenadas.
50
d1G, p1
G
d2G, p2
G
d3G, p3
G
d4G, p4
G
d5G, p5
G
d6G, p6
G
θ
Figura 3-10 – Esforços e deslocamentos no sistema global.
Com auxílio das figuras 3.8 e 3.10, e empregando relações trigonométricas
elementares se obtém as expressões abaixo que relacionam os esforços no sistema global
com os esforços no sistema local
1 1 2 4 4 5
2 1 2 5 4 5
3 3 6 6
cos sin cos sinsin cos sin cos
G L L G L L
G L L G L L
G L G L
p p p p p pp p p p p p
p p p p
θ θ θ θθ θ θ θ
= − = −= + = +
= =(3.91)
O grupo de equações (3.91) pode ser escrito na forma matricial ficando
−
−
=
L
L
L
L
L
L
G
G
G
G
G
G
pppppp
pppppp
6
5
4
3
2
1
6
5
4
3
2
1
1000000cossen0000sencos0000001000000cossen0000sencos
θθθθ
θθθθ
(3.92)
A equação (3.92) pode ser escrita de maneira mais compacta como
G T L=p R p (3.93)
onde se tem a matriz TR que relaciona os esforços no sistema global com os esforços no
sistema local e é chamada na análise matricial estrutural de matriz de transformação.
51
É fácil demonstrar (verificar) que a matriz de rotação é ortogonal, ou seja, sua
inversa é igual a sua transposta e ainda que os deslocamentos nos sistemas global e local
estão também relacionados pela matriz de rotação, ou seja
G T L=d R d (3.94)
ou aindaL G=d Rd (3.95)
Para o sistema local é sabido que os esforços e os deslocamentos nodais estão
relacionados pela matriz de rigidez no sistema local, dada pela equação
L L L=p k d (3.96)
que substituída na equação (3.93) e empregando a equação (3.95) fica
G T L L T L G= =p R k d R k Rd (3.97)
de onde se tira queG T L=k R k R (3.98)
é matriz de rigidez no elemento no sistema de coordenadas globais.
3.4 MATRIZ DE MASSA CONSISTENTE
A matriz de massa consistente é aquela obtida utilizando-se as mesmas funções de
interpolação para deslocamentos e acelerações ao longo da viga.
O método dos elementos finitos resolve por aproximação determinado problema
variacional, para tanto ele necessita empregar funções interpoladoras ou aproximadoras, as
quais são definidas para cada elemento.
O cálculo da matriz de massa consistente para elemento de barra empregando a
teoria de viga de Bernoulli-Euler pode ser encontrada em diversos livros. Como referência
bibliográfica, recomenda-se o trabalho de Pereira (2002) por apresentar de maneira bem
detalhada a dedução da matriz de massa consistente.
Como a matriz de massa consistente para um elemento de viga é bastante conhecida,
limita-se aqui a fornecer o resultado, qual seja
52
−−−−
−−
=
22
22
4220313022156013540
000021
3130422013540221560
002100
420
LLLLLL
NNLLLLLL
NN
mLM (3.99)
onde N = 140, e, m é a massa por unidade de comprimento da barra.
3.5 FORMULAÇÃO CO-ROTACIONAL
O desenvolvimento a seguir é baseado na apresentação de De Souza (2000). É
utilizado um sistema de coordenadas sem modos de corpo rígido, denominado sistema
básico, já mostrado na figura 3.6. A transformação entre os sistemas básico e local, ou entre
os sistemas básico e global é feita através da formulação co-rotacional.
Utilizando esse procedimento a formulação do elemento no sistema básico é
completamente independente da transformação, isto é, no sistema básico o elemento pode
ser formulado como linear e a não-linearidade geométrica é introduzida na transformação.
Para esse caso, a formulação pode representar arbitrariamente grandes movimentos de corpo
rígido, porém com pequenas deformações ao longo do elemento. Ainda assim, se os
membros estruturais forem divididos em pequenos elementos, a formulação co-rotacional
pode ser utilizada para a solução de problemas com deformações finitas.
3.5.1 CONFIGURAÇÃO INICIAL DA BARRA DE PÓRTICO
A figura 3.11 mostra um elemento de viga na configuração indeformada, cujas
coordenadas no sistema local são dadas por ( )LL yx , . Sejam 1Le e 2ˆ
Le os vetores unitários que
formam a base do espaço vetorial do sistema local. O vetor 1Le pode ser escrito como
1cosˆsen
LL IJ
L Lαα
= =
Xe (3.100)
com
53
J IIJ J I
J I
X XY Y
= = −
X X - X (3.101)
que representa a diferença entre as coordenadas globais dos nós J e I , e, L é o comprimento
inicial, ou indeformado, da barra, dado por
( )1/ 2TIJ IJ IJL = =X X X (3.102)
O subscrito IJ é empregado para representar a diferença entre duas grandezas
relacionadas com os nós I e J.
O outro vetor da base ortonormal do sistema local é representado por 2ˆLe e é obtido
através da expressão
2senˆ
cos
LL
L
αα
−=
e (3.103)
x G
Gy
XII
XJ
L JIX
Lx1e
x
J
L
2e L
y L
Figura 3.11 – Barra na configuração indeformada, em relação aos sistemas local e global.
3.5.2 CONFIGURAÇÃO DEFORMADA DA BARRA DE PÓRTICO
O elemento de barra mostrado na figura 3.11 possui dois nós I e J e seis graus de
liberdade tanto em relação ao sistema local como em relação ao sistema global. Cada nó
possui dois graus de liberdade de translação e um grau de liberdade de rotação.
54
O vetor dos deslocamentos em relação ao sistema global é então representado por
1
2
3
4
5
6
GI
GI I
GI IG
GJ J
GJ J
GJ
UdVd
dUdVd
d
γ γ
γγ
= = =
U
dU
(3.104)
onde IU e JU são os vetores com as duas translações, nas direções X e Y, e Iγ , Jγ são as
rotações em torno do eixo Z, com os subscritos representando os respectivos nós.
O elemento na configuração deformada é representado na figura 3.12, onde IU e JU
são os vetores que representam os deslocamentos dos nós I e J medidos em relação ao
sistema local. Sejam 1Be e 2ˆ
Be os vetores unitários que formam a base do espaço vetorial do
sistema básico. O vetor 1Be pode ser escrito como
1cosˆsen
BB IJ IJ
B lαα
+= =
X Ue (3.105)
IJX
Lx
J
I
IX JX
Ly
By
Bx
IU
JU
IJ IJ+X U
1ˆBe
2ˆBe
Lα
Bα
Gx
Gy
2ˆLe 1ˆLe
Figura 3.12 – Barra na configuração deformada, em relação aos sistemas básico, local e
global.
onde
55
IJ J I=U U -U (3.106)
representa a diferença entre as coordenadas locais dos deslocamentos dos nós J e I , e, l é o
comprimento deformado da barra, dado por
( ) ( )1/ 2T
IJ IJ IJ IJ IJ IJl = + = + + X U X U X U (3.107)
O outro vetor unitário do sistema básico é
2senˆ
cos
BB
B
αα
−=
e (3.108)
3.5.3 TRANSFORMAÇÃO DE DESLOCAMENTOS ENTRE SISTEMAS DE
COORDENADAS
A figura 3.13 mostra os deslocamentos do elemento de barra no sistema básico. De
acordo com esta figura, o deslocamento axial da barra é dado pela diferença entre o
comprimento deformado e o comprimento inicial do elemento, ou seja
Llud B −==1 (3.109)
As rotações mostradas na figura 3.13 em coordenadas no sistema básico são dadas
por
( ) BLII
Bd ααγθ −+==2 (3.110)
( ) BLJJ
Bd ααγθ −+==3 (3.111)
O ângulo Bα pode ser obtido através da equação
++
=IJIJ
IJIJB
UXVY
arctanα (3.112)
onde o numerador da expressão (3.112) corresponde à soma das ordenadas dos vetores IJX
e IJU , e o denominador corresponde a soma das abscissas dos mesmos vetores.
De Souza (2000) chama a atenção para o uso das expressões (3.110), (3.111) e
(3.112) em razão de que estas expressões somente são válidas para o caso em que °< 90Bα
devido a presença da função computacional arctan.
Esta restrição pode ser evitada utilizando-se o procedimento a seguir.
56
As equações (3.110) e (3.111) podem ser modificadas para
( )( ) ( ) IB
IBB
IBL
II βαβααβααγθ cossensencossensensen −=−=−+= (3.113)
( )( ) ( )cos cos cos cos cos sen senL B B B BI I I I Iθ γ α α β α α β α β= + − = − = + (3.114)
com LII αγβ += .
O valor de Iθ é então dado por
cos sen sen cosarctancos cos sen sen
B BI I
I B BI I
α β α βθα β α β
−= +
(3.115)
�������
���������������
Bx
ByIγ
Lα
Ju
L
Bα
l
������������
��������
������������
Iθ
BαJθ
Lx
Figura 3.13 – Barra na configuração deformada, com seus respectivos deslocamentos.
Procedendo da mesma maneira e considerando LJJ αγβ += , o valor de Jθ é dado
por
cos sen sen cosarctancos cos sen sen
B BJ J
J B BJ J
α β α βθα β α β
−= +
(3.116)
Notar que mesmo nos casos onde Bα é grande, Iθ e Jθ são em geral menor que 90°,
o que permite a utilização das equações (3.115) e (3.116) sem nenhum problema.
57
3.5.4 TRANSFORMAÇÃO DE FORÇAS
As forças no sistema global estão relacionadas com as forças no sistema básico
através das seguintes expressões (exatas), que são deduzidas a partir da figura 3.14BBBG Qpp αα sencos11 −−= (3.117)BBBG Qpp αα cossen12 +−= (3.118)
BG pp 23 = (3.119)
BBBG Qpp αα sencos14 += (3.120)BBBG Qpp αα cossen15 −= (3.121)
BG pp 36 = (3.122)
onde
lpp
lMM
QBB
JI 32 +=
+= (3.123)
é a força cortante nas extremidades do elemento, considerando o comprimento deformado
no sistema básico.
1Bp
Nó I Nó J
x xxy x
xxQ 1p
Gx p 4GQ Gx
p 2G
G
p 5G
p 1B
B
y G
BBy
Gy
BB
B
B
B
Figura 3-14 – Transformação de forças entre os sistemas global e básico.
Substituindo a equação (3.123) nas equações de número (3.117) a (3.122) e
reescrevendo na forma matricial, tem-se
58G T B=p T p (3.124)
onde
−−
−−
−−
=
1cossen0cossen
0cossen1cossen0sencos0sencos
llll
llllBBBB
BBBB
BBBB
αααα
αααααααα
T (3.125)
é a matriz de transformação de forças.
A relação tangencial entre os deslocamentos nos sistemas básico e global pode ser
calculada através da derivada dos deslocamentos no sistema básico Bd , dado pelas
equações (3.109), (3.110) e (3.111), em relação aos deslocamentos no sistema global Gd ,
tal que
BB G
Gδ δ∂=∂dd dd
(3.126)
onde Bδd e Gδd são as variações infinitesimais dos deslocamentos nos sistemas básico e
global respectivamente.
É importante observar que trabalhar com a transformação entre os sistemas básico e
global ou básico e local é muito semelhante, pois a diferença é que quando se trabalha com o
sistema básico e o sistema local é necessário fazer mais uma transformação de coordenadas,
do local para o global com o emprego da matriz de transformação dada por
cos sen 0 0 0 0sen cos 0 0 0 00 0 1 0 0 00 0 0 cos sen 00 0 0 sen cos 00 0 0 0 0 1
θ θθ θ
θ θθ θ
−
= −
R (3.127)
Para o cálculo da matriz B G∂ ∂d / d é necessário que as variações das quantidades que
definem as componentes dos deslocamentos no sistema básico, como lδ , δα e B1eδ , sejam
determinadas. A variação do comprimento deformado é obtida através da equação (ver
equação (3.107))
59
( ) ( ) ( )
( )
1/ 21 22
1
T TIJ IJ IJ IJ IJ IJ IJ
TIJ IJ IJ
l
l
δ δ
δ
− = + + + =
= +
X U X U X U U
X U U(3.128)
1( )B T GIJlδ δ δ= =e U r d (3.129)
onde, de acordo com as equações (3.104) e (3.106), tem-se
1 1ˆ ˆ[ ( ) 0 ( ) 0] cos sen 0 cos sen 0B T B T B B B Bα α α α = − = − − r e e (3.130)
A variação do vetor unitário ˆB1e é obtida a partir da equação (3.105)
1 12
1 1 1ˆ ˆ( ) ( )B BIJ IJ IJ IJl l
l l lδ δ δ δ δ= − + = −e U X U U e (3.131)
A variação do ângulo de rotação Bα pode ser determinada com o emprego da
equação (3.105) e da equação (3.108)
1 2cos senˆ ˆsen cos
B BB B B B
B B
α αδ δ δα δα
α α −
= = =
e e (3.132)
Em razão de ˆB1e e ˆB2e serem vetores unitários, multiplicando-se ambos os lados da
equação (3.132) por ˆ( )TB2e , chega-se a
2 1ˆ ˆ( )B B T Bδα δ= e e (3.133)
Substituindo a equação (3.131) na equação (3.133) e considerando que os vetores ˆB1e
e ˆB2e são ortogonais, tem-se
( )2 1 21 1 1ˆ ˆ ˆ( ) ( )B B T B B T G
IJ IJll l l
δα δ δ δ δ= − = =e U e e U s d (3.134)
com
2 2ˆ ˆ[ ( ) 0 ( ) 0] sen cos 0 sen cos 0B T B T B B B Bα α α α = − = − − s e e (3.135)
De posse das variações lδ e Bδα , a variação dos deslocamentos do sistema básico
podem ser calculados com o emprego das equações (3.109), (3.110) e (3.111), obtendo-se
60
01
1
B B G GI I I
BJ J J
u l
l
l
δ δδ δθ δγ δα δγ δ δ
δθ δγ δα δγ
= = − = + − = − −
r
d s d T d
s
(3.136)
Considerando as equações (3.104), (3.130) e (3.135), chega-se a
0 0 0 0 0 010 0 1 0 0 0
0 0 0 0 0 1
B
G
l
l
∂ = = + − ∂
−
rdT sd
s(3.137)
cos sen 0 cos sen 0sen cos sen cos1 0
sen cos sen cos0 1
B B B B
B B B B
B B B B
l l l l
l l l l
α α α α
α α α α
α α α α
− − = − − − −
T (3.138)
Notar que a equação (3.138) é a mesma equação (3.125).
Esta correspondência é uma extensão do Princípio da Contragradiência da análise
estrutural para o caso da não-linearidade geométrica. Este princípio pode ser derivado do
princípio do trabalho virtual, considerando que o trabalho efetuado pelas forças Bp junto
com os deslocamentos virtuais δ Bd no sistema básico, é igual ao trabalho realizado pelas
forças Gp junto com os deslocamentos virtuais δ Gd no sistema global. Empregando a
equação (3.136)
T T TB B G G G T Bδ δ δ= =d p d p d T p (3.139)
ou seja, como Gδd são arbitrários
G T B=p T p (3.140)
61
3.5.5 MATRIZ DE RIGIDEZ TANGENTE NO SISTEMA GLOBAL
A matriz de rigidez do elemento no sistema básico relaciona os incrementos de
deslocamentos com os incrementos de forças, ou sejaB B Bδ δ=p k d (3.141)
De acordo com De Souza (2000) a matriz de rigidez tangente Gk no sistema global
é obtida através da linearização da equação (3.124) e com o emprego das equações (3.136)
e (3.141), tal que
( )( )
G T B T B T B T B B T B
T B G G GGeo
δ δ δ δ δ δ
δ δ
= = + = + =
= + =
p T p T p T p T k d T pT k T k d k d
(3.142)
ondeG T B
Geo= +k T k T k (3.143)
é a matriz de rigidez no sistema global de coordenadas. O segundo termo do membro à
direita da equação (3-124) é a matriz de rigidez geométrica, dada porT
BGeo G
∂=∂Tk : pd
(3.144)
com o símbolo “:” representando uma contração, tal que3
1
T B T B Gr r Geo
rpδ δ δ
=
= =∑T p t k d (3.145)
com rt representando as linhas da matriz de transformação T.
Portanto, a matriz de rigidez geométrica é determinada empregando as variações rtδ
para cada uma das três linhas da matriz T e multiplicando o resultado obtido pela
correspondente força Brp no sistema básico.
Outra maneira de se obter a equação (3.143) é empregando a notação indicial
(considerando a regra do somatório implícito de Einstein). Para isto se parte da relação
( )BG Bki kG Bi ki k
ij k kiG G G Gj j j j
T pp T pk p Td d d d
∂∂ ∂ ∂= = = +∂ ∂ ∂ ∂
(3.146)
62
como
ij
BkiGeo kG
j
Tk pd∂
=∂
(3.147)
a equação (3.146) fica
ij
B BG k lij Geo ki B G
l j
p dk k Td d∂ ∂
= +∂ ∂
(3.148)
comoB
B kkl B
l
pkd∂
=∂
(3.149)
Bl
lj Gj
dTd∂
=∂
(3.150)
tem-se finalmente
ij
Gij Geo ki kl ljk k T k T= + (3.151)
onde as expressões (3.143) e (3.151) são idênticas.
Partindo da equação (3.136) e empregando as equações (3.130), (3.132), (3.134) e
(3.135), a seguinte relação pode ser obtida
11T T T B T G
lδ δ δα δ= = =t r s s s d (3.152)
Os segundo e terceiro termos T3
T2 tt δδ = são encontrados com o emprego da segunda
(ou terceira) linha da matriz T na equação (3.136) e o uso das equações (3.129) e (3.134)
( )
2 2 2
2
1 1 1 1 1
1
T T T T T B T G
T T G
ll l l l l
l
δ δ δ δ δα δ
δ
= − = − + = + =
= +
t s s s r s r d
r s s r d(3.153)
Substituindo as equações (3.152) e (3.153) em (3.145), tem-se
2 312
B BBT B T T T Gp pp
l lδ δ
+ = + +
T p s s r s s r d (3.154)
de onde se tira que a matriz de rigidez geométrica é
2 312
B BBT T T
Geop pp
l l + = + +
k s s r s s r (3.155)
podendo ser representada de forma mais compacta como
633
1
BGeo r r
rp
=
=∑k G (3.156)
onde2 2
2 2
2 2
2 2
0 00 0
0 0 0 0 0 01 10 00 0
0 0 0 0 0 0
T
s cs s cscs c cs c
l l s cs s cscs c cs c
− − − −
= = − − − −
1G s s (3.157)
e
2 3 2
1 T T
l = = + G G r s s r (3.158)
2 2 2 2
2 2 2 2
2 3 2 2 2 2 2
2 2 2 2
2 0 2 02 0 2 0
0 0 0 0 0 012 0 2 0
2 0 2 00 0 0 0 0 0
cs c s cs c sc s cs c s cs
l cs c s cs c sc s cs c s cs
− − − + − − + −
= = − + − −
− + − −
G G (3.159)
com Bc αcos= e Bs αsen= .
De Souza (2000) informa que de acordo com Kassimali a expressão (3.143) para a
matriz da rigidez geométrica foi obtida pela primeira vez por Oran em 1973 e aperfeiçoada
por Crisfield em 1990. Oran obteve a expressão da matriz de rigidez geométrica de forma
elegante e consistente para o caso bidimensional. Embora o trabalho de Powell de 1969
apresente uma formulação semelhante.
3.5.6 CARGA DISTRIBUÍDA EQUIVALENTE
O método da rigidez considera de forma direta apenas cargas atuando nos nós da
estrutura. As cargas distribuídas que atuam na estrutura devem ser consideradas como um
carregamento equivalente. Por essa razão, segundo Gere e Weaver (1981), as cargas são
divididas em dois tipos: cargas nodais e cargas que atuam nos membros.
64
A figura 3.15 mostra o elemento de barra com carregamentos distribuídos nas duas
direções x e y, com referência ao sistema básico de coordenadas.
qy
qx p1
B
p3Bp2
BB
B
Figura 3.15 – Barra com carregamentos distribuídos Bxq e B
yq no sistema básico.
Para se resolver o problema das cargas atuantes fora dos nós de uma estrutura será
necessário o emprego das funções de forma. Para a barra da figura 3.15 as funções de forma
( )N x são dadas pela equação (3.62), repetida aqui por conveniência
2 3 2 3
2 2
0 0( )
20
xlx
x x x xxl l l l
= − + − +
N (3.160)
O vetor de cargas distribuídas Bq é dado por
BxBy
qq =
Bq (3.161)
Empregando o princípio do trabalho virtual se tem que o trabalho realizado pelo
carregamento externo é igual ao trabalho realizado pelo carregamento equivalente, ou seja
eq extW Wδ δ= (3.162)
( ) ( )0
LT TB B B Beq dxδ δ= ∫d p u q (3.163)
como ( )B B Bxδ δ=u N d e ( ) ( ) ( )( )T T TB B B xδ δ=u d N , a equação (3.163) fica
( ) ( ) ( )0
LT TB B B T Beq x dxδ δ= ∫d p d N q (3.164)
como Bδd é arbitrário, tem-se
65
( )0
LB T Beq x dx= ∫p N q (3.165)
Aplicando as equações (3.160) e (3.161) na equação (3.165) e resolvendo a integral,
tem-se o vetor de cargas equivalentes no sistema básico, dado por
2
2
2
12
12
Bx
ByB
eq
By
q L
q L
q L
= −
p (3.166)
A equação (3.166) é um resultado bastante conhecido na análise estrutural linear.
Porém, é oportuno ressaltar que o comprimento L corresponde ao comprimento do elemento
de barra na configuração indeformada.
Para o caso linear, a transformação do carregamento equivalente do sistema básico
para o sistema local é realizado com o emprego da matriz de transformação dada pela
equação (3.83) e com a adição das reações de apoio do elemento no sistema básico com
sinal trocado, ou seja
1
21
32
43
5
6
1 0 01 10
20 1 0 01 0 0 0
1 102
0 0 1 0
Bx
L By
LB
LB
LB
L By
L
q Lp q Lp L L pp
pp
pp q L
L Lp
−− − = +
− − −
(3.167)
O segundo termo do membro à direita da equação (3.167) é representado por LRp e
corresponde aos esforços referentes às reações de apoio do sistema básico medidas no
sistema local.
Para o caso não linear o vetor das cargas equivalentes no sistema básico é facilmente
determinado com o auxílio da figura 3.16, que mostra a barra na configuração indeformada e
da figura 3.17, que mostra a barra na configuração deformada.
66
x L
x G
x L
y L
x
Lq
yG qy
L
Figura 3.16 – Cargas distribuídas na configuração indeformada.
Considera-se nesta formulação que a direção e a intensidade da resultante das cargas
distribuídas não se altera à medida que o elemento se deforma. Como exemplo deste tipo de
carregamento tem-se o peso próprio da estrutura.
By
xxq
yqB
BB
yG
x G
Bx
Figura 3.17 – Cargas distribuídas na configuração deformada.
A figura 3.18 mostra as componentes das resultantes das cargas distribuídas em
relação aos sistemas básico e global.
67
x
yq L
LxBx
yL
Lq
B
xq
LL
y LqB
G
x
x G
Figura 3.18 – Equilíbrio das cargas distribuídas.
Por imposição direta do equilíbrio, mostrado na figura 3.18, tem-se
( )cos senB L Lx x yq q qα α= + (3.168)
( )sen cosB L Ly x yq q qα α= − + (3.169)
ondeB Lα α α= − (3.170)
representa a diferença entre as inclinações da barra nos sistemas básico e local.
A equação (3.166) permanece válida para o caso não-linear, ou seja
2
2
2
12
12
Bx
ByB
eq
By
q L
q L
q L
= −
p (3.171)
Finalmente com o emprego da equação (3.171) se tem o vetor de cargas
equivalentes, para o caso não linear, dado por
( )
( )
( )
2
2
cos sen2
sen cos12
sen cos12
L Lx y
B L Leq x y
L Lx y
Lq q
Lq q
Lq q
α α
α α
α α
+ = − +
−
p (3.172)
Para a obtenção das reações de apoio do sistema básico no sistema local, para o caso
não linear, considera-se inicialmente a figura 3.19 que mostra o elemento de barra no
68
sistema básico com o carregamento distribuído nas duas direções e as respectivas reações de
apoio.
p1B
p2B
p3B
p1B p2
B+ p3B
l
p2B+ p3
B
l
B
Bqx
qy
qy
L
B
2
qy
2
B+
qxB LL+
I
J
Figura 3.19 – Equilíbrio de forças no sistema básico.
Para simplificar, na figura 3.19, o elemento é mostrado na configuração
indeformada, porém está sendo considerado o efeito não linear geométrico, sendo o
comprimento l do elemento função dos deslocamentos sofridos pelos nós.
O vetor de forças no sistema local é determinado com a imposição direta do
equilíbrio de forças nos nós I e J, conforme mostra a figura 3.20. Nesta figura se tem que
2 3
2
B BB
I yp p LQ q
l+
= − e 2 3
2
B BB
J yp p LQ q
l+
= + são as forças constantes no nó I e nó J,
respectivamente. Os resultados obtidos são
1
1 2
3
sen sencos cos sen2
B
L B B Bx y
B
pLp p q L q
l lp
α αα α α = − − − + − +
(3.173)
1
2 2
3
cos cossen sen cos2
B
L B B Bx y
B
pLp p q L q
l lp
α αα α α = − + − −
(3.174)
69
( )1
3 2
3
0 1 0 0
B
L B
B
pp p
p
= +
(3.175)
1
4 2
3
sen sencos sen2
B
L B By
B
pLp p q
l lp
α αα α = +
(3.176)
1
5 2
3
cos cossen cos2
B
L B By
B
pLp p q
l lp
α αα α = − − + −
(3.177)
( )1
6 2
3
0 0 1 0
B
L B
B
pp p
p
= +
(3.178)
p
Nó I Nó J
L+ qx
xy
I xxQ p
J pQ
pp p
yy
y
B
1B B
x
1
L
4L
2L
L
5L
1B
L
BB
L
B
x x L
x x
Figura 3.20 – Transformação de forças do sistema básico para o sistema local.
É conveniente reorganizar as equações de número (3.173) a (3.178) em uma única
equação matricial de tal forma que facilite tanto a notação como os cálculos que serão feitos
a seguir. Essa equação na forma matricial fica
70
1
21
32
43
5
6
sen sencos cos sen2
cos cossen sen cos2
0 1 0 0sen sencos
cos cossen
0 0 1
B Bx y
LB B
L x yB
LB
LB
yL
L
Lq L ql l
p Lq L qp l l pp
pp qp
l lpp
l l
α αα α α
α αα α α
α αα
α αα
− − − − + − − − = + − −
sen2
cos20
B
By
L
Lq
α
α
−
(3.179)
Na equação (3.179) a segunda parcela do membro à direita, será denotada por LRp ,
que corresponde às reações de apoio do sistema básico medidas em coordenadas locais.
Neste instante cabe uma digressão. A transformação de forças do sistema básico para
o sistema local está sendo feita com o emprego de um ângulo α que é obtido através da
diferença de dois outros ângulos (ver equação (3.170)). Por outro lado a transformação do
sistema local para o sistema global consiste simplesmente do emprego da matriz de rotação.
Com o objetivo de tornar o programa mais eficiente propõe-se a transformação direta do
sistema básico para o sistema global, cuja única mudança no que foi feito até agora é a
substituição do ângulo α pelo ângulo Bα . No mais a forma das equações é a mesma, apenas
com a troca das forças obtidas, inicialmente para o sistema local, e que agora são para o
sistema global.
A equação (3.179), agora para o sistema global, fica
1
21
32
43
5
6
sen sencos cos sen2
cos cossen sen
0 1 0sen sencos
cos cossen
0 0 1
B BB B B B B
x y
G B BB B
G xB
GB
G B BBB
G
B BGB
Lq L ql lp
q Lp l l pp
pp
pl lp
pl l
α αα α α
α αα
α αα
α αα
− − − − +
− −
= + − −
cos2
0
sen2
cos20
B B By
B By
B By
Lq
Lq
Lq
α α
α
α
− −
(3.180)
71
ou, G T B GR=p T p + p , onde
cos sen2
sen cos2
0
sen2
cos20
B B B Bx y
B B B Bx y
GR
B By
B By
Lq L q
Lq L q
Lq
Lq
α α
α α
α
α
− + − − = −
p (3.181)
A solução para a matriz de rigidez no sistema básico com a atuação das cargas
distribuídas pode ser obtida a partir das equaçõesB B B B
eq= −p k d p (3.182)
B Gδ δ=d T d (3.183)G T B G
R= +p T p p (3.184)
Tomando inicialmente a equação (3.184), a variação do vetor Gp é dada por
( )G T B G T B T B GR Rδ δ δ δ δ= + = + +p T p p T p T p p (3.185)
porém T B GGeoδ δ=T p k d , e com o emprego da equação (3.182), a equação (3.185) fica
( )G G T B B B GGeo eq Rδ δ δ δ= + − +p k d T k d p p (3.186)
G G T B B T B GGeo eq Rδ δ δ δ δ= + − +p k d T k d T p p (3.187)
empregando a equação (3.183), a equação (3.187) toma a formaG G T B G T B G
Geo eq Rδ δ δ δ δ= + − +p k d T k T d T p p (3.188)
( )G T B G T B GGeo eq Rδ δ δ δ= + − +p k T k T d T p p (3.189)
0G G G T B G
eq Rδ δ δ δ= − +p k d T p p (3.190)
onde
0G T B
Geo= +k k T k T (3.191)
é a matriz de rigidez em coordenadas globais, sem levar em conta os efeitos do
carregamento distribuído. Os termos referentes ao carregamento distribuído são dados por
72G G T B Gq eq Rδ δ δ= − +k d T p p (3.192)
ou, em função do vetor de reações de apoio no sistema local, porG G T B T Lq eq Rδ δ δ= − +k d T p R p (3.193)
Escrevendo dessa forma, o termo da matriz de rigidez referente ao carregamento
distribuído pode ser encontrado facilmente. Considerando inicialmente o vetor de cargas
equivalentes no sistema básico dado pela equação (3.172), a variação deste vetor pode ser
obtida diferenciando ambos os lados dessa expressão, cujo resultado é
( )( )( )
6 sen cos
cos sen12
cos sen
L Lx y
B L Leq x y
L Lx y
q qLL q q
L q q
α α
δ α α δα
α α
− + = − −
+
p (3.194)
ou
6sen 6coscos sen
12cos sen
LxB
eq Ly
qL L Lq
L L
α αδ α α δα
α α
− = − −
p (3.195)
É importante observar que como B Lα α α= − , tem-seB L Bδα δα δα δα= − = (3.196)
pois 0Lδα = , já que Lα é constante.
Com os resultados obtidos na equação (3.195), pode-se escreverB Leq aδ δα=p q (3.197)
onde
6sen 6coscos sen
12cos sen
L L LL L
α αα αα α
− = − −
a (3.198)
e Lq é o vetor de cargas distribuídas no sistema local, conforme definido na eq. (3.161).
Fazendo a substituição da carga distribuída no sistema básico pela carga distribuída
no sistema local, dadas pelas equações (3.168) e (3.169), a variação do vetor das reações de
apoio do sistema básico, medido em coordenadas globais, é obtida diferenciando ambos os
lados da expressão (3.181) modificada para o sistema local, ou seja
73
( ) ( )
( ) ( )
( ) ( )
( ) ( )
2 2
2 2
2 2
2 2
2cos sen cos sen2 2
cos -sen 2sen cos2 2
0
2sen cos cos sen2 2
cos sen 2sen2 cos2 2
L Lx y
L Lx y
LR
L Lx y
L Lx y
L Lq q
L Lq q
L Lq q
L Lq q
α α δα α α δα
α α δα α α δα
δα α δα α α δα
α α δα α α
− − + − − − + −
= − + − − + − −
p
0
δα
(3.199)
A equação (3.199) pode ser arrumada ficando
sen 2 cos 2cos 2 sen 2
0 0sen 2 cos 22
cos 2 sen 20 0
LxL
R Ly
qLq
α αα α
δ δαα αα α
− − −
= −
p (3.200)
Na equação (3.200) a matriz multiplicada por / 2L pode ser representada por b e o
vetor das cargas distribuídas no sistema local por Lq , com essas notações a equação ficaL LRδ δα=p bq (3.201)
Para o sistema global se temG T L T LR Rδ δ δα=p R p = R bq (3.202)
Levando as equações (3.197) e (3.201) até a equação (3.193), tem-seG G T L T Lq δ δα δα= − +k d T aq R bq (3.203)
( )G G T L T Lq δ δα= − +k d T aq R bq (3.204)
Empregando a expressão (3.134) que fornece a variação do ângulo entre o sistema
básico e o sistema global, a equação (3.204) fica
( )1G G T T L Gq lδ δ= −k d T a + R b q s d (3.205)
Como Gδd é arbitrário, tem-se
74
( )1G T T Lq l= +k -T a R b q s (3.206)
A equação (3.193) levada até a equação (3.190) fica
0G G G G G
qδ δ δ= +p k d k d (3.207)
que pode ser escrita como
( )0G G G G
qδ δ= +p k k d (3.208)
onde o termo entre parêntesis corresponde a matriz de rigidez do elemento procurada.
O termo dado pela equação (3.206) corresponde à contribuição do carregamento
distribuído na matriz de rigidez global do elemento.
Portanto, a matriz de rigidez é composta das parcelas referentes ao material, à não-
linearidade geométrica e ao carregamento distribuído ao longo do elemento.
Capítulo 4
Equation Section 4MÉTODOS COMPUTACIONAIS
Neste capítulo apresentam-se os métodos computacionais empregados nesta
dissertação.
4.1 MÉTODO DA RIGIDEZ DIRETA PARA ANÁLISE LINEAR
O método da rigidez direta possui este nome em razão de que a matriz de rigidez da
estrutura é montada a partir da contribuição direta das matrizes de rigidez dos elementos.
O outro método empregado na análise matricial das estruturas é o método da
flexibilidade (ou método das forças) que basicamente difere do método da rigidez pelo fato
de empregar forças como incógnitas enquanto que o da rigidez emprega os deslocamentos
dos nós como incógnitas.
Uma estrutura pode ser satisfatoriamente representada por uma associação de
elementos discretos (com suas propriedades físicas e geométricas), interligados de modo a
compor a estrutura como um todo.
O vetor das forças nodais P está relacionado com o vetor dos deslocamentos nodais
D, para o caso linear, através da expressão
=P KD (4.1)onde K representa a matriz de rigidez da estrutura em estudo.
Seja Gd o vetor dos deslocamentos do elemento, que estão relacionados com os
deslocamentos nodais através da equaçãoG
el=d A D (4.2)em que elA é uma matriz retangular, denominada matriz de incidência cinemática, composta
de zeros exceto nas posições correspondente aos graus de liberdade do elemento.
O carregamento externo referente ao vetor de deslocamentos da estrutura é extP .
Relacionando o carregamento externo com o correspondente deslocamento virtual de D
com a utilização do princípio do trabalho virtual, tem-se que o trabalho virtual externo é
dado por
76
T extextWδ δ= D P (4.3)
A energia de deformação virtual (trabalho virtual interno) pode ser determinada
através de
1
elnT
intel
Wδ δ=
=∑ d p (4.4)
De acordo com o princípio do deslocamento virtual as equações (4.3) e (4.4) são
iguais, ou seja
ext intW Wδ δ= (4.5)
1
elnT ext T
elδ δ
=
=∑D P d p (4.6)
Expandindo o segundo membro da equação (4.6) e com o emprego da equação (4.2),
obtém-se
( )
( ) ( )1 1
1 1
el el
el el
n nTT ext T
elel el
n nT T T T T T
el el elel el
δ δ δ
δ δ δ
= =
= =
= = =
= = +
∑ ∑
∑ ∑
D P d p A D p
D A p D A D A p(4.7)
como a matriz elA é constante, tem-se que 0elδ =A , portanto a equação (4.7) fica
1
elnT ext T T
elel
δ δ=
=∑D P D A p (4.8)
como Tδ D é arbitrário
1
elnext T
elel=
= ∑P A p (4.9)
A equação (4.9) é uma equação de equilíbrio que relaciona as cargas nodais extP
com as forças internas p .
A equação (4.1) para o elemento fica
=p kd (4.10)Levando as equações (4.2) e (4.10)até a expressão (4.9), tem-se
1 1
el eln next T T
el el elel el= =
= = =∑ ∑P A kd A kA D KD (4.11)
onde
77
1
elnTel el
el==∑K A kA (4.12)
é a matriz de rigidez de toda a estrutura.
Pré-multiplicando ambos os lados da equação (4.11) pela inversa da matriz de
rigidez da estrutura 1−K , tem-se1 1ext− −= =K P K KD D (4.13)
ou seja1 ext−=D K P (4.14)
Onde a equação (4.14) possibilita a obtenção dos deslocamentos da estrutura a partir
da inversa da matriz de rigidez.
A equação (4.14) também ilustra a filosofia do método dos deslocamentos que é
empregar os deslocamentos nodais como incógnitas principais do problema.
4.2 ANÁLISE MODAL
A equação do movimento para um sistema linear com vários graus de liberdade sem
amortecimento é dada por
( )ext t+ =&&MD KD P (4.15)
e quando se tem o amortecimento é
( )ext t+ + =&& &MD CD KD P (4.16)
A análise modal trata da resposta dinâmica do sistema quando solicitado pela carga
externa. ( )ext tP .
A determinação das freqüências naturais é um problema de autovalor e pode ser
resolvido da seguinte maneira:
a) A solução da equação (4.15) para o caso de vibrações livres ( )( )ext t =P 0 é dada
por cos( )tω θ= −D φ , onde φ é o vetor correspondente aos modos de vibração
da estrutura.
b) Substituindo a solução obtida anteriormente na equação (4.15) ajustada para
vibrações livres, tem-se
78
( )2ω− = 0K M φ (4.17)
c) A equação (4.17) admite a solução 0φ = . A fim de que esta equação admita
outras soluções é necessário impor que
( )2det 0ω− =K M (4.18)
que nada mais é do que um problema generalizado de autovalor. De posse dos
valores de 2ω , retornar-se à equação (4.17) e se determinam os valores de φ .
4.3 ANÁLISE ESTÁTICA NÃO-LINEAR
A análise estática não-linear neste trabalho é feita com o emprego do método de
Newton-Raphson. Este método foi publicado por Isaac Newton em 1687 e sistematizado por
Joseph Raphson em 1860, razão pela qual o método é conhecido com o nome dos dois
matemáticos, conforme relata Cunha (2000).
O método de Newton-Raphson é simplesmente um processo numérico para a
determinação de raízes (ou solução) de um sistema não-linear de equações.
Aplicado ao problema em estudo, o método pode ser entendido como um processo
iterativo em que a posição dos nós é ajustada a cada iteração, partindo-se de uma
configuração inicial arbitrária, até a convergência, quando se atinge a configuração final de
equilíbrio, dentro da precisão previamente estabelecida.
O método de Newton-Raphson combina duas idéias básicas bastante empregadas nos
procedimentos numéricos: linearização e iteração.
Na linearização, procura-se, em uma certa vizinhança, substituir um problema
complicado por uma aproximação linear que, na maioria das vezes, é bem mais simples de
ser resolvida. Essa aproximação linear pode ser obtida através do emprego dos primeiros
termos de uma série de Taylor para a função do problema. No caso do método de Newton-
Raphson, a linearização consiste na troca da curva por sua tangente.
Já um processo iterativo consiste na repetição sistemática de um procedimento até
que seja atingida a aproximação desejada.
79
4.3.1 ALGORITMO DE NEWTON-RAPHSON
De acordo com Bathe (1996), o principal requisito no equilíbrio dos elementos é
encontrar a solução que satisfaça a seguinte expressão
( ) 0ext −P P D = (4.19)
onde extP é o vetor das cargas externas e P é o vetor das forças internas, ambas aplicadas
nos nós. O vetor das forças internas P é função dos deslocamentos nodais, porém é uma
função não-linear que torna necessário o emprego de um método iterativo para resolver a
equação (4.19).
Desprezando momentaneamente para a análise estática os dois primeiros termos da
equação (4.16), tomando ( )tP constante ou variando “quasi” estaticamente e tomando-se
uma pequena variação na equação (4.1),tem-se
δ δ=P K D (4.20)
A figura 4.1 mostra a convergência obtida através do método de Newton-Raphson.
P
D
)Dn+1(Pi
P(i D)
P )(D1i
D(P )i 2
ni D(P )
Pext
Dn+1nD2D1D
Figura 4.1 – Convergência através do Método de Newton-Raphson.
A diferença entre o carregamento externo e o vetor de forças internas é dada por( ) ( )j jext∆ −P = P P (4.21)
80
A rigidez é definida como a derivada das forças internas em relação ao
deslocamento, ou seja
( )∂∂
P D= K
D(4.22)
A primeira iteração é obtida partindo da origem, ou seja
0 =D 0 (4.23)A força interna, a diferença entre o carregamento externo e a força interna, a rigidez
e o deslocamento são então dados respectivamente por
( ) ( )( )
( ) ( )
0
0
00
-10 0 0 0 0 0
1 0 0
ext ext
=
∆ − =
∂ ∂=
∂ ∂∆ ∆ ⇒ ∆ ∆
+ ∆
0 0
0
0
P D = P
P = P P P
P D PK =
D DK D = P D = K P
D = D D
(4.24)
A segunda iteração é feita a partir do cálculo da força interna no ponto obtido através
da última equação do grupo de equações (4.24) e em seguida são realizadas as mesmas
operações da iteração anterior, tem-se então
( )( )
1 1
11
-11 1 1 1 1
2 1 1
exti∆ −
∂∂
∆ ∆ ⇒ ∆ ∆+ ∆
P = P P D
P DK =
DK D = P D = K P
D = D D
(4.25)
O procedimento é repetido até que se tenha
( )( )
( )
-1
1
1 1 erro
extn i n
nn
n n n n n n
n n next
n i n
+
+ +
∆ −
∂∂
∆ ∆ ⇒ ∆ ∆+ ∆
∆ − ≤
P = P P D
P DK =
DK D = P D = K P
D = D DP = P P D
(4.26)
onde o erro é igual a precisão requerida ou desejada para o problema em estudo.
Há um modelo alternativo para o método de Newton-Raphson apresentado. Trata-se
do chamado método de Newton-Raphson modificado cuja diferença do método mostrado
81
está no fato de ser empregada sempre a mesma rigidez, eliminado, portanto, um passo no
procedimento descrito anteriormente.
Não há necessariamente uma vantagem geral em se empregar o método tradicional
ou o modificado. Em alguns casos, a convergência pode ser mais rápida para o primeiro
método e em outros casos pode ser mais rápida para o segundo, cabendo ao usuário fazer a
verificação de qual método é o mais adequado para o problema em questão. Neste trabalho
se optou pelo método tradicional em razão da grande não-linearidade geométrica das
estruturas analisadas, barras e cabos.
4.4 ANÁLISE DINÂMICA LINEAR
A integração no tempo da equação do movimento pode ser realizada através de
diversos métodos, dentre os quais citam-se:
a) Método da Diferença Central;
b) Método θ de Wilson;
c) Método de Houbolt;
d) Métodos de Newmark:
d1) Método da Aceleração Linear;
d2) Método da Aceleração Média.
Os métodos empregados nesta dissertação são os métodos de Newmark. Os métodos
de Newmark são métodos fundamentados na escolha da variação da aceleração. Estes
métodos foram desenvolvidos por N. M. Newmark em 1959 baseados nas seguintes
equações:
( ) ( )1 11i i i it tγ γ+ += + − ∆ + ∆ & & && &&D D D D (4.27)
( ) ( )2 21 1
1( )2i i i i it t tβ β+ +
= + ∆ + − ∆ + ∆ & && &&D D D D D (4.28)
Nas expressões (4.27) e (4.28) β e γ definem como a aceleração varia em um certo
intervalo de tempo, a estabilidade e a precisão características do método. Segundo Chopra
(1995) é usual o emprego de 21
=γ , e de acordo com as condições de precisão e estabilidade
para o método, β deve estar compreendido no intervalo 41
61
≤≤ β . Através das equações
82
(4.27) e (4.28) combinadas com a equação do movimento dada pela expressão (4.16), sendo
conhecidos ,i i&D D e i
&&D no tempo i, pode-se determinar 1 1,i i+ +&D D e 1i+
&&D no tempo i+1, ou
seja
1 1 1 1ext
i i i i+ + + ++ + =&& &MD CD P P (4.29)
Os métodos de Newmark são estáveis quando se tem que βγπ 2
12
1−
<∆
nTt . Para
o caso da aceleração linear 21
=γ e 61
=β , resultando que 551,0<∆
nTt , já para o método da
aceleração média 21
=γ e 41
=β , portanto ∞<∆
nTt , que torna o método incondicionalmente
estável.
4.4.1 MÉTODO DE NEWMARK PARA ANÁISE DINÂMICA LINEAR
Neste parágrafo são estudados os métodos de Newmark da aceleração média e o da
aceleração linear.
a) Método da Aceleração Média
A figura 4.2 ilustra esquematicamente a idéia básica do método da aceleração média.
De acordo com a figura se tem que o intervalo de tempo empregado para a avaliação da
convergência do método é dado por ii ttt −=∆ +1 , ou seja este é o valor do passo de tempo
utilizado para se fazer a integração no tempo.
A aceleração em um tempo arbitrário τ , (com 1i it tτ +< < ) é dada por
( ) ( )112 i iτ += +&& && &&D D D (4.30)
integrando a equação (4.30), obtém-se
( ) ( )12i i iττ += + +& & && &&D D D D (4.31)
( )1 12i i i it
+ +∆
= + +& & && &&D D D D (4.32)
83
i
i i+1t t t
i+1
D
D
D
Figura 4.2 –Método da Aceleração Média.
A equação (4.31) fornece a velocidade em um instante arbitrárioτ e a equação
(4.32) fornece a velocidade no final do intervalo de tempo t∆ . Integrando a equação (4.31)
chega-se ao deslocamento em um instante τ e no final do intervalo t∆ , dados por
( ) ( )2
14i i i iττ τ += + + +& && &&D D D D D (4.33)
( ) ( )2
1 14i i i i i
tt+ +
∆= + ∆ + +& && &&D D D D D (4.34)
b) Método da Aceleração Linear
A idéia básica do método da aceleração linear é ilustrada esquematicamente na
figura 4.3. De acordo com a figura se tem que o intervalo de tempo empregado para a
avaliação da convergência do método é também dado por ii ttt −=∆ +1 , ou seja este é o valor
do passo de tempo utilizado para se fazer a integração no tempo.
A aceleração em um tempo arbitrário τ é dada por
( ) ( )1i i itττ += +∆
&& && && &&D D D - D (4.35)
integrando a equação (4.35), chega-se a
( ) ( )2
12i i i itττ τ += +∆
& & && && &&D D + D D - D (4.36)
( )1 12i i i it
+ +∆
=& & && &&D D + D + D (4.37)
84
que fornecem a velocidade em um instante arbitrárioτ e em intervalo de tempo t∆ . É
importante observar que a equação (4.37) é idêntica à equação (4.32). Integrando as
equações (4.36) e (4.37), chega-se ao deslocamento em um instante τ e no final do intervalo
t∆ , respectivamente dados por
( ) ( )2 3
1 12 6i i i i itτ ττ τ + += + + + −
∆& && && &&D D D D D D (4.38)
( )21 1
1 16 3i i i i it t+ +
= + ∆ + ∆ +
& && &&D D D D D (4.39)
tt t i+1
i+1
i
DD
Di
Figura 4.3 –Método da Aceleração Linear.
c) Sistematização do Método
Chopra (1995) sugere os seguintes passos para o emprego dos métodos de Newmark:
I) Escolha do método a ser empregado, o da aceleração linear
==
61,
21 βγ
ou o da aceleração média
==
41,
21
βγ ;
II) Escolha, de acordo com as condições iniciais do problema em questão, das
aproximações iniciais: 0 0&P ,D e 0D , e determinação das matrizes de rigidez
( )K , massa ( )M e amortecimento ( )C ;
III) Escolha do passo de intervalo de tempo t∆ , que atenda aos critérios de
estabilidade do método, e de precisão da análise;
85
IV) Cálculo da aceleração inicial através da equação
( )10 0 0 0
−= − −&& &D M P CD KD (4.40)
V) Cálculo da rigidez equivalente ( )K com o emprego da expressão
( )21ˆ
t tγβ β∆ ∆
K = K + C + M (4.41)
VI) Cálculo dos termos a e b dados pelas equações
1t
γβ β
= +∆
a M C (4.42)
1 12 2
t γβ β
= + ∆ −
b M C (4.43)
VII) Sendo conhecido o vetor de forças externas e ext ext exti i+1 i∆ = −P P P , pode-se
definir a variação da força externa comoext
i iˆ
i i∆ = ∆ + +& &&P P aD bD (4.44)satisfazendo a equação
ˆ ˆi i∆ = ∆P K D (4.45)
VIII) Para cada passo de tempo se deve calcular-1ˆ ˆ
i i∆ = ∆D K P (4.46)
12i i i it
tγ γ γβ β β
∆ = ∆ − + ∆ − ∆ & & &&D D D D (4.47)
( )21 1 1
2i i i itt β ββ∆ = ∆ − −
∆∆&& & &&D D D D (4.48)
IX) De posse dos resultados das equações de números (4.46), (4.47) e (4.48),
calcula-se
1i i i+ = + ∆D D D (4.49)
1i i i+ = + ∆& & &D D D (4.50)
1i i i+ = + ∆&& && &&D D D (4.51)X) Calcula-se o novo valor de ext ext ext
i i+1 i∆ = −P P P , agora para o tempo i+1, e
com os resultados das expressões (4.49), (4.50) e (4.51), inicia-se
novamente o procedimento a partir da etapa VII até ser atingido o tempo
final.
86
4.5 ANÁLISE DINÂMICA NÃO-LINEAR
Neste trabalho, a análise dinâmica não-linear, é realizada empregando-se os métodos
de Newmark conjuntamente com o método de Newton-Raphson. O algoritmo de Newton-
Raphson é utilizado para a determinação da posição de equilíbrio levando em conta que a
rigidez não é mais constante e sim uma função dos deslocamentos nodais da estrutura, em
razão da não linearidade geométrica do elemento.
A seqüência desenvolvida no item 4.4.1 pode ser empregada com as devidas
adaptações para o caso não-linear, como mostrado a seguir.
4.5.1 MÉTODO DE NEWMARK PARA ANÁLISE DINÂMICA NÃO-LINEAR
Da mesma forma que para o caso linear deverá ser adotado um ponto de partida que
consiste na definição de 0 0&P ,D e 0D , de maneira que atenda as condições iniciais do
problema que estiver sendo estudado, assim como o método a ser empregado, o da
aceleração linear ou o da aceleração média. São necessários também as matrizes de rigidez
inicial, massa e amortecimento, que podem ser tomadas inicialmente iguais às do caso
linear. Portanto são conhecidos também a rigidez inicial ( )0K , massa ( )M e
amortecimento ( )C . Outro passo importante e independente dos demais é a escolha do
intervalo de incremento do tempo t∆ que deve garantir a estabilidade e precisão do método.
Em seguida se pode determinar a aceleração inicial dada pela equação
( )-10 0 0 0
ext= − −&& &D M P CD P (4.52)
onde 0P representa a força elástica (força interna) no instante (ou iteração) inicial.
O cálculo dos parâmetros a e b são feitos empregando as mesmas equações do caso
linear, ou seja, as equações (4.42) e (4.43).
O vetor da diferença das forças externas em dois instantes consecutivos é
1ext ext ext
i i i+∆ = −P P P . A variação da força externa equivalente é
ˆ exti i i i∆ = ∆ + +& &&P P aD bD (4.53)
87
Como a matriz de rigidez é dada pela derivada do vetor das forças internas em
relação ao vetor dos deslocamentos, ou seja pela equação (4.22), a rigidez tangente é dada
pela derivada da força equivalente em relação ao deslocamento e se escreve
( )21ˆ
i i t tγβ β
= + +∆ ∆
K K C M (4.54)
A figura 4.4 mostra a variação da força interna em função dos deslocamentos. As
seguintes considerações são os dados iniciais para o emprego do método de Newton-
Raphson com o objetivo de determinar a matriz de rigidez no tempo i(1)
i∆ = ∆R P (4.55)(0)
1i i+ =D D (4.56)( )0
i i=P P (4.57)ˆ ˆ
T i=K K (4.58)Para cada iteração com ,....3,2,1=j ,deve-se fazer
( ) -1 ( )ˆj jT∆ = ∆D K R (4.59)
( ) ( -1) ( )1 1j j j
i i+ += + ∆D D D (4.60)
( )( ) ( ) ( -1) ( )ˆj j j ji i i T i∆ = − + − ∆P P P K K D (4.61)
( 1) ( ) ( )j j ji
+∆ = ∆ −∆R R P (4.62)onde os passos estabelecidos pelas equações de (4.59) a (4.62) devem ser repetidos a cada
iteração, até ser atingida a precisão requerida.
(1)D
i n
D(i 1 )P
D(iP 2 )
D(P )
P
DDD Di+1
n(2)
)Dn+1(Pi
P(i D)ext
Pi+1^
Pi
Pi
iP
i i i
Figura 4.4 –Método de Newmark Não-Linear.
88
De posse de ˆiK e i∆P , determina-se i∆D e com o emprego das equações de número
(4.47) a (4.51) o problema é resolvido.
4.6 DESCRIÇÃO DO PROGRAMA COMPUTACIONAL
O programa utilizado para validação do elemento proposto foi elaborado na
plataforma Matlab.
Ele é composto por duas partes, uma chamada Programa Principal que processa os
dados fornecidos pela outra parte referente ao programa de entrada de dados, que possui a
forma de arquivo script (roteiro).
O programa de entrada de dados deve conter:
a) o número de nós, suas coordenadas, as restrições e massas nodais, e as cargas
atuantes nos nós da estrutura;
b) o número de elementos, a conectividade e os dados dos elementos
(compreendendo a área, o momento de inércia, o módulo de elasticidade do
material, a massa distribuída no elemento e os carregamentos distribuídos nas
direções longitudinal e transversal ao elemento);
c) os tipos de análise a serem realizadas pelo programa principal, ou seja, se será
feita uma análise estática, modal ou dinâmica da estrutura, ou ainda, se será feito
um estudo linear ou não-linear;
d) o intervalo de tempo e o tempo total de análise;
e) a tolerância (precisão requerida no problema) e um vetor descrevendo a variação
do carregamento no tempo.
De posse dos dados de entrada o programa principal é chamado e inicialmente
determina o número de graus de liberdade da estrutura, para em seguida montar o vetor de
carga.
Após esta etapa, é verificado pelo programa o tipo de análise a ser realizada, estática,
modal ou dinâmica.
Dentro das opções análise estática ou análise dinâmica é dada a opção de se fazer um
estudo linear ou não-linear geométrico.
89
No caso da análise estática, é utilizado o método de Newton-Raphson para a
determinação do vetor de deslocamentos e forças internas.
Quando do emprego do método de Newton-Raphson é chamada uma rotina para a
determinação do vetor de forças internas e da matriz de rigidez da estrutura, os quais são
obtidos somando-se a contribuição de cada elemento, de acordo com as equações (4.9) e
(4.12).
Para a análise linear considera-se a matriz de rigidez tradicional do elemento de
pórtico.
Caso seja feita a análise não linear, utiliza-se a formulação co-rotacional para
determinação do vetor de forças e matriz de rigidez do elemento.
Na montagem do vetor de forças internas e da matriz de rigidez de cada elemento é
empregado inicialmente o sistema básico, depois através de uma transformação de
coordenadas chega-se ao vetor de forças e a matriz de rigidez do elemento no sistema
global.
De posse do vetor de forças e da rigidez no sistema global soma-se a contribuição de
cada elemento e chega-se ao vetor de forças e a matriz de rigidez da estrutura.
É oportuno citar que a matriz de rigidez no sistema básico é obtida com a formulação
co-rotacional considerando a presença do carregamento distribuído.
Caso seja solicitada pelo usuário uma análise modal o programa principal também
empregará o método de Newton-Raphson estático, porém fornecerá alem dos vetores de
forças internas e deslocamentos a matriz de massa da estrutura, para que possa resolver o
problema de autovalor e determinar as freqüências modais e modos de vibração da estrutura.
A outra opção é o caso dinâmico, onde então o programa principal chamará a rotina
do método de Newmark, para a determinação da variação dos deslocamentos da estrutura no
tempo.
Para se obter o equilíbrio da força interna com o carregamento externo em cada
passo no tempo emprega-se o método de Newton-Raphson para o caso dinâmico.
Capítulo 5
Equation Section 5EXEMPLOS
Para verificação e comprovação do elemento proposto foi utilizado o programa
descrito no capítulo 4, item 4.6, fazendo-se a análise de diversos exemplos, os quais são
mostrados a seguir.
Os exemplos foram executados em computador Pentium III, com 384 MB de
memória Ram, HD de 20GB e processador de 1,7 GB.
5.1 TESTE DE VALIDAÇÃO PARA MATRIZ DE RIGIDEZ
Para se fazer uma verificação da dedução dos termos da matriz de rigidez da
estrutura, pode-se determinar numericamente as derivadas referentes aos termos da
matriz de rigidez, por diferenças finitas, obtendo-se assim uma matriz secante.
Para isto, considera-se um elemento de barra com comprimento indeformado de
1,20 m, base da seção igual a 20 cm, altura igual a 2 cm, módulo de elasticidade igual a121, 2 10× kgf/m² e sujeito a uma carga distribuída de 1000 kgf/m.
Para o cálculo das matrizes é imposto um vetor de deslocamentos igual a
2,2 1,3 5,2 3,5 6,3 4,8 TG = −d .
Calculando as derivadas referentes a matriz de rigidez analiticamente, obtém-se o
seguinte resultado9 8 5 9 8 53.3132 10 3.4348 10 -1.2800 10 -3.3132 10 -3.4348 10 -1.2800 108 9 4 8 9 43.4348 10 3.8282 10 6.4000 10 -3.4348 10 -3.8282 10 6.4000 105 4 5 5 4 5-1.2810 10 6.4000 10 5.3333 10 1.2810 10 6.4000 10 - 2.6667 10
-3.31gk
× × × × × ×
× × × × × ×
× × × × × ×= 9 8 5 9 8 532 10 -3.4348 10 1.2800 10 3.3132 10 3.4348 10 1.2800 108 9 4 8 9 4-3.4348 10 -3.8282 10 -6.4000 10 3.4348 10 3.8282 10 -6.4000 105 4 5 5 4 5-1.2790 10 6.4000 10 2.6667 10 1.2790 10 -6.4000 10 5.3333 10
× × × × × ×
× × × × × ×
× × × × × ×
(5.1)
Calculando estas derivadas numericamente, obtém-se9 8 5 9 8 53.3132 10 3.4348 10 -1.2798 10 -3.3132 10 -3.4348 10 -1.2808 108 9 4 8 9 43.4348 10 3.8282 10 6.3896 10 -3.4348 10 -3.8282 10 6.3896 105 4 5 5 4 5-1.2810 10 6.4000 10 5.3333 10 1.2810 10 6.4000 10 - 2.6667 10
-3numk
× × × × × ×
× × × × × ×
× × × × × ×= 9 8 5 9 8 5.3132 10 -3.4348 10 1.2798 10 3.3132 10 3.4348 10 1.2808 108 9 4 8 9 4-3.4348 10 -3.8282 10 -6.3896 10 3.4348 10 3.8282 10 -6.3896 105 4 5 5 4 5-1.2790 10 6.4000 10 2.6667 10 1.2790 10 -6.4000 10 5.3333 10
× × × × × ×
× × × × × ×
× × × × × ×
(5.2)
91O erro pode ser representado na forma de uma matriz, dado pela diferença
entre os resultados (5.1) e (5.2) dividido por (5.1), que resulta em-9 -7 -4 -8 -7 -4 6.1431 10 5.8950 10 -1.3208 10 -2.2641 10 -5.2110 10 6.1297 10-8 -7 -3 -7 -7 -3 3.4198 10 1.4507 10 -1.6222 10 -5.2110 10 -1.0404 10 -1.6222 10-7 - -2.0106 10 -3.8906 10
dif
× × × × × ×
× × × × × ×
× ×=
7 -8 -7 -7 -8 -1.0710 10 -2.0106 10 -3.8906 10 -1.0710 10-9 -7 -4 -8 -7 -4 6.1431 10 5.8950 10 -1.3208 10 -2.2641 10 -5.2110 10 6.1297 10-8 -7 -3 -7 3.4198 10 1.4507 10 -1.6222 10 -5.2110 10
× × × ×
× × × × × ×
× × × × -7 -3 -1.0404 10 -1.6222 10-7 -8 -8 -7 -8 -8 1.5081 10 -2.5262 10 -1.0710 10 1.5081 10 -2.5262 10 -1.0710 10
× ×
× × × × × ×
(5.3)
A listagem do programa empregado para o cálculo da matriz de rigidez através
das diferenças finitas se encontra no anexo 3 deste trabalho.
5.2 VIGA BI-APOIADA COM CARGA DISTRIBUÍDA
Este exemplo consiste na análise não linear estática de uma viga bi-apoiada
submetida a um carregamento transversal distribuído. A viga possui comprimento de 50
cm, base de 100 cm e altura de 1 cm, conforme ilustra esquematicamente a figura 5.1. O
módulo de elasticidade do material da viga é igual a 220.000 Mpa.
50 cm
W
v
E
y
=
=
220000
300
MPa
MPaσ
��������
1 cm100 cm
������� �������
Figura 5.1 – Viga bi-apoiada sujeita a carregamento uniformemente distribuído.
(De Souza/2000)
Este exemplo consta no trabalho de De Souza (2000), tendo sido examinado
anteriormente por Backlund (1976), que fez uma análise numérica considerando
diferentes combinações das não-linearidades (geométrica e física) e por Coulter e Miller
(1988), que também consideraram as duas não-linearidades. Uma solução analítica para
o problema foi apresentada por Timoshenko e Woinowsky-Krieger (1959), conforme
informa De Souza (2000).
92No presente estudo, faz-se uma análise da convergência considerando 1, 2, 4, 8
e 16 elementos para metade da viga, em função da simetria do problema. É empregada a
formulação co-rotacional com a utilização do carregamento distribuído aplicado ao longo
de toda a barra.
O carregamento é aplicado em 20 passos de carga, até atingir o valor total de 2
MN/m. Os deslocamentos obtidos no meio do vão para cada passo de carga, referentes
aos diferentes graus de discretização, assim como a solução analítica do problema são
mostrados na figura 5.2.
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.0160
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Deslocamento vertical no meio do vão(m)
Inte
nsid
ade
do c
arre
gam
ento
W (
MN
/m)
Nels = 1 Nels = 2 Nels = 4 Nels = 8 Nels = 16 Analitico - Timoshenko(1959)
Figura 5.2 – Curva Carga x Deslocamentos para diferentes níveis de discretização.
Nota-se também que devido à presença de dois apoios do segundo gênero, os
quais impedem o movimento das extremidades da viga, introduzindo um esforço de
tração, a rigidez da estrutura aumenta em função do carregamento, demonstrando a
característica não-linear do problema.
Cabe aqui lembrar que a não-linearidade não está somente associada a grandes
deslocamentos, estando principalmente relacionada com a falta de proporcionalidade
entre deslocamentos e deformações.
93Observa-se que com 8 elementos o valor obtido através da formulação
proposta já é bastante próximo do resultado analítico.
O tempo gasto para obtenção dos resultados pelo Matlab é de aproximadamente 6
segundos, mostrando a rapidez de processamento do programa apresentado.
5.3 COMPORTAMENTO ESTÁTICO DO CABO 1
Este exemplo consiste na análise estática de um cabo pré-tensionado, com
comprimento de 312,7m, submetido a um carregamento distribuído transversal de 46,11
N/m, correspondente a seu peso próprio, e a uma força de tração aplicada nas
extremidades e igual a 17794 N, conforme ilustra a figura 5.3. O módulo de elasticidade
do cabo é igual a 111,31 10× Pa, a área da seção transversal é igual a 45, 48 10−× m², e o
momento de inércia igual a 12 41,4 10 m−× .
312,7 m
PPq
Figura 5.3 – Cabo pré-tensionado com carregamento distribuído transversal.
O objetivo deste exemplo é mostrar a precisão do elemento proposto, realizando
uma análise da convergência para o valor do deslocamento no meio do vão do cabo. Faz-
se uma comparação entre os resultados obtidos considerando o peso próprio como
carregamento distribuído propriamente dito, e como cargas concentradas aplicadas nos
nós.
Foi realizada uma análise estática não-linear considerando-se 1, 2, 4, 8, 16 e 32
elementos.
A figura 5.4 mostra os resultados obtidos para o deslocamento no meio do vão,
em função do número de elementos empregados.
94Na figura mostra-se também o resultado analítico apresentado por Labegalini
et al (1992), cujo valor é dado por
0
0
cosh 12
cT wLyw T
= −
(5.4)
onde 0T é a componente horizontal da tração no cabo, w é o carregamento distribuído e
cL é o comprimento deformado do cabo, dado por
11
2 senh2
Apc
LL C
C
=
(5.5)
onde ApL é a distância entre os apoios do cabo e 1C é dado por
01
TCw
= (5.6)
Observando a figura 5.4, percebe-se que quando se utiliza o carregamento
distribuído com poucos elementos o resultado está muito longe da resposta correta.
Quando se emprega a carga concentrada a convergência é mais rápida. Com 8 elementos
as duas formulações já apresentam resultados praticamente iguais.
0 5 10 15 20 25 30 3530
31
32
33
34
35
36
37
38
39
Número de elementos
Fle
cha
no m
eio
do v
ão d
a vi
ga
Carga DistribuídaCarga ConcentradaAnalitico
Figura 5.4 – Análise de convergência da flecha do cabo
95Estes resultados são inesperados, pois a consideração do peso próprio como
carga distribuída deveria apresentar melhores resultados do que, considerando-se cargas
concentradas.
Esta discrepância é justificada pelo fato do comportamento do cabo ser
essencialmente axial, com uma rigidez a flexão praticamente desprezível, o que introduz
erros numéricos na solução do problema.
5.4 VIGA ENGASTADA E LIVRE
Este exemplo consiste na análise não-linear estática de uma viga engastada e livre
com um vão de 3 metros, base de 15 cm e altura de 2 cm, sujeita a ação de seu peso
próprio. O material da viga possui módulo de elasticidade de 2 MPa e peso específico de
25 N/m³.
Assim como no exemplo anterior, realiza-se uma análise da convergência para o
valor do deslocamento na extremidade da viga. Faz-se também, uma comparação entre
os resultados obtidos considerando o peso próprio como carregamento distribuído
propriamente dito, e como cargas concentradas aplicadas nos nós.
A discretização é realizada com 1, 2, 4, 8 e 16 elementos.
A figura 5.5 mostra os resultados obtidos.
Em função dos resultados obtidos, observa-se que este exemplo corresponde a
um problema altamente não linear, com deslocamentos de mesma ordem de grandeza do
comprimento da viga.
Pode-se observar que no caso da viga, os resultados obtidos com o emprego do
carregamento distribuído são melhores do que os resultados obtidos considerando-se
cargas concentradas, ao contrário do que obtido para o cabo do exemplo 5.3.
Este resultado serve como uma confirmação sobre os comentários feitos no final
da seção anterior. No problema em questão, qual seja de uma viga engastada e livre, o
efeito de flexão da estrutura é preponderante em relação ao efeito axial. Provavelmente,
por isso, neste caso a taxa de convergência referente ao carregamento distribuído seja
melhor do que à referente às cargas concentradas.
96
0 2 4 6 8 10 12 14 162.1
2.15
2.2
2.25
2.3
2.35
2.4
2.45
2.5
Número de elementos
Fle
cha
na e
xtre
mid
ade
do b
alan
ço d
a vi
ga (
m)
Carga DistribuídaCarga Concentrada
Figura 5.5 – Flecha na extremidade do balanço de uma viga engastada e livre.
5.5 COMPORTAMENTO ESTÁTICO DO CABO 2
Este exemplo consiste em uma nova análise do cabo apresentado na seção 5.3.
Nesta nova análise, aumenta-se a força de tração do cabo, e verifica-se a diminuição da
flecha. A força inicial de tração considerada é 17794T N= , a qual é posteriormente
multiplicada por 2, 3, 4, …,20. A distância entre apoios para o valor da força de tração
inicial T é igual a 304,8m, conforme ilustra a figura 5.6.
T T
304,8 m
Figura 5.6 – Cabo sujeito a esforço de tração nos apoios.
97Este estudo foi inicialmente apresentado por Ali e Ghaffar (1995) e por
Oliveira (2002), que empregou elementos isoparamétricos de cabos com quatro nós.
Oliveira (2002) obteve a resposta estática utilizando os elementos Lagrangeanos
atualizado e total, além de ter empregado também o Elemento de Catenária, por ele
desenvolvido.
É importante esclarecer que Oliveira (2002) obteve o mesmo resultado tanto para
a formulação Lagrangeana atualizada, como para a Lagrangeana total e para o Elemento
de catenária.
É registrada para cada valor da força de tração a flecha relativa no meio do vão
do cabo conforme mostrado na figura 5.7.
Conclui-se, através da figura 5.7, que com o aumento da força de tração no cabo,
a flecha obtida no meio do vão tende a reduzir. Isto demonstra o efeito da parcela de
rigidez da estrutura referente a matriz de rigidez geométrica, a qual aumenta em função
de aumento no esforço de tração.
Observa-se que os resultados obtidos pelo método apresentado são bem próximos
dos resultados encontrados por Oliveira (2002).
0 2 4 6 8 10 12 14 16 18 200
1
2
3
4
5
6
7
8
9
10
Incremento da força de tração no cabo (T/To)
(Fle
cha/
vão
do c
abo)
x100
Método apresentadoOliveira(2002)
Figura 5.7 – Redução da flecha com o aumento da força de tração.
98O tempo de processamento foi de 49 segundos. Este tempo corresponde a
todas as 20 iterações efetuadas para o acréscimo da força de tração.
5.6 PÓRTICO DE LEE
Este exemplo consiste no estudo estático não-linear do pórtico de Lee – um
exemplo clássico apresentado na literatura sobre análise estrutural geometricamente não-
linear.
O pórtico é constituído por duas barras, cada uma com comprimento de 120 cm e
seção transversal com 3 cm de base e 2 cm de altura. As barras são constituídas por
material com módulo de elasticidade igual a 70608 MPa. A figura 5.8 mostra o pórtico
de Lee.
De Souza (2000) efetuou o estudo do comportamento não-linear geométrico
considerando inicialmente o material linear elástico e em seguida o material como elasto-
plástico.
Neste trabalho é feito o estudo apenas para o caso não-linear geométrico com
material linear elástico. Neste trabalho foi empregada uma discretização correspondente
a 10 elementos por barra, no total de 20 elementos e 21 nós.
Este é o mesmo número de elementos empregados por Cichon (1983) conforme
mostra De Souza (2000).
E
E EH
y
=
=
=
70608
01
1020
MPa
MPa
.
σ
120
cm
P,v
120 cm
96 cm24 cm ������
2 cm3 cm
��������������
��������������
Figura 5.8 - Pórtico de Lee. (De Souza/2000)
99A análise realizada não considerou o carregamento distribuído (peso próprio
das barras), considerando-se somente a carga concentrada.
Os resultados encontrados mostram que a estrutura sofre grandes deslocamentos
comprovando a sua grande não-linearidade geométrica.
A figura 5.9 mostra os resultados obtidos com a formulação proposta e os
resultados obtidos por De Souza (2000), para o caso da análise não-linear geométrica
considerando o material linear elástico.
A partir da figura 5.9, conclui-se que os resultados obtidos utilizando-se o
método apresentado são bem próximos dos resultados encontrados por De Souza (2000).
Entretanto, deve-se destacar que a trajetória de equilíbrio apresentada foi interrompida
no primeiro ponto limite, em função da impossibilidade de se ultrapassar este ponto
utilizando-se o Método de Newton-Raphson. Na verdade, a trajetória de equilíbrio
completa deste problema é bem complexa, pois a estrutura apresenta os fenômenos de
snap-back e snap-through.
A figura 5.10 mostra a trajetória de equilíbrio completa. A figura apresentada foi
gerada por De Souza (2002).
O tempo de processamento foi de 7 segundos.
0 5 10 15 20 25 30 35 40 45 500
2
4
6
8
10
12
14
16
18
20
Deslocamento(cm)
Car
ga A
plic
ada
(kN
)
Método apresentadoDe souza
Figura 5.9 – Deslocamento vertical do ponto de aplicação da carga para o pórtico de Lee.
100
0 10 20 30 40 50 60 70 80 90 100-10
-5
0
5
10
15
20
linear elastic
elasto-plasticwith kinematichardening
Lee`s Frame
Displacement v (cm)
App
lied
Load
P (k
N)
proposed flex. formul. - 3 elmtsCichon (1983) - 10 elmts
Figura 5.10 – Trajetória completa de equilíbrio do pórtico de Lee (De Souza/2002)
5.7 BARRA COM UM GRAU DE LIBERDADE
Este exemplo corresponde a análise dinâmica linear de um sistema com um grau
de liberdade constituído por uma barra com seção transversal de 21 m , momento de
inércia de 41/12 m e módulo de elasticidade de 210N/m .
O único deslocamento permitido para a barra é o deslocamento horizontal no nó
2, conforme mostra a figura 5.11.
100 cm
P(t) NP(t)
1010 sen( t/0,6)~111
t0,6
Figura 5.11 – Barra com um grau de liberdade.
101A massa da barra é de 0,2533 kg e o coeficiente de amortecimento
0,1592 Ns/mc = . O passo de tempo empregado para a análise é de 0,1 s, em um tempo
total de 1 s e a força aplicada é do tipo senoidal.
É utilizado para a solução do problema o algoritmo de Newmark, método da
aceleração média, onde 0, 25β = e 0,5γ = .
A figura 5.12 mostra os deslocamentos horizontais do nó 2 determinados pelo
elemento de barra linear implementado e os resultados apresentados por Chopra (1995).
Este exemplo teve por objetivo apenas a validação do programa para sistemas
lineares amortecidos. Observa-se que os resultados obtidos pelo programa apresentado e
os teóricos de Chopra (1995) são coincidentes.
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-0.5
0
0.5
1
1.5
Intervalo de tempo
Des
loca
men
to h
oriz
onta
l no
nó 2
E lemento propostoChopra
Figura 5.12 – Deslocamento horizontal de barra sujeita a uma força senoidal.
5.8 TORRE CARREGADA NO TOPO
Chopra (1995) faz a análise dinâmica linear de uma torre com 600 pés de altura,
empregando 5 elementos.
A figura 5.13 mostra a torre e a carga aplicada na extremidade livre.
102
t
1000
P(t) kipsP
600 ft
Figura 5.13 – Torre sujeita a carga degrau.
Na extremidade de cada um dos quatro primeiros elementos é considerada uma
massa concentrada de 208,6 kip-s²/ft. No final do quinto elemento a massa é igual a
metade da massa dos demais nós, conforme indica a figura 5.14.
m
m
m
m
m/2
120
120
120
120
120
Figura 5.14 – Torre e massas concentradas.
A rigidez da seção da torre é 10 2EI=5,469 10 kip-ft× e o intervalo de tempo
utilizado na análise é de 0,1 s.
Uma extremidade da torre é considerada engastada e a outra livre. Na
extremidade livre é aplicada uma força degrau horizontal de intensidade igual a 1000
kips.
Os resultados encontrados pelo método apresentado e os obtidos por Chopra
(1995) encontram-se na figura 5.15.
Observa-se que os resultados obtidos são coincidentes.
103
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20
0.5
1
1.5
2
2.5
3
tempo(s)
Des
loca
men
to h
oriz
onta
l na
extr
emid
ade
do b
alan
ço(p
és)
Método apresentadoChopra(1995)
Figura 5.15 – Deslocamentos na extremidade da torre.
5.9 VIGA BIAPOIADA
Este exemplo consiste em uma análise dinâmica linear de uma viga bi-apoiada,
submetida a uma carga degrau aplicada no meio do vão.
A figura 5.16 mostra a viga e o carregamento atuante em seu ponto médio.
t
2846,72
P(t) N P(t)
50 cm 50 cm
Figura 5.16 – Viga bi-apoiada sujeita a carga degrau.
104A viga possui base de 15 cm, altura de 40 cm, comprimento total de 100 cm,
massa de 15 kg e o módulo de elasticidade do material da viga é 112,0685 10× N/m².
Este problema é resolvido analiticamente por Chopra (1995). A solução
apresentada corresponde aos deslocamentos verticais do ponto médio da barra.
A figura 5.17 mostra o resultado obtido para o deslocamento vertical da viga no
meio do vão. Nesta figura mostram-se tantos os resultados obtidos pelo programa
elaborado e aqueles encontrados com o emprego da resposta analítica apresentado por
Chopra (1995).
Não foi considerado o efeito do peso próprio da viga, ou seja, não está sendo
empregado o carregamento distribuído neste exemplo.
Pode-se observar que a viga apresenta pequenos deslocamentos em seu ponto
médio, e que os resultados obtidos numericamente aproximam satisfatoriamente a
solução analítica apresentada por Chopra (1195).
0 0.005 0.01 0.015 0.02 0.025 0.030
0.5
1
1.5
2
2.5
3
3.5
4
4.5x 10
-3 Desloc. Vertical no meio do vão
Tempo
Des
loca
men
to V
ertic
al (
cm)
Analítico Chopra Método apresentado
Figura 5.17 – Deslocamento vertical no ponto médio da viga bi-apoiada sujeita a carga
degrau.
105
5.10 VIGA BI-ENGASTADA
Neste exemplo faz-se uma análise dinâmica não-linear geométrica da viga bi-
engastada, sujeita a uma carga degrau concentrada no meio do vão, analisada por
Mondkar e Powell (1977), e Yang e Saigal (1984).
A viga possui comprimento total de 20 polegadas, seção transversal retangular,
com base de 1 polegada e altura de 1/8 de polegada. O módulo de elasticidade do
material da viga é 73 10 psiE = × , o coeficiente de Poisson 0ν = e a massa20,0026lb×s/inq = . O método de Newmark é empregado para 0, 25β = e 0,5γ = . A
força externa aplicada é do tipo degrau com valor de 640 lbf. A figura 5.18 mostra a viga
e o carregamento aplicado.
É importante ressaltar que as unidades estão sendo empregadas no sistema inglês
com o objetivo de se ter uma boa comparação com o resultado apresentado por Yang e
Saigal (1984).
10 int 10 in
640
P(t) lb P(t)
Figura 5.18 – Viga bi-engastada sujeita a carga degrau.
A figura 5.18 mostra os deslocamentos verticais do nó central, onde está sendo
aplicada a carga degrau de 640 lbf. A análise foi realizada com o emprego de 10
elementos.
A resposta obtida corresponde ao emprego de 10 elementos e com intervalo de
tempo de 200 sµ e um tempo total de análise de 50000 sµ .
Yang e Saigal (1984) realizaram o estudo discretizado da viga com o emprego de
3, 4 e 6 elementos, sendo mostrado na figura 5.19, o gráfico referente a 6 elementos. É
106importante destacar que o estudo desta viga foi realizado inicialmente por Mondkar e
Powell (1977).
Os resultados apresentados por Yang e Saigal (1984) e Mondkar e Powell (1977),
apresentam um erro em seus gráficos, pois ambos mostram a escala de tempo um décimo
menor que a real.
A confirmação da falha de ambos foi feita através da verificação do gráfico
apresentado por Mondkar e Powell (1977) para a análise linear da mesma viga que
também apresentou o mesmo problema no eixo do tempo.
O resultado apresentado na figura 5.18 referente à solução de Yang e Saigal
(1984) teve a escala de tempo corrigida (multiplicada por 10), a fim de proporcionar a
correta comparação com os resultados obtidos através do método proposto. Observa-se
que os resultados apresentam boa concordância, com pequenas diferenças provavelmente
causadas pelo pequeno grau de discretização utilizado, ou pelo fato do passo de tempo
considerado não ser pequeno o suficiente.
A figura 5.20 mostra o comportamento linear apresentado por Mondkar e Powell
(1977).
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.050
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
Tempo (s)
Des
loca
man
to V
ertic
al (
in)
Desloc. Vertical no meio do vão - Caso Não LInear
Elemento PropostoYang e Saigal
Figura 5.19 – Resposta dinâmica da viga bi-engastada.
107
O resultado dos deslocamentos do ponto médio da viga apresenta da mesma
forma que para o caso não-linear, uma diferença no tempo de 1/10.
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.050
2
4
6
8
10
12
Tempo
Des
loca
man
to V
ertic
al (
in)
Desloc. Vertical no meio do vão - Caso Linear
Formulação Co-rotacionalMondkar e Powell
Figura 5.20 – Resposta dinâmica linear da viga bi-engastada.
Como o exemplo da viga bi-apoiada com a carga degrau, para o caso linear já foi
verificado, observa-se que os artigos de Yang e Saigal (1984) e Mondkar e Powell
(1977) apresentam esta falha.
Com fins apenas ilustrativos, a figura 5.21 mostra tanto o resultado linear como o
não-linear para a análise dinâmica da viga estudada.
O tempo total de processamento foi de 58 segundos. Neste intervalo de tempo
foram realizadas tanto a análise linear como a análise não-linear da viga em estudo.
108
0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.050
2
4
6
8
10
12Desloc. Vertical no meio do vão
Tempo
Cas
o Li
near
x 1
- C
aso
Não
Lin
ear
x 10
Caso Linear Caso Não-linear
Figura 5.21 – Resposta dinâmica da viga bi-engastada, casos linear e não-linear.
5.11 ANÁLISE MODAL DO CABO
O último exemplo deste trabalho consiste na determinação das freqüências
naturais do cabo dos exemplos 5.3 e 5.5.
Neste exemplo utiliza-se o programa para a determinação numérica das
freqüências naturais do cabo e faz-se uma comparação entre os resultados obtidos através
do programa com valores teóricos, considerando a teoria de corda vibrante desenvolvida
no capítulo 2.
A análise foi realizada empregando-se 2, 4, 8, 16 e 32 elementos.
A expressão utilizada para a determinação dos valores teóricos foi
2n TfL m
= (5.7)
onde n é o modo de vibração, L o comprimento do cabo, T a força de tração no cabo e m
a massa por unidade de comprimento do cabo. Esta fórmula foi deduzida considerando-
109se que o cabo é inicialmente reto, que o esforço de tração é constante ao longo do
comprimento, e que a sua rigidez a flexão é nula.
Jordan (1985) apresenta uma expressão que leva em consideração a rigidez à
flexão do cabo no cálculo das freqüências naturais. A expressão é
2 2
2
12nn n EIf TL m L
π = +
(5.8)
onde além dos termos especificados na equação (5.7) tem-se E que corresponde ao
módulo de elasticidade do material do cabo e I que é o momento de inércia.
Como a inércia do cabo é pequena os resultados apresentados pelas equações
(5.7) e(5.8) são praticamente iguais.
A precisão dos resultados numéricos depende do grau de discretização da
estrutura.
Os resultados obtidos para os seis primeiros modos de vibração, com o emprego
de 32 elementos, assim como as diferenças entre os valores do programa e os teóricos
encontram-se na tabela 5.1.
Modo Elemento Proposto (Hz) Freqüência Natural (Hz)
(Analítico)
Diferença
1 0,3105 0,3111 -0,0004
2 0,6949 0,6222 0,0727
3 0,9845 0,9333 0,0512
4 1,2445 1,2361 0,0084
5 1,5556 1,5211 0,0345
6 1,8667 1,8413 0,0254
Tabela 5.1 – Valores das freqüências naturais do cabo.
As diferenças encontradas devem-se principalmente ao fato de que o valor
teórico não leva em conta a variação do esforço normal ao longo do seu comprimento.
A figura 5.22 mostra a convergência da freqüência do primeiro modo de
vibração, calculada pelo programa em função do número de elementos utilizados na
discretização.
110
0 5 10 15 20 25 30 350.31
0.32
0.33
0.34
0.35
0.36
0.37
0.38
Número de Elementos
Fre
quên
cia
natu
ral (
Hz)
P rimeiro modo de Vibração
Elemento PropostoAnalítico
Figura 5.22 – Convergência do primeiro modo de vibração.
A figura 5.23 mostra a convergência da freqüência do segundo modo de vibração,
calculada pelo programa em função do número de elementos utilizados na discretização.
0 5 10 15 20 25 30 350.5
1
1.5
2
2.5
3
3.5
Número de Elementos
Fre
quên
cia
natu
ral (
Hz)
Segundo modo de Vibração
Elemento PropostoAnalítico
Figura 5.23 – Convergência do segundo modo de vibração.
111
A figura 5.24 mostra a convergência da freqüência do terceiro modo de vibração,
calculada pelo programa em função do número de elementos utilizados na discretização.
0 5 10 15 20 25 30 350
2
4
6
8
10
12
Número de Elementos
Fre
quên
cia
natu
ral (
Hz)
Terceiro modo de Vibração
Elemento PropostoAnalítico
Figura 5.23 – Convergência do terceiro modo de vibração.
Capítulo 6
CONCLUSÕES E SUGESTÕES
6.1 CONCLUSÕES
Uma das grandes preocupações das engenharias civil e mecânica é a construção e/ou
fabricação de estruturas que atendam as suas finalidades de forma segura e econômica.
Para que se tenha segurança sobre o que se está produzindo é necessário um bom
domínio sobre a estrutura que se está projetando.
Quando se trata de estruturas mecânicas, tais como máquinas em geral, obras civis, e
cabos das linhas aéreas de transmissão, é importante ter-se conhecimento de como elas se
comportam diante das mais diversas solicitações a que estão sujeitas.
Embora várias estruturas possam ser analisadas considerando pequenas deformações
e carregamentos estáticos, existem casos onde as estruturas são sujeitas a grandes
deformações e/ou carregamentos dinâmicos.
Diante do exposto, justifica-se a necessidade e importância do trabalho aqui
apresentado.
Os resultados obtidos para os exemplos apresentados no texto comprovam a
eficiência do elemento proposto na solução destes tipos de problemas. Além disso, ficou
demonstrada a facilidade do elemento proposto para o tratamento de grandes rotações.
Na fase de implementação computacional, foi observada a facilidade da elaboração
do programa na plataforma Matlab. O Matlab apresenta diversas rotinas já implementadas, o
que ajudou bastante o desenvolvimento do programa.
De uma forma geral pode-se dizer que o objetivo principal desta dissertação foi
atingido, qual seja, o estudo da formulação e implementação computacional de um elemento
finito para análise dinâmica não-linear geométrica para aplicação em estruturas constituídas
de barras e estruturas formada por cabos, com desenvolvimento da teoria para consideração
de carregamento distribuído na formulação co-rotacional.
Também foram apresentados exemplos com aplicação a cabos de linhas aéreas de
transmissão de energia, os quais como foi visto, são submetidos a grandes deformações.
113
De uma maneira geral pode-se relacionar as seguintes conclusões:
a) O emprego da formulação co-rotacional permitiu tratar com simplicidade as
grandes deformações a que as estruturas analisadas estavam sujeitas. Neste
método cada elemento é tratado como linear em um sistema auxiliar
denominado básico e a não-linearidade é introduzida na transformação de
coordenadas entre este sistema básico e o sistema global.
b) Na resolução de sistemas de equações não-lineares foi empregado o método de
Newton-Raphson. Existem duas versões para o método, o tradicional e o
método de Newton-Raphson modificado.
O método tradicional, quando aplicado a problemas mecânicos, é caracterizado
por adotar uma matriz de rigidez que é atualizada a cada passo de carga. O
método modificado efetua as iterações empregando a mesma matriz de rigidez
para cada passo de carga.
Nesta dissertação foi empregado o método tradicional por se acreditar que para
as estruturas de interesse ele deva apresentar uma convergência mais rápida.
c) O método de Newmark foi empregado em razão da facilidade de
implementação computacional e da rapidez e precisão. Dentre os métodos de
Newmark foi utilizado o da aceleração média devido ao fato deste ser
incondicionalmente estável.
d) A grande dificuldade encontrada no trabalho foi a dedução e validação da
consideração do carregamento distribuído para a formulação co-rotacioanal.
É oportuno registrar que não se tem conhecimento de trabalhos sobre o assunto.
e) O programa elaborado pode ser perfeitamente empregado em diversas
situações práticas onde se tenha a não-linearidade geométrica suplantando os
efeitos da não-linearidade física, ou seja, problemas em que grandes
deslocamentos alteram basicamente a forma da estrutura ainda com
comportamento elástico.
f) A matemática envolvida na dissertação é bastante interessante, principalmente
no aspecto de aplicação do princípio dos trabalhos virtuais e a sua relação com
o cálculo das variações.
114
As freqüências naturais obtidas com o emprego do programa desenvolvido
neste trabalho apresentam diferença quando comparadas aos valores teóricos.
Tais diferenças são decorrente da não consideração pelas fórmulas
tradicionalmente empregadas nos textos, da rigidez à flexão dos cabos, assim
como da variação do esforço axial de tração do cabo ao longo do seu
comprimento.
g) Quando se trabalha com carregamento distribuído se obtém uma convergência
melhor que aquela apresentada pelo carregamento concentrado. No caso dos
cabos tal fato não ocorre, provavelmente em razão da pequena rigidez à flexão
destes.
h) Nos exemplos apresentados foi verificado que o programa elaborado
apresentou resultados excelentes quando comparados com outros existentes na
literatura.
i) Embora não apresentado no texto do trabalho, verificou-se também a
convergência quadrática dos deslocamentos em função da utilização do método
clássico de Newton-Raphson.
6.2 SUGESTÕES PARA TRABALHOS FUTUROS
Para trabalhos futuros são apresentadas as seguintes sugestões:
a) A consideração da teoria de viga de Timoshenko utilizando o método das
forças, e verificação da sua influência nos valores das altas freqüências naturais
de vibração dos cabos.
b) A representação do cabo como elemento sem rigidez a flexão (elemento de
cabo, propriamente dito).
c) A consideração do efeito da ação aleatória do vento nos cabos das linhas aéreas
de transmissão de energia elétrica;
d) O estudo da não-linearidade geométrica de barras com seção variável;
e) A implementação espacial do elemento apresentado nesta dissertação;
f) O estudo da não-linearidade física;
ANEXOS
Anexo 1 - Programa Principal para Análise de Pórticos
% Universidade Federal do Para - UFPa% Centro Tecnológico - CT% Departamento de Engenharia Mecânica - DEM% Curso de Mestrado em Engenharia Mecânica% Dissertação% Programa didático para Analise Estática, Dinâmica, Linear e Não Linear de Pórticos% Autores: Remo Magalhães de Souza ([email protected])% Gandhy Yeddo da Rocha Aranha Júnior ([email protected])% Versão 1.0 Data: 10/2003% programa principal para analise - AnalisePort.m% CARACTERÍSTICAS DO PROGRAMA:% Analise Estática e Dinâmica, Linear e Não Linear de Pórticos Planos% Formulação Co-rotacional para Matriz de Rigidez% Numeração dos gdl: inicialmente os gdl livres e em seguida os gdl com restrição% Nos resultados a numeração dos nós é a numeração tradicional% obs_1: gdl = graus de liberdade% obs_2: Este Programa emprega a matriz de massa concentrada e consistente% function [f,DT,RA,ti] = AnalisePort01(numNos, numElems, coords, restrs, cargasNos,...% massaNos, conect,dadosElem, tipoEst_Din,...% tipoLin_NL, lambda, delta_t, tempo, tol);% versão 03/10/2003 - unificação das analises estática dinâmica, linear e não-linear% Dados de entrada:% numNos = numero de nos% numElems = numero de elementos% coords = coordenadas dos nos% restrs = restrições nodais% cargasNos = cargas nodais% massaNos = massa nodais% conect = conectividade dos elementos (incidência nodal)% dadosElem = dados do elemento [E A I qx qy];% tipoEst_Din = 1 -> Analise Estática% = 2 -> Analise Dinâmico% tipoLin_NL = 1 -> Analise Linear% = 2 -> Analise Não-Linear% tol = tolerância% ndp = número de incrementos da força externa% Parâmetros de saída:% RA = vetor de forças nodais inclusive reações de apoio% DT = vetor de deslocamentos nodais% f = freqüências naturais da estrutura
123
% ti = vetor de tempo para construção dos gráficos% Parâmetros da estrutura% nGrLib = numero de graus de liberdade da estruturafunction [DT,RA,ti,freq, Modos] = AnalisePort01(numNos, numElems, coords, restrs,...
cargasNos,massaNos, conect, dadosElem, tipoEst_Din,...tipoLin_NL, lambda, delta_t, ndp, tol, tempo);
% especifica parâmetros da estruturaNGLN = 3; % numero de graus de liberdade por no% obtém os graus de liberdade da estruturanGrLibEst = NGLN * numNos;[ngdll, gdl] = obtemGrausLibEstrut(numNos, restrs);% monta vetor de carga% P = carga externafor n = 1 : numNos
for m = 1 : 3P(gdl(n,m),1) = cargasNos(n,m);
endend% Tipo da Analisefreq = 0;Modos = 0;switch tipoEst_Din
case 1 % Estatico% emprego do método de Newton-Raphson para o caso Estático[Pint, Desl] = NewtonRaphsonEstatico(nGrLibEst, numNos, numElems, coords, restrs,...
cargasNos,massaNos, conect, dadosElem,...tipoLin_NL,lambda, ndp, tol, ngdll, gdl, P);
ti = 0:1/ndp:ndp;case 2 % Modal% emprego do método de Newton-Raphson para o caso Estático[Pint, Desl, K] = NewtonRaphsonEstatico(nGrLibEst, numNos, numElems, coords,...
restrs, cargasNos,massaNos, conect,...dadosElem, tipoLin_NL,lambda, ndp, tol,...ngdll, gdl, P);
M = obtemMatMassaEstrutura(nGrLibEst, numNos, numElems, coords, massaNos,conect,..dadosElem, gdl);
ti = 0:lambda:lambda*ndp;[freq, Modos] = AnaliseModal(nGrLibEst, ngdll, K, M);case 3 % Dinâmico% análise dinâmica através do método de newmark[Pint, Desl, ti] = newmark (nGrLibEst, numNos, numElems, coords, restrs, cargasNos,...
massaNos, conect, dadosElem, tipoLin_NL,lambda,ndp,...tol, ngdll, gdl, P, delta_t, tempo)
end% transforma para a numeração tradicional de graus de liberdadeDr = zeros (nGrLibEst-ngdll, ndp);D1 = [Desl;
124
Dr];for no = 1 : numNos
for m = 1 : 3g = 3*(no-1) + m;RA(g,:) = Pint(gdl(no,m),:);DT(g,:) = D1(gdl(no,m),:);
endendreturn%---------------------------------------------------------------------------% análise dinâmica através do método defunction [Pinter, Desl, ti] = newmark(nGrLibEst, numNos, numElems, coords, restrs,...
cargasNos, massaNos, conect, dadosElem, tipoLin_NL,lambda, ndp, tol, ngdll, gdl, P, delta_t, tempo)
% monta a matriz de amortecimentoCt = zeros(nGrLibEst,nGrLibEst);% obtém a matriz de massa da estruturaM = obtemMatMassaEstrutura(nGrLibEst, numNos, numElems, coords, massaNos, conect,
dadosElem, gdl);Mt = M;M = Mt(1:ngdll, 1:ngdll);C = Ct(1:ngdll, 1:ngdll);P = P(1:ngdll,1);% dados para o problema 15.1 do livro do A.K.Chopra%M = [1 0;% 0 1];%C = [ 0 0;% 0 0];%P = [60.98;% -55.02];% kt = [2.895 0;% 0 104.4];% ngdll = 2;%Características do método da aceleração médiagama=1/2;beta=1/4;% Características do método da aceleração linear% gama=1/2;% beta=1/6;alfa = (gama/(beta*delta_t))*C + 1/(beta*delta_t^2)*M;% matriz de amortecimentoA=(1/(beta*delta_t))*M + (gama/beta)*C;B=(0.5/beta)*M + delta_t*(0.5*gama/beta - 1)*C;P0 = zeros(ngdll,1);u0 = zeros(ngdll,1);v0 = zeros(ngdll,1);a0 = zeros(ngdll,1);
125
ui=u0;vi=v0;ai=a0;nt = tempo/delta_t;Pi = P0;t0 = 0;for i=1:nt
t = t0 + delta_t;Pi1 = P*lambda(i);lambda_atual = lambda(i);delta_Pi = Pi1 - Pi;delta_Pci = delta_Pi + A*vi + B*ai;
% para o caso não linearswitch tipoLin_NLcase 1[Pint, delta_ui] = ObtemDeslocPintDinamico(nGrLibEst, numNos, numElems, coords,...
restrs, cargasNos,massaNos, conect,...dadosElem, tipoLin_NL, lambda_atual,ndp,..tol, ngdll, gdl, P, delta_t, delta_Pci,...ui, alfa)
case 2[Pint, delta_ui] = NewtonRaphsonDinamico(nGrLibEst, numNos, numElems, coords,...
restrs, cargasNos,massaNos, conect,...dadosElem, tipoLin_NL, lambda_atual, ndp,...tol, ngdll, gdl, P, delta_t, delta_Pci,ui, alfa)
enddelta_vi = gama/(beta*delta_t)*delta_ui - gama/(beta)*vi + delta_t*(1-0.5*gama/beta)*ai;delta_ai = 1/(beta*delta_t^2)*delta_ui - 1/(beta*delta_t)*vi - (0.5/beta)*ai;ui1 = ui + delta_ui;vi1 = vi + delta_vi;ai1 = ai + delta_ai;Pinter(:,i)= Pint;Desl(:,i)= ui1;V(:,i)= vi1;Ac(:,i)= ai1;Pi = Pi1;ui=ui1;vi=vi1;ai=ai1;ti(i,1)=t;t0 = t;i = i+1;endreturn%---------------------------------------------------------------------------% emprego do método de Newton-Raphson para o caso Estático
126
function [Pinter, Desl, K] = NewtonRaphsonEstatico(nGrLibEst, numNos, numElems,coords,restrs, cargasNos,massaNos, conect,...dadosElem, tipoLin_NL,lambda, ndp, tol,...ngdll, gdl, P);
D0 = zeros (ngdll,1);Pi = P (1:ngdll,1);numInt = 40;Dr = zeros (nGrLibEst-ngdll,1);Pr = P (ngdll+1:nGrLibEst);for n = 1: ndplambda_atual = lambda(n);Pe = lambda_atual*Pi;i = 0;erro = 1e20;while (erro > tol)D = [D0;
Dr];[Pint,K] = obtemForcaMatRigEstrutura(nGrLibEst, numNos, numElems, coords, conect,...
dadosElem, tipoLin_NL, lambda_atual, gdl, D);% particiona a matriz de rigidez e o vetor de forças internasKm = K (1:ngdll,1:ngdll);Pii = Pint (1:ngdll,1);dPm = Pe - Pii;dD = Km\dPm;Dl = D0 + dD;D0 = Dl;i = i + 1;erro = norm(dD)if (i > numInt)disp ('Numero máximo de iterações ultrapassado');breakendendfprintf ('\n numero de interações: %d', i); % armazena os deslocamentos referentes a todos os passos de carga em uma única matrizDesl(:,n) = D0;Pinter(:,n) = Pint;endreturn%---------------------------------------------------------------------------% emprego do método de Newton-Raphson para o caso Dinâmicofunction [Pint, delta_ui] = NewtonRaphsonDinamico(nGrLibEst, numNos, numElems,
coords,restrs, cargasNos, massaNos, conect,...dadosElem, tipoLin_NL, lambda_atual,...ndp,tol,ngdll, gdl, P, delta_t,...delta_Pci,ui, alfa);
uj_1 = ui;
127
delta_Rj = delta_Pci;TOL = 1e-8;err = 1;j = 1;numIters = 40;Dr = zeros (nGrLibEst-ngdll,1);D = [ui
Dr];while (err > TOL)[Pint,K] = obtemForcaMatRigEstrutura(nGrLibEst, numNos, numElems, coords,conect,...
dadosElem, tipoLin_NL, lambda_atual, gdl, D);fsj_1 = Pint(1:ngdll);ktj_1 = K(1:ngdll, 1:ngdll);kcj_1 = ktj_1 + alfa;delta_uj = kcj_1\delta_Rj;uj = uj_1 + delta_uj;D = [uj
Dr];[Pint,K] = obtemForcaMatRigEstrutura(nGrLibEst, numNos, numElems, coords,conect,...
dadosElem, tipoLin_NL, lambda_atual, gdl, D);fsj = Pint(1:ngdll);delta_fj = fsj - fsj_1 + alfa * delta_uj;deltaRj1 = delta_Rj - delta_fj;err = norm (delta_uj)j = j + 1delta_Rj = deltaRj1;D = [uj
Dr];uj_1 = uj;if (j > numIters)disp ('Newton-Raphson ultrapassou o num. máximo de iterações');break;endendj;ui1 = uj;delta_ui = ui1 - ui;return%---------------------------------------------------------------------------% emprego do método de Newton-Raphson para o caso Dinâmicofunction [Pint, delta_ui] = ObtemDeslocPintDinamico(nGrLibEst, numNos, numElems,
coords,restrs, cargasNos,massaNos, conect,...dadosElem, tipoLin_NL, lambda_atual, ndp,...tol, ngdll, gdl, P, delta_t, delta_Pci,ui, alfa)
Dr = zeros (nGrLibEst-ngdll,1);D = [ui
Dr];
128
[Pint,K] = obtemForcaMatRigEstrutura(nGrLibEst, numNos, numElems, coords,conect,...dadosElem, tipoLin_NL, lambda_atual, gdl, D);
kt = K(1:ngdll, 1:ngdll);% kt = [2.895 0;% 0 104.4];kc = kt + alfa;delta_ui = kc\delta_Pci;return%---------------------------------------------------------------------------% calcula o vetor de forcas resistentes e a matriz de rigidez global da estruturafunction [Pint,K] = obtemForcaMatRigEstrutura(nGrLibEst, numNos, numElems, coords,
conect,dadosElem, tipoLin_NL,lambda_atual, gdl, D);% zeros (numLinhas, numColunas)K = zeros (nGrLibEst, nGrLibEst);Pint = zeros (nGrLibEst,1);for el = 1: numElems% determina os nos do elementonoI = conect(el,1);noJ = conect(el,2);% determina os graus de liberdade do elementogrLibElem = obtemGrausLibElemento(noI, noJ, gdl);% extrai os deslocamentos das barras, do vetor de deslocamentos da estruturafor i = 1: 6
dg(i,1) = D(grLibElem(i));end% determina as coordenadas dos nos do elementocoordNoI = coords(noI,:); % copia a linha 'noI' da matriz coords na matriz coordNoIcoordNoJ = coords(noJ,:); % copia a linha 'noJ' da matriz coords na matriz coordNoJ% calcula a matriz de rigidez do elemento em coordenadas globais[pg,kg] = obtemForcaMatRigGlobElem(dadosElem(el,:), coordNoI, coordNoJ,
lambda_atual, dg, tipoLin_NL);% adiciona matrix de rigidez do elemento a matriz de rigidez da estruturafor l = 1: 6L = grLibElem(l);
for c = 1: 6C = grLibElem(c);K(L, C) = K(L, C) + kg(l,c);
endend% adiciona forcas internas do elemento ao vetor forcas resistentes da estruturafor l = 1: 6
L = grLibElem(l);Pint(L,1) = Pint(L,1) + pg(l,1);
endendreturn%---------------------------------------------------------------------------
129
% calcula o vetor de forcas resistentes e a matriz de rigidez global da estruturafunction M = obtemMatMassaEstrutura(nGrLibEst, numNos, numElems, coords, massaNos,
conect, dadosElem, gdl);% zeros (numLinhas, numColunas)M = zeros (nGrLibEst, nGrLibEst);% monta a matriz de massa nodalfor n = 1 : numNos
for m = 1 : 3i = gdl(n,m);M(i,i) = massaNos(n,m);
endend% monta a matriz de massa consistente da estruturafor el = 1: numElems% determina os nos do elementonoI = conect(el,1);noJ = conect(el,2);% determina os graus de liberdade do elementogrLibElem = obtemGrausLibElemento(noI, noJ, gdl);% determina as coordenadas dos nos do elementocoordNoI = coords(noI,:); % copia a linha 'noI' da matriz coords na matriz coordNoIcoordNoJ = coords(noJ,:); % copia a linha 'noJ' da matriz coords na matriz coordNoJ% calcula a matriz de rigidez do elemento em coordenadas globaismg = obtemMatMassaGlobElem(dadosElem(el,:), coordNoI, coordNoJ);% adiciona matrix de rigidez do elemento a matriz de rigidez da estruturafor l = 1: 6L = grLibElem(l);
for c = 1: 6C = grLibElem(c);M(L, C) = M(L, C) + mg(l,c);
endendendreturn%---------------------------------------------------------------------------% obtém os graus de liberdade da estruturafunction [ngdll, gdl] = obtemGrausLibEstrut(numNos, restrs);% ngdll = numero de graus de liberdade livre da estrutura% gdl = numero de graus de liberdade da estrutura% cria um vetor com as restricoes nodaiscon=0;for no = 1 : numNos
for i=1:3if (restrs(no,i) == 0)con = con + 1;gdl(no,i) = con;end
130
endendngdll=con;for no = 1 : numNos
for i = 1 : 3if (restrs(no,i) == 1)
con = con + 1;gdl(no,i) = con;
endend
endreturn%---------------------------------------------------------------------------% determina os graus de liberdade do elementofunction grLibElem = obtemGrausLibElemento(noI, noJ, gdl)grLibElem(1:3) = gdl(noI,:);grLibElem(4:6) = gdl(noJ,:);return%---------------------------------------------------------------------------% calcula a força e a matriz de rigidez do elemento em coordenadas globaisfunction [pg,kg] = obtemForcaMatRigGlobElem(dadosElem, coordNoI, coordNoJ,...
lambda_atual, dg, tipoLin_NL);switch tipoLin_NLcase 1[pg,kg] = obtemForcaMatRigGlobElemLinear(dadosElem, coordNoI, coordNoJ,
lambda_atual,dg);case 2[pg,kg] = obtemForcaMatRigGlobElemNL(dadosElem, coordNoI, coordNoJ,...
lambda_atual, dg);endreturn%---------------------------------------------------------------------------% calcula a força e a matriz de rigidez do elemento em coordenadas globaisfunction [pg,kg] = obtemForcaMatRigGlobElemLinear(dadosElem, coordNoI, coordNoJ,...
lambda_atual, dg);E = dadosElem(1);A = dadosElem(2);I = dadosElem(3);m = dadosElem(4);qx = dadosElem(5) * lambda_atual;qy = dadosElem(6) * lambda_atual;[L, cosAlfa, senAlfa] = compOrientElem(coordNoI, coordNoJ);EA = E * A;EI = E * I;% 1 2 3 4 5 6kl = [ EA/L 0 0 -EA/L 0 0; % 1
0 12*EI/L^3 6*EI/L^2 0 -12*EI/L^3 6*E*I/L^2; % 2
131
0 6*EI/L^2 4*EI/L 0 -6*EI/L^2 2*E*I/L; % 3-EA/L 0 0 EA/L 0 0; % 40 -12*EI/L^3 -6*EI/L^2 0 12*EI/L^3 -6*E*I/L^2; % 50 6*EI/L^2 2*EI/L 0 -6*EI/L^2 4*E*I/L]; % 6
% carga nodal equivalentepeql = [qx*L/2;
qy*L/2;qy*(L^2)/12;qx*L/2;qy*L/2;qy*(L^2)/12];
% transforma para coordenadas globaisR = obtemMatRotacao(cosAlfa, senAlfa);kg = R' * kl * R;peqg = R' * peql;pg = kg * dg - peqg;return%---------------------------------------------------------------------------% calcula a matriz de massa do elemento em coordenadas globaisfunction mg = obtemMatMassaGlobElem(dadosElem, coordNoI, coordNoJ);m = dadosElem(4);[L, cosAlfa, senAlfa] = compOrientElem(coordNoI, coordNoJ);ml = obtemMatMasConLocElem(m,L);R = obtemMatRotacao(cosAlfa, senAlfa);mg = R'*ml*R;return%---------------------------------------------------------------------------% calcula a matriz de massa consistente do elemento em coordenadas locaisfunction ml = obtemMatMasConLocElem(m,L);% 1 2 3 4 5 6ml =(m*L/420)* [ 140 0 0 70 0 0; % 1
0 156 22*L 0 54 -13*L; % 20 22*L 4*L^2 0 13*L -3*L^2; % 370 0 0 140 0 0; % 40 54 13*L 0 156 -22*L; % 50 -13*L -3*L^2 0 -22*L 4*L^2]; % 6
return%---------------------------------------------------------------------------% calcula o comprimento e orientação dos elementosfunction [L, cosAlfa, senAlfa] = compOrientElem(coordNoI, coordNoJ);XJ = coordNoJ';XI = coordNoI';XIJ = XJ - XI;L = norm (XIJ);e1L = XIJ/L;cosAlfa = e1L(1,1);senAlfa = e1L(2,1);
132
return%---------------------------------------------------------------------------% calcula a matriz de rotação do elemento (do sistema global para local)function R = obtemMatRotacao(cost, sint)c = cost;s = sint;R = [ c s 0 0 0 0;
-s c 0 0 0 0;0 0 1 0 0 0;0 0 0 c s 0;0 0 0 -s c 0;0 0 0 0 0 1];
return%---------------------------------------------------------------------------% calcula a matriz de rotacao do elemento (do sistema global para local)function [freq, Modos] = AnaliseModal(nGrLibEst, ngdll, K, M)Kll = K(1:ngdll,1:ngdll);Krl = K(ngdll+1:nGrLibEst,1:ngdll);Mll = M(1:ngdll,1:ngdll);Mrl = K(ngdll+1:nGrLibEst,1:ngdll);% determinacao dos autovalores e autovetores[Modos,W2] = eig(Kll,Mll);W = sqrt(W2);freq = diag(sqrt(W2))/(2*pi);return%---------------------------------------------------------------------------
133
Anexo 2 - Programa para Cálculo da Matriz de Rigidez do Elemento em Coordenadas
Globais
%---------------------------------------------------------------------------% calcula a força e a matriz de rigidez do elemento em coordenadas globais
function [pg,kg] = obtemForcaMatRigGlobElemNaoLinear(dadosElem, coordNoI,coordNoJ, lambda_atual, dg);
E = dadosElem(1);A = dadosElem(2);I = dadosElem(3);m = dadosElem(4);qx = dadosElem(5) * lambda_atual;qy = dadosElem(6) * lambda_atual;XJ = coordNoJ';XI = coordNoI';XIJ = XJ - XI;L = norm (XIJ);
e1L = XIJ/L;cosAlfaL = e1L(1,1);senAlfaL = e1L(2,1);UI = dg(1:2,1);gamaI = dg(3,1);UJ = dg(4:5,1);gamaJ = dg(6,1);UIJ = UJ - UI;l = norm(XIJ + UIJ);alfaL = atan2(senAlfaL,cosAlfaL);e1B = (XIJ+UIJ)/l;cosAlfaB = e1B(1,1);senAlfaB = e1B(2,1);alfaB = atan2(senAlfaB,cosAlfaB);db(1,1) = (1/(l+L))*(2*XIJ+UIJ)'*UIJ;db(2,1) = gamaI + alfaL - alfaB;db(3,1) = gamaJ + alfaL - alfaB;EA = E * A;EI = E * I;% calcula matriz de rigidez no sistgema básicokb = [EA/L 0 0;
0 4*EI/L 2*EI/L;0 2*EI/L 4*EI/L];
%cosAlfa = cosAlfaL*cosAlfaB + senAlfaL*senAlfaB;%senAlfa = -cosAlfaL*senAlfaB + senAlfaL*cosAlfaB;alfa = alfaB-alfaL;
134
cosAlfa = cos(alfa);senAlfa = sin(alfa);qxb = ( qx*cosAlfa + qy*senAlfa);qyb = (-qx*senAlfa + qy*cosAlfa);peqb = [qxb*L/2;
qyb*(L^2)/12;qyb*(L^2)/12];
pb0 = kb*db;pb = pb0 - peqb;% calcula matriz de transformaçãoc = cosAlfaB;s = senAlfaB;T = [-c -s 0 c s 0;
-s/l c/l 1 s/l -c/l 0;-s/l c/l 0 s/l -c/l 1];
% monta matriz de rigidez geométricag1 = (1/l)*[ s^2 -c*s 0 -s^2 c*s 0;
-c*s c^2 0 c*s -c^2 0;0 0 0 0 0 0;-s^2 c*s 0 s^2 -c*s 0;c*s -c^2 0 -c*s c^2 0;0 0 0 0 0 0];
g2 = (1/l^2)*[-2*c*s c^2-s^2 0 2*c*s -c^2+s^2 0;c^2-s^2 2*c*s 0 -c^2+s^2 -2*c*s 0;0 0 0 0 0 0;2*c*s -c^2+s^2 0 -2*c*s c^2-s^2 0;-c^2+s^2 -2*c*s 0 c^2-s^2 2*c*s 0;0 0 0 0 0 0];
g3 = g2;prg = [-qxb*l*cosAlfaB + qyb*(l/2)*senAlfaB;
-qxb*l*senAlfaB - qyb*(l/2)*cosAlfaB;0;
qyb*(l/2)*senAlfaB;-qyb*(l/2)*cosAlfaB;
0];pg = T'*pb + prg;kgeo = g1*pb(1) + g2*pb(2) + g3*pb(3);kg1 = T'*kb*T + kgeo;a = (L/12)*[-6*senAlfa 6*cosAlfa;
-L*cosAlfa -L*senAlfa;L*cosAlfa L*senAlfa];
s2a = 2*senAlfa*cosAlfa;c2a = (cosAlfa)^2 - (senAlfa)^2;b = (L/2)*[ s2a -c2a;
-c2a -s2a; 0 0;-s2a c2a;
135
c2a s2a;0 0];
q = [qx;qy];
r = [-c -s 0 c s 0];S = [s -c 0 -s c 0];R = obtemMatRotacao(cosAlfaL, senAlfaL);k0 = (-(1/l)*T'*a*q+(1/l)*R'*b*q)*S;kg = kg1 + k0;return%---------------------------------------------------------------------------% calcula a matriz de rotacao do elemento (do sistema global para local)function R = obtemMatRotacao(cost, sint)c = cost;s = sint;R = [ c s 0 0 0 0;
-s c 0 0 0 0;0 0 1 0 0 0;0 0 0 c s 0;0 0 0 -s c 0;0 0 0 0 0 1];
return
136
Anexo 3 - Programa para cálculo por diferenças finitas da matriz de rigidez
function matrignum2
% funcao para calculo da matriz de rigidez numericamenteE = 2e12;b = 0.20; %(m)h = 0.02; %(m)A = b*h; %(m^2)I = (b*h^3)/12; %(m^4)qx = 0; %(kgf/m)qy = -1000; %(kgf/m)
dadosElem(1,:) = [E A I qx qy];L = 1.2;
coordNoI = [0 0];coordNoJ = [L 0];
lambda_atual = 1;delta = 1e-10;
%dg1 = [0 0 0 0 0 0]';dg1 = [2.2 1.3 -5.2 3.5 6.3 4.8]'
[pg1,kg1] = obtemForcaMatRigGlobElemNaoLinear(dadosElem, coordNoI, coordNoJ,... lambda_atual, dg1);
for i=1:6
dg2 = dg1; dg2(i) = dg1(i) + delta
[pg2,kg2] = obtemForcaMatRigGlobElemNaoLinear(dadosElem, coordNoI, coordNoJ,... lambda_atual, dg2); deltaPg = (pg2 - pg1)/delta;
knum(:,i) = deltaPg;
%dg1 = dg2end;kg1knum(knum-kg1)./kg1
return
137
Anexo 4 - Programa para Entrada de Dados 1
% Universidade Federal do Para - UFPa% Centro Tecnológico - CT% Departamento de Engenharia Mecânica - DEM% Curso de Mestrado em Engenharia Mecânica% Dissertação: Exemplo Aplicado a Análise Dinâmica Não-linear de Vigas% Autores: Remo Magalhães de Souza ([email protected])% Gandhy Y R Aranha Júnior ([email protected])% Versão 1.0 Data: 10/2003% Arquivo exemplo para entrada de dados - inpYang1984vigaCompleta.m
function inpYang1984vigaCompleta% dados de entradanumElems = 10; % numero de elementos da estruturanumNos = numElems+1; % numero de nos da estruturadelta_t = 200e-6; % intervalo de acréscimo do tempotempo = 50000e-6; % tempo de análiseL = 20; % comprimento total da vigandp = tempo/delta_t; % número de passos de cargatol = 1e-12;% coordenadas dos nos% X Ydx = L/numElems;coords = zeros(numNos,2);for i = 1:numNos % determina coords.coords(i,1) = dx*(i-1); % coord. xcoords(i,2) = 0; % coord. yend% restricoes nodais (condições de apoio)restrs = zeros(numNos,3);NU = numElems+1; % último nó% Dx Dy Rzrestrs(1,:) = [1 1 1]; % no 1restrs(NU,:) = [1 1 1]; % no último nó% conectividade dos elementos% noI noJfor i = 1: numElemsconect(i,1) = i;conect(i,2) = i+1;end% cargas nodais% Fx Fy FznoCentral = numElems/2 + 1;cargasNos = zeros (numNos,3);
138
cargasNos(noCentral,:) = [0 640 0]; % no 6% Características geométricas e do materialrho = 0.0026 % lb-sec^2/in]b = 1;h = 1/8;A = b*h;I = b*h^3/12;% massas nodaismassaNos = zeros(numNos,3);% monta o vetor lambida (variação da força externa)i=1;lambda(1,1) = 0;for t = delta_t: delta_t: tempoi=i+1;lamb = 1;lambda(1,i) = lamb;end% propriedades da secaom=rho; % massa distribuída do elemento% Modulo de deformacao longitudinalE = 3e7;% carga distribuídaqx = 0;qy = 0;dadosElem = zeros(numElems,6);for i = 1:numElemsdadosElem(i,:) = [E A I m qx qy];endtipoEst_Din = 3;tipoLin_NL = 2;close all;% desenha estruturafigure(1);DesenhaEstrut2d(numNos, numElems, coords, conect);% chama o arquivo para analise estrutural[DT,RA,ti,freq, Modos] = AnalisePort01(numNos, numElems, coords, restrs, cargasNos,...
massaNos, conect, dadosElem, tipoEst_Din,...tipoLin_NL,lambda,delta_t,ndp, tol, tempo)
d = DT(3*noCentral-1,:)tipoEst_Din = 3;tipoLin_NL = 1;[DT,RA,ti,freq, Modos] = AnalisePort01(numNos, numElems, coords, restrs, cargasNos,...
massaNos, conect, dadosElem, tipoEst_Din,...tipoLin_NL,lambda,delta_t,ndp, tol, tempo)
d2 = DT(3*noCentral-1,:)d3 = DT(3*noCentral-1,:)J = delta_t: delta_t: tempo;
139
G = [0.00020442 0.00883440.0006177 0.0353910.0012332 0.0663380.0018531 0.106170.0024752 0.150450.0028974 0.194780.0033173 0.234670.0035395 0.279060.0043639 0.327730.0047838 0.367610.005006 0.4120.0054304 0.460770.0060503 0.500610.006668 0.5360.0074857 0.571340.0081034 0.606730.0089210 0.642070.0095365 0.673010.010154 0.708400.010372 0.743900.011385 0.77030.01197 0.739030.011956 0.712370.012141 0.681210.012521 0.641110.012501 0.601110.012876 0.552120.012854 0.507680.013434 0.467530.013816 0.431870.014198 0.396210.014776 0.351620.015156 0.311510.015540 0.280300.015922 0.244640.017109 0.217660.018504 0.20840.0195 0.199240.020088 0.176870.020875 0.149990.021866 0.131950.022679 0.15840.023297 0.193790.023712 0.224790.024128 0.255790.024352 0.304620.024575 0.34900
140
0.024997 0.393330.025223 0.446610.025239 0.477710.025463 0.526540.025885 0.570870.026508 0.615150.026528 0.655150.027147 0.694980.027759 0.717040.028545 0.690170.028932 0.66340.029318 0.636630.029705 0.609860.031098 0.596150.031887 0.573720.032273 0.546950.032856 0.511240.033238 0.475580.033618 0.435480.034202 0.404210.034387 0.373050.034567 0.333000.034746 0.292960.034733 0.266290.035118 0.235080.035497 0.194980.035884 0.168210.036069 0.137050.036457 0.114720.037239 0.0789580.038639 0.0785850.039646 0.0916500.040061 0.122650.040481 0.162540.040699 0.198040.041117 0.233480.041532 0.264480.041748 0.295530.042363 0.326480.043379 0.357320.043999 0.397160.044416 0.432600.045041 0.481320.045465 0.53010.045883 0.565540.046303 0.605430.046527 0.65426
141
0.046552 0.70314];figure(1);%desenha a estruturaDesenhaEstrut2d(numNos, numElems, coords, conect);hold on;figure(2);plot(J,d2,'k-*',J,d3,'b-o')xlabel ('Tempo')ylabel ('Deslocamanto Vertical (in)')title('Desloc. Vertical no meio do vão - Caso Linear')legend ('Formulação Co-rotacional','Mondkar e Powell ')gridzoom on;figure(3);plot(J,d,'k-*',G(:,1),G(:,2),'r-o')xlabel ('Tempo (s)')ylabel ('Deslocamanto Vertical (in)')title('Desloc. Vertical no meio do vão - Caso Não LInear')gridlegend ('Elemento Proposto','Yang e Saigal')figure(4);plot(J,d2,'k-*',J,10*d,'b-o')title('Desloc. Vertical no meio do vão')xlabel ('Tempo')ylabel ('Caso Linear x 1 - Caso Não Linear x 10')gridlegend ('Caso Linear','Caso Não-linear')zoom on;
142
Anexo 4 - Programa para Entrada de Dados 2
% Universidade Federal do Para - UFPa% Centro Tecnologico - CT% Departamento de Engenharia Mecânica - DEM% Curso de Mestrado em Engenharia Mecânica% Dissertação: Exemplo Aplicado a Análise Dinâmica Linear de Vigas% Programa Didatico para Analise Não Linear Estática de Porticos Planos% Autores: Remo Magalhaes de Souza ([email protected])% Gandhy Y R Aranha Júnior ([email protected])% Versão 1.0 Data: 10/2003% Arquivo exemplo para entrada de dados - inpVigaBiapoiadaChopra2.m
function inpVigaBiapoiadaChopra2
% dados de entrada
numElems = 30; % numero de elementos da estruturanumNos = numElems+1; % numero de nos da estruturadelta_t = 0.5e-3; % intervalo de acréscimo do tempotempo = 300e-3; % tempo de análiseL = 20; % comprimento total da viga em inndp = tempo/delta_t; % número de passos de carga
tol = 1e-12;
% coordenadas dos nos% X Ydx = L/numElems;coords = zeros(numNos,2); for i = 1:numNos % determina coords. coords(i,1) = dx*(i-1); % coord. x coords(i,2) = 0; % coord. y end
% restricoes nodais (condicoes de apoio)restrs = zeros(numNos,3);NU = numElems+1; % último nó% Dx Dy Rzrestrs(1,:) = [1 1 0]; % no 1restrs(NU,:) = [1 1 0]; % no último nó
% conectividade dos elementos
143
% noI noJfor i = 1: numElems conect(i,1) = i; conect(i,2) = i+1; end
% cargas nodais % Fx Fy FzP = 640; % Força em lbnoCentral = numElems/2 + 1;cargasNos = zeros (numNos,3);cargasNos(noCentral,:) = [0 P 0]; % no 6
% Características geométricas e do materialrho = 0.0026 % massa em lb.s^2/inb = 1; % base da viga em inh = b/8; % altura da viga em inA = b*h;I = b*h^3/12;
% massas nodaismassaNos = zeros(numNos,3);
% monta o vetor lambida (variação da força externa)i=1;lambda(1,1) = 0;for t = delta_t: delta_t: tempo i=i+1; lamb = 1; lambda(1,i) = lamb;end
% propriedades da secaom=rho/L; % massa distribuída do elemento% Modulo de deformacao longitudinalE = 3e7; % módulo de elasticidade do material em N/m^2% carga distribuídaqx = 0;qy = 0;
dadosElem = zeros(numElems,6);for i = 1:numElems dadosElem(i,:) = [E A I m qx qy];end
144
tipoEst_Din = 3;tipoLin_NL = 1;
close all;
% desenha estruturafigure(1);DesenhaEstrut2d(numNos, numElems, coords, conect);
% chama o arquivo para analise estrutural[DT,RA,ti,freq, Modos] = AnalisePort01(numNos, numElems, coords, restrs, cargasNos,... massaNos, conect, dadosElem, tipoEst_Din,... tipoLin_NL,lambda,delta_t,ndp, tol, tempo)
d = DT(3*noCentral-1,:)
% solucao analiticact = 2*P*L^3/((pi^4)*E*I);
J = delta_t: delta_t: tempo;
phi = [1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1];
soma = 0;for n =1:20 wn = ((n*pi/L)^2)*sqrt(E*I/m); soma = soma + (phi(n)/n^4)*(1-cos(wn*J))*sin(n*pi/2);endd2 = ct*soma;
figure(1);%desenha a estruturaDesenhaEstrut2d(numNos, numElems, coords, conect);
hold on;figure(2);plot(J,d2,'b-')title('Desloc. Vertical no meio do vão - Analítico Chopra')xlabel ('Tempo')ylabel ('Deslocamento Vertical')gridzoom on;
145
figure(3);plot(J,d,'k-*')title('Desloc. Vertical no meio do vão - Elemento Proposto')xlabel ('Tempo')ylabel ('Deslocamento Vertical')
grid
figure(4);plot(J,d2,'b-o',J,d,'k-*')title('Desloc. Vertical no meio do vão')xlabel ('Tempo')ylabel ('Deslocamento Vertical (in)')grid
legend ('Analítico Chopra','Elemento Proposto')zoom on;
